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$?

  • $\begingroup$ You should tell as if you want a closed form for the posterior or if you want to use some MCMC algorithm. $\endgroup$
    – niandra82
    Sep 3 '18 at 14:16
  • $\begingroup$ Closed form, I think. What would be the advantage of using some MCMC algorithm? $\endgroup$
    – Shum
    Sep 3 '18 at 14:29
  • $\begingroup$ 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 $\endgroup$
    – niandra82
    Sep 3 '18 at 14:36
  • $\begingroup$ 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). $\endgroup$
    – Shum
    Sep 3 '18 at 15:09
  • $\begingroup$ 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? $\endgroup$
    – 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.

  • $\begingroup$ 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/…. $\endgroup$
    – 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.