Part 0.1 — Stats 101

Goal: Summarize security data quickly and correctly (means, medians, quantiles, rates) and know when the “obvious” summary will mislead you.

Why summary statistics matter for security decisions

Every security decision rests on a number somebody summarized somewhere, from alert volumes that decide staffing levels to Mean-time-to-respond (MTTR) metrics that appear in board decks to breach-cost estimates that feed risk registers. If the summary is wrong or if it uses the wrong measure then the decision it supports is wrong too.

Security data is rarely well-behaved. Response times are right-skewed, loss amounts span orders of magnitude, and click rates cluster near zero. The standard “take the average” reflex can be actively dangerous in these cases. This notebook covers the summary statistics that work under skew and heavy tails, and shows when each one breaks down.

None of this requires advanced statistics. Instead it require discipline, picking the right summary for the shape of your data, attaching an uncertainty interval, and knowing which traps to avoid.

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
from decision_security.synth import make_rng, sample

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"


def annotated_hist(ax, data, lines, bins=60, color=DARK_BG, xlabel="", title=""):
    """Histogram with labeled vertical guide lines.

    lines: iterable of (label, value, color) tuples drawn as dashed verticals.
    """
    ax.hist(data, bins=bins, edgecolor="white", alpha=0.7, color=color)
    for label, val, lc in lines:
        ax.axvline(val, color=lc, linestyle="--", linewidth=1.5, label=label)
    ax.set_xlabel(xlabel)
    ax.set_ylabel("Count")
    if title:
        ax.set_title(title)
    ax.legend()
    return ax

Location: Mean vs Median

The mean and median both answer “what’s typical?” but they disagree when the data is skewed.

The arithmetic mean \(\bar x = \frac{1}{n}\sum x_i\) weights every observation equally, which makes it sensitive to outliers. Consider MTTR across incidents: most tickets close within a day, but a handful of complex investigations run for weeks. In the simulated sample below the mean is about 57 hours while the median sits at 20 hours. A staffing model built on the mean over-provisions; an SLA built on it looks easy to hit until the next tail event.

The median is the middle value when you sort the data. Half the observations fall above it, half below. It ignores how extreme the extremes are — whether your worst incident took 200 hours or 2,000, the median doesn’t move.

The trimmed mean — dropping the top and bottom \(k\)% before averaging — offers a middle ground. A 10% trimmed mean on MTTR discards the worst outliers while still reflecting more of the distribution than the median alone.

In security, almost everything is right-skewed: response times, dwell times, breach costs, vulnerability age (roughly lognormal). For these distributions:

• The median tells you the typical experience — what most incidents actually look like

• The mean tells you the budget expectation — what you’d spend per incident on average if you had to plan for the total

Use the median when the distribution is skewed or heavy-tailed. Use the mean when the distribution is roughly symmetric (rare) or when you genuinely need the total divided by count (e.g., average cost per user for budget allocation). A good habit: report both, plus the ratio. When mean/median > 2, you are deep in skew territory and the mean alone is misleading.

For time-to-event metrics like MTTR or mean-time-to-detect, even the median can mislead — survival curves handle censoring and reveal the full resolution trajectory (see Part 0.7).

Note: On lognormal data the median equals \(e^{\mu}\) and is a far more stable summary than the mean.

response_hours = sample('lognormal', 500, meanlog=3.0, sdlog=1.5, rng=make_rng(42))

mean_val = np.mean(response_hours)
median_val = np.median(response_hours)

print(f"Mean response time:   {mean_val:.1f} hours")
print(f"Median response time: {median_val:.1f} hours")
mean_median_ratio = mean_val / median_val
if mean_median_ratio > 2:
    print(f"\nThe mean is {mean_median_ratio:.1f}x the median — pulled by outliers.")
else:
    print(f"\nThe mean is {mean_median_ratio:.1f}x the median — the two agree closely.")

Mean response time:   57.3 hours
Median response time: 20.2 hours

The mean is 2.8x the median — pulled by outliers.

fig, ax = plt.subplots(figsize=(9, 4))
annotated_hist(
    ax, response_hours,
    lines=[(f"Mean = {mean_val:.0f} h", mean_val, ACCENT),
           (f"Median = {median_val:.0f} h", median_val, PRIMARY)],
    xlabel="Response Time (hours)",
    title="Incident Response Times: Mean vs Median",
)
plt.tight_layout()
plt.show()

Spread: Variance, Standard Deviation, IQR

