B.7 - Check of approximated formulae

get.prediction <- function(ns, N, p, E.pi1, sigma.pi1, E.pi2, sigma.pi2) {
   E.ps <- p
   E.fP <- E.pi1*E.ps + E.pi2*(1-E.ps)

   s.fP.R   <- sqrt( E.pi1*(1-E.pi1)*E.ps +  E.pi2*(1-E.pi2)*(1-E.ps) )/sqrt(ns)  
   s.fP.pi1 <- sigma.pi1*E.ps
   s.fP.pi2 <- sigma.pi2*(1-E.ps)
   s.ps     <- sqrt(p*(1-p)*(1-ns/N))/sqrt(ns)  
   s.fP.ps  <- s.ps*abs(E.pi1-E.pi2)           

   s.fP.stat <- sqrt(s.fP.R^2+s.fP.ps^2)
   s.fP.syst <- sqrt(s.fP.pi1^2+s.fP.pi2^2)
   s.fP = sqrt(s.fP.stat^2 + s.fP.syst^2)
   return(list(E.fP=E.fP, s.fP=s.fP))
} 

N  = 10^7  
p.v  <- 0:4 / 10
ns.v <- c(300, 1000, 3000, 10000)

E.pi1 = 0.978;  sigma.pi1 = 0.007
E.pi2 = 0.115;  sigma.pi2 = 0.022
for(case in 1:3) {
  if(case == 2) sigma.pi2 = sigma.pi1
  if(case == 3) E.pi2 = 1 - E.pi1
  cat(sprintf("\n [pi1 = %.3f+-%.3f;", E.pi1, sigma.pi1))
  cat(sprintf(" pi2 = %.3f+-%.3f]\n", E.pi2, sigma.pi2))
  for(i in 1:length(ns.v)) {
    ns <- ns.v[i]
    cat(sprintf(" ns = %d\n p: ", ns))
    for(j in 1:length(p.v))  cat(sprintf("   %.2f      ", p.v[j]))
    cat("\nfP: ")
    for(j in 1:length(p.v)) {
     p <- p.v[j]   
      pred <- get.prediction(ns, N, p, E.pi1, sigma.pi1, E.pi2, sigma.pi2)
     cat(sprintf("%.3f+-%.3f ", pred$E.fP, pred$s.fP))
    }
    cat("\n")
  }
}