Measure each build's HBAI distribution and single ages against its base

.github/workflows/pull_request.yaml

@@ -109,12 +109,26 @@ jobs:
         with:
           name: fit_summary.json
           path: fit_summary.json
+      # After the tests, and whatever they did: the measurement downloads the
+      # ONS projection, and an ONS outage must not skip the test suite.
+      - name: Measure the income distribution and single ages
+        if: ${{ !cancelled() }}
+        # A failed measurement fails the comparison (missing report), not
+        # the test job, which stays about the tests.
+        continue-on-error: true
+        run: uv run --frozen python -m policyengine_uk_data.utils.distribution_validation measure policyengine_uk_data/storage/enhanced_frs_2024_25.h5 --out distribution_report.json
+      - name: Save the distribution report
+        if: ${{ !cancelled() }}
+        uses: actions/upload-artifact@v7
+        with:
+          name: distribution_report.json
+          path: distribution_report.json
   # The no-worsening check (W3). The base is built the same way as the head:
   # TESTING, one clone and the head's lockfile, so both builds run one model
   # and environment and the comparison measures the data change. As in the
-  # test job, only the fit summary leaves it. No dataset or calibration log
-  # is uploaded, and nothing is cached: setup-uv restores uv's package cache
-  # and saves none.
+  # test job, only aggregate reports leave it: the fit summary and the
+  # distribution report. No dataset or calibration log is uploaded, and
+  # nothing is cached: setup-uv restores uv's package cache and saves none.
   base-build:
     name: Build the base for comparison
     runs-on: ubuntu-latest
@@ -164,6 +178,19 @@ jobs:
         with:
           name: base_fit_summary.json
           path: base_fit_summary.json
+      - name: Measure the base's income distribution and single ages
+        working-directory: base
+        # The head's measurement code, on the base's build.
+        run: |
+          mkdir -p policyengine_uk_data/targets/reference_data
+          cp ../head/policyengine_uk_data/utils/distribution_validation.py policyengine_uk_data/utils/
+          cp ../head/policyengine_uk_data/targets/reference_data/hbai_fye_2025.yaml policyengine_uk_data/targets/reference_data/
+          uv run --frozen python -m policyengine_uk_data.utils.distribution_validation measure policyengine_uk_data/storage/enhanced_frs_2024_25.h5 --out ../base_distribution_report.json
+      - name: Save the base distribution report
+        uses: actions/upload-artifact@v7
+        with:
+          name: base_distribution_report.json
+          path: base_distribution_report.json
   fit-comparison:
     name: Calibration fit against the base
     needs: [test, base-build]
@@ -196,6 +223,13 @@ jobs:
           pattern: "*fit_summary.json"
           path: logs
           merge-multiple: true
+      - name: Download the distribution reports
+        uses: actions/download-artifact@v8
+        continue-on-error: true
+        with:
+          pattern: "*distribution_report.json"
+          path: reports
+          merge-multiple: true
       - name: Compare the final fits
         id: compare
         env:
@@ -218,6 +252,18 @@ jobs:
             set -e
             { echo "$marker"; cat body.md 2>/dev/null || echo '## Calibration fit against the base build: FAIL'; echo; echo "[Workflow run]($RUN_URL)"; } > fit.md
           fi
+          if [ -f reports/base_distribution_report.json ] && [ -f reports/distribution_report.json ]; then
+            set +e
+            uv run --frozen python -m policyengine_uk_data.utils.distribution_validation compare \
+              reports/base_distribution_report.json reports/distribution_report.json --markdown distribution.md
+            distribution_status=$?
+            set -e
+            { echo; cat distribution.md 2>/dev/null || echo '## Income distribution and single ages: FAIL (the comparison did not finish)'; } >> fit.md
+            if [ "$distribution_status" != 0 ]; then echo "status=1" >> "$GITHUB_OUTPUT"; fi
+          else
+            { echo; echo '## Income distribution and single ages: FAIL'; echo; echo "The report is missing for one build, so the distributions could not be compared."; } >> fit.md
+            echo "status=1" >> "$GITHUB_OUTPUT"
+          fi
           cat fit.md >> "$GITHUB_STEP_SUMMARY"
       - name: Post or update the pull request comment
         if: ${{ !cancelled() && github.event.pull_request.head.repo.full_name == github.repository }}
