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

Statistics Fundamentals taught you to summarise data you already have. Probability goes the other way: it describes data you have not seen yet. How many defective parts will tomorrow's batch contain? How much will an insurance policy pay out next year? How long until the next server request arrives? Each of these is a **random variable**, and two numbers, its expected value and its variance, drive pricing, staffing, quality limits and model evaluation in every career this course serves.

## The core idea

A **random variable** is a rule that turns the outcome of a random process into a number. Roll a die and record the face: that number is a random variable, usually written with a capital letter, `X`. A particular value it takes is written in lower case, `x`.

A **discrete** random variable takes separate, countable values (0, 1, 2, ... defects). A **continuous** one can take any value in a range (a waiting time of 2.37 seconds). This lesson focuses on discrete ones; continuous variables come in lesson 3.

A discrete random variable is fully described by its **probability mass function (PMF)**: a list of each possible value and its probability. Two rules always hold: every probability is between 0 and 1, and the probabilities add up to exactly 1.

```text
X = face of a fair die
x      1    2    3    4    5    6
P(x)  1/6  1/6  1/6  1/6  1/6  1/6     (sum = 1)
```

## Expected value

The **expected value** (or mean) `E[X]` is the probability-weighted average of the values:

```text
E[X] = sum of x * P(x)
die: (1 + 2 + 3 + 4 + 5 + 6) / 6 = 3.5
```

You can never roll 3.5. The expected value is the long-run average over many repetitions, not a value you expect to see on one try.

## Variance and standard deviation

The **variance** measures how far values typically land from the mean, weighting each squared distance by its probability:

```text
Var(X) = sum of P(x) * (x - E[X])^2
die: [(1-3.5)^2 + (2-3.5)^2 + ... + (6-3.5)^2] / 6 = 17.5 / 6 = 35/12 ≈ 2.917
SD(X)  = sqrt(Var(X)) ≈ 1.708
```

A handy equivalent form is `Var(X) = E[X^2] - (E[X])^2`. It is quicker by hand but can lose precision on a computer when the two terms are large and nearly equal, so code usually uses the first form.

## Worked example: an insurance policy

A one-year policy pays $10,000 if a claim happens (probability 0.01) and $0 otherwise.

```text
x        0       10,000
P(x)   0.99       0.01

E[X]   = 0 * 0.99 + 10,000 * 0.01 = $100
Var(X) = 0.99 * (0 - 100)^2 + 0.01 * (10,000 - 100)^2
       = 9,900 + 980,100 = 990,000
SD(X)  ≈ $995
```

