Distributions and Sampling

A small incident assistant helps engineers with CI failures, rollbacks, access requests, and log triage. A launch review asks for 5,000 synthetic requests before traffic arrives. Each request contributes four measurements: it either resolves or escalates, takes one route, makes a nonnegative number of tool calls, and takes a positive number of seconds.

Now choose how to draw those four values. A wrong choice can make an average look reasonable while hiding retry bursts, producing impossible negative latency, or inventing a route that never appears in the taxonomy. Before trusting a simulator, name the claim behind each draw and the trace that could falsify it.

In Statistics and Uncertainty, 16 abuse labels out of 100 reviewed flags became a point estimate plus an interval. That asked how much a finite sample supports one unknown rate. Simulation reverses the direction: start with a process, draw new samples, and ask whether the resulting traces look like the system you intend to model.

A probability distribution specifies the allowed values of a random quantity and how probability is assigned to them. A launch review fails if a draw can be negative, invent a route, or hide a retry burst. Write down the allowed values first, then check the frequencies and tails the simulator actually produces.

The numbers below are planning assumptions. A real launch would estimate them from logged, reviewed traffic. The standard-library examples need Python 3.12 or later because they use Random.binomialvariate. Record the interpreter version along with the seed; distribution-sampling algorithms can change between versions.^{[1]Reference 1random: Generate pseudo-random numbershttps://docs.python.org/3/library/random.html}

Four measurements from one request

The first simulation model column is a claim about shape, not a default. If traces disagree, you revise the model instead of defending the convenient formula.

Discrete values versus continuous time

The first split is between discrete and continuous quantities.

For a discrete quantity you can ask for the probability of one exact value, such as P(tool_calls = 2). For a continuous response time, one infinitely precise instant has probability zero under the model. Ask for a range or a percentile instead: P(time > 15 seconds) or the p95 response time.

The allowed values are the distribution's support. Route labels need a finite set of names, tool calls need nonnegative integers, and response time needs nonnegative seconds. A measured zero may mean a duration was rounded below the clock's resolution; a negative duration is a different problem.

Before sampling, write down four things:

1. What values are allowed?
2. Which parameter sets the rate, center, or spread?
3. Which tail or threshold would hurt users or cost?
4. What observation would prove this model is a poor fit?

For example, ten tool calls may be rare under a baseline but common during retry loops. That observation should change the model, not just the seed.

PMF, PDF, and CDF

Use the kind of quantity to choose the probability object. For the binary resolution variable R, the two masses are P(R = 0) = 0.28 and P(R = 1) = 0.72. A probability mass function (PMF) lists those probabilities for each discrete value. The masses add to 1.

Response time needs different language. A probability density function (PDF) describes how densely probability is spread across possible times; area over an interval gives its probability. Its height isn't the probability of one exact time and may exceed 1. For instance, a uniform duration between 0 and 0.5 seconds has density 2 per second, but total area 2 × 0.5 = 1.

A cumulative distribution function (CDF) gives the chance of finishing by a deadline. Let T mean response time and t the chosen deadline:

Read the CDF as an operations question: F_T(10) is the chance a response finishes within 10 seconds, while 1 - F_T(15) is the chance it exceeds 15 seconds. For the smooth latency model below, p95 is the time where the CDF reaches 0.95. Discrete distributions have CDFs too, but those rise in steps and can jump past 0.95.^{[2]Reference 2Machine Learning: A Probabilistic Perspective.https://probml.github.io/pml-book/book0.html}

A chart or API that says "probability" is easy to misread. A PMF bar is a probability. A PDF curve needs an interval before it becomes one. A CDF's height is already accumulated probability.

The probability chapter treated a distribution as a function over a known sample space. Naming a response-time family still doesn't produce a duration. How does a computer turn one uniform random number into a simulated time?

Inverse CDF transform: how computers sample continuous variables

Start from uniform draws, $U \sim \text{Uniform}(0, 1)$. A pseudo-random generator supplies finite-precision numbers that approximate this mathematical input. It can also supply random bits and integers.

The inverse CDF transform, or inverse-transform sampling, maps uniform draws through a distribution's quantile function. If $F$ is strictly increasing over its support, this is its ordinary inverse. A continuous CDF can also have flat regions, so the general definition is $F^{-1}(u) = \inf\{t : F(t) \geq u\}$, the first point where the CDF reaches the chosen probability. This generalized inverse works for discrete distributions too.

With an ideal uniform draw and the exact quantile function, $T = F^{-1}(U)$ follows the target distribution:

The probability integral transform runs in the other direction: for a draw $T$ from a continuous CDF $F$, the value $F(T)$ is uniform. The two constructions are related, but they aren't the same operation. Numerical quantile implementations and finite-precision RNGs approximate the ideal sampling statement.SciPy's numerical inversion documentation describes this distinction in practice.

To see how that works analytically, consider an exponential service duration with completion rate $\mu$:

Because $1 - u$ is also distributed uniformly on $(0, 1)$, code can draw $u \sim \text{Uniform}(0, 1)$ and compute $s = -\ln(u) / \mu$. For a lognormal variable with median $e^\mu$ and log-spread $\sigma$, you invert standard normal noise: draw $u$, find the standard normal quantile $z = \Phi^{-1}(u)$, and exponentiate: $t = \exp(\mu + \sigma z)$.

View larger

Notice the geometry in the plot: where the CDF is steep (near the 6-second median), a wide swath of uniform draws projects into a narrow slice of time. Where the CDF flattens out (the right tail), a tiny 5% sliver of uniform draws ($u \in [0.95, 1.0]$) stretches out across all latencies beyond 14.8 seconds.

Before running the snippet below, predict the output: 5,000 inverted uniform draws should land near the theoretical 0.05-second exponential mean, the 6.0-second lognormal median, and the 14.83-second p95 threshold. Python's random() includes zero, while the normal quantile requires 0 < u < 1; redraw that endpoint instead of passing it into inv_cdf.

from math import ceil, exp, log, log1p
from random import Random
from statistics import NormalDist, median

rng = Random(42)
rate = 20.0
def open_uniform(rng: Random) -> float:
    u = rng.random()
    while u == 0.0:
        u = rng.random()
    return u  # random() excludes 1; redrawing zero makes this (0, 1).

uniform_draws = [open_uniform(rng) for _ in range(5_000)]
exp_samples = [-log1p(-u) / rate for u in uniform_draws]

mu = log(6.0)
sigma = 0.55
standard_normal = NormalDist(0, 1)
lognormal_samples = [
    exp(mu + sigma * standard_normal.inv_cdf(u)) for u in uniform_draws
]

print(f"exponential mean: {sum(exp_samples) / len(exp_samples):.3f}s (target: {1/rate:.3f}s)")
print(f"lognormal median: {median(lognormal_samples):.2f}s (target: 6.00s)")
p95_index = ceil(0.95 * len(lognormal_samples)) - 1  # Nearest-rank quantile.
print(f"lognormal p95:    {sorted(lognormal_samples)[p95_index]:.2f}s (target: 14.83s)")

exponential mean: 0.050s (target: 0.050s)
lognormal median: 6.03s (target: 6.00s)
lognormal p95:    14.85s (target: 14.83s)

Bernoulli: resolved or escalated

Begin with the smallest distribution. Let R = 1 if the assistant resolves a request without a human, and R = 0 if it escalates. The planning value is p = 0.72.

The expected value is:

Because R is encoded as 0 or 1, its long-run average is the resolution rate. Calculate that weighted average directly:

p_resolved = 0.72
outcomes = {"escalated": 0, "resolved": 1}
probabilities = {"escalated": 1 - p_resolved, "resolved": p_resolved}

expected_value = sum(
    outcomes[name] * probabilities[name] for name in outcomes
)

print(f"P(resolved):  {probabilities['resolved']:.2f}")
print(f"P(escalated): {probabilities['escalated']:.2f}")
print(f"E[R]:         {expected_value:.2f}")

assert expected_value == p_resolved

P(resolved):  0.72
P(escalated): 0.28
E[R]:         0.72

One request is now a coin with a known bias. The next question is what happens when you count many of those coins.

Binomial: resolutions in a 100-request batch

If you count resolutions across n = 100 independent requests that all share p = 0.72, the count X has a binomial distribution. A single request is Bernoulli; the sum of these 100 outcomes is binomial.^{[3]Reference 3NIST/SEMATECH e-Handbook: Binomial Distributionhttps://www.itl.nist.gov/div898/handbook/eda/section3/eda366i.htm}

Its expected count is np = 72, and its variance is np(1-p) = 20.16. The count can land anywhere from 0 to 100, even though 72 is the center of the model.

One particular arrangement of 72 resolutions and 28 escalations has probability 0.72^72 × 0.28^28. There are many arrangements with that same count. The binomial coefficient $\binom{100}{72}$ counts them, giving:

The binomial model is a sum of independent Bernoulli trials. If requests arrive in correlated incident waves, or if each request has a different resolution probability, this simple count model needs a richer design.

Before running the snippet, predict its shape: six rates should scatter around .72, not all equal .72. The code checks the formula, then draws six independent 100-request batches while holding p = 0.72 fixed. Random.binomialvariate is the standard-library binomial sampler.

from math import comb
from random import Random

n = 100
p = 0.72
k = 72
p_exact = comb(n, k) * (p**k) * ((1 - p) ** (n - k))

print(f"P(X={k}) = {p_exact:.3f}")
print(f"E[X] = {n * p:.0f}")
print(f"Var(X) = {n * p * (1 - p):.2f}")

rng = Random(12)
resolved_per_batch = [rng.binomialvariate(n, p) for _ in range(6)]
rates = [count / n for count in resolved_per_batch]

print("resolved counts:", resolved_per_batch)
print("batch rates:    ", [round(rate, 2) for rate in rates])
print(f"range:           {min(rates):.2f} to {max(rates):.2f}")

P(X=72) = 0.089
E[X] = 72
Var(X) = 20.16
resolved counts: [72, 70, 75, 69, 71, 78]
batch rates:     [0.72, 0.7, 0.75, 0.69, 0.71, 0.78]
range:           0.69 to 0.78

A confidence interval from the last chapter uses a finite sample to quantify uncertainty about an unknown rate. A distribution does the opposite: it assumes a process, then generates the samples you might see. The interval asks what the data support; the simulator asks what the assumption would produce.

Binary resolution is only one label. Real traffic is split across several routes.

Categorical: one route per request

An incident-assistant request isn't always binary. It may be exactly one of several route labels. For this planning mix:

The four probabilities sum to 1.00 because every request receives one route in this simplified taxonomy. A categorical distribution draws one of those labels according to the supplied weights.

One request gives one categorical outcome. If 500 requests are independent and use this same route mix, the four counts together follow a multinomial distribution. Their sum must be 500, so the counts aren't independent of one another. This sampler draws requests one at a time, then counts labels.

Predict before running it: the observed counts should be near 225, 100, 100, and 75, but they don't need to match those expected counts exactly. Sampling variation is the behavior we want to inspect.

from random import Random

intents = ["ci_failure", "rollback", "access_request", "log_triage"]
probabilities = [0.45, 0.20, 0.20, 0.15]
requests = 500

assert abs(sum(probabilities) - 1.0) < 1e-12

rng = Random(21)
sampled = rng.choices(intents, weights=probabilities, k=requests)
expected = {
    intent: int(requests * probability)
    for intent, probability in zip(intents, probabilities)
}
observed = {intent: sampled.count(intent) for intent in intents}

print("expected counts:", expected)
print("sampled counts: ", observed)
print("all labels known:", set(sampled).issubset(set(intents)))

expected counts: {'ci_failure': 225, 'rollback': 100, 'access_request': 100, 'log_triage': 75}
sampled counts:  {'ci_failure': 218, 'rollback': 105, 'access_request': 106, 'log_triage': 71}
all labels known: True

The sampled counts needn't equal the expectations. That's ordinary variation. An unknown label such as database_restore would be a different problem: the taxonomy or the sampler contract is wrong. choices accepts relative weights, so [45, 20, 20, 15] would encode the same distribution; our explicit probabilities make the planning assumptions easier to inspect.^{[1]Reference 1random: Generate pseudo-random numbershttps://docs.python.org/3/library/random.html}

Routes tell you what kind of request arrived. Tool calls tell you how expensive that request was to serve.

Poisson counts, then a bursty counterexample

Tool calls are counts. A simple baseline is a Poisson distribution with parameter lambda, written $\lambda$. If lambda = 2.2, the simulated agent makes 2.2 tool calls per request on average.

With a mean of 2.2, a Poisson model predicts about 26.8% of requests will make exactly two calls. To calculate that value, let C denote the count and k one nonnegative integer. In the formula below, k! means the factorial (2! = 2 × 1), and e is the base of natural logarithms:

With $\lambda = 2.2$ and $k = 2$:

The loop below prints a few more masses, including the expensive tail P(C > 5).

from math import exp, factorial

lam = 2.2

def poisson_probability(k: int) -> float:
    if k < 0:
        raise ValueError("k must be nonnegative")
    return lam**k * exp(-lam) / factorial(k)

for k in range(5):
    print(f"P(calls = {k}): {poisson_probability(k):.3f}")

share_above_5 = 1 - sum(poisson_probability(k) for k in range(6))
print(f"P(calls > 5): {share_above_5:.3f}")

try:
    poisson_probability(-1)
except ValueError as error:
    print(error)

P(calls = 0): 0.111
P(calls = 1): 0.244
P(calls = 2): 0.268
P(calls = 3): 0.197
P(calls = 4): 0.108
P(calls > 5): 0.025
k must be nonnegative

The Poisson model makes a stronger claim than "counts can't be negative." Its mean and variance are both $\lambda$.^{[4]Reference 4NIST/SEMATECH e-Handbook: Poisson Distributionhttps://www.itl.nist.gov/div898/handbook/eda/section3/eda366j.htm} A Poisson process produces such counts in a fixed exposure window when events arrive independently at a constant rate. Tool calls within a request aren't automatically such a process: the agent often decides its next call based on a previous result. Here Poisson is a candidate count distribution, not a fact implied by the data type.

An agent can break that assumption. A difficult rollback may trigger a deploy lookup, a runbook lookup, a retry, and a handoff together. Those bursts create overdispersion, where variance is much larger than the mean.

Both streams in the next fixture have a mean near 2.2. The mixed stream uses ordinary requests 88% of the time and hard requests 12% of the time. Its exact expected count is 0.88 × 1.2 + 0.12 × 10 = 2.256, close to the baseline's 2.2.

Predict what a good Poisson fit would show: a variance-to-mean ratio near 1 and a share above five calls near .025. The transparent product-based sampler below is restricted to small rates for this exercise. For general-purpose simulation, use a tested library sampler such as NumPy's Generator.poisson, not this loop at large rates.

from math import exp
from random import Random

def sample_poisson(rng: Random, lam: float) -> int:
    if not 0 <= lam <= 20:
        raise ValueError("this teaching sampler requires 0 <= lambda <= 20")
    threshold = exp(-lam)
    count = 0
    product = 1.0
    while True:
        product *= rng.random()
        if product <= threshold:
            return count
        count += 1

def mean(xs: list[float]) -> float:
    return sum(xs) / len(xs)

def pop_variance(xs: list[float]) -> float:
    center = mean(xs)
    return sum((x - center) ** 2 for x in xs) / len(xs)

rng = Random(31)
steady = [sample_poisson(rng, 2.2) for _ in range(5_000)]
bursty = []
for _ in range(5_000):
    hard_request = rng.random() < 0.12
    bursty.append(sample_poisson(rng, 10.0 if hard_request else 1.2))

def report(name: str, counts: list[int]) -> None:
    average = mean(counts)
    variance = pop_variance(counts)
    above_five = sum(count > 5 for count in counts) / len(counts)
    print(
        f"{name:>6}: mean={average:.2f} variance={variance:.2f} "
        f"variance/mean={variance / average:.2f} share>5={above_five:.3f}"
    )

report("steady", steady)
report("bursty", bursty)

steady: mean=2.23 variance=2.20 variance/mean=0.98 share>5=0.025
bursty: mean=2.29 variance=10.66 variance/mean=4.65 share>5=0.115

The mixed sample has both excess variance and a larger expensive-tail share. It doesn't identify a retry storm: this simulator chooses request types independently and has no temporal clustering. Overdispersion alone doesn't determine every tail probability, and a tiny sample can have an unstable variance estimate. Inspect the count frequencies and the actual threshold you care about. Also separate request types or incident waves before deciding whether one distribution should describe them all.

Counts are only half of the cost story. The other half is how long the request occupies a worker.

Lognormal response time

End-to-end response time can't be negative. Some requests can also take much longer than the typical request, forming a right tail. A lognormal distribution is one candidate for positive, right-skewed values; positivity alone doesn't establish that it's the right family.

A normal distribution is the familiar symmetric bell-shaped model, described by a mean and standard deviation. A variable is lognormal when its natural logarithm is normal. For this example, take the logarithm of the numerical duration expressed in seconds.

random.lognormvariate(mu, sigma) takes the mean and standard deviation of those log values, not of the durations. Set mu = log(6) for a six-second median and sigma = 0.55 for log-space spread. This gives a mean near 6.98 seconds. The conversion rules are:^{[5]Reference 5NIST/SEMATECH e-Handbook: Lognormal Distributionhttps://www.itl.nist.gov/div898/handbook/eda/section3/eda3669.htm}

To find F_T(10), measure how many log-space standard deviations 10 seconds lies above the center. This standardized distance is a z-score:

The normal CDF with mean zero and standard deviation one gives about 0.823 at this z-score. Thus roughly 82% finish within 10 seconds. Its 95th-percentile z-score is about 1.645, so the model's p95 is exp(log(6) + 0.55 × 1.645) ≈ 14.83 seconds. A finite sample needn't hit either value exactly.

Predict before sampling: the median should stay near 6 seconds, p95 should land near 14.8 seconds, and no draw should be nonpositive. The output below checks those three consequences.

⚠️ Common mistake: Passing 6.0 as mu does not make the median six seconds. It makes the median $e^{6} \approx 403$ seconds. NumPy's Generator.lognormal(mean, sigma) uses the same log-space mean.

from math import log
from random import Random
from statistics import NormalDist, fmean, median

def percentile(values: list[float], q: float) -> float:
    ordered = sorted(values)
    position = (len(ordered) - 1) * q / 100
    lower = int(position)
    upper = min(lower + 1, len(ordered) - 1)
    return ordered[lower] + (ordered[upper] - ordered[lower]) * (position - lower)

median_seconds = 6.0
sigma = 0.55
mu = log(median_seconds)

rng = Random(44)
response_seconds = [rng.lognormvariate(mu, sigma) for _ in range(5_000)]

print(f"target median seconds: {median_seconds:.2f}")
print(f"sample median seconds: {median(response_seconds):.2f}")
print(f"sample mean seconds:   {fmean(response_seconds):.2f}")
print(f"sample p95 seconds:    {percentile(response_seconds, 95):.2f}")
print(f"nonpositive values:    {sum(t <= 0 for t in response_seconds)}")
print(f"model P(T <= 10):      {NormalDist(mu, sigma).cdf(log(10)):.3f}")

target median seconds: 6.00
sample median seconds: 5.94
sample mean seconds:   6.94
sample p95 seconds:    14.80
nonpositive values:    0
model P(T <= 10):      0.823

The median describes a typical request. p95 tells you about a noticeably slow slice. Reporting only the mean would hide how long the slow responses take. These lognormal draws stayed positive. What happens if a familiar bell curve, whose support includes negative numbers, is used for a duration instead?

A normal model can go negative

The normal distribution's support extends across the whole real line. It can model signed quantities, such as a measurement error around zero. It can also approximate a positive quantity when the chance of a negative draw is negligible for the task, but that needs checking.

With mean 8 seconds and standard deviation 6 seconds, negative response times are far from negligible. Generate 5,000 durations and check support before accepting the plausible-looking mean or p95:

from random import Random
from statistics import fmean

def percentile(values: list[float], q: float) -> float:
    ordered = sorted(values)
    position = (len(ordered) - 1) * q / 100
    lower = int(position)
    upper = min(lower + 1, len(ordered) - 1)
    return ordered[lower] + (ordered[upper] - ordered[lower]) * (position - lower)

rng = Random(44)
bad_response_seconds = [rng.gauss(8.0, 6.0) for _ in range(5_000)]

print(f"minimum seconds:    {min(bad_response_seconds):.2f}")
print(f"share below zero:   {sum(t < 0 for t in bad_response_seconds) / len(bad_response_seconds):.3f}")
print(f"p95 seconds:        {percentile(bad_response_seconds, 95):.2f}")
print(f"sample mean:        {fmean(bad_response_seconds):.2f}")

assert any(t < 0 for t in bad_response_seconds)

minimum seconds:    -15.70
share below zero:   0.095
p95 seconds:        17.85
sample mean:        8.09

About 9.5% of this sample is physically impossible. Clipping every negative draw to zero would create a pile of zero-duration requests and change the distribution, not establish a fit. A positive family, a justified truncated model, or sampling representative empirical durations are alternatives to evaluate against actual traces.

Even a plausible distribution for isolated durations misses another source of latency: waiting behind other requests. First distinguish time spent doing work from time spent in a queue.

Exponential service time

An exponential distribution models a duration with a constant completion rate, conditional on not having finished yet. It's the waiting-time distribution between arrivals in a constant-rate Poisson process. Using it for service time is a separate modeling assumption.

For a toy single-server service, let S be service time and μ = 20 completions per second when the server is busy. The mean duration is 1 / 20 = 0.05 seconds. This 50 ms service is a separate queueing example, not a new parameter for the six-second agent response-time model.

Its density and CDF are:

At 20 requests per second, P(S > 0.10) = e^{-20 \times 0.10} \approx 0.135. The distribution is memoryless: after a request has already taken 50 ms, its remaining service-time distribution is still the same as at the start. This makes it a useful queueing baseline. Real tool work can depend on elapsed time and request state.^{[2]Reference 2Machine Learning: A Probabilistic Perspective.https://probml.github.io/pml-book/book0.html[6]Reference 6Pattern Recognition and Machine Learning.https://www.microsoft.com/en-us/research/publication/pattern-recognition-machine-learning/}

Before running the sampler, predict a mean near .05 seconds and a share near .135 above .10 seconds. Random.expovariate(mu) takes the rate μ; NumPy's Generator.exponential takes the mean duration as scale = 1 / μ. Passing 20 as the scale instead of 0.05 makes the mean 400 times larger. Although exact zero has probability zero in the mathematical model, a finite-precision sampler can return it, so the support check allows zero.^{[1]Reference 1random: Generate pseudo-random numbershttps://docs.python.org/3/library/random.html}

from math import exp
from random import Random
from statistics import fmean

rate_per_second = 20.0
scale_seconds = 1 / rate_per_second
rng = Random(66)
service_seconds = [rng.expovariate(rate_per_second) for _ in range(5_000)]

print(f"target mean seconds: {scale_seconds:.3f}")
print(f"sample mean seconds: {fmean(service_seconds):.3f}")
print(f"P(S > 0.10) formula: {exp(-rate_per_second * 0.10):.3f}")
print(f"sample share > 0.10: {sum(t > 0.10 for t in service_seconds) / len(service_seconds):.3f}")
print(f"min seconds: {min(service_seconds):.6f}")

assert all(t >= 0 for t in service_seconds)

target mean seconds: 0.050
sample mean seconds: 0.049
P(S > 0.10) formula: 0.135
sample share > 0.10: 0.132
min seconds: 0.000025

The sample follows the chosen model. It doesn't show that real service times are exponential. Repeated simulation can estimate a model's tail probability more precisely without improving the model's fit to reality.

More draws reduce simulation noise

Monte Carlo simulation means using repeated random samples to approximate a quantity you care about. Here the quantity is the resolution rate, estimated from repeated simulated batches.

For N = 100 independent Bernoulli requests with p = 0.72, the batch-rate standard deviation is sqrt(0.72 × 0.28 / 100) ≈ 0.045. With N = 10,000, it's about 0.0045. This is the earlier standard-error formula: sqrt(p(1-p)/N). Multiplying the sample size by 100 divides this uncertainty by 10, not by 100.

Predict the three standard deviations before running the cell: they should be roughly .045, .014, and .0045. This pattern shows convergence under the assumed process. Production evidence requires fitting or validating the assumption against real traces.

from random import Random

def mean(xs: list[float]) -> float:
    return sum(xs) / len(xs)

def pop_std(xs: list[float]) -> float:
    center = mean(xs)
    return (sum((x - center) ** 2 for x in xs) / len(xs)) ** 0.5

rng = Random(55)
p_resolved = 0.72

for batch_size in [100, 1_000, 10_000]:
    repeated_rates = [
        rng.binomialvariate(batch_size, p_resolved) / batch_size
        for _ in range(2_000)
    ]
    print(
        f"N={batch_size:>5}: mean={mean(repeated_rates):.3f} "
        f"sd_of_rates={pop_std(repeated_rates):.4f}"
    )

N=  100: mean=0.720 sd_of_rates=0.0457
N= 1000: mean=0.720 sd_of_rates=0.0140
N=10000: mean=0.720 sd_of_rates=0.0045

Simulation answers a conditional question: "What happens if these assumptions are the generating process?" It doesn't validate p = 0.72, lambda = 2.2, or the latency parameters. Those must come from data, review design, and production measurement.

The draw budget also depends on the question. Suppose a costly failure has assumed probability q = 0.0001 per independent request. A 5,000-request simulation expects only Nq = 0.5 failures. The probability of seeing none is (1-q)^N, about 61%. A run with zero failures would be common even under this model.

from math import ceil, log, log1p

q = 0.0001
draws = 5_000
miss_probability = (1 - q) ** draws
draws_for_detection = ceil(log(0.05) / log1p(-q))

print(f"expected failures: {draws * q:.2f}")
print(f"probability of zero failures: {miss_probability:.3f}")
print(f"draws for 95% chance of at least one failure: {draws_for_detection:,}")

expected failures: 0.50
probability of zero failures: 0.607
draws for 95% chance of at least one failure: 29,956

That last budget makes P(no failures) ≤ 0.05 under the assumed q. It gives a chance to detect the event, not a precise estimate of its frequency. Even that budget expects only about three failures. Before interpreting a quiet simulation, check whether it had enough draws to expose the tail you care about.

To investigate a surprising run, record enough information to replay it rather than hoping the next random draw looks better.

Recorded seeds replay a stream

Random(seed) creates a deterministic pseudo-random stream. In the same environment, reusing the seed and draw order recreates a run. Record the code, Python and library versions, parameters, and number of draws too. Sampler algorithms can change across upgrades, so matching seeds alone may produce different results.^{[1]Reference 1random: Generate pseudo-random numbershttps://docs.python.org/3/library/random.html}

A seed doesn't make observations independent, and it doesn't repair a wrong model. Independence is a modeling assumption; separate streams help avoid accidentally replaying the same sequence across replicas.

For parallel replications, create separate Random objects and record each stream ID. The snippet derives two IDs from one root label. NumPy's Generator.spawn offers a designed child-stream mechanism for parallel work; distinct streams avoid replaying identical pseudo-observations as if they were new evidence.^{[7]Reference 7Parallel random number generationhttps://numpy.org/doc/stable/reference/random/parallel.html}

from random import Random

root = "incident-sim-2026"
replica_one = Random(f"{root}:replica-1")
replica_two = Random(f"{root}:replica-2")
counts_one = [replica_one.binomialvariate(20, 0.72) for _ in range(4)]
counts_two = [replica_two.binomialvariate(20, 0.72) for _ in range(4)]

replay_one = Random(f"{root}:replica-1")
replay_two = Random(f"{root}:replica-2")
replayed_one = [replay_one.binomialvariate(20, 0.72) for _ in range(4)]
replayed_two = [replay_two.binomialvariate(20, 0.72) for _ in range(4)]

print("replica 1:", counts_one)
print("replica 2:", counts_two)
print("replay matches:", counts_one == replayed_one and counts_two == replayed_two)

assert counts_one == replayed_one
assert counts_two == replayed_two

replica 1: [16, 14, 16, 12]
replica 2: [14, 16, 13, 17]
replay matches: True

Replay is a reporting contract. It also helps inspect a more difficult sampling problem: a Bayesian posterior that can't be sampled directly.

Sample a posterior when direct draws aren't available

Bayesian Inference from Evidence updated a Beta(2, 2) prior with eight passes and two failures to obtain an exact Beta(10, 4) posterior. That direct update worked because the prior and likelihood were conjugate. A richer agent model might include several interacting skill parameters or task-difficulty effects without a convenient closed-form posterior.

Markov chain Monte Carlo (MCMC) addresses that harder case by constructing dependent parameter states with the desired posterior as their stationary distribution. Reaching that distribution requires a valid transition rule and adequate exploration; a chain trapped in one region can give misleading estimates. Unlike the independent Monte Carlo draws above, each new state starts from the previous state.^{[2]Reference 2Machine Learning: A Probabilistic Perspective.https://probml.github.io/pml-book/book0.html}

Metropolis-Hastings proposes a candidate state, then either moves to it or stays at the current state. Start with a concrete comparison: the Beta(10, 4) posterior's unnormalized density is about 2.642 times as high at 0.6 as at 0.5. With a symmetric proposal, accept the move to 0.6 with probability 1; accept the reverse move with probability 1 / 2.642 ≈ 0.379.

Here symmetric means the proposal is equally likely in either direction, as with a normal step centered on the current value. If $\pi$ denotes the target density and $p'$ a proposed value, this rule is:^{[8]Reference 8Monte Carlo Sampling Methods Using Markov Chains and Their Applicationshttps://www2.stat.duke.edu/~sschmid/Courses/Stat376/Papers/Basic/Hastings1970.pdf}

The unknown normalization constant in $\pi$ cancels in the ratio. The earlier Beta(10, 4) result lets us check the mechanics even though this posterior can be sampled directly with betavariate. Its density is proportional to $p^9(1-p)^3$ for 0 < p < 1. An asymmetric proposal needs an additional reverse-to-forward proposal-density ratio; don't reuse the simplified rule unchanged.

The next cell checks both directions with a fixed comparison value of 0.4. Predict why that value accepts the forward move but rejects the reverse move. This is one deterministic acceptance check; a complete sampler must draw a fresh comparison value at every iteration.

def unnormalized_posterior(probability: float) -> float:
    return probability**9 * (1 - probability) ** 3

forward_ratio = unnormalized_posterior(0.6) / unnormalized_posterior(0.5)
reverse_ratio = 1 / forward_ratio
comparison = 0.4

print(f"forward acceptance: {min(1.0, forward_ratio):.3f}")
print(f"reverse acceptance: {min(1.0, reverse_ratio):.3f}")
print(f"reverse accepted at u=0.4: {comparison < reverse_ratio}")

forward acceptance: 1.000
reverse acceptance: 0.379
reverse accepted at u=0.4: False

A rejection must remain in the chain as another copy of the current state. Keeping only accepted moves would sample a different process. This complete teaching sampler uses log densities to avoid underflow, rejects proposals outside [0, 1], and compares four starting points with the known posterior mean 10 / 14:

from math import log, log1p
from random import Random
from statistics import fmean

def log_target(p):
    if not 0 < p < 1:
        return float("-inf")
    return 9 * log(p) + 3 * log1p(-p)

def sample_chain(seed, start, warmup=2_000, kept=10_000):
    rng = Random(seed)
    current = start
    draws = []
    accepted = 0
    for step in range(warmup + kept):
        proposal = current + rng.gauss(0, 0.15)
        log_ratio = log_target(proposal) - log_target(current)
        # 1-random() lies in (0, 1], so log is defined.
        if log(1 - rng.random()) < min(0.0, log_ratio):
            current = proposal
            if step >= warmup:
                accepted += 1
        if step >= warmup:
            draws.append(current)  # include repeats after rejection
    return draws, accepted / kept

for seed, start in enumerate([0.1, 0.4, 0.7, 0.95], start=101):
    draws, acceptance = sample_chain(seed, start)
    center = fmean(draws)
    variance_sum = sum((p - center) ** 2 for p in draws)
    lag_one = sum((a - center) * (b - center)
                  for a, b in zip(draws[:-1], draws[1:])) / variance_sum
    print(f"start={start:.2f} mean={center:.3f} "
          f"acceptance={acceptance:.3f} lag1={lag_one:.3f}")
    assert abs(center - 10 / 14) < 0.02  # a loose known-target smoke check
print(f"exact posterior mean: {10 / 14:.3f}")

start=0.10 mean=0.708 acceptance=0.644 lag1=0.715
start=0.40 mean=0.712 acceptance=0.637 lag1=0.720
start=0.70 mean=0.716 acceptance=0.637 lag1=0.720
start=0.95 mean=0.711 acceptance=0.632 lag1=0.726
exact posterior mean: 0.714

Early warmup or burn-in states can reflect initialization rather than the posterior. Discarding them addresses initialization; assess later exploration with chain diagnostics. The printed lag-one correlation measures similarity between adjacent draws; it can stay large even when the sample mean looks right.

Autocorrelation reduces the information in many ordinary chains. Effective sample size estimates precision for a particular quantity relative to independent draws; chain length alone doesn't determine it. Established tools report rank-normalized split R-hat, which compares variation within and between chains, alongside effective sample sizes. Passing the loose mean assertion above is a code smoke check against a known answer, not a convergence proof.^{[9]Reference 9Posterior Analysishttps://mc-stan.org/docs/reference-manual/analysis.html}

Gibbs sampling instead updates one variable from its conditional distribution given the current values of the others. Exact conditional draws can still produce a strongly correlated chain. Neither method establishes that the underlying statistical model represents the real system.

Posterior-sampling diagnostics qualify uncertainty about model parameters. They still say nothing about queues, which is where isolated latency samples go wrong.

Waiting time when the server fills up

A simulation that samples latency in isolation ignores how requests interact. Suppose three requests arrive at 0, 10, and 20 ms, and each needs 50 ms on one server. The first starts immediately. The second waits until 50 ms, and the third until 100 ms. Their total times are therefore 50, 90, and 130 ms, even though each does the same amount of work.

That is queueing delay: time waiting to start, separate from service time. A probabilistic queueing model adds random arrivals and service durations to this same mechanism.

Little's Law states that the average number of requests in a stable system $L$ equals the arrival rate $\lambda$ times the average time $W$ a request spends in the system:

For an M/M/1 baseline, assume constant-rate Poisson arrivals, independent exponential service times, one server, first-come-first-served processing, unlimited waiting space, and no abandonment. The two Ms denote the memoryless arrival and service models. At a steady arrival rate of 15 requests per second and a service rate of 20 per second, the server is busy on average 15 / 20 = 75% of the time.

Write utilization as $\rho = \lambda / \mu$. A steady-state distribution requires $\lambda < \mu$. Here $\lambda$ counts arrivals per second; the earlier tool-count parameter 2.2 counted calls per request. The same Greek letter doesn't imply the same units.^{[10]Reference 10Introduction to Queueing Theory: 6.263/16.37, Lectures 5 and 6https://web.mit.edu/modiano/www/6.263/lec5-6.pdf}

Average total time in the system $W$ and average time waiting in the queue $W_q$ are:

As utilization approaches 100%, both quantities rise sharply. Consider a service with a mean service time of 50 ms ($\mu = 20$ requests/second):

The service time stays at 50 ms in both rows, but moving from 75% to 95% utilization lifts average total time from 0.20 to 1.00 seconds. Queueing, not slower work, caused the jump.

These are steady-state averages, not predictions for a three-request burst or a p99 latency. At or above λ = μ, this model has no finite steady-state mean waiting time; plugging those values into the formula is invalid. Little's Law alone doesn't cause the sharp increase: that shape comes from the M/M/1 assumptions. Real batching, parallel workers, deadlines, and dependent retries need a different queueing model.

Two ways a latency simulator can lie are already on the table: the wrong positive-time family, and a missing queue. The last diagnostic uses the CDF you already have.

Score a candidate with the CDF

Before using a simulator for planning, compare it with representative held-out traces, including difficult requests and timeouts. A quantile-quantile (Q-Q) plot compares observed quantiles against candidate quantiles; disagreement near the high end exposes a tail mismatch. A timeout recorded at its deadline is a censored observation, not necessarily a completed response at exactly that duration.

The empirical CDF of a fully observed sample at a deadline is the fraction of durations at or below it. If 160 of 200 responses finish within 10 seconds, and all 200 eventually finish with observed durations, it equals 0.80 there. A candidate CDF of 0.823 has a gap of 0.023 at that deadline.

Now suppose those 200 completions came from 250 requests. Another 50 requests were still unfinished when measurement stopped at a 15-second deadline. These are right-censored completion times: each is known to exceed 15 seconds, but its eventual duration is unknown.^{[11]Reference 11scipy.stats.CensoredDatahttps://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.CensoredData.html} The deadline-completion fraction at 10 seconds uses all 250 requests, not just the 200 that finished:

finished_by_10 = 160
finished_between_10_and_15 = 40
unfinished_at_15 = 50
completed = finished_by_10 + finished_between_10_and_15
requests = completed + unfinished_at_15

print(f"completed-only fraction by 10s: {finished_by_10 / completed:.2f}")
print(f"all-request fraction by 10s: {finished_by_10 / requests:.2f}")
print(f"unfinished fraction at 15s: {unfinished_at_15 / requests:.2f}")

completed-only fraction by 10s: 0.80
all-request fraction by 10s: 0.64
unfinished fraction at 15s: 0.20

Dropping the unfinished requests turns 0.64 into 0.80 and makes latency look better. Because 20% remain unfinished at 15 seconds, this cohort's completion-time p95 lies beyond that deadline; these logs can't identify its actual value. If instead you measure time until the client receives either a response or a timeout, 15 seconds is an observed timeout duration. That's a different quantity from time to completion. Use a censoring-aware fit for completion times, and don't treat the deadline as a completed response time in an ordinary KS test. The shared measurement deadline here is fixed in advance; shorter deadlines assigned to harder requests would require accounting for that selection rule too.

The Kolmogorov-Smirnov (KS) distance $D$ takes the largest such vertical gap over all deadlines. In the notation below, sup means the largest gap, F_n is the empirical CDF, and F is the candidate CDF:

On the same sample, a larger D means a larger worst-case CDF discrepancy. It isn't a verdict on every operational risk: two models can differ mainly in a rare tail that contributes little to D. Keep your explicit timeout and costly-tail checks.

The ordinary one-sample KS p-value assumes independent observations from a continuous distribution whose parameters were specified in advance. Fitting parameters on the same observations needs a calibration procedure that accounts for fitting, such as a parametric bootstrap. SciPy's goodness_of_fit implements this approach.^{[12]Reference 12scipy.stats.goodness_of_fithttps://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.goodness_of_fit.html} This example only calculates distances, not p-values.

This check specifies both candidates first, using the same lognormal planning values as the rest of the chapter, then scores them against 200 positive draws. The normal candidate gets the lognormal's mean and standard deviation, so the comparison tests shape beyond those summaries. The generating lognormal model has lower expected KS distance; either candidate can score lower on a finite sample.

from math import erf, exp, log, sqrt
from random import Random

def normal_cdf(x: float, mu: float, sigma: float) -> float:
    z = (x - mu) / sigma
    return 0.5 * (1.0 + erf(z / sqrt(2.0)))

def lognormal_cdf(t: float, mu: float, sigma: float) -> float:
    if t <= 0:
        return 0.0
    return normal_cdf(log(t), mu, sigma)

def ks_distance(samples: list[float], cdf) -> float:
    if not samples:
        raise ValueError("KS distance needs at least one observation")
    ordered = sorted(samples)
    n = len(ordered)
    distance = 0.0
    for index, value in enumerate(ordered, start=1):
        fitted = cdf(value)
        distance = max(
            distance,
            abs(index / n - fitted),
            abs((index - 1) / n - fitted),
        )
    return distance

mu = log(6.0)
sigma = 0.55
rng = Random(88)
observed_traces = [rng.lognormvariate(mu, sigma) for _ in range(200)]

normal_mean = exp(mu + sigma**2 / 2)
normal_std = sqrt((exp(sigma**2) - 1) * exp(2 * mu + sigma**2))

d_norm = ks_distance(
    observed_traces,
    lambda x: normal_cdf(x, normal_mean, normal_std),
)
d_log = ks_distance(
    observed_traces,
    lambda t: lognormal_cdf(t, mu, sigma),
)

print(f"normal D:    {d_norm:.3f}")
print(f"lognormal D: {d_log:.3f}")
print(f"normal mean: {normal_mean:.2f}s  std: {normal_std:.2f}s")

normal D:    0.135
lognormal D: 0.056
normal mean: 6.98s  std: 4.15s

The lognormal gap is smaller in this seeded run. A finite sample still gives D > 0. This check asks whether our comparison distinguishes the known generating shape from an alternative. Validate against actual service traces before using the model for planning.

Build a launch-simulation report

Now put the four request fields together. This is deliberately an independent-field baseline: it draws each field without conditioning on the others. It has no arrivals, worker limits, or queue, so it describes synthetic requests, not system capacity. Keep that limitation attached to the report.

from math import exp, log
from random import Random
from statistics import fmean, median

def sample_poisson(rng: Random, lam: float) -> int:
    if not 0 <= lam <= 20:
        raise ValueError("this teaching sampler requires 0 <= lambda <= 20")
    threshold = exp(-lam)
    count = 0
    product = 1.0
    while True:
        product *= rng.random()
        if product <= threshold:
            return count
        count += 1

def percentile(values: list[float], q: float) -> float:
    ordered = sorted(values)
    position = (len(ordered) - 1) * q / 100
    lower = int(position)
    upper = min(lower + 1, len(ordered) - 1)
    return ordered[lower] + (ordered[upper] - ordered[lower]) * (position - lower)

def simulate_incident_agent(seed: int, n_requests: int) -> dict[str, float]:
    if n_requests <= 0:
        raise ValueError("n_requests must be positive")

    rng = Random(seed)
    resolved = [rng.binomialvariate(1, 0.72) for _ in range(n_requests)]

    intents = ["ci_failure", "rollback", "access_request", "log_triage"]
    intent_probs = [0.45, 0.20, 0.20, 0.15]
    routed = rng.choices(intents, weights=intent_probs, k=n_requests)

    tool_calls = [sample_poisson(rng, 2.2) for _ in range(n_requests)]
    response_seconds = [
        rng.lognormvariate(log(6.0), 0.55) for _ in range(n_requests)
    ]

    return {
        "resolution_rate": fmean(resolved),
        "rollback_share": routed.count("rollback") / n_requests,
        "mean_tool_calls": fmean(tool_calls),
        "share_over_5_tools": sum(count > 5 for count in tool_calls) / n_requests,
        "median_seconds": median(response_seconds),
        "p95_seconds": percentile(response_seconds, 95),
        "nonpositive_times": float(sum(t <= 0 for t in response_seconds)),
    }

report = simulate_incident_agent(seed=7, n_requests=5_000)
for name, value in report.items():
    print(f"{name:>20}: {value:.3f}")

assert report["nonpositive_times"] == 0
assert report["p95_seconds"] > report["median_seconds"]

resolution_rate: 0.718
      rollback_share: 0.202
     mean_tool_calls: 2.206
  share_over_5_tools: 0.024
      median_seconds: 6.070
         p95_seconds: 15.688
   nonpositive_times: 0.000

Use that output as a simulator report. The rate, label mix, count tail, and latency tail are only as credible as their fitted parameters and assumption checks. Validate those assumptions against production data before making operational claims.

This fixture samples each field separately from one generator. It therefore treats resolution, route, tool count, and response time as independent once their parameters are fixed.

A shared random-number generator (RNG) makes the run replayable; it doesn't create or prove independence. Real rollback requests may use more tools and take longer. Preserve that relationship with conditional sampling, such as drawing route first, then drawing tool count and latency from models for that route.

To see what separate summaries miss, compare two four-request datasets. They have identical tool counts and identical durations, but pair those values differently. Predict how often a request is both expensive and slow:

tools = [1, 1, 6, 6]
aligned_seconds = [2, 2, 20, 20]
reversed_seconds = [20, 20, 2, 2]

def joint_share(counts, seconds):
    return sum(c > 5 and t > 15 for c, t in zip(counts, seconds, strict=True)) / len(counts)

assert sorted(aligned_seconds) == sorted(reversed_seconds)
print(f"aligned expensive-and-slow:  {joint_share(tools, aligned_seconds):.2f}")
print(f"reversed expensive-and-slow: {joint_share(tools, reversed_seconds):.2f}")

aligned expensive-and-slow:  0.50
reversed expensive-and-slow: 0.00

Every one-field histogram, mean, and percentile is identical between these datasets. Yet the joint threshold shares differ. Matching each marginal distribution separately doesn't validate the relationships among fields, which matter for timeouts and simultaneous resource demand.

Those checks apply to the simulator's fields. The same categorical math appears inside the model that writes the assistant's answer.

Next tokens: categorical draws, temperature scaling, and the Gumbel-Max trick

The incident-assistant labels form a categorical distribution because one label is selected from a finite list. An LLM uses the same family over a much larger vocabulary. At each step it produces a logit, an unnormalized score, for each candidate token. Softmax exponentiates those scores and divides by their sum, producing probabilities. A decoder may sample one token from that distribution.^{[13]Reference 13Deep Learning.https://www.deeplearningbook.org/}

Take a tiny continuation for the prompt The sky is. Four words are enough to see the mechanics. blue has the largest probability, while a sample can still choose clear; predict that distinction before running the cell.

from math import exp
from random import Random

tokens = ["blue", "clear", "cloudy", "green"]
logits = [2.2, 1.5, 0.2, -1.5]

def stable_softmax(values: list[float]) -> list[float]:
    maximum = max(values)
    shifted = [value - maximum for value in values]
    weights = [exp(value) for value in shifted]
    total = sum(weights)
    return [weight / total for weight in weights]

probabilities = stable_softmax(logits)
rng = Random(5)
chosen = rng.choices(tokens, weights=probabilities, k=1)[0]

for token, probability in zip(tokens, probabilities):
    print(f"{token:>8}: {probability:.3f}")
print("sampled token:", chosen)
print("probabilities sum to one:", abs(sum(probabilities) - 1.0) < 1e-12)

blue: 0.604
   clear: 0.300
  cloudy: 0.082
   green: 0.015
sampled token: clear
probabilities sum to one: True

blue is the mode, the most probable outcome, while this draw returned clear. Greedy decoding would choose blue; sampling can select any word with nonzero probability. Subtracting the largest logit before exponentiating prevents overflow for these finite inputs without changing the probability ratios. Displayed rounded probabilities may add to 1.001; the unrounded values still sum to one.

Reshaping entropy with temperature scaling

Before feeding logits into softmax, generation runtimes introduce a positive temperature parameter $T > 0$ that scales logits: $\tilde{z}_i = z_i / T$.

Dividing logits by temperature alters the distribution's entropy without changing which token has the highest score:

• Low temperature ($T < 1.0$, such as $T = 0.3$): Logit gaps stretch wider. Dividing by 0.3 multiplies differences by $3.33\times$. Here blue climbs from 60.4% to 91.1%, while clear still has about 8.8% probability. This remains sampling. As $T \to 0$, probability concentrates on the highest logit; with a unique maximum, the limit is greedy selection. Tied maxima share the limiting mass unless a separate tie-breaking rule selects one.
• Baseline temperature ($T = 1.0$): Softmax evaluates unscaled logits. That doesn't establish calibration; a model can assign confident probabilities to incorrect continuations.
• High temperature ($T > 1.0$, such as $T = 2.0$): Logit gaps shrink. Here blue drops to 44.8% and green rises to 7.1%. For a fixed finite vocabulary and finite logits, the distribution approaches uniform ($1/4 = 25\%$ here) as $T \to \infty$. At T=2, blue remains the most likely word. More varied continuations don't by themselves establish better answers or a particular hallucination rate.

GPU-parallel sampling with the Gumbel-Max trick

On GPU hardware serving thousands of concurrent requests, how does code draw a categorical sample from 128,000 candidate tokens?

One option is to build a cumulative probability array and locate a uniform draw in it. Prefix sums can run in parallel too. Another option is the Gumbel-Max trick, which expresses categorical sampling as independent noise plus an argmax. Runtimes choose among kernels based on hardware, vocabulary size, filtering, and reproducibility requirements.

Draw an independent standard Gumbel noise value for each vocabulary token from uniform random numbers:

Add that Gumbel noise directly to the temperature-scaled logits, then take the parallel argmax:

With independent standard Gumbel draws, token $k$ wins with its softmax probability in the mathematical model:Gumbel-Max derivation.

The transformation can run in parallel, and it also needs random-number generation, logarithms, and a reduction. Measure speed and memory-bandwidth use on the target workload. As checked on October 10, 2026, vLLM's PyTorch-native sampler uses an equivalent exponential-noise race, while eligible filtered CUDA paths can use FlashInfer rejection-sampling kernels. The mathematical distribution can agree while seeded outputs differ between implementations.

The hard argmax here selects a discrete token. It doesn't provide useful gradients through that chosen index. The Gumbel-Softmax training method replaces it with a smooth relaxation; this example uses the categorical sampler, not that relaxation.Gumbel-Softmax paper

View larger

In the bottom panel, word clear starts with a lower base logit (1.50) than blue (2.20). This draw gives clear a +2.25 Gumbel perturbation and blue a +0.77 perturbation, moving their scores to 3.75 and 2.97. clear therefore wins the argmax. Across thousands of draws, each candidate's empirical win frequency approaches its softmax probability.

The next snippet compares empirical frequencies with the softmax probabilities across 10,000 draws at two temperatures. It illustrates the equivalence rather than proving it. Redraw endpoints before evaluating the nested logarithms: 1 - random() can equal one when random() returns zero.

from math import exp, log
from random import Random

def stable_softmax(logits: list[float], temperature: float = 1.0) -> list[float]:
    scaled = [z / temperature for z in logits]
    max_z = max(scaled)
    weights = [exp(z - max_z) for z in scaled]
    total = sum(weights)
    return [w / total for w in weights]

def gumbel_sample(rng: Random, logits: list[float], temperature: float = 1.0) -> int:
    # Standard Gumbel noise: g = -log(-log(u)) where u ~ Uniform(0, 1)
    perturbed = []
    for z in logits:
        u = 1.0 - rng.random()
        while u == 1.0:
            u = 1.0 - rng.random()
        # Now 0 < u < 1, so both logarithms are defined.
        g = -log(-log(u))
        perturbed.append(z / temperature + g)
    return max(range(len(logits)), key=lambda i: perturbed[i])

tokens = ["blue", "clear", "cloudy", "green"]
logits = [2.2, 1.5, 0.2, -1.5]

rng = Random(42)
n_draws = 10_000

for temp in [1.0, 0.3]:
    expected_probs = stable_softmax(logits, temperature=temp)
    samples = [gumbel_sample(rng, logits, temperature=temp) for _ in range(n_draws)]
    observed_shares = [samples.count(i) / n_draws for i in range(len(tokens))]

    print(f"Temperature T={temp:.1f}:")
    for token, p_exp, p_obs in zip(tokens, expected_probs, observed_shares):
        print(f"  {token:>8}: softmax={p_exp:.3f}  gumbel_share={p_obs:.3f}")

Temperature T=1.0:
      blue: softmax=0.604  gumbel_share=0.610
     clear: softmax=0.300  gumbel_share=0.292
    cloudy: softmax=0.082  gumbel_share=0.084
     green: softmax=0.015  gumbel_share=0.014
Temperature T=0.3:
      blue: softmax=0.911  gumbel_share=0.902
     clear: softmax=0.088  gumbel_share=0.097
    cloudy: softmax=0.001  gumbel_share=0.001
     green: softmax=0.000  gumbel_share=0.000

Later chapters teach decoding controls in depth. A sampling decoder doesn't choose from raw logits directly. It draws from a probability distribution derived from those logits. A greedy decoder is different: it selects the highest-probability token instead of drawing a sample. For open-ended text generation, nucleus sampling chooses a dynamic high-probability token set and renormalizes it before sampling, rather than sampling from an unreliable long tail.^{[14]Reference 14The Curious Case of Neural Text Degeneration.https://arxiv.org/abs/1904.09751}

Keep the simulation honest

Before believing output from a random simulator, write down its contract:

The seed also belongs in the report. It makes one simulated run reproducible, but it doesn't make an incorrect model correct.

Practice: challenge a candidate launch model

A teammate proposes this incident-assistant simulator:

Write a review with four parts:

1. Name one parameter that needs evidence from logs or reviewed traffic.
2. Identify the impossible or missing state.
3. State one count diagnostic for the tool-call assumption.
4. State which latency summary should appear in the report.

Before trusting a simulator, build a launch-model worksheet from one seeded run. Record each assumed parameter, check labels and support, compare count variance with the mean, and report median, p95, and costly-tail rates. Use the worksheet to reject or revise assumptions before simulated capacity or cost numbers reach a design review.