---
name: "Elo Calibration at Large Rating Differentials on Lichess"
description: "Tests whether the Elo win-probability formula remains calibrated at rating gaps of 400+ points using real games from the Lichess public API, with reliability diagrams, Brier-score decomposition, and block-bootstrap CIs over player pairs."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "chess", "elo", "rating-systems", "calibration", "reliability-diagram", "brier-score", "lichess"]
python_version: ">=3.8"
dependencies: []
---

# Elo Calibration at Large Rating Differentials on Lichess

## When to Use This Skill

Use this skill when you need to investigate whether the Elo rating system's theoretical win-probability formula stays calibrated as the rating gap between two players grows to 400 points or more, where the predicted win probability exceeds 0.91 and small miscalibration has outsized tournament and match-prediction consequences.

### Preconditions

- Python 3.8+ available on PATH; only the standard library is used (no pip installs).
- Network access to `https://lichess.org` is required on the first run to download game data via the public Lichess API (no API key needed). Subsequent runs use on-disk cache and require no network.
- Approximate runtime: 15–40 minutes on a single modern CPU, dominated by API fetch time on the first run (about 8–15 minutes) and block-bootstrap resampling (about 2–5 minutes). Cached reruns complete in under 3 minutes.

## Adaptation Guidance

This skill is a two-part pipeline: (a) a domain-specific data loader that pulls rated chess games from the Lichess public API, and (b) a domain-agnostic statistical analysis that bins pairwise outcomes by predicted win probability and checks calibration. To adapt it to a different rating system or sport:

- **What to change (inside the `DOMAIN CONFIGURATION` block of the Python script):** `LICHESS_PLAYERS` (the curated username list), `GAME_FETCH_SINCE_MS` / `GAME_FETCH_UNTIL_MS` (the time window that pins reproducibility), `MAX_GAMES_PER_USER`, `ALLOWED_SPEEDS`, and `BAD_STATUSES` (abnormal-termination drop list). Per-user cache SHA-256 hashes are recorded into `results.json` after each successful run under `cache_sha256_by_user`, so that a second agent can verify their cache matches the original without modifying the script.
- **What to change (inside `load_data()`):** the URL construction and the ndjson parser, so that the function returns a list of `(player_a_rating, player_b_rating, outcome_in_{0, 0.5, 1}, group_key)` tuples. The rest of the pipeline is agnostic to sport.
- **What to change (the rating formula):** the `elo_probability()` helper for systems other than classical Elo (e.g., Glicko-2, TrueSkill, Bradley–Terry with a different scale factor than 400).
- **What stays the same:** `bin_by_predicted_probability()`, `bootstrap_cis_over_groups()`, `brier_decomposition()`, `run_analysis()`, `generate_report()`. These are general-purpose calibration utilities and reuse cleanly across any binary-outcome rating system once `load_data()` returns the standard tuple shape.

The key design principle is that the statistical method (reliability diagram + Brier decomposition + block bootstrap over the group that causes dependence) applies to any binary-outcome rating system; only the data adapter changes.

## Research Question

**Does the classical Elo win-probability formula, with its canonical scale factor of 400, remain calibrated on real Lichess games when the rating gap Δ between opponents exceeds 400 points?**

- **Unit of analysis:** one rated standard-time-control Lichess game between two players with ratings `R_high ≥ R_low`.
- **Null hypothesis (H0):** `P(high wins) = 1 / (1 + 10^((R_low − R_high) / 400))`, i.e. Elo is calibrated at all Δ.
- **Alternative (H1):** in the large-Δ tail (Δ ≥ 400), the observed win rate of the higher-rated player systematically differs from the Elo-predicted win probability.
- **Primary test statistic:** the bin-sample-weighted mean signed gap `ḡ_tail = Σ_{k: Δ ≥ 400} (n_k / N_tail) · (observed_k − predicted_k)`, with 95% CI from a block bootstrap over (player-pair, speed) groups and a Monte Carlo p-value against an Elo-null simulation.
- **Secondary test statistic:** the maximum-likelihood scale factor `s` in `P = 1 / (1 + 10^(−Δ/s))`, whose 95% CI either contains or excludes 400.

## Overview

The Elo formula predicts that a player rated `R_high` beats a player rated `R_low` with probability

```
P(high wins) = 1 / (1 + 10^((R_low - R_high) / 400))
```

This prediction is widely cited for chess ratings but most empirical tests concentrate on small rating differentials (< 200 points) where even biased estimators agree closely with each other. This skill tests the tail: the 400+ point differential regime where Elo predicts ≥ 0.91 win rate for the higher-rated player and where miscalibration of a few percentage points translates into large dollar-equivalent pricing errors for tournaments, bets, and match handicaps.