Knowing the center isn’t enough — you need to also know how wide the data spreads.

Variance \(s^2 = \frac{1}{n-1}\sum (x_i - \bar x)^2\) and its square root, the standard deviation, measure dispersion symmetrically around the mean. When the data is skewed, variance is inflated by the same tail observations that distort the mean. The interquartile range (IQR = \(Q_{75} - Q_{25}\)) is a robust alternative: it describes the spread of the middle 50% and is unaffected by extreme tails.

For security dashboards, pair the median with the IQR for a resilient center-and-spread summary. Reserve standard deviation for contexts where the data is approximately symmetric — for instance, login-attempt counts per hour during normal operations. If the std dev is larger than the median, that’s a strong signal that the mean-based summary is misleading.

Another useful measure is the coefficient of variation (CV = \(s / \bar x\)), which expresses spread relative to the mean. A CV above 1 signals high relative variability — common for breach costs and dwell times. When comparing variability across metrics with different units or scales (e.g., MTTR in hours vs loss in dollars), the CV provides a unit-free comparison that standard deviation alone cannot.

Note: The median absolute deviation (MAD) is another robust spread measure: \(\text{MAD} = \text{median}(|x_i - \text{median}(x)|)\). Scaled by 1.4826, it estimates standard deviation under normality but resists outlier contamination. It is a good alternative to IQR when you want a single-number robust spread.

var_val = np.var(response_hours, ddof=1)
std_val = np.std(response_hours, ddof=1)
q25, q75 = np.percentile(response_hours, [25, 75])
iqr_val = q75 - q25

print(f"Variance:  {var_val:,.1f}")
print(f"Std Dev:   {std_val:,.1f}")
print(f"IQR:       {iqr_val:.1f}  (Q25={q25:.1f}, Q75={q75:.1f})")
if std_val > median_val:
    print(f"\nStd dev ({std_val:.0f}) is larger than the median ({median_val:.0f}) — the distribution is skewed.")
else:
    print(f"\nStd dev ({std_val:.0f}) is smaller than the median ({median_val:.0f}) — moderate spread.")

Variance:  18,887.3
Std Dev:   137.4
IQR:       41.1  (Q25=7.4, Q75=48.5)

Std dev (137) is larger than the median (20) — the distribution is skewed.

fig, ax = plt.subplots(figsize=(9, 3))
bp = ax.boxplot(response_hours, orientation="horizontal", widths=0.6, patch_artist=True,
                boxprops=dict(facecolor=DARK_BG, alpha=0.3),
                medianprops=dict(color=PRIMARY, linewidth=2),
                meanprops=dict(marker="D", markerfacecolor=ACCENT, markeredgecolor="white", markersize=8),
                showmeans=True, showfliers=False)
jitter = make_rng(7).normal(1, 0.04, len(response_hours))
ax.scatter(response_hours, jitter, alpha=0.15, s=8, color=DARK_BG)
ax.set_xlabel("Response Time (hours)")
ax.set_title("Spread: Box Plot with IQR (box) and Mean (red diamond)")
ax.set_yticks([])
plt.tight_layout()
plt.show()

Quantiles as Guardrails

Quantiles answer “what value is exceeded only \((100-p)\)% of the time?” The 95th percentile (\(Q_{0.95}\)) of incident response time tells you the threshold below which 95% of incidents resolve. This maps directly to SLA language, budget reserves, and staffing ceilings.

A single central number hides tail risk. “Median response is 20 hours” sounds fine, but if p95 is 223 hours, 1 in 20 incidents takes over 9 days. SLAs are commitments about the worst acceptable case, not the typical case, and that’s a quantile, not a mean.

In practice, p50 (median) anchors “typical,” p90 flags the stress case, and p95 sets the worst-case budget line. Where you set a cyber-insurance retention depends on how much loss you can absorb; a retention near the p95 of your modeled annual loss means you insure only the catastrophic tail.

Why quantiles map to decisions:

• Risk appetite is set at a loss threshold (“we can absorb up to $X”) — that’s a percentile

• “95% of incidents resolved within 72 hours” directly defines a p95 target

• The gap between your current p95 and target p95 = the operational improvement needed

How to use them:

• SLA: “Our SLA is 48-hour response” \(\approx\) setting a target at ~p75 of current distribution

• Budget: “There is a 5% chance annual loss exceeds $2M” — that is the p95 of the loss distribution

• Staffing: If p90 concurrent incidents = 5, your shift model needs to handle that volume without queue buildup

