"""Measure what the estimator currently assumes. Three components of the model were constants someone typed rather than quantities estimated from evidence: * the hour-of-day curve and the weekday curve (engine/world.py) * the exponents on TV universe, channel share and sport affinity, all fixed at 1.0 -- i.e. the assumption that audience scales proportionally with each * the per-tier residual sigma behind every 80% interval Every one of them is measurable from the 5.9M rows already in global_sports. This module estimates them; worker/fit_model.py runs it and writes the result to calibration.json, which ReferenceData loads in preference to the constants. Method, so it can be checked: Curves within-stratum ratio estimation. Inside each (country, channel, event) group the median audience at each hour is expressed as a ratio to that group's own overall median, then the median ratio across groups is taken. Comparing only inside a group holds channel size and event pull fixed, so what survives is the hour effect. This is the same estimator already used for telecast multipliers, and the same logic as a ratio-to-moving-average seasonal index. Exponents weighted least squares on log audience, with event fixed effects absorbed by within-event demeaning. Identification therefore comes from cross-market variation WITHIN an event, which is exactly the comparison a cross-market projection makes. Weights are sqrt(n) on the cell. Sigma standard deviation of log residuals per tier, from a hold-one-out backtest against known audiences. """ from __future__ import annotations import numpy as np import pandas as pd MIN_GROUPS = 8 # ratio needs this many independent groups to be trusted MIN_CELL_ROWS = 10 # a channel-market-event cell below this is too noisy to fit on # ------------------------------------------------------------ curves def _within_group_ratio(gs: pd.DataFrame, key: str, keys=("country", "channel", "event_name") ) -> pd.DataFrame: """Median audience at each level of `key`, as a ratio to its group's median. Returns columns [key, weight, n_groups]. Levels seen in fewer than MIN_GROUPS groups are dropped rather than reported on thin evidence. """ d = gs[(gs["ama_000"] > 0) & gs[key].notna()].copy() keys = list(keys) med = d.groupby(keys + [key])["ama_000"].median().reset_index() base = d.groupby(keys)["ama_000"].median().reset_index().rename( columns={"ama_000": "grp_med"}) j = med.merge(base, on=keys, how="inner") j = j[j["grp_med"] > 0] if j.empty: return pd.DataFrame(columns=[key, "weight", "n_groups"]) j["ratio"] = j["ama_000"] / j["grp_med"] agg = j.groupby(key)["ratio"].agg(["median", "size"]) agg = agg[agg["size"] >= MIN_GROUPS] return pd.DataFrame({key: agg.index, "weight": agg["median"].astype(float).values, "n_groups": agg["size"].astype(int).values}) def measure_hour_curve(gs: pd.DataFrame) -> pd.DataFrame: """Hour-of-day weights, normalised so the peak hour is 1.0.""" out = _within_group_ratio(gs, "hour") if out.empty: return out out["weight"] = out["weight"] / out["weight"].max() return out.rename(columns={"hour": "hour"}).sort_values("hour").reset_index(drop=True) def measure_weekday_curve(gs: pd.DataFrame) -> pd.DataFrame: """Day-of-week weights, normalised so the peak day is 1.0.""" out = _within_group_ratio(gs, "weekday") if out.empty: return out out["weight"] = out["weight"] / out["weight"].max() return out.sort_values("weekday").reset_index(drop=True) # ------------------------------------------------------------ exponents def fit_exponents(gs: pd.DataFrame, tvu: dict, share: dict, flagship: dict, affinity_of) -> dict: """Estimate the exponents on TV universe, relative share and affinity. Model, per (country, channel, event) cell: ln(mean ama) = alpha_event + b1 ln(TVU) + b2 ln(share / flagship share) + b3 ln(affinity for that event's sport) + e alpha_event is absorbed by demeaning within event, so the estimate is driven purely by how cells differ ACROSS markets for the SAME event. That is the comparison a cross-market projection performs, which is why the coefficients are the right ones to carry into it. Returns the coefficients plus the diagnostics needed to judge them. """ X, Y, W, meta = [], [], [], [] cells = (gs[gs["ama_000"] > 0] .groupby(["country", "channel", "event_name"])["ama_000"] .agg(["mean", "size"]).reset_index()) cells = cells[cells["size"] >= MIN_CELL_ROWS] for r in cells.itertuples(index=False): t = tvu.get(r.country) s = share.get((r.country, r.channel)) fl = flagship.get(r.country) a = affinity_of(r.country, r.event_name) if not (t and s and fl and a and r.mean > 0): continue X.append([np.log(t), np.log(s / fl), np.log(a)]) Y.append(np.log(float(r.mean))) W.append(np.sqrt(float(r.size))) meta.append((r.event_name, r.country)) if len(X) < 60: return {"ok": False, "reason": f"only {len(X)} usable cells"} X = np.array(X); Y = np.array(Y); W = np.array(W) ev = np.array([m[0] for m in meta]) # Absorb the event fixed effect by weighted within-event demeaning. Xd, Yd, Wd = [], [], [] for e in np.unique(ev): m = ev == e if m.sum() < 3 or len({meta[i][1] for i in np.where(m)[0]}) < 3: continue # an event needs 3+ markets to inform a slope w = W[m] Xd.append(X[m] - np.average(X[m], axis=0, weights=w)) Yd.append(Y[m] - np.average(Y[m], weights=w)) Wd.append(w) if not Xd: return {"ok": False, "reason": "no event has 3+ markets"} Xd = np.vstack(Xd); Yd = np.concatenate(Yd); Wd = np.concatenate(Wd) beta, *_ = np.linalg.lstsq(Xd * Wd[:, None], Yd * Wd, rcond=None) resid = Yd - Xd @ beta sigma = float(np.sqrt(np.average(resid ** 2, weights=Wd))) ss_tot = float(np.sum(Wd * Yd ** 2)) r2 = 1 - float(np.sum(Wd * resid ** 2)) / ss_tot if ss_tot > 0 else float("nan") # Standard errors, so the panel can see which coefficients are real. n, k = Xd.shape xtx_inv = np.linalg.pinv((Xd * Wd[:, None]).T @ (Xd * Wd[:, None])) se = np.sqrt(np.diag(xtx_inv) * (np.sum(Wd * resid ** 2) / max(n - k, 1))) return {"ok": True, "tvu": float(beta[0]), "share": float(beta[1]), "affinity": float(beta[2]), "se_tvu": float(se[0]), "se_share": float(se[1]), "se_affinity": float(se[2]), "sigma": sigma, "r2_within": float(r2), "n_cells": int(n), "n_events": int(len(np.unique(ev)))} # ------------------------------------------------------------ sigma / calibration def residual_stats(pred: np.ndarray, truth: np.ndarray, tier: np.ndarray, shrink_n: float = 8.0) -> tuple[dict, dict]: """Per-tier log-residual sigma and a shrunk multiplicative correction. correction_tier = exp( w * mean(log truth/pred)_tier + (1-w) * global mean ), w = n / (n + shrink_n). That is a James-Stein / empirical-Bayes estimator: a tier with few backtest rows is pulled toward the global correction rather than trusting its own noisy mean. Without it a tier seen ten times would be calibrated as confidently as one seen ten thousand times. """ ok = np.isfinite(pred) & np.isfinite(truth) & (pred > 0) & (truth > 0) lr = np.log(truth[ok] / pred[ok]) tr = tier[ok] mu_g = float(np.mean(lr)) if len(lr) else 0.0 corr, sig = {}, {} for t in np.unique(tr): v = lr[tr == t] if len(v) < 8: continue w = len(v) / (len(v) + shrink_n) corr[int(t)] = float(np.exp(w * float(np.mean(v)) + (1 - w) * mu_g)) # Upper bound raised from 1.50: leave-one-market-out measured 2.45 for the # cross-market tier, and the old ceiling silently capped it -- reporting a # narrower interval than the evidence supports is the failure this whole # calibration exists to remove. sig[int(t)] = float(np.clip(np.std(v, ddof=1), 0.05, 3.00)) return corr, sig