Confidence Intervals

A single number from a sample — a mean, a conversion rate, a poll result — is almost certainly not exactly the population value. A confidence interval says how far off it could plausibly be. It's the difference between "support is 42%" and "support is 42% ± 3 points", and it's what makes a result reportable.

Intervals also answer the most common planning question in statistics: how big a sample do I need?

Estimate ± margin of error. Almost every interval has the shape estimate ± critical value × standard error. The standard error comes from Module 3; the critical value comes from a distribution and the confidence level.

What "95% confident" means. The method, used over and over on new samples, produces intervals that contain the true value 95% of the time. It does not mean there's a 95% probability that the true value is in this particular interval — the true value is fixed; it's the interval that varies.

Why the t distribution. If you knew the population σ, the critical value would come from the normal curve (1.96 for 95%). You almost never do, so you plug in the sample s — which is itself uncertain. Student's t distribution accounts for that: it has heavier tails than the normal, so its critical values are bigger, and the extra width shrinks as the sample grows. With n − 1 = 9 degrees of freedom the 95% value is 2.26; with 100 it's 1.98.

Width, confidence, and n. A wider interval is more likely to be right but less useful. You can narrow it by accepting lower confidence or by collecting more data — and because the standard error is σ/√n, halving the width costs four times the data.

Proportions need care near 0 and 1. The textbook (Wald) interval p̂ ± z·√(p̂(1−p̂)/n) misbehaves for small samples or rare events: with 1 success in 20 its lower end is negative. The Wilson interval fixes that and is what statisticians recommend.

• Some functions have no closed-form inverse. The t quantile is found by searching: bracket the answer, then bisect. Bisection halves the bracket every step, so 50 or so steps reach float64 precision — and it can't fail on a monotonic function, unlike Newton's method.
• Grow the bracket instead of guessing it. With one degree of freedom the 0.0005 quantile is −636.6; a fixed search range of ±100 would silently return the wrong answer.
• Return a small struct by value. Interval{Lo, Hi} with Width and Contains methods is clearer than two loose floats and costs nothing to copy.
• r.NormFloat64() draws from the standard normal, so mu + sigma*r.NormFloat64() simulates a measurement — useful for checking an interval's coverage, as below.
• math.Ceil for sample sizes. A sample size of 216.09 means 217 people: round up, then convert with int(...).

package main

import (
	"fmt"
	"math"
	"math/rand/v2"
)

func main() {
	r := rand.New(rand.NewPCG(4, 4))
	const mu, sigma, n, trials = 100.0, 15.0, 25, 1000

	z := math.Sqrt2 * math.Erfinv(0.95) // 1.96: the 0.975 quantile of the standard normal
	half := z * sigma / math.Sqrt(n)    // the margin of error, ±5.88

	hits := 0
	for i := range trials {
		sum := 0.0
		for range n {
			sum += mu + sigma*r.NormFloat64() // one person's score
		}
		xbar := sum / n
		ok := xbar-half <= mu && mu <= xbar+half
		if ok {
			hits++
		}
		if i < 5 {
			fmt.Printf("sample %d: mean %6.2f   interval [%6.2f, %6.2f]   contains 100? %v\n",
				i+1, xbar, xbar-half, xbar+half, ok)
		}
	}
	fmt.Printf("...\n%d of %d intervals contain the true mean (%.1f%%)\n", hits, trials, 100*float64(hits)/trials)
}

It prints:

Each sample gives a different interval, and about 95% of them catch the true mean of 100. Sample 5 only just does. Any one interval either contains 100 or doesn't; the 95% describes the method's long-run hit rate. (962 rather than exactly 950 is ordinary sampling variation; run 100,000 trials and it settles near 95.0%.)

• Saying "there's a 95% chance the mean is in this interval." The true mean isn't random. Say "we're 95% confident" and mean the method's track record.
• Using z* = 1.96 with a small sample and estimated σ. With n = 10, the right multiplier is 2.26; using 1.96 makes the interval about 13% too narrow and its real coverage about 92%.
• Using the Wald interval for rare events. Zero successes in 15 trials gives Wald the interval [0, 0] — "certain it's zero" — which is absurd. Use Wilson.
• Treating non-overlapping intervals as the only sign of a difference. Two 95% intervals can overlap a little while the difference between the means is still clearly non-zero. Test the difference directly (Module 5).

