Probability Distributions – Expected Value and the Tails You’ll See

Goal: Choose the right distribution family for your security data, compute expected values safely, and understand why tail choice changes everything.

Why distribution choice matters

Cyber risk quantification almost always decomposes into frequency x severity: how often does the bad thing happen, and how bad is it each time? The annual loss expectation (ALE) is the product of these two expected values (when frequency and severity are independent), but the shape of each distribution determines what happens at the 95th percentile and beyond – and that tail is where budgets break and executives lose sleep.

A normal model on breach costs will tell you the 95th-percentile loss is perhaps 2x the mean. A lognormal model might say 5x. A Pareto model might say 20x. The decision – how much insurance to buy, how much to invest in controls – changes by an order of magnitude depending on which family you choose.

Setup

We use decision-security’s synth module for sampling and mixture for combining distributions.

import numpy as np
import matplotlib.pyplot as plt

from decision_security.synth import make_rng, sample, mixture

rng = make_rng(42)

plt.rcParams.update({
    "font.family": "serif",
    "font.size": 10,
    "axes.labelsize": 11,
    "axes.titlesize": 12,
    "xtick.labelsize": 9,
    "ytick.labelsize": 9,
    "legend.fontsize": 9,
    "figure.dpi": 150,
    "axes.spines.top": False,
    "axes.spines.right": False,
})

PRIMARY = "#1A1A1A"
ACCENT = "#E74C3C"
DARK_BG = "#34495E"
LIGHT_GRAY = "#95A5A6"
MED_GRAY = "#7F8C8D"
VERY_LIGHT = "#BDC3C7"

Expected Value and its Limits

The expected value \(\mathbb{E}[X]\) is the probability-weighted average outcome. Linearity of expectation – \(\mathbb{E}[aX + bY] = a\mathbb{E}[X] + b\mathbb{E}[Y]\) – holds regardless of dependence and is the workhorse behind ALE calculations. The product rule for ALE (Wald’s identity) needs one more condition: the number of incidents must be independent of their sizes:

\[\mathbb{E}[\text{loss}] = \mathbb{E}[\text{freq}] \times \mathbb{E}[\text{sev}]\]

If average loss per incident is \(\mathbb{E}[X]\) and you expect \(\lambda\) incidents, the expected aggregate loss is \(\lambda \cdot \mathbb{E}[X]\).

But EV is a center measure. Two portfolios can have the same EV and wildly different 95th percentiles. In security, decisions are usually about surviving the bad tail, not optimizing the average, so we need to always pair EV with quantile estimates.

Linearity applies to expectations, not to quantiles or tail measures. \(p_{95}(X + Y) \neq p_{95}(X) + p_{95}(Y)\) unless the variables are comonotonic (perfect rank correlation). This is why Monte Carlo simulation (Part 0.3) matters – you can’t shortcut tail calculations with algebra alone.

Note: The variance of the sum does depend on correlation – which matters for Monte Carlo modeling in Part 0.3.

X = make_rng(42).normal(loc=100, scale=25, size=50_000)

Y = 2 * X + 5

ev_direct = Y.mean()
ev_linear = 2 * X.mean() + 5

print(f"E[2X + 5]  (direct)  = {ev_direct:,.2f}")
print(f"2*E[X] + 5 (linear)  = {ev_linear:,.2f}")
print(f"Difference           = {abs(ev_direct - ev_linear):.6f}")

E[2X + 5]  (direct)  = 204.97
2*E[X] + 5 (linear)  = 204.97
Difference           = 0.000000

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

axes[0].hist(X, bins=60, alpha=0.7, color=DARK_BG, edgecolor="white")
axes[0].axvline(X.mean(), color=ACCENT, linestyle="--", linewidth=2, label=f"E[X] = {X.mean():.1f}")
axes[0].set_title("X ~ Normal(100, 25)")
axes[0].set_xlabel("X")
axes[0].legend()

