y = 18  #data: 18 deer, priors: nu ~ Poisson(200), theta ~ Beta(2,38)
N = 5200
nu = numeric(N)
theta = numeric(N)
nu[1] = 200
theta[1] = 0.05
for(i in 1:(N-1)) {
  theta[i+1] = rbeta(1, y+2, nu[i]-y+38)
  nu[i+1] = y + rpois(1, 200*(1-theta[i+1]))
}
par(mfrow=c(2,1))
plot(nu,pch=".")
plot(theta,pch=".")

nu = nu[201:5200]
theta = theta[201:5200]
par(mfrow=c(1,2))
hist(nu)
hist(theta)

mean(nu); sd(nu); mean(theta); sd(theta)

acf(nu); acf(theta)

#priors: nu ~ Poisson(100), theta ~ Beta(2,38)
for(i in 1:(N-1)) {
  theta[i+1] = rbeta(1, y+2, nu[i]-y+38)
  nu[i+1] = y + rpois(1, 100*(1-theta[i+1]))
}
nu = nu[201:5200]; theta = theta[201:5200]; hist(nu); hist(theta)
mean(nu); sd(nu); mean(theta); sd(theta)

#priors: nu ~ Poisson(400), theta ~ Beta(2,38)
for(i in 1:(N-1)) {
  theta[i+1] = rbeta(1, y+2, nu[i]-y+38)
  nu[i+1] = y + rpois(1, 400*(1-theta[i+1]))
}
nu = nu[201:5200]; theta = theta[201:5200]; hist(nu); hist(theta)
mean(nu); sd(nu); mean(theta); sd(theta)

#priors: nu ~ Poisson(200), theta ~ Beta(1,79)
for(i in 1:(N-1)) {
  theta[i+1] = rbeta(1, y+1, nu[i]-y+79)
  nu[i+1] = y + rpois(1, 200*(1-theta[i+1]))
}
nu = nu[201:5200]; theta = theta[201:5200]; hist(nu); hist(theta)
mean(nu); sd(nu); mean(theta); sd(theta)

#priors: nu ~ Poisson(200), theta ~ Beta(1,1)
for(i in 1:(N-1)) {
  theta[i+1] = rbeta(1, y+1, nu[i]-y+1)
  nu[i+1] = y + rpois(1, 200*(1-theta[i+1]))
}
nu = nu[201:5200]; theta = theta[201:5200]; hist(nu); hist(theta)
mean(nu); sd(nu); mean(theta); sd(theta)