B.13 - Exact calculation of E$(p)$ and $\sigma(p)$ using a Beta prior (see Sec. , footnote )

nP = 201; ns = 1000
r0 = 1; s0 = 1
r1 = 409.1; s1 = 9.1
r2 =  25.1; s2 = 193.1

Nf   <- 0
sum.p  <- 0
sum.p2 <- 0
for (nPI in 0:nP) {
    for (nI in 0:ns) {
     l0 <-  ( lchoose(nI, nPI) + lchoose(ns-nI,nP-nPI) + lchoose(ns,nI)
          + lgamma(nPI+r1) + lgamma(nI-nPI+s1) - lgamma(r1+nI+s1)
          + lgamma (nP-nPI+r2) + lgamma(ns-nI-nP+nPI+s2) - lgamma(ns-nI+s2+r2) 
          + lgamma(s0+ns-nI) )  
     Nf     <- Nf     + exp( l0 +  lgamma(r0+nI)   - lgamma(r0+s0+ns) )
     sum.p  <- sum.p  + exp( l0 +  lgamma(r0+nI+1) - lgamma(r0+s0+ns+1) )
     sum.p2 <- sum.p2 + exp( l0 +  lgamma(r0+nI+2) - lgamma(r0+s0+ns+2) )
  }
}
E.p  <- sum.p/Nf
E.p2 <- sum.p2/Nf
cat(sprintf("nP=%d, ns=%d\n", nP, ns))
cat(sprintf("E(p)    : %.4f\n",E.p)) 
cat(sprintf("sigma(p): %.4f\n",sqrt(E.p2-E.p^2)))