Generating CDFs
A named distribution is usually the limit of a simple random process, not an arbitrary formula someone wrote down. In each of the nine mini-labs below, the group builds the machine first — draw, count, sort, exchange, whatever the story calls for — and only afterward checks the machine's own ECDF against the named distribution's CDF from scipy.stats. Success means the simulation's own staircase lines up with the textbook curve.
labslop/06_generating_cdfs/lab_00_generating_cdfs.ipynb. Every simulation below was re-run to confirm the reported quantiles; no external data required — every mini-lab generates its own.
- Think · 2 minPick any one of the nine distributions listed above by name. Before reading its story below: what physical or data-generating process would you guess produces it?
I've done the Think step — reveal Pair & Share
- Pair · 3 minCompare guesses, then jump to that section and see how close the actual story is.
- Share · 2 minAs a table, agree on which of the nine stories was the most surprising match to its distribution's name.
1 · A uniform distribution from a finite grid
Story: a machine chooses one of \(K\) equally spaced values. When \(K\) is small, each value has visible probability. As \(K\) grows, the probability of landing on any exact value goes to zero, but probabilities over intervals stay stable.
Algorithm: build the grid \(\{1/K, 2/K, \dots, K/K\}\), sample grid points with equal probability, and compare the resulting ECDF to stats.uniform(0,1).cdf.
draws = rng.integers(1, K + 1, size=B) / KAt \(K=5\) the ECDF is visibly a five-step staircase; by \(K=1{,}000\) the atom probability \(1/K = 0.001\) is small enough that the staircase is indistinguishable by eye from the smooth Uniform(0,1) CDF line. This is the same "grid gets finer" idea behind the PDF-as-limit-of-histogram puzzle from Class 06's lecture, run forward instead of backward.
2 · The normal distribution from many small additions
Story: many unrelated small shocks add together. After centering and scaling, the sum becomes normal — the Central Limit Theorem, built by hand instead of assumed.
Algorithm (uniform shocks): draw \(N\) values from Uniform\([-\sqrt3, \sqrt3]\) (mean 0, variance 1), average them, multiply by \(\sqrt N\), repeat \(B\) times.
shocks = rng.uniform(-a, a, size=(B, N)) # a = sqrt(3)
z_uniform = np.sqrt(N) * shocks.mean(axis=1)With \(B=30{,}000\) draws of \(N=40\) shocks each, the simulated quantiles already track the standard normal closely:
| quantile p | simulation | Normal(0,1) (scipy) |
|---|---|---|
| 0.10 | −1.277 | −1.282 |
| 0.25 | −0.684 | −0.674 |
| 0.50 | −0.007 | 0.000 |
| 0.75 | 0.678 | 0.674 |
| 0.90 | 1.289 | 1.282 |
The same convergence shows up starting from Bernoulli trials instead of uniform shocks — count successes in \(N\) yes/no trials, subtract the expected count \(Np\), divide by the standard deviation \(\sqrt{Np(1-p)}\) — the exact standardization behind the binomial-to-normal approximation.
3 · An exponential distribution from a failure clock
Story: a device is alive. In each tiny time interval of length \(\Delta\), it dies with probability \(\lambda \Delta\). As the time step shrinks, the lifetime distribution becomes exponential.
Algorithm: pick a death rate \(\lambda\) and a small step \(\Delta\) with \(\lambda\Delta \lt 1\); for each simulated device, draw the first period death occurs in (a geometric random variable), then convert periods to time by multiplying by \(\Delta\).
p_death = lam * Delta
periods_until_death = rng.geometric(p_death, size=B)
lifetimes = periods_until_death * Delta| quantile p | simulation (λ=0.8) | Exponential(scale=1/λ) |
|---|---|---|
| 0.10 | 0.14 | 0.132 |
| 0.25 | 0.36 | 0.360 |
| 0.50 | 0.88 | 0.866 |
| 0.75 | 1.74 | 1.733 |
| 0.90 | 2.88 | 2.878 |
4 · A Poisson distribution from soccer goals
Story: goals are rare in any tiny slice of a match, but there are many slices. Count the goals across a full 90-minute match.
Algorithm: split the match into many tiny intervals; in each, a goal occurs with probability \(\mu / \text{intervals}\); sum across intervals for the match total.
p_goal = mu / intervals
goals = rng.binomial(intervals, p_goal, size=B) # mu = 2.7, intervals = 900| goals ≤ | simulation P | Poisson(2.7) P |
|---|---|---|
| 0 | 0.0677 | 0.0672 |
| 1 | 0.2485 | 0.2487 |
| 2 | 0.4954 | 0.4936 |
| 3 | 0.7174 | 0.7141 |
| 4 | 0.8645 | 0.8629 |
This is exactly the Poisson-as-limit-of-binomial construction: a huge number of tiny, independent, low-probability trials (goal chances per second) collapses to a Poisson count over the whole match.
5 · A Pareto distribution from Zipf's law for cities
Story: rank cities from largest to smallest. Under Zipf's law, size is proportional to \(1/\text{rank}\). Choose a city uniformly by rank, and its size follows a Pareto distribution in the tail.
sampled_ranks = rng.integers(1, K + 1, size=B)
city_sizes = (K / sampled_ranks) ** alpha| quantile p | simulation | Pareto(b=1, scale=1) |
|---|---|---|
| 0.10 | 1.11 | 1.11 |
| 0.25 | 1.33 | 1.33 |
| 0.50 | 1.99 | 2.00 |
| 0.75 | 3.98 | 4.00 |
| 0.90 | 9.90 | 10.00 |
6 · A Boltzmann distribution from energy exchange
Story: many particles share a fixed pool of energy quanta and randomly pass energy around. After enough random exchanges, one particle's energy looks like a discrete exponential — the Boltzmann distribution of statistical mechanics.
Algorithm: assign \(Q\) quanta randomly to \(M\) particles; repeatedly pick a donor and receiver at random and move one quantum if the donor has any; after burn-in, sample particle energies periodically to reduce autocorrelation.
With \(M=500\) particles, \(Q=2{,}000\) quanta (mean energy 4 per particle), and 500,000 exchange steps: the matched Boltzmann parameter is \(\lambda = \log(1 + 1/\bar{E}) \approx 0.223\), and the quantiles land almost exactly on the fitted distribution:
| quantile p | simulation | Boltzmann(λ≈0.223) |
|---|---|---|
| 0.10 | 0 | 0 |
| 0.25 | 1 | 1 |
| 0.50 | 3 | 3 |
| 0.75 | 6 | 6 |
| 0.90 | 10 | 10 |
7 · A Gumbel distribution from maxima
Story: take the maximum of many independent draws. For many common light-tailed distributions, the centered maximum converges to a Gumbel distribution — the extreme-value analog of the CLT.
maxima = rng.exponential(scale=1, size=(B, n)).max(axis=1)
centered_maxima = maxima - np.log(n)| quantile p | simulation (n=1,000) | Gumbel (scipy) |
|---|---|---|
| 0.10 | −0.826 | −0.834 |
| 0.25 | −0.311 | −0.327 |
| 0.50 | 0.375 | 0.367 |
| 0.75 | 1.234 | 1.246 |
| 0.90 | 2.220 | 2.250 |
8 · A Beta distribution from order statistics
Story: sort a small sample of uniform draws. The location of the \(k\)th smallest value has a Beta distribution.
uniform_samples = rng.uniform(0, 1, size=(B, n))
kth_smallest = np.sort(uniform_samples, axis=1)[:, k - 1] # n=12, k=4| quantile p | simulation | Beta(4, 9) |
|---|---|---|
| 0.10 | 0.153 | 0.154 |
| 0.25 | 0.215 | 0.216 |
| 0.50 | 0.297 | 0.298 |
| 0.75 | 0.389 | 0.389 |
| 0.90 | 0.475 | 0.475 |
This is exactly Class 07's "order statistics & auctions" group-work exercise, given its closed form: the \(k\)th smallest of \(n\) uniforms is Beta(\(k\), \(n+1-k\)).
9 · A lognormal distribution from binomial asset pricing
Story: an asset price repeatedly moves up or down by small multiplicative factors. As the time step shrinks, the terminal price becomes lognormal — the discrete binomial-tree model converging to geometric Brownian motion.
u, d = np.exp(sigma*np.sqrt(dt)), np.exp(-sigma*np.sqrt(dt))
p_up = (np.exp(mu*dt) - d) / (u - d)
up_counts = rng.binomial(N, p_up, size=B)
final_prices = S0 * (u ** up_counts) * (d ** (N - up_counts))| quantile p | simulation | Lognormal (GBM) |
|---|---|---|
| 0.10 | $75.32 | $76.21 |
| 0.25 | $88.16 | $88.70 |
| 0.50 | $106.50 | $105.00 |
| 0.75 | $124.67 | $124.28 |
| 0.90 | $145.93 | $144.65 |
This is the same binomial-tree asset-pricing model from Class 04's lecture, run to its \(N=252\)-step (daily, one trading year) limit — the same construction, the same lognormal-limit result Black-Scholes assumes outright.
Every one of these nine stories starts from something crude and discrete (a grid, a coin flip, a random exchange) and ends at a smooth named distribution. What is the common mechanism — the one sentence that would explain all nine convergences at once, not just one of them?
Show one way to answer
Take many small, mostly-independent pieces of randomness and combine them (sum, count, max, or repeated exchange) — the combining operation itself erases the details of any one piece, leaving only a few aggregate quantities (a rate, a mean, a variance) to determine the limiting shape. That is the Law of Large Numbers and the Central Limit Theorem (and their extreme-value and combinatorial cousins), doing the same job nine different times.
FAQ
Why compare quantiles instead of just eyeballing the ECDF-vs-CDF plot?
A plot is the fastest way to see gross disagreement, but quantiles give an exact number to check — if the simulated 50th percentile and the theoretical median differ by more than sampling noise should allow, something in the algorithm is wrong, not just "a little off in the picture."
Why does the Boltzmann simulation need a burn-in period and thinning?
The random-exchange process starts from an arbitrary initial assignment of quanta to particles, which isn't yet a "typical" configuration — burn-in (discarding the first 100,000 of 500,000 steps) lets the process forget its starting point. Thinning (keeping only every 2,000th step afterward) reduces the correlation between consecutive snapshots, since two exchanges apart is not enough randomness for two snapshots to look independent.