#------------ 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 <- rho * r2 r2 ~ dgamma(1, 1e-6) rho ~ dgamma(1, 1e-6) } ", 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/T2, sqrt(x2)/T2) ) mu.rho <- ((x1+1)/T1)/((x2-1)/T2) sigma.rho <- sqrt(mu.rho*(T2/T1*(x1+2)/(x2-2)-mu.rho)) cat(sprintf(" rho = %.3f +- %.3f\n", mu.rho, sigma.rho))