MCMC Exercise 1

1. Make a 3-panel figure to illustrate Beta-Binomial conjugate updating. The first panel should show a Beta prior distribution, the second panel should show the likelihood function, and the third panel should show the posterior distribution. The posterior should be calculated using the conjugate relationship between the beta prior and the binomial likelihood. Data: 100 acorns are planted and 33 germinate.

# Beta Prior
a<-40; b<-20
x<-seq(from=0,to=1,len=1000)
priorF<-dbeta(x,shape1=a,shape2 = b)
par(mfrow=c(3,1),mar=c(3.0,3.5,1.5,1.5)+1)
plot(x,priorF,xlab='P',ylab='Density',type='l')

# Binomial Likelihood
likF<-dbinom(33,size=100,prob = x)
plot(x,likF,xlab='P',ylab='Density',type='l')

# Conjugate Beta Posterior
aPrime<-a+33; bPrime<-b+100-33
x<-seq(from=0,to=1,len=1000)
postF<-dbeta(x,shape1=aPrime,shape2 = bPrime)
plot(x,postF,xlab='P',ylab='Density',type='l')

2. Write code to sample from this same posterior distribution using the Metropolis algorithm (MCMC). Generate 10,000 samples from this posterior.

Likeilihood function:

likBinom<-function(p,s,n){
  lik<-dbinom(s,size=n,prob = p)
  # nllik<- -sum(dbinom(s,size=n,prob = parVec,log=TRUE))
  # cat("nllik= ",nllik,sep=" ",fill=T);cat(" ",sep=" ",fill=T)
  return(lik)
}

Prior function:

priorBeta<-function(p,a,b){
  prob<-dbeta(p,shape1=a,shape2 = b)
  # nllik<- -sum(dbinom(s,size=n,prob = parVec,log=TRUE))
  # cat("nllik= ",nllik,sep=" ",fill=T);cat(" ",sep=" ",fill=T)
  return(prob)
}

Proposal distribution:

proDist<-function(currentP,sd){
  p<-rnorm(n=1,mean=currentP,sd=sd)
  if(p>0 & p<1) return(p)
  proDist(currentP,sd)
}
  
proDist(0.2,0.2)

## [1] 0.2739407

Metropolis Algorithym

1. Intialize the chain at some initial value of p

2. Generate a proposed value of p

3. Calculate acceptabce probability r

4. Accept proposed value of p with probability r

5. Repeat n times

Coding this in R:

nIter<-10000
mySD<-0.2; a<-1;b<-1
p<-vector()
p[1]<-0.5
for(i in 1:nIter){
  proposedP<-proDist(p[i],mySD)
  r<-min(1, ( likBinom(proposedP,33,100)*priorBeta(proposedP,a,b) ) / ( likBinom(p[i],33,100) *priorBeta(p[i],a,b) ) )
  if(runif(1)<r) {
    p[i+1]<-proposedP
  } else {
    p[i+1]<-p[i]
  } 
}

hist(p[-1])

plot(1:length(p),p)

3. Use the R coda library to create a trace plot, density plot, autocorrelation plot, and an estimate of the 95% HDI from the chain produced in 2.

library(coda)
pMCMC<-as.mcmc(p)
densplot(pMCMC)

HPDinterval(pMCMC)

##          lower     upper
## var1 0.2490604 0.4354726
## attr(,"Probability")
## [1] 0.950005

autocorr.plot(pMCMC)

plot(pMCMC)