axes[1].hist(Y, bins=60, alpha=0.7, color=ACCENT, edgecolor="white")
axes[1].axvline(Y.mean(), color=ACCENT, linestyle="--", linewidth=2, label=f"E[2X+5] = {Y.mean():.1f}")
axes[1].axvline(2 * X.mean() + 5, color=PRIMARY, linestyle=":", linewidth=2, label=f"2\u00b7E[X]+5 = {2*X.mean()+5:.1f}")
axes[1].set_title("Y = 2X + 5")
axes[1].set_xlabel("Y")
axes[1].legend()

plt.tight_layout()
plt.show()

Frequency: Poisson vs Negative Binomial

The Poisson distribution models the count of events in a fixed window (incidents per quarter, breaches per year) when events occur independently at a constant rate \(\lambda\). Its key property: mean equals variance (\(\mathbb{E}[X] = \text{Var}(X) = \lambda\)).

In practice, security event counts often exhibit overdispersion – the variance exceeds the mean – because the underlying rate is not constant. Threat campaigns cluster; vulnerability disclosures come in waves. The negative binomial distribution adds a dispersion parameter \(k\), allowing variance to exceed the mean.

A practical test: compute the sample mean and variance of your monthly incident counts. If \(s^2 / \bar{x} > 1.5\), try negative binomial. If it is close to 1, Poisson suffices. Check for a trend first: a rate that rises or falls over the window pushes the ratio above 1.5 on its own, even without clustering (Part 3.4 tests this on real data). Short series fire the rule too: with 12 months of pure Poisson counts at 2.5 a month, the ratio exceeds 1.5 about 1 time in 8.

Note: Ransomware campaigns illustrate overdispersion well. In quiet months you see zero or one event; when a new variant spreads, you might see a cluster of five in a single week. The negative binomial captures this “bursty” pattern that Poisson cannot.

freq_pois   = sample('poisson', 10000, lam=2.5, rng=make_rng(11))
freq_negbin = sample('negbin',  10000, mean=2.5, k=3.0, rng=make_rng(12))

fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)

bins = np.arange(0, 20)

axes[0].hist(freq_pois, bins=bins, align='left', color=DARK_BG, edgecolor='white')
axes[0].set_title(f'Poisson (lam=2.5)\nmean={freq_pois.mean():.2f}, var={freq_pois.var():.2f}')
axes[0].set_xlabel('Events per year')
axes[0].set_ylabel('Count')

axes[1].hist(freq_negbin, bins=bins, align='left', color=ACCENT, edgecolor='white')
axes[1].set_title(f'NegBin (mean=2.5, k=3)\nmean={freq_negbin.mean():.2f}, var={freq_negbin.var():.2f}')
axes[1].set_xlabel('Events per year')

plt.tight_layout()
plt.show()

vmr_pois = freq_pois.var() / freq_pois.mean()
vmr_nb = freq_negbin.var() / freq_negbin.mean()
print(f"Poisson  p95 = {np.percentile(freq_pois, 95):.0f}   variance-to-mean ratio = {vmr_pois:.2f}")
print(f"NegBin   p95 = {np.percentile(freq_negbin, 95):.0f}   variance-to-mean ratio = {vmr_nb:.2f}")

Poisson  p95 = 5   variance-to-mean ratio = 1.00
NegBin   p95 = 7   variance-to-mean ratio = 1.83

Severity: Lognormal vs Pareto

Both are right-skewed with long tails, but the Pareto tail is much heavier.

Lognormal (soft tail): the default starting point in cyber risk. It arises naturally when losses are the product of many multiplicative factors (scope of breach \(\times\) records affected \(\times\) regulatory multiplier \(\times\) …). The lognormal is parameterized by \(\mu\) and \(\sigma\) on the log scale: the median is \(e^{\mu}\) and quantiles follow:

