Should updating one data point at a time or all change the posterior of a normal-inverse-gamma?

Question

I have implemented the normal inverse gamma distribution per section 3 of https://people.eecs.berkeley.edu/~jordan/courses/260-spring10/lectures/lecture5.pdf in some code. However, I've noticed something that I didn't expect. The posterior differs depending on how I choose to update.

So, imagine I have data with n=100. If I update the posterior by taking my prior, updating with n=10, using that posterior as my prior, and then continuing that process until I've updated with all of the 100 data points, my variance (i.e. inverse-gamma portion of the distribution) is different than if I just do one update with the entire set of 100 points once.

Am I wrong in thinking that the posterior should be the same independent of how I choose to update? Am I missing something in how the posterior is computed?

Edit: adding the posterior function for the inverse-gamma:

$\sigma^{2}|x \sim Inv-Ga(\alpha + \frac{n}{2}, \beta + \frac{1}{2}\sum(x_{i}-\bar{x})^2 + \frac{nn_{0}}{2(n + n_{0})}(\bar{x} - \mu_{0})^2)$

Parameter updating: $NIG(\frac{\nu\mu_{0}+n\bar{x}}{\nu+n}, \nu+n,\alpha + \frac{n}{2}, \beta + \frac{1}{2}\sum(x_{i}-\bar{x})^2 + \frac{n\nu}{2(n + \nu)}(\bar{x} - \mu_{0})^2)$

Here's our Scala implementation in case anyone is interested: https://github.com/udemy/statistics/blob/master/src/main/scala/com/udemy/statistics/distribution/NormalInverseGamma.scala

Answer 1

I suspect this might be a programming error; I've implemented the updates in R and gotten the expected result. Perhaps the code below, which I've tried to write for clarity rather than efficiency, when compared with the original Scala code, will help reveal the problem:

mu_update <- function(mu, v, x) (mu*v + sum(x)) / (v+length(x))

v_update <- function(v, n) v+n

alpha_update <- function(alpha, n) alpha+n/2

beta_update <- function(beta, v, mu, x) {
  n <- length(x)
  beta + (sum((x-mean(x))^2))/2 + (n*v/(2*(n+v)))*(mu-mean(x))^2  
}

update <- function(parms, x) {
  n <- length(x)
  list(alpha = alpha_update(parms$alpha, n),
       beta = beta_update(parms$beta, parms$v, parms$mu, x),
       mu = mu_update(parms$mu, parms$v, x),
       v = v_update(parms$v, n))
}

x1 <- rnorm(10)
x2 <- rnorm(10)
parms <- list(alpha=2, beta=2, mu=0, v=2)
parms_2step <- update(update(parms, x1), x2)
parms_1step <- update(parms, c(x1, x2))

Executing it and finding the difference of the two parameter lists yields:

> unlist(parms_1step) - unlist(parms_2step)
alpha  beta    mu     v 
    0     0     0     0 

indicating that the two updates generated the same final parameter estimates, as expected / hoped for.

A more extensive test involving 10 updates gives the same result:

parms <- list(alpha=2, beta=2, mu=0, v=2)
parms_2step <- parms
x <- list()
for (i in 1:10) {
   x[[i]] <- rnorm(10)
   parms_2step <- update(parms_2step, x[[i]])
}
parms_1step <- update(parms, unlist(x))

delta <- unlist(parms_2step) - unlist(parms_1step)

gives the following results:

> delta
       alpha         beta           mu            v 
0.000000e+00 7.105427e-15 0.000000e+00 0.000000e+00 

which is evidently due to roundoff error.

Also note that Scala, when dividing two integers, rounds down, so you have to make sure that expressions like n*v/(2*(n+v)) when executed will actually return the correct, floating point, number.

Answer 2

This is already pointed out by J. Bowman in his comments, but let me insist upon the fact that you cannot use this posterior today = prior tomorrow when looking solely at the parameter $\sigma^2$ because the integration over the other parameter $\mu$ turns the points of the sample into dependent random variables, which means that the actualisation principle cannot apply as such. More precisely, the distribution of $x_i$ given $x_1,\ldots,x_{i-1}$ and $\sigma^2$ is no longer a Normal independent from $x_1,\ldots,x_{i-1}$.