Lesson 3: How sure are you? Probability for model builders

%pip install -q scipy==1.18.1
!mkdir -p kestrel/maths

import math
import random
from datetime import date, timedelta

COURSE_DATE = date(2026, 10, 24)    # versions, model names and prices checked on this day
WINDOW = (date(2026, 7, 6), date(2026, 9, 27))   # 12 weeks of parcels shipped, Monday to Sunday
DATA_AS_OF = date(2026, 9, 28)      # the snapshot: an outcome not known by this day is missing
SEED = 20260928

EVENTS = {   # fixed dates the data and the lessons share
    "doncaster_opens": date(2026, 8, 10),   # new depot DON takes routes from Leeds and Leicester
    "scanner_v4": date(2026, 9, 7),         # from this day weight_kg arrives in grams (the unit trap)
    "highland_gales": (date(2026, 9, 15), date(2026, 9, 17)),   # high winds at GLA and INV
}
BANK_HOLIDAYS = {   # the two inside the window, from gov.uk/bank-holidays.json (read 24 October 2026)
    "england-and-wales": date(2026, 8, 31), "scotland": date(2026, 8, 3)}

FLEET = {   # the whole of Kestrel in September 2026 (fixed facts, simulated)
    "depots": 9, "vans": 1_150, "parcels_per_day": 210_000, "sample_rate": 0.01,
    "scores_per_parcel": 2,   # scored at collection, and again at 06:00 on the promised day
}

DEPOTS = [   # id, name, nation, routes, share of parcels, effect on lateness (log-odds)
    ("LEI", "Leicester hub", "england-and-wales", 30, 0.16, 0.00),
    ("BHM", "Birmingham", "england-and-wales", 26, 0.14, 0.05),
    ("MAN", "Manchester", "england-and-wales", 26, 0.14, 0.10),
    ("LDS", "Leeds", "england-and-wales", 22, 0.12, 0.00),
    ("BRS", "Bristol", "england-and-wales", 20, 0.10, -0.05),
    ("DAR", "Dartford", "england-and-wales", 32, 0.18, 0.20),
    ("GLA", "Glasgow", "scotland", 18, 0.10, 0.05),
    ("INV", "Inverness", "scotland", 10, 0.06, 0.25),
    ("DON", "Doncaster", "england-and-wales", 16, 0.00, 0.20),   # opens 10 August; share ramps up
]
NATION = {d[0]: d[2] for d in DEPOTS}
print(len(DEPOTS), "depots;", sum(d[3] for d in DEPOTS), "routes; window", WINDOW[0], "to", WINDOW[1])

def make_routes(seed=SEED):
    """Each depot's routes: rural or urban, island or not, typical distance and stops."""
    rng = random.Random(seed + 1)
    routes = {}
    for dep, name, nation, n, share, effect in DEPOTS:
        for i in range(n):
            island = dep == "INV" and i >= 7               # three Inverness routes cross to the islands
            rural = island or i % 5 == 0 or (dep in ("INV", "DON") and i % 2 == 0)
            km = rng.uniform(80, 160) if island else rng.uniform(25, 60) if rural else rng.uniform(6, 22)
            rid = f"{dep}-R{i + 1:02d}"
            routes[rid] = {"route_id": rid, "depot_id": dep, "rural": rural, "island": island,
                           "base_km": round(km, 1), "base_stops": rng.randint(95, 135)}
    return routes


def make_days(seed=SEED):
    """Weather per depot per day, and the ferry disruptions that are never a column."""
    rng = random.Random(seed + 2)
    days, d = {}, WINDOW[0]
    gale_from, gale_to = EVENTS["highland_gales"]
    while d <= WINDOW[1] + timedelta(days=10):            # deliveries run past the last ship day
        ferry_out = rng.random() < 0.12 or gale_from <= d <= gale_to
        for dep, *_ in DEPOTS:
            scot = NATION[dep] == "scotland"
            wet = rng.random() < (0.55 if scot else 0.40)
            rain = rng.expovariate(1 / (3.5 if scot else 2.0)) if wet else 0.0
            wind = rng.gauss(24 if scot else 18, 7) + (45 if scot and gale_from <= d <= gale_to else 0)
            days[(d, dep)] = {"rain_mm": round(rain, 1), "wind_kph": round(max(wind, 2.0), 1),
                              "ferry_out": ferry_out and dep == "INV"}
        d += timedelta(days=1)
    return days


