, ,

Probability and Distributions for Infra Telemetry (Infra to Data Science Series, Part 11)

Fit named distributions to your own latency and arrival data in Python, and see exactly where a clean fitted curve understates the tail that breaches your SLO.

Infra to Data Science Series · Part 11 of 26

A curve I had fitted to a month of latency told me that a request slower than 300 milliseconds should show up about nine times in a hundred thousand. My logs said it happened three hundred and twenty eight times in a hundred thousand, thirty six times more often than the model I trusted. Both numbers came from the same data. One came from assuming a clean shape, the other from counting. Naming the distribution your telemetry actually follows, and knowing exactly where that name stops being true, is the difference between a capacity plan that holds and one that pages you at 3 a.m.

Key takeaways: A log-normal fit to 500,000 request latencies landed at shape s 0.604 and scale 31.3 ms, an implied median of 31.3 against an actual 30.6. A Kolmogorov-Smirnov test on all 500,000 points returned p 7.22e-244 and rejected that fit, while the same test on a 1,000 point subsample returned p 0.377 and accepted it, so the test tracked sample size more than fit quality. The fitted curve put the chance of a latency past 200 ms at 0.107 percent while the data showed 0.764 percent, a sevenfold understatement, and at 300 ms the gap ran thirty six fold. Request arrivals matched a Poisson process, mean 42.08 per second against variance 41.8, near the equality Poisson demands.
Who this is for: An infrastructure engineer or SRE who has described a metric with the percentiles from Part 10 and now wants to name the shape underneath it. Terms on first use: a probability distribution is a rule that assigns likelihood to each possible value; log-normal means the logarithm of the value is bell shaped, common for latency and payload size; a Poisson process is a stream of independent events arriving at a steady average rate; a survival function, sf, is one minus the cumulative distribution, the chance of exceeding a value; goodness of fit is how closely a named distribution matches your data.

Name the shape before you trust a summary

Where the project stands: last part you described the latency export honestly, a median, tail percentiles and a bootstrap interval. This part names the shape those numbers came from, because once you can say latency is log-normal you get every percentile, every exceedance probability and a way to simulate more traffic from two parameters, instead of carrying a lookup table of quantiles. Latency, payload size and time between failures nearly always lean right, and log-normal is the first shape to try on a right leaning positive metric. Fit it, and fix the location at zero, because a request cannot take negative time and a free location parameter will happily slide the whole curve to invent one.

# tested on Python 3.10.12, numpy 2.2.6, scipy 1.15.3
import numpy as np
from scipy import stats
# lat: 500,000 request latencies in ms, the same monitoring export from Part 10
s, loc, scale = stats.lognorm.fit(lat, floc=0)   # fix location at 0, latency has a floor
print('median', round(np.median(lat), 1), 'mean', round(lat.mean(), 1))
print('lognorm s', round(s, 3), 'scale', round(scale, 1))
print('implied median', round(scale, 1))          # for a log-normal, median equals scale
median 30.6 mean 38.6
lognorm s 0.604 scale 31.3
implied median 31.3

Two parameters now stand in for the whole export. Shape s of 0.604 is the spread of the logarithm, and scale 31.3 ms is the median the fit implies, a hair above the actual 30.6 ms. That closeness is encouraging and it is also a trap, because a fit that nails the median can still miss the part of the curve that pages you. One production gotcha lives right here in the location parameter. Leave loc free and the fit sometimes settles on a small negative or positive shift that reads fine in a plot yet quietly breaks the survival function near zero and makes the tail probabilities wrong. Pass floc=0 whenever the metric has a physical floor, which for latency, size and counts is always. Distribution mechanics get built from scratch in the Data Science Series part on probability and distributions, which this series leans on rather than repeats.

