Sunday, February 8, 2015

R12. Normal Approximation to Binomial in R

According to the Central Limit Theorem, average of distributions such as Binomial can be approximated by the Normal for large n.


The function pnorm and dbinom are used. For the binomial, the probabilities are summed for a range near p*n.


For the normal, the pnorm at x = 0 is p = 0.5. Thus only the right side probability is used and multiplied by 2 since normal is symmetric about mean.

# ex12.R
p <- 0.4
n <- seq(100,2000,100)
s <- numeric(length(n))
nrm <- numeric(length(n))
print('Probability of successes n*(p-.01):n*(p+.01) in n: -> s')
for (j in 1:length(n)) {
  vals <- as.integer(c(n[j]*(p-.01), n[j]*(p+.01)))
  s[j] <- sum(dbinom(vals[1]:vals[2], size = n[j], prob = p ))
}
print("Normal approximation: -> nrm")
for (j in 1:length(n)) {
  mn <- n[j]*p
  std <- sqrt(n[j]*p*(1-p))
  vals <- as.integer(c(n[j]*(p-.01), n[j]*(p+.01)))
  nrm[j] <- 2*(pnorm(vals[2]+0.5, mean = mn, sd = std)-.5)
}
error <- (s-nrm)^2
plot(n, error, type='l', col = 'blue')
# [1] "Probability of successes n*(p-.01):n*(p+.01) in n: -> s"
# [1] "Normal approximation: -> nrm"

Output:

No comments:

Post a Comment