B.10 - JAGS model to infer $p$ (see Sec. )

#---- data and parameters
nr = 1000000
ns = 10000
nP = 2010
r0 = s0 = 1
r1 = 409.1; s1 = 9.1
r2 = 25.2;  s2 = 193.1 

#---- JAGS model ------------------------------
model = "tmp_model.bug"    # name of the model file ('temporary')
write("
model {
  nP ~ sum(nP.I, nP.NI)
  nP.I ~ dbin(pi1, n.I)
  nP.NI ~ dbin(pi2, n.NI)
  pi1 ~ dbeta(r1, s1)
  pi2 ~ dbeta(r2, s2)
  n.I ~ dbin(p, ns)
  n.NI <- ns - n.I
  p ~ dbeta(r0,s0)        
}
", model)

#---- call JAGS ---------------------------------------------
library(rjags)
data <- list(ns=ns, nP=nP, r0=r0, s0=s0, r1=r1, s1=s1, r2=r2, s2=s2)  
jm <- jags.model(model, data)
update(jm, 10000)
to.monitor <-  c('p', 'n.I')
chain <- coda.samples(jm, to.monitor, n.iter=nr)

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