A named shape earns its keep in three ways a quantile table cannot. You can extrapolate past the largest value you have on record, so a month of data can speak to a once a quarter event, with the honesty caveat this part keeps returning to. You can simulate, drawing a million synthetic requests to load test a queue or a retry budget before it ever meets real traffic. And you can compress, carrying a service latency profile as two numbers in a config file rather than a histogram in a database. Each of those is a lever the raw percentiles from Part 10 do not hand you, which is the whole reason to fit a distribution at all rather than stop at the summary.

Distributions that fit infrastructure signals

You do not need a zoo of distributions. Five shapes cover almost everything a running system emits, and each one carries a physical story that tells you when it applies. Time between events is exponential. Events counted in a window are Poisson. A right skewed positive quantity such as latency or file size is log-normal. Time to a hardware failure is Weibull, whose shape parameter says whether parts are dying young or wearing out. Rare enormous spikes belong to a heavy tail such as Pareto, where the variance can be so large it barely exists. Keep the table below as your reference artifact, one row per signal, the distribution it follows, the parameter that carries the meaning, and how a normal curve fails on it.

flowchart TD
  A[Operational signal] --> B{What does it measure}
  B -->|Time between events| C[Exponential]
  B -->|Events per window| D[Poisson]
  B -->|Latency or size, right skew| E[Log-normal]
  B -->|Time to failure| F[Weibull]
  C --> G[Confirm rate is steady]
  D --> H[Confirm mean equals variance]
  E --> I[Confirm the tail, not just the median]
  F --> I
Pick a candidate shape from what the signal physically measures, then confirm it with one check that matches that shape, never by eye alone.
Infrastructure signalNamed distributionParameter that mattersHow a normal curve fails
Time between requestsExponentialrate, events per secondinvents negative gaps
Requests per secondPoissonlambda, and mean equals variancemisses burst structure
Latency, payload sizeLog-normals, spread of the logmean minus 2 sd goes negative
Time between failuresWeibullshape k, wear in or wear outassumes a constant hazard
Rare huge spikesPareto, heavy tailtail indexvariance may not be finite

Two of those shapes deserve a second look because they carry a decision, not just a label. A Weibull shape parameter k below 1 means the failure rate falls with age, the infant mortality pattern where a fresh disk or a just deployed pod is likeliest to die early, while k above 1 means wear out, with risk climbing the longer a part stays in service, and those two readings call for opposite maintenance policies. A heavy tail such as Pareto is stranger still: once its tail index drops under 2 the variance is infinite in theory and useless in practice, so a standard deviation computed on spike sizes becomes a number that never settles no matter how much data you feed it. Naming that shape is how you stop quoting a summary statistic the data can never support, and it is a mistake I have watched cost a full afternoon of arguing about a chart.

A goodness of fit test that lies on big samples

Having fitted a shape, the textbook next move is a goodness of fit test, and the Kolmogorov-Smirnov test is the one every tutorial reaches for. It measures the largest gap between your data and the fitted curve and hands back a p-value. Run it on the full export and on a small subsample of the very same data.

from scipy import stats
D, p = stats.kstest(lat, 'lognorm', args=(s, 0, scale))
print('full N', len(lat), 'KS D', round(D, 5), 'p', format(p, '.2e'))
sub = np.random.default_rng(11).choice(lat, 1000, replace=False)
ss, _, sc = stats.lognorm.fit(sub, floc=0)
D2, p2 = stats.kstest(sub, 'lognorm', args=(ss, 0, sc))
print('subsample 1000 KS D', round(D2, 5), 'p', round(p2, 3))
full N 500000 KS D 0.02367 p 7.22e-244
subsample 1000 KS D 0.02866 p 0.377

