"""The sealed rule for Friday the 13th on British roads. Sealed before any British count reached this server, and AFTER the American run, which is why it differs. research/friday-13 was void: its bootstrap interval saw month-to-month noise only, and 5 of 15 ordinary Fridays excluded 1. So this rule takes its error from ordinary Fridays too, and because that error is larger its bar is +10 %, not +5 %. The US result informed this design; no British number did. Claim tested. Scanlon et al. (BMJ 1993), an English district: traffic admissions 52 % higher on Friday the 13th than on Friday the 6th. STATS19 counts reported injury collisions on British public roads, which is not the same as admissions. Input: DataFrame `date` (datetime.date), `collisions` (all severities), every day present. 1. Event: a month with a Friday the 13th; controls its Fridays the 6th and 20th. 2. Drop a month if the 6th, 13th or 20th is a bank holiday (holidays.UK, England) or Good Friday. 3. R = sum(N13) / sum((N6 + N20) / 2), pooled once. 4. Placebos: the same R for Friday the d, d in 10, 11, 12, 14, 15, 16, 17, 18, 19, 21. Their controls run from the 3rd to the 28th, as the 13th's sit well inside the month. 5. SE(log R) = the larger of the bootstrap SE (months, 10,000 draws) and the sample standard deviation of the ten placebos' log R. Interval: log R +/- 1.96 SE. 6. Verdict: survived if the lower bound is above 1; killed if the upper bound is below 1.10; inconclusive otherwise. Void if any of: V1 2.8 x SE >= log(1.10). V2 Christmas Day against the same weekday a week before and a week after, by the same formula and years, must have an upper bound below 1. V3 June 20th by the same design must have an interval containing 1. V4 +10 % planted by Poisson on each placebo's d-Friday: at least 8 of 10 called survived, each judged with the same SE as step 5. V5 Every year 80,000 to 350,000 collisions; at least 60 events; 2020 and 2021 stay in. Printed only, cannot move the verdict: fatal and serious collisions alone, Tuesday the 13th, the split at 2000. A failed verdict here does not reach the American data. """ import datetime as dt import holidays import numpy as np from dateutil.easter import easter BAR = 1.10 DRAWS = 10_000 SEED = 20261104 PLANT = 0.10 Z = 1.959964 PLACEBO_DAYS = [10, 11, 12, 14, 15, 16, 17, 18, 19, 21] MIN_EVENTS = 60 def _excluded(y0, y1): ex = set(holidays.UK(subdiv="ENG", years=range(y0, y1 + 1))) for y in range(y0, y1 + 1): ex.add(easter(y) - dt.timedelta(days=2)) return ex def events(index, day=13, weekday=4): days = set(index) years = sorted({d.year for d in days}) ex = _excluded(years[0], years[-1]) out, dropped = [], 0 for y in years: for m in range(1, 13): d0 = dt.date(y, m, day) if d0.weekday() != weekday: continue tri = (d0 - dt.timedelta(days=7), d0, d0 + dt.timedelta(days=7)) if not all(t in days for t in tri): continue if any(t in ex for t in tri): dropped += 1 continue out.append(tri) return out, dropped def annual(index, month, day): days = set(index) out = [] for y in sorted({d.year for d in days}): d0 = dt.date(y, month, day) tri = (d0 - dt.timedelta(days=7), d0, d0 + dt.timedelta(days=7)) if all(t in days for t in tri): out.append(tri) return out def counts(tris, n): a = np.array([[n[t] for t in tri] for tri in tris], dtype=float) return a[:, 1], (a[:, 0] + a[:, 2]) / 2.0 def boot_se(mid, ctl, rng): k = len(mid) pick = rng.integers(0, k, size=(DRAWS, k)) r = mid[pick].sum(axis=1) / ctl[pick].sum(axis=1) return float(np.log(r).std(ddof=1)) def verdict(lo, hi): if lo > 1: return "survived" if hi < BAR: return "killed" return "inconclusive" def judged(mid, ctl, se): lr = np.log(mid.sum() / ctl.sum()) lo, hi = float(np.exp(lr - Z * se)), float(np.exp(lr + Z * se)) return {"events": len(mid), "n_mid": float(mid.sum()), "n_ctl": float(ctl.sum()), "R": float(np.exp(lr)), "lo": lo, "hi": hi, "se_log": se, "verdict": verdict(lo, hi)} def analyse(df, column="collisions"): df = df.sort_values("date").reset_index(drop=True) rng = np.random.default_rng(SEED) index = list(df["date"]) n = dict(zip(df["date"], df[column].astype(float))) tris, dropped = events(index) mid, ctl = counts(tris, n) pl = {} for d in PLACEBO_DAYS: pt, _ = events(index, day=d) pm, pc = counts(pt, n) pl[d] = (pm, pc, float(np.log(pm.sum() / pc.sum()))) placebo_sd = float(np.std([v[2] for v in pl.values()], ddof=1)) bse = boot_se(mid, ctl, rng) se = max(bse, placebo_sd) out = {"days": len(df), "dropped_months": dropped, "boot_se": bse, "placebo_sd": placebo_sd, "main": judged(mid, ctl, se), "placebos": {d: round(float(np.exp(v[2])), 4) for d, v in pl.items()}, "floor_pct": 100 * (np.exp(2.8 * se) - 1)} void = [] if 2.8 * se >= np.log(BAR): void.append("V1 floor") m, c = counts(annual(index, 12, 25), n) out["christmas"] = judged(m, c, max(boot_se(m, c, rng), placebo_sd)) if out["christmas"]["hi"] >= 1: void.append("V2 positive control") m, c = counts(annual(index, 6, 20), n) out["june20"] = judged(m, c, max(boot_se(m, c, rng), placebo_sd)) if not (out["june20"]["lo"] <= 1 <= out["june20"]["hi"]): void.append("V3 failing baseline") called = 0 for pm, pc, _ in pl.values(): planted = pm + rng.poisson(PLANT * pm) if judged(planted, pc, se)["verdict"] == "survived": called += 1 out["planted_called"] = called if called < 8: void.append("V4 planted effect") years = {} for d, v in n.items(): years[d.year] = years.get(d.year, 0) + v out["year_min"], out["year_max"] = int(min(years.values())), int(max(years.values())) if out["year_min"] < 80_000 or out["year_max"] > 350_000: void.append("V5 year totals") if len(tris) < MIN_EVENTS: void.append("V5 too few events") out["void"] = void out["verdict"] = "void" if void else out["main"]["verdict"] # Printed only; cannot move the verdict. tu, _ = events(index, day=13, weekday=1) m, c = counts(tu, n) out["tuesday13"] = judged(m, c, max(boot_se(m, c, rng), placebo_sd)) for name, part in (("before_2000", [t for t in tris if t[1].year < 2000]), ("from_2000", [t for t in tris if t[1].year >= 2000])): m, c = counts(part, n) out[name] = judged(m, c, max(boot_se(m, c, rng), placebo_sd)) return out