3
$\begingroup$

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.

$\endgroup$
3

0

Your Answer

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