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)))