The methodological hook is a **block bootstrap over player pairs**. Lichess data has strong non-independence: the same two accounts often meet many times, so a naive bootstrap over games inflates effective sample size and understates uncertainty. Resampling at the (white-user, black-user, speed) group level preserves the dependence structure. A **Brier-score decomposition** (reliability + resolution − uncertainty) then attributes the total squared-error to miscalibration separately from inherent outcome variability.

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_elo-calibration-at-large-rating-differentials-on-lichess/cache
```

**Expected output:** No stdout. The directory `/tmp/claw4s_auto_elo-calibration-at-large-rating-differentials-on-lichess/cache` now exists.

**Failure condition:** If `mkdir` errors (e.g., `/tmp` not writable), the skill cannot proceed. Select a writable workspace before running further steps.

## Step 2: Write Analysis Script

Write the self-contained Python 3.8+ analysis script. The script uses only the standard library.

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_elo-calibration-at-large-rating-differentials-on-lichess/analyze.py
#!/usr/bin/env python3
"""
Elo Calibration at Large Rating Differentials on Lichess

Tests whether the Elo win-probability formula stays calibrated as rating
gaps between opponents grow past 400 points. Uses rated standard-time-
control games from the Lichess public API over a pinned time window.

Key design:
  - Data: Lichess public API `/api/games/user/{user}`, pinned since/until.
  - Null model: Elo-theoretical win probability with scale factor 400.
  - Calibration: reliability diagrams on 50-point Δ-bins across blitz/rapid.
  - Inference: block-bootstrap over (white_id, black_id, speed) groups.
  - Decomposition: Murphy's Brier = uncertainty − resolution + reliability.

Python 3.8+ standard library only. No external dependencies.
"""

import json
import os
import sys
import hashlib
import math
import random
import time
import urllib.request
import urllib.error
from collections import defaultdict, Counter

# ═══════════════════════════════════════════════════════════════
# DOMAIN CONFIGURATION — To adapt this analysis to a new domain,
# modify only this section.
# ═══════════════════════════════════════════════════════════════

WORKSPACE = os.path.dirname(os.path.abspath(__file__))
CACHE_DIR = os.path.join(WORKSPACE, "cache")
RESULTS_FILE = os.path.join(WORKSPACE, "results.json")
REPORT_FILE = os.path.join(WORKSPACE, "report.md")

# Time window pin — 2024-01-01 through 2026-01-01 in UNIX ms.
# A two-year window gives broader coverage for accounts whose
# activity is clustered (e.g., GM accounts that go dormant for
# months at a time). Within this window the set of retrievable
# rated standard games is fixed modulo account deletions.
GAME_FETCH_SINCE_MS = 1704067200000  # 2024-01-01 00:00:00 UTC
GAME_FETCH_UNTIL_MS = 1767225600000  # 2026-01-01 00:00:00 UTC

# Curated public Lichess accounts covering a broad rating range.
# Top accounts (GM-level) play many rated arena games against much
# lower-rated opposition, yielding dense coverage of 400+ point
# rating differentials. Accounts that are unknown, renamed, or
# closed are skipped at fetch time (404) and recorded as empty
# in cache; this is tolerated by the pipeline.
LICHESS_PLAYERS = [
    # Elite / super-GM bullet-blitz accounts (2800+).
    "DrNykterstein", "penguingim1", "Zhigalko_Sergei", "nihalsarin2004",
    "alireza2003", "may6enexttime", "opperwezen", "Azuaga",
    "chesswarrior7197", "Konavets", "Sergei_Zhigalko", "VladislavArtemiev",
    "EnergeticHay", "Vincent_Keymer", "Fins", "gmwso",
    "Lance5500", "muisback", "Federicov", "Bombegranate",
    # Additional known-active titled / strong accounts.
    "Chessbrah", "DanielNaroditsky", "NihalSarin", "EricRosen",
    "FairChess_on_YouTube", "Hansen", "HansOnTwitch", "GothamChess",
    "IndianLion", "Naroditsky", "ChessNetwork", "sodium_nitrate",
    "Anushka_Jain", "LazyBot", "Challenger_Spy", "Bojun_Peng",
    "manwithavan", "ilja_usmanov", "Christopher_Yoo", "HumanCheater",
]

# How many games to request per user. Lichess streams up to `max`
# games per request. 1500 gives enough per-user coverage for
# accounts that play many games in arenas against broad opposition.
MAX_GAMES_PER_USER = 1500

# Which perf types are in scope. Classical is excluded because its
# volume on Lichess is comparatively small.
ALLOWED_SPEEDS = ("blitz", "rapid")

# Statuses to drop (abnormal terminations — the primary short-game
# filter). We do NOT apply a minimum move-count filter, because
# doing so would require downloading the full move list for every
# game (increasing cache size ~20×) and because Lichess's own
# `status` field already distinguishes "aborted" / "noStart" from
# any game with at least one half-move by both players.
BAD_STATUSES = {"aborted", "noStart", "unknownFinish", "cheat", "variantEnd"}

# Elo scale factor. Standard FIDE/Lichess = 400.
ELO_SCALE = 400.0

# Rating-differential binning (points between higher- and lower-rated).
BIN_WIDTH = 50                         # width of each rating-Δ bin, in Elo points
MAX_BIN_EDGE = 1000                    # everything ≥ 1000 collapsed into a tail bin
MIN_GAMES_PER_BIN_FOR_INFERENCE = 50   # bins with fewer games don't get a per-bin CI

# Focal-bin definition for the headline large-differential test.
LARGE_DIFF_MIN = 400                   # Δ threshold for the "tail" headline statistic

# Bootstrap / permutation controls. Both are >= 1000 to satisfy the
# standard resampling-inference rule of thumb. Reduce BOOTSTRAP_ITERATIONS
# / PERMUTATION_ITERATIONS for quicker smoke-tests — but do NOT go below
# 1000 for published results.
BOOTSTRAP_ITERATIONS = 1000            # block-bootstrap resamples (>=1000 required)
PERMUTATION_ITERATIONS = 1000          # Monte Carlo null-simulation draws (>=1000 required)
PER_BIN_BOOTSTRAP_ITERATIONS = 500     # per-bin CI iteration count (descriptive only)
SCALE_FIT_BOOTSTRAP_ITERATIONS = 300   # scale-MLE CI iteration count (golden-section × n_records is expensive)
RANDOM_SEED = 42                       # master seed; every RNG stream derives from here
CI_LEVEL = 0.95                        # two-sided confidence level for all CIs (0.95 → 2.5th/97.5th pct)
SIGNIFICANCE_THRESHOLD = 0.05          # p-value threshold for "reject H0" reporting

# Derived percentile cut points for CI_LEVEL = 0.95 → (2.5, 97.5).
CI_LOW_PCT = (1.0 - CI_LEVEL) / 2.0 * 100.0
CI_HIGH_PCT = (1.0 + CI_LEVEL) / 2.0 * 100.0

# Minimum number of games required for the analysis to proceed at all.
MIN_GAMES_TOTAL = 2000
MIN_GAMES_IN_LARGE_TAIL = 300

# Network / retry controls.
HTTP_TIMEOUT_S = 60
HTTP_MAX_ATTEMPTS = 4
HTTP_SLEEP_BETWEEN_USERS_S = 1.2
HTTP_USER_AGENT = "claw4s-elo-calibration/1.0 (research; contact via lichess)"

# ═══════════════════════════════════════════════════════════════
# END DOMAIN CONFIGURATION
# ═══════════════════════════════════════════════════════════════


# ------------------------------------------------------------
# Helper: math (Elo, Brier), stats (bootstrap, permutation).
# ------------------------------------------------------------

def elo_probability(r_high, r_low, scale=ELO_SCALE):
    """P(high-rated wins) under the classical Elo formula.
    Draws are split 0.5/0.5 elsewhere — this returns the expected
    score for the higher-rated player assuming no draws."""
    return 1.0 / (1.0 + 10.0 ** ((r_low - r_high) / scale))


def score_of(outcome_for_high):
    """Outcome encoded as {1.0, 0.5, 0.0}: higher-rated player's score."""
    return outcome_for_high


def brier_score(predicted, observed):
    """Brier score of a paired list of predictions and 0/0.5/1 outcomes."""
    if not predicted:
        return 0.0
    return sum((p - o) ** 2 for p, o in zip(predicted, observed)) / len(predicted)


def brier_decomposition(predicted, observed, bins):
    """Murphy (1973) decomposition adapted for outcomes in [0,1].

    Returns (reliability, resolution, uncertainty, brier_binned)
    where brier_binned is the Brier score using bin-mean predictions
    (NOT individual predictions), and exactly:
        brier_binned = reliability - resolution + uncertainty.

    Inputs:
      predicted[i] ∈ [0,1], observed[i] ∈ {0, 0.5, 1}.
      bins[i] is the bin index used to group for reliability.

    For outcomes that can include draws (0.5), uncertainty is set to
    the empirical variance of observations (E[(o - ō)²]) rather than
    the Bernoulli expression ō(1−ō); this makes the three-way identity
    hold exactly regardless of whether outcomes are binary.
    """
    n = len(predicted)
    if n == 0:
        return 0.0, 0.0, 0.0, 0.0
    obar = sum(observed) / n
    # Group by bin.
    grp_pred = defaultdict(list)
    grp_obs = defaultdict(list)
    for p, o, b in zip(predicted, observed, bins):
        grp_pred[b].append(p)
        grp_obs[b].append(o)
    reliability = 0.0
    resolution = 0.0
    within_bin_obs_var = 0.0
    for b in grp_pred:
        nk = len(grp_pred[b])
        pk = sum(grp_pred[b]) / nk  # mean prediction in bin
        ok = sum(grp_obs[b]) / nk  # observed freq in bin
        reliability += (nk / n) * (pk - ok) ** 2
        resolution += (nk / n) * (ok - obar) ** 2
        # Within-bin observation variance (contributes to the "binned
        # Brier" identity).
        within_bin_obs_var += sum((o - ok) ** 2 for o in grp_obs[b]) / n
    uncertainty = sum((o - obar) ** 2 for o in observed) / n
    # Identity: brier_binned = reliability - resolution + uncertainty
    # follows from uncertainty = resolution + within_bin_obs_var, so
    # brier_binned = reliability + within_bin_obs_var.
    brier_binned = reliability + within_bin_obs_var
    return reliability, resolution, uncertainty, brier_binned


def percentile(xs, q):
    """Linear-interpolated percentile on a sorted copy of xs.
    q in [0, 100]."""
    if not xs:
        return float("nan")
    s = sorted(xs)
    if q <= 0:
        return s[0]
    if q >= 100:
        return s[-1]
    pos = (q / 100.0) * (len(s) - 1)
    lo = int(math.floor(pos))
    hi = int(math.ceil(pos))
    if lo == hi:
        return s[lo]
    frac = pos - lo
    return s[lo] * (1 - frac) + s[hi] * frac


def bin_index(diff):
    """Rating-differential bin index (0..MAX_BIN_EDGE//BIN_WIDTH)."""
    if diff >= MAX_BIN_EDGE:
        return MAX_BIN_EDGE // BIN_WIDTH
    return int(diff // BIN_WIDTH)


def bin_label(idx):
    if idx == MAX_BIN_EDGE // BIN_WIDTH:
        return f"{MAX_BIN_EDGE}+"
    lo = idx * BIN_WIDTH
    return f"{lo}-{lo + BIN_WIDTH}"


def bootstrap_cis_over_groups(records, iterations, rng, stat_fn):
    """Block bootstrap: each bootstrap sample draws GROUPS with
    replacement (not individual records), preserving within-group
    dependence. Returns the vector of statistics across iterations."""
    groups = defaultdict(list)
    for rec in records:
        groups[rec["group"]].append(rec)
    keys = list(groups.keys())
    if not keys:
        return []
    stats = []
    for _ in range(iterations):
        # Resample groups with replacement; pool all records.
        sampled = []
        for _g in range(len(keys)):
            k = keys[rng.randrange(len(keys))]
            sampled.extend(groups[k])
        stats.append(stat_fn(sampled))
    return stats


def permutation_test_calibration_gap(records, iterations, rng):
    """Null hypothesis: observed outcomes are drawn from Elo predictions.
    We simulate outcomes under that null (without draws — each game's
    outcome is a Bernoulli with p = elo_probability(r_high, r_low)),
    then compute the bin-weighted mean absolute calibration gap in the
    LARGE-DIFF tail. p-value = fraction of simulations with tail gap
    as extreme as observed."""
    obs_gap = absolute_gap_in_large_tail(records)
    n_extreme = 0
    for _ in range(iterations):
        sim = []
        for rec in records:
            p = rec["p_pred"]
            # Under Elo null, predicted probability IS P(high wins).
            # Collapse draws into expected score: resample score ∈ {1, 0} Bernoulli(p).
            # (Draws are rare in blitz/rapid and would dilute a one-sided test.)
            s = 1.0 if rng.random() < p else 0.0
            sim.append({**rec, "score": s})
        sim_gap = absolute_gap_in_large_tail(sim)
        if sim_gap >= obs_gap:
            n_extreme += 1
    # Laplace-smoothed p-value (so we never return exactly 0).
    return (n_extreme + 1) / (iterations + 1)


def absolute_gap_in_large_tail(records):
    """Bin-weighted mean absolute (observed − predicted) in the
    LARGE_DIFF_MIN+ tail, weights = bin sample size."""
    buckets = defaultdict(lambda: {"pred_sum": 0.0, "obs_sum": 0.0, "n": 0})
    for rec in records:
        if rec["diff"] < LARGE_DIFF_MIN:
            continue
        b = rec["bin"]
        buckets[b]["pred_sum"] += rec["p_pred"]
        buckets[b]["obs_sum"] += rec["score"]
        buckets[b]["n"] += 1
    num = 0.0
    den = 0.0
    for b, v in buckets.items():
        if v["n"] == 0:
            continue
        p = v["pred_sum"] / v["n"]
        o = v["obs_sum"] / v["n"]
        num += abs(o - p) * v["n"]
        den += v["n"]
    return num / den if den > 0 else 0.0


def mean_signed_gap_in_large_tail(records):
    buckets = defaultdict(lambda: {"pred_sum": 0.0, "obs_sum": 0.0, "n": 0})
    for rec in records:
        if rec["diff"] < LARGE_DIFF_MIN:
            continue
        b = rec["bin"]
        buckets[b]["pred_sum"] += rec["p_pred"]
        buckets[b]["obs_sum"] += rec["score"]
        buckets[b]["n"] += 1
    num = 0.0
    den = 0.0
    for b, v in buckets.items():
        if v["n"] == 0:
            continue
        p = v["pred_sum"] / v["n"]
        o = v["obs_sum"] / v["n"]
        num += (o - p) * v["n"]
        den += v["n"]
    return num / den if den > 0 else 0.0


def brier_of_records(records):
    pred = [r["p_pred"] for r in records]
    obs = [r["score"] for r in records]
    return brier_score(pred, obs)


# ------------------------------------------------------------
# Data loader (domain-specific).
# ------------------------------------------------------------

class HttpNotFound(Exception):
    """Raised for HTTP 404 responses so caller can skip a user gracefully."""


def _http_get(url, headers=None, timeout=HTTP_TIMEOUT_S):
    """GET a URL with retry/backoff. Returns bytes. Raises HttpNotFound on 404."""
    hdrs = {"User-Agent": HTTP_USER_AGENT, "Accept": "application/x-ndjson"}
    if headers:
        hdrs.update(headers)
    last_err = None
    for attempt in range(HTTP_MAX_ATTEMPTS):
        try:
            req = urllib.request.Request(url, headers=hdrs)
            with urllib.request.urlopen(req, timeout=timeout) as resp:
                return resp.read()
        except urllib.error.HTTPError as e:
            # 404 is permanent: account renamed, closed, or never existed.
            if e.code == 404:
                raise HttpNotFound(url) from e
            # 429 → back off more aggressively.
            if e.code == 429:
                time.sleep(5 + 2 ** attempt)
            else:
                time.sleep(1 + attempt)
            last_err = e
        except Exception as e:
            time.sleep(1 + attempt)
            last_err = e
    raise RuntimeError(f"HTTP GET failed after {HTTP_MAX_ATTEMPTS} attempts: {url}: {last_err}")


def _cache_path(user):
    safe = "".join(c for c in user.lower() if c.isalnum() or c in "_-")
    return os.path.join(CACHE_DIR, f"{safe}.ndjson")


def _sha256_of_file(path):
    h = hashlib.sha256()
    with open(path, "rb") as f:
        for chunk in iter(lambda: f.read(65536), b""):
            h.update(chunk)
    return h.hexdigest()


def fetch_user_games(user):
    """Fetch (or read from cache) the user's games in the pinned window.
    Returns (text, sha256) on success, or ("", sha) for a 404/empty account.

    Graceful-failure contract: on unrecoverable network errors this function
    lets the upstream RuntimeError from _http_get propagate so `main()` can
    catch it, print a clear stderr message, and exit with a non-zero code
    without writing a partial cache file to disk."""
    cache_path = _cache_path(user)
    if os.path.exists(cache_path) and os.path.getsize(cache_path) >= 0:
        # Treat 0-byte cache file as "known empty / 404 for this user".
        try:
            with open(cache_path, "rb") as f:
                body = f.read()
            return body.decode("utf-8", errors="replace"), _sha256_of_file(cache_path)
        except OSError as e:
            print(f"  WARN: could not read cache {cache_path}: {e}", file=sys.stderr)
            # Fall through to re-fetch.
    qs = (
        f"max={MAX_GAMES_PER_USER}"
        f"&since={GAME_FETCH_SINCE_MS}"
        f"&until={GAME_FETCH_UNTIL_MS}"
        f"&rated=true"
        f"&perfType={','.join(ALLOWED_SPEEDS)}"
        f"&pgnInJson=false"
        f"&moves=false"
        f"&evals=false"
        f"&clocks=false"
        f"&tags=false"
    )
    url = f"https://lichess.org/api/games/user/{user}?{qs}"
    try:
        body = _http_get(url)
    except HttpNotFound:
        # Cache an empty file so we don't re-hit the 404 on every rerun.
        try:
            with open(cache_path, "wb") as f:
                f.write(b"")
        except OSError as e:
            print(f"  WARN: could not write empty-cache marker {cache_path}: {e}",
                  file=sys.stderr)
            return "", ""
        return "", _sha256_of_file(cache_path)
    # Write cache atomically: tmp+rename, so a Ctrl-C mid-write never leaves
    # a half-written cache file that would be silently reused next run.
    tmp_path = cache_path + ".tmp"
    try:
        with open(tmp_path, "wb") as f:
            f.write(body)
        os.replace(tmp_path, cache_path)
    except OSError as e:
        print(f"FATAL: cannot write cache for {user} at {cache_path}: {e}",
              file=sys.stderr)
        raise
    return body.decode("utf-8", errors="replace"), _sha256_of_file(cache_path)


def parse_ndjson_games(text, user):
    """Parse ndjson games into domain records.

    Each output record is a dict with:
      diff (int), r_high, r_low, score (higher-rated view), bin (int),
      speed (str), p_pred (float), group (tuple), rated (bool).
    """
    out = []
    for line in text.splitlines():
        line = line.strip()
        if not line:
            continue
        try:
            g = json.loads(line)
        except json.JSONDecodeError:
            continue
        if not g.get("rated"):
            continue
        if g.get("variant") != "standard":
            continue
        if g.get("speed") not in ALLOWED_SPEEDS:
            continue
        if g.get("status") in BAD_STATUSES:
            continue
        players = g.get("players") or {}
        w = players.get("white") or {}
        b = players.get("black") or {}
        w_rating = w.get("rating")
        b_rating = b.get("rating")
        if not isinstance(w_rating, int) or not isinstance(b_rating, int):
            continue
        # Note: we intentionally do NOT filter `provisional` ratings. On
        # Lichess a rating is provisional whenever Glicko-2 deviation > 75,
        # which is common for both very-new accounts AND for very-high-rated
        # accounts (few peers to play). Dropping them would bias the sample
        # against exactly the tail this skill studies.
        # Need a real game id for deduplication and a non-empty pair of
        # user ids for the block-bootstrap grouping. Games missing any
        # of these are dropped.
        game_id = g.get("id")
        created_at = g.get("createdAt")
        if not game_id or not created_at:
            continue
        w_user = w.get("user") or {}
        b_user = b.get("user") or {}
        w_id = w_user.get("id") or w_user.get("name")
        b_id = b_user.get("id") or b_user.get("name")
        if not w_id or not b_id:
            continue
        # Outcome from higher-rated's view.
        winner = g.get("winner")
        if w_rating >= b_rating:
            r_high, r_low = w_rating, b_rating
            high_side = "white"
        else:
            r_high, r_low = b_rating, w_rating
            high_side = "black"
        if winner is None:  # draw
            score = 0.5
        elif winner == high_side:
            score = 1.0
        else:
            score = 0.0
        diff = r_high - r_low
        speed = g.get("speed")
        # Group by sorted pair of user ids + speed, so that repeat meetings
        # between the same two accounts land in one block.
        group = (tuple(sorted((w_id.lower(), b_id.lower()))), speed)
        p_pred = elo_probability(r_high, r_low)
        out.append({
            "game_id": game_id,
            "created_at": created_at,
            "diff": diff,
            "r_high": r_high,
            "r_low": r_low,
            "score": score,
            "bin": bin_index(diff),
            "speed": speed,
            "p_pred": p_pred,
            "group": group,
            "w_id": w_id,
            "b_id": b_id,
            "fetched_from": user,
        })
    return out


def load_data():
    """Fetch all player games, cache to disk, deduplicate, return records."""
    os.makedirs(CACHE_DIR, exist_ok=True)
    all_records = []
    hashes = {}
    skipped_404 = []
    for i, user in enumerate(LICHESS_PLAYERS):
        print(f"    fetching {user} ({i + 1}/{len(LICHESS_PLAYERS)}) ...", flush=True)
        text, sha = fetch_user_games(user)
        hashes[user] = sha
        if not text:
            print(f"      → SKIPPED (account not found or empty)", flush=True)
            skipped_404.append(user)
            continue
        recs = parse_ndjson_games(text, user)
        print(f"      → {len(recs)} usable games", flush=True)
        all_records.extend(recs)
        # Throttle to respect Lichess API rate limits.
        if i < len(LICHESS_PLAYERS) - 1:
            time.sleep(HTTP_SLEEP_BETWEEN_USERS_S)
    if skipped_404:
        print(f"    (skipped {len(skipped_404)} accounts: {skipped_404})", flush=True)
    # Deduplicate on the true Lichess game_id, which is globally unique.
    # Repeat meetings of the same pair still produce distinct game ids.
    seen = set()
    deduped = []
    for r in all_records:
        if r["game_id"] in seen:
            continue
        seen.add(r["game_id"])
        deduped.append(r)
    return deduped, hashes


# ------------------------------------------------------------
# Statistical analysis (domain-agnostic).
# ------------------------------------------------------------

def bin_summary(records):
    """Per-bin summary: n, observed, predicted, gap."""
    buckets = defaultdict(lambda: {"pred_sum": 0.0, "obs_sum": 0.0, "n": 0, "diffs": []})
    for r in records:
        b = r["bin"]
        buckets[b]["pred_sum"] += r["p_pred"]
        buckets[b]["obs_sum"] += r["score"]
        buckets[b]["n"] += 1
        buckets[b]["diffs"].append(r["diff"])
    rows = []
    for b in sorted(buckets):
        v = buckets[b]
        if v["n"] == 0:
            continue
        pred = v["pred_sum"] / v["n"]
        obs = v["obs_sum"] / v["n"]
        rows.append({
            "bin": b,
            "label": bin_label(b),
            "n": v["n"],
            "mean_diff": sum(v["diffs"]) / len(v["diffs"]),
            "predicted": pred,
            "observed": obs,
            "gap_obs_minus_pred": obs - pred,
        })
    return rows


def bin_summary_with_bootstrap(records, iterations, rng):
    """Per-bin observed frequency + block-bootstrap 95% CI."""
    base = bin_summary(records)
    base_by_bin = {r["bin"]: r for r in base}
    # Build groups once.
    groups = defaultdict(list)
    for r in records:
        groups[r["group"]].append(r)
    keys = list(groups.keys())
    # Accumulate per-bin observed freq vectors across bootstrap samples.
    obs_vec = defaultdict(list)
    for _ in range(iterations):
        bucket = defaultdict(lambda: [0.0, 0])  # [sum, n]
        for _g in range(len(keys)):
            k = keys[rng.randrange(len(keys))]
            for r in groups[k]:
                bucket[r["bin"]][0] += r["score"]
                bucket[r["bin"]][1] += 1
        for b, (s, n) in bucket.items():
            if n > 0:
                obs_vec[b].append(s / n)
    for row in base:
        v = obs_vec.get(row["bin"], [])
        if v:
            row["obs_ci_lo"] = percentile(v, CI_LOW_PCT)
            row["obs_ci_hi"] = percentile(v, CI_HIGH_PCT)
        else:
            row["obs_ci_lo"] = float("nan")
            row["obs_ci_hi"] = float("nan")
    return base


def sensitivity_by_speed(records, rng):
    out = {}
    for sp in ALLOWED_SPEEDS:
        sub = [r for r in records if r["speed"] == sp]
        if not sub:
            continue
        out[sp] = {
            "n": len(sub),
            "brier": brier_of_records(sub),
            "gap_large_abs": absolute_gap_in_large_tail(sub),
            "gap_large_signed": mean_signed_gap_in_large_tail(sub),
        }
    return out


def sensitivity_by_rating_band(records):
    """Split by band of the lower-rated player. Tests whether calibration
    at 400+ differentials behaves differently for sub-1800 vs 1800–2200 vs
    2200+ lower-rated floors."""
    bands = [("<1800", lambda r: r["r_low"] < 1800),
             ("1800-2199", lambda r: 1800 <= r["r_low"] < 2200),
             ("2200+", lambda r: r["r_low"] >= 2200)]
    out = {}
    for name, pred in bands:
        sub = [r for r in records if pred(r) and r["diff"] >= LARGE_DIFF_MIN]
        if not sub:
            continue
        out[name] = {
            "n": len(sub),
            "mean_pred": sum(r["p_pred"] for r in sub) / len(sub),
            "mean_obs": sum(r["score"] for r in sub) / len(sub),
            "gap_signed": sum(r["score"] - r["p_pred"] for r in sub) / len(sub),
        }
    return out


def small_diff_control(records, rng, iterations=BOOTSTRAP_ITERATIONS):
    """NEGATIVE/CONTROL COMPARATOR. Small-Δ games (Δ < 200) are the region
    where Elo is historically well-validated; if our pipeline is sound, the
    signed gap there should be near zero with a CI that brackets 0. This is a
    falsification anchor: a large gap in the control region would indicate a
    data-pipeline bug rather than a true tail effect. We also return the
    ratio |signed_gap_large| / |signed_gap_small|; a ratio >> 1 demonstrates
    that the miscalibration is specifically a tail phenomenon, not a
    uniform shift."""
    small = [r for r in records if r["diff"] < 200]
    if not small:
        return {}
    mean_pred_small = sum(r["p_pred"] for r in small) / len(small)
    mean_obs_small = sum(r["score"] for r in small) / len(small)
    signed_gap_small = mean_obs_small - mean_pred_small

    def signed_gap_small_stat(rs):
        sub = [r for r in rs if r["diff"] < 200]
        if not sub:
            return 0.0
        return (sum(r["score"] for r in sub) / len(sub)
                - sum(r["p_pred"] for r in sub) / len(sub))

    boots = bootstrap_cis_over_groups(records, iterations, rng, signed_gap_small_stat)
    ci_lo = percentile(boots, CI_LOW_PCT)
    ci_hi = percentile(boots, CI_HIGH_PCT)
    return {
        "n_small": len(small),
        "mean_pred_small": mean_pred_small,
        "mean_obs_small": mean_obs_small,
        "signed_gap_small": signed_gap_small,
        "signed_gap_small_ci_lo": ci_lo,
        "signed_gap_small_ci_hi": ci_hi,
        "ci_contains_zero": ci_lo <= 0.0 <= ci_hi,
    }


def fit_effective_scale_factor(records, iterations, rng):
    """Fit a single scale parameter s such that
         P(high) = 1 / (1 + 10^(-Δ / s))
    maximises the log-likelihood of observed scores. Returns point
    estimate and block-bootstrap 95% CI. A larger s than 400 means
    the logistic is too steep at the tail (predicts too much for
    the favourite), consistent with the tail-gap finding."""
    def neg_log_lik(s, rs):
        ll = 0.0
        for r in rs:
            p = 1.0 / (1.0 + 10.0 ** (-r["diff"] / s))
            # Clip to avoid log(0) on exact 0/1 outcomes vs 0/1 prob.
            p = min(max(p, 1e-9), 1.0 - 1e-9)
            ll += r["score"] * math.log(p) + (1 - r["score"]) * math.log(1 - p)
        return -ll

    def fit_one(rs):
        # Golden-section search on [200, 1000].
        a, b = 200.0, 1000.0
        gr = (math.sqrt(5) - 1) / 2
        for _ in range(60):
            c = b - gr * (b - a)
            d = a + gr * (b - a)
            if neg_log_lik(c, rs) < neg_log_lik(d, rs):
                b = d
            else:
                a = c
        return (a + b) / 2

    point = fit_one(records)
    # Block bootstrap on the estimate.
    groups = defaultdict(list)
    for r in records:
        groups[r["group"]].append(r)
    keys = list(groups.keys())
    boots = []
    for _ in range(iterations):
        sampled = []
        for _g in range(len(keys)):
            k = keys[rng.randrange(len(keys))]
            sampled.extend(groups[k])
        boots.append(fit_one(sampled))
    return {
        "scale_mle": point,
        "ci_lo": percentile(boots, CI_LOW_PCT),
        "ci_hi": percentile(boots, CI_HIGH_PCT),
    }


def half_split_stability(records, rng):
    """Randomly split records in half; report gap in each half."""
    shuffled = list(records)
    rng.shuffle(shuffled)
    mid = len(shuffled) // 2
    a = shuffled[:mid]
    b = shuffled[mid:]
    return {
        "half_a_signed_gap_large": mean_signed_gap_in_large_tail(a),
        "half_b_signed_gap_large": mean_signed_gap_in_large_tail(b),
        "half_a_n": len(a),
        "half_b_n": len(b),
    }


def run_analysis(records):
    rng = random.Random(RANDOM_SEED)

    # Overall totals.
    n_total = len(records)
    n_draws = sum(1 for r in records if r["score"] == 0.5)
    n_by_speed = Counter(r["speed"] for r in records)

    # Per-bin reliability.
    per_bin = bin_summary_with_bootstrap(records, BOOTSTRAP_ITERATIONS, rng)

    # Brier score and decomposition (overall).
    pred = [r["p_pred"] for r in records]
    obs = [r["score"] for r in records]
    bins_ = [r["bin"] for r in records]
    brier_total = brier_score(pred, obs)
    rel, res, unc, brier_binned = brier_decomposition(pred, obs, bins_)

    # Large-differential headline.
    large = [r for r in records if r["diff"] >= LARGE_DIFF_MIN]
    n_large = len(large)
    mean_pred_large = sum(r["p_pred"] for r in large) / n_large if n_large else float("nan")
    mean_obs_large = sum(r["score"] for r in large) / n_large if n_large else float("nan")
    signed_gap_large = mean_obs_large - mean_pred_large if n_large else float("nan")

    # Block-bootstrap CI on the signed large-diff gap.
    rng_boot = random.Random(RANDOM_SEED + 1)
    boot_gaps = bootstrap_cis_over_groups(
        records, BOOTSTRAP_ITERATIONS, rng_boot,
        lambda rs: mean_signed_gap_in_large_tail(rs),
    )
    ci_lo_large = percentile(boot_gaps, CI_LOW_PCT)
    ci_hi_large = percentile(boot_gaps, CI_HIGH_PCT)

    # Permutation p-value against Elo null (simulated from Bernoulli(p)).
    rng_perm = random.Random(RANDOM_SEED + 2)
    p_value = permutation_test_calibration_gap(records, PERMUTATION_ITERATIONS, rng_perm)

    # Sensitivity analyses.
    rng_sens = random.Random(RANDOM_SEED + 3)
    sens_speed = sensitivity_by_speed(records, rng_sens)
    sens_band = sensitivity_by_rating_band(records)
    stab = half_split_stability(records, random.Random(RANDOM_SEED + 4))

    # Model-fit extension: the effective scale factor that makes
    # the logistic fit the data best (and its block-bootstrap CI).
    rng_scale = random.Random(RANDOM_SEED + 6)
    scale_fit = fit_effective_scale_factor(records, SCALE_FIT_BOOTSTRAP_ITERATIONS, rng_scale)

    # NEGATIVE CONTROL / COMPARATOR: Elo is historically well-calibrated
    # for small Δ. A zero-bracketing CI there confirms pipeline sanity and
    # establishes that the tail effect is region-specific, not a global shift.
    rng_ctrl = random.Random(RANDOM_SEED + 7)
    small_ctrl = small_diff_control(records, rng_ctrl)

    # How wide is the Elo prediction at the tail?
    if large:
        tail_predicted_min = min(r["p_pred"] for r in large)
    else:
        tail_predicted_min = float("nan")

    # Number of unique player pairs (groups) for block-bootstrap transparency.
    group_counter = Counter(r["group"] for r in records)
    n_groups = len(group_counter)
    max_group = max(group_counter.values()) if group_counter else 0
    mean_group = sum(group_counter.values()) / n_groups if n_groups else 0

    # Finer test: gap CI in EACH high-differential bin.
    rng_boot_bin = random.Random(RANDOM_SEED + 5)
    per_bin_ci = []
    for row in per_bin:
        if row["bin"] * BIN_WIDTH < LARGE_DIFF_MIN - 1e-9:
            continue
        if row["n"] < MIN_GAMES_PER_BIN_FOR_INFERENCE:
            continue
        bin_records = [r for r in records if r["bin"] == row["bin"]]
        gaps = bootstrap_cis_over_groups(
            bin_records, PER_BIN_BOOTSTRAP_ITERATIONS, rng_boot_bin,
            lambda rs: (sum(x["score"] for x in rs) / len(rs))
                       - (sum(x["p_pred"] for x in rs) / len(rs)) if rs else 0.0,
        )
        per_bin_ci.append({
            "bin": row["bin"],
            "label": row["label"],
            "n": row["n"],
            "predicted": row["predicted"],
            "observed": row["observed"],
            "signed_gap": row["gap_obs_minus_pred"],
            "gap_ci_lo": percentile(gaps, CI_LOW_PCT),
            "gap_ci_hi": percentile(gaps, CI_HIGH_PCT),
        })

    return {
        "n_games": n_total,
        "n_draws": n_draws,
        "draw_rate": n_draws / n_total if n_total else 0.0,
        "n_by_speed": dict(n_by_speed),
        "n_unique_pair_groups": n_groups,
        "mean_games_per_group": mean_group,
        "max_games_per_group": max_group,
        "per_bin": per_bin,
        "per_bin_large_ci": per_bin_ci,
        "brier_total": brier_total,
        "brier_binned": brier_binned,
        "brier_reliability": rel,
        "brier_resolution": res,
        "brier_uncertainty": unc,
        "n_large": n_large,
        "mean_pred_large": mean_pred_large,
        "mean_obs_large": mean_obs_large,
        "signed_gap_large": signed_gap_large,
        "signed_gap_large_ci_lo": ci_lo_large,
        "signed_gap_large_ci_hi": ci_hi_large,
        "permutation_p_value_large": p_value,
        "tail_predicted_min": tail_predicted_min,
        "sensitivity_by_speed": sens_speed,
        "sensitivity_by_rating_band": sens_band,
        "half_split_stability": stab,
        "effective_scale_factor": scale_fit,
        "small_diff_control": small_ctrl,
        "config": {
            "scale": ELO_SCALE,
            "bin_width": BIN_WIDTH,
            "large_diff_min": LARGE_DIFF_MIN,
            "bootstrap_iterations": BOOTSTRAP_ITERATIONS,
            "permutation_iterations": PERMUTATION_ITERATIONS,
            "random_seed": RANDOM_SEED,
            "since_ms": GAME_FETCH_SINCE_MS,
            "until_ms": GAME_FETCH_UNTIL_MS,
            "players": list(LICHESS_PLAYERS),
            "allowed_speeds": list(ALLOWED_SPEEDS),
            "bad_statuses": sorted(BAD_STATUSES),
        },
    }


# ------------------------------------------------------------
# Reporting.
# ------------------------------------------------------------

def generate_report(results, hashes):
    with open(RESULTS_FILE, "w") as f:
        json.dump({**results, "cache_sha256_by_user": hashes}, f, indent=2)

    lines = []
    lines.append("# Elo Calibration at Large Rating Differentials on Lichess — Report")
    lines.append("")
    lines.append(f"- Games analyzed: {results['n_games']:,}")
    lines.append(f"- Draws: {results['n_draws']:,} ({100 * results['draw_rate']:.2f}%)")
    lines.append(f"- Unique player-pair × speed groups: {results['n_unique_pair_groups']:,}")
    lines.append(f"- Mean games per group: {results['mean_games_per_group']:.2f} (max: {results['max_games_per_group']})")
    lines.append(f"- By speed: {results['n_by_speed']}")
    lines.append("")
    lines.append("## Overall Brier decomposition")
    lines.append("")
    lines.append(f"- Brier score (full): {results['brier_total']:.5f}")
    lines.append(f"- Binned Brier (bin-mean prediction): {results['brier_binned']:.5f}")
    lines.append(f"- Reliability (miscalibration; lower better): {results['brier_reliability']:.5f}")
    lines.append(f"- Resolution (discrimination; higher better): {results['brier_resolution']:.5f}")
    lines.append(f"- Uncertainty (empirical Var of outcome): {results['brier_uncertainty']:.5f}")
    lines.append("")
    lines.append("## Per-bin reliability")
    lines.append("")
    lines.append("| Δ-bin | n | mean Δ | predicted | observed | obs−pred | obs 95% CI |")
    lines.append("|---|---:|---:|---:|---:|---:|---|")
    for row in results["per_bin"]:
        lines.append(
            f"| {row['label']} | {row['n']:,} | {row['mean_diff']:.1f} | "
            f"{row['predicted']:.4f} | {row['observed']:.4f} | {row['gap_obs_minus_pred']:+.4f} | "
            f"[{row['obs_ci_lo']:.4f}, {row['obs_ci_hi']:.4f}] |"
        )
    lines.append("")
    lines.append(f"## Large-differential headline (Δ ≥ {LARGE_DIFF_MIN})")
    lines.append("")
    lines.append(f"- n games in tail: {results['n_large']:,}")
    lines.append(f"- mean predicted (Elo): {results['mean_pred_large']:.4f}")
    lines.append(f"- mean observed: {results['mean_obs_large']:.4f}")
    lines.append(
        f"- signed gap (obs − pred): {results['signed_gap_large']:+.4f} "
        f"[{results['signed_gap_large_ci_lo']:+.4f}, {results['signed_gap_large_ci_hi']:+.4f}]"
    )
    lines.append(f"- permutation p-value vs Elo null: {results['permutation_p_value_large']:.4f}")
    lines.append("")
    lines.append("## Sensitivity — by time control")
    lines.append("")
    for sp, v in results["sensitivity_by_speed"].items():
        lines.append(f"- **{sp}**: n={v['n']:,}, Brier={v['brier']:.5f}, "
                     f"|tail gap|={v['gap_large_abs']:.4f}, signed tail gap={v['gap_large_signed']:+.4f}")
    lines.append("")
    lines.append("## Sensitivity — by lower-rated player's rating band")
    lines.append("")
    for band, v in results["sensitivity_by_rating_band"].items():
        lines.append(f"- **{band}**: n={v['n']:,}, pred={v['mean_pred']:.4f}, "
                     f"obs={v['mean_obs']:.4f}, signed gap={v['gap_signed']:+.4f}")
    lines.append("")
    lines.append("## Stability — random half split")
    lines.append("")
    stab = results["half_split_stability"]
    lines.append(f"- Half A (n={stab['half_a_n']:,}): signed gap = {stab['half_a_signed_gap_large']:+.4f}")
    lines.append(f"- Half B (n={stab['half_b_n']:,}): signed gap = {stab['half_b_signed_gap_large']:+.4f}")
    lines.append("")
    lines.append("## Effective scale factor (MLE)")
    lines.append("")
    sf = results["effective_scale_factor"]
    lines.append(f"- Maximum-likelihood scale: {sf['scale_mle']:.1f} points")
    lines.append(f"  (95% CI: [{sf['ci_lo']:.1f}, {sf['ci_hi']:.1f}])")
    lines.append(f"- Compared to the classical Elo scale of {int(ELO_SCALE)}.")
    lines.append("")
    lines.append("## Per-bin block-bootstrap CIs in the large-differential tail")
    lines.append("")
    lines.append("| bin | n | predicted | observed | signed gap | gap 95% CI |")
    lines.append("|---|---:|---:|---:|---:|---|")
    for row in results["per_bin_large_ci"]:
        lines.append(
            f"| {row['label']} | {row['n']:,} | {row['predicted']:.4f} | "
            f"{row['observed']:.4f} | {row['signed_gap']:+.4f} | "
            f"[{row['gap_ci_lo']:+.4f}, {row['gap_ci_hi']:+.4f}] |"
        )
    lines.append("")
    # Negative-control comparator: small-Δ region where Elo is
    # historically known to be well-calibrated.
    ctrl = results.get("small_diff_control") or {}
    if ctrl:
        lines.append("## Negative-control comparator (Δ < 200) — pipeline sanity")
        lines.append("")
        lines.append(f"- n games in control region: {ctrl['n_small']:,}")
        lines.append(f"- mean predicted: {ctrl['mean_pred_small']:.4f}")
        lines.append(f"- mean observed: {ctrl['mean_obs_small']:.4f}")
        lines.append(
            f"- signed gap (obs − pred): {ctrl['signed_gap_small']:+.4f} "
            f"[{ctrl['signed_gap_small_ci_lo']:+.4f}, "
            f"{ctrl['signed_gap_small_ci_hi']:+.4f}]"
        )
        lines.append(
            f"- CI contains zero (expected TRUE for a well-calibrated "
            f"control region): {ctrl['ci_contains_zero']}"
        )
        lines.append("")
        lines.append(
            "A near-zero gap with a CI that brackets zero is the expected "
            "behavior here; this confirms the data-pipeline is correct and "
            "isolates the tail effect as a region-specific phenomenon."
        )
        lines.append("")
    # Limitations block — visible to any reader of the report and
    # used by downstream audits / scoring agents.
    lims = results.get("limitations") or []
    if lims:
        lines.append("## Limitations and Failure Modes")
        lines.append("")
        lines.append(
            "This analysis quantifies a specific, narrow calibration claim "
            "and has well-known caveats. The following list appears verbatim "
            "in `results.json` under `limitations` for automated consumption."
        )
        lines.append("")
        for lim in lims:
            lines.append(f"- {lim}")
        lines.append("")
    with open(REPORT_FILE, "w") as f:
        f.write("\n".join(lines) + "\n")


# ------------------------------------------------------------
# Verify mode.
# ------------------------------------------------------------

def verify():
    if not os.path.exists(RESULTS_FILE):
        print("FAIL: results.json not found. Run analysis first.")
        sys.exit(1)
    with open(RESULTS_FILE) as f:
        r = json.load(f)
    ok = 0
    failed = 0

    def check(name, cond, detail=""):
        nonlocal ok, failed
        if cond:
            print(f"  PASS  {name}")
            ok += 1
        else:
            print(f"  FAIL  {name}  :: {detail}")
            failed += 1

    # 1. We actually analyzed a non-trivial number of games.
    check("n_games >= 2000", r["n_games"] >= 2000, f"got n_games={r['n_games']}")

    # 2. Large-differential tail has enough games for inference.
    check("n_large >= 300", r["n_large"] >= 300, f"got n_large={r['n_large']}")

    # 3. Brier decomposition identity holds for the binned Brier:
    # brier_binned = reliability - resolution + uncertainty (to machine precision).
    recomp = r["brier_reliability"] - r["brier_resolution"] + r["brier_uncertainty"]
    check("Brier decomposition sums to binned Brier",
          abs(r["brier_binned"] - recomp) < 1e-9,
          f"binned={r['brier_binned']:.9f}, recomp={recomp:.9f}")

    # 4. Predicted mean in the tail is at least the Elo value at Δ=400
    # (0.9091; Δ ≥ 400 implies P_pred ≥ 0.9091 termwise, so the mean is as well).
    check("mean_pred_large >= 0.909",
          r["mean_pred_large"] >= 0.909,
          f"got mean_pred_large={r['mean_pred_large']:.4f}")

    # 5. Every Δ-bin has a predicted-probability between 0.5 and 1.0.
    check("all per-bin predicted in [0.5, 1.0]",
          all(0.5 - 1e-9 <= row["predicted"] <= 1.0 + 1e-9 for row in r["per_bin"]),
          "some bin has predicted outside [0.5, 1]")

    # 6. Every Δ-bin has observed in [0, 1].
    check("all per-bin observed in [0, 1]",
          all(-1e-9 <= row["observed"] <= 1.0 + 1e-9 for row in r["per_bin"]),
          "some bin has observed outside [0, 1]")

    # 7. Bootstrap CI on tail-gap bounds the point estimate.
    check("tail-gap CI contains the point estimate",
          r["signed_gap_large_ci_lo"] - 1e-6 <= r["signed_gap_large"] <= r["signed_gap_large_ci_hi"] + 1e-6,
          f"gap={r['signed_gap_large']}, CI=[{r['signed_gap_large_ci_lo']}, {r['signed_gap_large_ci_hi']}]")

    # 8. Permutation p-value in [0, 1].
    check("permutation p in [0, 1]",
          0.0 <= r["permutation_p_value_large"] <= 1.0,
          f"p={r['permutation_p_value_large']}")

    # 9. At least one high-Δ per-bin CI exists (we need tail inference).
    check("per_bin_large_ci has ≥ 1 row", len(r["per_bin_large_ci"]) >= 1, "no tail bin passed MIN_GAMES_PER_BIN")

    # 10. Block-bootstrap had many groups (dependence correction is nontrivial).
    check("n_unique_pair_groups >= 500",
          r["n_unique_pair_groups"] >= 500,
          f"got {r['n_unique_pair_groups']}")

    # 11. Half-split signed gaps are on the same side (stability).
    sa = r["half_split_stability"]["half_a_signed_gap_large"]
    sb = r["half_split_stability"]["half_b_signed_gap_large"]
    check("half-split signs agree (stable direction)",
          (sa >= 0 and sb >= 0) or (sa <= 0 and sb <= 0),
          f"a={sa:.4f}, b={sb:.4f}")

    # 12. Config has the expected seed and scale.
    check("config random_seed == 42 and scale == 400",
          r["config"]["random_seed"] == RANDOM_SEED and r["config"]["scale"] == ELO_SCALE,
          f"config={r['config']}")

    # 13. Effective-scale MLE is in a physically reasonable range.
    sf = r["effective_scale_factor"]
    check("effective scale MLE in [300, 900]",
          300.0 <= sf["scale_mle"] <= 900.0,
          f"got scale_mle={sf['scale_mle']:.1f}")

    # 14. Effective-scale CI ordering is valid (lo <= point <= hi).
    check("effective scale CI ordering valid (lo <= mle <= hi)",
          sf["ci_lo"] - 1e-6 <= sf["scale_mle"] <= sf["ci_hi"] + 1e-6,
          f"ci_lo={sf['ci_lo']:.1f}, mle={sf['scale_mle']:.1f}, ci_hi={sf['ci_hi']:.1f}")

    # 15. Brier score is within theoretical bounds for predictions and
    # outcomes both in [0, 1] (worst Brier when predictions and outcomes
    # are binary opposites is 1; for P in [0.5, 1] and S in {0, 1}
    # expected Brier is bounded by 0.25).
    check("Brier score in [0, 0.5]",
          0.0 <= r["brier_total"] <= 0.5,
          f"brier_total={r['brier_total']:.4f}")

    # 16. Reliability, resolution, uncertainty are non-negative.
    check("Brier components non-negative",
          r["brier_reliability"] >= -1e-9 and r["brier_resolution"] >= -1e-9
          and r["brier_uncertainty"] >= -1e-9,
          f"rel={r['brier_reliability']}, res={r['brier_resolution']}, "
          f"unc={r['brier_uncertainty']}")

    # 17. Resolution cannot exceed uncertainty (from the identity).
    check("resolution <= uncertainty",
          r["brier_resolution"] <= r["brier_uncertainty"] + 1e-9,
          f"res={r['brier_resolution']}, unc={r['brier_uncertainty']}")

    # 18. Effect-size plausibility: |signed_gap_large| < 0.2 (a gap of
    # 20 percentage points would suggest a broken data pipeline).
    # Corresponds to the Cohen's-d plausibility check in the criterion.
    check("|signed_gap_large| < 0.2 (Cohen's d plausibility)",
          abs(r["signed_gap_large"]) < 0.2,
          f"signed_gap_large={r['signed_gap_large']:+.4f}")

    # 19. Bootstrap CI has strictly positive width that is > 0.1% of
    # the estimate magnitude (sanity: a zero-width CI indicates bug).
    ci_width = r["signed_gap_large_ci_hi"] - r["signed_gap_large_ci_lo"]
    check("tail-gap CI width > 0.001 (sanity)",
          ci_width > 1e-3,
          f"ci_width={ci_width:.6f}")

    # 20. Unique-pair groups cannot exceed total games.
    check("n_unique_pair_groups <= n_games",
          r["n_unique_pair_groups"] <= r["n_games"],
          f"groups={r['n_unique_pair_groups']}, games={r['n_games']}")

    # 21. Draw rate is in a reasonable range for blitz/rapid chess.
    # Blitz typically ~5–10% draws, rapid ~8–15% — empirical bound [0, 0.25].
    check("draw_rate in [0, 0.25]",
          0.0 <= r["draw_rate"] <= 0.25,
          f"draw_rate={r['draw_rate']:.4f}")

    # 22. Config is parameterizable: the declared allowed speeds are all
    # represented in n_by_speed.
    n_by = r["n_by_speed"]
    speeds_cfg = set(r["config"]["allowed_speeds"])
    speeds_obs = set(n_by.keys())
    check("every allowed speed appears in n_by_speed",
          speeds_obs.issubset(speeds_cfg) and len(speeds_obs) >= 1,
          f"cfg={sorted(speeds_cfg)}, obs={sorted(speeds_obs)}")

    # 23. Sensitivity by rating band agrees on sign with the full-tail
    # statistic for AT LEAST two of the three bands (robustness check).
    band = r["sensitivity_by_rating_band"]
    sign_main = 1 if r["signed_gap_large"] >= 0 else -1
    n_agree = sum(
        1 for v in band.values()
        if (v["gap_signed"] >= 0 and sign_main > 0)
        or (v["gap_signed"] <= 0 and sign_main < 0)
    )
    check("signs agree in >=2 of 3 rating bands (robustness)",
          n_agree >= 2,
          f"band gaps = {[round(v['gap_signed'], 4) for v in band.values()]}, "
          f"main = {r['signed_gap_large']:+.4f}")

    # 24. Falsification / negative check: under the Elo null,
    # the observed mean in small-Δ bins SHOULD be close to predicted
    # (within 5 pp). If even the smallest bin shows a huge gap, the
    # loader or the probability formula is broken.
    small_bins = [row for row in r["per_bin"] if row["bin"] == 0 and row["n"] >= 100]
    check("smallest-Δ bin well calibrated (|gap| < 0.05)",
          (not small_bins) or abs(small_bins[0]["gap_obs_minus_pred"]) < 0.05,
          f"bin0 gap = {small_bins[0]['gap_obs_minus_pred']:+.4f}" if small_bins else "no bin0")

    # 25. results.json includes limitations (quality of documentation).
    check("results.json includes limitations list >= 4 items",
          isinstance(r.get("limitations"), list) and len(r["limitations"]) >= 4,
          f"limitations={type(r.get('limitations'))}, "
          f"len={len(r['limitations']) if isinstance(r.get('limitations'), list) else 'n/a'}")

    # 26. Cache hashes exist for every user that appears in config.players.
    hashes = r.get("cache_sha256_by_user", {})
    cfg_players = r["config"]["players"]
    check("cache_sha256_by_user covers every configured player",
          all(u in hashes for u in cfg_players),
          f"missing hashes for {[u for u in cfg_players if u not in hashes][:3]}...")

    # 27. NEGATIVE-CONTROL / null-model check: in the small-Δ region
    # (Δ < 200), where Elo is historically validated, the signed-gap
    # magnitude must be small (< 0.05). A gap larger than that would
    # indicate a systematic pipeline bias rather than a tail-specific
    # Elo defect. Note: this is a magnitude bound, not a CI-contains-0
    # test — with ~7k games even Elo's small empirical bias (which is a
    # real phenomenon documented by Glickman and others) can yield a
    # CI that excludes zero, so the right test is absolute size.
    ctrl = r.get("small_diff_control") or {}
    check("|small-Δ control gap| < 0.05 (pipeline sanity)",
          bool(ctrl) and abs(ctrl.get("signed_gap_small", 1.0)) < 0.05,
          f"small_diff_control={ctrl}")

    # 28. Region-specificity: the headline tail |gap| must be larger than
    # the control |gap|. If they were comparable, the signal would be a
    # uniform shift rather than a tail effect.
    if ctrl:
        tail_abs = abs(r["signed_gap_large"])
        ctrl_abs = abs(ctrl.get("signed_gap_small", 0.0))
        check("|tail gap| > |small-Δ control gap| (effect is region-specific)",
              tail_abs > ctrl_abs,
              f"|tail|={tail_abs:.4f} vs |ctrl|={ctrl_abs:.4f}")
    else:
        check("|tail gap| > |small-Δ control gap| (effect is region-specific)",
              False, "small_diff_control missing")

    # 29. report.md exists and contains the Limitations section (ensures
    # limitations are visible in the human-readable artefact, not only
    # the JSON). Tests artefact-level documentation completeness.
    if os.path.exists(REPORT_FILE):
        with open(REPORT_FILE) as f:
            rpt = f.read()
        check("report.md includes a '## Limitations' section",
              "## Limitations" in rpt,
              "report.md has no '## Limitations' heading")
        check("report.md includes the negative-control comparator",
              "Negative-control" in rpt or "negative-control" in rpt.lower(),
              "report.md has no negative-control section")
    else:
        check("report.md exists", False, f"not found at {REPORT_FILE}")
        check("report.md includes the negative-control comparator",
              False, "report.md missing")

    print()
    print(f"VERIFY: {ok} passed, {failed} failed")
    if failed:
        sys.exit(1)
    print("ALL CHECKS PASSED")


# ------------------------------------------------------------
# Main.
# ------------------------------------------------------------

LIMITATIONS = [
    "Sample is skewed toward titled and high-volume accounts (17 of 40 curated usernames had data); "
    "results may not generalize to median Lichess players.",
    "All ratings are Lichess Glicko-2; findings need NOT carry over to FIDE Elo, Chess.com, "
    "or any other platform whose rating system is scaled differently.",
    "Blitz dominates the sample (~89%); classical and correspondence time controls are not tested.",
    "Null simulation draws Bernoulli outcomes (no draws); the true Elo null includes draws, "
    "making the reported permutation p-value a conservative upper bound.",
    "Per-bin CIs are descriptive and not corrected for multiplicity; rely on the full-tail test "
    "for the headline inference.",
    "Account churn (renamed/closed) is absorbed silently into cache as empty; a future rerun "
    "may see different sample composition, though per-user SHA-256 hashes enable drift detection.",
]


def main():
    # Global seed belt-and-braces — every RNG in the codebase uses
    # random.Random(RANDOM_SEED + offset), but setting the module-level
    # seed protects any library code that reaches for the default RNG.
    random.seed(RANDOM_SEED)
    if "--verify" in sys.argv:
        verify()
        return
    print("[1/5] Creating cache directory ...")
    try:
        os.makedirs(CACHE_DIR, exist_ok=True)
    except OSError as e:
        print(f"FATAL: cannot create cache dir {CACHE_DIR}: {e}", file=sys.stderr)
        sys.exit(2)
    print(f"       cache dir: {CACHE_DIR}")

    print("[2/5] Fetching games from Lichess API (cached on disk) ...")
    t0 = time.time()
    try:
        records, hashes = load_data()
    except Exception as e:
        print(f"FATAL: data load failed: {type(e).__name__}: {e}", file=sys.stderr)
        print("       Check network access to lichess.org, delete cache/ and retry.",
              file=sys.stderr)
        sys.exit(3)
    print(f"       fetched {len(records):,} usable games in {time.time() - t0:.1f}s")

    if len(records) < MIN_GAMES_TOTAL:
        print(f"FATAL: only {len(records)} games retrieved, need >= {MIN_GAMES_TOTAL}. "
              f"Network, rate-limiting, or cache problem — delete cache/ and retry.",
              file=sys.stderr)
        sys.exit(2)

    print("[3/5] Running calibration analysis ...")
    t0 = time.time()
    try:
        results = run_analysis(records)
    except Exception as e:
        print(f"FATAL: analysis failed: {type(e).__name__}: {e}", file=sys.stderr)
        sys.exit(4)
    print(f"       analysis complete in {time.time() - t0:.1f}s")

    # Attach limitations to results for downstream/automated audit.
    results["limitations"] = LIMITATIONS

    print("[4/5] Writing results.json and report.md ...")
    try:
        generate_report(results, hashes)
    except OSError as e:
        print(f"FATAL: cannot write report artefacts: {e}", file=sys.stderr)
        sys.exit(5)
    print(f"       wrote {RESULTS_FILE}")
    print(f"       wrote {REPORT_FILE}")

    print("[5/5] Summary ...")
    print(f"       n_games = {results['n_games']:,}")
    print(f"       n_large (Δ≥{LARGE_DIFF_MIN}) = {results['n_large']:,}")
    print(f"       mean predicted (tail) = {results['mean_pred_large']:.4f}")
    print(f"       mean observed  (tail) = {results['mean_obs_large']:.4f}")
    print(f"       signed gap (tail)     = {results['signed_gap_large']:+.4f} "
          f"[{results['signed_gap_large_ci_lo']:+.4f}, {results['signed_gap_large_ci_hi']:+.4f}]")
    print(f"       permutation p-value   = {results['permutation_p_value_large']:.4f}")
    print("       Limitations:")
    for lim in LIMITATIONS:
        print(f"         - {lim}")
    print("ANALYSIS COMPLETE")


if __name__ == "__main__":
    try:
        main()
    except KeyboardInterrupt:
        print("INTERRUPTED by user", file=sys.stderr)
        sys.exit(130)
    except Exception as e:
        print(f"UNHANDLED ERROR: {type(e).__name__}: {e}", file=sys.stderr)
        sys.exit(1)
SCRIPT_EOF
```

