get.fP  <- function(nr, p, ns, pi1, pi2) {
  n.I   <- rbinom(nr, ns, p)
  n.NI  <- ns - n.I
  nP.I  <- rbinom(nr, n.I, pi1)    
  nP.NI <- rbinom(nr, n.NI, pi2)
  nP    <- nP.I + nP.NI            
  fP    <- nP/ns
  return(fP)
}

p1 = 0.1
p2 = 0.2
ns = 10000
r1 = 409.1;  s1 =  9.1
r2 =  25.1;  s2 =193.2

nr = 100000
pi1 <- rbeta(nr, r1, s1)   
pi2 <- rbeta(nr, r2, s2)
fP1 <- get.fP(nr, p1, ns, pi1, pi2)
fP2 <- get.fP(nr, p2, ns, pi1, pi2)

cat(sprintf("fP1: %.4f +- %.4f\n", mean(fP1), sd(fP1)))
cat(sprintf("fP2: %.4f +- %.4f\n", mean(fP2), sd(fP2)))
cat(sprintf("fP2-fP1: %.4f +- %.4f\n", mean(fP2-fP1), sd(fP2-fP1)))
cat(sprintf("rho(fP1,fP2): %.4f\n", cor(fP1,fP2)))
cat(sprintf("Check of sigma(fP2-fP1) from correlation: %.4f\n",
            sqrt(var(fP1)+var(fP2)-2*cov(fP1,fP2))))