Level 2 · Statistics & Probability: Spread and Chance

Level 1 found the centre of a dataset. This level asks the next two questions. First: how spread out are the values around that centre? Two classes can both average 70 on a test, one with every student between 65 and 75 and the other with half the class failing; the mean alone cannot tell them apart. You will measure spread with the variance and the standard deviation, locate values with percentiles and quartiles, and describe any single value by how unusual it is (its z-score).

Second: what is chance? Spread is really about uncertainty, and the language of uncertainty is probability. You will count equally likely outcomes to get exact probabilities, use complements to turn "at least one" questions into easy ones, and then do what this site does best: simulate. A few lines with a seeded random generator estimate a probability by running an experiment thousands of times, and you can check the estimate against the exact answer.

Two links to the other paths. The z-score is exactly the standardisation used in the Machine Learning path to put features on a common scale before measuring distances. And the quartiles come from a sorted list, the same sort that gave the median; in the Algorithms path you saw why sorting once and indexing is cheaper than searching repeatedly.

Variance: the average squared distance

To measure spread, start from the deviations x - mean of Level 1. Their sum is always zero, so their average is useless. Two fixes: take absolute values (the mean absolute deviation) or square them. Squaring is the standard choice, and the average squared deviation is the variance:

\frac{1}{n}\sum_{i=1}^{n} (x_i - \bar{x})^2

This is the "From formula to Python" pattern of Level 1: the summand (x_i - \bar{x})^2 becomes (x - m) ** 2, the sum over i becomes for x in data, and the \frac{1}{n} becomes the division by len(data).

def mean(data):
    return sum(data) / len(data)

def variance(data):
    m = mean(data)
    return sum((x - m) ** 2 for x in data) / len(data)

calm = [68, 70, 72, 69, 71, 70]
wild = [40, 100, 55, 85, 70, 70]
for name, d in (("calm", calm), ("wild", wild)):
    print(name, "mean", mean(d), "variance", round(variance(d), 2))

calm mean 70.0 variance 1.67
wild mean 70.0 variance 375.0

Both classes average 70, but the variance of wild is several hundred times larger. Why squares rather than absolute values? Squares are smooth (no corner at zero), which makes the algebra and the calculus of later levels work; the least-squares line in Level 3 and the squared-error loss in the Machine Learning path both come from this choice. The price is that squaring exaggerates large deviations: a value 30 away contributes 900, nine times as much as a value 10 away.

There is a second formula that needs only one pass: the variance is the mean of the squares minus the square of the mean,

\frac{1}{n}\sum_{i=1}^{n} x_i^2 - \bar{x}^2

Keep a running sum s (for \sum x_i) and a running sum of squares s2 (for \sum x_i^2):

wild = [40, 100, 55, 85, 70, 70]
n = s = s2 = 0
for x in wild:
    n += 1
    s += x
    s2 += x * x
print(s2 / n - (s / n) ** 2)

375.0