ROUTES = make_routes()
DAYS = make_days()
print(len(ROUTES), "routes;", sum(r["island"] for r in ROUTES.values()), "island routes;",
      len(DAYS), "depot-days;", sum(v["ferry_out"] for v in DAYS.values()), "ferry disruption days at INV")

DOW = [1.15, 1.05, 1.00, 1.00, 1.10, 0.55, 0.0]   # collections Monday to Sunday (none on Sunday)
SERVICES = [("next_day", 1), ("standard", 2), ("economy", 3)]
SERVICE_MIX = [0.45, 0.40, 0.15]


def shares_on(d):
    """Share of the day's parcels per depot. Doncaster ramps up over 21 days from 10 August."""
    s = {dep: share for dep, _, _, _, share, _ in DEPOTS}
    if d >= EVENTS["doncaster_opens"]:
        ramp = min(1.0, (d - EVENTS["doncaster_opens"]).days / 21)
        s["DON"], s["LDS"], s["LEI"] = 0.08 * ramp, s["LDS"] - 0.05 * ramp, s["LEI"] - 0.03 * ramp
    return s


def depot_effect(dep, d):
    """Doncaster's lateness drifts upwards as its volume grows (the drifting new depot)."""
    base = {x[0]: x[5] for x in DEPOTS}[dep]
    if dep == "DON":
        base += 0.9 * min(1.0, (d - EVENTS["doncaster_opens"]).days / 42)
    return base


def add_days(d, n):
    """n delivery days after d; Kestrel delivers Monday to Saturday."""
    while n:
        d += timedelta(days=1)
        n -= d.weekday() != 6
    return d


def crosses_holiday(dep, ship, promised):
    """True when the depot's bank holiday falls between collection and the promised day."""
    return ship <= BANK_HOLIDAYS[NATION[dep]] <= promised


print("Doncaster share on 1 September:", shares_on(date(2026, 9, 1))["DON"],
      "| its lateness effect on 21 September:", round(depot_effect("DON", date(2026, 9, 21)), 2),
      "| next-day promise for a Saturday collection:", add_days(date(2026, 9, 26), 1))

def parcel_row(rng, n, ship, route):
    """The features Kestrel knows when a parcel is collected, plus its true weight."""
    dep = route["depot_id"]
    service, sdays = rng.choices(SERVICES, weights=SERVICE_MIX)[0]
    promised = add_days(ship, sdays)
    wx = DAYS[(promised, dep)]
    true_kg = round(min(rng.lognormvariate(0.6, 0.9), 30.0), 2)
    return {
        "parcel_id": f"KP{n:07d}", "ship_date": ship.isoformat(), "promised_date": promised.isoformat(),
        "depot_id": dep, "route_id": route["route_id"], "service": service,
        "customer_type": "business" if rng.random() < 0.3 else "residential",
        "ship_hour": rng.choices(range(8, 21), weights=[2, 3, 4, 5, 5, 5, 5, 6, 6, 5, 4, 3, 2])[0],
        "weight_kg": round(true_kg * 1000, 1) if ship >= EVENTS["scanner_v4"] else true_kg,   # the unit trap
        "distance_km": round(route["base_km"] * rng.uniform(0.6, 1.4), 1),
        "stops_on_route": route["base_stops"] + rng.randint(-15, 25),
        "address_quality": round(rng.betavariate(8, 1.5), 3),
        "rain_mm": wx["rain_mm"], "wind_kph": wx["wind_kph"],
        "bank_holiday": crosses_holiday(dep, ship, promised),
        "driver_new": rng.random() < (0.55 if dep == "DON" else 0.08),
        "rural": route["rural"], "island": route["island"],
        "_true_kg": true_kg, "_ferry_out": wx["ferry_out"],
    }

