02 - Discrete Random Variables: Bernoulli, Binomial, Poisson

You take a quiz of 20 multiple-choice questions, each with four options, and answer by guessing. How many correct answers do you expect? And with what probability do you hit at least ten, enough to pass? These are questions about the count of successes in repeated trials, and their answer has a name: the binomial distribution.

Here you derive it, compute it exactly (the mean is 5, the variance 15/4, the probability of reaching ten is about 1.4%) and confirm it with a million simulated quizzes. Then you see how the binomial becomes, when success is rare, a Poisson distribution, how large the error you make is, and how the geometric distribution describes the wait for the first success. You need to have read the first article of the series on axioms, independence and conditional probability; the rest I define here.

Giving a number to every outcome

A random variable is a rule that associates a number with every outcome of the experiment: in the quiz, the number of correct answers. It is called discrete if the possible values can be listed (0, 1, 2, …). Its complete description is the probability mass function (PMF), p(k) = P(X = k): a non-negative number for each value, with sum 1. Two summaries can be drawn from the PMF.

The expected value E[X] = Σ k p(k) is the mean of the values weighted by the probabilities, the "centre of mass" of the distribution. The variance Var(X) = E[(X − E[X])²] = E[X²] − (E[X])² measures how much the values scatter around the expected value; its square root is the standard deviation, in the same unit of measurement as X.

Two properties do almost all the work in the calculations that follow. The first is the linearity of expectation: E[X + Y] = E[X] + E[Y] always, even if X and Y depend on each other.

The second concerns the variance, and here caution is needed: Var(X + Y) = Var(X) + Var(Y) holds when X and Y are independent, not in general (in the article on covariance you see the term that is missing when they are not). The difference between the two rules is the reason why the mean of a sum can always be computed, while its variance requires knowing how the pieces relate to each other.

Bernoulli and binomial: counting successes

The simplest variable is the Bernoulli with parameter p: it is 1 ("success") with probability p and 0 with probability q = 1 − p. Its expected value is p and its variance pq, because X² = X and therefore E[X²] = p. A quiz question, answered at random, is a Bernoulli with p = 1/4. If you add up n independent Bernoullis with the same p you obtain the number of successes Y, which has a binomial distribution Bin(n, p).

The PMF is derived by counting: a particular sequence with k successes and n − k failures has probability p^kq^n−k by independence, and the sequences with k successes number C(n, k) = n! / (k!(n − k)!), whence P(Y = k) = C(n, k) p^kq^n−k.

The sum of n Bernoullis gives the moments at once, without computing series: by linearity E[Y] = np, and by independence Var(Y) = npq.

The program checks these formulas in the most direct way: it builds the PMF with exact fractions (for the quiz, n = 20 and p = 1/4), checks that it sums to 1, computes the mean and variance from the definition, and compares the probability of reaching the threshold with scipy.stats.binom and with a million simulated quizzes. At the end it also looks for the threshold beyond which guessing costs less than one chance in a thousand.

from fractions import Fraction as F
from math import comb
import numpy as np
from scipy.stats import binom

# A quiz of 20 multiple-choice questions with 4 options each, answered at random: X = number of correct answers
n, p = 20, F(1, 4)
q = 1 - p
pmf = [comb(n, k) * p**k * q**(n - k) for k in range(n + 1)]     # exact fractions
assert sum(pmf) == 1

media = sum(k * pk for k, pk in enumerate(pmf))
secondo = sum(k * k * pk for k, pk in enumerate(pmf))
var = secondo - media**2
print("E[X] =", media, "(np =", n * p, ")  Var(X) =", var, "(npq =", n * p * q, ")")
print("mode:", max(range(n + 1), key=lambda k: pmf[k]), " P(X=5) =", f"{float(pmf[5]):.4f}")

soglia = 10                                  # 'pass' if at least 10 out of 20
esatto = sum(pmf[soglia:])
print(f"P(X >= {soglia}) exact = {float(esatto):.6f}")
print(f"P(X >= {soglia}) scipy = {binom.sf(soglia - 1, n, 0.25):.6f}")