@@ -246,5 +292,5 @@ jobs:
       - name: Fail on a fit regression
         if: steps.compare.outputs.status != '0'
         run: |
-          echo "::error::The head build's calibration fit is worse than the base's; see the job summary."
+          echo "::error::The calibration fit or distribution comparison failed; see the job summary."
           exit 1

changelog.d/distribution-checks.added.md

@@ -0,0 +1 @@
+Measure each pull request build's HBAI income distribution (medians, means, Gini and relative low income by group against DWP's HBAI FYE 2025 estimates and intervals) and single-age population (against the ONS 2024-based projection), and report them against the base build in the pull request comment. They flag but do not fail yet.

policyengine_uk_data/targets/reference_data/hbai_fye_2025.yaml

@@ -0,0 +1,27 @@
+# Households below average income, FYE 2025 (DWP, 26 March 2026), from the
+# published data tables:
+# https://assets.publishing.service.gov.uk/media/69c2902a7e02b81c0d1c7762/hbai-data_tables-2025.zip
+# uncertainty-hbai-2024-25-tables.ods, tables 8a (medians, means, Gini) and
+# 8b (low income), read on 15 October 2026. Estimates and 95% confidence
+# intervals for the United Kingdom, FRS 2024/25, which policyengine-uk models
+# as period 2024. Incomes are equivalised, in £ a week at 2024/25 prices;
+# rates and Gini coefficients are per cent.
+schema_version: 1
+source: DWP Households below average income FYE 2025
+period: 2024
+measures:
+  median_bhc: {estimate: 719.481062, low: 703.567322, high: 735.575892, table: 8a}
+  median_ahc: {estimate: 623.144059, low: 605.209749, high: 637.071794, table: 8a}
+  mean_bhc: {estimate: 854.979494, low: 836.367495, high: 870.351587, table: 8a}
+  mean_ahc: {estimate: 749.82155, low: 732.329379, high: 764.496483, table: 8a}
+  gini_bhc: {estimate: 32.736632, low: 31.889881, high: 33.978308, table: 8a}
+  gini_ahc: {estimate: 36.780956, low: 35.897706, high: 38.293054, table: 8a}
+  # Percentage of individuals below 60% of the relative median.
+  relative_low_income_bhc_all: {estimate: 15.911572, low: 14.446526, high: 17.637657, table: 8b}
+  relative_low_income_bhc_children: {estimate: 20.519181, low: 17.356404, high: 23.955585, table: 8b}
+  relative_low_income_bhc_working_age: {estimate: 14.281508, low: 12.908347, high: 16.018814, table: 8b}
+  relative_low_income_bhc_pensioners: {estimate: 15.895071, low: 13.690455, high: 17.991754, table: 8b}
+  relative_low_income_ahc_all: {estimate: 19.640844, low: 17.9939, high: 21.271606, table: 8b}
+  relative_low_income_ahc_children: {estimate: 27.41076, low: 23.49369, high: 30.644085, table: 8b}
+  relative_low_income_ahc_working_age: {estimate: 18.579734, low: 17.012601, high: 20.343162, table: 8b}
+  relative_low_income_ahc_pensioners: {estimate: 13.866371, low: 11.618244, high: 15.898207, table: 8b}

policyengine_uk_data/tests/test_distribution_validation.py

