What is the correct posterior when data are sufficient statistics?

Question

Say you have N observations that are iid.

$$ \forall i, \quad p(X_i=x_i|\mu,\sigma,I) = \frac{1}{\sqrt{2\pi}\sigma} \exp\left(-\frac{1}{2\sigma^2}(x_i-\mu)^2\right)$$ then $$ p(x_1,\dots,x_N|\mu,\sigma,I) = \frac{1}{(\sqrt{2\pi}\sigma)^N} \exp\left(-\frac{N}{2\sigma^2}[(\bar{x}-\mu)^2 + \bar{\sigma}^2]\right)$$ with $\bar{x}$ the sample mean and $\bar{\sigma}^2$ the sample variance. Say instead of the observations, you are given their sufficient statistics: e.g. sample mean, sample variance and number of points. You want to estimate the true mean, variance and number of points. What is the correct posterior (or likelihood) for this case, assuming that the data follow a gaussian distribution? So rigorously speaking, $$p(\bar{x},\bar{\sigma}^2,N|I)=p(x_1,\dots,x_N|\mu,\sigma,I) J(\bar{x},\bar{\sigma}^2,N)$$ where $J$ is the jacobian of the transformation from the $x_i$ to the sufficent statistics. But I can drop the Jacobian when I write my likelihood function of the sufficient statistics, because the Jacobian does not depend on the parameters to be estimated. Right?

Now assume you are only given sample mean and standard deviation. Is there any way to estimate the number of points involved? I feel like not but I have difficulties expressing why.

Answer 1

When N is known, you can use the fact that the sample mean and variance are independent (conditionally on $\mu, \sigma$) and have known distributions. You can justify this approach because you already known that the sample mean and variance are sufficient statistics and the likelihood has to factor this way (your comment about the Jacobian is correct as well). Or like @jpillow says, you could also derive the result by starting with the full likelihood and charging through the algebra. You'll end up at the same place

But the sample mean and variance aren't sufficient for $\mu, \sigma$ if $N$ is unknown. It isn't be hard to show mathematically. Start with the likelihood you derived earlier, replacing $N$ with $N(\boldsymbol{x})$ to make the dependence explicit :

$$ p(x_1,\dots,x_{N(\boldsymbol{x})}|N=N(\boldsymbol{x}),\mu,\sigma,I) = \frac{1}{(\sqrt{2\pi}\sigma)^{N(\boldsymbol{x})}} \exp\left(-\frac{N(\boldsymbol{x})}{2\sigma^2}[(\bar{x}-\mu)^2 + \bar{\sigma}^2]\right)$$

Now if $\bar x, \bar\sigma$ were sufficient, the above would have to factor as $h(\boldsymbol{x})g_{\mu, \sigma}(\bar \mu, \bar\sigma)$ where $h$ doesn't involve $\mu, \sigma$ and $g$ only depends on the data through the sufficient statistics. Because of the leading term $\frac{1}{(\sqrt{2\pi}\sigma)^{N(\boldsymbol{x})}}$ you can't get there.

I suppose the intuition is that if I give you a sample mean and variance but not the sample size you can't be sure if it came from 5 observations or 5 million...

Answer 2

If you don't know the value of $N$, then we simply integrate it out, the same as any other nuisance parameter. We know that $N$ is discrete, and as long as your standard deviation is not zero, then you know that you have at least $2$ observations, so $N\geq 2$. So we write out the posterior as:

$$p(\mu,\sigma|\overline{x},s^2,I)=\sum_{N=2}^{\infty}p(\mu,\sigma,N|\overline{x},s^2,I)$$ $$=\sum_{N=2}^{\infty}p(N|\overline{x},s^2,I)p(\mu,\sigma|N,\overline{x},s^2,I)=\frac{\sum_{N=2}^{\infty}p(N|\overline{x},s^2,I)p(\mu,\sigma|N,\overline{x},s^2,I)}{\sum_{N=2}^{\infty}p(N|\overline{x},s^2,I)}$$

The last form is written to show that we basically have a weighted average of $p(\mu,\sigma|N,\overline{x},s^2,I)$ for different $N$. This term is the standard posterior distribution you get in most cases, so I won't discuss it here. I will focus on the first term, which is the weights. Now we can re-arrange it using Bayes theorem to get:

$$p(N|\overline{x},s^2,I)\propto p(N|I)p(\overline{x},s^2|N,I)$$