Tip: The “five-number summary” (min, Q25, median, Q75, max) plus p95 should be your default descriptive output for any security metric. It fits on one line and immediately reveals skew, spread, and tail risk.

p50 = np.percentile(response_hours, 50)
p90 = np.percentile(response_hours, 90)
p95 = np.percentile(response_hours, 95)

print(f"p50 (median): {p50:.1f} hours")
print(f"p90:          {p90:.1f} hours")
print(f"p95:          {p95:.1f} hours")

fig, ax = plt.subplots(figsize=(9, 4))
annotated_hist(
    ax, response_hours,
    lines=[(f"p50 = {p50:.0f} h", p50, PRIMARY),
           (f"p90 = {p90:.0f} h", p90, MED_GRAY),
           (f"p95 = {p95:.0f} h", p95, ACCENT)],
    xlabel="Response Time (hours)",
    title="Incident Response Time Distribution with Quantile Guardrails",
)
plt.tight_layout()
plt.show()

p50 (median): 20.2 hours
p90:          122.9 hours
p95:          223.0 hours

Rates and Proportions

Phishing click rates, patch compliance percentages, and detection rates are all proportions. A point estimate is not enough; you need a confidence interval, especially with small samples.

A single campaign’s click rate is noisy. If 16 out of 200 people clicked (8%), is your organisation’s true susceptibility 8%? Or could it be 5% on a good month and 12% on a bad one? Without a CI, you cannot tell whether a drop from 8% to 5% next quarter is a real improvement or sampling noise. The CI gives you the range of rates consistent with the data. Next quarter’s rate is noisy too, so to judge a change, compare the two rates with a test that accounts for both, not one rate against the other’s band.

The naive confidence interval \(\hat p \pm z \sqrt{\hat p(1-\hat p)/n}\) (Wald interval) fails badly when \(\hat p\) is near 0 or 1 — exactly where security proportions tend to live. A campaign with 2 clicks out of 500 emails gives \(\hat p = 0.004\); the Wald interval can go negative.

The Wilson score interval corrects this by inverting the score test, producing valid intervals even at extreme proportions and small samples. Always use Wilson (or Agresti–Coull) for click rates, false-positive rates, and compliance percentages.

# 12 monthly campaigns with a fixed true click rate of 8%: month-to-month variation is pure noise
n_sent = 200
n_campaigns = 12
clicks = sample('binomial', n_campaigns, n=n_sent, p=0.08, rng=make_rng(42))

click_rates = clicks / n_sent
print("Monthly click counts:", clicks)
print("Click rates:         ", np.round(click_rates, 3))
print(f"Overall mean rate:    {np.mean(click_rates):.3f}")

Monthly click counts: [19 15 20 18 11 24 19 19 12 15 15 22]
Click rates:          [0.095 0.075 0.1   0.09  0.055 0.12  0.095 0.095 0.06  0.075 0.075 0.11 ]
Overall mean rate:    0.087

# Wilson CI for the pooled rate
total_clicks = int(np.sum(clicks))
total_sent = n_sent * len(clicks)
p_hat = total_clicks / total_sent

z = scipy.stats.norm.ppf(0.975)  # 95% CI

denominator = 1 + z**2 / total_sent
centre = (p_hat + z**2 / (2 * total_sent)) / denominator
margin = (z / denominator) * np.sqrt(p_hat * (1 - p_hat) / total_sent + z**2 / (4 * total_sent**2))

wilson_lo = centre - margin
wilson_hi = centre + margin

print(f"Observed rate:    {p_hat:.4f}")
print(f"Wilson 95% CI:    [{wilson_lo:.4f}, {wilson_hi:.4f}]")
print(f"\nAssumes all {total_sent:,} sends are i.i.d.; months that are genuinely riskier would widen this.")

Observed rate:    0.0871
Wilson 95% CI:    [0.0765, 0.0990]

Assumes all 2,400 sends are i.i.d.; months that are genuinely riskier would widen this.

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

# Left: click rates per campaign. The band is where ONE 200-email campaign lands 95% of
# the time if the rate never changes; the pooled Wilson CI is for the overall rate and far narrower.
camp_lo, camp_hi = scipy.stats.binom.ppf([0.025, 0.975], n_sent, p_hat) / n_sent
n_inside = int(np.sum((click_rates >= camp_lo) & (click_rates <= camp_hi)))
print(f"Per-campaign 95% band: [{camp_lo:.1%}, {camp_hi:.1%}]  ({n_inside} of {len(click_rates)} campaigns inside)")

