Probability Distributions

Statistics is mostly the art of asking "how surprising is this?" — and you can't answer that without a model of what unsurprising looks like. A probability distribution is that model. The binomial describes counts of successes, the Poisson describes counts of rare events, and the normal describes measurements that are the sum of many small influences: heights, test scores, measurement error, and — as Module 3 shows — averages of almost anything.

Every p-value, confidence interval, and regression standard error in the rest of this course comes from one of these curves or a close relative.

Discrete versus continuous. A discrete distribution gives each whole number a probability: P(X = 3) for "three heads". A continuous one gives a density — a height on a curve — and probabilities are areas: P(1.8 m < height < 1.9 m). The chance of any single exact value is zero.

PMF, PDF, and CDF. The probability mass function (PMF) gives P(X = k) for a discrete variable; the probability density function (PDF) gives the curve's height for a continuous one. The cumulative distribution function (CDF) gives P(X ≤ x) for both, and it's the one you'll use most: every "what fraction falls below…" question is a CDF.

The quantile function runs the CDF backwards: given a probability p, which value has fraction p below it? The 0.975 quantile of the standard normal is 1.96, which is where the famous "±1.96" in 95% confidence intervals comes from.

Mean and spread of a distribution. A binomial with N trials and success probability P has mean N·P and variance N·P·(1 − P). A Poisson's mean and variance are both λ. A normal is completely described by its mean μ and standard deviation σ.

The normal approximation. For large N, a binomial looks like a normal curve with the same mean and standard deviation. Because the binomial's bar for k spans k − ½ to k + ½, adding that half (the continuity correction) makes the approximation much closer.

• The math package has what you need. math.Erfc and math.Erfinv give the normal CDF and quantile exactly; math.Lgamma gives log(n!) as Lgamma(n+1) without ever computing n!.
• Work in log space for products of many factors. 171! already overflows a float64 to +Inf, and 0.5²⁰⁰⁰ underflows to 0. Add logarithms, then math.Exp once at the end.
• math.Lgamma returns two values — the log of |Γ(x)| and its sign — so write lg, _ := math.Lgamma(x).
• Use math.Log1p(-p) for log(1 − p) when p is tiny; plain math.Log(1-p) loses digits because 1 − p rounds.
• Use Erfc, not 1 − something, in the tails. Φ(−9) is about 1.1×10⁻¹⁹; 1 − Φ(9) computed the naive way is exactly 0.
• Methods on value types suit small structs like Normal{Mu, Sigma}: StandardNormal.CDF(1.96) reads like the maths.

package main

import (
	"fmt"
	"math"
)

// binomPMF is C(n,k)·pᵏ·(1−p)ⁿ⁻ᵏ, computed in log space.
func binomPMF(n, k int, p float64) float64 {
	a, _ := math.Lgamma(float64(n + 1))
	b, _ := math.Lgamma(float64(k + 1))
	c, _ := math.Lgamma(float64(n - k + 1))
	return math.Exp(a - b - c + float64(k)*math.Log(p) + float64(n-k)*math.Log1p(-p))
}

// normalCDF is P(X <= x) for a normal with mean mu and sd sigma.
func normalCDF(x, mu, sigma float64) float64 {
	return 0.5 * math.Erfc(-(x-mu)/(sigma*math.Sqrt2))
}

func main() {
	n, p := 100, 0.5 // a fair coin, flipped 100 times

	exact := 0.0 // P(60 or more heads), adding up the bars
	for k := 60; k <= n; k++ {
		exact += binomPMF(n, k, p)
	}

	mu, sigma := float64(n)*p, math.Sqrt(float64(n)*p*(1-p)) // 50 and 5
	plain := 1 - normalCDF(60, mu, sigma)                    // the curve from 60 up
	corrected := 1 - normalCDF(59.5, mu, sigma)              // the bar for 60 starts at 59.5

	fmt.Printf("exact binomial          P(X >= 60) = %.5f\n", exact)
	fmt.Printf("normal, no correction   P(X >= 60) = %.5f\n", plain)
	fmt.Printf("normal, continuity fix  P(X >= 60) = %.5f\n", corrected)
}

It prints:

