Hypothesis Testing

Version B of a page converted better. The new drug group recovered faster. The die came up six a lot. Each time the question is the same: is this a real effect, or the noise you'd expect from a sample this size? Hypothesis tests give a disciplined answer, and they're the backbone of A/B testing, clinical trials, quality control, and the coefficient tables every regression prints.

They're also the most misused tool in statistics, so this module spends as much time on how not to fool yourself as on the formulas.

The null hypothesis. Start by assuming nothing is going on: the means are equal, the die is fair, the variables are independent. That's H₀. The alternative, H₁, is what you suspect.

The test statistic measures how far the data is from what H₀ predicts, in standardized units. A t statistic is "difference ÷ its standard error"; a chi-square statistic adds up squared gaps between observed and expected counts.

The p-value is the probability of a statistic at least as extreme as yours if H₀ were true. A small p-value means the data would be surprising under H₀. It is not the probability that H₀ is true, and 1 − p is not the probability that your effect is real.

α and the two kinds of error. Decide a threshold α (often 0.05) before looking. Rejecting a true H₀ is a Type I error (a false alarm, rate α); failing to reject a false H₀ is a Type II error (a miss). A test's power is its chance of catching a real effect, and it grows with sample size and effect size.

Choosing a t-test. One sample against a fixed value: one-sample t. Two independent groups: Welch's t, which doesn't assume equal variances and is the safe default. The same subjects measured twice: the paired t-test, which works on the differences and is far more sensitive because it removes person-to-person variation.

Counts get chi-square tests. Goodness of fit compares observed counts with expected ones (is the die fair?). The test of independence asks whether two categorical variables are related (does the choice of plan depend on the region?).

Many tests, many false alarms. At α = 0.05, twenty tests of pure noise produce about one "discovery" on average. Bonferroni divides α by the number of tests (strict); Benjamini–Hochberg controls the false discovery rate — the expected share of false positives among your discoveries — and keeps more power.

Effect size. A p-value mixes the size of an effect with the size of the sample. Cohen's d — the difference in means in standard-deviation units — reports the size alone.

• Compute small tail probabilities directly. 2*tCDF(-math.Abs(t), df) keeps a p-value of 1e-12 accurate; 2*(1-tCDF(math.Abs(t), df)) subtracts two numbers near 1 and can round it to 0.
• Return a struct, not three floats. Result{Statistic, DF, P} keeps the numbers together and lets you add methods like Significant(alpha).
• Sort indices, not values, when results must come back in input order: build idx := []int{0, 1, …} and sort it with slices.SortFunc comparing ps[a] and ps[b].
• Don't modify the caller's slice. Adjusting p-values in place would silently change data the caller still needs; allocate a result slice.
• Integer counts, float expectations. Observed counts are []int; expected counts usually aren't whole numbers. Convert with float64(o) before subtracting.

package main

import (
	"fmt"
	"math"
)

// zTest compares two group means when the population sd is known (sigma),
// returning the z statistic, the two-sided p-value, and Cohen's d.
func zTest(diff, sigma float64, nPerGroup int) (z, p, d float64) {
	se := sigma * math.Sqrt(2/float64(nPerGroup)) // sd of the difference in means
	z = diff / se
	p = math.Erfc(math.Abs(z) / math.Sqrt2) // = 2·P(Z > |z|)
	d = diff / sigma
	return z, p, d
}

func main() {
	// Two made-up studies of a test whose scores have sd 15.
	studies := []struct {
		name string
		diff float64 // difference between the group means
		n    int     // people per group
	}{
		{"huge study, tiny effect ", 0.5, 100000},
		{"small study, large effect", 8, 20},
	}
	for _, s := range studies {
		z, p, d := zTest(s.diff, 15, s.n)
		fmt.Printf("%s  diff %4.1f  n %6d  z %5.2f  p %.2g  d %.2f  significant at 0.05? %v\n",
			s.name, s.diff, s.n, z, p, d, p < 0.05)
	}
}

It prints:

The first study's difference — half a point on a test with a standard deviation of 15 — is wildly "significant" because 100,000 people per group can detect almost anything. It is also almost certainly meaningless (d = 0.03). The second shows a sizeable effect (d = 0.53) that misses 0.05 because 20 people per group is too few to be sure. Report the effect size and an interval, not just the p-value.

• Reading p as the probability the null is true. p = 0.03 means "data this extreme would happen 3% of the time if H₀ were true", nothing more.
• Treating p = 0.051 and p = 0.049 as different worlds. They're the same evidence. α is a decision rule, not a law of nature.
• Testing until something is significant. Peeking at an A/B test daily and stopping when p < 0.05, or trying twenty outcome measures, guarantees false discoveries. Fix the plan in advance or correct for multiple tests.
• Using an independent-samples test on paired data. Before/after measurements on the same people need the paired test; the unpaired test throws away the pairing and usually misses the effect.
• Concluding "no effect" from a non-significant result. Absence of evidence isn't evidence of absence — especially with small samples. Look at the confidence interval: if it's wide, you simply don't know.