\[Q_p = \exp(\mu + \sigma \, z_p)\]

where \(z_p\) is the standard normal quantile. Most breach cost data fits here. The log-transformed values are roughly normal, so the distribution has a well-behaved mean.

Pareto (hard tail): models the “power law” regime where a small number of events account for most of the total loss. Empirically, the largest cyber losses (nation-state attacks, massive data breaches) exhibit Pareto-like behavior. The key diagnostic is the tail index \(\alpha\): when \(\alpha \leq 2\), variance is infinite and sample means are unstable.

When to choose which: Use lognormal as the default for severity. Escalate to Pareto (or a generalized Pareto tail splice) when you have evidence of extreme tail weight – for instance, when the top 5% of losses account for more than 50% of the total. Another diagnostic: \(p_{99}/p_{95}\) much larger than \(p_{95}/p_{90}\). Both flag a heavy tail, not the Pareto family itself (a lognormal with a large \(\sigma\) passes them too). To tell the two apart, check whether the log-log survival plot is a straight line, or run a likelihood-ratio test of Pareto against lognormal.

In the matched-median example below, both distributions share a $10K median, yet moving from lognormal to Pareto (\(\alpha = 1.1\)) lifts the p95 by roughly 1.6x and the p99 by roughly 3.4x — and the log-log survival plot shows the gap widening the further you go into the tail.

meanlog, sdlog = np.log(10_000), 1.0
alpha = 1.1
xm_pareto = 10_000 / 2 ** (1 / alpha)   # pick xm so the Pareto median matches the lognormal median ($10K)

sev_logn   = sample('lognormal', 100_000, meanlog=meanlog, sdlog=sdlog, rng=make_rng(1))
sev_pareto = sample('pareto',    100_000, alpha=alpha, xm=xm_pareto, rng=make_rng(2))

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

axes[0].hist(np.log10(sev_logn), bins=60, color=DARK_BG, edgecolor='white', alpha=0.6, label='Lognormal')
axes[0].hist(np.log10(sev_pareto), bins=60, color=ACCENT, edgecolor='white', alpha=0.6, label='Pareto')
axes[0].axvline(np.log10(10_000), color=PRIMARY, ls='--', lw=1.5, label='shared median $10K')
axes[0].set_title('Matched medians, different tails')
axes[0].set_xlabel('log10(loss $)')
axes[0].set_ylabel('Count')
axes[0].legend(fontsize=8)

# A Pareto tail is a straight line on log-log axes
for label, arr, color in [('Lognormal', sev_logn, DARK_BG), ('Pareto', sev_pareto, ACCENT)]:
    xs = np.sort(arr)
    ccdf = 1.0 - np.arange(1, len(xs) + 1) / len(xs)
    axes[1].loglog(xs, ccdf, color=color, lw=1.5, label=label)
axes[1].set_title('Survival P(Loss > x): Pareto tail dominates')
axes[1].set_xlabel('loss $ (log)')
axes[1].set_ylabel('P(Loss > x) (log)')
axes[1].legend(fontsize=8)

plt.tight_layout()
plt.show()

for label, arr in [('Lognormal', sev_logn), ('Pareto', sev_pareto)]:
    print(f"{label:10s}  median=${np.median(arr):>14,.0f}   "
          f"p95=${np.percentile(arr, 95):>14,.0f}   "
          f"p99=${np.percentile(arr, 99):>14,.0f}")

p95_ratio = np.percentile(sev_pareto, 95) / np.percentile(sev_logn, 95)
p99_ratio = np.percentile(sev_pareto, 99) / np.percentile(sev_logn, 99)
print(f"\nPareto / lognormal:  p95 {p95_ratio:.1f}x   p99 {p99_ratio:.1f}x")

Lognormal   median=$         9,946   p95=$        51,478   p99=$       101,905
Pareto      median=$         9,990   p95=$        81,925   p99=$       349,632