@@ -0,0 +1,245 @@
+"""Single-age and HBAI distribution checks (utils/distribution_validation.py)."""
+
+import numpy as np
+import pytest
+from hypothesis import given, settings
+from hypothesis import strategies as st
+
+from policyengine_uk_data.utils.distribution_validation import (
+    SCHEMA,
+    TOP_AGE_BAND,
+    compare,
+    hbai_measures,
+    load_hbai_reference,
+    render,
+    single_ages,
+    summarise_single_ages,
+    weighted_gini,
+    weighted_median,
+)
+
+samples = st.lists(
+    st.tuples(st.floats(-1e4, 1e5, allow_nan=False), st.integers(1, 5)),
+    min_size=1,
+    max_size=40,
+)
+
+
+def _expand(pairs):
+    return np.repeat([v for v, _ in pairs], [w for _, w in pairs]).astype(float)
+
+
+@settings(max_examples=300, deadline=None)
+@given(pairs=samples)
+def test_weighted_median_equals_the_median_of_the_expanded_sample(pairs):
+    """Differential: integer weights are repeated records."""
+    values = np.array([v for v, _ in pairs])
+    weights = np.array([w for _, w in pairs], float)
+    expanded = np.sort(_expand(pairs))
+    # The lower median of the expanded sample.
+    lower = expanded[(len(expanded) - 1) // 2]
+    upper = expanded[len(expanded) // 2]
+    assert lower <= weighted_median(values, weights) <= upper
+    # Scaling every weight leaves it unchanged.
+    assert weighted_median(values, weights * 7.5) == weighted_median(values, weights)
+
+
+@settings(max_examples=300, deadline=None)
+@given(
+    pairs=st.lists(
+        st.tuples(st.floats(0.01, 1e5), st.integers(1, 5)), min_size=2, max_size=40
+    )
+)
+def test_weighted_gini_equals_the_mean_difference_gini_of_the_expanded_sample(pairs):
+    expanded = _expand(pairs)
+    n = len(expanded)
+    mean_difference = np.abs(expanded[:, None] - expanded[None, :]).sum() / (n * n)
+    reference = mean_difference / (2 * expanded.mean())
+    values = np.array([v for v, _ in pairs])
+    weights = np.array([w for _, w in pairs], float)
+    assert weighted_gini(values, weights) == pytest.approx(reference, abs=1e-9)
+
+
+def test_gini_of_equal_incomes_is_zero():
+    assert weighted_gini([5.0, 5.0, 5.0], [1, 2, 3]) == pytest.approx(0)
+
+
+def test_hbai_measures_on_a_known_population():
+    income = np.array([100.0, 200.0, 300.0, 400.0, 1000.0])
+    weights = np.ones(5)
+    child = np.array([1, 0, 0, 0, 0], bool)
+    pensioner = np.array([0, 0, 0, 0, 1], bool)
+    m = hbai_measures(income, income - 50, weights, child, pensioner)
+    assert m["median_bhc"] == 300.0
+    assert m["mean_bhc"] == 400.0
+    # Below 60% of 300: 100 only, a child.
+    assert m["relative_low_income_bhc_all"] == pytest.approx(20.0)
+    assert m["relative_low_income_bhc_children"] == pytest.approx(100.0)
+    assert m["relative_low_income_bhc_working_age"] == pytest.approx(0.0)
+    assert m["relative_low_income_bhc_pensioners"] == pytest.approx(0.0)
+
+
+def test_the_hbai_reference_is_the_published_table():
+    """Spot values from DWP HBAI FYE 2025, tables 8a and 8b."""
+    ref = load_hbai_reference()
+    assert ref["period"] == 2024
+    m = ref["measures"]
+    assert m["median_bhc"]["estimate"] == pytest.approx(719.481062)
+    assert m["median_ahc"]["estimate"] == pytest.approx(623.144059)
+    assert m["relative_low_income_ahc_children"]["estimate"] == pytest.approx(27.41076)
+    assert len(m) == 14
+
+
+def test_single_ages_band_the_oldest_and_split_by_sex():
+    age = np.array([0, 0, 79, 80, 85, 90])
+    male = np.array([1, 0, 1, 1, 0, 1], bool)
+    weights = np.array([2.0, 3.0, 1.0, 1.0, 4.0, 5.0])
+    reference = {
+        (s, a): 1.0 for s in ("male", "female") for a in range(TOP_AGE_BAND + 1)
+    }
+    rows = single_ages(age, male, weights, reference)
+    assert rows["male"]["0"]["value"] == 2.0 and rows["female"]["0"]["value"] == 3.0
+    assert rows["male"]["80"]["value"] == 6.0 and rows["female"]["80"]["value"] == 4.0
+    summary = summarise_single_ages(rows)
+    assert summary["error_80_plus"] == {"male": 5.0, "female": 3.0}
+
+
+def _report(**hbai_values):
+    ref = load_hbai_reference()["measures"]
+    hbai = {}
+    for name, r in ref.items():
+        value = hbai_values.get(name, r["estimate"])
+        hbai[name] = {
+            "value": value,
+            "reference": r["estimate"],
+            "low": r["low"],
+            "high": r["high"],
+            "inside": r["low"] <= value <= r["high"],
+        }
+    single = {
+        "mean_abs_error_0_79": 0.2,
+        "max_abs_error_0_79": 0.9,
+        "max_abs_error_15_19": 0.6,
+        "error_80_plus": {"male": 0.0, "female": 0.0},
+    }
+    return {
+        "schema": SCHEMA,
+        "periods": {"hbai": 2024, "single_age": 2025},
+        "hbai": hbai,
+        "summary": {"single_age": single},
+    }
+
+
+def test_compare_flags_a_measure_that_leaves_its_interval_but_does_not_fail():
+    result = compare(_report(), _report(median_bhc=760.0))
+    assert result["status"] == "pass"
+    assert result["flags"] == [
+        "HBAI median_bhc leaves its 95% interval: 719.48 -> 760.00 (HBAI 719.48, 703.57 to 735.58)."
+    ]
+    text = render(result)
+    assert "760.00 (outside)" in text and "do not fail yet" in text
+
+
+def test_compare_with_itself_passes_with_no_flags():
+    result = compare(_report(), _report())
+    assert result["status"] == "pass" and result["flags"] == []
+    assert all(r["further_from_reference"] == 0 for r in result["rows"])
+
+
+@pytest.mark.parametrize(
+    "change, message",
+    [
+        (lambda r: r.update(schema="other"), "malformed"),
+        (lambda r: r["periods"].update(hbai=2025), "different periods"),
+        (lambda r: r["hbai"]["median_bhc"].update(value=float("nan")), "not finite"),
+    ],
+)
+def test_compare_fails_on_a_broken_or_mismatched_report(change, message):
+    head = _report()
+    change(head)
+    result = compare(_report(), head)
+    assert result["status"] == "fail"
+    assert any(message in f for f in result["failures"])
+
+
+def test_negative_incomes_before_housing_costs_are_reset_to_zero_as_hbai_does():
+    """Sol review of #600: HBAI resets negative BHC to zero and derives AHC from
+    the adjusted figure, so AHC can still be negative."""
+    from policyengine_uk_data.utils.distribution_validation import (
+        WEEKS_PER_YEAR,
+        hbai_incomes,
+    )
+
+    weekly_bhc = np.array([-100.0, 100.0, 200.0, 300.0, 400.0])
+    annual_bhc = weekly_bhc * WEEKS_PER_YEAR
+    annual_ahc = (weekly_bhc - 50) * WEEKS_PER_YEAR
+    bhc, ahc = hbai_incomes(annual_bhc, annual_ahc, np.ones(5), np.ones(5))
+    np.testing.assert_allclose(bhc, [0, 100, 200, 300, 400])
+    np.testing.assert_allclose(ahc, [-50, 50, 150, 250, 350])
+    m = hbai_measures(bhc, ahc, np.ones(5), np.zeros(5, bool), np.zeros(5, bool))
+    assert m["mean_bhc"] == pytest.approx(200) and m["mean_ahc"] == pytest.approx(150)
+    assert m["gini_bhc"] == pytest.approx(40.0)
+    assert m["gini_ahc"] == pytest.approx(53.333333, rel=1e-6)
+    # No children or pensioners: those groups have no rate, not a NaN.
+    assert m["relative_low_income_bhc_children"] is None
+    assert m["relative_low_income_ahc_pensioners"] is None
+    # Median 200, line 120: incomes 0 and 100 are below it.
+    assert m["relative_low_income_bhc_working_age"] == pytest.approx(40.0)
+
+
+def test_an_incomplete_or_non_finite_report_fails():
+    """Sol review of #600: every reference measure and age summary must be there."""
+    head = _report()
+    del head["hbai"]["median_bhc"]
+    result = compare(_report(), head)
+    assert result["status"] == "fail"
+    assert result["failures"] == ["The head report has no median_bhc."]
+    head = _report()
+    head["summary"]["single_age"]["mean_abs_error_0_79"] = float("nan")
+    result = compare(_report(), head)
+    assert result["status"] == "fail" and "not finite" in result["failures"][0]
+    head = _report()
+    name = next(iter(head["hbai"]))
+    head["hbai"][name]["value"] = None  # an empty group
+    result = compare(_report(), head)
+    assert result["status"] == "fail"
+    assert f"HBAI {name} is missing or not finite on one side." in result["failures"]
+
+
+def test_the_workflow_fails_when_a_distribution_report_is_missing():
+    from pathlib import Path
+
+    text = (
+        Path(__file__).parents[2] / ".github" / "workflows" / "pull_request.yaml"
+    ).read_text()
+    block = text[text.index("The report is missing for one build") :]
+    assert 'echo "status=1" >> "$GITHUB_OUTPUT"' in block.split("\n          fi", 1)[0]
+    assert "## Income distribution and single ages: FAIL" in text
+
+
+def test_the_measurement_runs_after_the_tests_and_cannot_skip_them():
+    """Delta reviews of #600: an ONS outage must not skip the tests or the
+    fit summary, nor turn the test job red."""
+    from pathlib import Path
+
+    import yaml
+
+    workflow = Path(__file__).parents[2] / ".github" / "workflows" / "pull_request.yaml"
+    steps = yaml.safe_load(workflow.read_text())["jobs"]["test"]["steps"]
+    names = [step.get("name") for step in steps]
+    measure = names.index("Measure the income distribution and single ages")
+    assert names.index("Run tests") < measure
+    assert names.index("Save the fit summary") < measure
+    assert steps[measure]["if"] == "${{ !cancelled() }}"
+    assert steps[measure]["continue-on-error"] is True
+    # The fit summary runs after the tests too, whatever they did.
+    summarise = names.index("Summarise the calibration fit")
+    assert names.index("Run tests") < summarise
+    assert steps[summarise]["if"] == "${{ !cancelled() }}"
+    assert steps[names.index("Save the fit summary")]["if"] == "${{ !cancelled() }}"
+    # On the base, the distribution measurement runs after the fit summary.
+    base = yaml.safe_load(workflow.read_text())["jobs"]["base-build"]["steps"]
+    base_names = [step.get("name") for step in base]
+    assert base_names.index("Save the base fit summary") < base_names.index(
+        "Measure the base's income distribution and single ages"
+    )

policyengine_uk_data/tests/test_fit_report.py

@@ -248,7 +248,10 @@ def test_the_commands_summarise_and_exit_one_on_a_failure(tmp_path):
 
 
 def test_the_pull_request_workflow_uploads_only_summaries():
-    """The datasets are licensed microdata, and raw calibration logs stay on the runner."""
+    """The datasets are licensed microdata, and raw calibration logs stay on the runner.
+
+    Only aggregate reports leave a job: the fit summaries and the distribution reports.
+    """
     from pathlib import Path
 
     import yaml
@@ -261,7 +264,12 @@ def test_the_pull_request_workflow_uploads_only_summaries():
         for step in job.get("steps", [])
         if step.get("uses", "").startswith("actions/upload-artifact")
     )
-    assert uploads == ["base_fit_summary.json", "fit_summary.json"]
+    assert uploads == [
+        "base_distribution_report.json",
+        "base_fit_summary.json",
+        "distribution_report.json",
+        "fit_summary.json",
+    ]
     for job in jobs.values():
         for step in job.get("steps", []):
             assert not step.get("uses", "").startswith("actions/cache"), step