`panel6pt2.fn` <-
function(ni=5000,nb=1000,nthin=2,nc=2,nz=225){
# This R script fits Model Mh in WinBUGS using data augmentation.
# (see Panel 6.2)
# The data are the breeding bird survey data for the species
# richness estimation problem
library("R2WinBUGS")

start.time = Sys.time()

nx<-c(15, 7, 5, 2, 5, 6, 3, 2, 1, 2, 
   4, 4, 0, 1, 0, 2, 3, 1, 1, 1,
   1, 0, 0, 0, 1, 1, 0, 1, 0, 1,
   0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
J<-length(nx)
nind<-sum(nx)
y<-rep(1:50,nx)
nz<-250
y<-c(y,rep(0,nz))

sink("winbugsmodel.txt")
cat("
model {
  psi~dunif(0,1)
  mu.p~dnorm(0,0.001)
  tau~dgamma(.001,.001)
  sigma<-sqrt(1/tau)
  for(i in 1:(nind+nz)){
    z[i]~dbin(psi,1)
    eta[i]~dnorm(mu.p,tau)     
    logit(p[i])<- eta[i]
    mu[i]<-p[i]*z[i]
    y[i] ~ dbin(mu[i],J)
}
N<-sum(z[1:(nind+nz)])
}
",fill=TRUE)
sink()

data <- list ("y","nind","nz","J")
zst<-c(rep(1,nind),rbinom(nz,1,.25))
inits <- function(){
  list (mu.p=rnorm(1,-2,.2),tau=runif(1,1,2),z=zst)
}
parameters <- c("N","mu.p","sigma")
out <- bugs (data, inits, parameters, "winbugsmodel.txt", n.thin=nthin,n.chains=nc, n.burnin=nb,n.iter=ni,debug=FALSE)

end.time = Sys.time()
 elapsed.time = difftime(end.time, start.time, units='mins')
 cat(paste(paste('Posterior computed in ', elapsed.time, sep=''), ' minutes\n', sep=''))





out


}