The Era of High-Performance Computing

In today’s world of advanced computing, various algorithms exist that, while not achieving absolute precision, can get remarkably close. For practical purposes, this level of accuracy is often sufficient. In R, we can derive our final probability using:

> dbinom(4500, size = 10000, prob = 0.5)

[1] 1.422495e-24

The dbinom() function, unlike our earlier formula, employs an ingenious algorithm created by Catherine Loader, a Stanford statistician, in 2000 during her tenure at Bell Labs. This method retains 15 pre-computed values for specific parameters and calculates the distribution based on these values, maintaining precision up to 18 digits.

Historical Context: The Good Old Days

Before the advent of high-performance computing, statisticians used tables for binomial probabilities with smaller n values. For larger n values (typically n > 20), they relied on the normal approximation of the binomial distribution to estimate probability densities. While distributions like the binomial and Poisson distributions are defined for discrete integer values, the normal distribution is continuous, permitting broader calculation options.

To visualize the discrete nature of the binomial distribution, let’s plot it as a continuous function:

library(ggplot2)

ggplot() +

xlim(40, 60) +

geom_function(fun = function(x) dbinom(x, 100, 0.5), color = "blue") +

labs(x = "Successes", y = "Density") +

theme_minimal()

Normal Approximation of the Binomial Distribution

Next, we can overlay a normal distribution with a mean of 50 and a standard deviation of 5 onto the previous plot:

ggplot() +

xlim(40, 60) +

geom_function(fun = function(x) dbinom(x, 100, 0.5), color = "blue") +

geom_function(fun = function(x) dnorm(x, 50, 5), color = "red") +

labs(x = "Successes", y = "Density") +

theme_minimal()

The resulting graph shows that the density from the normal distribution aligns closely with that of the binomial distribution for relevant integer values.

For a binomial distribution B(n, p), the expected number of successful events is np, and thus np serves as the mean. Each independent event (like a coin toss) has a variance of p(1-p), leading to a total variance of np(1-p) and a standard deviation of √(np(1-p)). Consequently, we can approximate B(n, p) by the normal distribution N(np, √(np(1-p))). The diagram above illustrates that B(100, 0.5) ~ N(50, 5).

This approximation is fairly accurate unless n is quite small or the event is heavily biased (i.e., p diverges significantly from 0.5). To explore this, we can create a function that compares binomial and normal densities:

test_norm <- function(size, prob, range) {

df <- data.frame(

n = range,

binom = dbinom(range, size = size, prob = prob),

norm = dnorm(range, mean = size * prob, sd = sqrt(size * prob * (1 - prob)))

)

ggplot(data = df, aes(x = n)) +

geom_line(y = df$binom, color = "red") +

geom_line(y = df$norm, color = "blue") +

theme_minimal()

}

Let’s evaluate this for 5, 10, 25, and 50 coin tosses using a fair coin:

library(patchwork)

for (n in c(5, 10, 50, 100)) {

assign(paste0("p", n), test_norm(n, 0.5, (n/2 - 0.3*n):(n/2 + 0.3*n)))

}

p5 + p10 + p25 + p50

Even with just 10 tosses, the red and blue lines nearly overlap. Now, let’s analyze biased coins with 100 tosses:

for (p in c(0.1, 0.2, 0.3, 0.4)) {

assign(paste0("p", p), test_norm(100, p, (100*p - 0.6*100*p):(100*p + 0.6*100*p)))

}

p0.1 + p0.2 + p0.3 + p0.4

Even for coins biased at p = 0.2, distinguishing the normal approximation from the binomial distribution remains challenging.

Practical Applications of Statistical Concepts

A notable question posed during a 1991 Cambridge University mathematics entrance exam was:

"A fair coin is flipped 10,000 times. Each head scores 1 point, while each tail deducts 1 point. Estimate the probability of achieving a final score greater than 100."

To solve this, we need at least 100 more heads than tails across the 10,000 flips, meaning the number of heads must exceed 5,050. Thus, we seek P(x > 5050) in the binomial distribution B(10000, 0.5).

We can approximate B(10000, 0.5) as N(5000, 50). Therefore, we need to determine the area under the normal curve of N(5000, 50) for x ≥ 5051. Applying continuity correction, we seek P(x > 5050.5):

We calculate the z-score, which is +1.01, as we are looking for the area under the curve to the right of 50.5/50 = 1.01 standard deviations above the mean. Using tables, scientific calculators, or R functions, we find the upper-tail p-value for this z-score:

> pnorm(1.01, lower.tail = FALSE)

[1] 0.1562476

Let’s verify this with the R function for calculating the p-value for a binomial distribution:

> pbinom(5050, 10000, 0.5, lower.tail = FALSE)

[1] 0.1562476

The results perfectly align, revealing that the probability of scoring over 100 is approximately 15.62%.

For those interested in an alternative perspective using Python, Christian Janiake has provided a gist here.

What are your thoughts on this exploration of statistical fundamentals? Feel free to share your comments!