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:

Fit cloud

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:

Reconstructed vs true points.

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

Posterior distribution from the true and reconstructed data points.

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}$:

Reconstructed points - quadratic fit

Posterior distributions - quadratic fit

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

Fancy 5D fit

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}$.

The posteriors match in 5D

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!)

Posterior matching with biased points.



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.