# What is the correct posterior when data are sufficient statistics?

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.

• $N$ is not know? Could you please confirm? – suncoolsu Jul 15 '11 at 12:05
• not in the last question. – yannick Jul 15 '11 at 15:56

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

• great! almost there! now how do you write the likelihood if you're given the sample mean and variance? Cause you have to compute the Jacobian this time, no? Or you simply assume it's a gaussian? – yannick Jul 15 '11 at 16:08
• @yannick I'm not sure what you mean; are you looking for the likelihood in the case of unknown $N$? When you put down the likelihood it's conditional on any unknown parameters, including $N$, so the above expression should work. But you have more fundamental problems in that the data (that is, the sample mean and variance) carry no information about the sample size. You need more information about $N$ coming in through either your prior on $N$, or from some other data source. – JMS Jul 15 '11 at 16:32
• precisely, if the data don't give you N, you have no hope of getting any reliable estimate if your prior is uninformative. But you still want an estimate of the true mean and standard deviation. So how should you write the likelihood to estimate these? Just a one-dimensional gaussian, like I've seen a lot of times? But then you are making an assumption on the value of N(x)! – yannick Jul 15 '11 at 17:37
• @yannick I'm still not sure I understand your point. You can specify the likelihood $p(\bar\mu, \bar\sigma^2 | \mu, \sigma, N)$ based on the sampling distribution of $\bar\mu, \bar\sigma^2$ since you've assumed that $x_i\sim N(\mu, \sigma^2)$. What you want is the marginal posterior for $\mu, \sigma^2$. Taking $D$ as the data and slightly abusing notation, this would involve computing $\sum_{i=1}^{n=\infty} p(\mu, \sigma^2, N=i|D) = \sum_{i=1}^{n=\infty} p(\mu, \sigma^2| N=i, D)p(N=i|D)$. – JMS Jul 15 '11 at 19:21
• But the data have no information about $N$ so $p(N=i|D) = p(N=i)$ and you're just integrating over your prior. So yes, you need an informative prior on $N$. But any posterior inferences are going to be so heavily driven by your prior on $N$ that the results are basically useless unless you can (very strongly) justify your choice of prior. – JMS Jul 15 '11 at 19:24

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.

• that's a neat contribution! so one can estimate N if one uses a normal-inverse gamma prior for the mean and variance parameters. However, I feel that your estimate is strongly affected by the choice of the prior, and in that sense is not very robust. Basically, your data sets $(N-1)/N\sigma$ and your prior sets $\sigma$. So you can make inferences on $N$, but they don't come from the data, they're only due to the prior. – yannick Aug 20 '11 at 16:32
• @yannick - yes the prior information matters - but this hardly calls for a "non-robust" label - this gives the impression that it must be wrong. If the prior information does matter, then it should only be called "non-robust" if that prior information is incorrectly specified. If you have prior information, why would you want your inference to be independent of that information? That would be just as crazy as ignoring the data. And remember, my explanation was an attempt at intuitively understanding why the posterior was not equal to the prior. – probabilityislogic Aug 21 '11 at 21:03
• ...(cont'd)... there may be other features in the posterior which make more accurate predictions. – probabilityislogic Aug 21 '11 at 21:04

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.

• i think OP's problem is more complicated than usual because $N$ is also unknown. I am not sure if the question makes sense. $N$ must be known, right? – suncoolsu Jul 15 '11 at 12:05
• hmm not sure because the prefactor does not have the right exponent. What's "I" by the way? And still, if you change from the $X_i$ to the sufficient statistics, you have to compute a Jacobian. – yannick Jul 15 '11 at 16:10
• Sorry, that was just laziness--$I$ is the $N\times N$ identity matrix; this is another way of writing multivariate normal. It's the same as writing $2 \pi \sigma^2$ to the power $N/2$. – jpillow Jul 15 '11 at 17:01
• And no--no need to compute the Jacobian here (unless you want to for some reason). The likelihood here is just the product of the likelihoods from each observation. With a little effort you can collect the terms in the product of likelihoods and rewrite them in terms of the sufficient statistics. (Note, as JMS says, you DO need to know $N$). – jpillow Jul 15 '11 at 17:13