Monte Carlo simulation of posterior distribution

Question

I have a question regarding Bayesian estimation and whether a method such as MCMC is necessary. In the Bayesian setting we have that

$$\pi_{\Theta|X}(\theta|x)\propto f_{X|\Theta}(x|\theta)\,\pi_{\Theta}(\theta)$$

where $\pi_{\Theta}(\theta)$ is the prior density of parameters; $f_{X|\Theta}(x|\theta)$ is the density of the data given the parameters $\Theta=\theta$; and $\pi_{\Theta|X}(\theta|x)$ is the posterior density of interest.

Now, let's assume we know our prior density and are able to calculate the likelihood function of the data given a set of parameters. Further to this, assume we are merely interested in inferring the set of parameters that maximises the posterior density i.e. we want

$$\DeclareMathOperator*{\argmax}{argmax} \argmax_\theta\,\pi_{\Theta|X}(\theta|x)$$

If we follow the following procedure:

1. Randomly sample (Monte Carlo style) from the prior density such that we get a set of simulated $\theta$;
2. Evaluate the likelihood function at each of these simulated parameter values; and
3. Multiply the likelihood and the prior density.

Does taking the mode of this simulated posterior give us the set of parameters that maximise the posterior density? Because we are only interested in the set of parameters that maximise the posterior density, is MCMC unnecessary?

Furthermore, is it reasonable to state that the reason for MCMC is for situations when we would like to evaluate the full posterior density such that we can calculate some other function of the parameters (e.g. $\mathbb{E}[\theta]$) and not just the mode?

Answer 1

MCMC is not needed to find the posterior mode. I think the approach you propose is correct in principle, but unfeasible for all but the simplest tasks. What you suggest is a particularly naive optimization scheme for $\theta$, the problem is that

a) If the data (likelihood) provides non-negligible information, the posterior density will actually be very close to zero on most of the parameter values admitted by the prior. It is hard to guess the range where posterior has non-negligible density.

b) The posterior mode can be arbitrarily narrow, requiring you to run a lot of samples to have a chance of hitting close to it.

Of course this all gets exponentially worse as the dimension of $\theta$ grows.

Basically any classical optimization method (simulated annealing, gradient descent, ...) is likely to outperform your approach. In fact the Stan probabilistic language has a very fast optimization mode (with L-BFGS).

Note however that the posterior mode is a problematic quantity as it does change under reparametrization - unlike maximum likelihood or full Bayes.

You correctly state that MCMC is useful to evaluate the full posterior density (and expectations over it). In fact it is the only known practical way to get full posterior density in non-trivial models.

The following question also has some related thoughts: Bayesian posterior: is multiplying likelihood by prior (rather than simulation) an acceptable approach?

Answer 2

What you have described is essentially rejection sampling without the rejection. Since there is no guarantee that high density regions of the prior will align with high density regions of the posterior, if you stop sampling after generating $N$ variates, there's no guarantee you've sampled anywhere near the mode(s). If you instead reject poor variates ($f_{X|\Theta} < U$ for random uniform $U\in[0,1]$), your $N$ variants will be exact, but it might take a prohibitive amount of time if your prior is a bad estimate of the posterior. If your likelihood is log-concave, you can use adaptive rejection sampling for low dimensional $\theta$, but the curse of dimensionality will halt this method for large parameter spaces. To avoid these two sources of slowdown, MCMC methods are used. Sampling is generally used to evaluate expectations, but it is preferrable to use exact methods, if applicable.

Now if all you want is the MAP estimator, just maximize the (log of the) posterior likelihood with any non-linear optimizer. This will likely be more efficient than sampling anyway.