A fair coin lands 60 or more heads about 2.8% of the time — unusual, but not wildly so. The plain normal approximation underestimates that by a fifth; with the continuity correction it's within 0.0003 of the exact answer. That 2.8% is exactly the kind of number Module 5 will call a p-value.

• Reading a density as a probability. A normal with σ = 0.2 has a peak of almost 2. Only areas are probabilities.
• Using the normal approximation for small samples or extreme P. The rule of thumb is N·P ≥ 10 and N·(1 − P) ≥ 10. Below that, compute the binomial exactly — it's cheap.
• Forgetting independence. The binomial assumes each trial is independent with the same P. Survey answers from the same household, or defects from the same faulty machine, aren't.
• Using the Poisson when events cluster. If the variance of your counts is much larger than their mean, events aren't independent, and the Poisson will understate how variable they are.

Part A is the discrete distributions — binomial and Poisson PMFs and CDFs, in log space. Part B is the normal: density, CDF, quantile, and the area between two points. Part C compares the normal approximation with the exact binomial.

go test ./02-probability-distributions/

Two tests use numbers too big for direct computation — Binomial{2000, 0.5}.PMF(1000) and Poisson{230}.PMF(250) — so a version that multiplies factorials will fail them.

The files for this module

Everything below is also in the course download, in 02-probability-distributions/. Work in exercises.go; the tests decide when you’re done.

// Package dist is Module 2: probability distributions.
//
// A distribution is a rule that says how likely each outcome is. Discrete
// ones (Binomial, Poisson) give a probability to each whole number; the
// Normal is continuous, so it gives a density, and probabilities come from
// areas under the curve.
package dist

// Binomial counts successes in N independent trials, each a success with
// probability P: heads in N coin flips, defective parts in a batch,
// sign-ups among N visitors.
type Binomial struct {
	N int
	P float64
}

// Poisson counts events in a fixed stretch of time or space when they occur
// independently at an average rate Lambda: calls per hour, typos per page.
type Poisson struct {
	Lambda float64
}

// Normal is the bell curve with mean Mu and standard deviation Sigma.
type Normal struct {
	Mu, Sigma float64
}

// StandardNormal has mean 0 and standard deviation 1. Its values are
// z-scores.
var StandardNormal = Normal{Mu: 0, Sigma: 1}

//go:build !solutions

package dist

// ============================================================
// Part A — discrete distributions
// ============================================================

// PMF returns P(X = k) = C(N, k)·Pᵏ·(1−P)^(N−k), and 0 when k < 0 or k > N.
// Compute it in log space with math.Lgamma: C(2000, 1000) ≈ 10⁶⁰⁰ overflows
// a float64 and 0.5²⁰⁰⁰ underflows to 0, but their logarithms are ordinary
// numbers. Handle P = 0 and P = 1 separately
// (log 0 is −Inf, and 0·−Inf is NaN). Binomial{10, 0.5}.PMF(5) ≈ 0.2461.
func (b Binomial) PMF(k int) float64 {
	// TODO
	return 0
}

// CDF returns P(X <= k): the sum of PMF(0) … PMF(k). 0 for k < 0, 1 for
// k >= N. Target: O(k).
func (b Binomial) CDF(k int) float64 {
	// TODO
	return 0
}

// Mean returns N·P, the expected number of successes.
func (b Binomial) Mean() float64 {
	// TODO
	return 0
}

// Variance returns N·P·(1−P).
func (b Binomial) Variance() float64 {
	// TODO
	return 0
}

// PMF returns P(X = k) = λᵏ·e^(−λ) / k!, and 0 for k < 0. Log space again:
// exp(k·ln λ − λ − lgamma(k+1)). Poisson{3}.PMF(2) ≈ 0.2240.
func (p Poisson) PMF(k int) float64 {
	// TODO
	return 0
}

// CDF returns P(X <= k), or 0 for k < 0.
func (p Poisson) CDF(k int) float64 {
	// TODO
	return 0
}

// ============================================================
// Part B — the normal distribution
// ============================================================

// PDF returns the density at x: exp(−z²/2) / (σ·√(2π)), z = (x−μ)/σ.
// A density is not a probability — it can exceed 1 when σ is small.
func (d Normal) PDF(x float64) float64 {
	// TODO
	return 0
}

