rm(list=ls()); set.seed(20210406);
# Discussed in class on Tuesday, April 6
library(MASS)
# This function takes as input the hyperparameters for the
# Normal-Inverse-Wishart conjugate prior distribution to
# a multivariate normal model, and returns a random sample
# (indep draws) from the joint posterior of the mean
# vector mu (mu.sim) and covarance matrix Sigma (Sigma.sim)
Post.sim.MVN <- function(Sim.size, y, mu0, k0, nu0, Lambda0) { T <- Sim.size n <- dim(y)[1]; d <- dim(y)[2]; ybar <- apply(y, 2, mean) S <- (n-1) * cov(y) k.n <- k0 + n mu.n <- (k0*mu0 + n*ybar) / k.n nu.n <- nu0 + n Lambda.n <- Lambda0 + S + k0*n/k.n * outer(ybar-mu0, ybar-mu0) Lambda.n.inv <- solve(Lambda.n); rm(Lambda.n); mu.sim <- matrix(rep(NA, T*d), T, d) Sigma.sim <- matrix(rep(NA, T*d^2), T, d^2) for(t in 1:T) { Sigma <- solve(rWishart(1, nu.n, Lambda.n.inv)[,,1]) mu <- mvrnorm(1, mu.n, 1/k.n*Sigma) mu.sim[t,] <- mu; rm(mu); Sigma.sim[t,] <- c(Sigma); rm(Sigma); } list(mu.sim=mu.sim, Sigma.sim=Sigma.sim) } # Read in data file <- "http://www2.stat.duke.edu/~pdh10/FCBS/Exercises/agehw.dat" y <- read.table(file=file, header=T); rm(file); # This data file contains the ages for 100 randomly selected # opposite-sex married couples. # Our inferential goal is to summarize the posterior distributions # of mu_h (mean height of husbands), mu_w (mean height of wives), # and the correlation coefficient rho # This is Exercise 7.4 from Hoff (2009) # Prior mean, prior sample size for normal prior on mean vector # Prior scale matrix, prior df for inverse Wishart prior on Sigma mu0 <- c(45,40); k0 <- 1; nu0 <- 3; Lambda0 <- diag(2); # Now the simulations foo <- Post.sim.MVN(Sim.size=1000, y, mu0, k0, nu0, Lambda0) mu.sim <- foo$mu.sim; Sigma.sim <- foo$Sigma.sim; rm(foo); plot(mu.sim, xlab="mu_h", ylab="mu_w", main="Joint posterior of (mu_h, mu_w)" ) rho.sim <- Sigma.sim[,2] / sqrt(Sigma.sim[,1] * Sigma.sim[,4]) hist(rho.sim, breaks=30, freq=F); lines(density(rho.sim), col="red"); Results <- apply(mu.sim, 2, quantile, prob=c(.025,.25,.5,.75,.975)) Results <- t(Results); rownames(Results) <- c("ageh", "agew"); Results