p <- 1- pnorm(4) cat("Parameter to be estimated",p,"\n") n <- 10000000 x <- rnorm(n) k <- length(x[x>4]) cat("Count",k,"naive estimate",k/n,"error",p-k/n,"; s.d.",sqrt(p*(1-p)/n),"\n") kk <- length(x[abs(x)>4]) cat("Count",kk,"improved estimate",kk/(2*n),"error",p-kk/(2*n), "; s.d.",sqrt(2*p*(1-2*p)/(4*n)),"\n")