### This script simulates data under the model described in Royle and Nichols (2004) and fits that model in WinBUGS. Thanks to Yuichi Yamaura <yamaura@affrc.go.jp> and colleague Tatsuya Amano for providing this script.

par(mfrow=c(3,3))
lambda <- 2			# globally averaged local abundance across sites
nsite <- 100		# number of sites
N <- rpois(nsite,lambda)# abundance of each site generated by Poisson distribution
hist(N)
r <- 0.8			# detection probability of each individual

p <- 1-(1-r)^N		# probability of detecting at least one individuals per visit
plot(p~n,ylim=c(0,1))

v <- 5			# number of visits
y <- rbinom(nsite,v,p)	# number of detections (composed of 1/0) for each site
hist(y)
plot(y~N,ylim=c(0,max(y)))

### hierarchical model
library(R2WinBUGS)
sink("model4.1.txt")
cat("
	model{
		for(i in 1:nsite){
			N[i] ~ dpois(lambda)
			p[i] <- 1-pow(1-r,N[i])
			y[i] ~ dbin(p[i],v)
		}
	lambda ~ dlnorm(0,0.001)
	r ~ dunif(0,1)
	}
",fill=TRUE)
sink()

data <- list("y","nsite","v")
inits <- function()list(N=rpois(nsite,10),lambda=rlnorm(1),r=runif(1))
parameters <- c("lambda","r")
out <- bugs(data,inits,parameters,"model4.1.txt",
		n.chain=3,n.burnin=5000,n.iter=10000,debug=T)
out