This lab will provide an opportunity to explore MCMC code. Turn in your results to D2L under MCMC_Lab. You are welcome to work together if you’d like, but each person should turn in their own document.
Answer the following questions in this R Markdown document. Note: for this assignment you can supress much of the code by using the echo=F option.
Use the code below for this question
######### First Gibbs Sampler
set.seed(09182017)
### simulate data
num.obs <- 100
mu.true <- 0
sigmasq.true <- 1
y <- rnorm(num.obs,mu.true,sigmasq.true)
mean.y <- mean(y)
var.y <- var(y)
library(LearnBayes) # for rigamma
### initialize vectors and set starting values and priors
num.sims <- 10000
Phi <- matrix(0,nrow=num.sims,ncol=2)
Phi[1,1] <- 0 # initialize theta
Phi[1,2] <- 1 # initialize (sigmasq)
mu.0 <- 0
tausq.0 <- 1
nu.0 <- 1
sigmasq.0 <- 10
for (i in 2:num.sims){
# sample theta from full conditional
mu.n <- (mu.0 / tausq.0 + num.obs * mean.y / Phi[(i-1),2]) / (1 / tausq.0 + num.obs / Phi[(i-1),2] )
tausq.n <- 1 / (1/tausq.0 + num.obs / Phi[(i-1),2])
Phi[i,1] <- rnorm(1,mu.n,sqrt(tausq.n))
# sample (1/sigma.sq) from full conditional
nu.n <- nu.0 + num.obs
sigmasq.n.theta <- 1/nu.n*(nu.0*sigmasq.0 + sum((y - Phi[i,1])^2))
Phi[i,2] <- rigamma(1,nu.n/2,nu.n*sigmasq.n.theta/2)
}
# plot joint posterior
plot(Phi[1:5,1],1/Phi[1:5,2],xlim=range(Phi[,1]),ylim=range(1/Phi[,2]),pch=c('1','2','3','4','5'),cex=.8,
ylab=expression(sigma[2]), xlab = expression(theta), main='Joint Posterior',sub='first 5 samples')
plot(Phi[1:10,1],1/Phi[1:10,2],xlim=range(Phi[,1]),ylim=range(1/Phi[,2]),pch=as.character(1:15),cex=.8,
ylab=expression(sigma[2]), xlab = expression(theta), main='Joint Posterior',sub='first 10 samples')
plot(Phi[1:100,1],1/Phi[1:100,2],xlim=range(Phi[,1]),ylim=range(1/Phi[,2]),pch=16,col=rgb(0,0,0,1),cex=.8,
ylab=expression(sigma[2]), xlab = expression(theta), main='Joint Posterior',sub='first 100 samples')
plot(Phi[,1],1/Phi[,2],xlim=range(Phi[,1]),ylim=range(1/Phi[,2]),pch=16,col=rgb(0,0,0,.25),cex=.8,
ylab=expression(sigma[2]), xlab = expression(theta), main='Joint Posterior',sub='all samples')
points(0,1,pch='X',col='red',cex=2)
# plot marginal posterior of theta
hist(Phi[,1],xlab=expression(theta),main=expression('Marginal Posterior of ' ~ theta),probability=T)
abline(v=mu.true,col='red',lwd=2)
# plot marginal posterior of sigmasq
hist(Phi[,2],xlab=expression(sigma[2]),main=expression('Marginal Posterior of ' ~ sigma[2]),probability=T)
abline(v=sigmasq.true^2,col='red',lwd=2)
# plot trace plots
plot(Phi[,1],type='l',ylab=expression(theta), main=expression('Trace plot for ' ~ theta))
abline(h=mu.true,lwd=2,col='red')
plot(Phi[,2],type='l',ylab=expression(sigma[2]), main=expression('Trace plot for ' ~ sigma[2]))
abline(h=sigmasq.true^2,lwd=2,col='red')Choose and explore a few different starting points for \(\phi\). Summarize your findings. How does impact your graphical displays?
Set the starting point for \(\phi\) at (10,10). Approximately how many iterations do you think it takes for the sampler to “find” the posterior distribution?
Interpret the prior values specified in the code. In other words, how can we interpret these values. Do these values seem reasonable?
Change your simulation settings such that mu.true <- 100. Run the code and assess the findings. How influential was the prior?
With this setting, (mu.true <- 100), choose and defend new priors for \(\theta\) and \(\sigma^2\). Include graphical displays that indicate the location of your posterior distribution and the true values.