N = 10**6
rng = np.random.default_rng(11)
campione = rng.binomial(n, 0.25, size=N)     # N students guessing at random
freq = np.mean(campione >= soglia)
pe = float(esatto)
tol = 3 * (pe * (1 - pe) / N) ** 0.5
print(f"Monte Carlo N={N}: {freq:.6f}; tolerance 3 sigma {tol:.6f}; deviation {abs(freq - pe):.6f}")
print(f"sample mean {campione.mean():.4f}, sample variance {campione.var():.4f}")
assert abs(freq - pe) < tol

# Threshold below which guessing "by chance" has probability under 1 in 1000
k_min = next(k for k in range(n + 1) if sum(pmf[k:]) < F(1, 1000))
print(f"threshold: P(X >= {k_min}) = {float(sum(pmf[k_min:])):.5f} < 0.001, while P(X >= {k_min - 1}) = {float(sum(pmf[k_min - 1:])):.5f}")

# Output:
# E[X] = 5 (np = 5 )  Var(X) = 15/4 (npq = 15/4 )
# mode: 5  P(X=5) = 0.2023
# P(X >= 10) exact = 0.013864
# P(X >= 10) scipy = 0.013864
# Monte Carlo N=1000000: 0.013781; tolerance 3 sigma 0.000351; deviation 0.000083
# sample mean 5.0024, sample variance 3.7576
# threshold: P(X >= 12) = 0.00094 < 0.001, while P(X >= 11) = 0.00394

The mean 5 and the variance 15/4 come out exact from the definition and coincide with np and npq, so the proof by a sum of Bernoullis and the direct calculation say the same thing. The probability of getting at least 10 correct answers is 0.013864 in both calculations, exact and SciPy's, and the simulation gives 0.013781, with a deviation of 0.000083 against a tolerance of 0.000351 (three standard errors).

The practical result is clear-cut: guessing yields a pass in fewer than one case in seventy, and with 12 correct answers the probability drops below one in a thousand (0.00094). Note also the mode, the most probable value, which is 5 and is worth only 20%: the distribution is much wider than a peak, the standard deviation is about 1.94. The sample variance 3.758 is close to 15/4 = 3.75, as expected.

Poisson: the binomial when success is rare

When n is large and p is small, computing C(n, k) becomes awkward and, for many phenomena, not even n and p separately are known: all you know is the mean number of events, λ = np. The Poisson distribution with parameter λ has PMF P(Z = k) = e^−λ λ^k / k! for k = 0, 1, 2, …, and has mean and variance both equal to λ.

It is the limit of Bin(n, λ/n) as n tends to infinity: for fixed k, C(n, k) (λ/n)^k tends to λ^k/k!, and (1 − λ/n)^n−k tends to e^−λ. The equality between mean and variance is also a diagnostic test: count data with a variance much larger than the mean are not Poisson, and a different model is needed.

"Tends to" does not say how close the Poisson is for a concrete n. The natural measure is the total variation distance between two distributions, half the sum of the absolute differences of the PMFs: it is the largest difference the two laws can give to the probability of an event. For a sum of independent Bernoullis with probabilities p_i, Le Cam's inequality (1960) bounds it above by Σ p_i², which for n equal trials with p = λ/n is λ²/n.

The program computes the exact distance with λ = 3 for four values of n and checks the bound; then it applies the Poisson to a hypothetical use case, built for the example: an airline sells 105 tickets for 100 seats, and each passenger shows up with probability 0.95 independently of the others.

import numpy as np
from scipy.stats import binom, poisson

lam = 3.0
print(" n      TV(Bin(n,lam/n), Poisson(lam))   bound lam^2/n        ratio")
for n in (10, 100, 1000, 10000):
    k = np.arange(0, 200)
    tv = 0.5 * np.abs(binom.pmf(k, n, lam / n) - poisson.pmf(k, lam)).sum()
    bound = lam**2 / n                               # Le Cam inequality: sum p_i^2 = n (lam/n)^2
    assert tv <= bound
    print(f"{n:6d}  {tv:.6f}                       {bound:.6f}            {tv / bound:.3f}")

