B.1 - Monte Carlo evaluation of Eqs. () and ()

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))))