Part A is the three t-tests, all built on the given tCDF. Part B is the two chi-square tests for counts. Part C adjusts p-values for multiple testing and measures effect size.

go test ./05-hypothesis-testing/

Every expected value comes from SciPy (ttest_1samp, ttest_ind(equal_var=False), ttest_rel, chisquare, chi2_contingency(correction=False)) or statsmodels (multipletests), so passing means your numbers match the standard tools to nine digits.

The files for this module

Everything below is also in the course download, in 05-hypothesis-testing/. Work in exercises.go; the tests decide when you’re done.

// Package hypothesis is Module 5: hypothesis tests.
//
// Every test here follows one recipe: assume nothing is going on (the null
// hypothesis), compute a statistic measuring how far the data is from that,
// and turn it into a p-value — the probability of data at least that
// extreme if the null were true.
package hypothesis

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

// Result is the outcome of a test.
type Result struct {
	Statistic float64 // t or χ²
	DF        float64 // degrees of freedom
	P         float64 // p-value (two-sided for t tests)
}

// Significant reports whether the p-value is below alpha.
func (r Result) Significant(alpha float64) bool { return r.P < alpha }

// The helpers below are given. tCDF is what you built in Module 4;
// chiSquareCDF is P(X <= x) for a chi-square with k degrees of freedom,
// built on the incomplete gamma function in internal/special.

func tCDF(t, df float64) float64        { return special.TCDF(t, df) }
func chiSquareCDF(x, k float64) float64 { return special.ChiSquareCDF(x, k) }

// 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))
}

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

//go:build !solutions

package hypothesis

// ============================================================
// Part A — t tests
// ============================================================

// OneSampleT tests whether the population mean equals mu0:
// t = (x̄ − mu0) / (s/√n), df = n − 1, and a two-sided p-value
// 2·P(T <= −|t|). Matches SciPy's ttest_1samp.
func OneSampleT(xs []float64, mu0 float64) Result {
	// TODO
	return Result{}
}

// WelchT tests whether two independent groups have the same mean, without
// assuming equal variances (Welch's t-test — the safe default):
//
//	va, vb = s²/n for each group
//	t  = (x̄a − x̄b) / √(va + vb)
//	df = (va + vb)² / (va²/(na−1) + vb²/(nb−1))
//
// Matches SciPy's ttest_ind(a, b, equal_var=False).
func WelchT(a, b []float64) Result {
	// TODO
	return Result{}
}

// PairedT tests paired measurements (the same subjects, before and after):
// it's OneSampleT on the differences a[i] − b[i] with mu0 = 0. Matches
// SciPy's ttest_rel(a, b).
func PairedT(a, b []float64) Result {
	// TODO
	return Result{}
}

// ============================================================
// Part B — chi-square tests for counts
// ============================================================

// ChiSquareGOF is the goodness-of-fit test: do the observed counts match
// the expected counts? χ² = Σ (O − E)² / E, df = k − 1, and
// p = 1 − chiSquareCDF(χ², df). Matches SciPy's chisquare.
func ChiSquareGOF(observed []int, expected []float64) Result {
	// TODO
	return Result{}
}

// ChiSquareIndependence tests whether the rows and columns of a contingency
// table are independent. Each cell's expected count is
// rowTotal·colTotal / grandTotal; χ² = Σ (O − E)² / E over every cell,
// df = (rows − 1)(cols − 1). Matches SciPy's
// chi2_contingency(table, correction=False).
func ChiSquareIndependence(table [][]int) Result {
	// TODO
	return Result{}
}

// ============================================================
// Part C — many tests, and how big the effect is
// ============================================================

// Bonferroni adjusts p-values for running m tests: each becomes
// min(1, p·m). Compare the adjusted values with alpha as usual. Returns a
// new slice in the same order.
func Bonferroni(ps []float64) []float64 {
	// TODO
	return nil
}

// BenjaminiHochberg adjusts p-values to control the false discovery rate.
// Sort the p-values ascending; the k-th smallest (k = 1…m) becomes
// p·m/k; then, from the largest down, replace each with the minimum of
// itself and the one after it, and cap at 1. Return them in the ORIGINAL
// order. Matches statsmodels' multipletests(method="fdr_bh").
func BenjaminiHochberg(ps []float64) []float64 {
	// TODO
	return nil
}

// CohenD returns the standardized difference between two group means,
// (x̄a − x̄b) / s_pooled, with
// s_pooled = √(((na−1)·s²a + (nb−1)·s²b) / (na + nb − 2)).
// Roughly: 0.2 is small, 0.5 medium, 0.8 large.
func CohenD(a, b []float64) float64 {
	// TODO
	return 0
}