Part A builds Student's t distribution: its CDF from the given incomplete beta function, and its quantile by bisection. Part B builds intervals for a mean and a sample-size planner. Part C compares the Wald and Wilson intervals for a proportion.

go test ./04-confidence-intervals/

Get TCDF right first — TQuantile and TInterval both depend on it. The TQuantile(0.0005, 1) case only passes if your bracket keeps growing.

The files for this module

Everything below is also in the course download, in 04-confidence-intervals/. Work in exercises.go; the tests decide when you’re done.

// Package confint is Module 4: confidence intervals.
//
// An interval estimate says "the population mean is probably between 96.3
// and 104.1" instead of "it's 100.2" — and says how probably. This module
// builds the t distribution, then intervals for means and proportions.
package confint

import (
	"math"

	"github.com/bozcode/go-statistics-and-regression/internal/special"
)

// Interval is a range [Lo, Hi] that, at a stated confidence level, contains
// an unknown population value.
type Interval struct {
	Lo, Hi float64
}

// Width returns Hi − Lo.
func (iv Interval) Width() float64 { return iv.Hi - iv.Lo }

// Contains reports whether x lies inside the interval.
func (iv Interval) Contains(x float64) bool { return iv.Lo <= x && x <= iv.Hi }

// regIncBeta is the regularized incomplete beta function I_x(a, b). It's
// numerical machinery (a continued fraction), so it's given; TCDF is built
// on it.
func regIncBeta(a, b, x float64) float64 { return special.RegIncBeta(a, b, x) }

// The helpers below are given — you wrote them in Modules 1 and 2.

// normalQuantile returns the z with P(Z <= z) = p for a standard normal.
func normalQuantile(p float64) float64 { return special.NormalQuantile(p) }

// mean returns the arithmetic mean of xs.
func mean(xs []float64) float64 {
	sum := 0.0
	for _, x := range xs {
		sum += x
	}
	return sum / float64(len(xs))
}

// stdDev returns the sample standard deviation of xs (divides by n − 1).
func stdDev(xs []float64) float64 {
	m, ss := mean(xs), 0.0
	for _, x := range xs {
		ss += (x - m) * (x - m)
	}
	return math.Sqrt(ss / float64(len(xs)-1))
}

//go:build !solutions

package confint

// ============================================================
// Part A — Student's t distribution
// ============================================================

// TCDF returns P(T <= t) for Student's t distribution with df degrees of
// freedom. The upper tail beyond |t| is ½·I_x(df/2, ½) with
// x = df / (df + t²), using the given regIncBeta; the curve is symmetric, so
// for t >= 0 return 1 − tail, otherwise tail. TCDF(2.0, 10) ≈ 0.9633.
func TCDF(t, df float64) float64 {
	// TODO
	return 0
}

// TQuantile returns the t with TCDF(t, df) = p, for 0 < p < 1. There's no
// closed form, so search: start with [−1, 1], double each end until it
// brackets p, then bisect until the bracket is narrower than 1e-12.
// TQuantile(0.975, 10) ≈ 2.2281.
func TQuantile(p, df float64) float64 {
	// TODO
	return 0
}

// ============================================================
// Part B — intervals for a mean
// ============================================================

// ZInterval returns xbar ± z*·σ/√n, for the rare case where the population
// σ is known. z* = normalQuantile((1 + level)/2): 1.96 for level 0.95.
func ZInterval(xbar, sigma float64, n int, level float64) Interval {
	// TODO
	return Interval{}
}

// TInterval returns x̄ ± t*·s/√n with t* = TQuantile((1 + level)/2, n − 1):
// the everyday interval for a mean when σ is estimated from the sample.
func TInterval(xs []float64, level float64) Interval {
	// TODO
	return Interval{}
}

// SampleSize returns the smallest n whose z-interval has a margin of error
// (half-width) of at most margin: n = ⌈(z*·σ / margin)²⌉.
// SampleSize(15, 2, 0.95) = 217.
func SampleSize(sigma, margin, level float64) int {
	// TODO
	return 0
}

