B.3 - Monte Carlo estimate of $f(n_P\,\vert\,N=10^5,n_s=10^4, p=0.1)$ reported in Fig. 

N = 100000; ns  = 10000; p = 0.1
N.I  <- as.integer(p*N)
N.NI <- N - N.I
E.pi1 = 0.978; E.pi2 = 0.115

nr =100000
n.I  <- rhyper(nr, m=N.I, n=N.NI, k=ns)
n.NI <- ns - n.I

nP.I  <- rbinom(nr, n.I, E.pi1)
nP.NI <- rbinom(nr, n.NI, E.pi2)
nP <- nP.I + nP.NI
cat(sprintf("mean = %.0f,  sigma = %.0f\n", mean(nP),  sd(nP) ))
hist(nP, nc=90, col='cyan', freq=FALSE)