# Sequential Bayesian updating of mean and variance of normal distribution

I am trying to write some code to learn the parameters of a normal distribution. I am new to this, and I have patched together the equations from various sources, which may be part of the problem. In particular, I seem to manage to update the prior with a batch of data, though I am having some difficulties interpreting the constant involved. However, I struggle to move from that to sequential updating, probably due to the same probelm with the constant.

Take a data vector $$x = (x_i, ... , x_n)$$, with $$x_i \sim \mathcal{N}(\mu,\lambda^{-1})$$, assuming exchangeable $$x_i$$, and where $$\lambda \triangleq \frac{1}{\sigma^2}$$ is the precision. I want to learn $$p(\mu,\lambda | x) = p(x | \mu , \lambda)p(\mu,\lambda)$$. Exploiting that $$p(\mu,\lambda) = p(\mu | \lambda) p(\lambda)$$, I can use the following normal-gamma prior: $$p(\mu,\lambda) = \mathcal{N}(\mu|\mu_0,(n_0\lambda)^{-1}) Gam(\lambda|\alpha_0,\beta_0),$$ where $$n_0$$ should be some sort of scaling constant.

I should then be able to use the following equations to learn the model parameters (see here): $$\alpha = \alpha_0 + \frac{n}{2}\\ \beta = \beta_0 + \frac{1}{2}\sum_{i=1}^{n}(x_i - \mu)^2 = \beta_0 + \frac{1}{2}\sum_{i=1}^{n}(x_i - \overline{x})^2 + \frac{1}{2}\frac{n n_0}{n + n_0}(\overline{x} - \mu_0)^2) \\ \lambda = \frac{\alpha}{\beta} \\ E[\mu|x,\lambda] =\frac{n\lambda}{n\lambda + n_0\lambda} \overline{x} + \frac{n_0\lambda}{n\lambda + n_0\lambda}\mu_0 = \mu_0 + \frac{n\lambda}{n\lambda + n_0\lambda} (\overline{x} - \mu_0),$$ where $$\overline{x} = \frac{1}{n}\sum_i x_i$$ is the sample mean. This seems to work on simulated data, but the results are very sensitive to my choice of $$n_0$$. The best results obtain for $$n_0 = \frac{\alpha_0}{\beta_0} var(x)$$. Hence my first question:

Q1: Is $$n_0 = \frac{\alpha_0}{\beta_0} var(x)$$ the correct way of fixing this constant'? Can it be referred to as a constant at all, since it is clearly dependent on the data at hand?

Next, I would like to move on to sequential sampling, i.e. updating my prior with one observation at a time. The equation for $$\alpha$$ obtains trivially by setting $$n=1$$. For $$\beta$$, the middle term should drop out, since $$x_i = \overline{x}$$ for a single observation. My trouble is, once more, with $$n_0$$. For a single observation, the sample variance is not defined. This suggests that my solution above may have been wrong, or at least not generally applicable. I am also very confused of whether this constant' should now be given one fixed value, or whether it should rather be adapted observation by observation.

Q2: how can I update the different parameters observation by observation? In particular, I am unclear what to do with the parameter $$n_0$$ if I set $$n=1$$.

The $$n_0$$ parameter is a parameter of a prior distribution. It is not estimated but assumed a priori. It serves as a "pseudocount" of the number of observations seen a priori and translates to "how strongly we believe" in the prior (Murphy, 2007). Notice this part of the equation

$$\frac{n\lambda}{n\lambda + n_0\lambda} \overline{x} + \frac{n_0\lambda}{n\lambda + n_0\lambda}\mu_0$$

where we take the weighted average between the observed $$\overline{x}$$ and the prior $$\mu_0$$. The weights depend on the number of samples you observed $$n$$ and $$n_0$$ "imaginary" samples from the prior. The more "imaginary" samples, the more weight goes to the prior.

If you want to update sequentially $$n=1$$ for each step. To choose the $$n_0$$ you need to consider how many samples would you consider to be enough to overwhelm the prior. Keep in mind that this is just intuition on how to think about it, because how strong the prior is depends on how much information it contains as compared to the data, so it is not only about the sample size.

Check "Conjugate Bayesian analysis of the Gaussian distribution" by Kevin P. Murphy for a more detailed description.

• OK, thank you, that makes a lot of sense. It seems like I had made an additional mistake in my previous simulations. Now changing $n_0$ yields the expected results. However, in terms of sequential updating I am still missing something: putting $n=1$ in the equation for the mean above gives a constant weight to any new observation, and there thus seems to be excess learning. It seems one would need to augment $n_0$ at every iteration, to reflect that this is now based on a richer observation set. Is that correct, or is there some other way to adjust the equation? Aug 10, 2022 at 15:40
• @MaxMontana no, why? What is excess about each observation having the same weight? It's equivalent to using all the data at once, in such a case why would you weigh different observations differently? If you learned something about the data, it should make you more certain about it rather than re-weighting the prior to “discard” what you've learned.
– Tim
Aug 10, 2022 at 18:03
• Then something else must be wrong with my equations. I get over-learning in the sense that I get recency effects in sequential updating. This directly violates exchangeability, which I know to hold true in my simulated data. Using the same data as for the batch updating, I thus get markedly different results. This happens with alpha[i] <- alpha[i-1] + 1/2, beta[i] <- beta[i-1] + (n0)/(2*(1 + n0)) * (x[i] - mu0[i-1])^2, lambda[i] <- alpha[i]/beta[i], and mu0[i] <- mu0[i-1] + (lambda[i])/(lambda[i] + n0 * lambda[i] ) * (x[i] - mu0[i-1]). Any idea where I am going wrong? Aug 10, 2022 at 19:01
• @MaxMontana you need to have a bug in your code: it should give exactly the same result if you run it with all the observations vs if you update sequentially. Example of code that works: gist.github.com/twolodzko/5704deaa0617df13895557df59c165f7
– Tim
Aug 11, 2022 at 8:17

Thank you Tim for all the useful input and references!

Augmenting $$n_0$$, as I had suspected, is indeed what did the trick. With each observation, one additional "pseudocount" is thus added to the prior, making our belief in the prior increasingly strong. The notes of Kevin Murphy referenced by Tim above mention this case at the bottom of page 9. The equations thus become: $$\alpha_i = \alpha_{i-1} + \frac{1}{2}\\ n_i = n_{i-1} + 1\\ \beta_i = \beta_{i-1} + \frac{1}{2}\frac{n_i}{n_i+1}(x_i - \mu_{i-1})^2 \\ \mu_i = \mu_{i-1} + \frac{1}{n_i + 1}(x_i - \mu_{i-1}) .$$

• The intuition of "making it increasingly strong" is incorrect. Notice that when updating the prior using all the data at once you update $n_0 = n_0 + n$, while sequentially you update $n_i = n_{i-1} +1$, so the update is the same in both cases. At each step, you update a prior and get a posterior that serves as a new prior, but the "original" prior is still $n_0$. You are not changing your prior beliefs, but learning from the data.
– Tim
Aug 11, 2022 at 10:33
• Hi Tim, that is fair enough. The issue was that the part $n_p = n_0 + n$ has been missing from my initial equation system, which is what caused the problem in sequential updating (in batch updating this problem did not show up because it was a one-time process). "Learning from the data" is indeed a better description of the process than "making the belief in the prior increasingly strong", but it remains that the key was to increase $n_0$ at every iteration to obtain the posterior, which then serves of the new prior. Aug 11, 2022 at 11:33