Binomial & Poisson Distributions

Many questions in statistics are about counts. How many questions will a student get right by pure guessing? How many sixes in an over? How many orders will a bakery receive in the next 10 minutes? How many parcels in a batch of 1,000 will arrive damaged? Two distributions answer most such questions: the binomial and the Poisson.

In this lesson you will learn when each one applies, compute their probabilities by hand from the formulas, check them with scipy.stats, and confirm them by simulation.

The binomial setting

A count X follows a binomial distribution when all four of these conditions hold:

1. There is a fixed number of trials, n.
2. Each trial has two outcomes, called "success" and "failure".
3. The probability of success, p, is the same on every trial.
4. The trials are independent.

Guessing on a multiple-choice test fits well: 10 questions (n = 10), each answered right or wrong, each with a 1 in 4 chance (p = 0.25) if you pick randomly, and one guess does not affect another.

The probability of exactly k successes is:

P(X = k) = C(n, k) x p^k x (1 - p)^(n - k)

mean     = n x p
variance = n x p x (1 - p)

Why? Any one particular order of k successes and n - k failures has probability p^k x (1 - p)^(n - k), by the multiplication rule for independent events. There are C(n, k) such orders (the number of ways to choose which k trials are the successes), so you multiply by C(n, k).

Binomial probabilities by hand and with SciPy

What is the chance of getting exactly 5 out of 10 right by guessing?

import math
from scipy import stats

# Guessing all 10 questions of an MCQ quiz, 4 options each
n, p = 10, 0.25
k = 5
by_hand = math.comb(n, k) * p**k * (1 - p)**(n - k)
print(f"P(exactly {k} right) by hand: {by_hand:.4f}")
print(f"stats.binom.pmf:             {stats.binom.pmf(k, n, p):.4f}")
print("Mean n*p =", n * p, "  SD =", round(math.sqrt(n * p * (1 - p)), 3))

P(exactly 5 right) by hand: 0.0584
stats.binom.pmf:             0.0584
Mean n*p = 2.5   SD = 1.369

By hand: C(10, 5) = 252, and 252 x 0.25^5 x 0.75^5 = 0.0584. A guesser expects 2.5 right answers on average. In SciPy, pmf stands for probability mass function: the probability of one exact value.

For "at most" and "at least" questions, use the cumulative distribution function cdf(k) = P(X ≤ k) and the survival function sf(k) = P(X > k) = 1 - cdf(k). Be careful with the boundary: "at least 5" means X ≥ 5, which is X > 4, so you need sf(4).

from scipy import stats

quiz = stats.binom(10, 0.25)
for k in range(11):
    print(f"{k:>2} right: {quiz.pmf(k):.4f}  {'#' * round(quiz.pmf(k) * 100)}")
print("P(at most 2)  =", round(quiz.cdf(2), 4))
print("P(at least 5) =", round(quiz.sf(4), 4))

 0 right: 0.0563  ######
 1 right: 0.1877  ###################
 2 right: 0.2816  ############################
 3 right: 0.2503  #########################
 4 right: 0.1460  ###############
 5 right: 0.0584  ######
 6 right: 0.0162  ##
 7 right: 0.0031  
 8 right: 0.0004  
 9 right: 0.0000  
10 right: 0.0000  
P(at most 2)  = 0.5256
P(at least 5) = 0.0781

stats.binom(10, 0.25) creates a "frozen" distribution, so you do not have to repeat n and p in every call. The whole table adds up to 1. A guesser has under an 8% chance of scoring 5 or more, and getting all 10 right has a probability of about 1 in a million (0.25^10). This is how a teacher can tell that a student who scores 9 out of 10 was very probably not guessing.

Checking with a simulation

rng.binomial(n, p, size) simulates many binomial counts at once. Here, 100,000 students guess their way through the test.

import numpy as np

rng = np.random.default_rng(4)
right = rng.binomial(n=10, p=0.25, size=100_000)
print("Simulated P(at least 5):", (right >= 5).mean())
print("Simulated mean:", right.mean(), " variance:", round(right.var(), 3))

Simulated P(at least 5): 0.07769
Simulated mean: 2.49806  variance: 1.866

The simulation agrees with the formulas: P(at least 5) close to 0.0781, mean close to 2.5 and variance close to n x p x (1 - p) = 1.875.

The Poisson distribution

Sometimes there is no fixed number of trials. You count how many times something happens in a fixed stretch of time or space: orders per 10 minutes, calls to a helpline per hour, potholes per kilometre of road. If events happen independently and at a steady average rate λ (lambda), the count follows a Poisson distribution:

P(X = k) = e^(-λ) x λ^k / k!        k = 0, 1, 2, ...

mean = λ        variance = λ

The fact that the mean and the variance are equal is a quick check: if your counts have a variance much larger than their mean, the Poisson model probably does not fit.

