B.5 - Example of JAGS inference of rates and their ratio for the model of Fig. 

#------------  Data  --------------------------------------------
x1 = 3; T1=3
x2 = 6; T2=6
nr = 1e5

#-------------  JAGS model --------------------------------------
library(rjags)
model = "tmp_model.bug"    # name of the model file ('temporary')
write("
model {
  x1 ~ dpois(lambda1)
  x2 ~ dpois(lambda2)
  lambda1 <- r1 * T1
  lambda2 <- r2 * T2
  r1 ~ dgamma(1, 1e-6) 
  r2 ~ dgamma(1, 1e-6)
  rho <- r1/r2
}
", model)

#------------  JAGS call via rjags ------------------------------
data <- list(x1=x1, T1=T1, x2=x2, T2=T2)
jm <- jags.model(model, data)
update(jm, 100)
to.monitor <-  c('r1', 'r2', 'rho')
chain <- coda.samples(jm, to.monitor, n.iter=nr)

#------------  Results -----------------------------------------
print(summary(chain))
plot(chain, col='blue')

cat(sprintf("Exact: \n"))
cat(sprintf("   r1 =  %.3f +- %.3f\n", (x1+1)/T1, sqrt(x1+1)/T1))
cat(sprintf("   r2 =  %.3f +- %.3f\n", (x2+1)/T2, sqrt(x2+1)/T2))
mu.rho    <-  ((x1+1)/T1)/(x2/T2)
sigma.rho <- sqrt(mu.rho*(T2/T1*(x1+2)/(x2-1)-mu.rho)) 
cat(sprintf("  rho =  %.3f +- %.3f\n", mu.rho, sigma.rho))