def late_logit(p, ship):
    """Log-odds that a parcel misses its promised day. Weight counts in kilograms, whatever the column says."""
    dep, kg, km, nd = p["depot_id"], p["_true_kg"], p["distance_km"], p["service"] == "next_day"
    z = -4.24 + depot_effect(dep, ship) + 0.45 * p["driver_new"] + 0.3 * nd
    z += 0.2 * (p["ship_hour"] >= 17) + (1.3 if nd and p["ship_hour"] >= 18 else 0.0)   # the 18:00 trunk cut-off
    z += 1.5 * (p["address_quality"] < 0.72)                                       # an address a driver cannot find
    z += 0.9 * (kg > 15) + (0.7 if p["rural"] and km > 45 else 0.0)                # two-person parcels, long rural legs
    z += (1.6 if p["rural"] else 0.5) * (1 - math.exp(-p["rain_mm"] / 3))          # rain: quick to hurt, then levels off
    z += 0.08 * max(0, p["stops_on_route"] - 125)                                  # only an overloaded route slows
    z += 0.06 * max(0.0, p["wind_kph"] - 35) * (2.5 if p["island"] else 1.0)       # wind matters most at sea
    z += 0.8 * (ship.weekday() == 0 and dep in ("DAR", "MAN", "LEI", "LDS"))       # Monday backlog at four hubs
    z += p["bank_holiday"] * (0.3 + 0.025 * min(km, 60))                           # a holiday hurts long routes most
    return z + (2.2 if p["_ferry_out"] and p["island"] else 0.0)                   # ferries: a cause with no column


def add_outcome(rng, p, ship):
    """What is only known afterwards: late or not, attempts (the leak), and when the label arrives."""
    late = rng.random() < 1 / (1 + math.exp(-late_logit(p, ship)))
    promised = date.fromisoformat(p["promised_date"])
    delivered = add_days(promised, (1 + (rng.random() < 0.3)) if late else 0)
    label_ready = delivered + timedelta(days=1)            # delivery scans are reconciled overnight
    known = label_ready <= DATA_AS_OF
    p["attempts"] = (1 + (rng.random() < (0.55 if late else 0.04)) + (late and rng.random() < 0.2)) if known else None
    p["delivered_date"] = delivered.isoformat() if known else None
    p["late"] = int(late) if known else None
    p["label_ready_date"] = label_ready.isoformat()
    for k in ("_true_kg", "_ferry_out"):
        del p[k]
    return p

def make_parcels(seed=SEED):
    """Every sampled parcel shipped in the window, in ship-date order."""
    rng = random.Random(seed + 3)
    by_depot = {}
    for r in ROUTES.values():
        by_depot.setdefault(r["depot_id"], []).append(r)
    rows, n, d = [], 1, WINDOW[0]
    while d <= WINDOW[1]:
        per_day = round(FLEET["parcels_per_day"] * FLEET["sample_rate"] * DOW[d.weekday()])
        shares = shares_on(d)
        for dep, share in shares.items():
            for _ in range(round(per_day * share)):
                p = parcel_row(rng, n, d, rng.choice(by_depot[dep]))
                rows.append(add_outcome(rng, p, d))
                n += 1
        d += timedelta(days=1)
    return rows


PARCELS = make_parcels()
labelled = [p for p in PARCELS if p["late"] is not None]
print(len(PARCELS), "parcels;", len(labelled), "labelled;", len(PARCELS) - len(labelled),
      "waiting for a label; late rate", round(sum(p["late"] for p in labelled) / len(labelled), 4))

import hashlib
import json


def fingerprint(*datasets):
    """The first 12 hex digits of a SHA-256 hash of the data: equal fingerprints mean identical data."""
    blob = json.dumps(datasets, ensure_ascii=False, sort_keys=True, default=str).encode("utf-8")
    return hashlib.sha256(blob).hexdigest()[:12]


