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))))(Note how the `random' sequences of values of pi1 and pi2 are generated before the calls to get.fP(). This is crucial in order to get the correlation among fP1 and fP2 discussed in Sec. . If these two sequences are defined, each time, inside the function, or they a generated before each call to the function, the correlation will disappear. This way of generating the events is consequence of our model assumptions, as stressed in Sec. .)