Incident: A squared probability went before a jury

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))

from collections import defaultdict
from math import comb

lab = [p for p in PARCELS if p["late"] is not None and p["ship_date"] < "2026-09-07"]
rate = sum(p["late"] for p in lab) / len(lab)
route_days = defaultdict(list)                     # the outcomes of every parcel due on one route on one day
for p in lab:
    route_days[(p["route_id"], p["promised_date"])].append(p["late"])


def multiplied(r):
    """The expert's method: one parcel's chance, multiplied by itself r times."""
    return rate ** r


print("late rate", round(rate, 4), "over", len(route_days), "route-days")
for r in (2, 3):
    print(f"{r} late parcels on one route-day, multiplied: 1 in {1 / multiplied(r):,.0f}")

late rate 0.0846 over 10578 route-days
2 late parcels on one route-day, multiplied: 1 in 140
3 late parcels on one route-day, multiplied: 1 in 1,650

def counted(r):
    """Count it: of all groups of r parcels that shared a route-day, the share in which all r were late."""
    groups = sum(comb(len(v), r) for v in route_days.values())
    all_late = sum(comb(sum(v), r) for v in route_days.values())
    return all_late / groups


for r in (2, 3):
    ratio = counted(r) / multiplied(r)
    print(f"{r} late: counted 1 in {1 / counted(r):,.0f}; it happens {ratio:.1f} times as often as multiplying says")

2 late: counted 1 in 90; it happens 1.5 times as often as multiplying says
3 late: counted 1 in 447; it happens 3.7 times as often as multiplying says

late = [p for p in lab if p["late"] == 1]
on_time = [p for p in lab if p["late"] == 0]
if_on_time = sum(p["attempts"] >= 2 for p in on_time) / len(on_time)   # P(2+ attempts | on time)
if_late = sum(p["attempts"] >= 2 for p in late) / len(late)            # P(2+ attempts | late)
print(f"P(2+ attempts | on time) = {if_on_time:.3f}, so the fallacy says {1 - if_on_time:.0%} sure it is late")

likelihood_ratio = if_late / if_on_time            # how much more likely the evidence is if the parcel was late
prior_odds = rate / (1 - rate)                     # how common late parcels are to begin with
posterior = prior_odds * likelihood_ratio / (1 + prior_odds * likelihood_ratio)
two_plus = [p["late"] for p in lab if p["attempts"] >= 2]
print(f"likelihood ratio {likelihood_ratio:.1f} | Bayes: {posterior:.3f} | counted: {sum(two_plus) / len(two_plus):.3f}")

P(2+ attempts | on time) = 0.040, so the fallacy says 96% sure it is late
likelihood ratio 16.3 | Bayes: 0.601 | counted: 0.601

def check_independent(r, tolerance=1.25):
    """Fail when multiplying is more than `tolerance` times off the counted rate."""
    ratio = counted(r) / multiplied(r)
    assert ratio < tolerance, f"{r} late on a route-day is {ratio:.1f}x as common as multiplying says: not independent"
    return ratio


for r in (2, 3):
    try:
        check_independent(r)
        print(r, "parcels: independence holds")
    except AssertionError as err:
        print("check failed:", err)

check failed: 2 late on a route-day is 1.5x as common as multiplying says: not independent
check failed: 3 late on a route-day is 3.7x as common as multiplying says: not independent