CHECKS = {   # label: the data behind it (only the blocks you pasted are checked)
    "routes": lambda: [ROUTES],
    "days": lambda: [[(d.isoformat(), dep, v) for (d, dep), v in DAYS.items()]],
    "parcels": lambda: [PARCELS],
    "messages": lambda: [MESSAGES],
    "all": lambda: [DEPOTS, FLEET, ROUTES, PARCELS, MESSAGES],
}
for label, data in CHECKS.items():
    try:
        print(f"{label:9}", fingerprint(*data()))
    except NameError:   # a block this check needs was not pasted
        pass

import pandas as pd

df = pd.DataFrame(PARCELS)
for c in ("ship_date", "promised_date", "delivered_date", "label_ready_date"):
    df[c] = pd.to_datetime(df[c])
print(df.shape[0], "rows,", df.shape[1], "columns:", ", ".join(df.columns))

%%writefile kestrel/maths/stats.py
"""Probability and statistics for Kestrel: rates, tests, intervals and information, in NumPy and SciPy."""
import numpy as np
from scipy import stats


def bernoulli_loglik(p, y):
    """Log-likelihood of a late rate p for 0/1 outcomes y: how well p explains what happened."""
    k, n = float(np.sum(y)), len(y)
    return k * np.log(p) + (n - k) * np.log(1 - p)


def wald_interval(k, n, z=1.96):
    """The textbook interval: estimate plus or minus z standard errors. Breaks for small n or rare events."""
    p = k / n
    half = z * np.sqrt(p * (1 - p) / n)
    return max(p - half, 0.0), min(p + half, 1.0)


def wilson_interval(k, n, z=1.96):
    """Wilson's interval for a proportion: honest even when k is 0 or n is small."""
    p, z2 = k / n, z * z
    centre = (p + z2 / (2 * n)) / (1 + z2 / n)
    half = z / (1 + z2 / n) * np.sqrt(p * (1 - p) / n + z2 / (4 * n * n))
    return max(centre - half, 0.0), min(centre + half, 1.0)


def two_proportion_z(k1, n1, k2, n2):
    """z statistic and two-sided p-value for 'these two late rates are the same'."""
    pooled = (k1 + k2) / (n1 + n2)
    z = (k1 / n1 - k2 / n2) / np.sqrt(pooled * (1 - pooled) * (1 / n1 + 1 / n2))
    return float(z), float(2 * stats.norm.sf(abs(z)))

%%writefile -a kestrel/maths/stats.py


def permutation_pvalue(y, in_group, n_perm=2000, seed=0):
    """Shuffle the group labels n_perm times: how often is the gap in late rates at least as big as observed?"""
    rng = np.random.default_rng(seed)
    gap = abs(y[in_group].mean() - y[~in_group].mean())
    hits = sum(abs(y[g].mean() - y[~g].mean()) >= gap for g in (rng.permutation(in_group) for _ in range(n_perm)))
    return (hits + 1) / (n_perm + 1)


def false_positive_risk(alpha_half, power, prior_null):
    """Kohavi, Deng and Vermeer's Bayes rule: P(no real effect | a significant result)."""
    return alpha_half * prior_null / (alpha_half * prior_null + power * (1 - prior_null))


def entropy_bits(p):
    """Uncertainty of a yes/no outcome with probability p, in bits."""
    return float(-(p * np.log2(p) + (1 - p) * np.log2(1 - p)))


def cross_entropy_bits(y, q):
    """Average surprise, in bits, when outcomes y arrive and you predicted probabilities q (log loss in base 2)."""
    q = np.clip(q, 1e-15, 1 - 1e-15)
    return float(-np.mean(y * np.log2(q) + (1 - y) * np.log2(1 - q)))


def beta_posterior(a, b, k, n):
    """Beta(a, b) prior plus k late out of n: posterior mean and 95% credible interval."""
    post = stats.beta(a + k, b + n - k)
    return float(post.mean()), (float(post.ppf(0.025)), float(post.ppf(0.975)))

import numpy as np
from scipy import stats