// CDF returns P(X <= x) = ½·erfc(−z/√2), z = (x−μ)/σ. Use math.Erfc, not
// 1 − ½·erfc: it stays accurate far out in the tails.
// StandardNormal.CDF(1.96) ≈ 0.9750.
func (d Normal) CDF(x float64) float64 {
	// TODO
	return 0
}

// Quantile returns the x with CDF(x) = p, the inverse of CDF:
// μ + σ·√2·erfinv(2p − 1). Quantile(0.975) of the standard normal ≈ 1.96.
// p = 0 and p = 1 give −Inf and +Inf; p outside [0, 1] gives NaN.
func (d Normal) Quantile(p float64) float64 {
	// TODO
	return 0
}

// Between returns P(a <= X <= b) = CDF(b) − CDF(a).
func (d Normal) Between(a, b float64) float64 {
	// TODO
	return 0
}

// ============================================================
// Part C — the normal approximation
// ============================================================

// NormalApprox returns the normal curve with the binomial's mean and
// standard deviation: μ = N·P, σ = √(N·P·(1−P)). Good when N·P and
// N·(1−P) are both at least about 10.
func NormalApprox(b Binomial) Normal {
	// TODO
	return Normal{}
}

// ApproxCDF estimates the binomial P(X <= k) with the normal approximation
// and a continuity correction: the bar for k covers k−½ to k+½, so use
// NormalApprox(b).CDF(k + 0.5).
func ApproxCDF(b Binomial, k int) float64 {
	// TODO
	return 0
}

package dist

import (
	"math"
	"testing"
)

// near reports whether a and b agree to about 9 significant digits (or 1e-9
// absolutely, for values near zero). Expected values come from SciPy.
func near(a, b float64) bool {
	return math.Abs(a-b) <= 1e-9*math.Max(1, math.Abs(b))
}

// relNear is near for very small numbers: 9 significant digits, no floor.
func relNear(a, b float64) bool {
	return math.Abs(a-b) <= 1e-9*math.Abs(b)
}

func TestBinomialPMF(t *testing.T) {
	for _, tc := range []struct {
		b    Binomial
		k    int
		want float64
	}{
		{Binomial{10, 0.5}, 0, 0.0009765625},
		{Binomial{10, 0.5}, 5, 0.24609375},
		{Binomial{20, 0.15}, 3, 0.24282889614926756},
		{Binomial{2000, 0.5}, 1000, 0.017839011145854313}, // C(2000,1000) overflows, 0.5^2000 underflows
		{Binomial{5, 0}, 0, 1},
		{Binomial{5, 1}, 5, 1},
		{Binomial{5, 0.3}, 6, 0},
		{Binomial{5, 0.3}, -1, 0},
	} {
		if got := tc.b.PMF(tc.k); !relNear(got, tc.want) && !(tc.want == 0 && got == 0) {
			t.Errorf("%+v.PMF(%d) = %v, want %v", tc.b, tc.k, got, tc.want)
		}
	}
}

func TestBinomialCDF(t *testing.T) {
	for _, tc := range []struct {
		b    Binomial
		k    int
		want float64
	}{
		{Binomial{20, 0.15}, 3, 0.6477251741567033},
		{Binomial{100, 0.5}, 59, 0.9715560331795096},
		{Binomial{12, 0.3}, 7, 0.99051062887},
		{Binomial{12, 0.3}, -1, 0},
		{Binomial{12, 0.3}, 12, 1},
	} {
		if got := tc.b.CDF(tc.k); !near(got, tc.want) {
			t.Errorf("%+v.CDF(%d) = %v, want %v", tc.b, tc.k, got, tc.want)
		}
	}
}

func TestBinomialMoments(t *testing.T) {
	b := Binomial{N: 40, P: 0.25}
	if got := b.Mean(); !near(got, 10) {
		t.Errorf("Mean = %v, want 10", got)
	}
	if got := b.Variance(); !near(got, 7.5) {
		t.Errorf("Variance = %v, want 7.5", got)
	}
}