package hypothesis

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 check(t *testing.T, name string, got, want Result) {
	t.Helper()
	if !near(got.Statistic, want.Statistic) || !near(got.DF, want.DF) || !near(got.P, want.P) {
		t.Errorf("%s = %+v, want %+v", name, got, want)
	}
}

// bottle fills in ml (made-up data); the label says 500
var fills = []float64{498.2, 501.1, 497.6, 499.3, 498.8, 500.2, 497.9, 499.5, 498.4, 499.0, 497.2, 499.8}

// seconds to finish a task on two page designs (made-up data)
var designA = []float64{34.1, 29.8, 41.2, 37.5, 30.9, 36.4, 33.0, 39.8}
var designB = []float64{28.4, 31.2, 26.9, 30.1, 27.5, 33.8, 29.0, 25.6, 31.7, 28.8}

func TestOneSampleT(t *testing.T) {
	check(t, "OneSampleT(fills, 500)", OneSampleT(fills, 500),
		Result{Statistic: -3.3033597986209187, DF: 11, P: 0.007035921490324698})
	r := OneSampleT([]float64{4.9, 5.1, 5.0, 4.8, 5.2}, 5)
	if math.Abs(r.Statistic) > 1e-9 || !near(r.P, 1) {
		t.Errorf("a sample centred exactly on mu0: got %+v, want t = 0, p = 1", r)
	}
}

func TestWelchT(t *testing.T) {
	check(t, "WelchT(designA, designB)", WelchT(designA, designB),
		Result{Statistic: 3.6755093164374473, DF: 10.89200578962794, P: 0.00371382564279196})
	// swapping the groups flips the sign and nothing else
	r := WelchT(designB, designA)
	if !near(r.Statistic, -3.6755093164374473) || !near(r.P, 0.00371382564279196) {
		t.Errorf("WelchT(designB, designA) = %+v", r)
	}
}

func TestPairedT(t *testing.T) {
	before := []float64{142, 138, 150, 145, 160, 155, 139, 148}
	after := []float64{136, 135, 144, 146, 151, 149, 137, 141}
	check(t, "PairedT(before, after)", PairedT(before, after),
		Result{Statistic: 4.203727949252743, DF: 7, P: 0.0040169099978379225})
}

func TestChiSquareGOF(t *testing.T) {
	// is this die fair? 100 rolls
	fair := []float64{100.0 / 6, 100.0 / 6, 100.0 / 6, 100.0 / 6, 100.0 / 6, 100.0 / 6}
	check(t, "ChiSquareGOF(die)", ChiSquareGOF([]int{16, 18, 16, 14, 12, 24}, fair),
		Result{Statistic: 5.12, DF: 5, P: 0.4014115932460956})
	expected := []float64{40, 20, 30, 40, 50, 20}
	check(t, "ChiSquareGOF(uneven)", ChiSquareGOF([]int{30, 14, 34, 45, 57, 20}, expected),
		Result{Statistic: 6.438333333333333, DF: 5, P: 0.26587161209871035})
}

func TestChiSquareIndependence(t *testing.T) {
	check(t, "2×3 table", ChiSquareIndependence([][]int{{42, 33, 25}, {18, 37, 45}}),
		Result{Statistic: 15.542857142857141, DF: 2, P: 0.00042161053562412366})
	check(t, "2×2 table", ChiSquareIndependence([][]int{{30, 10}, {20, 25}}),
		Result{Statistic: 8.16269841269841, DF: 1, P: 0.004276070303601571})
}

var ps = []float64{0.001, 0.008, 0.039, 0.041, 0.042, 0.06, 0.074, 0.205, 0.212, 0.216}

func TestBonferroni(t *testing.T) {
	want := []float64{0.01, 0.08, 0.39, 0.41, 0.42, 0.6, 0.74, 1, 1, 1}
	got := Bonferroni(ps)
	if len(got) != len(want) {
		t.Fatalf("Bonferroni returned %d values, want %d", len(got), len(want))
	}
	for i := range want {
		if !near(got[i], want[i]) {
			t.Fatalf("Bonferroni = %v, want %v", got, want)
		}
	}
	if ps[0] != 0.001 {
		t.Error("Bonferroni modified its input")
	}
}

func TestBenjaminiHochberg(t *testing.T) {
	// deliberately unsorted: results must come back in input order
	in := []float64{0.042, 0.205, 0.001, 0.074, 0.039, 0.216, 0.008, 0.06, 0.212, 0.041}
	want := []float64{0.084, 0.216, 0.01, 0.10571428571428572, 0.084, 0.216, 0.04, 0.1, 0.216, 0.084}
	got := BenjaminiHochberg(in)
	if len(got) != len(want) {
		t.Fatalf("BenjaminiHochberg returned %d values, want %d", len(got), len(want))
	}
	for i := range want {
		if !near(got[i], want[i]) {
			t.Fatalf("BenjaminiHochberg = %v, want %v", got, want)
		}
	}
	if in[0] != 0.042 {
		t.Error("BenjaminiHochberg modified its input")
	}
}

