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$?
 A: 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.
A: 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.
