Level 4 · Statistics & Probability: Randomness and Sampling

Level 2 computed probabilities of single events. This level gives randomness a structure. A random variable is a number whose value depends on chance (the score on a die, the number of emails in an hour, the height of a person picked at random), and its distribution says how likely each value is. A few distributions turn up again and again because the situations that produce them are common: the Bernoulli and binomial for yes/no trials, the Poisson for counts of rare events, and the normal for anything built from many small effects. You will compute each one from its formula and then simulate it, so the formula and the behaviour always sit side by side.

The second half of the level turns the question around. In practice you do not know the distribution; you have one sample from it. How much would your mean have changed with a different sample? The law of large numbers, sampling distributions and the standard error answer that, and two simulation methods, Monte Carlo and the bootstrap, let a computer answer it for almost any quantity. These ideas carry all of Level 5: a confidence interval is a statement about a sampling distribution.

Links to the other paths: the expected value is the "average loss" that the Machine Learning path minimises, and Monte Carlo simulation shares an idea with the randomised heuristics of the Algorithms path (simulated annealing, genetic algorithms): give up certainty in exchange for a method that works when exact computation is out of reach.

Random variables and expected value

A game costs £1 to play. You roll a die; on 4 or 5 you get £1 back, on 6 you get £3, otherwise nothing. The amount you get back is a random variable X with a small distribution: P(X = 0) = 3/6, P(X = 1) = 2/6, P(X = 3) = 1/6.

The expected value E[X] is the probability-weighted average of the values:

The sum runs over the possible values x of X. With the distribution stored as a dictionary from value to probability, it becomes sum(x * p for x, p in dist.items()), the first line of the code below.

It is the mean you would see over a very long run of plays. Compute it exactly, then play 100,000 times:

import random
from fractions import Fraction

dist = {0: Fraction(3, 6), 1: Fraction(2, 6), 3: Fraction(1, 6)}
ev = sum(x * p for x, p in dist.items())
print("E[X] =", ev, " expected profit per game =", ev - 1)

payout = {1: 0, 2: 0, 3: 0, 4: 1, 5: 1, 6: 3}
rng = random.Random(2)
n = 100000
profit = sum(payout[rng.randint(1, 6)] - 1 for _ in range(n))
print("simulated profit per game:", round(profit / n, 4))

E[X] = 5/6  expected profit per game = -1/6
simulated profit per game: -0.1656

On average you lose about 17p a game. No single game loses 17p (you lose £1, break even or win £2), which is the same point as in Level 1: the expected value describes the long run, not a typical outcome.

The variance of a random variable is the expected squared distance from its mean \mu = E[X], and its square root is the standard deviation. They are the population versions of the formulas in Level 2, with probabilities as weights in place of \frac{1}{n}:

from fractions import Fraction
import math

dist = {0: Fraction(3, 6), 1: Fraction(2, 6), 3: Fraction(1, 6)}
mu = sum(x * p for x, p in dist.items())
var = sum((x - mu) ** 2 * p for x, p in dist.items())
print("mean", mu, "variance", var, "sd", round(math.sqrt(var), 3))

mean 5/6 variance 41/36 sd 1.067

One rule makes expected values easy to work with: linearity. The expected value of a sum is the sum of the expected values, E[X + Y] = E[X] + E[Y], always, even when X and Y are dependent. Ten games lose on average ten times 1/6 of a pound. Variances add too, but only when the variables are independent.

Common mistakes: forgetting that the probabilities must add up to 1 (check it in code). Averaging the values without weighting them by their probabilities. Adding variances of dependent quantities. And reading a positive expected value as "safe": a game with a high average can still ruin you if a single loss is larger than you can afford.

Practice: What Is a Scratch Card Worth?, How Many Rainy Spells?

Bernoulli and binomial

A Bernoulli random variable is a single yes/no trial: 1 with probability p, 0 otherwise. Its mean is p and its variance p(1 - p). The binomial distribution counts the successes in n independent Bernoulli trials with the same p: the number of heads in 10 tosses, the number of free throws scored out of 10.

The probability of exactly k successes is

where \binom{n}{k} (read "n choose k", also written C(n, k), and computed by math.comb(n, k)) counts the ways to choose which k of the n trials succeed, and each such arrangement has probability p^k (1 - p)^{n - k}. A player scores 70% of free throws; how many will she score out of 10?