Look at the two D values first. They are almost identical, 0.02367 against 0.02866, so both samples sit about the same distance from the fitted curve. Yet the p-value swings from 7.22e-244, a flat rejection, to 0.377, a comfortable pass. Only the sample size changed. With half a million points a test can resolve a mismatch far too small to matter, so it rejects a fit that is, for every practical purpose, excellent. That contradicts the way significance testing is usually taught, where a tiny p-value means a real problem. On large operational samples a KS p-value measures how many rows you have, not whether the shape is right. A second, quieter fault sits underneath: because you estimated s and scale from the same data you then tested, the p-value is not even valid without a Lilliefors or Monte Carlo correction. Judge the fit by a quantile-quantile plot and by the exceedance error you care about, not by this p-value.

A quantile-quantile plot is the honest replacement, and it needs no p-value at all. Sort your data, sort an equal count of draws from the fitted distribution, and plot one against the other; a good fit lands on the diagonal and every departure is visible in the units you already reason in. On this latency the points hug the line through the body and bow away above the 95th percentile, which is the same tail miss the exceedance numbers are about to make concrete, seen as a shape instead of a p-value. That picture survives a large sample where the KS number does not, because your eye is measuring distance in milliseconds rather than counting rows, and a departure that matters looks large while one that does not stays small.

Verdict: To judge a distribution fit on real telemetry, my pick is a quantile-quantile plot paired with the exceedance error at the thresholds you alert on, since both speak in the units of the decision. What to avoid is gating the fit on a KS p-value when the sample runs past a few thousand points, because it will reject a fit you should keep and it is statistically invalid once parameters come from the same data.

Arrivals, counts and the Poisson process

Latency is one half of the story, request arrivals are the other, and they follow a different shape you can name. If requests arrive independently at a steady average rate, the gaps between them are exponential and the count in any fixed window is Poisson. That pairing is worth knowing because a Poisson count has one signature you can check in a line: its mean equals its variance. Build per second counts from exponential inter arrival gaps and read the two numbers.

rng = np.random.default_rng(11)
rate = 42.0                                    # requests per second on one service
inter = rng.exponential(1.0 / rate, 300_000)   # inter arrival gaps in seconds
t = np.cumsum(inter)
counts = np.bincount(np.floor(t).astype(int))[:-1]   # arrivals per one second bin
print('counts per second: mean', round(counts.mean(), 2), 'var', round(counts.var(), 2))
counts per second: mean 42.08 var 41.8

Mean 42.08 and variance 41.8 sit almost on top of each other, which is the Poisson signature and a genuine capacity signal. When variance runs far above the mean your traffic is overdispersed, bursty rather than smooth, and a plain Poisson queue model will undersize you. Now try to fit the count the obvious way, by asking the Poisson object to fit itself, and watch it fail.

lam = stats.poisson.fit(counts)                # the reflex call
AttributeError: 'poisson_gen' object has no attribute 'fit'

Discrete distributions in scipy do not carry a fit method, only the continuous ones do, so the call dies with that AttributeError. There is no need for it anyway. A Poisson rate has a one line maximum likelihood estimate, the sample mean, so estimate it directly.

lam = counts.mean()                            # MLE for a Poisson rate is the sample mean
print('lambda', round(lam, 2))                 # lambda 42.08

Overdispersion is not a curiosity, it is a capacity number. Queueing results that assume Poisson arrivals, the kind that turn a target utilisation into a p99 wait, turn optimistic the moment real traffic clumps, because bursts build queues a smooth arrival rate never would. When I checked one gateway its per minute variance ran about 3 times its mean, so I stopped sizing headroom from the average rate and started sizing from the 99th percentile of the per minute count, which was roughly 1.7 times that mean. That single change moved a planned node count from 6 to 9, and it was the whole distance between a quiet peak hour and a paged one. Check mean against variance before you trust any arrival based sizing, because the average rate is the most flattering number in the room.

Tail risk your fitted curve understates

Here is what a fitted distribution is actually for: answering a question the raw data answers slowly, such as how often a request will breach an SLO. The survival function gives it in one call. Ask the fitted log-normal for the chance of exceeding 100, 200 and 300 ms, and put its answer next to the empirical count from the data.

