# How to score predictions in test set taking into account the full predictive posterior distribution?

I have three predictive models (regressions) which parameters are estimated by Markov Chain Monte Carlo. Predictions are made over a test set of size $N$. Since I compare the models under different settings, I would like to summarize the results in a table.

I have considered using Means Squared Error (MSE) or $R^2$ but I wonder if I could be fully bayesian and use all the posterior sampled by MCMC, and not just the means of the predictions.

For instance, MSE would be something like: $$\text{MSE} = \frac{1}{N} \sum_{i=1}^N (y_i^{true} - \mathbb{E}[y_i^{MC}])^2$$

where $\mathbb{E}[y_i^{MC}]$ is the average of the $y_i$ samples.

However, it seems reasonable to me that a tight posterior around the good prediction should have a better score than a flatter posterior. And yet, I don't find any reference to this kind of method.

Questions:

• What is the/a correct bayesian way to do it? Does it have a name?

• What is the name for techniques that check the accuracy of posteriors-based predictions in test sets where we know the true outputs? 

 : I say that because I've realized that a lot of bayesian literature talks about assessing the quality of a prediction (e.g.: the outcome of the elections) before we know the real outcome.

Edit:

What about taking the mean of the mean squared errors? That is averaging the squared errors for the samples at every test point, and then averaging the results through all the points.

\begin{align} \frac{1}{N}\sum_{i=1}^N\left[\frac{1}{M}\sum_{j=1}^M\left(y_i^{true}- y_i^{(j)}\right)^2\right] \end{align} where $M$ is the number of samples from the MCMC and $y_i^{(j)}$ is the prediction, for the input $x_i$, using the $j$-th sample of the model parameters $\boldsymbol{\theta^{(j)}}$.

• This looks like you should be use the posterior variance - this is what you would be using in the second formula. Beware of using highly auto-correlated MCMC chains to do this though, as it will make the variance seem smaller. – probabilityislogic Feb 28 '15 at 22:27
• Sorry! There was a mistake in the third formula. I meant the squared difference between the posterior prediction and the true value. Do you mean that I should report both the MSE and the posterior variance? – alberto Feb 28 '15 at 22:35
• The average posterior variance for each data point is the posterior expectation of what you have called MSE. – probabilityislogic Mar 1 '15 at 0:08
• I just found this survey projecteuclid.org/download/pdfview_1/euclid.ssu/1356628931 and discovered the "L-criterion" (page 198, eq 111). Isn't it closer (or exactly) to what I am doing in my edit? – alberto Mar 1 '15 at 9:10

You want the posterior predictive distribution, i.e. for a test vector $\tilde{y}$ and training vector $y$

$$p(\tilde{y}|y) = \int p(\tilde{y}|\theta) p(\theta|y) d\theta).$$

As you point out, if this distribution is tighter around the test data $\tilde{y}$, then it will give a higher density.

If you only have samples from the posterior, then you can approximate this via

$$p(\tilde{y}|y) \approx \frac{1}{M} \sum_{j=1}^M p(\tilde{y}|\theta^{(j)})$$

where $\theta^{(j)} \sim p(\theta|y)$.

• I was re-reading Gelman's Bayesian Data Analysis and I think it's exactly that ! So, to be clear $\bar{y}$ would be the observed / true / target predictions, wouldn't they? – alberto Mar 2 '15 at 21:36
• For the record, I see that if we factorized $\tilde{y}$ into the different $\tilde{y_i}$ and take the logarithmic version of that is called "log pointwise predictive density". – alberto Mar 3 '15 at 8:17
• Yes $\tilde{y}$ is the true values. Depending on the model structure, you may be able to assume the $\tilde{y}_i$ are conditionally independent and thus factorize. – jaradniemi Mar 3 '15 at 16:09