import math
import random

n, p = 10, 0.7
pmf = [math.comb(n, k) * p ** k * (1 - p) ** (n - k) for k in range(n + 1)]

rng = random.Random(4)
trials = 10000
counts = [0] * (n + 1)
for _ in range(trials):
    k = sum(rng.random() < p for _ in range(n))
    counts[k] += 1

for k in range(4, n + 1):
    print(f"{k:2d}  exact {pmf[k]:.3f}  simulated {counts[k] / trials:.3f}")
print("mean:", round(sum(k * pmf[k] for k in range(n + 1)), 6), " n*p =", n * p)

 4  exact 0.037  simulated 0.038
 5  exact 0.103  simulated 0.103
 6  exact 0.200  simulated 0.206
 7  exact 0.267  simulated 0.261
 8  exact 0.233  simulated 0.232
 9  exact 0.121  simulated 0.121
10  exact 0.028  simulated 0.029
mean: 7.0  n*p = 7.0

(Values below 4 have tiny probabilities and are left out of the printout.) The simulation matches the formula to about two decimal places, as the 1/sqrt(n) rule from Level 2 predicts for 10,000 trials. The most likely count is 7, but it happens only about a quarter of the time: a 70% shooter scores 6 or fewer out of 10 in about 35% of games. When a commentator says a player "has gone cold" after scoring 5 of 10, remember that the binomial gives this about one game in ten.

The mean of a binomial is np (linearity: n Bernoullis with mean p each) and its variance is np(1 - p) (independent trials, so the variances add).

Common mistakes: using the binomial when the trials are not independent (drawing cards without replacement changes p after every draw) or when p changes between trials. Forgetting the C(n, k) factor and giving the probability of one particular sequence. And writing p ** k * (1 - p) ** n - k, which subtracts k after the power; the brackets around n - k matter.

Practice: Accept the Batch?, How Many Tickets for the Ferry?

Poisson: counting rare events

Some counts have no fixed number of trials: emails arriving in an hour, typing errors on a page, goals in a football match. Events happen independently at a steady average rate, and you count how many fall in a fixed window. If the average count is λ (lambda), the number of events follows the Poisson distribution:

Where does it come from? Split the hour into 1000 tiny slots, each with a small chance λ/1000 of holding an email. The count is then binomial with a huge n and a tiny p, and as the slots get finer the binomial turns into the Poisson:

import math

lam = 3
for k in range(8):
    poisson = math.exp(-lam) * lam ** k / math.factorial(k)
    binom = math.comb(1000, k) * 0.003 ** k * 0.997 ** (1000 - k)
    print(k, round(poisson, 4), round(binom, 4))

0 0.0498 0.0496
1 0.1494 0.1491
2 0.224 0.2242
3 0.224 0.2244
4 0.168 0.1683
5 0.1008 0.1009
6 0.0504 0.0503
7 0.0216 0.0215

The two columns agree to about three decimal places. You can also simulate the arrivals themselves. When events happen at rate λ, the gaps between them follow an exponential distribution, which random.expovariate(lam) samples. Count how many arrivals fit into one hour, for many hours:

import math
import random

rng = random.Random(9)
lam, hours = 3, 20000
counts = []
for _ in range(hours):
    t, k = rng.expovariate(lam), 0
    while t < 1:
        k += 1
        t += rng.expovariate(lam)
    counts.append(k)

m = sum(counts) / hours
v = sum((c - m) ** 2 for c in counts) / hours
print("mean", round(m, 3), "variance", round(v, 3))
print("P(0 emails)", counts.count(0) / hours, " exact", round(math.exp(-3), 4))

mean 3.018 variance 2.981
P(0 emails) 0.04775  exact 0.0498

A signature of the Poisson is that its mean equals its variance (both λ). That gives a quick check on real count data: if the variance is much larger than the mean, the events are clumping together (emails arriving in bursts), and a Poisson model will underestimate how often extreme counts occur.

Common mistakes: using the Poisson when events are not independent (a busy shop gets queues, a sick classroom spreads a cold). Using a rate for the wrong window: 3 emails an hour means λ = 1.5 for half an hour. And computing lam ** k / math.factorial(k) for very large k in floats, which can overflow; for big counts, work with logarithms (Level 6 does this for likelihoods).

Practice: Messages at the Help Desk, Bursty or Steady?

The normal distribution