months = np.arange(1, len(click_rates) + 1)
axes[0].bar(months, click_rates, color=DARK_BG, edgecolor="white")
axes[0].axhline(p_hat, color=ACCENT, linestyle="--", linewidth=1.5, label=f"Overall rate = {p_hat:.1%}")
axes[0].axhspan(camp_lo, camp_hi, alpha=0.15, color=ACCENT, label="95% band for one campaign")
axes[0].set_xlabel("Campaign")
axes[0].set_ylabel("Click Rate")
axes[0].set_title("Phishing Click Rate per Campaign")
axes[0].legend(fontsize=9)

# Right: distribution of click counts
vals, counts = np.unique(clicks, return_counts=True)
axes[1].bar(vals, counts, color=DARK_BG, edgecolor="white")
axes[1].set_xlabel("Clicks per Campaign (out of 200)")
axes[1].set_ylabel("Frequency")
axes[1].set_title("Distribution of Click Counts")

plt.tight_layout()
plt.show()

Per-campaign 95% band: [5.0%, 13.0%]  (12 of 12 campaigns inside)

In the left panel, every campaign bar that falls inside the shaded band (the range one 200-email campaign hits 95% of the time if the click rate never changes) is consistent with the same underlying click rate — the month-to-month variation is noise, not signal. A bar outside the band would warrant investigation.

Poisson Counts

Security incidents per time period — alerts per shift, breaches per quarter, phishing campaigns per month — are counts. The Poisson distribution models these: events arrive independently at a constant average rate \(\lambda\).

For count data the Poisson confidence interval applies. With \(k\) observed events the exact interval uses the chi-squared distribution; for \(k \geq 20\) the normal approximation \(k \pm z\sqrt{k}\) is adequate.

If you observed 8 incidents this month, the 95% CI on the rate is roughly [3.5, 15.8]. But next month’s count is noisy too, so to judge a change, compare the two counts directly: if the rate never changed, a split at least as uneven as 8 vs 12 happens about half the time (p ≈ 0.50).

Note: A SOC that observes 8 critical incidents in January and 12 in February might panic about a 50% increase. But comparing the two months directly gives p ≈ 0.50 — February’s 12 is well within the noise.

from scipy.stats import binomtest, chi2, poisson

monthly_incidents = sample('poisson', 24, lam=3.2, rng=make_rng(42))

lambda_hat = np.mean(monthly_incidents)
n_months = len(monthly_incidents)
total_events = int(np.sum(monthly_incidents))

# Exact (gamma / chi-square) 95% CI for the MEAN rate, based on the total count.
rate_lo = chi2.ppf(0.025, 2 * total_events) / 2 / n_months
rate_hi = chi2.ppf(0.975, 2 * (total_events + 1)) / 2 / n_months

# Single-month PREDICTIVE band: where one month's count is expected to fall,
# given the estimated rate.
pred_lo, pred_hi = poisson.ppf([0.025, 0.975], lambda_hat)

print(f"Observed monthly incidents: {monthly_incidents}")
print(f"Estimated mean rate (lambda):     {lambda_hat:.2f}  ({total_events} events over {n_months} months)")
print(f"95% CI of the mean rate:          [{rate_lo:.2f}, {rate_hi:.2f}]  (exact, chi-square)")
print(f"95% single-month predictive band: [{pred_lo:.0f}, {pred_hi:.0f}] incidents")

# The "8 then 12" example from the text.
k = 8
k_lo = chi2.ppf(0.025, 2 * k) / 2
k_hi = chi2.ppf(0.975, 2 * (k + 1)) / 2
print(f"\nExact 95% CI on the rate for a month with {k} events: [{k_lo:.1f}, {k_hi:.1f}]")
# Both months are noisy, so test the two counts against each other: under a constant
# rate, each of the k + next_month events falls in either month with probability 1/2.
next_month = 12
p_change = binomtest(next_month, k + next_month, 0.5).pvalue
if p_change < 0.05:
    print(f"{k} vs {next_month}: p = {p_change:.2f} — more than noise alone explains.")
else:
    print(f"{k} vs {next_month}: p = {p_change:.2f} — consistent with noise.")

