#!/usr/bin/env python3 """ metrics.py — retrieval and ranking metrics, derived rather than imported. Used by Project 2 (ANN recall), Project 8 (recommender quality) and Projects 9/10 (simulator and A/B analysis). Every function here is short on purpose: you should be able to read the formula off the code, because in an interview or a design review you will be asked to define NDCG from memory and "I call sklearn" is not an answer. A warning that governs all of them: **a quality metric is only meaningful against a stated ground truth and a stated k.** "Recall 0.95" is not a result. "recall@10 = 0.95 against exact cosine ground truth on 100k SIFT-1M vectors, efSearch=64" is a result. python3 metrics.py demo """ from __future__ import annotations import math from collections import Counter from typing import Iterable, Sequence # --------------------------------------------------------------------------- # Set metrics -- do the returned items overlap the correct items at all? # --------------------------------------------------------------------------- def recall_at_k(retrieved: Sequence[int], relevant: Iterable[int], k: int) -> float: """|retrieved[:k] ∩ relevant| / min(k, |relevant|). For ANN evaluation `relevant` is the exact top-k from brute force, so the denominator is k and this is the standard "ANN recall". For a recommender, `relevant` is the set of items the user actually engaged with, |relevant| can be < k or > k, and the min() keeps the metric in [0,1] rather than punishing you for a user who only clicked one thing. """ rel = set(relevant) if not rel or k <= 0: return 0.0 hits = sum(1 for d in retrieved[:k] if d in rel) return hits / min(k, len(rel)) def precision_at_k(retrieved: Sequence[int], relevant: Iterable[int], k: int) -> float: """Fraction of the k slots you spent that were worth spending. Recall asks "did you find the good stuff"; precision asks "did you waste the user's screen". A recommender that returns 10 items of which 2 are good has precision@10 = 0.2 no matter how many good items existed. """ rel = set(relevant) if k <= 0: return 0.0 return sum(1 for d in retrieved[:k] if d in rel) / k # --------------------------------------------------------------------------- # Rank-aware metrics -- does position matter? # --------------------------------------------------------------------------- def dcg(gains: Sequence[float]) -> float: """sum_i gain_i / log2(i + 2), i zero-indexed. The log2(i+2) discount is a *modelling choice*, not a law: it encodes the assumption that a user's probability of examining position i decays logarithmically. If your product is an infinite scroll feed rather than ten blue links, that assumption is wrong and you should say so before quoting the number. """ return sum(g / math.log2(i + 2) for i, g in enumerate(gains)) def ndcg_at_k(retrieved: Sequence[int], gains: dict[int, float], k: int) -> float: """DCG of your ranking / DCG of the best possible ranking. Normalising by the ideal makes NDCG comparable across queries with different numbers of relevant items -- which is exactly what raw DCG is not, and the reason you almost never report raw DCG. """ got = [gains.get(d, 0.0) for d in retrieved[:k]] ideal = sorted(gains.values(), reverse=True)[:k] idcg = dcg(ideal) return dcg(got) / idcg if idcg > 0 else 0.0 def mrr(retrieved: Sequence[int], relevant: Iterable[int]) -> float: """1 / (rank of first relevant item), 0 if none. MRR only looks at the first hit, which makes it the right metric for "one correct answer" tasks (a lookup, a navigational query) and the wrong metric for a feed, where the 2nd through 10th items are most of the value. """ rel = set(relevant) for i, d in enumerate(retrieved): if d in rel: return 1.0 / (i + 1) return 0.0 # --------------------------------------------------------------------------- # Catalogue-level metrics -- what the system does across ALL users # --------------------------------------------------------------------------- def coverage(all_recommendations: Sequence[Sequence[int]], catalogue_size: int) -> float: """Fraction of the catalogue that was shown to at least one user. Coverage is the metric that catches the failure NDCG rewards: a recommender that shows the same 50 popular articles to everyone can score well on accuracy and still be worthless, because it is a bestseller list with extra steps. """ shown = set() for rec in all_recommendations: shown.update(rec) return len(shown) / catalogue_size if catalogue_size else 0.0 def novelty(all_recommendations: Sequence[Sequence[int]], popularity: dict[int, int]) -> float: """Mean self-information -log2(p(item)) of recommended items. An item shown to 1% of users carries -log2(0.01) = 6.6 bits; an item shown to everyone carries 0 bits. Averaging over recommendations gives "how surprising is this feed", in bits. Higher is more novel, and more novel is not automatically better -- it trades against accuracy, which is why you report both. """ total = sum(popularity.values()) if not total: return 0.0 vals = [] for rec in all_recommendations: for d in rec: p = popularity.get(d, 0) / total if p > 0: vals.append(-math.log2(p)) return sum(vals) / len(vals) if vals else 0.0 def intra_list_diversity(items: Sequence[int], sim) -> float: """1 - mean pairwise similarity within one result list. `sim(a, b)` returns similarity in [0,1]. This is the metric that catches "ten articles about the same news story", which is invisible to every accuracy metric because all ten are genuinely relevant. """ n = len(items) if n < 2: return 0.0 tot, cnt = 0.0, 0 for i in range(n): for j in range(i + 1, n): tot += sim(items[i], items[j]) cnt += 1 return 1.0 - tot / cnt def gini(all_recommendations: Sequence[Sequence[int]]) -> float: """Gini coefficient of the exposure distribution: 0 = perfectly even, 1 = all exposure on one item. The concentration number to put next to coverage; coverage says how many items were shown at all, Gini says how unevenly the impressions were split among them.""" counts = Counter() for rec in all_recommendations: counts.update(rec) xs = sorted(counts.values()) n = len(xs) if n == 0: return 0.0 cum = sum((2 * (i + 1) - n - 1) * x for i, x in enumerate(xs)) return cum / (n * sum(xs)) # --------------------------------------------------------------------------- # Statistics for Project 10 # --------------------------------------------------------------------------- def welch_t(a: Sequence[float], b: Sequence[float]) -> tuple[float, float]: """Welch's t statistic and degrees of freedom (unequal variances). Student's t assumes equal variance in the two arms. Treatment arms in a recommender routinely have different variance -- a more aggressive ranker produces a wider spread of session lengths -- so Welch is the default and Student is the special case. """ na, nb = len(a), len(b) ma, mb = sum(a) / na, sum(b) / nb va = sum((x - ma) ** 2 for x in a) / (na - 1) vb = sum((x - mb) ** 2 for x in b) / (nb - 1) se2 = va / na + vb / nb t = (mb - ma) / math.sqrt(se2) df = se2 ** 2 / ((va / na) ** 2 / (na - 1) + (vb / nb) ** 2 / (nb - 1)) return t, df def required_n_per_arm(baseline_std: float, mde_abs: float, alpha: float = 0.05, power: float = 0.8) -> int: """n per arm = 2 * (z_{1-a/2} + z_{power})^2 * sigma^2 / delta^2. At the conventional alpha=0.05 / power=0.8 the bracket is (1.96 + 0.8416)^2 = 7.849, which is where the folklore "16 sigma^2 / delta^2" comes from (2 * 7.849 ~ 15.7). Memorise the derivation, not the 16. The practical consequence: halving your minimum detectable effect QUADRUPLES the sample you need. That is the single fact that kills most proposed experiments, and you should compute it before building the variant, not after. """ z_a = 1.959963984540054 # Phi^-1(0.975) z_b = 0.8416212335729143 # Phi^-1(0.80) return math.ceil(2 * (z_a + z_b) ** 2 * baseline_std ** 2 / mde_abs ** 2) def srm_chi2(observed: Sequence[int], expected_ratio: Sequence[float]) -> float: """Chi-square statistic for sample-ratio mismatch. You assigned 50/50 and got 50.4/49.6 on 400k users: is that chance, or is your bucketing broken? With 1 degree of freedom, chi2 > 10.83 is p < 0.001, the conventional SRM alarm threshold. An SRM invalidates the experiment outright -- it means the arms are not comparable populations, so no amount of clever analysis on the metric will save it. """ n = sum(observed) total_r = sum(expected_ratio) exp = [n * r / total_r for r in expected_ratio] return sum((o - e) ** 2 / e for o, e in zip(observed, exp)) def _demo() -> int: print("== rank metrics ==") truth = [7, 3, 11, 2, 9] # exact top-5 from brute force approx = [7, 3, 42, 2, 88, 9, 11, 5, 6, 1] # what the ANN index returned for k in (1, 5, 10): print(f" recall@{k:<2} = {recall_at_k(approx, truth, k):.3f} " f"precision@{k:<2} = {precision_at_k(approx, truth, k):.3f}") gains = {7: 3.0, 3: 2.0, 11: 2.0, 2: 1.0, 9: 1.0} print(f" ndcg@5 = {ndcg_at_k(approx, gains, 5):.4f}") print(f" ndcg@10 = {ndcg_at_k(approx, gains, 10):.4f}") print(f" mrr = {mrr(approx, truth):.4f}") print(" note ndcg@10 > ndcg@5: two relevant items sat at ranks 6 and 7.") print(" A metric reported without its k is not a metric.") print("\n== catalogue metrics: accuracy is not enough ==") catalogue = 1000 popular_only = [[1, 2, 3, 4, 5] for _ in range(200)] personalised = [[(u * 7 + i) % catalogue for i in range(5)] for u in range(200)] pop = Counter() for r in popular_only: pop.update(r) for name, recs in (("bestseller list", popular_only), ("personalised ", personalised)): p = Counter() for r in recs: p.update(r) print(f" {name}: coverage={coverage(recs, catalogue):.3f} " f"gini={gini(recs):.3f} novelty={novelty(recs, p):.2f} bits") print("\n== experiment sizing ==") for mde in (0.05, 0.025, 0.0125): n = required_n_per_arm(baseline_std=0.5, mde_abs=mde) print(f" sigma=0.5, MDE={mde:<7} -> n={n:>9,} per arm") print(" Halving the MDE quadruples n. Check the ratios above: they are 4x.") print("\n== sample ratio mismatch ==") for obs in ([200000, 200000], [201000, 199000], [202000, 198000]): c = srm_chi2(obs, [0.5, 0.5]) flag = "ALARM (p<0.001)" if c > 10.83 else "ok" print(f" {obs} -> chi2={c:8.3f} {flag}") t, df = welch_t([0.30] * 40 + [0.34] * 60, [0.33] * 40 + [0.37] * 60) print(f"\n== welch t = {t:.3f}, df = {df:.1f} ==") return 0 if __name__ == "__main__": import sys raise SystemExit(_demo() if len(sys.argv) > 1 and sys.argv[1] == "demo" else print(__doc__))