import math
from scipy import stats

# A Pune bakery gets 3 online orders per 10 minutes on average (invented)
lam = 3
for k in range(8):
    by_hand = math.exp(-lam) * lam**k / math.factorial(k)
    print(f"P({k} orders) = {by_hand:.4f}   scipy {stats.poisson.pmf(k, lam):.4f}")
print("P(more than 5) =", round(stats.poisson.sf(5, lam), 4))
print("Mean and variance:", stats.poisson.mean(lam), stats.poisson.var(lam))

P(0 orders) = 0.0498   scipy 0.0498
P(1 orders) = 0.1494   scipy 0.1494
P(2 orders) = 0.2240   scipy 0.2240
P(3 orders) = 0.2240   scipy 0.2240
P(4 orders) = 0.1680   scipy 0.1680
P(5 orders) = 0.1008   scipy 0.1008
P(6 orders) = 0.0504   scipy 0.0504
P(7 orders) = 0.0216   scipy 0.0216
P(more than 5) = 0.0839
Mean and variance: 3.0 3.0

The bakery gets no orders in a 10-minute window only about 5% of the time, and more than 5 orders about 8% of the time. If the owner plans for up to 5 orders per 10 minutes, the kitchen will be overloaded in roughly one window in twelve. The rate must match the window: 3 orders per 10 minutes is λ = 18 per hour.

Poisson as an approximation to the binomial

When n is large and p is small, the binomial is very close to a Poisson with λ = n x p. Before computers this saved a lot of work, and it still explains why rare events in large groups (defects, accidents, typing errors per page) often follow a Poisson pattern.

from scipy import stats

# 1,000 parcels, each damaged with probability 0.002 (invented)
n, p = 1000, 0.002
for k in range(6):
    print(f"k = {k}: binomial {stats.binom.pmf(k, n, p):.4f}   Poisson {stats.poisson.pmf(k, n * p):.4f}")

k = 0: binomial 0.1351   Poisson 0.1353
k = 1: binomial 0.2707   Poisson 0.2707
k = 2: binomial 0.2709   Poisson 0.2707
k = 3: binomial 0.1806   Poisson 0.1804
k = 4: binomial 0.0902   Poisson 0.0902
k = 5: binomial 0.0360   Poisson 0.0361

With λ = 1000 x 0.002 = 2, the two columns agree to about three decimal places.

Seeing the shapes

Plot both distributions as bar charts. Counts are discrete, so bars (not a smooth curve) are the honest picture.

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

k = np.arange(0, 11)
fig, axes = plt.subplots(1, 2, figsize=(9, 3.2), sharey=True)
axes[0].bar(k, stats.binom.pmf(k, 10, 0.25), color="#7c3aed")
axes[0].set_title("Binomial: n = 10, p = 0.25")
axes[1].bar(k, stats.poisson.pmf(k, 3), color="#0ea5e9")
axes[1].set_title("Poisson: lambda = 3")
for ax in axes:
    ax.set_xlabel("k")
axes[0].set_ylabel("P(X = k)")
plt.show()

Both are right-skewed when the mean is small, because counts cannot go below zero but can stretch upwards. As n x p (or λ) grows, both shapes become more symmetric and bell-shaped, which leads to the normal distribution.

Common mistakes

• Using the binomial when trials are not independent. Drawing cards without replacement changes p each time; for small populations that breaks the binomial.
• Off-by-one errors with sf. P(X ≥ 5) is sf(4), not sf(5).
• Mismatched rates. λ must be for the same window you ask about.
• Forcing a Poisson model onto clustered events. Wedding-season orders or festival traffic come in bursts; the variance will be much larger than the mean.

Key points

• Binomial: fixed n independent trials with the same success probability p; P(X = k) = C(n, k) p^k (1 - p)^(n - k).
• Binomial mean = np and variance = np(1 - p).
• Poisson: counts of independent events at an average rate λ; P(X = k) = e^(-λ) λ^k / k!, with mean = variance = λ.
• In SciPy: pmf for exactly k, cdf(k) for at most k, sf(k) for more than k.
• For large n and small p, binomial(n, p) is close to Poisson(λ = np).

import math
from scipy import stats

n, p = 30, 0.08
print("Expected sixes:", round(n * p, 2))
print("P(no sixes) by hand:", round(math.comb(n, 0) * p**0 * (1 - p)**n, 4))

sixes = stats.binom(n, p)
print("P(no sixes):        ", round(sixes.pmf(0), 4))
print("P(at least 4):      ", round(sixes.sf(3), 4))
print("Poisson, at least 4:", round(stats.poisson.sf(3, n * p), 4))

Expected sixes: 2.4
P(no sixes) by hand: 0.082
P(no sixes):         0.082
P(at least 4):       0.2158
Poisson, at least 4: 0.2213