The insurer's expected cost is $100 per policy, but the standard deviation is almost ten times larger. A single policy is very risky; a portfolio of thousands is not, because (as you will see in module 2) averages are far less variable than single values. This is the basic logic of actuarial pricing, and it is the opening topic of the Society of Actuaries' Exam P (Probability). Exam details change, so read the current syllabus at [soa.org](https://www.soa.org/education/exam-req/edu-exam-p-detail/). Checked on 2026-10-23.

## Rules you will use constantly

For constants `a` and `b` and random variables `X` and `Y`:

- `E[aX + b] = a E[X] + b`: shifting and scaling move the mean the same way.
- `Var(aX + b) = a^2 Var(X)`: adding a constant does not change spread; scaling by `a` multiplies variance by `a^2`.
- `E[X + Y] = E[X] + E[Y]`, always, even if `X` and `Y` are related.
- `Var(X + Y) = Var(X) + Var(Y)` only when `X` and `Y` are **independent**.

Example: the sum of two independent dice has mean 3.5 + 3.5 = 7 and variance 35/12 + 35/12 ≈ 5.83.

## In Python

A dictionary from value to probability is a natural PMF:

```python
die = {face: 1 / 6 for face in range(1, 7)}
mean = sum(x * p for x, p in die.items())
print(mean)   # 3.5
```

In data work you will also meet `numpy` arrays of values and probabilities, and `scipy.stats` objects such as `scipy.stats.binom(n, p).mean()` that compute these for named distributions.

## Common mistakes

- **Forgetting to check the PMF.** If the probabilities do not add to 1 (within rounding), every result built on them is wrong. Validate first.
- **Adding standard deviations.** Variances of independent variables add; standard deviations do not. Two independent dice have SD `sqrt(5.83) ≈ 2.42`, not `1.71 + 1.71`.
- **Reading E[X] as "most likely".** The most likely value is the mode. For the insurance policy the mode is $0, while the mean is $100.
- **Confusing population and sample variance.** A PMF describes the whole population, so there is no `n - 1` here.

## Recap

- A random variable maps random outcomes to numbers; a PMF lists values and probabilities that sum to 1.
- `E[X]` is the long-run average; `Var(X)` is the probability-weighted squared distance from it.
- Means always add; variances add only for independent variables.

## Try it

Open the **Do** tab and write the variance and validity check for a PMF stored as a Python dictionary.


Describe an uncertain quantity with a random variable and compute its expected value and variance.

Random variables

Expected Value and Variance of Discrete Random Variables

Mean (expected value) of a discrete random variable | AP Statistics | Khan Academy

## Why this matters

Writing a PMF by hand works for a die. Real problems have thousands of possible values: the number of defective chips in a lot of 5,000, the number of support tickets in an hour, the number of mutations in a stretch of DNA. Most of these follow a small set of **named distributions** with known formulas. If you can recognise which one fits, you get the probabilities, mean and variance for free.

## The binomial distribution

Use the **binomial** when all of these hold:

1. A fixed number of trials, `n`.
2. Each trial is a success or a failure.
3. Every trial has the same success probability, `p`.
4. Trials are independent.

Then `X` = number of successes has

```text
P(X = k) = C(n, k) * p^k * (1 - p)^(n - k)      k = 0, 1, ..., n
E[X] = n p        Var(X) = n p (1 - p)
```

`C(n, k)` ("n choose k") counts the ways to place `k` successes among `n` trials. In Python it is `math.comb(n, k)`.

### Worked example: defects

A machine makes parts with a 2% defect rate. In a sample of 50 parts:

```text
P(no defects)   = 0.98^50 ≈ 0.364
P(at most 2)    = P(0) + P(1) + P(2) ≈ 0.922
E[X] = 50 * 0.02 = 1 defect        SD = sqrt(50 * 0.02 * 0.98) ≈ 0.99
```

So about one batch in three has no defects at all, and a sample with 4 or more defects (probability about 1.8%) is a strong signal that the process has drifted. This is how acceptance sampling and lab QC limits are set.

```python
from math import comb

def binomial_pmf(k, n, p):
    return comb(n, k) * p**k * (1 - p) ** (n - k)

print(sum(binomial_pmf(k, 50, 0.02) for k in range(3)))   # 0.9216
```

## The Poisson distribution

Use the **Poisson** for counts of events in a fixed interval of time, space or volume when events happen independently at a steady average rate `λ` (lambda): calls per hour, typos per page, bacterial colonies per plate, photons per detector window.

```text
P(X = k) = e^(-λ) * λ^k / k!        k = 0, 1, 2, ...  (no upper limit)
E[X] = λ          Var(X) = λ
```

### Worked example: a help desk

Calls arrive at an average of 3 per hour.

```text
P(no calls in an hour)  = e^-3 ≈ 0.050
P(at most 4 calls)      = P(0) + ... + P(4) ≈ 0.815
P(5 or more calls)      = 1 - 0.815 ≈ 0.185
```

If one person can handle 4 calls an hour, they will be overloaded in roughly one hour out of five. That is a staffing decision made with one formula.

Rates scale with the interval: 3 calls per hour means `λ = 1.5` for a half-hour window and `λ = 24` for an 8-hour shift.

### Poisson as a binomial limit

When `n` is large and `p` is small, a binomial is almost exactly Poisson with `λ = n p`. For `n = 1000` and `p = 0.002`, `P(X = 2)` is 0.2709 from the binomial and 0.2707 from the Poisson with `λ = 2`. This is why rare-event counts (claims, failures, mutations) are so often modelled as Poisson.

## A check you can always run

Mean equals variance for a Poisson. If your count data has a variance much larger than its mean (common with web traffic and insurance claims), the data is **overdispersed** and a plain Poisson model will understate risk. Analysts then reach for the negative binomial distribution, which you will meet in regression courses.

## Other discrete distributions to know by name

- **Bernoulli**: a single trial (binomial with `n = 1`).
- **Geometric**: number of trials until the first success; mean `1/p`.
- **Negative binomial**: trials until the `r`-th success, or overdispersed counts.
- **Hypergeometric**: successes when sampling **without** replacement from a small finite lot.

## In libraries

`numpy` draws samples (`rng.binomial(50, 0.02, size=1000)`, `rng.poisson(3, size=1000)`), and `scipy.stats` gives exact values (`scipy.stats.binom.pmf(2, 50, 0.02)`, `scipy.stats.poisson.cdf(4, 3)`). Knowing the formulas lets you sanity-check both.

## Common mistakes

- **Using the binomial without independence or a constant p.** Drawing 10 cards from a deck without replacement is hypergeometric, not binomial.
- **Forgetting to rescale λ.** A rate per hour does not apply directly to a 15-minute window.
- **"At least" errors.** `P(X ≥ 5) = 1 - P(X ≤ 4)`, not `1 - P(X ≤ 5)`.
- **Assuming Poisson because the data is a count.** Check that the variance is close to the mean.

## Recap

- Binomial: fixed `n` independent yes/no trials with constant `p`; mean `np`, variance `np(1-p)`.
- Poisson: counts in an interval at rate `λ`; mean and variance both `λ`.
- Large `n`, small `p`: binomial ≈ Poisson with `λ = np`.

## Try it

In the **Do** tab, write the Poisson PMF and CDF and use them to answer a staffing question.


Model counts of successes with the binomial distribution and counts of events in an interval with the Poisson distribution.

Discrete distributions

The Binomial Distribution: Crash Course Statistics #15

An Introduction to the Poisson Distribution

## Why this matters

Lengths, concentrations, temperatures, lifetimes and waiting times are not counts: they can take any value in a range. You cannot list their probabilities one by one, so continuous random variables are described by a **density** instead. Three continuous distributions cover a huge share of everyday work: the **exponential** for waiting times, the **normal** for measurements, and the **gamma** for totals of waiting times and skewed positive quantities.

## Density and CDF

A **probability density function (PDF)** `f(x)` is a curve whose **area** gives probability:

```text
P(a ≤ X ≤ b) = area under f(x) between a and b
total area under f(x) = 1
```

Two consequences surprise people:

- `P(X = 2.5)` is exactly 0 for a continuous variable. Only intervals have probability.
- A density value can be larger than 1 (a very narrow curve is very tall). Density is not probability.

The **cumulative distribution function (CDF)** `F(x) = P(X ≤ x)` is usually what you compute with, because

```text
P(a < X ≤ b) = F(b) - F(a)
P(X > x)     = 1 - F(x)
```

## The exponential distribution

The time between events in a Poisson process (lesson 2) is **exponential** with rate `λ`:

```text
f(x) = λ e^(-λx)        F(x) = 1 - e^(-λx)        x ≥ 0
E[X] = 1/λ              SD(X) = 1/λ
```

**Worked example.** Calls arrive at 3 per hour, so `λ = 3` per hour and the mean gap is 20 minutes. The chance of waiting more than 30 minutes (0.5 h) for the next call is `e^(-3 × 0.5) = e^(-1.5) ≈ 0.223`.

**Worked example.** A component has a mean life of 1,000 hours (`λ = 0.001` per hour). `P(fails before 500 h) = 1 - e^(-0.5) ≈ 0.393`.

The exponential is **memoryless**: `P(X > s + t | X > s) = P(X > t)`. A part that has survived 800 hours is, under this model, as good as new. That is realistic for random shocks and electronic failures but wrong for parts that wear out, which is why reliability engineers also use the Weibull distribution.

## The normal distribution

You met the normal in Statistics Fundamentals. Here is the tool for computing with it exactly. For `X ~ Normal(μ, σ)`, standardise with `z = (x - μ) / σ`, then use the standard normal CDF, which Python can build from the error function:

```text
Φ(z) = 0.5 * (1 + erf(z / sqrt(2)))
```

**Worked example.** Test scores are Normal(100, 15). The share between 85 and 130 is `Φ(2) - Φ(-1) ≈ 0.977 - 0.159 = 0.819`.

```python
import math
from statistics import NormalDist

def normal_cdf(x, mu, sigma):
    return 0.5 * (1 + math.erf((x - mu) / (sigma * math.sqrt(2))))

print(normal_cdf(130, 100, 15) - normal_cdf(85, 100, 15))   # 0.8186
print(NormalDist(100, 15).cdf(130) - NormalDist(100, 15).cdf(85))   # same
```

`statistics.NormalDist` also has `inv_cdf` for the reverse question: "what score marks the top 5%?" (`NormalDist(100, 15).inv_cdf(0.95) ≈ 124.7`).

## The gamma distribution

Add up `k` independent exponential waiting times, each with rate `λ`, and the total follows a **gamma** distribution with shape `k` and rate `λ`:

```text
E[X] = k / λ          Var(X) = k / λ^2
```

The time until the 3rd call at 3 calls per hour has mean `3/3 = 1` hour and SD `sqrt(3)/3 ≈ 0.58` hours. The gamma is also a flexible model for any positive, right-skewed quantity: claim sizes, rainfall totals, reaction times. It allows non-integer shapes too, and with `k = 1` it is the exponential.

## In libraries

`scipy.stats` provides all three with the same interface (`pdf`, `cdf`, `ppf` for the inverse CDF, `rvs` for random draws):

```python
from scipy import stats
stats.expon(scale=1/3).sf(0.5)        # P(X > 0.5) = 0.223; scale is 1/rate
stats.norm(100, 15).cdf(130)
stats.gamma(a=3, scale=1/3).mean()    # 1.0
```

Watch the parameterisation: SciPy and NumPy use **scale = 1/rate** for the exponential and gamma. Mixing up rate and scale is the most common error in this area.

## Common mistakes

- **Treating density as probability.** Only areas are probabilities.
- **Rate vs scale.** `rng.exponential(3)` in NumPy has mean 3, not rate 3.
- **Using the exponential for wear-out.** Memorylessness rarely fits mechanical parts.
- **Assuming everything is normal.** Waiting times and claim sizes are skewed; look at a histogram first.

## Recap

- Continuous variables use a density; probabilities are areas, computed with the CDF.
- Exponential: waiting time between Poisson events, mean `1/λ`, memoryless.
- Normal: standardise, then use `Φ` (via `math.erf` or `NormalDist`).
- Gamma: sum of `k` exponentials, mean `k/λ`, variance `k/λ²`.

## Try it

In the **Do** tab, write the exponential and normal CDFs and use them for a waiting-time and a measurement question.


Use densities and CDFs to work with exponential, normal and gamma distributions for measurements and waiting times.

Continuous distributions

Continuous Probability Distributions - Basic Introduction

The Exponential Distribution: Time Between Poisson Events

Random variables, expectation and variance, and the discrete and continuous distributions (binomial, Poisson, exponential, normal, gamma) that model real counts and waiting times.