fig, ax = plt.subplots(figsize=(9, 4))
months = np.arange(1, n_months + 1)
ax.bar(months, monthly_incidents, color=DARK_BG, edgecolor="white")
ax.axhline(lambda_hat, color=ACCENT, linestyle="--", linewidth=1.5, label=f"Mean rate = {lambda_hat:.1f}")
ax.fill_between(months, rate_lo, rate_hi, alpha=0.2, color=ACCENT,
                label=f"95% CI on mean rate [{rate_lo:.1f}, {rate_hi:.1f}]")
ax.axhline(pred_lo, color=LIGHT_GRAY, linestyle=":", linewidth=1, label=f"95% single-month band [{pred_lo:.0f}, {pred_hi:.0f}]")
ax.axhline(pred_hi, color=LIGHT_GRAY, linestyle=":", linewidth=1)
ax.set_xlabel("Month")
ax.set_ylabel("Incident Count")
ax.set_title("Monthly Security Incidents")
ax.legend()
plt.tight_layout()
plt.show()

Observed monthly incidents: [4 4 5 1 7 1 4 2 2 6 3 1 6 2 5 0 2 2 3 4 0 2 2 3]
Estimated mean rate (lambda):     2.96  (71 events over 24 months)
95% CI of the mean rate:          [2.31, 3.73]  (exact, chi-square)
95% single-month predictive band: [0, 7] incidents

Exact 95% CI on the rate for a month with 8 events: [3.5, 15.8]
8 vs 12: p = 0.50 — consistent with noise.

Five-Number Summary + p95

The five-number summary (min, Q25, median, Q75, max) plus the 95th percentile should be your default descriptive output for any security metric. It fits on one line and immediately reveals skew (is the median much closer to Q25 than Q75?), spread (how wide is the IQR?), and tail risk (how far is p95 from the median?).

def five_number_plus(data, label="value"):
    """Return a five-number summary plus p95."""
    keys = ["min", "Q1", "median", "Q3", "p95", "max"]
    vals = [
        np.min(data),
        np.percentile(data, 25),
        np.median(data),
        np.percentile(data, 75),
        np.percentile(data, 95),
        np.max(data),
    ]
    print(f"Five-number summary + p95 for {label}:")
    for k, v in zip(keys, vals):
        print(f"  {k:>6s}: {v:>10.1f}")
    return dict(zip(keys, vals))

_ = five_number_plus(response_hours, label="response time (hours)")

Five-number summary + p95 for response time (hours):
     min:        0.4
      Q1:        7.4
  median:       20.2
      Q3:       48.5
     p95:      223.0
     max:     1588.9

Putting It Together: A Descriptive Workflow

When you encounter a new security dataset, follow this sequence:

1. Compute the five-number summary plus p95

2. Check the mean/median ratio to assess skew

3. Choose median + IQR or mean + SD based on skew

4. For proportions, compute Wilson CIs

5. For counts, compute Poisson CIs

This five-step routine takes under a minute and gives you a defensible, honest summary before any modeling begins.

A well-chosen summary statistic paired with an appropriate uncertainty interval is worth more than a complex model built on a poorly understood dataset.

Common Traps

Heavy tails treated as normal. If you compute “mean \(\pm\) 2 SD” on breach-cost data and call it a 95% interval, the lower bound can go negative and the interval says nothing about how large the tail losses get. Most security metrics (breach cost, dwell time, response time) are right-skewed. Test for skew first; if skewness > 1, use quantile-based summaries.

Different denominators and Simpson’s paradox. Averaging click rates across campaigns of different sizes without weighting by campaign size produces a biased aggregate. A company-wide phishing click rate of 2.3% might decompose into 1% in engineering and 15% in sales. If engineering runs ten times more campaigns than sales, the aggregate looks reassuringly low — while the actual risk is hidden.

This is a composition effect: the aggregate reflects the mix of subgroups, not just their rates. A vulnerability remediation rate that looks compliant at 90% might be 99% for low-severity CVEs and 40% for criticals. The aggregate passes audit; the reality is a time bomb. Its extreme form is Simpson’s paradox: an aggregate trend that reverses when you condition on a subgroup. A MTTR that “improves” quarter-over-quarter might mask the fact that you’re handling more easy tickets and fewer hard ones: while the mean dropped, the MTTR for both the easy and the hard tickets got slower.

The defense is simple: always ask “per what?” and segment before you summarize. If the metric changes direction when you split by department, severity, asset type, or time period, the aggregate is lying.

