`panel9pt1.fn` <-
function(ni=6000,nb=2000,nt=1,nc=3){
# this script fits the model given in Panel 9.1 to the crossbill data from
# the Swiss bird survey.  The model assumes p = 1.  The survey consists of
# 3 replicate samples and only the 2nd sample is used here (since we assume
# that p = 1).

library("R2WinBUGS")

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

psi~dunif(0,1)
for(i in 1:(nyear-1)){
gamma[i]~dunif(0,1)
phi[i]~dunif(0,1)
}

for(i in 1:nsite){
z[i,1]~dbern(psi)
for(t in 2:nyear){
muZ[i,t]<- z[i,t-1]*phi[t-1] + (1-z[i,t-1])*gamma[t-1]
z[i,t]~dbern(muZ[i,t])
}
}

lambda[1]<-psivec[2]/psivec[1]
lambda[2]<-psivec[3]/psivec[2]
lambda[3]<-psivec[4]/psivec[3]

psivec[1]<-psi
for(t in 2:nyear){
psivec[t]<-psivec[t-1]*phi[t-1] + (1-psivec[t-1])*gamma[t-1]
turnover[t-1]<- ((1 - psivec[t-1]) * gamma[t-1])/( (1 - psivec[t-1]) * gamma[t-1] + phi[t-1]*psivec[t-1])
}
}
",fill=TRUE)
sink()

data<-crossbill.data
y<-data$data
y[y>1]<-1

nsite<-dim(y)[1]
nrep<-3
nyear<-4
# extract middle sample
z<-y[,2,]
data <- list ( "z","nsite","nyear")
gst<-runif(3)
inits <- function()
  list (gamma=gst)
parameters <- c("gamma","phi","turnover","psivec","lambda")
out <- bugs (data, inits, parameters, "model.txt", n.thin=nt,n.chains=nc, n.burnin=nb,n.iter=ni,debug=FALSE)
return(out)



}