`panel3pt8bayes.fn` <-
function(niter=2000,nburn=1000,nthin=1,nc=3){
# fits occupancy model with detection and occurrence covariates
# using WinBUGS. Uses Swiss MHB willow tit data
# See Section 3.6.2 in the book.

library("R2WinBUGS")
source("utilfns.Rd")

data<- read.table("wtmatrix.csv",header=TRUE,sep=",",na.strings=c("NA"))

y<-as.matrix(data[,c("y.1","y.2","y.3")])
M<-nrow(y)
J<- ncol(y)

length<-as.vector(data[,"length"])
forest<-as.vector(scale(data[,"forest"],center=TRUE))
elev<-as.vector(scale(data[,"elev"],center=TRUE))
elev2<-elev*elev
date<-as.matrix(data[,c("day.1","day.2","day.3")])

mdate<-mean(date,na.rm=TRUE)
sddate<-sqrt(var(date[1:length(date)],na.rm=TRUE))
date<-(date-mdate)/sddate
date2<-date*date

dur<-as.matrix(data[,c("dur.1","dur.2","dur.3")])
intensity<-dur/length
mint<-mean(intensity,na.rm=TRUE)
sdint<-sqrt(var(intensity[1:length(intensity)],na.rm=TRUE))
intensity<-(intensity-mint)/sdint
length<- (length-mean(length))/sqrt(var(length))


date[is.na(y)]<-date2[is.na(y)]<-intensity[is.na(y)]<-0

data <- list ( "y","M","J","forest","elev","elev2","date","date2","intensity")

inits <- function()
  list ( z=rbinom(M,1,.4),psi0=runif(1),b1=rnorm(1),b2=rnorm(1),
b3=rnorm(1),p0=runif(1),a1=rnorm(1),a2=rnorm(1),a3=rnorm(1))

parameters <- c("b0","b1","b2","b3","a0","a1","a2","a3","p0","psi0")

sink("model.txt")
cat("

model { 
                               # prior distributions
        p0~dunif(0,1)
        a0<-log(p0/(1-p0))
        a1~dnorm(0,.001)
        a2~dnorm(0,.001)
        a3~dnorm(0,.001)

        psi0 ~ dunif(0,1)      
        b0 <-log(psi0/(1-psi0))
        b1 ~ dnorm(0,.001) 
        b2 ~ dnorm(0,.001) 
        b3 ~ dnorm(0,.001) 

  for(i in 1:M){

     z[i] ~ dbin(psi[i],1) # STATE MODEL
     logit(psi[i]) <- b0 + b1*elev[i] + b2*elev2[i] + b3*forest[i]
     
        for(t in 1:J){
           mu[i,t]<-a0 + a1*date[i,t] + a2*date2[i,t] + a3*intensity[i,t]
           p[i,t]<-exp(mu[i,t])/(1+exp(mu[i,t]))

           muy[i,t]<-z[i]*p[i,t]     #
           y[i,t] ~ dbin(muy[i,t],1) # OBSERVATION MODEL
}
}


      }
",fill=TRUE)
sink()


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

out






}