Pareto / lognormal:  p95 1.6x   p99 3.4x

Expert Elicitation: Lognormal from Two Quantiles

When historical data is sparse, you can calibrate a lognormal from expert-elicited quantiles. Experts rarely think in \(\mu_{\log}\) and \(\sigma_{\log}\). They say things like: “median loss is about $50K, and it could go as high as $2M in a bad case.”

Ask the subject-matter expert: “What is the median loss?” and “What is a reasonable worst case – the value exceeded only 5% of the time?” These give you two equations:

\[\text{median} = e^{\mu}, \quad Q_{0.95} = \exp(\mu + 1.645\,\sigma)\]

Solving:

\[\mu = \ln(\text{median}), \quad \sigma = \frac{\ln(Q_{0.95}) - \ln(\text{median})}{1.645}\]

(Note: FAIR itself uses a modified Beta-PERT fitted from min, most likely, and max — a different elicitation format. The lognormal approach here is an alternative when you want a heavier tail or prefer quantile-based calibration.) The from_quantiles option below takes two quantile-value pairs and back-solves for the lognormal parameters automatically.

Tip: Train your experts to think in orders of magnitude. “Is the worst case 3x the median or 30x?” maps directly to \(\sigma\) – a ratio of 3 implies \(\sigma \approx 0.67\); a ratio of 30 implies \(\sigma \approx 2.07\).

sev_expert = sample(
    'lognormal', 10000,
    from_quantiles=True,
    x1=50_000,  p1=0.5,
    x2=2_000_000, p2=0.95,
    rng=make_rng(42),
)

sigma_implied = (np.log(2_000_000) - np.log(50_000)) / 1.645
p99_implied = 50_000 * np.exp(2.326 * sigma_implied)
print(f"p95/median ratio: {2_000_000/50_000:.0f}x  ->  implied sigma (log-scale) = {sigma_implied:.2f}")
print(f"Implied p99 loss: ${p99_implied:,.0f}  (sample p99 = ${np.percentile(sev_expert, 99):,.0f})")

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Clipped at p95, otherwise the tail squashes the bulk
axes[0].hist(sev_expert / 1e3, bins=80, color=DARK_BG, edgecolor='white',
             range=(0, np.percentile(sev_expert, 95) / 1e3))
axes[0].axvline(np.median(sev_expert) / 1e3, color=PRIMARY, ls='--',
                label=f'median = USD {np.median(sev_expert)/1e3:,.0f}K')
axes[0].axvline(np.percentile(sev_expert, 95) / 1e3, color=ACCENT, ls='--',
                label=f'p95 = USD {np.percentile(sev_expert, 95)/1e3:,.0f}K')
axes[0].set_xlabel('Loss (USD K)')
axes[0].set_ylabel('Count')
axes[0].set_title('Bulk of the distribution (up to p95)')
axes[0].legend(fontsize=9)

axes[1].hist(np.log10(sev_expert), bins=60, color=DARK_BG, edgecolor='white')
for pct, col, ls, lbl in [(50, PRIMARY, '--', 'median'), (95, ACCENT, '--', 'p95'), (99, ACCENT, ':', 'p99')]:
    val = np.percentile(sev_expert, pct)
    axes[1].axvline(np.log10(val), color=col, ls=ls,
                    label=f'{lbl} = USD {val/1e3:,.0f}K')
axes[1].set_xlabel('log10(Loss)')
axes[1].set_ylabel('Count')
axes[1].set_title('Full range on log scale — the tail is real')
axes[1].legend(fontsize=9)

plt.tight_layout()
plt.show()

p95/median ratio: 40x  ->  implied sigma (log-scale) = 2.24
Implied p99 loss: $9,209,872  (sample p99 = $9,851,897)

Mixtures: Bimodal Breach Data

