2
$\begingroup$

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

|cite|improve this question
$\endgroup$
  • 2
    $\begingroup$ You are not wrong in thinking that the posteriors should be the same; "yesterday's posterior is today's prior" is your guide. If you want further help, posting the math you used as the basis for your implementation would be the next step. $\endgroup$ – jbowman Sep 7 '18 at 1:19
  • 2
    $\begingroup$ You also need the update for $\mu$, because after you've updated with, say, 10 observations, and are starting to update with the next 10 observations, $\mu_0$ in the above expression is different; it's not the original $\mu_0$ any more, but instead the posterior mean for $\mu$ from the update you've just completed. $\endgroup$ – jbowman Sep 7 '18 at 2:15
  • 1
    $\begingroup$ But you are updating that $\mu_0$ at every step in your repeated-update path to the final procedure? $\endgroup$ – jbowman Sep 7 '18 at 2:39
  • 1
    $\begingroup$ ... and $n_0$ is updated (by $+10$) after every $10$ observations and $\bar{x}$ is the sample mean of just the $10$ observations? (Sort of going down the list of places where mistakes might be made...) $\endgroup$ – jbowman Sep 7 '18 at 2:48
  • 1
    $\begingroup$ You should probably write out the update of each parameter that is updated, e.g., how $\alpha$ is updated ($alpha' \leftarrow \alpha + n/2$), how $\beta$ is updated, how $n_0$ is updated, how $\mu_0$ is updated, ... That will make it easier to figure out which ones aren't being updated in a way that's consistent between the two procedures. $\endgroup$ – jbowman Sep 7 '18 at 3:10
1
$\begingroup$

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.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Thanks for putting this together! This helped confirm a hunch that I had. The likelihood that gave me the inconsistent posteriors was not from a normal distribution. However, when I tested my implementation how you test yours, i.e. by drawing from a normal and updating using those values, I do get consistent posteriors. This makes sense to me. I just need to make sure my likelihoods are in fact drawn from normal distributions. Thanks again! This was very helpful. $\endgroup$ – Robert Sep 8 '18 at 17:58
  • $\begingroup$ Actually, my hunch was wrong. So, I get differences even using your r code. I don't get them consistently, but I do get them using your code as is if I run it a few times. And, the differences are greater and more likely the fewer draws from the Normal, viz. the smaller n in the update. Does that seem right? $\endgroup$ – Robert Sep 8 '18 at 18:21
  • 1
    $\begingroup$ Huh, I didn't get any. Let me make some more extensive tests. $\endgroup$ – jbowman Sep 8 '18 at 19:22
  • 1
    $\begingroup$ I didn't get any differences that weren't small enough to be attributable to roundoff error even with 10 updates of 10 each (see added code section in the answer above). Can you give me the parameters you were initializing with as well as the sample sizes you were using? $\endgroup$ – jbowman Sep 8 '18 at 19:35
  • $\begingroup$ Yeah, you're right. Sorry, I should have paid attention to the magnitude of differences. And your final thought was spot on. I had originally specified $\nu$ as a Long since the size of the sample would always be an integer. The consequence, though, as you rightly pointed out, was that divisions by $\nu$ resulted in rounded numbers. Thanks again for all of your help! $\endgroup$ – Robert Sep 8 '18 at 23:32
1
$\begingroup$

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}$.

|cite|improve this answer
$\endgroup$

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.