`panel10pt2.fn` <-
function(ni=3000,nb=2000,nthin=1,nchain=3,nz=160){
#
# Jolly-Seber Model w/Individual Heterogeneity -- This R script fits the 
# occupancy formulation of the Jolly-Seber model induced by data 
# augmentation to the dipper data, and allows for  individual 
# heterogeneity in survival probability. This script produced the 
# results in Table 10.3. The Markov chains mix slowly, so a 
# long run is necessary to reduce Monte Carlo error. Note also that
# results in Table 10.3 were run with nz=100 which probably was not high enough


library("R2WinBUGS")

y<-dipper.data$y
# set pre-capture NA to 0 -- this is important
y[is.na(y)]<-0
y<-rbind(y,matrix(rep(0,ncol(y)*nz),ncol=ncol(y),byrow=TRUE))
M<-nrow(y)
T<-ncol(y)
z<-y



sink("jsmodel.txt")
cat("
model {

for(j in 1:(T-1)){
phi0[j]~dunif(0,1)
lphi[j]<-log(phi0[j]/(1-phi0[j]))
}

ptmp[1]<- 1
ptmp[T]<- 1
for(i in 2:(T-1)){
ptmp[i]~dunif(0,1)
}

sigma~dunif(0,10)
tau<-1/(sigma*sigma)

for(j in 1:M){
eta[j]~dnorm(0,tau)
for(i in 1:(T-1)){
logit(phi[j,i])<-lphi[i]  + eta[j]
}
}


# one parameterization is to give each gamma a uniform(0,1) prior
# the other way is to put the dirichlet prior on the multinomial equivalent

a[1]<- 1.4/14
a[2]<- 1.4/13
a[3]<- 1.4/12
a[4]<- 1.4/11
a[5]<- 1.4/10
a[6]<- 1.4/9
a[7]<- 1.4/8
for(i in 1:T){
b[i]<-(1.4-a[i])
gamma[i]~dbeta(a[i],b[i])
lgamma[i]<- log(gamma[i]/(1-gamma[i]))
gammaSA[i]<- cprob[i]/psi.imp
}

cprob[1]<- gamma[1]
cprob[2]<- gamma[2]*(1-gamma[1])
cprob[3]<- gamma[3]*(1-gamma[2])*(1-gamma[1])
cprob[4]<- gamma[4]*(1-gamma[3])*(1-gamma[2])*(1-gamma[1])
cprob[5]<- gamma[5]*(1-gamma[4])*(1-gamma[3])*(1-gamma[2])*(1-gamma[1])
cprob[6]<- gamma[6]*(1-gamma[5])*(1-gamma[4])*(1-gamma[3])*(1-gamma[2])*(1-gamma[1])
cprob[7]<- gamma[7]*(1-  gamma[6])*(1-gamma[5])*(1-gamma[4])*(1-gamma[3])*(1-gamma[2])*(1-gamma[1])
psi.imp<- sum(cprob[1:7])

for(i in 1:M){
z[i,1]~dbern(gamma[1])
mu[i]<-z[i,1]*ptmp[1]
y[i,1]~dbern(mu[i])
recruitable[i,1]<-1

for(j in 2:T){
 survived[i,j]<- phi[i,j-1]*z[i,j-1]
 recruitable[i,j]<- recruitable[i,(j-1)]*(1-z[i,j-1])
muz[i,j]<-  survived[i,j] + gamma[j]*recruitable[i,j] 
 z[i,j]~dbern(muz[i,j])
 muy[i,j]<-z[i,j]*ptmp[j]
 y[i,j]~dbern(muy[i,j])
}
} 


for(j in 1:M){
recruit[j,1]<-z[j,1]
for(i in 2:T){
recruit[j,i]<-z[j,i-1]*(1-z[j,i])
}
}

for(i in 1:T){
N[i]<-sum(z[1:M,i])
B[i]<-sum(recruit[1:M,i])
}

for(i in 1:M){
Nind[i]<-sum(z[i,1:T])
Nalive[i]<-1-equals(Nind[i],0)
}

Nsuper<-sum(Nalive[1:M])

}
",fill=TRUE)
sink()    

data <- list ("y","M","T")
inits <- function(){
  list ( phi0=runif(T-1,.2,.9),z=z,sigma=1)
}

parameters <- c("ptmp","phi0","N","Nsuper","gamma","gammaSA","cprob")

out <- bugs (data, inits, parameters, "jsmodel.txt", n.thin=nthin,n.chains=nchain, n.burnin=nb,n.iter=ni,debug=TRUE)


return(out)


}