// ============================================================
// Part C — intervals for a proportion
// ============================================================

// WaldInterval returns p̂ ± z*·√(p̂(1−p̂)/n) with p̂ = successes/n — the
// textbook interval. Don't clamp it to [0, 1]: the test checks that it
// misbehaves near 0, which is the point of WilsonInterval.
func WaldInterval(successes, n int, level float64) Interval {
	// TODO
	return Interval{}
}

// WilsonInterval returns the Wilson score interval, which behaves well even
// for small n and proportions near 0 or 1:
//
//	centre = (p̂ + z²/2n) / (1 + z²/n)
//	half   = z/(1 + z²/n) · √(p̂(1−p̂)/n + z²/4n²)
func WilsonInterval(successes, n int, level float64) Interval {
	// TODO
	return Interval{}
}

package confint

import (
	"math"
	"testing"
)

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

func intervalNear(got, want Interval) bool {
	return near(got.Lo, want.Lo) && near(got.Hi, want.Hi)
}

func TestTCDF(t *testing.T) {
	for _, tc := range []struct{ t, df, want float64 }{
		{2.0, 10, 0.9633059826146299},
		{-1.5, 4, 0.104},
		{0.3, 1, 0.5927735790777423},
		{-6.2, 2.5, 0.007067944574802309},
		{3.1, 120, 0.9987941426721137},
		{0, 7, 0.5},
	} {
		if got := TCDF(tc.t, tc.df); !near(got, tc.want) {
			t.Errorf("TCDF(%v, %v) = %v, want %v", tc.t, tc.df, got, tc.want)
		}
	}
}

func TestTQuantile(t *testing.T) {
	for _, tc := range []struct{ p, df, want float64 }{
		{0.975, 10, 2.228138851986274},
		{0.025, 10, -2.2281388519862753},
		{0.995, 3, 5.840909309733355},
		{0.0005, 1, -636.6192487687196}, // one degree of freedom has very heavy tails
		{0.9, 50, 1.2987136941948096},
	} {
		if got := TQuantile(tc.p, tc.df); math.Abs(got-tc.want) > 1e-6*math.Max(1, math.Abs(tc.want)) {
			t.Errorf("TQuantile(%v, %v) = %v, want %v", tc.p, tc.df, got, tc.want)
		}
	}
	if got := TQuantile(0.5, 12); math.Abs(got) > 1e-6 {
		t.Errorf("TQuantile(0.5, 12) = %v, want 0", got)
	}
}

func TestZInterval(t *testing.T) {
	got := ZInterval(103.2, 15, 36, 0.95)
	want := Interval{98.30009003864987, 108.09990996135014}
	if !intervalNear(got, want) {
		t.Errorf("ZInterval(103.2, 15, 36, 0.95) = %+v, want %+v", got, want)
	}
}

var temps = []float64{98.6, 101.2, 97.9, 103.5, 99.1, 100.4, 102.8, 96.7, 101.9, 99.8}

func TestTInterval(t *testing.T) {
	for _, tc := range []struct {
		level float64
		want  Interval
	}{
		{0.95, Interval{98.62547303332073, 101.75452696667926}},
		{0.99, Interval{97.94238633915057, 102.43761366084942}},
	} {
		if got := TInterval(temps, tc.level); math.Abs(got.Lo-tc.want.Lo) > 1e-6 || math.Abs(got.Hi-tc.want.Hi) > 1e-6 {
			t.Errorf("TInterval(temps, %v) = %+v, want %+v", tc.level, got, tc.want)
		}
	}
}

func TestSampleSize(t *testing.T) {
	if got := SampleSize(15, 2, 0.95); got != 217 {
		t.Errorf("SampleSize(15, 2, 0.95) = %d, want 217", got)
	}
	if got := SampleSize(8, 1, 0.99); got != 425 {
		t.Errorf("SampleSize(8, 1, 0.99) = %d, want 425", got)
	}
}

