MCMC Exercise 2:

1. Simulated data from a transect from closed forest into a canopy gap. We plant 10 seeds in each quadrat and then simulate survival of 0.1 in the forest and 0.8 in the canopy gap. We will randomly place the gqp-forest bourndary somewhere near the mid-point of the transect.

aPrime<-a+33; bPrime<-b+100-33
x<-seq(from=0,to=1,len=1000)
postF<-dbeta(x,shape1=aPrime,shape2 = bPrime)

library(MCMCpack)

Loading required package: coda
Loading required package: MASS
##
## Markov Chain Monte Carlo Package (MCMCpack)
## Copyright (C) 2003-2017 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
##
## Support provided by the U.S. National Science Foundation
## (Grants SES-0350646 and SES-0350613)
##

1. Make a figure similiar to Fig. 6.4: Beta-Binomial conjugagte updating.

You plant 100 acorns. 33 germinate. What is your posterior estimate for the transition probability?

First, a flat prior using the beta distribution:

# 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. Find the posterior distribution using MCMC

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)

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])

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

plot(cars)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).