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. We'll start with those allowed values, then check the frequencies, tails, and dependencies that a simulator produces.

The numbers below are planning assumptions, not measured production results. 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:

$$F_T(t) = P(T \leq t)$$

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. Here you pick a named family, sample from it, and check whether traces match the family.

Inverse CDF transform: how computers sample continuous variables

How does a random-number generator turn deterministic code into smooth, continuous response times? A computer's core RNG primitive only knows how to produce uniform numbers between zero and one: $U \sim \text{Uniform}(0, 1)$.

The mathematical bridge from uniform numbers to any continuous distribution is the inverse CDF transform (also known as the probability integral transform). Because a continuous CDF $F(t) = P(T \leq t)$ is strictly increasing from 0 to 1, its inverse $F^{-1}(u)$ exists. If you draw $U \sim \text{Uniform}(0, 1)$, the transformed variable $T = F^{-1}(U)$ follows the target distribution $F$ exactly:

$$P(F^{-1}(U) \leq t) = P(U \leq F(t)) = F(t)$$

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

$$F(s) = 1 - e^{-\mu s} = u \implies e^{-\mu s} = 1 - u \implies s = -\frac{\ln(1 - u)}{\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)$.

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 very close to the theoretical 0.05-second exponential mean, the 6.0-second lognormal median, and the 14.83-second p95 threshold.

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:

$$E[R] = 0 \times 0.28 + 1 \times 0.72 = 0.72$$

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

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:

$$P(X=72) = \binom{100}{72}(0.72)^{72}(0.28)^{28} \approx 0.089$$

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.

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.

What is the difference between p = 0.72 and a sampled batch rate of 0.68?

Answer

p = 0.72 is the assumed long-run resolution probability in the simulator. 0.68 is one finite batch generated under that assumption. A different batch can land above or below 0.72.

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.

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:

$$P(C=k) = \frac{\lambda^k e^{-\lambda}}{k!}$$

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

$$P(C=2) = \frac{2.2^2 e^{-2.2}}{2} = \frac{4.84 \times e^{-2.2}}{2} \approx 0.268$$

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

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.

The mixed sample has both excess variance and a larger expensive-tail share. 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.

Why isn't a Poisson distribution automatically correct for agent tool calls just because tool calls are counts?

Answer

Poisson is a useful first count model, but it assumes the variance matches the mean. Agent retries and difficult request clusters can create much larger variance and a heavier expensive tail. Check the observed count shape before relying on the baseline.

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}

$$\text{median}(T) = e^\mu \qquad E[T] = e^{\mu + \sigma^2/2}$$

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:

$$z = \frac{\log 10 - \log 6}{0.55} = \frac{\log(10/6)}{0.55} \approx 0.93$$

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.

The median describes a typical request. p95 tells you about a noticeably slow slice. Reporting only the mean would hide the operational question customers feel: how long do the slow responses take?

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:

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:

$$f_S(s) = \mu e^{-\mu s}, \qquad F_S(s) = P(S \leq s) = 1 - e^{-\mu s}, \qquad s \geq 0$$

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. That makes it a useful queueing baseline, not a claim that real tool work forgets its elapsed time.^{[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}

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. That pattern is convergence under the assumption, not evidence that the assumption came from production.

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.

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. The seed alone isn't a promise that every distribution sampler produces identical results across upgrades.^{[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}

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 a sequence of dependent parameter states whose long-run distribution matches the desired posterior. 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}

$$a(p,p')=\min\left(1,\frac{\pi(p')}{\pi(p)}\right)$$

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.

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:

Early warmup or burn-in states can reflect initialization rather than the posterior. Discarding them doesn't guarantee that later states explore well. 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:

$$L = \lambda W$$

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:

$$W = \frac{1}{\mu - \lambda} \qquad W_q = \frac{\lambda}{\mu(\mu - \lambda)}$$

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 at a deadline is simply the fraction of recorded durations at or below it. If 160 of 200 completed responses finish within 10 seconds, it equals 0.80 there. A candidate CDF of 0.823 has a gap of 0.023 at that deadline.

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:

$$D = \sup_t \left| F_n(t) - F(t) \right|$$

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.^{[11]Reference 11scipy.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. We expect the generating lognormal model to fit better on average, though a finite sample doesn't guarantee that ordering.

The lognormal gap is smaller in this seeded run. A finite sample still doesn't give D = 0. This is a check that our comparison can distinguish the known generating shape from an alternative, not validation against actual service traces.

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.

Use that output as a simulator report, not a production claim. The rate, label mix, count tail, and latency tail are only as credible as their fitted parameters and assumption checks.

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:

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.^{[12]Reference 12Deep Learning.https://www.deeplearningbook.org/}

Take a tiny continuation for the prompt Deploy is. Four tokens are enough to see the mechanics. approved has the largest probability, but a sample can still choose pending; predict that distinction before running the cell.

approved is the mode, the most probable outcome, but this draw returned pending. Greedy decoding would choose approved; sampling can select any token 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$. Exponentiating those stretched gaps concentrates almost all mass on the top token (approved climbs from 60.4% to 91.1%), squashing competitor probabilities. As $T \to 0$, categorical sampling converges to deterministic greedy argmax.
• Baseline temperature ($T = 1.0$): Softmax evaluates raw logits directly, preserving the model's calibrated uncertainty.
• High temperature ($T > 1.0$, such as $T = 2.0$): Logit gaps shrink toward zero. The probabilities flatten toward a uniform distribution ($1/4 = 25\%$), making generation exploratory and diverse, but increasing the risk of low-probability hallucinated tokens like coupon.

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?

Computing a cumulative probability array (CDF) and performing sequential binary search requires expensive prefix sums and branches that cause GPU thread divergence. Instead, production runtimes use the Gumbel-Max trick.

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

$$g_i = -\ln(-\ln(u_i)), \qquad u_i \sim \text{Uniform}(0, 1)$$

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

$$i^* = \operatorname{argmax}_i \left( \frac{z_i}{T} + g_i \right)$$

A celebrated result in extreme value theory proves that the probability of token $k$ achieving the highest perturbed score equals its softmax probability exactly:

$$P\left(\operatorname{argmax}_i \left( \frac{z_i}{T} + g_i \right) = k\right) = \frac{\exp(z_k / T)}{\sum_j \exp(z_j / T)}$$

This reparameterization turns sampling into a parallel vector addition and an argmax reduction, running at native GPU memory bandwidth.

In the bottom panel, notice how token pending had a lower base logit (1.50) than approved (2.20). But because pending drew a lucky Gumbel perturbation (+2.25 versus +0.77), its perturbed score reached 3.75, winning the argmax. Across thousands of draws, each token wins with an empirical frequency that mirrors its softmax probability.

The next snippet verifies that equivalence across 10,000 Gumbel-Max draws at both baseline and sharp temperatures:

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.^{[13]Reference 13The 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.

What is a strong answer to the practice review?

Answer

The resolution probability needs evidence from a representative set of requests. The categorical model is missing log_triage, so it can't represent every stated route. Tool-call counts should be checked for overdispersion by comparing variance with mean and reporting the share above a costly threshold. A normal response-time model can create negative seconds when spread is large; use a positive baseline and report median plus p95.

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.