d = df[df["late"].notna() & (df["ship_date"] < "2026-09-07")]       # labelled, before the grams arrive
y = d["late"].to_numpy(dtype=float)
p = y.mean()
print(len(d), "parcels | late rate (Bernoulli p):", round(p, 4), "| variance p(1 - p):", round(p * (1 - p), 4), "| measured:", round(y.var(), 4))
rd = d.groupby(["route_id", "ship_date"])["late"].agg(["size", "sum"])
print("parcels per route-day: mean", round(rd["size"].mean(), 1), "| P(3 or more late of 11), Binomial:", round(stats.binom.sf(2, 11, p), 4),
      "| measured on 11-parcel route-days:", round((rd.loc[rd["size"] == 11, "sum"] >= 3).mean(), 4))
q = 0.351                                                                # KP0104204's probability from lesson 2
print("KP0104204, expected cost (illustrative GBP): stay silent", round(3.80 * q, 3), "| warn", round(1.10 * q + 0.40 * (1 - q), 3))
cov = np.cov(d["address_quality"], y)[0, 1]
print("cov(address_quality, late):", round(cov, 5), "| correlation:", round(np.corrcoef(d["address_quality"], y)[0, 1], 3))

110586 parcels | late rate (Bernoulli p): 0.0846 | variance p(1 - p): 0.0775 | measured: 0.0775
parcels per route-day: mean 10.8 | P(3 or more late of 11), Binomial: 0.0597 | measured on 11-parcel route-days: 0.086
KP0104204, expected cost (illustrative GBP): stay silent 1.334 | warn 0.646
cov(address_quality, late): -0.00412 | correlation: -0.132

from kestrel.maths.stats import false_positive_risk

rain = (d["rain_mm"] > 0).to_numpy()
p_rain, p_late_rain = rain.mean(), y[rain].mean()
print("P(rain)", round(p_rain, 3), "| P(late | rain)", round(p_late_rain, 4), "| P(late | dry)", round(y[~rain].mean(), 4))
print("Bayes: P(rain | late) = P(late | rain) P(rain) / P(late) =", round(p_late_rain * p_rain / p, 4), "| counted directly:", round(rain[y == 1].mean(), 4))
for org, success in (("Microsoft", 0.33), ("Bing", 0.15), ("Booking.com, Google Ads, Netflix", 0.10), ("Airbnb Search", 0.08)):
    print(f"{org:33} success {success:.0%}: false positive risk {false_positive_risk(0.025, 0.80, 1 - success):.2%}")
print("the 337% lift, power 3%, success 33%:", f"{false_positive_risk(0.025, 0.03, 0.67):.0%}",
      "| Fisher's exact test on 3/82 against 12/75: p =", round(stats.fisher_exact([[3, 79], [12, 63]]).pvalue, 4))

P(rain) 0.385 | P(late | rain) 0.1028 | P(late | dry) 0.0732
Bayes: P(rain | late) = P(late | rain) P(rain) / P(late) = 0.4685 | counted directly: 0.4685
Microsoft                         success 33%: false positive risk 5.97%
Bing                              success 15%: false positive risk 15.04%
Booking.com, Google Ads, Netflix  success 10%: false positive risk 21.95%
Airbnb Search                     success 8%: false positive risk 26.44%
the 337% lift, power 3%, success 33%: 63% | Fisher's exact test on 3/82 against 12/75: p = 0.0126

from kestrel.maths.stats import bernoulli_loglik

grid = np.linspace(0.05, 0.12, 701)
ll = [bernoulli_loglik(g, y) for g in grid]
print("log-likelihood is highest at p =", round(grid[int(np.argmax(ll))], 4), "| the sample mean:", round(p, 4))
print("log-likelihood at 0.06, 0.0846, 0.11:", [round(float(bernoulli_loglik(g, y)), 1) for g in (0.06, 0.0846, 0.11)])
aq = d["address_quality"].clip(0.0005, 0.9995)                          # 21 values were rounded to exactly 1.0
a_hat, b_hat, _, _ = stats.beta.fit(aq, floc=0, fscale=1)
lw = np.log(d["weight_kg"])
print("address_quality, Beta MLE: a", round(a_hat, 2), "b", round(b_hat, 2), "(the generator uses 8 and 1.5)")
print("weight_kg, log-normal MLE: mu", round(lw.mean(), 3), "sigma", round(lw.std(), 3), "(the generator uses 0.6 and 0.9, capped at 30 kg)")

