How do we reuse and store the output of an MCMC?

Question

For many Bayesian models, the posterior distribution is intractable... a solution is then to sample points from this unknow distribution with a Markov Chain Monte Carlo (MCMC). But at the end, how do we store the result of our sampling process, in other words, how do we predict a new value given the fact that we have only points that follow some distribution. My guess is that we have to approximate the distribution with for example a Mixture of Gaussians, we obtain this way an analytical expression of the distribution. Therefore we are able to predict a new output value by summing over parameters that follow this posterior distribution...

Another solution would be to store all the samples we have and then to compute the sum over all those samples ? But in this case we have to store millions of sample in case of complicated distributions and it doesn't look right to me ...

So the question is simple, which solution is employed ?

Thank you so much for your help !

Answer 1

I'll frame this in the context of linear regression.

Say your model is a relatively simple one. Something like

$$ \begin{align} \beta_0 &\sim \mathcal{N}(0,1)\\ \beta_1 &\sim \mathcal{N}(0,1)\\ y_i &\sim \mathcal{N}(\beta_0 + \beta_1 x_i, \sigma) \end{align}$$

Here, $\sigma$ is known. You then stumble upon a new input $x_n$ and would like to predict $y_n$.

There are a few ways to, as you say, how do we predict a new value given the fact that we have only points that follow some distribution. I'll outline only one here and assume we are working in a language like R.

Some pseudocode for how to compute this is as follows:

mean_samples<- beta_0_samples + beta_1_samples*x_n

prediction = mean(mean_samples)

And this totally makes sense in the context of Bayesianism. The mean of the likelihood is a random variable. It has expectation, variance, etc. Why not just approximate the mean of the distribution of the mean of the likelihood and use that as a prediction? You note that we have to store the samples in order to do this sort of computation, and you're right. But we don't usually need millions of samples. A few thousand might do fine, especially with new methods which find the typical set very quickly.

There are other ways to do this (e.g. MAP, median, etc), but the mean of the posterior samples is easiest to understand.

Answer 2

Expanding the answer by Demetri Pananos, recall that what we estimate is the posterior distribution of the parameters

$$ p(\theta|X) = \frac{ p(X|\theta) \; p(\theta) }{p(X)} $$

So we are not making the predictions at this stage. Unlike with point estimation, we end up in here with estimates of distributions for the parameters. If we had point estimates, to make prediction from the model, we would plug-in the estimated parameters and the data to our model (a function $f$ of data $X$ and parameters $\hat\theta$) and return the outputs as our prediction

$$ \hat y = f(X; \hat\theta) $$

Since we have distributions of the parameters, we plug-in the distributions and obtain posterior predictive distribution of the predicted values.

As you noticed, in many cases the posterior distributions are intracable and instead of finding the distributions, we use MCMC to obtain samples from those distributions. If we have large enough number of samples, we can treat empirical statistics from those samples as estimates from the posterior distribution, for example, to estimate expected value of $\theta$ you would take mean of the samples from the posterior distribution of $\theta$

$$ E[\theta|X] \approx n^{-1} \sum_{i=1}^n \hat\theta_i $$

where $\hat\theta_1,\hat\theta_2,\dots,\hat\theta_n$ are $n$ samples from the posterior distribution. To obtain posterior predictive distribution, you would take the samples of the parameters from the posterior distribution and plug-in them to the model function to obtain samples from the posterior predictive distribution

$$ \hat y_i = f(X; \hat\theta_i) $$

Alternatively, if you want to make predictions using some other data, that was not used for training, say data from the test set $X_\text{test}$, you plug it in as the same way

$$ \hat y_{\text{test},i} = f(X_\text{test}; \hat\theta_i) $$

Now, when you have those samples, you can estimate all the statistics from the posterior predictive distribution the same way as you would from the posterior distribution. To obtain point estimates, you can take things such as mean, median, or mode of the posterior predictive distribution, you can obtain interval estimates etc.

So answering your question: having the MCMC samples from the posterior distribution of the parameters enables us to calculate all the statistics of interest about the posterior distribution, to visualize it (plot histograms, or kernel density estimates from the samples), and make predictions.