# A (hypothetical) use case: 105 tickets sold for 100 seats, each passenger shows up with prob. 0.95
vendute, posti, p_assente = 105, 100, 0.05
assenti = binom(vendute, p_assente)
esatto = assenti.cdf(4)                              # overbooking <=> no-shows <= 4  (shown up >= 101)
approx = poisson.cdf(4, vendute * p_assente)         # Poisson with lam = 5.25 on the number of NO-SHOWS
print(f"P(overbooking) exact = {esatto:.5f}  Poisson(5.25) = {approx:.5f}  difference {abs(esatto - approx):.5f}")
print(f"Le Cam bound on the distance: {vendute * p_assente**2:.5f}")

rng = np.random.default_rng(5)
N = 10**6
presenti = rng.binomial(vendute, 1 - p_assente, size=N)
freq = np.mean(presenti > posti)
tol = 3 * (esatto * (1 - esatto) / N) ** 0.5
print(f"Monte Carlo N={N}: {freq:.5f}; tolerance {tol:.5f}; deviation {abs(freq - esatto):.5f}")
assert abs(freq - esatto) < tol

# Output:
#  n      TV(Bin(n,lam/n), Poisson(lam))   bound lam^2/n        ratio
#     10  0.086409                       0.900000            0.096
#    100  0.007607                       0.090000            0.085
#   1000  0.000752                       0.009000            0.084
#  10000  0.000075                       0.000900            0.083
# P(overbooking) exact = 0.39243  Poisson(5.25) = 0.39777  difference 0.00534
# Le Cam bound on the distance: 0.26250
# Monte Carlo N=1000000: 0.39187; tolerance 0.00146; deviation 0.00057

The exact distance falls like 1/n, from 0.086 with n = 10 to 0.000075 with n = 10,000, and is always below the bound: in the four cases measured, Le Cam's bound is between ten and twelve times wider than the true value (ratio 0.083-0.096). So it serves to guarantee a maximum error, not to estimate it.

In the tickets case the Poisson is applied to the number of no-shows, not to that of the passengers present: the no-shows are the rare event, with mean 105 · 0.05 = 5.25, and there is overbooking when the no-shows are at most 4. The exact probability is 0.39243, the Poisson(5.25) gives 0.39777, and the simulation 0.39187. The approximation is off by 0.005 while Le Cam's bound, 0.2625, is fifty times wider than the deviation observed on the requested probability: the guarantee is prudent, not precise. The number is hypothetical and illustrates the method, it does not describe a real airline.

Geometric: waiting for the first success

If you repeat independent trials until you get the first success, the number X of trials has a geometric distribution: P(X = k) = q^k−1p, with mean 1/p and variance q/p². Here "attempts" counts the last one too, the success; other texts count only the failures and start from zero, so the mean becomes q/p: always check which convention the text or the library uses.

The characteristic property is memorylessness: P(X > m + n | X > m) = P(X > n), because the probability of exceeding m trials without success is q^m and the ratio q^m+n/q^m is qⁿ. The example is the die: you roll until a 6 comes up, so p = 1/6. The program checks the equality with fractions for three pairs (m, n), sums the series for the mean, and simulates a million sequences.

from fractions import Fraction as F
import numpy as np

# I roll a die until a 6 comes up: X = number of rolls (X = 1, 2, 3, ...)
p = F(1, 6); q = 1 - p
coda = lambda m: q**m                                # P(X > m): the first m rolls fail
for m, n in ((2, 3), (5, 4), (10, 1)):
    cond = coda(m + n) / coda(m)                     # P(X > m+n | X > m)
    print(f"m={m:2d} n={n}: P(X>m+n | X>m) = {cond}  =  P(X>n) = {coda(n)}  -> {cond == coda(n)}")

media = 1 / p; var = q / p**2
somma = sum(k * q**(k - 1) * p for k in range(1, 400))   # partial sum of the series for E[X]
print(f"E[X] = 1/p = {media}, Var(X) = q/p^2 = {var}; series truncated at 399 terms: {float(somma):.6f}")