This is the one-pass pattern from the Algorithms path, but be careful with it. When the values are huge and close together (say times in milliseconds since 1970), s2 / n and (s / n) ** 2 are two enormous, nearly equal numbers, and subtracting them loses most of the digits. The two-pass version (compute the mean first, then the squared deviations) is numerically safer, and streaming code uses a careful update rule (Welford's method) instead.

Common mistakes: forgetting to square, so the deviations cancel to zero. Squaring the sum, (\sum x_i)^2 or sum(xs) ** 2, instead of summing the squares, \sum x_i^2 or sum(x * x for x in xs). Reading the variance in the data's units: the variance of scores is in "points squared", which is why the next section takes its square root.

Practice: Price Bands for the Map

Standard deviation, and why n − 1

The standard deviation is the square root of the variance. It is back in the data's own units, and it reads as "a typical distance from the mean".

import math

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

calm = [68, 70, 72, 69, 71, 70]
wild = [40, 100, 55, 85, 70, 70]
print(round(sd(calm), 2), round(sd(wild), 2))

1.29 19.36

So a typical calm score is about 1.3 points from 70 and a typical wild score about 19 points away. For data with a single hump and no strong skew, a useful rule of thumb is that about two thirds of the values lie within one standard deviation of the mean and about 95% within two. For skewed data (like the bus lateness of Level 1) the rule fails, and the quartiles of the next section describe the data better.

Now the puzzle every statistics course has to explain. Calculators and libraries offer two versions: divide by n or divide by n - 1. Which is right?

If your data is the whole population you care about (every student in this class), divide by n. That is the population variance.
If your data is a sample from a larger population, and you want to estimate the population's variance, divide by n - 1. That is the sample variance.

\sigma^2 = \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2

s^2 = \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2

The usual names are \sigma^2 (sigma squared) for the population variance and s^2 for the sample variance; their square roots \sigma and s are the two standard deviations. In code they share the sum ss = sum((x - m) ** 2 for x in xs) and differ only in the last step, ss / n against ss / (n - 1).

The reason is that the deviations are measured from the sample's own mean, and the sample mean is, by construction, the point that makes the sum of squared deviations as small as possible. The true mean is somewhere else, so deviations from it would be larger. Dividing by n therefore underestimates the spread on average, and dividing by the slightly smaller n - 1 corrects that exactly. You do not have to take this on trust: simulate it. Rolling a fair die has a known variance, 35/12 ≈ 2.917. Take thousands of small samples of three rolls and average both versions:

import random
rng = random.Random(1)
faces = [1, 2, 3, 4, 5, 6]
m = sum(faces) / 6
true_var = sum((f - m) ** 2 for f in faces) / 6
print("true variance of one roll:", round(true_var, 4))

n, trials = 3, 20000
by_n = by_n1 = 0
for _ in range(trials):
    sample = [rng.choice(faces) for _ in range(n)]
    sm = sum(sample) / n
    ss = sum((x - sm) ** 2 for x in sample)
    by_n += ss / n
    by_n1 += ss / (n - 1)
print("average of ss/n:    ", round(by_n / trials, 3))
print("average of ss/(n-1):", round(by_n1 / trials, 3))

true variance of one roll: 2.9167
average of ss/n:     1.93
average of ss/(n-1): 2.894

Dividing by n is systematically low: on average it gives about two thirds of the truth, which is exactly (n - 1) / n for n = 3. Dividing by n - 1 lands close to 2.917, and the small remaining gap is simulation noise that shrinks with more trials. An estimator whose average equals the truth is called unbiased; Level 6 returns to this idea. For large samples the two versions barely differ (dividing by 999 or 1000 hardly matters); for small samples the difference is real.

Python's statistics module has both, and they agree with the formulas:

import statistics
wild = [40, 100, 55, 85, 70, 70]
print(statistics.pvariance(wild), statistics.variance(wild))
print(round(statistics.pstdev(wild), 3), round(statistics.stdev(wild), 3))

375 450
19.365 21.213

The p stands for population (divide by n); the plain names are the sample versions (divide by n - 1).

Common mistakes: using the population formula on a small sample and underestimating the spread. Taking the square root of each deviation instead of the variance. Expecting the rule of thumb (two thirds within one standard deviation) to hold for skewed data. And calling statistics.variance on a single value: with n - 1 = 0 there is no spread to estimate, and it raises an error.

Practice: Whole Batch or Sample?, One Spread from Many Class Reports, One Number for Every Day

Percentiles and quartiles

The median splits sorted data in half. Percentiles generalise it: the 90th percentile is the value below which 90% of the data lies. The quartiles are the 25th, 50th and 75th percentiles, called Q1, Q2 (the median) and Q3.

There is no single agreed formula for a percentile of a finite list, because the exact position usually falls between two values. A common and simple choice: the p-th percentile sits at position p / 100 * (n - 1) in the sorted list, and if that position is not a whole number, interpolate linearly between its two neighbours.

def percentile(data, p):
    s = sorted(data)
    pos = p / 100 * (len(s) - 1)
    lo = int(pos)
    frac = pos - lo
    if lo + 1 < len(s):
        return s[lo] + frac * (s[lo + 1] - s[lo])
    return s[lo]

late = [3, 0, 5, 2, 0, 4, 1, 0, 3, 27, 2, 0, 6, 3, 1]
for p in (0, 25, 50, 75, 90, 100):
    print(p, percentile(late, p))

The 0th and 100th percentiles are the minimum and the maximum, and the 50th is the median, as they should be. Python's statistics.quantiles with method="inclusive" uses the same rule; its default method uses a slightly different position formula, and spreadsheets have several more. On large datasets the methods agree closely; on fifteen values they differ noticeably (below, Q3 is 3.5 by one method and 4.0 by the other), so say which one you used.

import statistics
late = [3, 0, 5, 2, 0, 4, 1, 0, 3, 27, 2, 0, 6, 3, 1]
print(statistics.quantiles(late, n=4, method="inclusive"))
print(statistics.quantiles(late, n=4))

[0.5, 2.0, 3.5]
[0.0, 2.0, 4.0]

The interquartile range (IQR) is Q3 − Q1: the width of the middle half of the data. Unlike the range and the standard deviation, it ignores the top and bottom quarters entirely, so a single outlier cannot move it. That makes it the natural measure of spread to pair with the median. It also gives a widely used rule for flagging outliers: anything more than 1.5 IQRs below Q1 or above Q3.

def percentile(data, p):
    s = sorted(data)
    pos = p / 100 * (len(s) - 1)
    lo = int(pos)
    return s[lo] + (pos - lo) * (s[min(lo + 1, len(s) - 1)] - s[lo])

late = [3, 0, 5, 2, 0, 4, 1, 0, 3, 27, 2, 0, 6, 3, 1]
q1, q3 = percentile(late, 25), percentile(late, 75)
iqr = q3 - q1
low_fence, high_fence = q1 - 1.5 * iqr, q3 + 1.5 * iqr
print("IQR", iqr, "fences", low_fence, high_fence)
print("flagged:", [x for x in late if x < low_fence or x > high_fence])

IQR 3.0 fences -4.0 8.0
flagged: [27]

The 27-minute day is flagged, and nothing else. The rule does not say the value is wrong, only that it deserves a look. The five numbers minimum, Q1, median, Q3 and maximum (the five-number summary) are what a box plot draws.

Common mistakes: computing percentiles on unsorted data. Assuming every tool uses the same percentile method. Treating a flagged outlier as an error to delete; it may be the most important observation (the day the bus broke down). And forgetting that "in the 90th percentile" means better than 90% of the group, not a score of 90%.

Practice: Percentiles of the Fun Run, The Numbers Behind a Box Plot, Percentile Ranks for Every Candidate

z-scores: how unusual is a value?

A score of 78 means little on its own. It is excellent on a test where the class averaged 60 with a standard deviation of 8, and mediocre where the class averaged 85. The z-score turns a value into "how many standard deviations above (or below) the mean":

z = \frac{x - \bar{x}}{s}

where \bar{x} is the group's mean and s its standard deviation.

Aisha scored 78 in maths and 82 in history. Which is the more impressive result?

import math

def mean_sd(data):
    m = sum(data) / len(data)
    return m, math.sqrt(sum((x - m) ** 2 for x in data) / len(data))

maths   = [55, 62, 78, 58, 70, 49, 66, 60, 57, 65]
history = [80, 85, 82, 90, 78, 88, 84, 79, 86, 88]
for name, data, score in (("maths", maths, 78), ("history", history, 82)):
    m, s = mean_sd(data)
    print(f"{name}: mean {m:.1f}, sd {s:.2f}, z = {(score - m) / s:.2f}")

maths: mean 62.0, sd 7.80, z = 2.05
history: mean 84.0, sd 3.92, z = -0.51

Her 78 in maths is just over two standard deviations above the class mean, the best in the class. Her 82 in history, the higher raw mark, is half a standard deviation below that class's mean. The z-score removes the units and the difficulty of each test, so the two numbers become comparable.

Converting a whole dataset to z-scores always gives a new dataset with mean 0 and standard deviation 1: subtracting the mean moves the centre to 0, and dividing by the standard deviation rescales the spread to 1. The shape does not change: a skewed dataset stays skewed.

That last point is the main way z-scores mislead. "More than 2 standard deviations from the mean is unusual" is based on the bell-shaped (normal) distribution of Level 4. For skewed data such as waiting times, values with z above 2 are far more common than 5%, and a z-score near 0 may not be a typical value at all.

Common mistakes: mixing population and sample standard deviations between groups you want to compare (use the same one throughout). Computing a z-score against the wrong group's mean. Reading a z-score as a percentage or a probability; it is a distance, measured in standard deviations.

Practice: Her Strongest Subject, Peeling Off the Glitches

Normalisation

Normalising data means putting it on a standard scale. The two common recipes are standardisation (convert to z-scores: mean 0, standard deviation 1) and min-max scaling (map the smallest value to 0 and the largest to 1).

import math

def standardise(data):
    m = sum(data) / len(data)
    s = math.sqrt(sum((x - m) ** 2 for x in data) / len(data))
    return [(x - m) / s for x in data]

def min_max(data):
    lo, hi = min(data), max(data)
    return [(x - lo) / (hi - lo) for x in data]

heights = [150, 160, 165, 170, 185]
print([round(z, 2) for z in standardise(heights)])
print([round(v, 2) for v in min_max(heights)])

[-1.38, -0.52, -0.09, 0.35, 1.64]
[0.0, 0.29, 0.43, 0.57, 1.0]

Why bother? Because many methods compare values on different scales. Suppose you describe people by height in centimetres and yearly income in pounds and ask who is most similar to whom, by distance. The income numbers are thousands of times larger, so the distance is decided almost entirely by income, and height is ignored. This is exactly the problem that feature scaling solves in the Machine Learning path, before nearest-neighbour methods measure distances:

import math

people = {"A": (160, 30000), "B": (190, 31000), "C": (161, 36000), "D": (175, 60000)}

def dist(p, q):
    return math.dist(people[p], people[q])

print("raw:    A-B", round(dist("A", "B")), " A-C", round(dist("A", "C")))

cols = list(zip(*people.values()))
lo = [min(c) for c in cols]
hi = [max(c) for c in cols]
scaled = {k: tuple((v - l) / (h - l) for v, l, h in zip(vals, lo, hi)) for k, vals in people.items()}
print("scaled: A-B", round(math.dist(scaled["A"], scaled["B"]), 2),
      " A-C", round(math.dist(scaled["A"], scaled["C"]), 2))

raw:    A-B 1000  A-C 6000
scaled: A-B 1.0  A-C 0.2

On the raw numbers, A looks much closer to B: the 30 cm height difference barely registers, and the distance is essentially the £1000 income gap. After scaling each column to 0–1 (heights by their range of 30 cm, incomes by their range of £30000), the 30 cm difference counts fully and A is much closer to C. Neither answer is "true": the choice of scale is a decision about what similarity should mean, and it should be made on purpose.

Which recipe? Min-max scaling keeps everything in [0, 1] but depends on the two most extreme values, so one outlier squashes all the others together. Standardisation is less sensitive to a single extreme value and is the usual default. With strong outliers, some people scale by the median and the IQR instead.

Common mistakes: dividing by zero when every value is the same (the range and standard deviation are 0; there is nothing to scale). Computing the mean and standard deviation on all the data, including the test data, before a machine-learning experiment; the ML path explains why the scaling must be learnt from the training data only. And forgetting to apply the same transformation to new values later: a new person must be scaled with the old minimum and maximum.

Practice: Put the Flat Features on One Scale, Who Is Most Alike Depends on the Scale

Probability by counting equally likely outcomes

When every outcome of an experiment is equally likely, the probability of an event A is a count:

P(A) = \frac{\text{number of outcomes in } A}{\text{total number of outcomes}}

Rolling two dice has 36 equally likely outcomes, the pairs (first, second). itertools.product lists them all, and fractions.Fraction keeps the answer exact:

from itertools import product
from fractions import Fraction

outcomes = list(product(range(1, 7), repeat=2))
sevens = [o for o in outcomes if sum(o) == 7]
print(len(outcomes), sevens)
print(Fraction(len(sevens), len(outcomes)))

36 [(1, 6), (2, 5), (3, 4), (4, 3), (5, 2), (6, 1)]
1/6

Counting the outcomes for every total gives the whole probability distribution of the sum:

from itertools import product

counts = {}
for a, b in product(range(1, 7), repeat=2):
    counts[a + b] = counts.get(a + b, 0) + 1
for total in range(2, 13):
    print(f"{total:2d} {counts[total]}/36 {'#' * counts[total]}")

 2 1/36 #
 3 2/36 ##
 4 3/36 ###
 5 4/36 ####
 6 5/36 #####
 7 6/36 ######
 8 5/36 #####
 9 4/36 ####
10 3/36 ###
11 2/36 ##
12 1/36 #

Press Run to see the same distribution as a bar chart:

from itertools import product

totals = range(2, 13)
ways = [sum(a + b == t for a, b in product(range(1, 7), repeat=2)) for t in totals]
plot(list(totals), [w / 36 for w in ways], kind="bar", xlabel="total of two dice", ylabel="probability")

The histogram is a triangle: 7 is six times as likely as 2, because six pairs make 7 and only (1, 1) makes 2. This is the crucial condition in the counting rule. The eleven totals 2 to 12 are not equally likely, so "P(total is 12) = 1/11" is wrong; the equally likely outcomes are the 36 pairs. Whenever you count, first make sure the things you are counting really are equally likely.

For two events that cannot happen together (disjoint events), probabilities add: P(total is 2 or 12) = 1/36 + 1/36 = 1/18. For events that can overlap, adding counts the overlap twice; counting the outcomes directly avoids the mistake.

Dice are rolled one after the other, so order matters. When you draw a group instead (a team, a hand of cards, tickets from a box), the equally likely outcomes are the sets of items, and order does not matter. The number of ways to choose k items out of n is the binomial coefficient \binom{n}{k} = \frac{n!}{k!\,(n-k)!}, math.comb(n, k) in Python; itertools.combinations lists the sets themselves:

from itertools import combinations
import math

people = ["Ana", "Ben", "Cy", "Dee", "Eli"]
teams = list(combinations(people, 2))
print(len(teams), math.comb(5, 2))
print(teams[:4])

10 10
[('Ana', 'Ben'), ('Ana', 'Cy'), ('Ana', 'Dee'), ('Ana', 'Eli')]

To count the sets with a given make-up, choose each part separately and multiply. A box holds 4 red and 6 blue tickets and you draw 3 without putting them back. The draws with exactly 2 red choose 2 of the 4 red and 1 of the 6 blue:

from fractions import Fraction
import math

good = math.comb(4, 2) * math.comb(6, 1)
print(good, math.comb(10, 3), Fraction(good, math.comb(10, 3)))

36 120 3/10

Common mistakes: counting outcomes that are not equally likely (totals instead of pairs; "it either happens or it does not, so it is 50%"). Treating (1, 2) and (2, 1) as one outcome: the dice are different objects even if they look the same. Mixing ordered and unordered counting: if you count the drawn group as a set (comb), count every outcome as a set too. Using floats and getting 0.16666666666666666 when the exact answer 1/6 is easier to check.

Practice: Two Unusual Dice, A Good Opening Hand

Complements: the easy way to "at least one"

The complement of an event is everything else: "not A". Exactly one of them happens, so

P(\text{not } A) = 1 - P(A)

This turns many hard questions into easy ones, especially questions with "at least one". What is the probability of at least one six in four rolls of a die? "At least one" includes one six, two sixes, three and four, in many arrangements. Its complement is a single case: no six at all, which is (5/6) for each roll, and the rolls do not affect each other, so the chances multiply:

from fractions import Fraction
from itertools import product

no_six = Fraction(5, 6) ** 4
print("P(no six)       =", no_six)
print("P(at least one) =", 1 - no_six, "≈", round(float(1 - no_six), 4))

rolls = list(product(range(1, 7), repeat=4))
print("by counting:     ", Fraction(sum(6 in r for r in rolls), len(rolls)))

P(no six)       = 625/1296
P(at least one) = 671/1296 ≈ 0.5177
by counting:      671/1296

Counting all 1296 outcomes agrees with the complement. So a gambler betting on "at least one six in four rolls" wins slightly more often than not: a bet that was a real money-maker in the seventeenth century.

The same trick solves the famous birthday problem: in a group of k people, how likely is it that at least two share a birthday (ignoring 29 February, and assuming all 365 days equally likely)? The complement is "all birthdays different": the second person avoids 1 day, the third avoids 2, and so on.

def p_shared(k):
    p_all_different = 1.0
    for i in range(k):
        p_all_different *= (365 - i) / 365
    return 1 - p_all_different

for k in (10, 23, 40, 60):
    print(k, round(p_shared(k), 3))

The loop computes a product, \prod_{i=0}^{k-1} \frac{365 - i}{365}, and products follow the same pattern as sums: math.prod((365 - i) / 365 for i in range(k)) is the one-line version. Press Run to see how fast the chance climbs:

import math

ks = range(1, 61)
plot(list(ks), [1 - math.prod((365 - i) / 365 for i in range(k)) for k in ks],
     xlabel="people in the room", ylabel="P(shared birthday)")

With just 23 people the chance of a shared birthday is already above one half. Intuition fails here because it thinks about my birthday matching someone's, while the question is about any pair, and 23 people form 253 pairs.

Common mistakes: adding instead of multiplying for "and" across independent rolls (1/6 + 1/6 + 1/6 + 1/6 = 2/3 is wrong: with seven rolls that reasoning would give more than 1). Taking the complement of "at least one" to be "exactly one" or "all of them"; it is "none". Multiplying probabilities of events that influence each other; Level 3 handles that case with conditional probability.

Practice: How Many Seeds to Plant?, First to Roll a Winner

Simulation with `random`

Counting works when the outcomes are few and equally likely. For everything else there is simulation: run the random experiment many times in code and count how often the event happens. The fraction of runs in which it happened estimates its probability.

Python's random module provides the randomness. It is a pseudo-random generator: a deterministic algorithm whose output looks random. Seeding it fixes its starting point, so a simulation gives the same numbers every time it runs, which is what you want for checking, debugging and sharing results. A private generator, random.Random(seed), is better than seeding the global one, because other code cannot disturb its sequence.

import random
rng = random.Random(42)
print(rng.random())            # a float in [0, 1)
print(rng.randint(1, 6))       # a whole number from 1 to 6, both included
print(rng.choice(["H", "T"]))  # one element of a list
rng2 = random.Random(42)
print(rng2.random())           # same seed, same sequence

0.6394267984578837
1
T
0.6394267984578837

Now estimate the probability of at least one six in four rolls, which you know exactly is 671/1296 ≈ 0.5177, with more and more trials:

import random

exact = 1 - (5 / 6) ** 4
for n in (100, 1000, 10000, 100000):
    rng = random.Random(7)
    hits = 0
    for _ in range(n):
        if any(rng.randint(1, 6) == 6 for _ in range(4)):
            hits += 1
    est = hits / n
    print(f"{n:>6} trials: estimate {est:.4f}, error {abs(est - exact):.4f}")

   100 trials: estimate 0.4700, error 0.0477
  1000 trials: estimate 0.5110, error 0.0067
 10000 trials: estimate 0.5169, error 0.0008
100000 trials: estimate 0.5173, error 0.0005

With 100 trials the estimate is off by almost 0.05; with 100,000 it is right to about three decimal places. The error of a simulation shrinks roughly like 1/\sqrt{n}: a hundred times more trials buys about ten times more accuracy. Level 4 shows why (it is the law of large numbers, and the standard deviation of an average), and Level 5 turns it into a confidence interval. For now, the working rule is: an estimate from n trials is only good to roughly 1 / sqrt(n), so 10,000 trials give about two reliable decimal places.

Simulation shines where intuition fails. In the Monty Hall game there are three doors, a car behind one. You pick a door; the host, who knows where the car is, opens another door with no car behind it, and offers you the chance to switch to the remaining closed door. Most people feel it makes no difference. Play it ten thousand times:

import random

rng = random.Random(3)
stay = switch = 0
trials = 10000
for _ in range(trials):
    car = rng.randrange(3)
    pick = rng.randrange(3)
    opened = next(d for d in range(3) if d != pick and d != car)
    other = 3 - pick - opened      # the doors are 0, 1, 2, which add up to 3
    stay += pick == car
    switch += other == car
print("stay wins:  ", stay / trials)
print("switch wins:", switch / trials)

stay wins:   0.3344
switch wins: 0.6656

Switching wins about two thirds of the time. The simulation does not explain why, but it tells you what to explain: staying wins only when your first pick was right (probability 1/3), and in every other case the host's forced choice leaves the car behind the other door. Estimate first, then find the exact reasoning and check that the two agree. That habit, simulate and then compute, runs through every level of this path.

Common mistakes: forgetting to seed, so results change on every run and cannot be checked. Re-seeding inside the loop, so every trial is identical. Using randint(0, 6) for a die (it includes both ends, so that gives seven faces). Trusting three decimal places from a thousand trials. And simulating a different game from the one asked: if the host opened a door at random (and sometimes revealed the car), the answer would change, so the code must follow the rules exactly.

Practice: Simulating a Lucky Streak, Nobody Draws Their Own Name, The Game Show with Many Doors

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.

Checklist

You can now:

Compute the variance and standard deviation by hand, in two passes or one, and say what the standard deviation means in the data's units.
Explain when to divide by n and when by n - 1, and demonstrate the bias of dividing by n with a simulation.
Compute percentiles by sorting and interpolating, find the quartiles and the IQR, and flag outliers with the 1.5 × IQR rule.
Convert values to z-scores to compare results from different groups, and say why z-scores mislead for skewed data.
Normalise data by standardisation or min-max scaling, and explain why unscaled features distort distances.
Compute exact probabilities by counting equally likely outcomes with itertools.product and Fraction.
Use complements to answer "at least one" questions, including the birthday problem.
Estimate a probability by a seeded simulation, compare it with the exact answer, and judge its accuracy from the number of trials.

Spread and Chance