1 In this exercise the posterior (target) distribution is

\[ \pi(\theta|y_{1:n}) \propto \prod_{i=1}^n p_N(y_i|\theta_1+\theta_2^2,\sigma^2), \] where \(\theta=(\theta_1,\theta_2)\).

posterior = function(thetas){
  sum(dnorm(y,thetas[1]+thetas[2]^2,sigy,log=TRUE))+
  dnorm(thetas[1],0,sigth,log=TRUE)+
  dnorm(thetas[2],0,sigth,log=TRUE)
}

2 Simulating some data

set.seed(31415)
n     = 20
sigy  = 2
sigth = 1
y     = rnorm(n,1,sigy)
y
##  [1]  4.2940258 -1.2239743 -2.0970912 -0.4967569  4.0528222 -1.2758508
##  [7]  2.5250727  2.7833250 -1.2507839  3.4323365 -0.4404296  0.8419455
## [13]  3.6728787  5.3518124 -0.2399929 -1.5464641 -2.5428390  3.8169362
## [19]  0.8109975 -0.3364036

3 Computing (on a grid) the joint and marginal posterior distributions

N   = 200
th1 = seq(-3,3,length=N)
th2 = seq(-3,3,length=N)
post = matrix(0,N,N)
for (i in 1:N)
  for (j in 1:N)
    post[i,j] = posterior(c(th1[i],th2[j]))
post = exp(post-max(post))    
    
post.th1 = apply(post,1,sum)
post.th1 = post.th1/sum(post.th1)/(th1[2]-th1[1])
sum(post.th1*(th1[2]-th1[1]))
## [1] 1
post.th2 = apply(post,2,sum)
post.th2 = post.th2/sum(post.th2)/(th2[2]-th2[1])
sum(post.th2*(th2[2]-th2[1]))
## [1] 1

4 Plots

par(mfrow=c(1,3))
image(th1,th2,post,col=cm.colors(60),main="",xlab=expression(theta[1]),ylab=expression(theta[2]))
contour(th1,th2,post,add=TRUE,drawlabels=FALSE)
plot(th1,post.th1,type="l",xlab=expression(theta[1]),ylab="Density")
plot(th2,post.th2,type="l",xlab=expression(theta[2]),ylab="Density")

5 Sampling importance resampling (SIR)

time.SIR = system.time({
set.seed(54321)
N  = 30000
N  = 2400000
N1 = 20000

theta1.d = runif(N,-4,3)
theta2.d = runif(N,-3,3)
w = rep(0,N)
for (i in 1:N)
w[i] = posterior(c(theta1.d[i],theta2.d[i]))
w = exp(w-max(w))
ind = sample(1:N,size=N1,replace=TRUE,prob=w)
theta1.dp = theta1.d[ind]
theta2.dp = theta2.d[ind]
})
time.SIR
##    user  system elapsed 
##  18.978   0.170  19.165
par(mfrow=c(1,3))
plot(theta1.dp,theta2.dp,xlab=expression(theta[1]),ylab=expression(theta[2]),col=grey(0.8),pch=16,xlim=range(th1,theta1.dp),ylim=range(th2,theta2.dp))
contour(th1,th2,post,drawlabels=FALSE,add=TRUE)
hist(theta1.dp,prob=TRUE,xlab=expression(theta[1]),main="")
lines(th1,post.th1,lwd=2)
hist(theta2.dp,prob=TRUE,xlab=expression(theta[2]),main="")
lines(th2,post.th2,lwd=2)

6 RANDOM-WALK METROPOLIS-HASTINGS ALGORITHM

time.rw = system.time({
set.seed(2345)
theta1=-2
theta2=0
M0 = 730000
M = 20000
niter = M0+M
draws = matrix(0,niter,2)
for (iter in 1:niter){
theta1s = rnorm(1,theta1,0.1)
theta2s = rnorm(1,theta2,0.1)
logaccept = min(0,posterior(c(theta1s,theta2s))-posterior(c(theta1,theta2)))
if (log(runif(1))<logaccept){
  theta1 = theta1s
  theta2 = theta2s
}
draws[iter,] = c(theta1,theta2)
}
})
par(mfrow=c(2,3))
plot(theta1.dp,theta2.dp,xlab=expression(theta[1]),ylab=expression(theta[2]),col=grey(0.8),pch=16,xlim=range(th1,theta1.dp),ylim=range(th2,theta2.dp))
contour(th1,th2,post,drawlabels=FALSE,add=TRUE)
title("SIR")
hist(theta1.dp,prob=TRUE,xlab=expression(theta[1]),main="")
lines(th1,post.th1,lwd=2)
hist(theta2.dp,prob=TRUE,xlab=expression(theta[2]),main="")
lines(th2,post.th2,lwd=2)