rng = np.random.default_rng(3)
N = 10**6
x = rng.geometric(1/6, size=N)                       # numpy counts the attempts, success included
sigma_media = float(var) ** 0.5 / N ** 0.5
print(f"Monte Carlo: mean {x.mean():.4f} (expected {float(media):.0f} +/- {3*sigma_media:.4f}), variance {x.var():.2f} (expected {float(var):.0f})")
assert abs(x.mean() - float(media)) < 3 * sigma_media
# Conditional: among the runs with X > 5, the distribution of X-5 is again geometric
resto = x[x > 5] - 5
print(f"after 5 failures: {len(resto)} runs, mean number of rolls still needed {resto.mean():.4f}")

# Output:
# m= 2 n=3: P(X>m+n | X>m) = 125/216  =  P(X>n) = 125/216  -> True
# m= 5 n=4: P(X>m+n | X>m) = 625/1296  =  P(X>n) = 625/1296  -> True
# m=10 n=1: P(X>m+n | X>m) = 5/6  =  P(X>n) = 5/6  -> True
# E[X] = 1/p = 6, Var(X) = q/p^2 = 30; series truncated at 399 terms: 6.000000
# Monte Carlo: mean 6.0037 (expected 6 +/- 0.0164), variance 29.97 (expected 30)
# after 5 failures: 402794 runs, mean number of rolls still needed 5.9900

The three equalities are exact (for example 625/1296 on both sides), the series truncated at 399 terms gives 6.000000, and the simulation finds mean 6.0037 against an expected 6 ± 0.0164 and variance 29.97 against 30.

The last line shows the memory that is missing: among the sequences in which the first five rolls have failed, you still need on average 5.99 rolls, almost the same 6 as a fresh sequence, and not one fewer as the idea that the 6 is "now due" would suggest. It is the correction to the gambler's reasoning: the die does not remember past rolls. The continuous analogue of this property is the exponential distribution, which you meet in the next article.

When the model does not hold

The formulas in this article require assumptions, and assumptions are the weak point of every application. The binomial presupposes independent trials with the same probability: if the passengers travel in groups, or if the quiz questions are not equally difficult, the real variance differs from npq and the distribution is no longer binomial. The Poisson presupposes rare and independent events: counts with a variance much larger than the mean point to clustering (events in bunches) and call for more flexible families, such as the negative binomial, which I do not cover here.

This article also leaves out the moment generating function, the hypergeometric (draws without replacement) and the estimation of parameters from data: these are topics of inference, where the estimation of p from an observed frequency takes up the calculation of mean and variance done here. For the behaviour of the binomial when n grows with p fixed, the right result is the central limit theorem, in article 05.

Frequently asked questions

What is the difference between binomial and Poisson?

The binomial counts successes in a fixed number n of trials, with an upper limit equal to n. The Poisson has no upper limit and depends on a single parameter, the mean; it approximates the binomial well when n is large and p is small.

Why do mean and variance coincide in a Poisson?

Because the Poisson is the limit of a binomial with p tending to zero: the variance npq tends to np = λ, while the mean stays np = λ. If in real data the variance far exceeds the mean, the Poisson model is not adequate.

Does the geometric start from 0 or from 1?

It depends on the convention. Counting the attempts, success included, it starts from 1 and has mean 1/p; counting only the failures it starts from 0 and has mean (1 − p)/p. NumPy uses the first, like the program above.

What you can do now

You can describe a discrete variable with its PMF and compute its mean and variance from the definition, recognize when a count is binomial and when it is reasonable to approximate it with a Poisson, and quantify the error of that approximation instead of trusting a rule of thumb. You can check a formula in two ways, with exact fractions and with a fixed-seed simulation and declared tolerance. The next article moves everything to the continuous case: density, distribution function, and the uniform, exponential and normal distributions.

The same series also includes Exploratory Data Analysis (EDA) with Python. The full list is on the series page, Probability and Statistics.