Sampling

You almost never get the whole population. You get a sample, compute a number from it, and want to know how far that number could be from the truth. The answer comes from the sampling distribution: how the statistic would vary if you could repeat the sampling many times. This module builds that idea by brute force — draw thousands of samples, look at the results — so that the formulas in Modules 4 and 5 stop being magic.

Two results come out of it. The central limit theorem says sample means are approximately normal whatever the data looks like. And the bootstrap lets you estimate the variability of almost any statistic from the one sample you have.

A statistic is a random variable. Draw a different sample and you get a different mean. The distribution of those means across all possible samples is the sampling distribution of the mean.

The standard error. The sampling distribution of the mean is centred on the population mean, and its standard deviation — called the standard error — is σ / √n. Quadrupling the sample size halves the error. That square root is why polls of 1,000 people are useful and polls of 10,000 aren't ten times better.

The central limit theorem. For a large enough n, the sample mean is approximately normal no matter how skewed or lumpy the population is. "Large enough" is often 30 for moderately skewed data and more for wildly skewed data. This is why the normal curve from Module 2 shows up in nearly every test and interval.

The bootstrap. You can't resample the population, but you can resample your sample: draw n values with replacement from your n data points, compute the statistic, and repeat thousands of times. The spread of those results estimates the statistic's sampling variability. It works for medians, correlations, and other statistics that have no tidy formula.

Permutation tests. To ask whether two groups differ, assume they don't: if group labels are meaningless, shuffling them shouldn't matter. Shuffle thousands of times and see how often a difference as large as the real one appears by chance. That fraction is a p-value, computed without any distribution tables.

• Use math/rand/v2. rand.New(rand.NewPCG(seed1, seed2)) gives a generator you control; r.IntN(n) draws an int in [0, n), r.Float64() one in [0, 1), and r.Shuffle(n, swap) shuffles in place.
• Seed everything you simulate. With a fixed seed, the same program prints the same numbers every run, which makes results reproducible and tests possible. The package-level functions (rand.IntN) are randomly seeded — fine for games, wrong for analysis.
• Pass the generator in. Functions that take a *rand.Rand parameter, rather than reaching for a global, can be tested and reproduced, and two simulations can't interfere with each other.
• for range 10000 (Go 1.22+) loops a fixed number of times without an unused index variable.
• Functions are values. Bootstrap(r, data, 5000, mean) passes the mean function itself; any func([]float64) float64 works, including a closure.
• slices.Concat(a, b) builds a new pooled slice, so shuffling it never disturbs the caller's data.

package main

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

// histogram prints one bar per bin; each # is 100 of the 10,000 results.
func histogram(title string, bins map[float64]int, lo, hi, step float64) {
	fmt.Println(title)
	for x := lo; x <= hi+1e-9; x += step {
		n := bins[math.Round(x*10)/10]
		fmt.Printf("  %4.1f %-36s %d\n", x, strings.Repeat("#", n/100), n)
	}
}

func main() {
	r := rand.New(rand.NewPCG(1, 2)) // fixed seed: the same output every run
	const reps = 10000

	one := map[float64]int{} // a single die: flat
	for range reps {
		one[float64(r.IntN(6)+1)]++
	}

	avg := map[float64]int{} // the average of 10 dice, rounded to the nearest 0.5
	for range reps {
		sum := 0
		for range 10 {
			sum += r.IntN(6) + 1
		}
		avg[math.Round(float64(sum)/10*2)/2]++
	}

	histogram("one die:", one, 1, 6, 1)
	histogram("average of 10 dice:", avg, 1, 6, 0.5)
}

It prints:

A single die is flat — every face equally likely. The average of ten dice is already a bell, centred on 3.5 and far narrower: the standard deviation of one die is 1.71, and of a 10-die average 1.71 / √10 = 0.54. Nothing about dice is normal; averaging made it so. (Other seeds give slightly different counts; the shape doesn't change.)

• Confusing the standard deviation with the standard error. The sd describes individual values; the SE describes the mean. A sample of 10,000 has a tiny SE but the same sd as a sample of 100.
• Bootstrapping without replacement. Drawing n of n values without replacement just reorders the data — every resample has the same mean.
• Expecting the bootstrap to fix a bad sample. It estimates the variability of your statistic around your sample's value. If the sample is biased (only customers who answered a survey), resampling it faithfully reproduces the bias.
• Reporting p = 0 from a simulation. Ten thousand shuffles with no extreme result means p is below about 1/10,000, not zero. The +1 in (count + 1) / (reps + 1) keeps that honest.

Part A simulates sampling distributions and checks the standard error formula against them. Part B is the bootstrap: resampling, applying any statistic, and a percentile interval. Part C is a permutation test for a difference in means.

go test ./03-sampling/

The tests check statistical properties with generous tolerances — the mean of 4,000 sample means should be close to the population mean, for example — so any correct implementation passes, whatever order it draws its random numbers in.

The files for this module

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

// Package sampling is Module 3: sampling distributions, the central limit
// theorem, and the bootstrap.
//
// Everything here is simulation, so every function takes a *rand.Rand from
// math/rand/v2. Build it from a fixed seed with NewRand and a run repeats
// exactly: same seed, same numbers, same answer.
package sampling

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

// NewRand returns a random-number generator seeded with seed. The same seed
// always produces the same stream of numbers.
func NewRand(seed uint64) *rand.Rand {
	return rand.New(rand.NewPCG(seed, seed^0x9e3779b97f4a7c15))
}

// The helpers below are given — you wrote them in Module 1.

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

// quantile returns the q-quantile of xs by linear interpolation.
func quantile(xs []float64, q float64) float64 {
	s := slices.Clone(xs)
	slices.Sort(s)
	h := float64(len(s)-1) * q
	lo := int(math.Floor(h))
	if lo >= len(s)-1 {
		return s[len(s)-1]
	}
	return s[lo] + (h-float64(lo))*(s[lo+1]-s[lo])
}

//go:build !solutions

package sampling

import "math/rand/v2"

// ============================================================
// Part A — sampling distributions
// ============================================================

// Sample draws n values from pop at random, with replacement: each draw is
// pop[r.IntN(len(pop))]. Target: O(n).
func Sample(r *rand.Rand, pop []float64, n int) []float64 {
	// TODO
	return nil
}

// SampleMeans draws reps samples of size n from pop and returns each
// sample's mean — a simulated sampling distribution of the mean. By the
// central limit theorem it looks normal for large n, whatever pop looks
// like. Target: O(reps·n).
func SampleMeans(r *rand.Rand, pop []float64, n, reps int) []float64 {
	// TODO
	return nil
}

// StandardError returns the standard deviation of a sample mean, sd / √n,
// where sd is the standard deviation of the individual values.
func StandardError(sd float64, n int) float64 {
	// TODO
	return 0
}

// ============================================================
// Part B — the bootstrap
// ============================================================

// Resample returns a bootstrap sample: len(data) draws from data with
// replacement. Some values appear twice, some not at all.
func Resample(r *rand.Rand, data []float64) []float64 {
	// TODO
	return nil
}

// Bootstrap resamples data reps times and returns stat of each resample —
// an estimate of stat's sampling distribution, using only the one sample
// you have. stat can be any function: mean, median, a correlation…
func Bootstrap(r *rand.Rand, data []float64, reps int, stat func([]float64) float64) []float64 {
	// TODO
	return nil
}

// PercentileInterval returns the middle fraction level of boot: its
// (1−level)/2 and (1+level)/2 quantiles (use the given quantile helper).
// With level 0.95 that's the 2.5th and 97.5th percentiles — a 95% bootstrap
// confidence interval.
func PercentileInterval(boot []float64, level float64) (lo, hi float64) {
	// TODO
	return 0, 0
}

// ============================================================
// Part C — a permutation test
// ============================================================

// PermutationTest asks whether groups a and b could plausibly come from the
// same population. The observed statistic is |mean(a) − mean(b)|. Pool both
// groups, then reps times: shuffle the pool (r.Shuffle), call the first
// len(a) values group a and the rest group b, and count how often the
// shuffled |difference| is at least the observed one. Return the p-value
// (count + 1) / (reps + 1) — the +1s count the observed split itself, so
// the p-value is never exactly 0.
func PermutationTest(r *rand.Rand, a, b []float64, reps int) float64 {
	// TODO
	return 0
}

package sampling

import (
	"math"
	"slices"
	"testing"
)

// These tests check statistical properties with generous tolerances, so any
// correct implementation passes whatever order it draws random numbers in.

var die = []float64{1, 2, 3, 4, 5, 6}

// skewed is a lopsided population: mean 4.4, standard deviation (÷N) 5.6249.
var skewed = []float64{1, 1, 1, 1, 2, 2, 3, 5, 8, 20}

func TestSample(t *testing.T) {
	s := Sample(NewRand(1), die, 60000)
	if len(s) != 60000 {
		t.Fatalf("Sample returned %d values, want 60000", len(s))
	}
	counts := map[float64]int{}
	for _, x := range s {
		counts[x]++
	}
	for _, face := range die {
		if f := float64(counts[face]) / 60000; math.Abs(f-1.0/6) > 0.01 {
			t.Errorf("face %v drawn %.4f of the time, want about 1/6", face, f)
		}
	}
	if len(counts) != 6 {
		t.Errorf("Sample drew values outside the population: %v", counts)
	}
	if !slices.Equal(Sample(NewRand(7), die, 20), Sample(NewRand(7), die, 20)) {
		t.Error("the same seed should give the same sample")
	}
}

func TestSampleMeans(t *testing.T) {
	means := SampleMeans(NewRand(2), skewed, 30, 4000)
	if len(means) != 4000 {
		t.Fatalf("SampleMeans returned %d values, want 4000", len(means))
	}
	if m := mean(means); math.Abs(m-4.4) > 0.1 {
		t.Errorf("mean of the sample means = %.3f, want about 4.4 (the population mean)", m)
	}
	// the spread of the means is the standard error: 5.6249 / √30 = 1.027
	if sd := stdDev(means); math.Abs(sd-1.027)/1.027 > 0.06 {
		t.Errorf("sd of the sample means = %.3f, want about 1.027", sd)
	}
}

func TestStandardError(t *testing.T) {
	if got := StandardError(5.624944444170093, 30); math.Abs(got-1.0269696522617728) > 1e-12 {
		t.Errorf("StandardError(5.6249, 30) = %v, want 1.0269696522617728", got)
	}
	if got := StandardError(12, 9); got != 4 {
		t.Errorf("StandardError(12, 9) = %v, want 4", got)
	}
}

func TestResample(t *testing.T) {
	data := []float64{3, 1, 4, 1, 5, 9, 2, 6}
	r := NewRand(3)
	sawRepeat := false
	for range 50 {
		s := Resample(r, data)
		if len(s) != len(data) {
			t.Fatalf("Resample returned %d values, want %d", len(s), len(data))
		}
		for _, x := range s {
			if !slices.Contains(data, x) {
				t.Fatalf("Resample produced %v, which isn't in the data", x)
			}
		}
		c := slices.Clone(s)
		slices.Sort(c)
		if len(slices.Compact(c)) < 7 { // fewer distinct values than the data has
			sawRepeat = true
		}
	}
	if !sawRepeat {
		t.Error("50 resamples never repeated a value — is Resample drawing with replacement?")
	}
}

func TestBootstrap(t *testing.T) {
	data := []float64{12.1, 9.8, 11.4, 10.2, 13.7, 8.9, 10.6, 12.8, 11.1, 9.5}
	boot := Bootstrap(NewRand(4), data, 3000, mean)
	if len(boot) != 3000 {
		t.Fatalf("Bootstrap returned %d values, want 3000", len(boot))
	}
	if m := mean(boot); math.Abs(m-mean(data)) > 0.05 {
		t.Errorf("mean of bootstrap means = %.3f, want about %.3f", m, mean(data))
	}
	// the bootstrap sd of the mean approximates s/√n = 0.4466 (a bit less: ÷n, not ÷(n−1))
	if sd := stdDev(boot); sd < 0.36 || sd > 0.50 {
		t.Errorf("bootstrap sd of the mean = %.3f, want about 0.42", sd)
	}
	counted := 0
	Bootstrap(NewRand(5), data, 7, func(xs []float64) float64 { counted++; return 0 })
	if counted != 7 {
		t.Errorf("stat was called %d times for 7 reps", counted)
	}
}

func TestPercentileInterval(t *testing.T) {
	boot := make([]float64, 101)
	for i := range boot {
		boot[100-i] = float64(i) // 0..100, in reverse order
	}
	if lo, hi := PercentileInterval(boot, 0.90); math.Abs(lo-5) > 1e-9 || math.Abs(hi-95) > 1e-9 {
		t.Errorf("PercentileInterval(0..100, 0.90) = (%v, %v), want (5, 95)", lo, hi)
	}
	if lo, hi := PercentileInterval(boot, 0.95); math.Abs(lo-2.5) > 1e-9 || math.Abs(hi-97.5) > 1e-9 {
		t.Errorf("PercentileInterval(0..100, 0.95) = (%v, %v), want (2.5, 97.5)", lo, hi)
	}
}

func TestPermutationTest(t *testing.T) {
	a := []float64{5.1, 4.9, 5.3, 5.0, 4.8, 5.2}
	b := []float64{5.0, 5.2, 4.9, 5.1, 5.3, 4.8} // the same values, reordered
	if p := PermutationTest(NewRand(6), a, b, 2000); p < 0.5 {
		t.Errorf("identical groups: p = %.4f, want a large p-value", p)
	}
	low := []float64{1, 2, 3, 4, 5, 6, 7, 8}
	high := []float64{11, 12, 13, 14, 15, 16, 17, 18}
	p := PermutationTest(NewRand(8), low, high, 5000)
	if p <= 0 || p > 0.01 {
		t.Errorf("completely separated groups: p = %.5f, want 0 < p <= 0.01", p)
	}
	if len(low) != 8 || low[0] != 1 || high[0] != 11 {
		t.Error("PermutationTest modified its inputs — shuffle a pooled copy")
	}
}

//go:build solutions

package sampling

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

// ============================================================
// Part A — sampling distributions
// ============================================================

func Sample(r *rand.Rand, pop []float64, n int) []float64 {
	out := make([]float64, n)
	for i := range out {
		out[i] = pop[r.IntN(len(pop))]
	}
	return out
}

func SampleMeans(r *rand.Rand, pop []float64, n, reps int) []float64 {
	out := make([]float64, reps)
	for i := range out {
		out[i] = mean(Sample(r, pop, n))
	}
	return out
}

func StandardError(sd float64, n int) float64 {
	return sd / math.Sqrt(float64(n))
}

// ============================================================
// Part B — the bootstrap
// ============================================================

func Resample(r *rand.Rand, data []float64) []float64 {
	return Sample(r, data, len(data))
}

func Bootstrap(r *rand.Rand, data []float64, reps int, stat func([]float64) float64) []float64 {
	out := make([]float64, reps)
	for i := range out {
		out[i] = stat(Resample(r, data))
	}
	return out
}

func PercentileInterval(boot []float64, level float64) (lo, hi float64) {
	return quantile(boot, (1-level)/2), quantile(boot, (1+level)/2)
}

// ============================================================
// Part C — a permutation test
// ============================================================

func PermutationTest(r *rand.Rand, a, b []float64, reps int) float64 {
	observed := math.Abs(mean(a) - mean(b))
	pool := slices.Concat(a, b)
	count := 0
	for range reps {
		r.Shuffle(len(pool), func(i, j int) { pool[i], pool[j] = pool[j], pool[i] })
		if math.Abs(mean(pool[:len(a)])-mean(pool[len(a):])) >= observed {
			count++
		}
	}
	return float64(count+1) / float64(reps+1)
}

• Bootstrap the median of a skewed sample and compare the interval's width with the mean's.
• Make the CLT fail: sample from a population with one enormous value (say 10,000 ones and a single 1,000,000) and watch how large n has to get before the means look normal.
• Practice: shuffle-test whether your weekday and weekend commute times really differ.

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