func TestPoisson(t *testing.T) {
	for _, tc := range []struct {
		p    Poisson
		k    int
		want float64
	}{
		{Poisson{3}, 0, 0.049787068367863944},
		{Poisson{3}, 2, 0.22404180765538775},
		{Poisson{230}, 250, 0.010830332030310759}, // 230^250 overflows; its log doesn't
		{Poisson{3}, -2, 0},
	} {
		if got := tc.p.PMF(tc.k); !relNear(got, tc.want) && !(tc.want == 0 && got == 0) {
			t.Errorf("%+v.PMF(%d) = %v, want %v", tc.p, tc.k, got, tc.want)
		}
	}
	if got := (Poisson{3}).CDF(4); !near(got, 0.8152632445237722) {
		t.Errorf("Poisson{3}.CDF(4) = %v, want 0.8152632445237722", got)
	}
	if got := (Poisson{7.5}).CDF(10); !near(got, 0.8622379834283881) {
		t.Errorf("Poisson{7.5}.CDF(10) = %v, want 0.8622379834283881", got)
	}
}

func TestNormalPDF(t *testing.T) {
	if got := StandardNormal.PDF(0); !near(got, 0.3989422804014327) {
		t.Errorf("StandardNormal.PDF(0) = %v, want 0.3989422804014327", got)
	}
	if got := (Normal{100, 15}).PDF(108); !near(got, 0.023070259545128195) {
		t.Errorf("Normal{100,15}.PDF(108) = %v, want 0.023070259545128195", got)
	}
	// a density, not a probability: it can be bigger than 1
	if got := (Normal{10, 0.2}).PDF(10); !near(got, 1.994711402007164) {
		t.Errorf("Normal{10,0.2}.PDF(10) = %v, want 1.994711402007164", got)
	}
}

func TestNormalCDF(t *testing.T) {
	for _, tc := range []struct {
		d    Normal
		x    float64
		want float64
	}{
		{StandardNormal, 1.96, 0.9750021048517795},
		{StandardNormal, -1, 0.15865525393145707},
		{Normal{100, 15}, 130, 0.9772498680518208},
		{StandardNormal, 0, 0.5},
	} {
		if got := tc.d.CDF(tc.x); !near(got, tc.want) {
			t.Errorf("%+v.CDF(%v) = %v, want %v", tc.d, tc.x, got, tc.want)
		}
	}
	// far tail: 1 − Φ(9) would round to 0; erfc keeps the digits
	if got := StandardNormal.CDF(-9); !relNear(got, 1.1285884059538324e-19) {
		t.Errorf("StandardNormal.CDF(-9) = %v, want 1.1285884059538324e-19", got)
	}
}

func TestNormalQuantile(t *testing.T) {
	for _, tc := range []struct {
		d       Normal
		p, want float64
	}{
		{StandardNormal, 0.975, 1.959963984540054},
		{StandardNormal, 0.05, -1.6448536269514729},
		{Normal{100, 15}, 0.9, 119.22327348316901},
		{StandardNormal, 0.5, 0},
	} {
		if got := tc.d.Quantile(tc.p); !near(got, tc.want) {
			t.Errorf("%+v.Quantile(%v) = %v, want %v", tc.d, tc.p, got, tc.want)
		}
	}
	// round trip
	d := Normal{-3, 2.5}
	for _, x := range []float64{-9, -3, 0.4, 5} {
		if got := d.Quantile(d.CDF(x)); !near(got, x) {
			t.Errorf("Quantile(CDF(%v)) = %v", x, got)
		}
	}
	if got := StandardNormal.Quantile(1.5); !math.IsNaN(got) {
		t.Errorf("Quantile(1.5) = %v, want NaN", got)
	}
}

func TestBetween(t *testing.T) {
	// the 68–95–99.7 rule
	for _, tc := range []struct{ k, want float64 }{
		{1, 0.6826894921370859}, {2, 0.9544997361036416}, {3, 0.9973002039367398},
	} {
		if got := StandardNormal.Between(-tc.k, tc.k); !near(got, tc.want) {
			t.Errorf("Between(-%v, %v) = %v, want %v", tc.k, tc.k, got, tc.want)
		}
	}
	if got := (Normal{100, 15}).Between(85, 115); !near(got, 0.6826894921370859) {
		t.Errorf("Normal{100,15}.Between(85, 115) = %v, want 0.6827", got)
	}
}

