# Blood coagulation time example from Section 11.6 of Gelman et al
rm(list=ls()); set.seed(20210312);
chain.length <- 2000 # Key in data J <- 4 yA <- c(62, 60, 63, 59) yB <- c(63, 67, 71, 64, 65, 66) yC <- c(68, 66, 71, 67, 68, 68) yD <- c(56, 62, 60, 61, 63, 64, 63, 59) n <- c(length(yA), length(yB), length(yC), length(yD)) ybar <- c(mean(yA), mean(yB), mean(yC), mean(yD)) s <- c(sd(yA), sd(yB), sd(yC), sd(yD)) rm(yA, yB, yC, yD) # Gibbs update functions theta.update <- function(mu, sigma, tau, J, n, ybar) { V.theta <- 1 / (1/tau^2 + n/sigma^2) theta.hat <- V.theta * (mu/tau^2 + n*ybar/sigma^2) rnorm(J, mean=theta.hat, sd=sqrt(V.theta)) } mu.update <- function(theta, tau, J) { mu.hat <- mean(theta) rnorm(1, mean=mu.hat, sd=tau/sqrt(J)) } sigma.update <- function(theta, n, ybar, s) { sigma2.hat <- sum((n-1)*s^2 + n*(ybar-theta)^2) / sum(n) sigma2 <- sum(n) * sigma2.hat / rchisq(1, df=sum(n)) sqrt(sigma2) } tau.update <- function(J, theta, mu) { tau2.hat <- sum((theta-mu)^2) / (J-1) tau2 <- (J-1) * tau2.hat / rchisq(1, df=J-1) sqrt(tau2) } # Putting it together # Inputs: chain length, data, starting values build.chain <- function(chain.length, J, n, y, s, theta0, mu0, sigma0, tau0) { T <- chain.length theta.chain <- matrix(NA, T, J) mu.chain <- rep(NA, T); sigma.chain <- rep(NA, T); tau.chain <- rep(NA, T); theta <- theta0; mu <- mu0; sigma <- sigma0; tau <- tau0; for(t in 1:T) { theta <- theta.update(mu, sigma, tau, J, n, ybar) mu <- mu.update(theta, tau, J) sigma <- sigma.update(theta, n, ybar, s) tau <- tau.update(J, theta, mu) theta.chain[t,] <- theta; mu.chain[t] <- mu; sigma.chain[t] <- sigma; tau.chain[t] <- tau; } list(theta.chain=theta.chain, mu.chain=mu.chain, sigma.chain=sigma.chain, tau.chain=tau.chain) } # Run a single chain # Initial values theta0 <- ybar; mu0 <- mean(theta0); sigma0 <- sqrt(mean(s^2)); tau0 <- sd(ybar); chain <- build.chain(chain.length, J, n, y, s, theta0, mu0, sigma0, tau0) theta.chain <- chain$theta.chain; mu.chain <- chain$mu.chain; sigma.chain <- chain$sigma.chain; tau.chain <- chain$tau.chain; rm(chain) Results <- matrix(NA, J+3, 5) probs <- c(.025, .25, .50, .75, .975) for(j in 1:J) { Results[j, ] <- quantile(theta.chain[,j], probs=probs) } Results[J+1,] <- quantile(mu.chain, probs=probs) Results[J+2, ] <- quantile(sigma.chain, probs=probs) Results[J+3, ] <- quantile(tau.chain, probs=probs) rownames(Results) <- c(paste("theta",1:J,sep=""), "mu", "sigma", "tau") colnames(Results) <- paste(probs*100, "pct", sep="") round(Results, 1)