# Combining bootstrap and maximum a posteriori estimation

I've recently read interesting paper Uncertainty in Neural Networks: Approximately Bayesian Ensembling by Pearce et al (2020), who suggested algorithm for approximating the posterior distribution by taking sample from the prior $$\theta_{acc} \sim p(\theta)$$ and finding MAP estimate using

$$\theta_\text{MAP}' = \underset{\theta}{\operatorname{arg\,max}} \; p(x | \theta) \; p_{acc}(\theta)$$

where $$p_{acc}(\theta)$$ is the prior distribution centred at the $$\theta_{acc}$$ value. This is repeated many times, and we treat the ensemble of such estimates as samples from the posterior distribution.

Unfortunately, the above approach underestimates the variability of the parameters. What if instead we varied the data, rather then priors? We could take bootstrap sample $$x^*$$ from the data $$x$$ and then estimate

$$\theta_\text{MAP}^* = \underset{\theta}{\operatorname{arg\,max}} \; p(x^* | \theta) \; p(\theta)$$

again, repeating this many times, and ensembling the results. My questions:

• Do you recall anyone describing such approach? So far, the closest I found is the paper by Laird and Louis (1985) about combining empirical Bayes with bootstrap.
• What could be the potential problems with this approach? Why could it, or could not work?

Doing this for a trivial example of estimating mean of normal distribution with known variance gives pretty nice results.

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sp

np.random.seed(42)

n = 15
# true parameters
μ = 5
σ = 2.7

x = sp.norm(μ, σ).rvs(n)

prior_μ = 0
prior_σ = 10
prior_dist = sp.norm(prior_μ, prior_σ)

def posterior(x, prior_μ, prior_σ):
n = len(x)
σ2 = σ ** 2
prior_σ2 = prior_σ ** 2

post_σ2 = 1 / (1/prior_σ2 + n/σ2)
post_μ = post_σ2 * (prior_μ/prior_σ2 + np.sum(x)/σ2)

return float(post_μ), np.sqrt(post_σ2)

post_μ, post_σ = posterior(x, prior_μ, prior_σ)
post_dist = sp.norm(post_μ, post_σ)

R = 5000
results = []

for _ in range(R):
idx = np.random.choice(n, n, replace=True)
results.append(posterior(x[idx], prior_μ, prior_σ))

results = np.vstack(results)

np.mean(results[:, 0]), np.std(results[:, 0])
## (4.9949681628669085, 0.6582492996849297)

# vs true posterior
post_μ, post_σ
## (5.003623405947105, 0.6954491092861294)


Trying to partially answer my own question, there are some examples of using regular bootstrap to approximate Bayesian inference e.g. Laird and Louis (1987) used empirical Bayes together with bootstrap to obtain confidence intervals, Djuric (1997) used MAP with bootstrap for model selection, and to improve MCMC sampling (1999), Osband and Van Roy (2018) used bootstrap together with Thompson sampling for multiarmed bandits problem, while Osband et al (2018) used it in reinforcement learning scenario (see figure below), and Everitt (2017) and Picchini and Everitt (2018) (see also this talk) described using bootstrap in Approximate Bayesian Computation.

Unfortunately, neither of those references described just using bootstrap with MAP, but apparently there are some promising results in related areas.