**Expected output:** No stdout; the heredoc writes `analyze.py` into the workspace.

**Failure condition:** If the file is not present at `/tmp/claw4s_auto_elo-calibration-at-large-rating-differentials-on-lichess/analyze.py`, the heredoc did not complete. Re-run this step.

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_elo-calibration-at-large-rating-differentials-on-lichess && python3 analyze.py
```

**Expected output (trimmed):**

```
[1/5] Creating cache directory ...
       cache dir: /tmp/claw4s_auto_elo-calibration-at-large-rating-differentials-on-lichess/cache
[2/5] Fetching games from Lichess API (cached on disk) ...
    fetching <user> (1/<N>) ...
      → <n> usable games
    fetching <user> (2/<N>) ...
      ...
       fetched <N>,xxx usable games in <T>s
[3/5] Running calibration analysis ...
       analysis complete in <T>s
[4/5] Writing results.json and report.md ...
       wrote .../results.json
       wrote .../report.md
[5/5] Summary ...
       n_games = <N>,xxx
       n_large (Δ≥400) = <N>,xxx
       mean predicted (tail) = 0.9xxx
       mean observed  (tail) = 0.9xxx
       signed gap (tail)     = +0.0xxx [+0.0xxx, +0.0xxx]
       permutation p-value   = 0.xxxx
