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