Heights, measurement errors and exam totals often have the same bell-shaped histogram: most values near the mean, fewer and fewer further out, and symmetric on both sides. That shape is the normal (or Gaussian) distribution. It has two parameters, the mean μ and the standard deviation σ, and its density is

A density is not a probability: probabilities are areas under it. The area within k standard deviations of the mean is the same for every normal distribution, and math.erf computes it:

import math
for k in (1, 2, 3):
    print(k, round(math.erf(k / math.sqrt(2)), 4))

1 0.6827
2 0.9545
3 0.9973

This is the 68–95–99.7 rule: about 68% of values lie within one standard deviation of the mean, 95% within two and 99.7% within three. It is where the rule of thumb in Level 2 came from. Check it on simulated heights with mean 170 cm and standard deviation 8 cm, drawn with random.gauss:

import random

rng = random.Random(12)
heights = [rng.gauss(170, 8) for _ in range(50000)]
for k in (1, 2, 3):
    inside = sum(abs(h - 170) < k * 8 for h in heights) / len(heights)
    print(k, round(inside, 4))

from statistics import NormalDist
print("P(height > 186) =", round(1 - NormalDist(170, 8).cdf(186), 4))
print("simulated:       ", round(sum(h > 186 for h in heights) / len(heights), 4))

1 0.6832
2 0.9541
3 0.9973
P(height > 186) = 0.0228
simulated:        0.0229

186 cm is exactly two standard deviations above the mean, so about 2.3% of people are taller (half of the 4.6% outside two standard deviations). statistics.NormalDist gives the areas directly, and it is the right tool once you understand what it computes.

The density formula turns into one line of Python with math.exp and math.sqrt. Press Run to lay it over a histogram of the simulated heights. Each bar is a 2 cm slice, so its share of the people should be close to 2 f(x) at the slice's centre:

import math
import random

