Model uncertainty with random variables and distributions, then estimate, test and update beliefs from data, with Python simulations you run in the browser.

Probability & Statistical Inference

## Why this matters

Two theorems explain why statistics works at all. The **law of large numbers (LLN)** says averages settle down to the true mean as you collect more data. The **central limit theorem (CLT)** says the error in that average is approximately normal, with a size you can calculate. Every confidence interval, t-test, A/B test and insurance reserve you will build leans on them. In this lesson you will check both by simulation rather than taking them on trust.

## The law of large numbers

Let `X1, X2, ..., Xn` be independent draws from the same distribution with mean `μ`. The sample mean

```text
X̄n = (X1 + X2 + ... + Xn) / n
```

gets closer and closer to `μ` as `n` grows. Flip a fair coin 10 times and you might see 70% heads; flip it 10,000 times and you will be very close to 50%.

The LLN is about the **average**, not about individual outcomes. It does not say that tails becomes "due" after a run of heads (the gambler's fallacy). Early streaks are not corrected; they are diluted by the many flips that follow.

## How fast? The standard error

If each draw has standard deviation `σ`, the sample mean has standard deviation

```text
SE(X̄n) = σ / sqrt(n)
```

This **standard error** shrinks with the square root of `n`. To halve the error you need four times as much data; to cut it by 10 you need 100 times as much. That single fact explains why survey sample sizes, lab replicates and simulation run counts look the way they do.

## The central limit theorem

The CLT adds the shape. For large enough `n`, whatever the original distribution (as long as it has a finite variance):

```text
X̄n is approximately Normal(μ, σ / sqrt(n))
```

The original data can be skewed, lumpy or discrete. The averages still come out bell-shaped. That is why the normal distribution appears in places where the raw data is not normal at all.

## Worked example: averaging skewed waiting times

Waiting times between calls are exponential with mean 1 minute (and SD 1 minute). A single waiting time is very skewed: most are short, a few are long. Now average 30 of them, and repeat that 5,000 times:

```python
import numpy as np

rng = np.random.default_rng(42)          # fixed seed: same numbers every run
draws = rng.exponential(scale=1.0, size=(5000, 30))
means = draws.mean(axis=1)               # one mean per row

print(means.mean())                      # about 1.00 (LLN)
print(means.std())                       # about 0.18 = 1 / sqrt(30) (standard error)
```

A histogram of `means` looks close to a normal curve centred on 1 with SD about 0.18, even though a histogram of single draws looks like a steep slide. Try `n = 2` and the histogram of means is still clearly skewed; by `n = 30` it is nearly symmetric. A quick numeric measure is **skewness** (0 for a symmetric shape): an exponential has skewness 2, and the mean of `n` draws has skewness about `2 / sqrt(n)`.

### Reading the code

- `np.random.default_rng(42)` creates a generator with a fixed seed. Always seed simulations you want others to reproduce.
- `size=(5000, 30)` makes a table: 5,000 experiments (rows) of 30 draws (columns).
- `mean(axis=1)` averages across each row.

## How large is "large enough"?

There is no single number. The rule of thumb `n ≥ 30` works for moderately skewed data. Symmetric data needs far less; heavily skewed data (insurance claims, income, file sizes) can need hundreds. With proportions, a common rule is at least 10 expected successes and 10 expected failures. When in doubt, simulate, exactly as above.

The CLT also fails when the variance is infinite. Some heavy-tailed distributions (such as the Cauchy) have averages that never settle down. These are rare in practice but do appear in some finance and physics data.

## Common mistakes

- **Applying the CLT to individual values.** The CLT is about averages and sums. A single claim size stays skewed no matter how many claims you collect.
- **Gambler's fallacy.** The LLN does not balance out short runs.
- **Confusing SD and SE.** The SD describes individual values; the SE describes the uncertainty in the mean. Report the one your reader needs.
- **Ignoring dependence.** Both theorems assume independent draws. Daily stock returns or time-series sensor readings are correlated, and the effective sample size is smaller than `n`.

## Recap

- LLN: the sample mean converges to the true mean.
- Standard error: `σ / sqrt(n)`; four times the data halves the error.
- CLT: sample means are approximately normal for large `n`, whatever the original shape.
- Simulate when you are unsure whether `n` is large enough.

## Try it

In the **Do** tab, simulate thousands of sample means with NumPy and check both theorems numerically.


See by simulation why sample means converge to the true mean and why their distribution becomes normal.

LLN and CLT

But what is the Central Limit Theorem?

Law of large numbers | Probability and Statistics | Khan Academy

## Why this matters

Many real probability questions have no tidy formula. What is the chance a project with 12 uncertain tasks finishes on time? How likely is a portfolio to lose more than 10% in a year? What share of manufactured assemblies will be out of tolerance when five parts each vary a little? **Monte Carlo simulation** answers all of them the same way: build a model of the random process, run it thousands of times, and count. It is a daily tool for data scientists, actuaries and operations research analysts, and it is also the best way to check your own formulas.

## The recipe

1. **Model** one run of the process with random draws.
2. **Repeat** it `N` times with a seeded random generator.
3. **Count** the runs where the event happens.
4. **Estimate** `p̂ = count / N`.
5. **Report the error**: `SE = sqrt(p̂ (1 - p̂) / N)`.

Step 5 is what separates a simulation from a guess. By the CLT (previous lesson), `p̂` is roughly normal around the true `p`, so `p̂ ± 2 SE` is an approximate 95% range.

## Worked example: the birthday problem

In a room of 23 people, what is the probability that at least two share a birthday (ignoring leap years and assuming all 365 days are equally likely)? The exact answer is about 0.507, which surprises most people. Here is a simulation using only the standard library:

```python
import random

def shared_birthday(n_people, rng):
    days = [rng.randrange(365) for _ in range(n_people)]
    return len(set(days)) < n_people       # a repeat shrinks the set

rng = random.Random(2026)                  # seeded: reproducible
trials = 20_000
hits = sum(shared_birthday(23, rng) for _ in range(trials))
p_hat = hits / trials
se = (p_hat * (1 - p_hat) / trials) ** 0.5
print(p_hat, se)                           # about 0.507, SE about 0.0035
```

With 20,000 trials the standard error is about 0.0035, so the estimate is reliable to roughly ±0.007. You can check it against the exact formula:

```text
P(no shared birthday) = (365/365) × (364/365) × ... × (343/365) ≈ 0.493
P(shared)             ≈ 1 - 0.493 = 0.507
```

## Worked example: estimating π

Throw random points into the unit square. The share that lands inside the quarter circle `x² + y² ≤ 1` is `π/4`:

```python
rng = random.Random(1)
N = 100_000
inside = sum(rng.random() ** 2 + rng.random() ** 2 <= 1 for _ in range(N))
print(4 * inside / N)       # about 3.14
```

This is a toy, but the same idea (estimate an area or integral by random sampling) is how physicists, engineers and Bayesian statisticians compute quantities that have no closed form.

## Worked example: tolerance stacking

An assembly is three parts end to end, each Normal(mean 10.0 mm, SD 0.05 mm), and it must be at most 30.12 mm long. With NumPy:

```python
import numpy as np
rng = np.random.default_rng(0)
total = rng.normal(10.0, 0.05, size=(100_000, 3)).sum(axis=1)
print((total > 30.12).mean())      # about 0.083
```

Here a formula exists (the sum is normal with SD `0.05 × sqrt(3) ≈ 0.087`), and the simulation agrees with it. But if one part were uniform, or the parts were correlated, the simulation would need only a one-line change while the formula would need new math.

## Choosing N

The error falls as `1/sqrt(N)`. For a probability near 0.5:

| N | SE |
|---|---|
| 1,000 | 0.016 |
| 10,000 | 0.005 |
| 100,000 | 0.0016 |

Rare events need more runs. To pin down a 1-in-1,000 probability to within about ±10% of its value (a 95% range), you need around 400,000 runs. Analysts use variance-reduction tricks (importance sampling, stratified sampling) for very rare events.

## Common mistakes

- **No seed.** Results change on every run and nobody can reproduce your figure.
- **No error bar.** "The probability is 0.507" from 200 trials (SE ≈ 0.035) is overconfident.
- **Modelling the wrong process.** A simulation is only as good as its assumptions. Birthdays are not quite uniform across the year; say so.
- **Reusing the same random numbers by accident**, such as reseeding inside the loop, which makes every trial identical.

## Recap

- Monte Carlo: model one run, repeat `N` times, count, divide.
- Always seed the generator and report `SE = sqrt(p̂(1 - p̂)/N)`.
- Error shrinks with `1/sqrt(N)`; rare events need many runs.
- Use simulation to answer questions with no formula and to check the ones that have a formula.

## Try it

In the **Do** tab, simulate the birthday problem and compare your estimate with the exact answer.


Estimate a probability by simulating a random process many times and report how accurate the estimate is.

Monte Carlo simulation

What is Monte Carlo Simulation?

6. Monte Carlo Simulation

## Why this matters

In Statistics Fundamentals you built a confidence interval for a mean with a formula. Formulas exist for a few statistics (means, proportions, regression slopes) under some assumptions. But what is the margin of error on a **median** repair time, a **90th percentile** response time, a ratio of two means, or a model's accuracy? The **bootstrap** gives a confidence interval for almost any statistic, using only your data and a computer. It is one of the most practical ideas in modern statistics.

## The core idea

If you could, you would draw many fresh samples from the population, compute the statistic on each, and look at how much it varies. You cannot, because you have one sample. The bootstrap's trick is to treat **your sample as a stand-in for the population** and draw new samples from it:

1. Take your data of size `n`.
2. Draw a **resample** of size `n` from it **with replacement** (some values appear twice, some not at all).
3. Compute the statistic on the resample.
4. Repeat `B` times (typically 2,000 to 10,000).
5. The spread of the `B` results estimates the sampling variability. For a 95% **percentile interval**, take the 2.5th and 97.5th percentiles of them.

"With replacement" is essential. Without it, every resample would be the original data shuffled, and the statistic would never change.

## Worked example: median repair time

A maintenance team logged 12 repair times (hours). The data is skewed by two long jobs, so the median is the better summary.

```python
import numpy as np

times = np.array([1.2, 0.8, 2.5, 1.1, 0.9, 6.0, 1.4, 1.7, 0.7, 2.1, 9.5, 1.3])
print(np.median(times))                    # 1.35

rng = np.random.default_rng(7)
boot = []
for _ in range(5000):
    resample = rng.choice(times, size=len(times), replace=True)
    boot.append(np.median(resample))

lo, hi = np.percentile(boot, [2.5, 97.5])
print(lo, hi)                              # roughly 1.0 to 2.3
```

Report it as: "median repair time 1.35 hours (95% bootstrap CI about 1.0 to 2.3 hours)". The interval is not symmetric around 1.35, and that is a feature: it reflects the skew in the data, which a plus-or-minus formula would hide.

## Bootstrapping other statistics

Swap `np.median` for any function of the data:

- `np.mean`, to check the textbook interval.
- `lambda x: np.percentile(x, 90)` for a 90th-percentile latency.
- A difference between two groups: resample **each group separately**, then subtract.
- A model's accuracy on a test set: resample the test rows.

`scipy.stats.bootstrap(...)` implements this (with a better "BCa" interval by default), and you will see it in data-science code. Writing the loop yourself once makes it clear what the library does.

## When the bootstrap works, and when it does not

The bootstrap assumes your sample represents the population and that observations are **independent**. It works well for means, medians, percentiles in the middle of the data, and most smooth statistics, given a reasonable `n` (say 20 or more).

It struggles with:

- **Very small samples.** Twelve values give a rough interval; five values give a misleading one.
- **Extremes.** The maximum of a resample can never exceed the sample maximum, so bootstrapping the max or the 99.9th percentile underestimates uncertainty.
- **Dependent data.** Time series and clustered data (many readings per patient, per machine) need block or cluster resampling, or the interval will be too narrow.

## How many resamples?

`B` controls only the **Monte Carlo** error of the interval endpoints, not the underlying uncertainty from your sample size. 1,000 resamples are fine for a quick look; use 5,000 to 10,000 for anything you publish. Always seed the generator so the interval is reproducible.

## Common mistakes

- **Resampling without replacement.** Every resample becomes identical to the data.
- **Resampling a different size.** Each resample must have the same `n` as the original.
- **Thinking more resamples means more certainty.** Only more real data narrows the interval.
- **Ignoring structure.** Resample whole clusters (patients, batches) when observations come in groups.

## Recap

- Resample your data with replacement, recompute the statistic, repeat thousands of times.
- The 2.5th and 97.5th percentiles of the results form a 95% percentile interval.
- Works for medians, percentiles, ratios and model metrics, not just means.
- Needs independent data and a reasonable sample; struggles with extremes.

## Try it

In the **Do** tab, write a reusable bootstrap function and use it for a median and a mean.


Build a percentile bootstrap confidence interval for a median or mean by resampling your data.

The bootstrap

Bootstrapping Main Ideas!!!

Bootstrap confidence intervals - explained

Why averages settle down (law of large numbers) and become normal (central limit theorem), and how to answer probability questions by simulation and the bootstrap.