plot(draws,xlab=expression(theta[1]),ylab=expression(theta[2]),col=grey(0.8),pch=16,
xlim=range(th1,draws[,1]),ylim=range(th2,draws[,2]))
contour(th1,th2,post,drawlabels=FALSE,add=TRUE)
title("RW-MH")
hist(draws[,1],prob=TRUE,xlab=expression(theta[1]),main="")
lines(th1,post.th1,lwd=2)
hist(draws[,2],prob=TRUE,xlab=expression(theta[2]),main="")
lines(th2,post.th2,lwd=2)

# MCMC with adapting strategies
#install.packages("adaptMCMC")
library(adaptMCMC)
## Loading required package: parallel
## Loading required package: coda
## Loading required package: Matrix
time.adMCMC = system.time({
    samp <- MCMC(posterior,n=niter,init=c(-2,0),scale=c(0.01,0.01),adapt=TRUE,acc.rate=0.15)
})
##   generate 750000 samples
samp = samp$samples[(M0+1):niter,]

7 Trace plots and autocorrelations

par(mfrow=c(2,4))
ts.plot(draws[,1],xlab="Iterations",ylab=expression(theta[1]),main="RW-MH")
ts.plot(draws[,2],xlab="Iterations",ylab=expression(theta[2]),main="RW-MH")
acf(draws[,1],lag.max=100,main="RW-MH")
acf(draws[,2],lag.max=100,main="RW-MH")
ts.plot(samp[,1],xlab="Interations",ylab=expression(theta[1]),main="With adaption")
ts.plot(samp[,2],xlab="Interations",ylab=expression(theta[2]),main="With adaption")
acf(samp[,1],lag.max=100,main="with adaption")
acf(samp[,2],lag.max=100,main="with adaption")

8 Comparison

par(mfrow=c(3,3),mai=c(0.6,0.6,0.2,0.2))
plot(theta1.dp,theta2.dp,xlab=expression(theta[1]),ylab=expression(theta[2]),col=grey(0.8),pch=16,xlim=range(th1,theta1.dp),ylim=range(th2,theta2.dp))
contour(th1,th2,post,drawlabels=FALSE,add=TRUE)
title(paste("SIR: N=",N," - N1=",N1,sep=""))
legend("topright",legend=paste(round(time.SIR[1],3),"secs",sep=""),bty="n")
hist(theta1.dp,prob=TRUE,xlab=expression(theta[1]),main="")
lines(th1,post.th1,lwd=2)
hist(theta2.dp,prob=TRUE,xlab=expression(theta[2]),main="")
lines(th2,post.th2,lwd=2)

plot(draws,xlab=expression(theta[1]),ylab=expression(theta[2]),col=grey(0.8),pch=16,
xlim=range(th1,draws[,1]),ylim=range(th2,draws[,2]))
contour(th1,th2,post,drawlabels=FALSE,add=TRUE)
title(paste("RW-MH: M0=",M0," - M=",M,sep=""))
legend("topright",legend=paste(round(time.rw[1],3),"secs",sep=""),bty="n")
hist(draws[,1],prob=TRUE,xlab=expression(theta[1]),main="")
lines(th1,post.th1,lwd=2)
hist(draws[,2],prob=TRUE,xlab=expression(theta[2]),main="")
lines(th2,post.th2,lwd=2)

plot(samp,xlab=expression(theta[1]),ylab=expression(theta[2]),col=grey(0.8),pch=16,
xlim=range(th1,draws[,1]),ylim=range(th2,draws[,2]))
contour(th1,th2,post,drawlabels=FALSE,add=TRUE)
title(paste("ad-MCMC: M0=",M0," - M=",M,sep=""))
legend("topright",legend=paste(round(time.adMCMC[1],3),"secs",sep=""),bty="n")
hist(samp[,1],prob=TRUE,xlab=expression(theta[1]),main="")
lines(th1,post.th1,lwd=2)
hist(samp[,2],prob=TRUE,xlab=expression(theta[2]),main="")
lines(th2,post.th2,lwd=2)