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

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.

• If you can calculate the posterior probabilities without using MCMC (so to save them), then why would you like to use it? – Tim Sep 1 '16 at 9:50
• Because I need MCMC to gain efficiency. If I just laid a grid over the parameter space and calculated unnormalized posterior probabilities for all resulting parameter values, I would waste many time on low-probability regions. Being able to get unnormalized posterior probability values for a given parameter value is a precondition for using MCMC. I need not be able to solve the posterior analytically. So I could take all saved probability values, divide them by their sum and the result would be an approximation of my posterior. – akraf Sep 1 '16 at 9:57
• @Tim: what he means is that for computing the acceptance probability of the proposed move, you evaluate the posterior at the current state and at the proposed state. If you keep these posterior values for every state reached then the OP thinks you can derive the whole posterior but that is not the case, at least I have never seen a theorem that proves that. By looking at the distribution of the states reached markov theory shows that you get a sample from the posterior ' in the end' – user83346 Sep 1 '16 at 10:39
• @fcop yes, I understand that and I think we are saying the same thing but in different words :) – Tim Sep 1 '16 at 10:41

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.
• Thank you for your extensive comments, please allow me some clarifying questions! I do not understand what you mean by "recycle" in your point 1 and how that is preventing the use of the unnormalized posterior values. To point 2: If the "approximation to the whole posterior based on the observed density values at some parameter values is quite a challenge", why is it less so if using just the frequencies of the samples resulting from the MCMC process? – akraf Nov 22 '16 at 10:14
• To point 3: Targeting $\pi(\theta|D)^{1/T}$ with $T>1$ is a common way of "tempering" the posterior, i.e. "flattening its peaks" to facilitate the mixing of the MCMC chains, where the flattening is stronger, the higher $T$ is. Might the approach you suggested be a way to recover the original untempered distribution $\pi(\theta|D)$, given samples from the tempered distribution $\pi(\theta|D)^{1/T}$? – akraf Nov 22 '16 at 10:16

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!