Real breach-cost data often looks bimodal: a large cluster of routine incidents (credential stuffing, single-endpoint malware, small data exposure, contained quickly) and a small cluster of catastrophic events (mass data exfiltration, ransomware with extended downtime, major litigation, regulatory action). No single lognormal fits both modes well.

A mixture model – say, 90% lognormal(\(\mu_1, \sigma_1\)) + 10% lognormal(\(\mu_2, \sigma_2\)) – captures this structure. The mixture inflates the tail beyond what the routine component alone would produce. In practice, even a two-component mixture substantially improves p95 estimates compared to a single fitted distribution.

A mixture model combines two (or more) component distributions with mixing weights.

rng = make_rng(42)

mixed = mixture(
    components=[
        ('lognormal', {'meanlog': 9, 'sdlog': 0.8}),
        ('lognormal', {'meanlog': 14, 'sdlog': 1.0}),
    ],
    weights=[0.9, 0.1],
    size=10000,
    rng=rng,
)

fig, ax = plt.subplots(figsize=(10, 4))
ax.hist(np.log10(mixed), bins=80, color=DARK_BG, edgecolor="white", alpha=0.8)
ax.set_xlabel(r"Log$_{10}$ Loss (\$)")
ax.set_ylabel("Count")
ax.set_title("Mixture: 90% Routine + 10% Catastrophic")
ax.axvline(np.log10(np.median(mixed)), color=ACCENT, linestyle="--",
           linewidth=1.5, label=f"Median = ${np.median(mixed):,.0f}")
ax.axvline(np.log10(np.percentile(mixed, 95)), color=PRIMARY, linestyle=":",
           linewidth=1.5, label=f"p95 = ${np.percentile(mixed, 95):,.0f}")
ax.legend()
plt.tight_layout()
plt.show()

Model Selection Cheat Sheet

Data type	Default family	Upgrade when…
Event counts (stable rate)	Poisson	Overdispersion (\(s^2/\bar{x} > 1.5\), after ruling out a trend) –> negative binomial
Event counts (many zeros)	Zero-inflated Poisson	Zero-inflated negative binomial
Loss per event	Lognormal	Extreme tail weight –> Pareto or GPD splice
Proportions	Beta	N/A (use Wilson CI for inference)
Bimodal losses	Lognormal mixture	Formal mixture model with EM

Beta vs Wilson — two answers to the same question. The Beta distribution is the Bayesian model for a proportion: with a \(\text{Beta}(\alpha, \beta)\) prior, observing \(s\) successes in \(n\) trials gives a \(\text{Beta}(\alpha + s,\ \beta + n - s)\) posterior you can read credible intervals off directly. The Wilson score interval is the frequentist counterpart – a closed-form confidence interval that needs no prior and stays well-behaved at extreme proportions and small \(n\). Reach for Beta when you have prior information worth encoding (or want a full posterior); reach for Wilson when you just need a quick, prior-free interval on an observed rate.

Common Traps

1. Normal on skewed data. The normal distribution allows negative values and has thin tails. Applying it to loss data truncates the tail and produces nonsensical negative-loss scenarios. Always use a right-skewed distribution (lognormal, gamma, Pareto) for positive-only, skewed severity quantities.

2. Estimating EV from too few tail samples. With 30 observations you can estimate the median reliably but the 95th percentile poorly. With 1,000 draws, your p99 estimate is based on roughly 10 observations. If you have fewer than ~100 data points, treat tail estimates as highly uncertain and report wide credible intervals. For stable tail estimates from a simulation, use at least 10,000 draws; more draws reduce simulation noise, not the uncertainty from fitting on few data points – and even then, treat extreme quantiles (p99.9) with skepticism.

3. Ignoring zero-inflation. Many security count series have excess zeros (months with no events of interest or incidents of a particular type). A Poisson with a low \(\lambda\) already handles some zero mass, and a negative binomial puts extra mass on zero, but if the zero mass is larger than they predict, standard Poisson or negative binomial will underfit the zero mass and overfit the positive counts. Test for zero-inflation and use ZIP/ZINB models if warranted, or a mixture with a point mass at zero.

