Sequential Bayesian updating of Normal-gamma model

585019-sequential-updating-normal-gama.r

# See:
# https://people.eecs.berkeley.edu/%7Ejordan/courses/260-spring10/lectures/lecture5.pdf
# https://www.cs.ubc.ca/%7Emurphyk/Papers/bayesGauss.pdf

post <- function(x, mu0, kappa0, alpha0, beta0) {
  n <- length(x)
  xbar <- mean(x)
  sse <- sum((x - xbar)^2)
  mu <- ((kappa0 * mu0) + (n * xbar)) / (kappa0 + n)
  kappa <- kappa0 + n
  alpha <- alpha0 + n/2
  beta <- beta0 + sse/2 + (kappa0 * n * (xbar - mu0)^2)/(2 * (kappa0 + n))
  return(c(mu, kappa, alpha, beta))
}

set.seed(42)

mu0 <- 0
kappa0 <- 1
alpha0 <- 1
beta0 <- 1 

x <- rnorm(10, 5, 16)

cat("All-at-once:\n")
print(post(x, mu0, kappa0, alpha0, beta0))

cat("\nSequential updating:\n")
for (i in 1:length(x)) {
  params <- post(x[i], mu0, kappa0, alpha0, beta0)
  print(params)
  mu0 <- params[1]
  kappa0 <- params[2]
  alpha0 <- params[3]
  beta0 <- params[4]
}