⚙️ Setup & data¶
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
AMBER="#d97706"; TEAL="#0d9488"; INK="#1a2138"; GRID="#e6e9f2"; PINK="#db2777"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
"axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
"axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
try: df = pd.read_csv("../../data/it_server_traffic.csv")
except FileNotFoundError: df = pd.read_csv(BASE+"it_server_traffic.csv")
print("loaded:", df.shape)
df.head()
loaded: (1000, 4)
| timestamp | server_id | api_endpoint | requests_per_second | |
|---|---|---|---|---|
| 0 | 2026-06-01 12:00:00 | SRV_NODE_A | /v1/checkout | 11 |
| 1 | 2026-06-01 12:00:01 | SRV_NODE_A | /v1/users | 8 |
| 2 | 2026-06-01 12:00:02 | SRV_NODE_B | /v1/users | 8 |
| 3 | 2026-06-01 12:00:03 | SRV_NODE_A | /v1/products | 11 |
| 4 | 2026-06-01 12:00:04 | SRV_NODE_B | /v1/users | 10 |
x = df["requests_per_second"]
lam = x.mean()
print(f"observations : {len(x):,}")
print(f"mean (lambda hat) : {x.mean():.3f}")
print(f"variance : {x.var():.3f}")
print(f"-> mean ~ variance, the Poisson signature ({x.var()/x.mean():.3f} ratio)")
observations : 1,000 mean (lambda hat) : 12.007 variance : 11.691 -> mean ~ variance, the Poisson signature (0.974 ratio)
fig,ax=plt.subplots(figsize=(7,3.4))
kmax=int(x.max())+1; k=np.arange(0,kmax+1)
ax.hist(x,bins=np.arange(-0.5,kmax+1.5),density=True,color=AMBER,alpha=0.7,label="observed requests/sec")
ax.plot(k, stats.poisson.pmf(k,lam),"o-",color=TEAL,ms=4,label=f"Poisson(lambda={lam:.1f})")
ax.set_xlabel("requests per second"); ax.set_ylabel("probability"); ax.set_title("Requests per second: observed vs Poisson"); ax.legend()
plt.tight_layout(); plt.show()
# Formal goodness-of-fit: chi-square on observed vs expected counts.
# Rule of thumb: merge adjacent bins so every EXPECTED count is at least 5, then
# chi-square = sum((obs-exp)**2 / exp), with df = (bins - 1 - number of fitted parameters).
def chisq_gof(counts, pmf, n_params):
N = len(counts); kmax = int(counts.max())
obs = np.array([(counts == k).sum() for k in range(kmax + 1)], float)
exp = np.array([N * pmf(k) for k in range(kmax + 1)], float)
obs = np.append(obs, 0.0); exp = np.append(exp, N - exp.sum()) # tail bin for k > kmax
mo, me, co, ce = [], [], 0.0, 0.0
for o, e in zip(obs, exp):
co += o; ce += e
if ce >= 5:
mo.append(co); me.append(ce); co = ce = 0.0
if ce > 0 and mo:
mo[-1] += co; me[-1] += ce
mo, me = np.array(mo), np.array(me)
chi2 = ((mo - me) ** 2 / me).sum(); dof = len(mo) - 1 - n_params
return chi2, dof, stats.chi2.sf(chi2, dof)
chi2, dof, pval = chisq_gof(x.values, lambda k: stats.poisson.pmf(k, lam), 1)
print(f"Chi-square goodness-of-fit: chi2 = {chi2:.2f}, df = {dof}, p = {pval:.3f} -> "
+ ("cannot reject the fit (Poisson is a good model)" if pval > 0.05 else "reject: not Poisson"))
Chi-square goodness-of-fit: chi2 = 22.06, df = 17, p = 0.182 -> cannot reject the fit (Poisson is a good model)
Beyond the overlay, a chi-square goodness-of-fit test scores the Poisson against the observed per-second counts (rare bins merged to keep expected counts at least 5): chi-square = 22.1 on 17 df, p = 0.18, so we cannot reject it. The mean-equals-variance signature plus a passing GoF test makes the Poisson a defensible model.
Mean and variance are both about 12, the Poisson signature, and the PMF tracks the observed counts closely. A single rate lambda ≈ 12 requests/second describes the whole traffic distribution.
for thr in [15,18,20,22]:
print(f"P(X > {thr}) = {stats.poisson.sf(thr,lam):.4f} (observed {(x>thr).mean():.4f})")
print(f"\n95th percentile of load = {stats.poisson.ppf(0.95,lam):.0f} requests/sec")
print(f"99th percentile of load = {stats.poisson.ppf(0.99,lam):.0f} requests/sec")
P(X > 15) = 0.1561 (observed 0.1510) P(X > 18) = 0.0376 (observed 0.0300) P(X > 20) = 0.0117 (observed 0.0160) P(X > 22) = 0.0031 (observed 0.0050) 95th percentile of load = 18 requests/sec 99th percentile of load = 21 requests/sec
A breach above 20 requests/second happens about 1.2% of the time under the Poisson, P(X > 20) ≈ 0.012, close to the 1.6% seen in the data. The 99th percentile load is about 21 requests/second, the number to provision against. The tail, not the average, drives capacity.
n, p = 100000, lam/100000 # many users, each with a tiny per-second hit probability
k=np.arange(0,int(x.max())+1)
binom = stats.binom.pmf(k,n,p); pois = stats.poisson.pmf(k,lam)
print(f"max |Binomial(100000, lambda/100000) - Poisson(lambda)| = {np.abs(binom-pois).max():.6f}")
rng=np.random.default_rng(50); sim=rng.poisson(lam,size=200_000)
print(f"simulated mean {sim.mean():.2f}, var {sim.var():.2f}, P(>20) {(sim>20).mean():.4f}")
max |Binomial(100000, lambda/100000) - Poisson(lambda)| = 0.000007 simulated mean 12.02, var 11.97, P(>20) 0.0118
A binomial with 100,000 users each hitting at probability λ/100,000 is indistinguishable from the Poisson, differing by less than a millionth. That is the Poisson story: a vast number of independent chances, each unlikely, collapsing to a single rate λ that is both the mean and the variance of the count.