Bayesian Preview: Flipping the Probability Tree

Raiffa (1968) introduced one of the most intuitive explanations of Bayes’ theorem by showing how to mechanically “flip” a probability tree. The forward tree starts with the state of the world (intrusion or benign) and branches to the observable evidence (alert fires or not). The flipped tree reverses the order: it starts with the evidence you actually observe (alert fired) and branches to the posterior probability of each state.

The flipping procedure is mechanical: compute each path’s joint probability as the product of its branch probabilities, sum the paths that lead to the same evidence outcome, and divide each joint by that sum. The result is the posterior (no formula memorization required).

Security scenario: An IDS alert fires. You need to decide if this is a real intrusion or a false positive.

• Prior: 2% of events are real intrusions – the base rate across all traffic, before the IDS fires. (This is not the fraction of alerts that are real; that is the posterior we compute below.)

• Detection: If it’s a real intrusion, the IDS fires 90% of the time (sensitivity).

• False alarm: If it’s benign, the IDS still fires 5% of the time (false positive rate).

The forward tree gives: \(P(\text{alert} \mid \text{intrusion}) = 0.90\), \(P(\text{alert} \mid \text{benign}) = 0.05\). Flipping the tree gives: \(P(\text{intrusion} \mid \text{alert}) = ?\)

Most analysts intuitively guess 50–80%. The Bayesian answer is approximately 27% – because the 98% base rate of benign traffic, even with only a 5% false alarm rate, generates far more false alerts than the 2% intrusion rate generates true ones.

This calculation has direct operational implications. It explains why even high-quality detection systems produce alert volumes dominated by false positives, why SOC teams develop “alert fatigue,” and why base-rate awareness is the single most important statistical concept for anyone triaging security events.

p_intrusion = 0.02       # base rate over all events, not over alerts
p_benign = 1 - p_intrusion

p_alert_given_intrusion = 0.90   # sensitivity
p_alert_given_benign = 0.05      # false positive rate

p_intrusion_and_alert = p_intrusion * p_alert_given_intrusion
p_benign_and_alert = p_benign * p_alert_given_benign

p_alert = p_intrusion_and_alert + p_benign_and_alert

p_intrusion_given_alert = p_intrusion_and_alert / p_alert

print("=== Forward tree ===")
print(f"  P(intrusion AND alert) = {p_intrusion:.2f} x {p_alert_given_intrusion:.2f} = {p_intrusion_and_alert:.4f}")
print(f"  P(benign AND alert)    = {p_benign:.2f} x {p_alert_given_benign:.2f} = {p_benign_and_alert:.4f}")
print(f"  P(alert)               = {p_alert:.4f}")
print()
print("=== Flipped tree ===")
print(f"  P(intrusion | alert)   = {p_intrusion_and_alert:.4f} / {p_alert:.4f} = {p_intrusion_given_alert:.1%}")
print(f"  P(benign | alert)      = {1 - p_intrusion_given_alert:.1%}")
print()
if p_intrusion_given_alert < 0.5:
    print(f"Most alerts are false positives: the {p_intrusion:.0%} base rate outweighs {p_alert_given_intrusion:.0%} sensitivity.")
else:
    print(f"Most alerts are real: an FPR of {p_alert_given_benign:.1%} overcomes the {p_intrusion:.0%} base rate.")

=== Forward tree ===
  P(intrusion AND alert) = 0.02 x 0.90 = 0.0180
  P(benign AND alert)    = 0.98 x 0.05 = 0.0490
  P(alert)               = 0.0670

=== Flipped tree ===
  P(intrusion | alert)   = 0.0180 / 0.0670 = 26.9%
  P(benign | alert)      = 73.1%

Most alerts are false positives: the 2% base rate outweighs 90% sensitivity.