# Bayesian update for a univariate normal distribution with unknown mean and variance

Suppose I have some random process $X$ which is emitting values which follow a normal distribution:

$$X \sim N(μ, σ^2)$$

Both $μ$ and $σ$ are unknown, so I want to model each of them with their own distribution which I will update every time I observe a new value.

How can I do this?

For $μ$ it seems obvious that I should model it with its own normal distribution: $μ \sim N(μ_μ, σ_μ^2)$. For $σ^2$ it's not clear what distribution I should use - my googling so far suggests that inverse-gamma would make the math work-out nicely but it's not clear to me that it even makes sense to use two independent distributions for $μ$ and $σ^2$.

So my question is: what mathematical model should someone use in this situation (or, if there's a choice, what are the options), and how exactly does one calculate the posterior parameters of the model given the prior parameters and an observation $x$?

• You should tell as if you want a closed form for the posterior or if you want to use some MCMC algorithm. Sep 3 '18 at 14:16
• Closed form, I think. What would be the advantage of using some MCMC algorithm?
– Shum
Sep 3 '18 at 14:29
• If you are willing to use MCMC algorithms, than you can chose any prior distribution for $\mu$ and $\sigma2$. While if you want a joint closed posterior closed form you must use an "inversegamma-normal" distribution for $(\mu,\sigma^2)$ (a special case of the normal inverse wishart distribution). For details have a look at cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf Sep 3 '18 at 14:36
• Thanks! I followed that paper and tried implementing a normal-gamma distribution. I initialized it with arbitrary parameter values then trained it using a thousand random samples from N(μ = 2, σ^2 = 9). The result was that the mean converged on the correct value (ie. 2), but the variance converged on zero (?). Am I doing something wrong, or have I misunderstood what this is for? I expected that my normal-gamma distribution would converge on N(2, 9).
– Shum
Sep 3 '18 at 15:09
• It also worries me that the parameters of the normal-gamma distribution don't converge. κ, α, and β just keep getting bigger and bigger with each update and I'm worried that after a billion updates I'll starting running into floating-point weirdness. Is there some technique people use to avoid this?
– Shum
Sep 3 '18 at 15:12

If $$X_i \stackrel{ind}{\sim} N(\mu,\sigma^2)$$ where $E[X_i] = \mu$ and $Var[X_i] = \sigma^2$, then the fully conjugate prior for an unknown mean $\mu$ and variance $\sigma^2$ is $$\mu|\sigma^2 \sim N(m,\sigma^2/k) \qquad \sigma^2 \sim \mbox{Inv-}\chi^2(v,s^2)$$ where $\mbox{Inv-}\chi^2(v,s^2)$ is the scaled inverse-chi-squared distribution with mean $vs^2/(v-2)$ for $v>2$ and variance $2v^2s^4/[(v-2)^2(v-4)$ for $v>4$ which is equivalent to $IG(v/2,vs^2/2)$, an inverse gamma distribution.

The posterior under this model and prior is $$\mu|\sigma^2,x_1,\ldots,x_n \sim N(m',\sigma^2/k') \qquad \sigma^2|x_1,\ldots,x_n \sim \mbox{Inv-}\chi^2(v',(s')^2)$$ with $$\begin{array}{rl} k' &= k+n \\ m' &= [km+n\overline{x}]/k' \\ v' &= v+n \\ v'(s')^2 &= vs^2 + (n-1)S^2 + kn(\overline{x}-m)^2/k' \end{array}$$ where $\overline{x} = \frac{1}{n}\sum_{i=1}^n x_i$ is the sample mean and $S^2 = \frac{1}{n-1}\sum_{i=1}^n (x_i-\overline{x})^2$ is the sample variance.

• Thanks! I tried this but I get the exact same problem as I got before: the variance converges on zero rather than the actual variance of the normal distribution I'm trying to learn. Does that make sense? I posted my code here if it helps explain what I'm doing: play.rust-lang.org/….
– Shum
Sep 4 '18 at 4:16

A bit late, but here is an implementation based on the formula stated on wikipedia (https://en.wikipedia.org/wiki/Conjugate_prior#When_likelihood_function_is_a_continuous_distribution):

import random

u0 = 0.0
v = 0
b = 1.

for i in range(10000):
x = random.gauss(2., 1.5)

b = b + 0.5 * v / (v + 1) * (x - u0)**2
u0 = (v*u0 + x) / (v + 1)
v += 1

print("mean:", u0, "sd:", (2*b*(v+1) / (v*v))**0.5)


Make sure to first update b, then u0 and then v.