B.6 - Monte Carlo estimate of $f_P$ JAGS from R via rjags (Sec. )

#--------------------------  JAGS model ------------------------------
model.name = "tmp_model.bug"  
write("
model {
  n.I ~ dbin(p, ns)
  n.NI <- ns - n.I 
  nP.I ~ dbin(pi1, n.I)
  nP.NI ~ dbin(pi2, n.NI)
  pi1 ~ dbeta(r1, s1)
  pi2 ~ dbeta(r2, s2)
  nP ~ sum(nP.I, nP.NI)
  fP <- nP / ns 
}
", model)

#--------------------------- call JAGS ---------------------------------
library(rjags)
data <- list(ns=1000, p=0.1, r1=409.1, s1=9.1, r2=25.2, s2=193.1)  
jm <- jags.model(model, data)
chain <- coda.samples(jm, c('n.I', 'fP'), n.iter=10000)

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