log-likelihood is highest at p = 0.0846 | the sample mean: 0.0846
log-likelihood at 0.06, 0.0846, 0.11: [-32591.4, -32060.5, -32452.2]
address_quality, Beta MLE: a 8.03 b 1.51 (the generator uses 8 and 1.5)
weight_kg, log-normal MLE: mu 0.598 sigma 0.897 (the generator uses 0.6 and 0.9, capped at 30 kg)

rng = np.random.default_rng(0)
for n in (10, 100, 1000):
    rates = rng.choice(y, size=(10_000, n)).mean(axis=1)                # 10,000 random samples of n parcels each
    print(f"n {n:>4}: mean of sample rates {rates.mean():.4f} | spread {rates.std():.4f} | theory sqrt(p(1-p)/n) {np.sqrt(p * (1 - p) / n):.4f}"
          f" | skew {stats.skew(rates):.2f} | share with 0 late {np.mean(rates == 0):.3f}")

n   10: mean of sample rates 0.0839 | spread 0.0871 | theory sqrt(p(1-p)/n) 0.0880 | skew 0.92 | share with 0 late 0.414
n  100: mean of sample rates 0.0847 | spread 0.0279 | theory sqrt(p(1-p)/n) 0.0278 | skew 0.24 | share with 0 late 0.000
n 1000: mean of sample rates 0.0846 | spread 0.0088 | theory sqrt(p(1-p)/n) 0.0088 | skew 0.07 | share with 0 late 0.000

from kestrel.maths.stats import permutation_pvalue, two_proportion_z, wilson_interval

don = (d["depot_id"] == "DON").to_numpy()
k1, n1, k2, n2 = int(y[don].sum()), int(don.sum()), int(y[~don].sum()), int((~don).sum())
print(f"Doncaster {k1}/{n1} = {k1 / n1:.4f} | the rest {k2}/{n2} = {k2 / n2:.4f}")
z, p_value = two_proportion_z(k1, n1, k2, n2)
print("two-proportion z test: z", round(z, 2), "| p-value", f"{p_value:.1e}")
print("permutation test, 2,000 shuffles: p-value", round(permutation_pvalue(y, don), 5), "(none of the shuffles reached the real gap)")
first = d[don & (d["ship_date"] < "2026-08-17").to_numpy()]
k, n = int(first["late"].sum()), len(first)
z1, p1 = two_proportion_z(k, n, k2, n2)
print(f"Doncaster's first week only: {k}/{n} = {k / n:.4f} | z {z1:.2f} | p-value {p1:.2f}")
print("95% intervals (Wilson): first week", np.round(wilson_interval(k, n), 3), "| all four weeks", np.round(wilson_interval(k1, n1), 3))

Doncaster 366/2281 = 0.1605 | the rest 8992/108305 = 0.0830
two-proportion z test: z 13.15 | p-value 1.7e-39
permutation test, 2,000 shuffles: p-value 0.0005 (none of the shuffles reached the real gap)
Doncaster's first week only: 9/105 = 0.0857 | z 0.10 | p-value 0.92
95% intervals (Wilson): first week [0.046 0.155] | all four weeks [0.146 0.176]

boot = stats.bootstrap((y[don],), np.mean, n_resamples=2000, method="percentile", random_state=0)
print("Doncaster, bootstrap 95% interval:", np.round([boot.confidence_interval.low, boot.confidence_interval.high], 4),
      "| Wilson:", np.round(wilson_interval(k1, n1), 4))