From @JMS comment, he seems to think that $p(\overline{x},s^2|N,I)$ does not depend on $N$, so that we are basically weighting by the prior $p(N|I)$. This is incorrect, unless we have very little prior information about $\mu$ and $\sigma$. I will show this with a specific example, followed by my attempt at explaining why this is the case. Now we can further decompose this as:

$$p(\overline{x},s^2|N,I)=\iint p(\overline{x},s^2,\mu,\sigma|N,I)\,d\mu\, d\sigma$$ $$=\iint p(\overline{x},s^2|\mu,\sigma,N,I)p(\mu,\sigma|N,I)\,d\mu\,d\sigma$$ $$=\iint p(\overline{x}|\mu,\sigma,N,I)p(s^2|\mu,\sigma,N,I)p(\mu,\sigma|N,I)\,d\mu\,d\sigma$$ The last line follows from the know fact that for normal with known mean and variance, sample mean is independent of sample standard deviation. We also have:

$$p(\overline{x}|\mu,\sigma,N,I)=\sqrt{\frac{N}{2\pi\sigma^2}}\exp\left(-\frac{N}{2\sigma^{2}}[\overline{x}-\mu]^2\right)$$ $$\frac{Ns^2}{\sigma^2}\sim \chi^2_{N-1}\implies p(s^2|\mu,\sigma,N,I)=\frac{1}{\Gamma[\frac{N-1}{2}]}\left(\frac{Ns^2}{2}\right)^{(N-1)/2}\sigma^{-(N-1)}\exp\left(-\frac{Ns^2}{2\sigma^2}\right)$$

Now if I was to specialise to the prior $p(\mu,\sigma|N,I)\propto \sigma^{-\alpha-1}\exp\left(-\frac{\beta}{\sigma^2}\right)$ we get:

$$p(\overline{x},s^2|N,I)\propto\frac{\Gamma[\frac{N-1}{2}+\alpha]}{\Gamma[\frac{N-1}{2}]}\left(\frac{Ns^2}{2}\right)^{(N-1)/2}\left(\frac{Ns^2}{2}+\beta\right)^{-(N-1+\alpha)/2}$$

Where factors independent of $N$ have been absorbed into the "$\propto$" symbol. Note that this still depends on $N$ unless we use the Jeffrey's prior $\alpha=\beta=0$. Hence, if you know something about the parameters of the normal distribution, you can use this information to estimate $N$ from $s^2$ and $\overline{x}$. Now where is the intuition behind this? I don't really know, but I'll have a go, as I was initially quite surprised myself at this result, but on a bit of reflection it does make sense.

We know from the chi-square pivotal quantity, that we can roughly expect $Ns^2\approx (N-1)\sigma^2$. Now if we use an informative prior for $\sigma^2$ like I did then we have $\sigma^2\approx \frac{\alpha}{2\beta}$ (this is the prior mean for $\sigma^2$). But we can then use this to solve for $N$, and we get $$N\approx\left(1-\frac{2s^2\beta}{\alpha}\right)^{-1}=\left(1-\frac{s^2}{\hat{\sigma}^{2}_{PRIOR}}\right)^{-1}$$

Now why doesn't this work for the Jeffreys prior? Well, when we use the Jeffreys prior, the prior estimate is undefined so the above equation is arbitrary for Jeffreys prior. Additionally, we essentially estimate $\sigma^2\approx s^2$ with the Jeffreys prior and so we already have $Ns^2\approx (N-1)\sigma^2$ independently of the value of $N$ - hence we cannot use this to help us estimate $N$. You can see this also in the form of $p(\overline{x},s^2|N,I)$ as it essentially has the heuristic form:

$$\frac{\text{normalisation with }\mu,\sigma\text{ estimated from N obs and prior information}}{\text{normalisation with }\mu,\sigma\text{ estimated from N obs}}$$

So if you are in some sort of sequential situation, you can get estimates for $N$ after the receiving the first mean and variance.

Answer 3

You can get there more easily if you just collect terms:

$P(\{X_i\}|\mu,\sigma^2) = \prod_{i=1}^N \frac{1}{\sqrt{2\pi \sigma^2}} \exp (-\frac{(x_i-\mu)^2}{2\sigma^2}) = \frac{1}{\sqrt{|2\pi \sigma^2 I|}}\exp (-\frac{\sum (x_i-\mu)^2}{2\sigma^2})$

If you massage the terms in that exponent (start by multiplying out the quadratics), you should be able to get to the desired expression. (This is probably easier than computing the Jacobian explicitly). Actually, it's probably easier if you plug in explicit expressions for sufficient statistics $\bar \mu$ and $\bar \sigma$ and work backward to the expression here.