B.11 - Inferring the proportions of infected in two different populations (see Sec. )

model = "tmp_model.bug" 
write("
model {
  n.I.1 ~ dbin(p1, ns1)
  p1 ~ dbeta(r0, s0)     
  n.NI.1 <- ns1 - n.I.1 
  nP.I.1 ~ dbin(pi1, n.I.1)
  nP.NI.1 ~ dbin(pi2, n.NI.1)
  nP.1 ~ sum(nP.I.1, nP.NI.1)

  n.I.2 ~ dbin(p2, ns2)
  p2 ~ dbeta(r0, r0)    
  n.NI.2 <- ns2 - n.I.2 
  nP.I.2 ~ dbin(pi1, n.I.2)
  nP.NI.2 ~ dbin(pi2, n.NI.2)
  nP.2 ~ sum(nP.I.2, nP.NI.2)

  Dp <- p2 - p1

  pi1 ~ dbeta(r1, s1)
  pi2 ~ dbeta(r2, s2)
  }
", model)

library(rjags)
set.seed(193)
nr = 1000000
ns1 = 10000
ns2 = 10000   # they could be different
nP.1 = 2000
nP.2 = 2200
r0 = s0 = 1   # flat prior
r1 = 409.1; s1 = 9.1
r2 = 25.2;  s2 = 193.1

data <- list(ns1=ns1, ns2=ns2, nP.1=nP.1, nP.2=nP.2,
             r0=r0, s0=s0, r1=r1, s1=s1, r2=r2, s2=s2)  
jm <- jags.model(model, data)
update(jm, 10000)
to.monitor <-  c('p1', 'p2', 'Dp')
chain <- coda.samples(jm, to.monitor, n.iter=nr)
print(summary(chain))