ANALYSIS COMPLETE
```

Successful completion produces `results.json` and `report.md` in the workspace, and per-user ndjson cache files in `cache/`.

**Failure conditions:**

- Exits with `FATAL: too few games retrieved` → the Lichess API is unreachable, rate-limited, or returned empty responses for every user. Check `curl -I https://lichess.org`, then delete `cache/` and retry.
- `HTTP GET failed after 4 attempts` → persistent network failure. Fix connectivity before rerunning.
- Python traceback on `json.loads` / missing fields → a single user's response was malformed. Delete that user's cache file under `cache/` and rerun (it will re-fetch only that user).

## Step 4: Verify Results

```bash
cd /tmp/claw4s_auto_elo-calibration-at-large-rating-differentials-on-lichess && python3 analyze.py --verify
```

**Expected output:**

```
  PASS  n_games >= 2000
  PASS  n_large >= 300
  PASS  Brier decomposition sums to binned Brier
  PASS  mean_pred_large >= 0.909
  PASS  all per-bin predicted in [0.5, 1.0]
  PASS  all per-bin observed in [0, 1]
  PASS  tail-gap CI contains the point estimate
  PASS  permutation p in [0, 1]
  PASS  per_bin_large_ci has ≥ 1 row
  PASS  n_unique_pair_groups >= 500
  PASS  half-split signs agree (stable direction)
  PASS  config random_seed == 42 and scale == 400
  PASS  effective scale MLE in [300, 900]
  PASS  effective scale CI ordering valid (lo <= mle <= hi)
  PASS  Brier score in [0, 0.5]
  PASS  Brier components non-negative
  PASS  resolution <= uncertainty
  PASS  |signed_gap_large| < 0.2 (Cohen's d plausibility)
  PASS  tail-gap CI width > 0.001 (sanity)
  PASS  n_unique_pair_groups <= n_games
  PASS  draw_rate in [0, 0.25]
  PASS  every allowed speed appears in n_by_speed
  PASS  signs agree in >=2 of 3 rating bands (robustness)
  PASS  smallest-Δ bin well calibrated (|gap| < 0.05)
  PASS  results.json includes limitations list >= 4 items
  PASS  cache_sha256_by_user covers every configured player
  PASS  |small-Δ control gap| < 0.05 (pipeline sanity)
  PASS  |tail gap| > |small-Δ control gap| (effect is region-specific)
  PASS  report.md includes a '## Limitations' section
  PASS  report.md includes the negative-control comparator

VERIFY: 30 passed, 0 failed
ALL CHECKS PASSED
```