rng = random.Random(12)
counts = [0] * 30
for _ in range(50000):
    b = int((rng.gauss(170, 8) - 140) // 2)
    if 0 <= b < 30:
        counts[b] += 1
centres = [141 + 2 * b for b in range(30)]
plot(centres, [c / 50000 for c in counts], kind="bar", label="simulated", xlabel="height (cm)", ylabel="share")
plot(centres, [2 * math.exp(-(x - 170) ** 2 / (2 * 8 ** 2)) / (8 * math.sqrt(2 * math.pi)) for x in centres], label="2 f(x)")

Why is the normal so common? Because a quantity built by adding up many small independent effects is approximately normal, whatever the effects look like individually. That is the central limit theorem, which you will see happen in the section on sampling distributions. It also explains when the normal fails: incomes and waiting times are shaped by multiplying effects or by a hard floor at zero, and they are skewed, not bell-shaped.

Common mistakes: assuming data is normal without looking at a histogram. Reading the density f(x) as the probability of exactly x (for a continuous variable that probability is 0; ask about intervals). Using the 68–95–99.7 rule for skewed data. And passing the variance where the standard deviation is expected: random.gauss(170, 64) gives people 170 ± 64 cm.

Practice: Ordering Jackets by Height, When Does the Bell Curve Fit the Bars?

The law of large numbers

Every simulation so far has relied on one fact: the average of many independent draws gets close to the expected value, and closer as the number of draws grows. That is the law of large numbers. Watch it for a fair coin:

import random

rng = random.Random(21)
heads = 0
checkpoints = {10, 100, 1000, 10000, 100000}
for i in range(1, 100001):
    heads += rng.random() < 0.5
    if i in checkpoints:
        print(f"{i:>6} tosses: share {heads / i:.4f}, heads minus half {heads - i / 2:+.0f}")

    10 tosses: share 0.4000, heads minus half -1
   100 tosses: share 0.5100, heads minus half +1
  1000 tosses: share 0.4990, heads minus half -1
 10000 tosses: share 0.5042, heads minus half +42
100000 tosses: share 0.4986, heads minus half -143

The share of heads closes in on 0.5. But look at the last column: the number of heads minus half the tosses does not shrink at all; after 100,000 tosses it is 143 away. The law of large numbers does not say that the coin corrects its past. It says that a surplus of a few dozen or a few hundred heads becomes a tiny fraction of a huge total. Past imbalances are diluted, not reversed.

The share of heads is a running mean of 0/1 values, the recurrence \bar{x}_k = \bar{x}_{k-1} + (x_k - \bar{x}_{k-1})/k from Level 1. Press Run to watch it settle (the horizontal axis is logarithmic, so each step to the right is ten times more tosses):

import random

rng = random.Random(21)
heads, ns, shares = 0, [], []
for i in range(1, 100001):
    heads += rng.random() < 0.5
    if i < 50 or i % 50 == 0:
        ns.append(i)
        shares.append(heads / i)
plot(ns, shares, logx=True, label="share of heads", xlabel="tosses", ylabel="share")
plot([1, 100000], [0.5, 0.5], label="0.5")

That is exactly why the gambler's fallacy is a fallacy. After five reds in a row at roulette, black is not "due": the wheel has no memory, and the next spin is as likely to be red as ever. You can check that in a simulation by looking only at the tosses that follow five heads in a row:

import random

rng = random.Random(22)
tosses = [rng.random() < 0.5 for _ in range(200000)]
after_streak = [tosses[i] for i in range(5, len(tosses)) if all(tosses[i - 5:i])]
print(len(after_streak), "tosses after five heads; share of heads:",
      round(sum(after_streak) / len(after_streak), 3))

6164 tosses after five heads; share of heads: 0.491

The rate at which the average settles down is the 1 / sqrt(n) rule you have been using. The standard deviation of the share of heads in n tosses is \sqrt{p(1 - p)/n}, which is 0.05 for 100 tosses and 0.0016 for 100,000. The next section explains where that comes from.

Common mistakes: expecting the count (rather than the proportion) to converge. Believing that a streak makes the opposite outcome more likely. And applying the law to a handful of observations: it is a statement about the long run, and says little about ten tosses.

Practice: Is a Six Due?

Sampling distributions and the standard error

A sample is a handful of values drawn from a much larger population, and a statistic (the sample mean, the sample median, a proportion) is computed from it. Draw a different sample and you get a different value. The distribution of the statistic over all the samples you could have drawn is its sampling distribution, and its standard deviation is the standard error. The standard error is the answer to "how much would my estimate wobble if I repeated the study?".

You can only see a sampling distribution by simulation (or with mathematics). Take a skewed population, waiting times with mean 4 minutes and standard deviation 4 (exponential), draw 5000 samples of size n, and look at the 5000 sample means:

import math
import random

rng = random.Random(30)
for n in (1, 4, 16, 64):
    means = []
    for _ in range(5000):
        sample = [rng.expovariate(1 / 4) for _ in range(n)]
        means.append(sum(sample) / n)
    m = sum(means) / len(means)
    se = math.sqrt(sum((x - m) ** 2 for x in means) / len(means))
    print(f"n={n:2d}  mean of means {m:.2f}  sd of means {se:.3f}  sigma/sqrt(n) {4 / math.sqrt(n):.3f}")

n= 1  mean of means 4.01  sd of means 4.150  sigma/sqrt(n) 4.000
n= 4  mean of means 3.96  sd of means 2.002  sigma/sqrt(n) 2.000
n=16  mean of means 3.99  sd of means 0.985  sigma/sqrt(n) 1.000
n=64  mean of means 4.00  sd of means 0.505  sigma/sqrt(n) 0.500

Two facts appear. The sample means are centred on the true mean 4 (the sample mean is unbiased). And their spread is

where \sigma is the standard deviation of the population and n the sample size: quadruple the sample size, halve the standard error. That is the formula behind the 1 / sqrt(n) rule, and it is why precision is expensive: one more decimal place of accuracy costs a hundred times more data.

The third fact is about shape. Here are histograms of single waiting times and of means of 16:

import random

rng = random.Random(31)

def histogram(values, top):
    counts = {}
    for v in values:
        b = min(int(v), top)
        counts[b] = counts.get(b, 0) + 1
    for b in range(top + 1):
        label = f"{b}+" if b == top else str(b)
        print(f"{label:>4}  {'#' * (counts.get(b, 0) // 100)}".rstrip())

for n in (1, 16):
    means = [sum(rng.expovariate(1 / 4) for _ in range(n)) / n for _ in range(5000)]
    print("means of", n)
    histogram(means, 10)

means of 1
   0  ##########
   1  ########
   2  ######
   3  #####
   4  ###
   5  ###
   6  ##
   7  #
   8  #
   9  #
 10+  ####
means of 16
   0
   1
   2  #######
   3  ##################
   4  ###############
   5  ######
   6  #
   7
   8
   9
 10+

Press Run to draw the two histograms as charts, with half-minute bins and each bar showing a share of the 5000 values (values above 12 minutes are left off):

import random

rng = random.Random(31)
centres = [b / 2 + 0.25 for b in range(24)]
for n in (1, 16):
    means = [sum(rng.expovariate(1 / 4) for _ in range(n)) / n for _ in range(5000)]
    shares = [sum(c - 0.25 <= m < c + 0.25 for m in means) / 5000 for c in centres]
    plot(centres, shares, kind="bar", new=True, title=f"means of {n}", xlabel="minutes", ylabel="share")

Single waiting times are strongly skewed: most are short, with a long tail to the right. Means of 16 are bunched around 4 and nearly symmetric, already close to a bell shape. This is the central limit theorem: means (and sums) of many independent values are approximately normal, whatever the shape of the population, and the approximation improves as n grows. It is why the normal distribution is used for inference about means in Level 5 even when the data itself is far from normal.

Common mistakes: confusing the standard deviation of the data (how much individuals vary, which does not shrink with more data) with the standard error of the mean (how much the estimate varies, which does). Assuming the central limit theorem kicks in for tiny samples from very skewed populations; with n = 4 the means of these waiting times are still clearly skewed (change 16 to 4 in the histogram code to see it). And forgetting that σ / sqrt(n) needs independent draws: surveying 100 people from ten households is not 100 independent observations.

Practice: A Thousand Small Surveys, Asking Half the Village, How Sure Is the Average Temperature?

Monte Carlo simulation

Monte Carlo methods answer a question about a probability or an average by simulating it many times, as you have been doing since Level 2. What this level adds is the ability to say how accurate the answer is: a Monte Carlo estimate is a sample mean, so its standard error is σ / sqrt(n).

The classic first example estimates π. Throw random points into the unit square; the fraction that lands inside the quarter circle x² + y² < 1 estimates its area, π/4:

import math
import random

rng = random.Random(2)
inside = 0
for i in range(1, 100001):
    x, y = rng.random(), rng.random()
    if x * x + y * y < 1:
        inside += 1
    if i in (100, 1000, 10000, 100000):
        est = 4 * inside / i
        se = 4 * math.sqrt((math.pi / 4) * (1 - math.pi / 4) / i)
        print(f"{i:>6}  estimate {est:.4f}  error {abs(est - math.pi):.4f}  standard error {se:.4f}")

   100  estimate 3.0000  error 0.1416  standard error 0.1642
  1000  estimate 3.1280  error 0.0136  standard error 0.0519
 10000  estimate 3.1356  error 0.0060  standard error 0.0164
100000  estimate 3.1436  error 0.0020  standard error 0.0052

Each point is a Bernoulli trial with p = π/4, so the standard error of the estimate is 4 · sqrt(p(1 - p) / n), and the actual errors are of that size. (A single run can be luckier or unluckier than the standard error; it is a typical size, not a guarantee.) Getting π right to six decimal places this way would take about ten trillion points, which is why nobody computes π like this. Monte Carlo is valuable not because it is fast, but because it works when there is no formula.

Here is a question with a formula you might not know: how many rolls of a die does it take, on average, to see all six faces? Simulate first:

import random

rng = random.Random(6)

def rolls_to_see_all():
    seen = set()
    rolls = 0
    while len(seen) < 6:
        seen.add(rng.randint(1, 6))
        rolls += 1
    return rolls

trials = 20000
results = [rolls_to_see_all() for _ in range(trials)]
print("simulated mean:", sum(results) / trials)
print("longest run:   ", max(results))
print("exact:         ", round(sum(6 / k for k in range(1, 7)), 2))

simulated mean: 14.67535
longest run:    67
exact:          14.7

The exact answer comes from linearity of expectation. Once you have seen k faces, each roll shows a new one with probability (6 - k) / 6, so you wait on average 6 / (6 - k) rolls; adding those waits for k = 0, …, 5 gives 6/6 + 6/5 + 6/4 + 6/3 + 6/2 + 6/1 = 14.7. In symbols that is \sum_{k=0}^{5} \frac{6}{6 - k}, and the code writes the same sum with the index counted the other way: sum(6 / k for k in range(1, 7)). The simulation agrees, and it also tells you something the formula does not: the spread. Sometimes it takes more than 50 rolls.

The randomised heuristics of the Algorithms path make a similar trade. The difference is that a Monte Carlo estimate comes with a measurable error: a small, known imprecision in exchange for a method that is simple and always applicable.

Common mistakes: reporting more digits than the standard error supports. Reusing the same random numbers for two estimates you then compare as if they were independent (or, when you want a fair comparison, forgetting to reuse them). And long loops in the browser: a million iterations of pure Python takes seconds, so start with 10,000 and increase only if the standard error is too large.

Practice: Only the Digits You Can Trust, Caught Without an Umbrella, The Queue at the Coffee Cart

The bootstrap

The sampling distribution needs many samples from the population, and in real life you have exactly one. The bootstrap is a clever substitute: treat your sample as if it were the population, and draw many new samples from it, each of the same size, with replacement (so a value can be picked twice or not at all). The spread of a statistic across these resamples estimates its standard error.

Here are the fifteen bus lateness values from Level 1. How precise is their median?

import math
import random

def median(data):
    s = sorted(data)
    n = len(s)
    return s[n // 2] if n % 2 else (s[n // 2 - 1] + s[n // 2]) / 2

def sd(values):
    m = sum(values) / len(values)
    return math.sqrt(sum((v - m) ** 2 for v in values) / len(values))

late = [3, 0, 5, 2, 0, 4, 1, 0, 3, 27, 2, 0, 6, 3, 1]
rng = random.Random(10)
medians, means = [], []
for _ in range(10000):
    resample = [rng.choice(late) for _ in late]
    medians.append(median(resample))
    means.append(sum(resample) / len(resample))
medians.sort()
means.sort()

print("median", median(late), " bootstrap SE", round(sd(medians), 2),
      " middle 95%:", medians[250], "to", medians[9749])
print("mean  ", round(sum(late) / 15, 2), " bootstrap SE", round(sd(means), 2),
      " middle 95%:", round(means[250], 2), "to", round(means[9749], 2))

median 2  bootstrap SE 0.87  middle 95%: 0 to 4
mean   3.8  bootstrap SE 1.67  middle 95%: 1.4 to 7.53

The median is 2 minutes, and across the resamples it mostly stays between 0 and 4. The mean is 3.8, but its resamples range from about 1.4 to 7.5, and the range is lopsided: resamples that happen to include the 27 twice or three times drag the mean far up. The bootstrap has measured, in numbers, the Level 1 observation that the mean is fragile when there is an outlier.

For the mean there is a formula to check against: the standard error is s / sqrt(n), with s the sample standard deviation. For the lateness data that is about 1.73, close to the bootstrap's 1.67. For the median, and for most statistics you will invent yourself, there is no simple formula, and that is where the bootstrap earns its place. The "middle 95%" of the resampled values is a first confidence interval, the subject of Level 5.

Common mistakes: resampling without replacement (you just get the same sample back, shuffled). Resamples of a different size from the original. Expecting the bootstrap to rescue a tiny or unrepresentative sample: it can only reuse the values you have, and it never produces a value that was not in the data. And bootstrapping values that are not independent, such as consecutive days of a trend, where each resample should keep neighbouring days together.

Practice: Which Summary Is Steadiest?, Did the New Layout Help?

Practice

These problems are required on this path at this level; the path adds one problem from each of the level's topics, picked for you and kept.

• What Is a Scratch Card Worth?
• Accept the Batch?
• How Many Tickets for the Ferry?
• A Thousand Small Surveys
• Only the Digits You Can Trust
• Which Summary Is Steadiest?

Checklist

You can now:

• Describe a random variable by its distribution, and compute its expected value and variance exactly.
• Use linearity of expectation, and know when variances add.
• Compute binomial probabilities with math.comb, and simulate the trials to check them.
• Explain where the Poisson distribution comes from, compute its probabilities, and use "mean equals variance" as a check on count data.
• Use the normal distribution and the 68–95–99.7 rule, and check it with random.gauss and statistics.NormalDist.
• State the law of large numbers and explain why it does not support the gambler's fallacy.
• Simulate a sampling distribution, compute the standard error σ / sqrt(n), and see the central limit theorem at work.
• Use Monte Carlo simulation to answer questions without a formula, and quote the result with an honest standard error.
• Bootstrap the standard error of any statistic, and say when the bootstrap cannot help.