B.2 - Monte Carlo evaluation of Eq. ()

ns = 10000
ps = 0.1
n.I = ps * ns
n.NI = (1 - ps) * ns
r1=409.1; s1=9.1
r2=25.1; s2=193.2
nr =100000
pi1 <- rbeta(nr, r1, s1)
pi2 <- rbeta(nr, r2, s2)
nP.I  <- rbinom(nr, n.I, pi1)
nP.NI <- rbinom(nr, n.NI, pi2)
nP <- nP.I + nP.NI
hist(nP, nc=100, col='cyan', freq=FALSE) 
cat(sprintf("nP: mean = %.1f, sd = %.1f\n",mean(nP),sd(nP)))