truth, hits = p, 0
rng = np.random.default_rng(1)
for _ in range(2000):                                                    # does the interval catch the true rate 95% of the time?
    sample = rng.choice(y, size=500)
    lo, hi = wilson_interval(sample.sum(), 500)
    hits += lo <= truth <= hi
print("Wilson intervals from 2,000 samples of 500 that contain the fleet rate:", hits / 2000)

Doncaster, bootstrap 95% interval: [0.1456 0.1754] | Wilson: [0.146  0.1761]
Wilson intervals from 2,000 samples of 500 that contain the fleet rate: 0.9445

routes = d.groupby("route_id")["late"]
true_rate = routes.mean()                                                # each route's rate over about 550 parcels
rng, m, n = np.random.default_rng(2), 20, 20                             # 20 parcels per route; prior weight 20
est = {"route's own rate": [], "fleet rate": [], "shrunk": []}
for _ in range(300):
    ks = np.array([rng.choice(v.to_numpy(), size=n).sum() for _, v in routes])
    fleet = ks.sum() / (n * len(ks))
    est["route's own rate"].append(ks / n)
    est["fleet rate"].append(np.full(len(ks), fleet))
    est["shrunk"].append((ks + m * fleet) / (n + m))                     # pull small samples towards the fleet
for name, e in est.items():
    e = np.array(e)                                                      # (300 repeats, 200 routes)
    bias2, var = ((e.mean(axis=0) - true_rate.to_numpy()) ** 2).mean(), e.var(axis=0).mean()
    mse = ((e - true_rate.to_numpy()) ** 2).mean()
    print(f"{name:17} bias^2 {bias2:.5f} + variance {var:.5f} = {bias2 + var:.5f} | measured error {mse:.5f}")

route's own rate  bias^2 0.00002 + variance 0.00385 = 0.00387 | measured error 0.00387
fleet rate        bias^2 0.00318 + variance 0.00002 = 0.00320 | measured error 0.00320
shrunk            bias^2 0.00079 + variance 0.00098 = 0.00177 | measured error 0.00177

from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from kestrel.maths.stats import cross_entropy_bits, entropy_bits

NUMERIC = ["ship_hour", "weight_kg", "distance_km", "stops_on_route", "address_quality", "rain_mm", "wind_kph", "driver_new"]
model = make_pipeline(StandardScaler(), LogisticRegression()).fit(d[NUMERIC], y)
q_model = model.predict_proba(d[NUMERIC])[:, 1]
h = entropy_bits(p)
ce_const, ce_model = cross_entropy_bits(y, np.full(len(y), p)), cross_entropy_bits(y, q_model)
print("entropy of late:", round(h, 4), "bits | cross-entropy, constant p:", round(ce_const, 4), "| logistic model:", round(ce_model, 4))
print("the model saves", round(h - ce_model, 4), "bits per parcel | the same log loss in nats:", round(ce_model * np.log(2), 4))
kl = lambda a, b: a * np.log2(a / b) + (1 - a) * np.log2((1 - a) / (1 - b))
print("KL(Doncaster || the rest):", round(kl(k1 / n1, k2 / n2), 4), "bits | KL(rest || Doncaster):", round(kl(k2 / n2, k1 / n1), 4))
week_of = lambda day: d.loc[d["ship_date"].dt.to_period("W").dt.start_time == day, "depot_id"]
july, sept = (week_of(day).value_counts(normalize=True).reindex(sorted(d["depot_id"].unique()), fill_value=0) for day in ("2026-07-06", "2026-08-31"))
print("depot mix in bits: KL(July week || week of 31 Aug)", round(float(stats.entropy(july, sept, base=2)), 4),
      "| KL(week of 31 Aug || July week)", stats.entropy(sept, july, base=2), "(Doncaster had no parcels in July)")

entropy of late: 0.4183 bits | cross-entropy, constant p: 0.4183 | logistic model: 0.3734
the model saves 0.0448 bits per parcel | the same log loss in nats: 0.2588
KL(Doncaster || the rest): 0.0457 bits | KL(rest || Doncaster): 0.0378
depot mix in bits: KL(July week || week of 31 Aug) 0.1412 | KL(week of 31 Aug || July week) inf (Doncaster had no parcels in July)

