`panel6pt5.fn` <-
function(ni=6000,nb=1000,nthin=3,nc=4,nz=125){
#This function fits a "double-observer" model with an individual
# covariate influencing detection probability.  Analysis proceeds using
# data augmentation in WinBUGS. The data are the aerial waterfowl
# survey data (U.S. Fish and Wildlife Service) and the individual 
# covariate is group size (number of ducks in the group). The two
# observers were front-seat (pilot) and back-seat (observer)

library("R2WinBUGS")

data<-source("malldata.R")$value
#0,1  1-7
#1,0  8-14
#1,1  15-21
x01<-apply(data[,1:7],2,sum)
x10<-apply(data[,8:14],2,sum)
x11<-apply(data[,15:21],2,sum)
X<-rbind(x01,x10,x11)

# this mess of code creates the multinomial trial representation of the data
ncap<-matrix(0,nrow=sum(X),ncol=4)
rid<-c( 1:sum(X[1,]), (sum(X[1,])+1):(sum(X[1,]) + sum(X[2,])), (sum(X[1,])+sum(X[2,])+1):nrow(ncap))
cid<-c( rep(1,sum(X[1,])),rep(2,sum(X[2,])),rep(3,sum(X[3,])) )
ncap[cbind(rid,cid)]<-1
gs<-c( rep(1:7,X[1,]), rep(1:7,X[2,]),  rep(1:7,X[3,]))
nind<-length(gs)

# augment the encounter histories with histories of the last type 
# which are those not seen by either observer (0,0,0,1).
# the covariate groupsize is augmented with missing values
ncap<-rbind(ncap,matrix(rep(c(0,0,0,1),nz),ncol=4,byrow=TRUE))
gs<-c(gs,rep(NA,nz))

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

psi~dunif(0,1)
logmu.gs~dnorm(0,.001)
mu.gs<-exp(logmu.gs)
mu1.p~dunif(-10,10)
mu2.p~dunif(-10,10)
beta~dunif(-10,10)

for(i in 1:(nind+nz)){
  z[i]~dbin(psi,1)
  gs[i]~dpois(mu.gs)
  logit(p1[i])<- mu1.p + beta*(1+gs[i])
  logit(p2[i])<- mu2.p + beta*(1+gs[i])
  
  cp1[i]<- (1-p1[i])*p2[i]
  cp2[i]<- p1[i]*(1-p2[i])
  cp3[i]<- p1[i]*p2[i]
  cp4[i]<- (1-p1[i])*(1-p2[i])

  mu[i,1]<-z[i]*cp1[i]
  mu[i,2]<-z[i]*cp2[i]
  mu[i,3]<-z[i]*cp3[i]
  mu[i,4]<-z[i]*cp4[i] + (1-z[i])

  ncap[i,1:4]~dmulti(mu[i,1:4],1)
}
 
logit(p1bar)<-mu1.p + beta
logit(p2bar)<-mu2.p + beta

Ng<-sum(z[1:(nind+nz)])
for(i in 1:(nind+nz)){
  tmp[i]<- 1 + gs[i]
  tmp2[i]<-tmp[i]*z[i]
}
Nind<-sum(tmp2[1:(nind+nz)])
D<-Nind/544.5
}
",fill=TRUE)
sink()

data <- list ("ncap","nind","nz","gs")
inits <- function(){
  list (logmu.gs=0,mu1.p=0,mu2.p=0,beta=rnorm(1))

}

parameters <- c("Ng","Nind","beta","mu1.p","mu2.p","mu.gs","p1bar","p2bar","D","psi")

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






out


}