# Recreating data variance from the posterior distribution

## Recreating data variance from the posterior distribution

Take a set of data points $$(x, y)$$ with (Gaussian) uncertainties $$\sigma_y$$ on the $$y$$ coordinate; they are modeled as $$y \sim f(x; \alpha) + \epsilon$$ with a function $$f$$ depending on some parameters $$\alpha$$ and including Gaussian noise $$\epsilon$$ with variance $$\sigma_y^2$$.

This results in some posterior distribution for the parameters $$\alpha \sim p(\alpha | y)$$. Suppose then we are given $$p(\alpha | y)$$ and the values $$x$$, but not $$y$$ nor $$\sigma_y$$. We can draw a "fit cloud" by sampling from $$p(\alpha | y)$$ and drawing a version of $$f(x; \alpha)$$ for each sample:

In this example $$f(x; \alpha) = a x + b$$ is a straight line (with parameters $$\alpha = (a, b)$$) and the errors are constant, $$\sigma_y \equiv 1$$. The assumption here is that we know about the black cloud, while we do not have access to the red points.

The question then is: given $$x$$ and $$p(\alpha | y)$$, can we make a reasonable guess $$(\bar{y}, \bar{\sigma}_y)$$ for the values $$y$$ and their variances $$\sigma_y^2$$ such that the new set of points, $$(x, \bar{y})$$ with $$\bar{\sigma}_y$$, can be fit with the same model as before resulting in the same posterior $$p(\alpha | y)$$?

### Constant model

The simplest example would be with a constant model $$f(x; c) = c$$: then, our fit would amount to finding the average of the $$y$$s, and the "fit cloud" would have an expected standard deviation of $$\sigma_y / \sqrt{N}$$, where $$N$$ is the number of data points.

So in this context, we could estimate the errors $$\sigma_y$$ as

$$\bar{\sigma}_y \approx \text{std}_{p(\alpha|y)}(f(x, \alpha)) \sqrt{N} = \text{std}(c_i) \sqrt{N},$$

and the $$y$$s as

$$\bar{y} \approx \text{mean}_{p(\alpha|y)}(f(x, \alpha) = \text{mean}(c_i),$$

where $$c_i$$ are sample points from the posterior distribution of the parameter $$c$$, while with notation such as $$\text{mean}_{p(\alpha|y)}f$$ I mean an average of the values $$f$$ takes on while we vary $$\alpha$$ according to the distribution $$p(\alpha|y)$$.

### The two-parameter case

Multiplying by $$\sqrt{N}$$ does not seem to generalize, though. In the aforementioned linear fit scenario, heuristically I have found that I need to multiply by roughly $$\sqrt{N / 2}$$ to reproduce the fit results. This is what that looks like:

To clarify, here I'm taking $$\bar{y}$$ to be the mean of $$f(x; \alpha)$$ for a given $$x$$ and varying $$\alpha$$ according to $$p(\alpha | y)$$, and $$\bar{\sigma}_y$$ to be its standard deviation. The results roughly match: this is what $$p(\alpha | y)$$ looks like

#### Explicit question

Is this a known problem in statistics? Is the scaling factor an approximation of some exact result? Is it true that, in general, the scaling is $$\sqrt{N/d}$$ with $$d$$ being the number of the fit parameters?

### Robustness checks

I have verified that the scaling relation $$\sqrt{N/d}$$ is robust to changes in $$N$$, changes in the range of the $$x$$ points (e.g. making them span [0, 30] instead of [0, 3], as well as the $$x$$ points for the reconstruction not being the same as the ones originally used for the fit.

### Increasing the dimensionality

I have tried increasing the model dimensionality to 3 ($$y \sim ax^2/2 + bx + c$$), and things seem to hold up if I multiply the errors by $$\sqrt{N/3}$$:

I also tried with an even more complex, nonpolynomial $$d=5$$ fit: $$y \sim A \cos(2 \pi f x - \phi) + bx + c$$

The prescription $$\sqrt{N/d}$$ still works! Here are the posteriors obtained by fitting the reconstructed data points with errors multiplied by $$\sqrt{N/5}$$.

### Related terminology

Reading up on the issue I found that (probably) the thing I'm looking for is called a "posterior predictive distribution", and it would be written as

$$p( \bar{y} | y) = \int p(\bar{y}|\alpha) p(\alpha | y) \text{d}\alpha$$

The issue with this way of expressing it seems to be that I would need to know $$p(\bar{y}|\alpha)$$, a Gaussian with mean $$f(x; \alpha)$$ and variance $$\sigma_y^2$$ - but I don't know the latter!

### The effect on Bayesian evidence

I'm performing these fits with the nested sampling package Ultranest, therefore I also have information on the evidence for these fits. In all cases, the evidence $$\log Z$$ for the model $$f$$ I get with $$(x, y, \sigma_y)$$ is much lower than the one I get with $$(x, \bar{y}, \bar{\sigma}_y)$$. For example, in the last case (with the cosine) the evidence is $$\log Z = -104.6 \pm 0.5$$ with the original data, $$\log \bar{Z} = -20.9 \pm 0.3$$ with the reconstructed data. I figure this is expected: the reconstructed data is by construction a much better fit to the model.

I can artificially make "dirtier" reconstructed data by computing the standard deviation as the $$\bar{\sigma}_y = \text{std}_{p(\alpha | y)}(f(x, \alpha))$$ and then adding a normally distributed offset to each data point: $$\bar{y} = \text{mean}_{p(\alpha | y)}(f(x, \alpha)) + \mathcal{N}(0, \bar{\sigma}_y)$$.

This results in a much closer value for the evidence, though still not compatible within the error bounds: $$\log \bar{Z} = -100.3 \pm 0.3$$. Simultaneously, though, the posteriors stop being so perfectly compatible: their variances remain correct, but a bias is introduced (as expected: it was put there by hand!)

• Related: stats.stackexchange.com/questions/145724/… but there is no proper answer there! Feb 7 at 14:46
• Tongue in cheek - don't Bayesians disavow variance in the data? Feb 7 at 17:35
• It would be much easier without it for sure! 😁 Feb 7 at 22:47