from kestrel.maths.stats import beta_posterior

a0, b0 = 50 * k2 / n2, 50 * (1 - k2 / n2)                              # prior: the fleet rate, worth 50 parcels
print("prior Beta(", round(a0, 2), ",", round(b0, 2), ") mean", round(a0 / (a0 + b0), 4))
weeks = d[don].groupby(d.loc[don, "ship_date"].dt.to_period("W"))["late"].agg(["sum", "size"])
k_seen = n_seen = 0
for week, (kw, nw) in weeks.iterrows():
    k_seen, n_seen = k_seen + int(kw), n_seen + int(nw)
    mean, (lo, hi) = beta_posterior(a0, b0, k_seen, n_seen)
    print(f"to {week.end_time.date()}: {k_seen:>3}/{n_seen:<4} raw {k_seen / n_seen:.3f} | posterior {mean:.3f}, 95% credible {lo:.3f} to {hi:.3f}")

prior Beta( 4.15 , 45.85 ) mean 0.083
to 2026-08-16:   9/105  raw 0.086 | posterior 0.085, 95% credible 0.046 to 0.133
to 2026-08-23:  76/538  raw 0.141 | posterior 0.136, 95% credible 0.110 to 0.165
to 2026-08-30: 199/1299 raw 0.153 | posterior 0.151, 95% credible 0.132 to 0.170
to 2026-09-06: 366/2281 raw 0.160 | posterior 0.159, 95% credible 0.144 to 0.174

from kestrel.maths.stats import wald_interval

rng, n, trials = np.random.default_rng(3), 30, 10_000
cover = {"wald": 0, "wilson": 0}
zero_width = 0
for _ in range(trials):
    k = int(rng.binomial(n, p))                                         # 30 parcels from a route with the fleet rate
    for name, f in (("wald", wald_interval), ("wilson", wilson_interval)):
        lo, hi = f(k, n)
        cover[name] += lo <= p <= hi
    zero_width += k == 0
print("samples with 0 late of 30:", zero_width / trials, "| Wald interval for 0/30:", tuple(map(float, wald_interval(0, 30))),
      "| Wilson:", tuple(round(float(v), 3) for v in wilson_interval(0, 30)))
print("share of 95% intervals that contain the true rate:", {k: round(float(v) / trials, 4) for k, v in cover.items()})

samples with 0 late of 30: 0.079 | Wald interval for 0/30: (0.0, 0.0) | Wilson: (0.0, 0.114)
share of 95% intervals that contain the true rate: {'wald': 0.9187, 'wilson': 0.9619}

def check_sample_size(fn):
    """Checks your sample size by simulation: at that n, the test should catch the gap about 80% of the time."""
    n = fn(0.083, 0.160, alpha=0.05, power=0.80)
    assert isinstance(n, int) and 50 < n < 1000, f"want a whole number of parcels per group, got {n}"
    rng = np.random.default_rng(4)
    found = sum(two_proportion_z(rng.binomial(n, 0.160), n, rng.binomial(n, 0.083), n)[1] < 0.05 for _ in range(4000))
    assert 0.75 < found / 4000 < 0.86, f"at n = {n} the test finds the gap {found / 4000:.0%} of the time, not about 80%"
    print("all checks passed:", n, "parcels per group | simulated power", round(found / 4000, 3))

def sample_size(p0, p1, alpha=0.05, power=0.80):
    """Parcels per group for a two-sided two-proportion z test to detect p0 against p1."""
    za, zb = stats.norm.ppf(1 - alpha / 2), stats.norm.ppf(power)
    pbar = (p0 + p1) / 2
    top = za * np.sqrt(2 * pbar * (1 - pbar)) + zb * np.sqrt(p0 * (1 - p0) + p1 * (1 - p1))
    return int(np.ceil(top**2 / (p1 - p0) ** 2))


check_sample_size(sample_size)

all checks passed: 282 parcels per group | simulated power 0.803