func TestCohenD(t *testing.T) {
	if got := CohenD(designA, designB); !near(got, 1.8438942692909126) {
		t.Errorf("CohenD(designA, designB) = %v, want 1.8438942692909126", got)
	}
	if got := CohenD(designB, designA); !near(got, -1.8438942692909126) {
		t.Errorf("CohenD(designB, designA) = %v, want -1.8438942692909126", got)
	}
}

//go:build solutions

package hypothesis

import (
	"math"
	"slices"
)

// ============================================================
// Part A — t tests
// ============================================================

// twoSidedT is the two-sided p-value for a t statistic: 2·P(T <= −|t|).
// Computing the lower tail directly keeps tiny p-values accurate; 1 − tCDF
// would round them to 0.
func twoSidedT(t, df float64) float64 {
	return 2 * tCDF(-math.Abs(t), df)
}

func OneSampleT(xs []float64, mu0 float64) Result {
	n := float64(len(xs))
	t := (mean(xs) - mu0) / math.Sqrt(variance(xs)/n)
	return Result{Statistic: t, DF: n - 1, P: twoSidedT(t, n-1)}
}

func WelchT(a, b []float64) Result {
	na, nb := float64(len(a)), float64(len(b))
	va, vb := variance(a)/na, variance(b)/nb
	t := (mean(a) - mean(b)) / math.Sqrt(va+vb)
	df := (va + vb) * (va + vb) / (va*va/(na-1) + vb*vb/(nb-1))
	return Result{Statistic: t, DF: df, P: twoSidedT(t, df)}
}

func PairedT(a, b []float64) Result {
	d := make([]float64, len(a))
	for i := range a {
		d[i] = a[i] - b[i]
	}
	return OneSampleT(d, 0)
}

// ============================================================
// Part B — chi-square tests for counts
// ============================================================

func ChiSquareGOF(observed []int, expected []float64) Result {
	chi := 0.0
	for i, o := range observed {
		d := float64(o) - expected[i]
		chi += d * d / expected[i]
	}
	df := float64(len(observed) - 1)
	return Result{Statistic: chi, DF: df, P: 1 - chiSquareCDF(chi, df)}
}

func ChiSquareIndependence(table [][]int) Result {
	rows := make([]float64, len(table))
	cols := make([]float64, len(table[0]))
	total := 0.0
	for i, row := range table {
		for j, c := range row {
			rows[i] += float64(c)
			cols[j] += float64(c)
			total += float64(c)
		}
	}
	chi := 0.0
	for i, row := range table {
		for j, c := range row {
			e := rows[i] * cols[j] / total
			chi += (float64(c) - e) * (float64(c) - e) / e
		}
	}
	df := float64((len(rows) - 1) * (len(cols) - 1))
	return Result{Statistic: chi, DF: df, P: 1 - chiSquareCDF(chi, df)}
}

// ============================================================
// Part C — many tests, and how big the effect is
// ============================================================

func Bonferroni(ps []float64) []float64 {
	m := float64(len(ps))
	out := make([]float64, len(ps))
	for i, p := range ps {
		out[i] = math.Min(1, p*m)
	}
	return out
}

func BenjaminiHochberg(ps []float64) []float64 {
	m := len(ps)
	idx := make([]int, m)
	for i := range idx {
		idx[i] = i
	}
	slices.SortFunc(idx, func(a, b int) int {
		switch {
		case ps[a] < ps[b]:
			return -1
		case ps[a] > ps[b]:
			return 1
		}
		return 0
	})
	out := make([]float64, m)
	running := 1.0
	for k := m; k >= 1; k-- {
		i := idx[k-1]
		running = math.Min(running, ps[i]*float64(m)/float64(k))
		out[i] = running
	}
	return out
}

func CohenD(a, b []float64) float64 {
	na, nb := float64(len(a)), float64(len(b))
	pooled := math.Sqrt(((na-1)*variance(a) + (nb-1)*variance(b)) / (na + nb - 2))
	return (mean(a) - mean(b)) / pooled
}

• Add PooledT(a, b) — Student's equal-variance t-test, df = na + nb − 2 — and compare its p-value with Welch's when one group is much more variable.
• Simulate power: draw two normal samples of size n with a true difference of 0.5 sd, run WelchT 2,000 times, and count rejections. How big must n be for 80% power?
• Practice: run a goodness-of-fit test on the last digits of numbers in a spreadsheet you trust and one you don't.

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