func TestWaldInterval(t *testing.T) {
	if got, want := WaldInterval(42, 120, 0.95), (Interval{0.2646608783681289, 0.43533912163187105}); !intervalNear(got, want) {
		t.Errorf("WaldInterval(42, 120, 0.95) = %+v, want %+v", got, want)
	}
	// 1 success in 20: the Wald interval dips below zero — an impossible proportion
	got := WaldInterval(1, 20, 0.95)
	if !near(got.Hi, 0.14551682940272123) || !near(got.Lo, -0.04551682940272123) {
		t.Errorf("WaldInterval(1, 20, 0.95) = %+v, want {-0.0455 0.1455} (unclamped)", got)
	}
}

func TestWilsonInterval(t *testing.T) {
	for _, tc := range []struct {
		s, n  int
		level float64
		want  Interval
	}{
		{42, 120, 0.95, Interval{0.27051899590285405, 0.4387867540227571}},
		{1, 20, 0.95, Interval{0.008881448800795402, 0.2361311934467421}},
		{0, 15, 0.90, Interval{0, 0.1528077046101403}},
	} {
		if got := WilsonInterval(tc.s, tc.n, tc.level); !intervalNear(got, tc.want) {
			t.Errorf("WilsonInterval(%d, %d, %v) = %+v, want %+v", tc.s, tc.n, tc.level, got, tc.want)
		}
	}
}

//go:build solutions

package confint

import "math"

// ============================================================
// Part A — Student's t distribution
// ============================================================

func TCDF(t, df float64) float64 {
	// Near t = 0, df/(df+t²) rounds to exactly 1 and the tail loses its
	// digits, so use the mirror identity I_x(a,b) = 1 − I_(1−x)(b,a) there.
	var tail float64
	if t*t < df {
		tail = 0.5 * (1 - regIncBeta(0.5, df/2, t*t/(df+t*t)))
	} else {
		tail = 0.5 * regIncBeta(df/2, 0.5, df/(df+t*t))
	}
	if t >= 0 {
		return 1 - tail
	}
	return tail
}

func TQuantile(p, df float64) float64 {
	lo, hi := -1.0, 1.0
	for TCDF(lo, df) > p {
		lo *= 2
	}
	for TCDF(hi, df) < p {
		hi *= 2
	}
	for hi-lo > 1e-12*math.Max(1, math.Abs(lo)) {
		mid := (lo + hi) / 2
		if TCDF(mid, df) < p {
			lo = mid
		} else {
			hi = mid
		}
	}
	return (lo + hi) / 2
}

// ============================================================
// Part B — intervals for a mean
// ============================================================

func ZInterval(xbar, sigma float64, n int, level float64) Interval {
	z := normalQuantile((1 + level) / 2)
	half := z * sigma / math.Sqrt(float64(n))
	return Interval{xbar - half, xbar + half}
}

func TInterval(xs []float64, level float64) Interval {
	n := len(xs)
	t := TQuantile((1+level)/2, float64(n-1))
	half := t * stdDev(xs) / math.Sqrt(float64(n))
	m := mean(xs)
	return Interval{m - half, m + half}
}

func SampleSize(sigma, margin, level float64) int {
	z := normalQuantile((1 + level) / 2)
	return int(math.Ceil(math.Pow(z*sigma/margin, 2)))
}

// ============================================================
// Part C — intervals for a proportion
// ============================================================

func WaldInterval(successes, n int, level float64) Interval {
	p := float64(successes) / float64(n)
	z := normalQuantile((1 + level) / 2)
	half := z * math.Sqrt(p*(1-p)/float64(n))
	return Interval{p - half, p + half}
}

func WilsonInterval(successes, n int, level float64) Interval {
	p := float64(successes) / float64(n)
	nf := float64(n)
	z := normalQuantile((1 + level) / 2)
	denom := 1 + z*z/nf
	centre := (p + z*z/(2*nf)) / denom
	half := z / denom * math.Sqrt(p*(1-p)/nf+z*z/(4*nf*nf))
	return Interval{centre - half, centre + half}
}

• Write TwoMeanInterval(a, b, level) for the difference of two means, using the Welch standard error √(s₁²/n₁ + s₂²/n₂) (degrees of freedom are in Module 5).
• Simulate the coverage of the Wald interval for n = 20 and p = 0.05. It's far below the advertised 95%; Wilson's is close.
• Practice: a poll finds 540 of 1,200 voters in favour. Compute both intervals and the sample size needed for a ±2-point margin.

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 5, so git log shows your own progress through the course.