func TestNormalApprox(t *testing.T) {
	n := NormalApprox(Binomial{100, 0.5})
	if !near(n.Mu, 50) || !near(n.Sigma, 5) {
		t.Fatalf("NormalApprox(Binomial{100, 0.5}) = %+v, want {50 5}", n)
	}
	// exact 0.97156; with the continuity correction 0.97128; without it 0.96407
	if got := ApproxCDF(Binomial{100, 0.5}, 59); !near(got, 0.9712834401839981) {
		t.Errorf("ApproxCDF(Binomial{100, 0.5}, 59) = %v, want 0.9712834401839981", got)
	}
	if got := ApproxCDF(Binomial{12, 0.3}, 7); !near(got, 0.99299036144302) {
		t.Errorf("ApproxCDF(Binomial{12, 0.3}, 7) = %v, want 0.99299036144302", got)
	}
}

//go:build solutions

package dist

import "math"

// ============================================================
// Part A — discrete distributions
// ============================================================

func (b Binomial) PMF(k int) float64 {
	if k < 0 || k > b.N {
		return 0
	}
	switch b.P {
	case 0:
		return boolProb(k == 0)
	case 1:
		return boolProb(k == b.N)
	}
	return math.Exp(logChoose(b.N, k) + float64(k)*math.Log(b.P) + float64(b.N-k)*math.Log1p(-b.P))
}

func boolProb(ok bool) float64 {
	if ok {
		return 1
	}
	return 0
}

// logChoose is log C(n, k).
func logChoose(n, k int) float64 {
	a, _ := math.Lgamma(float64(n + 1))
	b, _ := math.Lgamma(float64(k + 1))
	c, _ := math.Lgamma(float64(n - k + 1))
	return a - b - c
}

func (b Binomial) CDF(k int) float64 {
	if k < 0 {
		return 0
	}
	if k >= b.N {
		return 1
	}
	sum := 0.0
	for i := 0; i <= k; i++ {
		sum += b.PMF(i)
	}
	return sum
}

func (b Binomial) Mean() float64 { return float64(b.N) * b.P }

func (b Binomial) Variance() float64 { return float64(b.N) * b.P * (1 - b.P) }

func (p Poisson) PMF(k int) float64 {
	if k < 0 {
		return 0
	}
	lf, _ := math.Lgamma(float64(k + 1))
	return math.Exp(float64(k)*math.Log(p.Lambda) - p.Lambda - lf)
}

func (p Poisson) CDF(k int) float64 {
	sum := 0.0
	for i := 0; i <= k; i++ {
		sum += p.PMF(i)
	}
	return sum
}

// ============================================================
// Part B — the normal distribution
// ============================================================

func (d Normal) PDF(x float64) float64 {
	z := (x - d.Mu) / d.Sigma
	return math.Exp(-z*z/2) / (d.Sigma * math.Sqrt(2*math.Pi))
}

func (d Normal) CDF(x float64) float64 {
	return 0.5 * math.Erfc(-(x-d.Mu)/(d.Sigma*math.Sqrt2))
}

func (d Normal) Quantile(p float64) float64 {
	return d.Mu + d.Sigma*math.Sqrt2*math.Erfinv(2*p-1)
}

func (d Normal) Between(a, b float64) float64 {
	return d.CDF(b) - d.CDF(a)
}

// ============================================================
// Part C — the normal approximation
// ============================================================

func NormalApprox(b Binomial) Normal {
	return Normal{Mu: b.Mean(), Sigma: math.Sqrt(b.Variance())}
}

func ApproxCDF(b Binomial, k int) float64 {
	return NormalApprox(b).CDF(float64(k) + 0.5)
}

• Add an Exponential{Rate} distribution — the waiting time between Poisson events — with CDF(x) = 1 − e^(−rate·x).
• Plot the binomial bars and the approximating normal curve for N = 10, 30 and 100 (print a text histogram), and watch them converge.
• Practice: if a website converts 3% of visitors, how likely is a day with no sign-ups out of 200 visitors? (Binomial, and Poisson with λ = 6 — compare.)

Four questions. Answer each one in your head — or out loud — before you open the answer.

Commit your work before you move on to module 3, so git log shows your own progress through the course.