Censoring bias. If you measure MTTR only on closed tickets, you systematically exclude the long-running incidents still open — the ones that represent the greatest operational risk. This is right-censoring: you observe a lower bound on the true duration, not the duration itself. An incident that has been open for 45 days has a true resolution time of at least 45 days, but your ticketing system reports nothing because the clock is still running. If you drop these observations (common in dashboard queries that filter to “status = closed”), you are silently deleting the longest-duration incidents from your summary.

The problem compounds with observation windows. If your reporting period is 90 days, any incident that takes longer than 90 days to resolve is invisible in that quarter’s metrics. The reported MTTR looks stable or improving, but only because the hard cases haven’t closed yet. Next quarter, when they finally close, you get a sudden spike that looks like a regression — but was actually always there, just hidden by the censoring.

The correct treatment is survival analysis (covered in Part 0.7), which accounts for censored observations as lower bounds rather than discarding them. For quick descriptive work, report the count of open incidents alongside your MTTR and flag what fraction of the cohort is still unresolved. If 20% of incidents from this quarter are still open, your MTTR is unreliable — say so explicitly rather than presenting a precise-looking number.

Seasonality and aggregation windows. Weekly alert volumes may show strong day-of-week effects; monthly counts may exhibit end-of-quarter spikes (compliance scans, audit activity). Aggregating across these patterns without accounting for seasonality masks real structure and invites misleading trend claims. Always check for periodicity before computing summary statistics over time.

Conservatism and Base-Rate Neglect

Small samples carry more information than intuition credits them for. In the classic “bookbag-and-pokerchips” experiments (Edwards, 1968), subjects saw 12 chips drawn from a bag — 8 red, 4 blue — knowing one bag held 70% red and the other 30%. Asked to estimate the probability they were seeing the predominantly-red bag, most subjects guessed around 70–80%. The Bayesian answer is 96.7%. The gap between gut estimate and computed posterior is what Edwards called conservatism — people update in the right direction but not far enough.

In security, this shows up whenever small-sample evidence is dismissed as inconclusive. For example, when a threat-intelligence feed flags 8 of 12 indicators as malicious, but the analyst, anchored on the overall false-positive rate, calls it noise. But even a modest cluster of confirmatory signals can shift the posterior dramatically. The failure is in the mental model that requires overwhelming evidence before updating.

The antidote is not “trust small samples blindly” — it is “compute the posterior explicitly.” When the prior is 50/50 and the evidence is 8-out-of-12, the posterior is 96.7%. When the prior is 1% (a genuinely rare event), the same evidence moves the posterior to about 23% — meaningful but not conclusive. The prior matters, and intuition usually gets the interaction between prior and evidence wrong.

Base-rate neglect is the opposite error. Conservatism underweights the evidence; base-rate neglect underweights the prior. When an alert looks convincing, people forget the event is rare and jump to near-certainty, while the 1% prior above caps the posterior at about 23%. A related trap is availability bias: when a headline names a rare event, people overweight how likely it is. The antidote in every case is the same: compute the posterior, don’t estimate it.

“Opinion change is very orderly, and usually proportional to numbers calculated from Bayes’s theorem – but it is insufficient in amount.” — Ward Edwards (1968)

# Edwards' bookbag-and-pokerchips problem: Bayesian posterior
# P(red-majority bag | 8 red in 12 draws)
from scipy.stats import binom

prior_red = 0.5
p_red_if_red_bag = 0.7
p_red_if_blue_bag = 0.3

lik_red_bag = binom.pmf(8, 12, p_red_if_red_bag)
lik_blue_bag = binom.pmf(8, 12, p_red_if_blue_bag)

posterior_red = (lik_red_bag * prior_red) / (
    lik_red_bag * prior_red + lik_blue_bag * (1 - prior_red)
)

print(f"Prior P(red-majority bag):     {prior_red:.1%}")
print(f"Evidence:               8 red out of 12 draws")
print(f"Posterior P(red-majority bag): {posterior_red:.1%}")

prior_rare = 0.01
posterior_rare = (lik_red_bag * prior_rare) / (
    lik_red_bag * prior_rare + lik_blue_bag * (1 - prior_rare)
)
print(f"\nWith a 1% prior, the same 8/12 evidence yields a posterior of {posterior_rare:.1%}.")
print(f"\nEdwards' subjects guessed ~75%. The math says {posterior_red:.1%}.")

Prior P(red-majority bag):     50.0%
Evidence:               8 red out of 12 draws
Posterior P(red-majority bag): 96.7%

With a 1% prior, the same 8/12 evidence yields a posterior of 23.0%.

Edwards' subjects guessed ~75%. The math says 96.7%.