PDFs, CDFs, and inverse CDFs with scipy.stats — and a synthetic motor portfolio with realistic, long-tailed claims.
scipy.stats “frozen” distributions consistentlyCompanion notebook: lesson-07-continuous-losses.ipynb — worked examples and self-checking exercises.
scipy.stats ToolkitOne consistent interface for over a hundred distributions.
import numpy as np
from scipy import stats
N = stats.norm(loc=100, scale=15) # "freeze" the parameters
N.pdf(100) # 0.026596 density at 100
N.cdf(130) # 0.97725 P(X ≤ 130)
N.sf(130) # 0.02275 P(X > 130) = 1 - cdf
N.ppf(0.975) # 129.3995 inverse CDF: the 97.5% quantile
N.mean(), N.std() # (100.0, 15.0)
rng = np.random.default_rng(11)
N.rvs(size=3, random_state=rng) # [100.51 120.4 118.37]Use sf(x) rather than 1 - cdf(x) for far-tail probabilities: it keeps precision where the CDF rounds to 1.
norm for lognorm, nothing else changesGenerator as random_state — Lesson 6 reproducibility still appliesN.cdf([90, 130])pdf(x)Density: relative likelihood near x. Areas under it are probabilities.
For discrete distributions: pmf(k) = P(X = k)
cdf(x)P(X ≤ x). Rises from 0 to 1.
sf(x) = 1 − CDF: the survival / exceedance probability
ppf(p)The inverse CDF: the value x with P(X ≤ x) = p.
A quantile. VaR at 99.5% is ppf(0.995).
“What’s the chance a claim exceeds 50,000?” → sf(50_000). “What claim size is exceeded only 1 time in 200?” → ppf(0.995). Most actuarial questions are one of these two.
Getting the parameters right is most of the battle.
E = stats.expon(scale=5000) # scale = mean = 1/λ
E.mean() # 5000.0
E.cdf(5000) # 0.632121 1 - e^{-1}
E.ppf(0.5) # 3465.74 median = mean × ln 2
E.ppf(0.99) # 23025.85scipy uses scale = 1/λ, not the rate λ. stats.expon(5000) without the keyword sets loc — a shift — instead. Always name your parameters.
scale# If ln X ~ Normal(mu, sigma):
# s = sigma, scale = exp(mu)
sigma = 1.0
mu = np.log(5000) - sigma**2 / 2 # chosen so the mean is 5000
L = stats.lognorm(s=sigma, scale=np.exp(mu))
L.mean() # 5000.0
L.median() # 3032.65 median < mean: right-skewed
L.std() # 6554.16
L.ppf(0.995) # 39855.77| Math | scipy.stats.lognorm |
|---|---|
| σ (sd of ln X) | s=sigma |
| μ (mean of ln X) | scale=np.exp(mu) |
| median | = scale |
| mean | exp(μ + σ²/2) |
The single most common scipy bug in actuarial work: passing μ as scale, or σ² as s. Check L.mean() and L.median() against what you intended.
m, s = 42, 14 # underlying normal
lo, hi = 18, 85 # licensing and cut-off ages
# Bounds must be in STANDARD units: (bound - m) / s
a, b = (lo - m) / s, (hi - m) / s # (-1.7143, 3.0714)
T = stats.truncnorm(a, b, loc=m, scale=s)
T.mean(), T.std() # (43.2923, 12.6417)
T.cdf(25) # 0.0723 under-25s
T.ppf([0, 0.5, 1]) # [18. 42.74 85. ]
stats.truncnorm(18, 85, loc=42, scale=14).mean()
# 294.77 — bounds read as 18 and 85 standard deviations!m — trimming the left tail more than the right pulls it up to 43.29truncnorm’s a and b are in standard deviations from loc, not in the units of the data. Pass raw ages and you get nonsense, silently.
If U ~ Uniform(0, 1), then F−1(U) has CDF F. So any distribution with a PPF can be simulated from uniforms:
rng = np.random.default_rng(5)
u = rng.random(5)
# [0.805 0.8079 0.5153 0.2858 0.0539]
E.ppf(u)
# [8173.85 8249.76 3621.39 1682.97 277.2 ]This is essentially what .rvs() does for you.
Why learn it? It links a simulated loss to a percentile. Feed the same uniforms into two models and you compare them scenario by scenario — the foundation of stress testing and of dependence modeling with copulas (Course 3).
High u → large loss. Sort your uniforms and you’ve sorted your losses.
Exponential and Lognormal (σ = 1), both with mean 5,000. Probability a claim exceeds k × the mean:
| k | Exponential | Lognormal | Ratio |
|---|---|---|---|
| 2 | 1.35 × 10⁻¹ | 1.16 × 10⁻¹ | 0.9 |
| 5 | 6.74 × 10⁻³ | 1.75 × 10⁻² | 2.6 |
| 10 | 4.54 × 10⁻⁵ | 2.54 × 10⁻³ | 56 |
| 20 | 2.06 × 10⁻⁹ | 2.36 × 10⁻⁴ | 115,000 |
Skewness: exponential 2.0, lognormal 6.18.
Mean excess e(d) = E[X − d | X > d], from 1,000,000 lognormal draws:
| Threshold d | Exponential | Lognormal |
|---|---|---|
| 5,000 | 5,000 | 6,201 |
| 20,000 | 5,000 | 11,746 |
| 50,000 | 5,000 | 19,533 |
A mean excess that grows with the threshold is the signature of a heavy tail: the bigger a claim already is, the bigger you should expect it to get. That’s why excess-of-loss layers are hard to price.
Build a synthetic motor portfolio of 10,000 drivers — each with an age, an annual mileage, a claim count, and claim amounts after a deductible and a policy limit.
One seeded generator for everything, and draw in the order shown on the next slide — then your numbers will match the answer table.
rng = np.random.default_rng(2026)
n = 10_000
age_dist = stats.truncnorm((18 - 42) / 14, (85 - 42) / 14,
loc=42, scale=14)
ages = np.floor(age_dist.rvs(size=n, random_state=rng)).astype(int)
mileage = stats.lognorm(s=0.5, scale=12_000).rvs(size=n, random_state=rng)
base_freq = 0.08
age_factor = np.where(ages < 25, 1.8, np.where(ages >= 70, 1.3, 1.0))
freq = base_freq * (mileage / 12_000) ** 0.5 * age_factor
claims = rng.poisson(freq) # claim count per driver
sev_dist = stats.lognorm(s=1.2, scale=4_000)
deductible, limit = 500, 50_000
ground_up = sev_dist.rvs(size=claims.sum(), random_state=rng)
paid = np.minimum(np.maximum(ground_up - deductible, 0), limit)| Check | Value |
|---|---|
| Mean age · under-25 share | 43.06 · 6.55% |
| Median · mean mileage | 12,000 · 13,652 |
| Mean frequency | 0.0876 |
| Total claims | 852 |
| Claims per driver 0 / 1 / 2 / 3 / 4 | 9,189 / 773 / 36 / 1 / 1 |
| Below deductible · hit the limit | 34 · 15 |
| Mean paid per claim | 7,132.35 |
| Loss cost per driver | 607.68 |
E[min(X, 50,500)] − E[min(X, 500)] using sev_dist.expect(...) → 7,114.09. Within a standard error?ages: check the hard edges at 18 and 85Move the generator into src/foundations/simulation.py as a function:
def simulate_portfolio(n: int, rng: np.random.Generator,
deductible: float = 500,
limit: float = 50_000) -> dict:
"""Simulate drivers and their paid claims.
Returns a dict with arrays: ages, mileage,
claims (count per driver), paid (per claim).
"""
...This function is the frequency–severity engine of the next lesson’s capstone. Commit it.
scipy.stats distributions and name every parameterpdf, cdf, sf, and ppftruncnorm its bounds in standard unitsnp.minimum/np.maximumsimulate_portfolio() lives in src/foundations/Modules, docstrings, and pytest — then a Monte Carlo library that computes a portfolio’s VaR and TVaR.