When approximating a posterior using MCMC, why don't we save the posterior probabilities but use the parameter value frequencies afterwards?

Question

I'm currently estimating parameters of a model defined by several ordinary differential equations (ODEs). I try this with a bayesian approach by approximating the posterior distribution of the parameters given some data using Markov Chain Monte Carlo (MCMC).

A MCMC sampler generates a chain of parameter values where it uses the (unnormalized) posterior probability of a certain parameter value to decide (stochastically) whether it will add that value to the chain or add the previous value again. But, it seems to be the practice that the actual posterior probabilities do not need to be saved, rather is a n-dimensional histogram of the resulting parameter values generated and summary statistics like highest density regions (HDRs) of a parameter posterior ditribution is calculated from this histogram. At least that's what I think i learned from Kruschkes tutorial book on bayesian inference.

My Question: Wouldn't it be more straightforward to save the posterior probabilities of the sampled parameter values along with these and approximate the posterior distribution from these values and not from the frequencies of parameter values in the MCMC chain? The problem of the burn-in phase would not arise as the sampler would initially still sample low probability regions more often than they would "deserve" by their posterior probabilities but it would not be anymore the problem of giving unduly high probability values to these.

Answer 1

This is an interesting question, with different issues:

1. MCMC algorithms do not always recycle the computation of the posterior density at all proposed values, but some variance reduction techniques like Rao-Blackwellisation do. For instance, in a 1996 Biometrika paper with George Casella, we propose to use all simulated values, $\theta_i$ $(i=1,\ldots,T)$, accepted or not, by introducing weights $\omega_i$ that turn the average$$\sum_{i=1}^T \omega_ih(\theta_i)\big/\sum_{i=1}^T \omega_i$$into an almost unbiased estimator. (The almost being due to the normalisation by the sum of the weights.)
2. MCMC is often used on problems of large (parameter) dimension. Proposing an approximation to the whole posterior based on the observed density values at some parameter values is quite a challenge, including the issue of the normalising constant mentioned in Tim's answer and comments. One can imagine an approach that is a mix of non-parametric kernel estimation (as in e.g. krigging) and regression, but the experts I discussed with about this solution [a few years ago] were quite skeptical. The issue is that the resulting estimator remains non-parametric and hence "enjoys" non-parametric convergence speeds that are slower than Monte Carlo convergence speeds, the worse the larger the dimension.
3. Another potential use of the availability of the posterior values $\pi(\theta|\mathcal{D})$ is to weight each simulated value by its associated posterior, as in $$\frac{1}{T}\sum_{t=1}^T h(\theta_t) \pi(\theta_t|\mathcal{D})$$ Unfortunately, this creates a bias as the simulated values are already simulated from the posterior: $$\mathbb{E}[h(\theta_t) \pi(\theta_t|\mathcal{D})]=\int h(\theta) h(\theta_t) \pi(\theta_t|\mathcal{D})^2\text{d}\theta$$ Even without a normalisation issue, those simulations should thus be targeting $\pi(\theta|\mathcal{D})^{1/2}$ and use a weight proportional to $\pi(\theta|\mathcal{D})^{1/2}$ but I do not know of results advocating this switch of target. As you mention in the comments, this is connected with tempering in that all simulations produced within a simulated tempering cycle can be recycled for Monte Carlo (integration) purposes this way. A numerical issue, however, is to handle several importance functions of the form $\pi(\theta)^{1/T}$ with missing normalising constants.

Answer 2

As you correctly noticed, the probabilities we are dealing with are unnormalized. Basically, we use MCMC to compute the normalizing factor in Bayes theorem. We cannot use the probabilities because they are unnormalized. The procedure that you suggest: to save the unnormalized probabilities and then divide them by their sum is incorrect.

Let me show it to you by example. Imagine that you used Monte Carlo to draw ten values from Bernoulli distribution parametrized by $p=0.9$, they are as follows:

1 0 1 1 1 1 1 1 1 1

you also have corresponding probabilities:

0.9 0.1 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9

In this case the probabilities are normalized, but dividing them by their sum (that by axioms of probability is equal to unity) should not change anything. Unfortunatelly, using your procedure it does change the results to:

> f/sum(f)
 [1] 0.10975610 0.01219512 0.10975610 0.10975610 0.10975610 0.10975610 0.10975610 0.10975610 0.10975610 0.10975610

Why is that? The answer is simple, in your sample each saved "probability" f appears with probability f, so you are weighting the probabilities by themselves!