## Source:   chap7.txt (with modifications) in
## http://www.stat.colostate.edu/computationalstatistics/datasets/computationalstatistics.zip

#__________________________________________________________________

# MCMC:  Gibbs sampler
# Fur Seal Pup Capture-Recapture analysis, Example 7.6, page 212 in
# Givens & Hoeting, 2nd edition

# Key idea:  apply a Gibbs sampler to estimate the population size and
# sampling probabilities in a capture-recapture study.

#_______________________________________________________________________
#DATA
num.caught<-c(30,22,29,26,31,32,35) #number caught each time
                    #num.caught is called "c" in the book
I<-7               #total number of capture occasions
r.param<-84        #total number of unique animals caught
#_______________________________________________________________________
# Gibbs Sampler
# Implements the Gibbs sampler described on page 213-214

log.target<-function(exp.u,alpha,I) {
#Inputs:  exp.u = vector of (exp(u1),exp(u2))
#         alpha = vector of capture probabilities of length I
#         I = number of capture occasions
#
#Output:  log of the target distribution for U1,U2 (see equation (7.27) on page 229)
  I*(lgamma(sum(exp.u))-lgamma(exp.u[1])-lgamma(exp.u[2])) +
    exp.u[1]*sum(log(alpha))+exp.u[2]*sum(log(1-alpha))-sum(exp.u)/1000+
    log(exp.u[1])+log(exp.u[2])
}

#Gibbs set-up
num.its<-10000      #Number of MCMC iterations
N<-rep(0,num.its)   #Vector of N=total sample size realizations
N[1]<-sample(84:500,1)
u<-matrix(0,num.its,2)  #Matrix of u realizations, column 1= u1,
            #column 2= u2
u[1,]<-log(sample(2:3000,2))
alpha<-matrix(0,num.its,I)  #Vector of alpha realizations
set.seed(1)

for (i in 2:num.its) {
   #Equation (7.26) on page 229
   alpha[i,]<-rbeta(I,num.caught+exp(u[i-1,1]),N[i-1]-num.caught+exp(u[i-1,2]))

   #Equation (7.16) on page 217
   N[i]<-rnbinom(1,r.param,1-prod(1-alpha[i,]))+r.param

   #Implements Metropolis-Hastings random walk within Gibbs (eq (7.27) on page 229)
   u[i,]<-u[i-1,]+rnorm(2,c(0,0),c(.085,.085))
   R<-exp(log.target(exp(u[i,]),alpha[i,],I)-
    log.target(exp(u[i-1,]),alpha[i,],I))
   if (R<1)
    if(rbinom(1,1,R)==0)    u[i,]<-u[i-1,]
}

#_______________________________________________________________________
#Summary statistics of Gibbs sampler
#Note:  you may wish to remove iterations as a burn-in period
rbind(round(c("mean.N"=mean(N),"sd.N"=sd(N),quantile(N,c(.025,.975))),2))

mean.alpha=apply(alpha,2,mean)
sd.alpha=apply(alpha,2,sd)
d=apply(alpha,2,quantile,c(.025,.975))
round(rbind(mean.alpha,sd.alpha,d),2)  #alpha summary stats
#_______________________________________________________________________
#Plots of the output

# Note:  you may wish to change the burn-in period if you change the
# number of iterations
burn.in=1:1000

# sample path plot
par(mfrow=c(1,1))
plot(N[-burn.in],type="l")

par(mfrow=c(1,1))                     # sample path plot by parts
plot(N[seq(1001,num.its,by=5)],type="l")

# cusum plot for mean of N
temp <- N[-burn.in]
ntemp <- length(temp)
mtemp <- mean(temp)
temp <- cumsum(temp-mtemp)
par(mfrow=c(1,1))
plot(temp, type="l")

# autocorrelation plot
temp <- N[-burn.in]
acf(temp, lag.max=200, type="correlation", xlab="N");

#Similar to Figure 7.6 on page 214
boxplot(split(apply(alpha[-burn.in,],1,mean),N[-burn.in]),
  ylab="mean capture probability over all captures",xlab="N")

#Alpha plot
boxplot(data.frame(alpha[-burn.in,]),names=1:7,xlab="capture occasion",
  ylab="capture probability")

#Similar to Figure 7.7 on page 215
hist(N[-burn.in],probability=T,xlab="N")