B.14 - Approximated $f(p)$ from the first four moments of the distribution (see Sec. , Eq. ())

library("PearsonDS") # (package needs to be installed)
pause <- function() { cat ("\n >> press Enter to continue\n"); scan() }
nP = 201; ns = 1000
r0 = 1; s0 = 1
r1 = 409.1; s1 = 9.1
r2 =  25.1; s2 = 193.1
Nf = sum.p = sum.p2 = sum.p3 = sum.p4 = 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) )
     sum.p3 <- sum.p3 + exp( l0 +  lgamma(r0+nI+3) - lgamma(r0+s0+ns+3) )
     sum.p4 <- sum.p4 + exp( l0 +  lgamma(r0+nI+4) - lgamma(r0+s0+ns+4) )
  }
}
E.p  <- sum.p/Nf
E.p2 <- sum.p2/Nf
s2.p <- E.p2-E.p^2
s.p  <- sqrt(s2.p)
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))) 
E.p3 <- sum.p3/Nf
Skew <- ( E.p3 - 3*E.p2*E.p + 2*E.p^3  )  / s.p^3   
cat(sprintf("E(p^3) = %.3e; Skewness = %.3f \n",E.p3, Skew))
E.p4 <- sum.p4/Nf
Kurt <-  ( E.p4 - 4*E.p3*E.p + 6*E.p2*E.p^2 - 3*E.p^4 )   / s.p^4  
cat(sprintf("E(p^4) = %.3e; Kurtosis = %.3f\n", E.p4, Kurt))

moments <- c(mean = E.p,variance = s2.p,skewness = Skew, kurtosis = Kurt)
curve(dpearson(x, moments=moments), max(0.001,E.p-4*s.p), min(0.999,E.p+4*s.p),
      lwd=2, col='blue', xlab='p', ylab='f(p)')
pr <- rpearson(100000, moments = moments)
hist(pr , nc=100, freq=FALSE, col='cyan', add=TRUE)
cat(sprintf("\n\n mean    : %.4f\n", mean(pr))) 
cat(sprintf(" sigma: %.4f\n",sd(pr)))