for slo in [100, 200, 300]:
    emp = (lat > slo).mean() * 100
    fit = stats.lognorm.sf(slo, s, 0, scale) * 100    # sf is 1 minus cdf, the tail
    print('P(lat over', slo, 'ms) empirical', round(emp, 3), 'fitted', round(fit, 3))
P(lat over 100 ms) empirical 3.2 fitted 2.729
P(lat over 200 ms) empirical 0.764 fitted 0.107
P(lat over 300 ms) empirical 0.328 fitted 0.009

At 100 ms the fit is close, 2.729 against 3.2 percent, because that threshold still sits in the bulk the curve was fitted to. Push into the tail and the fit collapses. At 200 ms the data breaches seven times as often as the model predicts, and at 300 ms thirty six times as often. Nothing is wrong with the log-normal, it fit the middle beautifully, which is exactly the danger. Real latency is a mixture, a fast bulk from warm caches plus a slow component from garbage collection pauses, cold starts and retries, and a single clean curve smooths that second hump away. For capacity and SLO work, read the deep tail off the empirical data or fit a mixture, and never size a timeout from a single fitted distribution alone. Latency distributions are precisely what tracing tools plot when the AI Engineering Series argues for observability and tracing, where the shape of the tail, not a single average, tells you which calls are hurting.

Empirical tail versus fitted log-normalChance a latency exceeds each threshold, percent0123100 ms3.22.73200 ms0.7640.1077x300 ms0.3280.00936xempiricalfitted
The two bars match in the bulk at 100 ms and split apart in the tail, where the fitted curve falls to a ninth then a thirty sixth of the real breach rate. That gap is the slow component a single distribution cannot see.
ThresholdEmpirical breachFitted log-normalUnderstatement
100 ms3.2 percent2.729 percent1.2 times
200 ms0.764 percent0.107 percent7.1 times
300 ms0.328 percent0.009 percent36 times
War story: I sized a downstream timeout from a log-normal I had fitted to a month of latency. The curve said a breach past 300 ms would hit about 1 request in 11,000, so a 300 ms timeout felt generous and I signed off. Production breached at roughly 1 in 300, thirty times heavier, because the fit had erased a garbage collection hump the raw data plainly showed. We took a wave of timeout errors during the next deploy, and I spent a Saturday moving every SLO calculation from the fitted curve back to the empirical tail. The distribution was not wrong, my decision to trust it past the range it described was.

Fit a distribution to one signal this week

Do one thing before the next part. Pull a real latency series you own, fit a log-normal with floc=0, and compute the chance of exceeding your SLO two ways, from the fitted survival function and by counting the data. I expect the two to agree in the bulk and split in the tail, and that split is the single most useful thing this part can hand you. Then bin the arrivals per second and compare the mean to the variance, because if they diverge your traffic is burstier than a Poisson queue assumes. If you would rather practice on shared data first, the Numenta Anomaly Benchmark carries real AWS server metrics you can fit the same way. With the shapes of your telemetry named and their limits known, the next part turns to machine learning fundamentals framed for infra engineers, where these distributions become the features a model learns from.

Infra to Data Science Series · Part 11 of 26
« Previous: Part 10  |  Guide  |  Next: Part 12 »

References

About The Author


Discover more from Journal of Intelligent Infrastructure

Subscribe to get the latest posts sent to your email.

Leave a Reply

Your email address will not be published. Required fields are marked *

Architect’s Toolkit

About the Author

Dr. Pranay Jha is a Cloud and AI Consultant with 18+ years of experience in hybrid cloud, virtualization, and enterprise infrastructure transformation. He specializes in VMware technologies, multi-cloud strategy, and Generative AI solutions. He holds a PhD in Computer Applications with research focused on Cloud and AI, has published multiple research papers, and has been a VMware vExpert since 2016 and a VMUG Community Leader.

Discover more from Journal of Intelligent Infrastructure

Subscribe now to keep reading and get access to the full archive.

Continue reading