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=".")

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

mean(nu[201:5200]); sd(nu[201:5200]); mean(theta[201:5200]);sd(theta[201:5200])

acf(nu[201:5200]); acf(theta[201:5200])

#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]))
}
hist(nu[201:5200]); hist(theta[201:5200])
mean(nu[201:5200]); sd(nu[201:5200]); mean(theta[201:5200]);sd(theta[201:5200])

#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]))
}
hist(nu[201:5200]); hist(theta[201:5200])
mean(nu[201:5200]); sd(nu[201:5200]); mean(theta[201:5200]);sd(theta[201:5200])

#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]))
}
hist(nu[201:5200]); hist(theta[201:5200])
mean(nu[201:5200]); sd(nu[201:5200]); mean(theta[201:5200]);sd(theta[201:5200])

#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]))
}
hist(nu[201:5200]); hist(theta[201:5200])
mean(nu[201:5200]); sd(nu[201:5200]); mean(theta[201:5200]);sd(theta[201:5200])


#compare mixing to Metropolis-Hastings, sampling theta using a normal proposal
# explore convergence diagnostics (compare to original acf plots first)
accept.theta = 0
for(i in 1:(N-1)) {
  theta.star = rnorm(1, theta[i], .05)
  if(theta.star >= 0 & theta.star <= 1) {
    # alpha = (theta.star^(y+1))*(1-theta.star)^(nu[i]-y+37)/((theta[i]^(y+1))*(1-theta[i])^(nu[i]-y+37))
    alpha = (theta.star/theta[i])^(y+1)*((1-theta.star)/(1-theta[i]))^(nu[i]-y+37)
    if(runif(1) < alpha) {
      theta[i+1] = theta.star
      accept.theta = accept.theta + 1
    }
  }
  else theta[i+1] = theta[i]
  nu[i+1] = y + rpois(1, 200*(1-theta[i+1]))
}
par(mfrow=c(2,1))
plot(nu,pch="."); plot(theta,pch=".")
accept.theta/N

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

#try larger step sizes
accept.theta = 0
for(i in 1:(N-1)) {
  theta.star = rnorm(1, theta[i], .06)
  if(theta.star >= 0 & theta.star <= 1) {
    alpha = (theta.star/theta[i])^(y+1)*((1-theta.star)/(1-theta[i]))^(nu[i]-y+37)
    if(runif(1) < alpha) {
      theta[i+1] = theta.star
      accept.theta = accept.theta + 1
    }
  }
  else theta[i+1] = theta[i]
  nu[i+1] = y + rpois(1, 200*(1-theta[i+1]))
}
par(mfrow=c(2,1))
plot(nu,pch="."); plot(theta,pch=".")
accept.theta/N

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

#try even larger
  theta.star = rnorm(1, theta[i], .07)
  theta.star = rnorm(1, theta[i], .15)
  theta.star = rnorm(1, theta[i], .4)
  theta.star = rnorm(1, theta[i], 1)
  theta.star = rnorm(1, theta[i], 20)

#try smaller
  theta.star = rnorm(1, theta[i], .03)
  theta.star = rnorm(1, theta[i], .01)
  theta.star = rnorm(1, theta[i], .001)