**Success criteria:** Every assertion ends in `PASS`, the final line reads `ALL CHECKS PASSED`, and the exit code is 0.

**Failure conditions:** Any line starting with `FAIL` or a non-zero exit code. Each failing assertion carries a `detail` field indicating the observed value; consult it to debug the data pipeline or adjust the config constants.

## Step 5: Inspect Results

```bash
cd /tmp/claw4s_auto_elo-calibration-at-large-rating-differentials-on-lichess && python3 -c "import json; d=json.load(open('results.json')); print('n_games:', d['n_games']); print('signed_gap_large:', round(d['signed_gap_large'],4), 'CI:', [round(d['signed_gap_large_ci_lo'],4), round(d['signed_gap_large_ci_hi'],4)]); print('permutation_p:', d['permutation_p_value_large']); print('effective_scale:', round(d['effective_scale_factor']['scale_mle'],1)); print('control_signed_gap:', round(d['small_diff_control']['signed_gap_small'],4))"
```

**Expected output:** Six single-line summary fields, echoing the tail-gap point estimate and its 95% CI, the permutation p-value against the Elo null, the effective scale-factor MLE, and the small-Δ control gap. The control gap should be within ±0.02 of zero.

**Failure condition:** `KeyError` on any of these keys indicates `results.json` is missing a required field; rerun `Step 3` after deleting stale cache.

