r.pi1 = 409.1; s.pi1 = 9.1 r.pi2 = 25.1; s.pi2 = 193.1 p = 0.1 n = 100000 pi1 <- rbeta(n, r.pi1, s.pi1) pi2 <- rbeta(n, r.pi2, s.pi2) P.Inf.Pos.i <- pi1*p/(pi1*p + pi2*(1-p)) P.NoInf.Neg.i <- (1-pi2)*(1-p) / ((1-pi1)*p + (1-pi2)*(1-p)) P.Inf.Pos <- mean(P.Inf.Pos.i) P.NoInf.Neg <- mean(P.NoInf.Neg.i) cat(sprintf("Integral (by MC): P(Inf|Pos) = %.4f; P(NoInf|Neg) = %.4f \n", P.Inf.Pos, P.NoInf.Neg)) E.pi1 = r.pi1 / (r.pi1 + s.pi1) E.pi2 = r.pi2 / (r.pi2 + s.pi2) cat(sprintf("Using E(..): P(Inf|Pos) = %.4f; P(NoInf|Neg) = %.4f \n", E.pi1*p/(E.pi1*p + E.pi2*(1-p)), (1-E.pi2)*(1-p) / ((1-E.pi1)*p + (1-E.pi2)*(1-p))))