B.5 - Monte Carlo estimate of $f_P$ using R functions, as described in Sec.  (see Fig.)

p  = 0.1
ns = 1000
r1 = 409.1;  s1 =  9.1
r2 =  25.1;  s2 =193.2

nr = 10000
n.I   <- rbinom(nr, ns, p)       # 1. 
n.NI  <- ns - n.I
pi1   <- rbeta(nr, r1, s1)       # 2.
pi2   <- rbeta(nr, r2, s2)
nP.I  <- rbinom(nr, n.I, pi1)    # 3.
nP.NI <- rbinom(nr, n.NI, pi2)
nP    <- nP.I + nP.NI            # 4.

fP <- nP/ns
cat(sprintf("fP: %.3f +- %.3f\n", mean(fP), sd(fP)))
hist(fP, col='cyan')
# barplot(table(fP), col='cyan')  # alternative, for small ns and p