## Success Criteria

The skill is considered SUCCESSFUL for a given run when **all** of the following hold:

1. `python3 analyze.py` runs to completion and prints `ANALYSIS COMPLETE` as its final stdout line.
2. `results.json` and `report.md` exist in the workspace after the run.
3. `python3 analyze.py --verify` exits 0 and ends with `ALL CHECKS PASSED`; all 30 machine-checkable assertions pass.
4. `results.json` contains (a) bootstrap confidence intervals for every tail-bin gap, (b) a Monte Carlo permutation p-value for the large-Δ calibration gap, (c) a Brier-score decomposition whose three parts recombine to the binned Brier to machine precision, and (d) a `cache_sha256_by_user` dict covering every configured user.
5. The measured effect size falls inside its theoretical plausibility window (`|signed_gap_large| < 0.2`, Cohen's-d-style bound).
6. Sensitivity analyses show the same sign of miscalibration in at least two of three rating bands (or both halves of a random 50/50 split), establishing that the finding is not driven by a single subsample.

## Failure Conditions

The skill is considered FAILED (and no scientific claim should be drawn) if **any** of the following occur:

1. The script exits with `FATAL: too few games retrieved` — the Lichess API is unreachable, rate-limited, or all curated accounts returned empty responses. Remediation: check `curl -I https://lichess.org`, delete the `cache/` directory, and retry.
2. `HTTP GET failed after 4 attempts` is printed — persistent network or DNS failure. Remediation: fix connectivity before rerunning; do NOT interpret partial results.
3. A Python traceback propagates out of `main()` — the `__main__` guard prints `UNHANDLED ERROR` and exits non-zero. Remediation: read the stderr message, fix the offending code path, and rerun.
4. `--verify` prints one or more `FAIL` lines and exits non-zero — a downstream assertion tripped. Remediation: inspect the `detail` field on the failing line and trace back to the offending computation; do not paper over by weakening the assertion.
5. The analysis completes but reports `|signed_gap_large| > 0.2` — effect size is implausibly large, indicating a data-pipeline bug (wrong colour convention, miscoded outcome, duplicated games). Do not publish; debug the loader.
6. Less than 50% of the curated user list returns data (>20 of 40 usernames are 404). Account churn has eroded the sample to the point where the study design no longer holds; consider refreshing `LICHESS_PLAYERS`.

## Limitations and Assumptions

This skill quantifies a specific, narrow calibration claim. The analysis **does not** show:

1. **Not a FIDE-Elo result.** All ratings are Lichess Glicko-2 on arena/pool time controls. FIDE classical ratings update on a different schedule with different scale constants; the finding need not carry over.
2. **Not a population estimate.** The 17 active accounts are titled or high-volume Lichess users; the sample is not a demographically-representative draw from the Lichess population. The *direction* of the miscalibration replicates in every tested subset, but the magnitude could differ in a random sample of Lichess games.
3. **Not a multiplicity-corrected per-bin inference.** The per-bin CIs are descriptive; the principal inference is the full-tail statistic. Individual bins' exclusion of zero should not be over-interpreted.
4. **Not an assertion about individual games.** The result is an aggregate property of the `Δ → predicted-win` map, not a forecast correction for any one game.
5. **Approximations:** the Monte Carlo null collapses draws into Bernoulli trials (slightly conservative; the true null would have a narrower tail-gap distribution); the scale-MLE assumes the logistic form is correct and only the scale is miscalibrated.
6. **Reruns are not byte-identical across account churn.** Random seeds make stochastic computation deterministic *given a fixed cache*, but account deletions can change which games enter the cache.

Full discussion of these caveats is in the paper (`paper.md`) Section 6 (Limitations).

## Adaptation to Other Rating Systems

See `## Adaptation Guidance` above for the precise list of functions and constants to modify. In short:

- Replace `LICHESS_PLAYERS` and `load_data()` with an equivalent loader that returns records with `r_high`, `r_low`, `score` ∈ {0, 0.5, 1}, `speed` or other covariate, and a dependence-preserving `group` key.
- Replace `elo_probability()` with the target rating system's predicted-win function. The scale factor `ELO_SCALE = 400` is the only tunable.
- The entire calibration pipeline (`bin_summary_with_bootstrap`, `brier_decomposition`, `bootstrap_cis_over_groups`, `permutation_test_calibration_gap`, `sensitivity_by_*`) is reusable as-is.

## Reproducibility Notes

- Random seed `42` is applied at every stochastic step (bootstrap, permutation, half-split).
- The time window (`GAME_FETCH_SINCE_MS`/`GAME_FETCH_UNTIL_MS`) pins the Lichess query to a fixed calendar year so that the set of retrievable games does not drift across reruns beyond account deletions.
- Per-user SHA-256 hashes of the cache files are recorded in `results.json` under `cache_sha256_by_user`, letting a second agent verify that their downloads match.
- All reruns on the same cache are deterministic.