######### 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') # compute posterior mean and quantiles colMeans(Phi) apply(Phi,2,quantile,probs=c(.025,.975))