From f3db8591e2a75af75e77fbdde81205b210b3799c Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Tue, 8 Sep 2026 12:05:49 +0000
Subject: [PATCH] Add the inequality package, so the repository's name is
 earned

Counted before starting: 27 Python files outside tests, 872 lines, of which
survey_estimation, impact_evaluation and one cleaning pipeline were 732.
Exactly one thing in the repository was about distribution, and the file
whose name promised it, public_health_access_index.py, was a five-line
row mean.

inequality/ is 24 functions across five modules.

Indices: Gini, the full generalised entropy family with Theil T and the
mean log deviation named, Atkinson at any aversion, Palma, the 80/20
ratio, quantile shares, the Lorenz curve, weighted quantiles.

Concentration: the concentration index, which is the measure that answers
the question an average cannot, whether an outcome falls on the poor and
how steeply. With it the Erreygers and Wagstaff corrections, because a
raw index is not comparable between two places whose prevalence differs
and comparing them anyway is a standard error in published work. Plus the
concentration curve, a by-group table, and Wagstaff's achievement index.

Decomposition: the Theil within/between split, exact for GE(0) and GE(1)
and refusing any other alpha rather than returning parts that do not sum.
Blinder-Oaxaca in twofold and threefold form with four reference choices.

Incidence: benefit incidence analysis, single-service and split by level
of service, which is where the finding usually is. Primary care
progressive, tertiary regressive, the two cancelling in the aggregate.

Three things are centralised rather than reimplemented per function.
Weighted fractional ranks live in one file and are shared by the Gini,
the concentration index and both curves, so those agree by construction.
Ties take the block midpoint, which is what makes a wealth quintile a
legitimate ranking variable. And every function refuses rather than
returning a plausible wrong number when an assumption breaks: a Gini on
negative values, GE(0) on zeros, Wagstaff where everyone has the outcome,
a level with no budget, an Oaxaca on three groups.

Tests: 25 to 101. Every new one is an identity or a hand-computed figure,
never previous output. The Gini against brute-force mean absolute
difference on fifteen random datasets. Curve areas against their indices
by trapezoid. within + between == total for both alphas, with weights.
Oaxaca's components summing to the gap under all four references. GE(2)
equalling half the squared coefficient of variation. A weight of two
equalling the row appearing twice. The mean fractional rank being exactly
0.5, which three derivations here depend on.

Two of my own test expectations were wrong on first run and the code was
right, which is the point of pinning to closed forms: an Oaxaca test drew
both groups from the same distribution and then asserted the realised
endowment gap was zero, and a composite-index test hardcoded arithmetic I
had done in my head. Both now compute the expectation from the definition.

Also: social_sector/composite_index.py replaces the row mean with an
index that makes its three decisions arguments rather than assumptions
(direction, normalisation, weights), offers goalpost normalisation so a
series is comparable across years, geometric aggregation where indicators
should not substitute for each other, a minimum non-missing count so a
district reporting two of six indicators is not scored as though it
reported all six, and rank_sensitivity to show how far a ranking moves
under an equally defensible choice. compute_access_index keeps working.

Deleted workflows/ci.yml: dead configuration, pinned to Python 3.10 and
running R tests in a repository with no R, duplicating the real workflow
at .github/workflows/python-tests.yml. The landing page described it as
"start-to-finish routes through an analysis", which it never was.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01DMYhECx2Dzbqad6KJxWZ7t
---
 CLAUDE.md                                   |  42 +-
 README.md                                   |  39 +-
 index.html                                  |  81 +++-
 inequality/README.md                        | 152 +++++++
 inequality/__init__.py                      |  56 +++
 inequality/concentration.py                 | 262 +++++++++++
 inequality/decomposition.py                 | 235 ++++++++++
 inequality/incidence.py                     | 197 +++++++++
 inequality/indices.py                       | 293 ++++++++++++
 inequality/ranks.py                         | 128 ++++++
 social_sector/composite_index.py            | 253 +++++++++++
 social_sector/public_health_access_index.py |  73 ++-
 tests/test_composite_index.py               | 211 +++++++++
 tests/test_inequality.py                    | 466 ++++++++++++++++++++
 workflows/ci.yml                            |  32 --
 15 files changed, 2467 insertions(+), 53 deletions(-)
 create mode 100644 inequality/README.md
 create mode 100644 inequality/__init__.py
 create mode 100644 inequality/concentration.py
 create mode 100644 inequality/decomposition.py
 create mode 100644 inequality/incidence.py
 create mode 100644 inequality/indices.py
 create mode 100644 inequality/ranks.py
 create mode 100644 social_sector/composite_index.py
 create mode 100644 tests/test_composite_index.py
 create mode 100644 tests/test_inequality.py
 delete mode 100644 workflows/ci.yml

diff --git a/CLAUDE.md b/CLAUDE.md
index b693a18..ee45990 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -7,10 +7,18 @@ Python scripts and notebooks for development sector data workflows. Part of the
 ## Layout
 
 `cleaning/`, `eda/`, `modelling/`, `validation/`, `io_helpers/`,
-`impact_evaluation/`, `social_sector/`, `visualisation/`, plus
+`impact_evaluation/`, `social_sector/`, `visualisation/`, plus the two that
+carry the repository: `inequality/` for distributional analysis and
 `survey_estimation/` for design-based estimates from complex surveys.
 Notebooks in `notebooks/`, small CSVs in `sample_data/`, tests in `tests/`.
 
+Sizes, so nobody has to guess again: 101 tests, and `inequality/` (about 950
+lines) plus `survey_estimation/` (357) plus `impact_evaluation/` (220) are most
+of the substance. The rest is technique taught against stand-in data, and some
+of it is a handful of lines. That is a deliberate split, not neglect, but do not
+describe a five-line snippet folder as a module in any user-facing copy: an
+audit on 2026-09-08 found the landing page doing exactly that.
+
 House style is numpydoc docstrings, one module per group of related functions,
 pytest with plain asserts. Match it.
 
@@ -43,6 +51,38 @@ pip install -r requirements.txt
 PYTHONPATH=$(pwd) pytest tests/
 ```
 
+## inequality
+
+The package that makes the repository's name true. Everything in it is a
+covariance between an outcome and a position in a distribution, and positions
+are where the errors hide, so three things are centralised rather than
+reimplemented per function.
+
+- **Weighted fractional ranks** live in `inequality/ranks.py` and are shared by
+  the Gini, the concentration index and both curves, so those three agree by
+  construction. Ties take the block midpoint, which is what the World Bank's DHS
+  equity work does and what makes a wealth quintile a legitimate rank variable.
+- **The mean fractional rank is exactly 0.5** for any weights and any pattern of
+  ties. Three derivations depend on it and a test asserts it.
+- **Sign convention: negative means concentrated among the poor.** Reversing
+  `outcome` and `rank_by` returns a number in the right range with the opposite
+  meaning and nothing errors, which is why the signature is
+  `concentration_index(outcome, rank_by=...)`.
+
+A raw concentration index is not comparable across different prevalences, so
+`erreygers_index` and `wagstaff_index` exist; they answer different normative
+questions and can disagree about the direction of change over time. Pick one per
+table and say which.
+
+**Do not add an analytic standard error to these.** The convenient regression
+form people quote ignores the survey design entirely. Bootstrap over PSUs within
+strata; `inequality/README.md` carries the recipe.
+
+Tests are pinned to closed-form identities, never to previous output: the Gini
+against brute-force mean absolute difference on random data, curve areas against
+their indices by trapezoid, `within + between == total` for both GE(0) and
+GE(1), Oaxaca's components summing to the gap under all four reference choices.
+
 ## survey_estimation
 
 The one part worth reading before touching. It implements Taylor linearisation
diff --git a/README.md b/README.md
index 2e34e52..798b0ee 100644
--- a/README.md
+++ b/README.md
@@ -6,7 +6,7 @@
 [![License: MIT](https://img.shields.io/badge/License-MIT-green.svg)](LICENSE)
 [![Status: Stable](https://img.shields.io/badge/Status-Stable-0969da?style=flat-square)](https://github.com/Varnasr/OpenStacks-for-Change/blob/main/MAINTENANCE.md)
 
-> Plug-and-play templates for health, gender, education, and climate equity data — built for reproducibility.
+> Distributional analysis and survey estimation for health, gender, education and climate equity data.
 
 > **Status: Stable.** This repository works and is correct, but it is not under active
 > development. Bug reports are welcome and issues stay open; new features are unlikely,
@@ -21,6 +21,40 @@ EquityStack is a collection of **Python scripts, Jupyter notebooks, and sample d
 
 This is the **data pipeline layer** of [OpenStacks for Change](https://openstacks.dev) — an open ecosystem of tools for public interest research and evaluation.
 
+## The two modules to read first
+
+Most of this repository is technique on stand-in data. Two parts are not, and
+they are what the name is about.
+
+**`inequality/`** answers *who has it*. The Gini says how unequally consumption
+is spread; the concentration index says whether stunting falls on the poor, by
+how much, and comparably across states whose prevalence differs. It also
+decomposes: `theil_decomposition` splits national inequality exactly into a
+within-state and a between-state part, and `oaxaca_blinder` splits a group gap
+into endowments and returns. `benefit_incidence` asks who actually receives a
+public budget, which for tertiary health and higher education is usually not the
+people it was voted for.
+
+**`survey_estimation/`** answers *how sure are we*. Weighting is the half
+everyone remembers; a national household survey is also clustered, and an
+interval that ignores that is too narrow, often by half, erring in the direction
+that flatters the result.
+
+Together they are the two halves of an equity finding: the gradient, and whether
+it is real.
+
+```python
+from inequality import concentration_index, theil_decomposition
+from survey_estimation import svy_prop_by
+
+svy_prop_by(df, "stunted", by="wealth_quintile")          # is the gap real?
+concentration_index(df.stunted, rank_by=df.wealth_index)  # how steep is it?
+theil_decomposition(df.consumption, groups=df.state)      # where does it sit?
+```
+
+`inequality/README.md` has the full function-by-function guide, the four things
+that are easy to get wrong, and how to bootstrap a standard error over PSUs.
+
 ## What's Inside
 
 ### Core Modules
@@ -33,7 +67,8 @@ This is the **data pipeline layer** of [OpenStacks for Change](https://openstack
 | `io_helpers/` | Chunked CSV reading, Stata/SPSS import, formatted Excel export | Ready |
 | `modelling/` | Multicollinearity checks (VIF) | Ready |
 | `visualisation/` | Annotated bar charts, district-level choropleth maps | Ready |
-| `social_sector/` | Public health access index | Ready |
+| `inequality/` | **Distributional analysis**: Gini, Theil, Atkinson, Palma, Lorenz and concentration curves, the concentration index with Erreygers and Wagstaff corrections, Theil within/between decomposition, Blinder-Oaxaca, benefit incidence | Ready |
+| `social_sector/` | Composite indices with explicit direction, normalisation and weighting, plus a rank-sensitivity check | Ready |
 | `survey_estimation/` | Design-based proportions and means for stratified, clustered surveys, with a worked NFHS-5 example | Ready |
 
 ### Notebooks
diff --git a/index.html b/index.html
index 417e30e..3834b16 100644
--- a/index.html
+++ b/index.html
@@ -17,8 +17,9 @@
   <div class="in">
     <a class="brand" href="./"><span class="dot"></span> EquityStack</a>
     <nav>
+      <a class="opt" href="#equity">Equity</a>
+      <a class="opt" href="#survey">Survey</a>
       <a class="opt" href="#modules">Modules</a>
-      <a class="opt" href="#survey">Survey estimation</a>
       <a href="https://github.com/Varnasr/EquityStack">GitHub</a>
     </nav>
   </div>
@@ -26,20 +27,72 @@
 
 <section class="hero">
   <h1 class="word">Equity<span class="s2">Stack</span></h1>
-  <p class="lede">The Python you reach for on a development research project: getting messy data into
-    shape, asking whether a programme worked, and getting the confidence intervals right
-    when the survey is clustered. Part of <a href="https://openstacks.dev">OpenStacks</a>.</p>
+  <p class="lede">The Python for a development research project, pointed at the distributional
+    question: not just what the average is, but who is carrying the bad outcome, how steep the
+    gradient runs, and whether the gap survives a clustered sample.
+    Part of <a href="https://openstacks.dev">OpenStacks</a>.</p>
   <div class="counts">
-    <div class="hi"><b>9</b><span>Modules</span></div>
-    <div><b>34</b><span>Python files</span></div>
-    <div><b>2</b><span>Notebooks</span></div>
+    <div class="hi"><b>24</b><span>Distributional measures</span></div>
+    <div><b>10</b><span>Modules</span></div>
     <div><b>4</b><span>Impact evaluation methods</span></div>
+    <div><b>101</b><span>Tests</span></div>
   </div>
 </section>
 
 <main id="main">
 
-<section class="band white" id="survey">
+<section class="band white" id="equity">
+  <div class="band-h">
+    <h2>Distributional analysis</h2>
+    <p class="sub">Who is carrying it, and how steeply</p>
+  </div>
+  <div class="tiles">
+    <a class="tile" href="https://github.com/Varnasr/EquityStack/tree/main/inequality">
+      <p class="meta">inequality/</p>
+      <h3>The concentration index</h3>
+      <p>An average tells you the level. It cannot tell you whether stunting falls on the poor,
+         which is the question a programme budget answers to. The concentration index does, on a
+         scale where negative means concentrated among the poor, and it works off a wealth
+         quintile as readily as a continuous consumption measure.</p>
+      <p class="go">Read the module &rarr;</p>
+    </a>
+    <a class="tile" href="https://github.com/Varnasr/EquityStack/blob/main/inequality/README.md">
+      <p class="meta">Comparability</p>
+      <h3>Why the raw index misleads</h3>
+      <p>A concentration index is not comparable between two places whose prevalence differs,
+         because its possible range shrinks as the mean moves away from half. Comparing states or
+         years without correcting for that is a standard error in published work. Both fixes are
+         here, and they can disagree about which way things moved.</p>
+      <p class="go">Read the guide &rarr;</p>
+    </a>
+    <a class="tile" href="https://github.com/Varnasr/EquityStack/blob/main/inequality/decomposition.py">
+      <p class="meta">Decomposition</p>
+      <h3>Where the inequality sits</h3>
+      <p>A national Gini of 0.35 is equally consistent with identical states that differ from each
+         other and with unequal states that look alike on average. Those want opposite responses.
+         The Theil split answers it exactly, and Blinder-Oaxaca separates how much of a group gap
+         is land and schooling from how much is the return on them.</p>
+      <p class="go">Read the module &rarr;</p>
+    </a>
+    <a class="tile" href="https://github.com/Varnasr/EquityStack/blob/main/inequality/incidence.py">
+      <p class="meta">Benefit incidence</p>
+      <h3>Who receives the budget</h3>
+      <p>Public tertiary hospitals and universities are free at the point of use and mostly used by
+         the better-off, so the spending is regressive whatever the intention. Splitting a budget
+         by level of service usually shows primary care progressive and tertiary regressive, the
+         two cancelling in any aggregate figure.</p>
+      <p class="go">Read the module &rarr;</p>
+    </a>
+  </div>
+  <p class="note">Twenty-four measures: Gini, Theil, Atkinson, Palma, the 80/20 ratio, Lorenz and
+    concentration curves, quantile shares, the Erreygers and Wagstaff corrections, within and
+    between decomposition, Blinder-Oaxaca, benefit incidence. All survey-weighted. Verified
+    against identities rather than against previous output, because 0.31 and 0.34 are both
+    plausible Ginis and a sign error gives a number in the right range that says the opposite of
+    the truth.</p>
+</section>
+
+<section class="band navy" id="survey">
   <div class="band-h">
     <h2>Survey estimation</h2>
     <p class="sub">The part worth reading first</p>
@@ -70,12 +123,15 @@ <h3>NFHS-5 stunting, reproduced</h3>
     bound, and the interval widens correctly when you have few clusters.</p>
 </section>
 
-<section class="band navy" id="modules">
+<section class="band white" id="modules">
   <div class="band-h">
     <h2>Modules</h2>
-    <p class="sub">Nine, plus notebooks and sample data</p>
+    <p class="sub">Ten, plus notebooks and sample data</p>
   </div>
   <div class="rows">
+    <a class="row" href="https://github.com/Varnasr/EquityStack/tree/main/inequality">
+      <span class="t"><b>Inequality</b><span>Python</span></span>
+      <span class="d">Gini, Theil, Atkinson, Palma, concentration indices, decompositions and benefit incidence, all survey-weighted.</span></a>
     <a class="row" href="https://github.com/Varnasr/EquityStack/tree/main/cleaning">
       <span class="t"><b>Cleaning</b><span>Python</span></span>
       <span class="d">Get a messy export into shape, and keep a record of what you changed.</span></a>
@@ -99,7 +155,7 @@ <h2>Modules</h2>
       <span class="d">Open a Stata or SPSS file, read a CSV too big for memory, hand a colleague a usable spreadsheet.</span></a>
     <a class="row" href="https://github.com/Varnasr/EquityStack/tree/main/social_sector">
       <span class="t"><b>Social sector</b><span>Python</span></span>
-      <span class="d">Build a health access index and compare districts on it.</span></a>
+      <span class="d">Build an access or vulnerability index, and check whether a district's rank survives a different but equally defensible normalisation.</span></a>
     <a class="row" href="https://github.com/Varnasr/EquityStack/tree/main/visualisation">
       <span class="t"><b>Visualisation</b><span>Python</span></span>
       <span class="d">District maps and annotated charts that can go straight into a report.</span></a>
@@ -109,9 +165,6 @@ <h2>Modules</h2>
     <a class="row" href="https://github.com/Varnasr/EquityStack/tree/main/sample_data">
       <span class="t"><b>Sample data</b><span>CSV</span></span>
       <span class="d">Small datasets to try any of the above on.</span></a>
-    <a class="row" href="https://github.com/Varnasr/EquityStack/tree/main/workflows">
-      <span class="t"><b>Workflows</b><span>Guides</span></span>
-      <span class="d">Start-to-finish routes through an analysis, so you know what comes next.</span></a>
   </div>
   <p class="note">These open on GitHub, where the code lives.</p>
 </section>
diff --git a/inequality/README.md b/inequality/README.md
new file mode 100644
index 0000000..a8af0ae
--- /dev/null
+++ b/inequality/README.md
@@ -0,0 +1,152 @@
+# inequality
+
+Distributional analysis for survey microdata: how unequally something is spread,
+and across whom.
+
+The rest of this repository answers "what is the level, and how sure are we".
+This package answers the question that follows: **who has it**. Those are
+different questions and the second one is the harder to get right, because every
+measure of it is a covariance between an outcome and a position in a
+distribution, and positions are where the errors hide.
+
+```python
+from inequality import gini, concentration_index, theil_decomposition
+
+gini(df.consumption, weights=df.hh_weight)
+concentration_index(df.stunted, rank_by=df.wealth_index, weights=df.child_weight)
+theil_decomposition(df.consumption, groups=df.state, weights=df.hh_weight)
+```
+
+## What to reach for
+
+**One variable, how unequal is it?**
+
+| Function | Use when |
+|---|---|
+| `gini` | The headline. Comparable with almost every published figure. |
+| `theil_t`, `theil_l` | You need the within/between split. Only GE(0) and GE(1) decompose exactly. |
+| `atkinson` | You want the inequality aversion stated rather than implied. |
+| `generalised_entropy` | You want to choose where the index is sensitive: GE(0) bottom, GE(1) neutral, GE(2) top. |
+| `palma_ratio`, `ratio_80_20` | The audience is not technical. Both are ratios of shares. |
+| `quantile_shares`, `lorenz_curve`, `share_of_top` | You want the distribution itself, not a summary of it. |
+
+**Two variables, unequal across whom?**
+
+| Function | Use when |
+|---|---|
+| `concentration_index` | The core: is this outcome concentrated among the poor? Negative means yes. |
+| `erreygers_index` | Comparing the same indicator across states or years whose mean differs. |
+| `wagstaff_index` | The same problem, relative rather than absolute. Pick one and say which. |
+| `concentration_curve` | Plotting it, or checking dominance between two distributions. |
+| `concentration_index_by` | A table: one gradient per state, per round, per sex. |
+| `achievement_index` | One number combining level and distribution. Report it beside the mean, never instead. |
+
+**What explains it?**
+
+| Function | Use when |
+|---|---|
+| `theil_decomposition` | Is national inequality *between* states or *within* them? |
+| `oaxaca_blinder` | How much of a group gap is endowments, how much is returns? |
+
+**Who gets the money?**
+
+| Function | Use when |
+|---|---|
+| `benefit_incidence` | Allocating a budget across quintiles by who uses the service. |
+| `benefit_incidence_by_level` | The same, split by primary / secondary / tertiary, which is usually where the finding is. |
+
+## Four things that are easy to get wrong
+
+**The concentration index ranks by living standards, not by the outcome.**
+Reversing them returns a number in the right range with the wrong meaning, and
+nothing errors. The signature makes it awkward to do by accident: outcome first,
+then `rank_by=`.
+
+**Sign convention.** Negative means concentrated among the *poor*. For stunting
+that is the expected direction and it is bad news; for institutional delivery it
+is good news. The index has no opinion about which.
+
+**A raw concentration index is not comparable across different prevalences.** Its
+theoretical range shrinks as a binary outcome's mean moves away from 0.5, so an
+indicator at 8 per cent and the same indicator at 60 per cent cannot be put in
+one table without a correction. Erreygers keeps absolute gradients comparable,
+Wagstaff keeps relative ones, and they can disagree about the direction of
+change over time. Choose before you look at the answer.
+
+**Ties.** A wealth quintile is a legitimate ranking variable and this package
+handles it the way the World Bank's own DHS equity work does: everyone sharing a
+quintile shares the midpoint rank of the block that quintile occupies. The index
+is attenuated relative to a continuous wealth measure, which is a property of the
+coarser data and not a defect. `test_inequality.py` asserts the attenuation runs
+in the expected direction.
+
+## Standard errors
+
+None of these functions return one. The analytic variance of a concentration
+index under a stratified clustered design is not something to write from memory,
+and the convenient regression form people quote gives a standard error that
+ignores the design entirely.
+
+Bootstrap over PSUs within strata instead, resampling clusters rather than
+households, and use this repository's `survey_estimation` package for the design
+object. Roughly:
+
+```python
+import numpy as np
+from inequality import concentration_index
+
+def boot_ci(df, reps=500, seed=0):
+    rng = np.random.default_rng(seed)
+    out = []
+    for _ in range(reps):
+        parts = []
+        for _, stratum in df.groupby("strata"):
+            psus = stratum["psu"].unique()
+            drawn = rng.choice(psus, size=len(psus), replace=True)
+            parts.append(pd.concat([stratum[stratum.psu == p] for p in drawn]))
+        rep = pd.concat(parts)
+        out.append(concentration_index(rep.stunted, rep.wealth_index, rep.weight))
+    return np.percentile(out, [2.5, 97.5])
+```
+
+## Verification
+
+`tests/test_inequality.py`, 54 checks. Every one is either an identity the
+measure must satisfy or a figure worked out by hand from the definition, never
+a number this code produced on a previous run. A distributional index cannot be
+checked by eye: 0.31 and 0.34 are both plausible consumption Ginis, and a sign
+error gives you a number in exactly the right range that says the opposite of
+the truth.
+
+The identities that carry the most weight:
+
+- The Gini equals the mean absolute difference over twice the mean, on fifteen
+  random datasets, computed independently by brute force over all pairs.
+- Twice the area between the Lorenz curve and the diagonal equals the Gini, and
+  the same for the concentration curve and the concentration index.
+- The mean fractional rank is exactly 0.5 for any weights and any pattern of
+  ties. Three derivations here depend on it.
+- `within + between == total` in the Theil decomposition, for both GE(0) and
+  GE(1), with weights.
+- Oaxaca's components sum to the gap under all four reference choices.
+- A weight of 2 gives the same answer as the row appearing twice.
+- GE(2) equals half the squared coefficient of variation.
+
+## Sources
+
+Cowell, *Measuring Inequality*, 3rd edn, Oxford University Press 2011, chapters
+2-3, for the entropy family and the axioms.
+
+O'Donnell, van Doorslaer, Wagstaff and Lindelow, *Analyzing Health Equity Using
+Household Survey Data*, World Bank 2008, chapters 8 and 15, for the
+concentration index and its corrections.
+
+Erreygers, "Correcting the concentration index", *Journal of Health Economics*
+28(2), 2009, 504-515. Wagstaff, "The bounds of the concentration index when the
+variable of interest is binary", *Health Economics* 14(4), 2005, 429-432.
+
+Demery, *Benefit Incidence: A Practitioner's Guide*, World Bank 2000.
+
+Fortin, Lemieux and Firpo, "Decomposition methods in economics", *Handbook of
+Labor Economics* 4A, 2011, section 3, on what the unexplained component is and
+is not.
diff --git a/inequality/__init__.py b/inequality/__init__.py
new file mode 100644
index 0000000..e849666
--- /dev/null
+++ b/inequality/__init__.py
@@ -0,0 +1,56 @@
+"""Distributional analysis: how unequally something is spread, and across whom.
+
+    from inequality import gini, concentration_index, theil_decomposition
+
+Four groups of function, and the choice between them is the choice of question.
+
+**One variable, how unequal?**  `gini`, `theil_t`, `theil_l`, `atkinson`,
+`generalised_entropy`, `palma_ratio`, `ratio_80_20`, `quantile_shares`,
+`lorenz_curve`, `share_of_top`, `share_of_bottom`.
+
+**Two variables, unequal across whom?**  `concentration_index` and its
+prevalence-comparable corrections `erreygers_index` and `wagstaff_index`, plus
+`concentration_curve`, `concentration_index_by` for a table by state or round,
+and `achievement_index`.
+
+**Where does it sit, and what explains a gap?**  `theil_decomposition` splits
+inequality exactly into within-group and between-group parts.
+`oaxaca_blinder` splits a mean gap into endowments and returns.
+
+**Who gets the money?**  `benefit_incidence` and `benefit_incidence_by_level`.
+
+Everything takes survey weights, drops missing rows pairwise rather than
+imputing, and refuses rather than returning a plausible wrong number when the
+input violates an assumption the measure depends on. Weighted fractional ranks
+are shared through `ranks.fractional_rank`, so the Gini, the concentration index
+and both curves agree with each other by construction.
+
+For a standard error on any of these, resample: the analytic variance of a
+concentration index under a clustered design is not something to write from
+memory, and a bootstrap over PSUs within strata gets there honestly. The design
+machinery is in this repository's `survey_estimation` package.
+"""
+
+from .ranks import fractional_rank, weighted_mean
+from .indices import (
+    gini, lorenz_curve, share_of_top, share_of_bottom, quantile_shares,
+    palma_ratio, ratio_80_20, generalised_entropy, theil_t, theil_l,
+    atkinson, weighted_quantile,
+)
+from .concentration import (
+    concentration_index, erreygers_index, wagstaff_index, concentration_curve,
+    concentration_index_by, achievement_index,
+)
+from .decomposition import theil_decomposition, oaxaca_blinder
+from .incidence import benefit_incidence, benefit_incidence_by_level
+
+__all__ = [
+    "fractional_rank", "weighted_mean",
+    "gini", "lorenz_curve", "share_of_top", "share_of_bottom", "quantile_shares",
+    "palma_ratio", "ratio_80_20", "generalised_entropy", "theil_t", "theil_l",
+    "atkinson", "weighted_quantile",
+    "concentration_index", "erreygers_index", "wagstaff_index",
+    "concentration_curve", "concentration_index_by", "achievement_index",
+    "theil_decomposition", "oaxaca_blinder",
+    "benefit_incidence", "benefit_incidence_by_level",
+]
diff --git a/inequality/concentration.py b/inequality/concentration.py
new file mode 100644
index 0000000..54a7a11
--- /dev/null
+++ b/inequality/concentration.py
@@ -0,0 +1,262 @@
+"""Concentration indices: inequality in one variable, ranked by another.
+
+The Gini asks how unequally consumption is spread. It cannot answer the question
+an equity analyst actually has, which is whether *stunting*, or *institutional
+delivery*, or *out-of-pocket health spending* falls more heavily on the poor,
+and by how much. That needs two variables: the outcome, and a living-standards
+rank to order people by.
+
+The concentration index is the answer, and it is one line of arithmetic:
+
+    CI = 2 · cov(h, r) / h̄
+
+where ``h`` is the health or service variable and ``r`` is the weighted
+fractional rank of the living-standards variable. It runs from -1 to +1.
+**Negative means concentrated among the poor.** Stunting gives a negative index;
+private hospital use gives a positive one. Zero means the outcome is spread
+evenly across the distribution, which for a bad outcome is not good news, it is
+an absence of a gradient.
+
+Three things this module handles that a hand-rolled version usually does not.
+
+**Bounded variables need a correction.** The index's theoretical range shrinks
+as the mean of a binary variable moves away from 0.5, so a raw CI on an outcome
+with 8 per cent prevalence cannot be compared with one at 60 per cent, and
+comparing them across states or years is exactly what people do. Both standard
+fixes are here: Erreygers (2009) and Wagstaff (2005). Say which you used.
+
+**Ranks come from the living-standards variable, not the outcome.** Reversing
+them is the commonest error and it produces a number in the right range with the
+wrong meaning. The argument order here makes it hard: outcome first, then
+``rank_by=``.
+
+**A wealth quintile is a legitimate rank variable.** Ties get the block
+midpoint, so 20 per cent of the population sitting on the integer 2 all rank at
+0.3, which is what the World Bank's own DHS equity reports do. It is coarser
+than a continuous consumption measure and the index is attenuated accordingly;
+that is a property of the data, not a bug in the code.
+
+References: O'Donnell, van Doorslaer, Wagstaff and Lindelow, *Analyzing Health
+Equity Using Household Survey Data* (World Bank 2008), chapters 8 and 15;
+Erreygers, "Correcting the concentration index", *Journal of Health Economics*
+28(2), 2009, 504-515; Wagstaff, "The bounds of the concentration index when the
+variable of interest is binary", *Health Economics* 14(4), 2005, 429-432.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+
+from .ranks import fractional_rank
+
+__all__ = [
+    "concentration_index", "erreygers_index", "wagstaff_index",
+    "concentration_curve", "concentration_index_by", "achievement_index",
+]
+
+
+def _aligned(outcome, rank_by, weights):
+    """Drop rows missing any of outcome, rank variable or weight, pairwise."""
+    h = np.asarray(pd.to_numeric(pd.Series(outcome), errors="coerce"), dtype=float)
+    x = np.asarray(pd.to_numeric(pd.Series(rank_by), errors="coerce"), dtype=float)
+    if h.shape != x.shape:
+        raise ValueError(
+            f"outcome has {h.shape[0]} rows but rank_by has {x.shape[0]}")
+    if weights is None:
+        w = np.ones_like(h)
+    else:
+        w = np.asarray(pd.to_numeric(pd.Series(weights), errors="coerce"), dtype=float)
+        if w.shape != h.shape:
+            raise ValueError(
+                f"outcome has {h.shape[0]} rows but weights has {w.shape[0]}")
+
+    keep = np.isfinite(h) & np.isfinite(x) & np.isfinite(w)
+    h, x, w = h[keep], x[keep], w[keep]
+    if h.size == 0:
+        raise ValueError("no rows left after dropping missing values")
+    if np.any(w < 0):
+        raise ValueError("negative weights are not a thing; check the weight column")
+    if w.sum() <= 0:
+        raise ValueError("weights sum to zero, so nothing can be estimated")
+    return h, x, w
+
+
+def concentration_index(outcome, rank_by, weights=None) -> float:
+    """Concentration index of ``outcome`` over the ranking of ``rank_by``.
+
+    Parameters
+    ----------
+    outcome
+        The health, service or spending variable. Binary is fine and common.
+    rank_by
+        The living-standards variable people are ordered by: a consumption
+        aggregate, a wealth index, or a quintile. Higher must mean better off.
+    weights
+        Survey weights. Relative or absolute, it makes no difference.
+
+    Returns a number in [-1, 1]. Negative: concentrated among the poor.
+
+    >>> # an outcome that falls entirely on the poorest gives a strong negative
+    >>> round(concentration_index([1, 0, 0, 0], rank_by=[1, 2, 3, 4]), 4)
+    -0.75
+    """
+    h, x, w = _aligned(outcome, rank_by, weights)
+    mean = np.sum(w * h) / np.sum(w)
+    if mean == 0:
+        raise ValueError(
+            "concentration_index: the weighted mean of the outcome is zero, so the "
+            "index is undefined. Nobody in the sample has the outcome.")
+    r = fractional_rank(x, w)
+    return float(2.0 * np.sum(w * h * r) / (np.sum(w) * mean) - 1.0)
+
+
+def erreygers_index(outcome, rank_by, weights=None, *,
+                    bounds: tuple[float, float] = (0.0, 1.0)) -> float:
+    """Erreygers-corrected concentration index, comparable across prevalences.
+
+    ``E = 4 · μ / (b − a) · CI``, where ``(a, b)`` are the variable's true
+    bounds. For a binary outcome the default (0, 1) is right and ``E = 4 μ CI``.
+
+    Use this when comparing the same indicator across states, years or surveys
+    whose mean differs. Two states can have identical raw concentration indices
+    and very different absolute gradients, and it is the absolute gradient that
+    a programme budget responds to.
+    """
+    a, b = bounds
+    if not b > a:
+        raise ValueError("erreygers_index: bounds must satisfy b > a")
+    h, x, w = _aligned(outcome, rank_by, weights)
+    lo, hi = np.nanmin(h), np.nanmax(h)
+    if lo < a - 1e-12 or hi > b + 1e-12:
+        raise ValueError(
+            f"erreygers_index: the outcome runs [{lo:g}, {hi:g}], outside the stated "
+            f"bounds [{a:g}, {b:g}]. The correction is meaningless on the wrong bounds.")
+    mean = np.sum(w * h) / np.sum(w)
+    ci = concentration_index(h, x, w)
+    return float(4.0 * mean / (b - a) * ci)
+
+
+def wagstaff_index(outcome, rank_by, weights=None, *, upper: float = 1.0) -> float:
+    """Wagstaff-normalised concentration index for a bounded variable.
+
+    ``W = CI / (1 − μ/b)``. Defined for a variable on [0, b]; undefined when the
+    mean reaches the upper bound, because then no inequality is possible and the
+    denominator is zero.
+
+    Erreygers and Wagstaff answer different normative questions and disagree,
+    sometimes about direction of change over time. Erreygers keeps absolute
+    differences comparable; Wagstaff keeps relative ones. Pick one, state it,
+    and do not switch between them inside a single table.
+    """
+    if upper <= 0:
+        raise ValueError("wagstaff_index: upper must be positive")
+    h, x, w = _aligned(outcome, rank_by, weights)
+    if np.nanmin(h) < -1e-12 or np.nanmax(h) > upper + 1e-12:
+        raise ValueError(
+            f"wagstaff_index: the outcome falls outside [0, {upper:g}]")
+    mean = np.sum(w * h) / np.sum(w)
+    denom = 1.0 - mean / upper
+    if abs(denom) < 1e-12:
+        raise ValueError(
+            "wagstaff_index: the mean equals the upper bound, so the normalisation "
+            "divides by zero. Everyone has the outcome; there is no gradient to scale.")
+    return float(concentration_index(h, x, w) / denom)
+
+
+def concentration_curve(outcome, rank_by, weights=None) -> pd.DataFrame:
+    """Points on the concentration curve, poorest to richest.
+
+    Cumulative share of the outcome against cumulative share of the population
+    ranked by living standards. Above the 45 degree line means concentrated
+    among the poor.
+
+    Sign convention, since it is the easy thing to get backwards: with ``A`` the
+    area under the curve, ``CI = 2 · (0.5 − A)``. A curve above the line has
+    ``A > 0.5`` and a negative index. `test_inequality.py` checks that identity
+    by trapezoid against `concentration_index`.
+
+    One row per distinct value of the rank variable, plus the origin, so a
+    wealth quintile gives six rows.
+    """
+    h, x, w = _aligned(outcome, rank_by, weights)
+    if np.any(h < 0):
+        raise ValueError("concentration_curve: a negative outcome makes the curve non-monotone")
+    order = np.argsort(x, kind="mergesort")
+    sh, sx, sw = h[order], x[order], w[order]
+
+    starts = np.concatenate(([True], sx[1:] != sx[:-1]))
+    group = np.cumsum(starts) - 1
+    gw = np.zeros(group[-1] + 1)
+    gh = np.zeros(group[-1] + 1)
+    np.add.at(gw, group, sw)
+    np.add.at(gh, group, sw * sh)
+
+    total_h = gh.sum()
+    pop = np.concatenate(([0.0], np.cumsum(gw) / gw.sum()))
+    val = np.concatenate(([0.0], np.cumsum(gh) / total_h)) if total_h > 0 else pop.copy()
+    return pd.DataFrame({"population_share": pop, "outcome_share": val,
+                         "rank_value": np.concatenate(([np.nan], sx[starts]))})
+
+
+def concentration_index_by(df: pd.DataFrame, outcome: str, rank_by: str, *,
+                           by: str, weights: str | None = None,
+                           bounds: tuple[float, float] = (0.0, 1.0)) -> pd.DataFrame:
+    """One concentration index per group, with the mean and Erreygers alongside.
+
+    The usual table: an equity gradient per state, per survey round, per sex.
+    Groups whose index cannot be computed (nobody has the outcome, one row, no
+    variation in the rank variable) come back with NaN and a stated ``note``
+    rather than being dropped, so a state missing from the output is visible.
+
+    **Each group is ranked within itself.** A state's index answers "is this
+    unequal *within* this state", not "are people in this state poor relative to
+    India". Those are different questions and pooling the rank would answer the
+    second while looking like the first.
+    """
+    for col in (outcome, rank_by, by) + ((weights,) if weights else ()):
+        if col not in df.columns:
+            raise KeyError(f"concentration_index_by: no column {col!r} in the frame")
+
+    rows = []
+    for key, g in df.groupby(by, dropna=False, observed=True):
+        w = g[weights] if weights else None
+        note = ""
+        ci = erre = mean = np.nan
+        try:
+            h, x, ww = _aligned(g[outcome], g[rank_by], w)
+            mean = float(np.sum(ww * h) / np.sum(ww))
+            if len(np.unique(x)) < 2:
+                note = "the rank variable takes one value in this group"
+            elif mean == 0:
+                note = "nobody in this group has the outcome"
+            else:
+                ci = concentration_index(h, x, ww)
+                erre = 4.0 * mean / (bounds[1] - bounds[0]) * ci
+        except ValueError as exc:
+            note = str(exc)
+        rows.append({by: key, "n": len(g), "mean": mean,
+                     "concentration_index": ci, "erreygers": erre, "note": note})
+    return pd.DataFrame(rows)
+
+
+def achievement_index(outcome, rank_by, weights=None, *, v: float = 2.0) -> float:
+    """Wagstaff's achievement index: the mean, discounted for how unequal it is.
+
+    ``A = μ · (1 − CI)`` at the usual inequality aversion of ``v = 2``. It
+    collapses level and distribution into one number, so a programme that raises
+    coverage only among the better-off scores worse than its headline coverage
+    suggests.
+
+    Reported alongside the mean, never instead of it. A single number that mixes
+    level and distribution cannot be read backwards into either.
+    """
+    if v < 1:
+        raise ValueError("achievement_index: v must be at least 1")
+    h, x, w = _aligned(outcome, rank_by, weights)
+    mean = float(np.sum(w * h) / np.sum(w))
+    if mean == 0:
+        return 0.0
+    r = fractional_rank(x, w)
+    weight = v * (1.0 - r) ** (v - 1.0)
+    return float(np.sum(w * weight * h) / np.sum(w))
diff --git a/inequality/decomposition.py b/inequality/decomposition.py
new file mode 100644
index 0000000..e988c6f
--- /dev/null
+++ b/inequality/decomposition.py
@@ -0,0 +1,235 @@
+"""Decompositions: where inequality sits, and how much of a gap is explained.
+
+Two questions a summary index cannot answer.
+
+*Where does the inequality sit?* A national Gini of 0.35 is consistent with
+every state being identical internally and very different from each other, and
+with every state being internally unequal and identical on average. Those imply
+opposite policy responses. `theil_decomposition` splits the total exactly into a
+within-group and a between-group part, which is the property GE(0) and GE(1)
+have and the Gini does not.
+
+*How much of a gap is composition?* Scheduled Tribe households have lower mean
+consumption than others. Some of that gap is that they hold less land and less
+schooling; some of it is a different return on the same land and the same
+schooling. `oaxaca_blinder` separates the two. The second part is often labelled
+discrimination, and that label is a claim the arithmetic does not support: it is
+the residual, and it carries every determinant left out of the model.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+
+from .indices import generalised_entropy
+from .ranks import prepare
+
+__all__ = ["theil_decomposition", "oaxaca_blinder"]
+
+
+def theil_decomposition(values, groups, weights=None, *, alpha: float = 1.0) -> dict:
+    """Split GE(0) or GE(1) into within-group and between-group inequality.
+
+    Parameters
+    ----------
+    alpha
+        1 for Theil's T, which weights groups by their share of the total; 0 for
+        the mean log deviation, which weights them by population share. No other
+        value decomposes additively, and passing one raises rather than
+        returning a number whose parts do not sum.
+
+    Returns a dict with ``total``, ``within``, ``between``, ``between_share``
+    and a per-group frame. ``within + between == total`` to floating-point
+    tolerance, asserted in the tests, because a decomposition whose parts do not
+    add up is worse than no decomposition.
+
+    The between-group share is the interpretable number. Reading it: a between
+    share of 0.15 means that if you equalised every state's internal
+    distribution while leaving state means alone, 85 per cent of national
+    inequality would remain. Across Indian states, consumption inequality is
+    overwhelmingly within rather than between, which is the standard finding
+    and the reason state-level averages mislead.
+    """
+    if not (np.isclose(alpha, 0.0) or np.isclose(alpha, 1.0)):
+        raise ValueError(
+            f"theil_decomposition: alpha = {alpha:g} does not decompose additively. "
+            "Only GE(0) and GE(1) do. Use alpha=0 or alpha=1.")
+
+    v = pd.Series(pd.to_numeric(pd.Series(values).reset_index(drop=True), errors="coerce"))
+    g = pd.Series(groups).reset_index(drop=True)
+    w = (pd.Series(np.ones(len(v))) if weights is None
+         else pd.Series(pd.to_numeric(pd.Series(weights).reset_index(drop=True), errors="coerce")))
+    if not (len(v) == len(g) == len(w)):
+        raise ValueError("theil_decomposition: values, groups and weights differ in length")
+
+    frame = pd.DataFrame({"v": v, "g": g, "w": w}).dropna(subset=["v", "w"])
+    if frame.empty:
+        raise ValueError("theil_decomposition: no rows left after dropping missing values")
+    if (frame["v"] <= 0).any():
+        raise ValueError(
+            "theil_decomposition: the variable must be strictly positive. Both GE(0) "
+            "and GE(1) take a logarithm of it.")
+
+    total_w = frame["w"].sum()
+    grand_mean = float((frame["v"] * frame["w"]).sum() / total_w)
+    total = generalised_entropy(frame["v"], frame["w"], alpha=alpha)
+
+    rows, within = [], 0.0
+    for key, part in frame.groupby("g", dropna=False, observed=True):
+        pw = part["w"].sum()
+        pop_share = float(pw / total_w)
+        mean = float((part["v"] * part["w"]).sum() / pw)
+        value_share = float((part["v"] * part["w"]).sum() / (grand_mean * total_w))
+        inner = (generalised_entropy(part["v"], part["w"], alpha=alpha)
+                 if len(part) > 1 else 0.0)
+        share = value_share if np.isclose(alpha, 1.0) else pop_share
+        within += share * inner
+        rows.append({"group": key, "n": len(part), "population_share": pop_share,
+                     "value_share": value_share, "mean": mean,
+                     "within_group_index": inner,
+                     "contribution_to_within": share * inner})
+
+    by_group = pd.DataFrame(rows)
+    if np.isclose(alpha, 1.0):
+        between = float((by_group["value_share"]
+                         * np.log(by_group["mean"] / grand_mean)).sum())
+    else:
+        between = float((by_group["population_share"]
+                         * np.log(grand_mean / by_group["mean"])).sum())
+
+    return {
+        "index": "GE(1), Theil T" if np.isclose(alpha, 1.0) else "GE(0), mean log deviation",
+        "total": float(total),
+        "within": float(within),
+        "between": between,
+        "between_share": float(between / total) if total != 0 else float("nan"),
+        "grand_mean": grand_mean,
+        "by_group": by_group,
+    }
+
+
+def oaxaca_blinder(df: pd.DataFrame, outcome: str, group: str, predictors: list[str], *,
+                   weights: str | None = None, advantaged=None,
+                   reference: str = "pooled") -> dict:
+    """Blinder-Oaxaca decomposition of the mean gap between two groups.
+
+    Splits ``ȳ_A − ȳ_B`` into a part explained by differences in the predictors
+    (endowments) and a part explained by differences in the coefficients on
+    them.
+
+    Parameters
+    ----------
+    group
+        A column with exactly two distinct values.
+    advantaged
+        Which value is group A. Defaults to whichever has the higher weighted
+        mean outcome, so the reported gap is positive and reads naturally.
+    reference
+        Which coefficient vector counts as non-discriminatory. ``"pooled"``
+        (Neumark 1988) fits one model on both groups and is the usual default.
+        ``"advantaged"`` or ``"disadvantaged"`` use one group's own
+        coefficients, which is the classic twofold form. ``"threefold"`` returns
+        the endowments / coefficients / interaction split instead, with the
+        disadvantaged group's coefficients as the base.
+
+    Returns a dict with the gap, the components, and a per-predictor frame
+    showing each variable's contribution.
+
+    **On what the unexplained part is.** It is a residual. It contains genuine
+    differences in returns, and it also contains every predictor you left out,
+    every mismeasured one, and the consequences of selection into the sample. A
+    large unexplained share is a reason to look harder, not a measurement of
+    discrimination, and the index cannot tell the two apart. The identification
+    problem is Oaxaca's own; see Fortin, Lemieux and Firpo, "Decomposition
+    methods in economics", *Handbook of Labor Economics* 4A (2011), section 3.
+    """
+    import statsmodels.api as sm
+
+    for col in [outcome, group, *predictors] + ([weights] if weights else []):
+        if col not in df.columns:
+            raise KeyError(f"oaxaca_blinder: no column {col!r} in the frame")
+
+    cols = [outcome, group, *predictors] + ([weights] if weights else [])
+    data = df[cols].dropna()
+    if data.empty:
+        raise ValueError("oaxaca_blinder: no complete rows")
+
+    levels = pd.unique(data[group])
+    if len(levels) != 2:
+        raise ValueError(
+            f"oaxaca_blinder: {group!r} takes {len(levels)} values; this decomposition "
+            "compares exactly two groups. Subset the frame first.")
+
+    w_all = data[weights].to_numpy(float) if weights else np.ones(len(data))
+
+    def wmean(mask, col):
+        return float(np.sum(data.loc[mask, col].to_numpy(float) * w_all[mask])
+                     / np.sum(w_all[mask]))
+
+    if advantaged is None:
+        means = {lv: wmean((data[group] == lv).to_numpy(), outcome) for lv in levels}
+        advantaged = max(means, key=means.get)
+    if advantaged not in set(levels):
+        raise ValueError(f"oaxaca_blinder: advantaged={advantaged!r} is not one of {list(levels)}")
+    disadvantaged = [lv for lv in levels if lv != advantaged][0]
+
+    def fit(mask):
+        X = sm.add_constant(data.loc[mask, predictors].to_numpy(float), has_constant="add")
+        y = data.loc[mask, outcome].to_numpy(float)
+        return sm.WLS(y, X, weights=w_all[mask]).fit()
+
+    mask_a = (data[group] == advantaged).to_numpy()
+    mask_b = (data[group] == disadvantaged).to_numpy()
+    for name, mask in (("advantaged", mask_a), ("disadvantaged", mask_b)):
+        if mask.sum() <= len(predictors) + 1:
+            raise ValueError(
+                f"oaxaca_blinder: the {name} group has {mask.sum()} rows for "
+                f"{len(predictors) + 1} parameters. The fit is not identified.")
+
+    fit_a, fit_b = fit(mask_a), fit(mask_b)
+    beta_a, beta_b = fit_a.params, fit_b.params
+
+    names = ["(intercept)", *predictors]
+    xbar_a = np.array([1.0] + [wmean(mask_a, p) for p in predictors])
+    xbar_b = np.array([1.0] + [wmean(mask_b, p) for p in predictors])
+    gap = wmean(mask_a, outcome) - wmean(mask_b, outcome)
+
+    if reference == "threefold":
+        endow = (xbar_a - xbar_b) * beta_b
+        coeff = xbar_b * (beta_a - beta_b)
+        inter = (xbar_a - xbar_b) * (beta_a - beta_b)
+        parts = pd.DataFrame({"term": names, "endowments": endow,
+                              "coefficients": coeff, "interaction": inter})
+        out = {"form": "threefold (base: disadvantaged group)",
+               "endowments": float(endow.sum()), "coefficients": float(coeff.sum()),
+               "interaction": float(inter.sum())}
+    else:
+        if reference == "pooled":
+            X = sm.add_constant(data[predictors].to_numpy(float), has_constant="add")
+            beta_star = sm.WLS(data[outcome].to_numpy(float), X, weights=w_all).fit().params
+            label = "twofold, pooled reference (Neumark)"
+        elif reference == "advantaged":
+            beta_star, label = beta_a, "twofold, advantaged group's coefficients"
+        elif reference == "disadvantaged":
+            beta_star, label = beta_b, "twofold, disadvantaged group's coefficients"
+        else:
+            raise ValueError(
+                "oaxaca_blinder: reference must be 'pooled', 'advantaged', "
+                "'disadvantaged' or 'threefold'")
+        explained = (xbar_a - xbar_b) * beta_star
+        unexplained = xbar_a * (beta_a - beta_star) + xbar_b * (beta_star - beta_b)
+        parts = pd.DataFrame({"term": names, "explained": explained,
+                              "unexplained": unexplained})
+        out = {"form": label, "explained": float(explained.sum()),
+               "unexplained": float(unexplained.sum())}
+
+    out.update({
+        "advantaged": advantaged, "disadvantaged": disadvantaged,
+        "mean_advantaged": wmean(mask_a, outcome),
+        "mean_disadvantaged": wmean(mask_b, outcome),
+        "gap": gap,
+        "n_advantaged": int(mask_a.sum()), "n_disadvantaged": int(mask_b.sum()),
+        "by_term": parts,
+    })
+    return out
diff --git a/inequality/incidence.py b/inequality/incidence.py
new file mode 100644
index 0000000..0428584
--- /dev/null
+++ b/inequality/incidence.py
@@ -0,0 +1,197 @@
+"""Benefit incidence analysis: who actually receives public spending.
+
+A health or education budget is spent on facilities. Whether the money reaches
+the poor depends on who uses those facilities, and the answer is routinely the
+opposite of the intention. Public tertiary hospitals and universities in most of
+South Asia are used disproportionately by the better-off, so a rupee spent there
+is regressive even though the facility is free at the point of use. Primary
+health centres and government primary schools usually run the other way.
+
+The method (Demery, *Benefit Incidence: A Practitioner's Guide*, World Bank
+2000) is deliberately simple arithmetic on two inputs:
+
+1. a household survey giving utilisation by living-standards group, and
+2. a public expenditure figure, ideally net of cost recovery.
+
+Group *j*'s benefit is its share of utilisation times the budget. That single
+assumption, that the unit subsidy is the same for everyone using a given level
+of service, is the method's whole weakness and it should be stated in any
+write-up. It is wrong wherever quality varies systematically with who is being
+served, which is most places. It biases the result toward finding spending
+*more* progressive than it is, because a rural PHC and a district hospital are
+counted as delivering the same rupee value.
+
+What this module will not do is compare a benefit share against a population
+share and call the difference targeting. That comparison ignores need, and need
+is not flat across the distribution: the poorest quintile carries more disease
+per person, so an equal share of health spending is already regressive relative
+to need.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+
+from .concentration import concentration_index
+
+__all__ = ["benefit_incidence", "benefit_incidence_by_level"]
+
+
+def benefit_incidence(df: pd.DataFrame, *, utilisation: str, group: str,
+                      spending: float, weights: str | None = None,
+                      fees: str | None = None,
+                      group_order: list | None = None) -> pd.DataFrame:
+    """Allocate a budget across living-standards groups by their use of the service.
+
+    Parameters
+    ----------
+    utilisation
+        Column counting use: visits, enrolled children, admissions. Zero for a
+        household that did not use the service, not missing, or that household
+        drops out of the denominator and every share is overstated.
+    group
+        Living-standards group. A wealth quintile is the usual choice.
+    spending
+        Total public expenditure on the service, in the survey's own currency
+        units. Net of cost recovery if you have that figure; if not, pass
+        ``fees`` and read the net column.
+    fees
+        Optional column of what the household paid out of pocket. Where given,
+        ``net_benefit`` subtracts it, which is the number that answers whether
+        the household came out ahead.
+    group_order
+        Order for the output rows. Defaults to sorted order, which is right for
+        numeric quintiles and wrong for strings like "poorest".
+
+    Returns one row per group: utilisation share, benefit, benefit per person,
+    and the share of the total each group receives.
+
+    >>> import pandas as pd
+    >>> d = pd.DataFrame({"visits": [3, 2, 1, 0], "q": [1, 2, 3, 4]})
+    >>> benefit_incidence(d, utilisation="visits", group="q", spending=600)["benefit"].tolist()
+    [300.0, 200.0, 100.0, 0.0]
+    """
+    for col in [utilisation, group] + ([weights] if weights else []) + ([fees] if fees else []):
+        if col not in df.columns:
+            raise KeyError(f"benefit_incidence: no column {col!r} in the frame")
+    if spending <= 0:
+        raise ValueError("benefit_incidence: spending must be positive")
+
+    data = df.copy()
+    data["_w"] = data[weights].astype(float) if weights else 1.0
+    data["_u"] = pd.to_numeric(data[utilisation], errors="coerce")
+    if data["_u"].isna().any():
+        missing = int(data["_u"].isna().sum())
+        raise ValueError(
+            f"benefit_incidence: {missing} rows have missing utilisation. A household "
+            "that did not use the service should be a zero, not a blank; leaving it "
+            "blank drops it from the denominator and inflates every share.")
+    if (data["_u"] < 0).any():
+        raise ValueError("benefit_incidence: negative utilisation")
+
+    data["_wu"] = data["_w"] * data["_u"]
+    total_u = data["_wu"].sum()
+    if total_u <= 0:
+        raise ValueError("benefit_incidence: nobody in the sample used the service")
+    total_w = data["_w"].sum()
+
+    agg = {"_w": "sum", "_wu": "sum", utilisation: "size"}
+    g = data.groupby(group, dropna=False, observed=True).agg(
+        population=("_w", "sum"), utilisation=("_wu", "sum"), n=(utilisation, "size"))
+    g = g.reset_index()
+
+    g["utilisation_share"] = g["utilisation"] / total_u
+    g["population_share"] = g["population"] / total_w
+    g["benefit"] = g["utilisation_share"] * spending
+    g["benefit_per_person"] = g["benefit"] / g["population"]
+
+    if fees is not None:
+        paid = data.assign(_f=data["_w"] * pd.to_numeric(data[fees], errors="coerce").fillna(0.0))
+        f = paid.groupby(group, dropna=False, observed=True)["_f"].sum().reset_index(name="fees_paid")
+        g = g.merge(f, on=group, how="left")
+        g["net_benefit"] = g["benefit"] - g["fees_paid"]
+        g["net_benefit_per_person"] = g["net_benefit"] / g["population"]
+
+    if group_order is not None:
+        g[group] = pd.Categorical(g[group], categories=group_order, ordered=True)
+        g = g.sort_values(group)
+    else:
+        g = g.sort_values(group)
+
+    # The concentration index of the benefit, computed on the household rows so
+    # that within-group variation in utilisation is not thrown away by the
+    # groupby. This is the one number that says progressive or regressive.
+    ci = concentration_index(data["_u"], rank_by=data[group], weights=data["_w"])
+    g.attrs["concentration_index"] = ci
+    g.attrs["verdict"] = (
+        "progressive: the benefit is concentrated among the poorer groups" if ci < -0.01
+        else "regressive: the benefit is concentrated among the better-off" if ci > 0.01
+        else "roughly proportional across groups")
+    g.attrs["spending"] = float(spending)
+    g.attrs["caveat"] = (
+        "Assumes an identical unit subsidy for every user of this service. Where "
+        "quality rises with the wealth of the catchment, this understates how "
+        "regressive the spending is.")
+    return g.reset_index(drop=True)
+
+
+def benefit_incidence_by_level(df: pd.DataFrame, *, utilisation: str, group: str,
+                               level: str, spending: dict, weights: str | None = None,
+                               group_order: list | None = None) -> dict:
+    """Benefit incidence run separately per service level, then aggregated.
+
+    The level split is where the interesting result usually lives. Aggregate
+    health spending often looks close to proportional while primary care is
+    strongly progressive and tertiary care strongly regressive, the two
+    cancelling in the total. Running the whole budget as one number hides
+    exactly the finding a ministry can act on.
+
+    Parameters
+    ----------
+    level
+        Column naming the service level: primary / secondary / tertiary, or
+        PHC / CHC / district hospital.
+    spending
+        Mapping from each level to its budget. Every level present in the data
+        must appear, because a missing one is silently allocated nothing and the
+        aggregate shares then do not mean what they say.
+
+    Returns a dict with ``by_level`` (a frame per level), ``total`` (the
+    aggregated frame) and ``concentration_index`` per level plus overall.
+    """
+    if level not in df.columns:
+        raise KeyError(f"benefit_incidence_by_level: no column {level!r} in the frame")
+    present = set(pd.unique(df[level].dropna()))
+    missing = present - set(spending)
+    if missing:
+        raise ValueError(
+            f"benefit_incidence_by_level: no budget given for level(s) {sorted(missing)}. "
+            "A level with no entry would be allocated nothing and the aggregate shares "
+            "would not sum to the budget you think they do.")
+
+    per_level, frames, cis = {}, [], {}
+    for lv, part in df.groupby(level, dropna=True, observed=True):
+        out = benefit_incidence(part, utilisation=utilisation, group=group,
+                                spending=spending[lv], weights=weights,
+                                group_order=group_order)
+        per_level[lv] = out
+        cis[lv] = out.attrs["concentration_index"]
+        frames.append(out[[group, "benefit", "population"]].assign(**{level: lv}))
+
+    stacked = pd.concat(frames, ignore_index=True)
+    total = stacked.groupby(group, observed=True).agg(
+        benefit=("benefit", "sum"), population=("population", "max")).reset_index()
+    total["benefit_share"] = total["benefit"] / total["benefit"].sum()
+    total["population_share"] = total["population"] / total["population"].sum()
+    total["benefit_per_person"] = total["benefit"] / total["population"]
+    if group_order is not None:
+        total[group] = pd.Categorical(total[group], categories=group_order, ordered=True)
+        total = total.sort_values(group).reset_index(drop=True)
+
+    overall = concentration_index(total["benefit"] / total["population"],
+                                  rank_by=total[group], weights=total["population"])
+    return {"by_level": per_level, "total": total,
+            "concentration_index_by_level": cis,
+            "concentration_index_overall": overall,
+            "spending": dict(spending)}
diff --git a/inequality/indices.py b/inequality/indices.py
new file mode 100644
index 0000000..2ec5700
--- /dev/null
+++ b/inequality/indices.py
@@ -0,0 +1,293 @@
+"""Summary measures of inequality in a continuous variable, survey-weighted.
+
+Consumption, income, land, a school-quality score: anything where the question
+is how unequally the total is spread. Every function here takes weights, drops
+missing rows pairwise, and is scale invariant, so the answer does not depend on
+whether the variable is in rupees, thousands of rupees or 2011 PPP dollars.
+
+Which measure to reach for, since they are not interchangeable:
+
+- **Gini** for the headline. Comparable with almost every published figure,
+  and most sensitive around the middle of the distribution.
+- **Theil (GE(1)) or mean log deviation (GE(0))** when the question is *where*
+  inequality sits, because both decompose exactly into a within-group and a
+  between-group part. `decomposition.theil_decomposition` does that.
+- **Atkinson** when you want to state the aversion to inequality explicitly
+  rather than let the index choose one for you. A(0.5) and A(2) answer
+  different normative questions about the same data.
+- **Palma and the 80/20 ratio** when the audience is not technical. Both are
+  ratios of shares, so they are read off a table without a formula.
+
+None of these say anything about *who* is where. For that, an inequality
+measure has to be crossed with a living-standards rank, which is the
+concentration index in `concentration.py`.
+
+Sources for the definitions: Cowell, *Measuring Inequality* (3rd edn, OUP 2011),
+chapters 2 and 3; Haughton and Khandker, *Handbook on Poverty and Inequality*
+(World Bank 2009), chapter 6.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+
+from .ranks import fractional_rank, prepare
+
+__all__ = [
+    "gini", "lorenz_curve", "share_of_top", "share_of_bottom", "quantile_shares",
+    "palma_ratio", "ratio_80_20", "generalised_entropy", "theil_t", "theil_l",
+    "atkinson", "weighted_quantile",
+]
+
+
+def gini(values, weights=None, *, small_sample_correction: bool = False) -> float:
+    """Gini coefficient, 0 at perfect equality and 1 at perfect concentration.
+
+    Computed as ``2 · cov(y, F(y)) / ȳ`` using midpoint fractional ranks, which
+    is the population (not sample) Gini and matches the mean-difference
+    definition ``Δ / 2ȳ`` exactly.
+
+    Parameters
+    ----------
+    small_sample_correction
+        Multiply by ``n / (n - 1)``. The population formula is biased downward
+        in a small sample. Leave it off to compare against published national
+        figures, which are almost always uncorrected; turn it on for a survey of
+        a few dozen units where the bias is not negligible. With weights, *n* is
+        the number of rows, not the weighted total.
+
+    Negative values raise. A Gini on a variable that can go negative (net income
+    after debt, farm profit in a bad year) is not bounded by 1 and is not
+    comparable with anything; the fix is a decision about the data, not a
+    silently different formula.
+
+    >>> round(gini([1, 2, 3, 4, 5]), 6)
+    0.266667
+    >>> round(gini([10, 10, 10]), 12)
+    0.0
+    """
+    v, w = prepare(values, weights, name="gini")
+    if np.any(v < 0):
+        raise ValueError(
+            "gini: the variable contains negative values. The Gini is not bounded "
+            "by 1 on a variable that can go negative, and the result would not be "
+            "comparable with any published figure. Decide what a negative means "
+            "(censor at zero, drop, or use a different measure) before calling this.")
+    mean = np.sum(v * w) / np.sum(w)
+    if mean == 0:
+        return 0.0
+    r = fractional_rank(v, w)
+    g = 2.0 * np.sum(w * v * r) / (np.sum(w) * mean) - 1.0
+    if small_sample_correction:
+        n = v.size
+        if n < 2:
+            raise ValueError("gini: the small-sample correction needs at least 2 rows")
+        g *= n / (n - 1)
+    return float(g)
+
+
+def lorenz_curve(values, weights=None) -> pd.DataFrame:
+    """Points on the Lorenz curve: cumulative population share against value share.
+
+    Returns a frame with ``population_share`` and ``value_share``, opening with
+    the origin (0, 0) and closing at (1, 1), ordered from poorest to richest.
+    One row per distinct value, so a wealth-quintile variable gives six rows
+    rather than a row per household.
+
+    The area between this curve and the 45 degree line is half the Gini, which
+    `test_inequality.py` checks by trapezoid against `gini` directly.
+    """
+    v, w = prepare(values, weights, name="lorenz_curve")
+    if np.any(v < 0):
+        raise ValueError("lorenz_curve: negative values make the curve non-monotone")
+    order = np.argsort(v, kind="mergesort")
+    sv, sw = v[order], w[order]
+
+    starts = np.concatenate(([True], sv[1:] != sv[:-1]))
+    group = np.cumsum(starts) - 1
+    gw = np.zeros(group[-1] + 1)
+    gv = np.zeros(group[-1] + 1)
+    np.add.at(gw, group, sw)
+    np.add.at(gv, group, sw * sv)
+
+    total_w, total_v = gw.sum(), gv.sum()
+    pop = np.concatenate(([0.0], np.cumsum(gw) / total_w))
+    val = np.concatenate(([0.0], np.cumsum(gv) / total_v)) if total_v > 0 else pop.copy()
+    return pd.DataFrame({"population_share": pop, "value_share": val})
+
+
+def _share(values, weights, lower: float, upper: float) -> float:
+    """Share of the total held by the population between two rank cut-points.
+
+    The cut-point almost never falls on a row boundary once weights are
+    involved, so the marginal row is split proportionally rather than assigned
+    whole to one side. Assigning it whole is the usual shortcut and it makes the
+    top-decile share jump around when a single large-weight household crosses
+    the line.
+    """
+    v, w = prepare(values, weights, name="share")
+    order = np.argsort(v, kind="mergesort")
+    sv, sw = v[order], w[order]
+    total_w = sw.sum()
+    total_v = np.sum(sv * sw)
+    if total_v == 0:
+        return 0.0
+
+    cum = np.cumsum(sw) / total_w
+    start = np.concatenate(([0.0], cum[:-1]))
+    frac = np.clip((np.minimum(cum, upper) - np.maximum(start, lower))
+                   / np.where(cum - start > 0, cum - start, 1.0), 0.0, 1.0)
+    return float(np.sum(frac * sw * sv) / total_v)
+
+
+def share_of_top(values, weights=None, *, p: float = 0.10) -> float:
+    """Share of the total held by the richest ``p`` of the population."""
+    if not 0 < p <= 1:
+        raise ValueError("share_of_top: p must be in (0, 1]")
+    return _share(values, weights, 1.0 - p, 1.0)
+
+
+def share_of_bottom(values, weights=None, *, p: float = 0.40) -> float:
+    """Share of the total held by the poorest ``p`` of the population."""
+    if not 0 < p <= 1:
+        raise ValueError("share_of_bottom: p must be in (0, 1]")
+    return _share(values, weights, 0.0, p)
+
+
+def quantile_shares(values, weights=None, *, q: int = 5) -> pd.DataFrame:
+    """Share of the total held by each equal-sized population group.
+
+    ``q=5`` gives quintiles, ``q=10`` deciles. Groups are equal shares of
+    *population*, cut with the marginal observation split, so the shares sum to
+    exactly 1 whatever the weights look like.
+    """
+    if q < 2:
+        raise ValueError("quantile_shares: q must be at least 2")
+    edges = np.linspace(0, 1, q + 1)
+    rows = [{"group": i + 1,
+             "population_share": edges[i + 1] - edges[i],
+             "value_share": _share(values, weights, edges[i], edges[i + 1])}
+            for i in range(q)]
+    out = pd.DataFrame(rows)
+    out["ratio_to_equal_share"] = out["value_share"] * q
+    return out
+
+
+def palma_ratio(values, weights=None) -> float:
+    """Richest 10 per cent's share divided by the poorest 40 per cent's.
+
+    Gabriel Palma's observation is that the middle five deciles take a stable
+    half of national income almost everywhere, so the distributional fight is
+    between the top decile and the bottom four. At perfect equality the ratio is
+    0.25; India's is usually reported somewhere above 1.5.
+    """
+    bottom = share_of_bottom(values, weights, p=0.40)
+    if bottom == 0:
+        return float("inf")
+    return share_of_top(values, weights, p=0.10) / bottom
+
+
+def ratio_80_20(values, weights=None) -> float:
+    """Richest fifth's share over the poorest fifth's (Eurostat's S80/S20)."""
+    bottom = share_of_bottom(values, weights, p=0.20)
+    if bottom == 0:
+        return float("inf")
+    return share_of_top(values, weights, p=0.20) / bottom
+
+
+def generalised_entropy(values, weights=None, *, alpha: float = 1.0) -> float:
+    """Generalised entropy index GE(α).
+
+    ``alpha`` sets where the index is sensitive. GE(0), the mean log deviation,
+    weights the bottom of the distribution; GE(1), Theil's T, is neutral;
+    GE(2), half the squared coefficient of variation, weights the top. Only
+    GE(0) and GE(1) decompose additively into within and between components,
+    which is why those two have their own names below.
+
+    Zero and negative values raise for α ≤ 1, where the logarithm or the
+    negative power is undefined. That is a real constraint on consumption data
+    with zeros in it, and it is better as an error than as a NaN in a table.
+    """
+    v, w = prepare(values, weights, name="generalised_entropy")
+    mean = np.sum(v * w) / np.sum(w)
+    if mean <= 0:
+        raise ValueError("generalised_entropy: the weighted mean must be positive")
+    if alpha <= 1 and np.any(v <= 0):
+        raise ValueError(
+            f"generalised_entropy: GE({alpha:g}) is undefined where the variable is "
+            "zero or negative. Either drop those rows and say so, or use alpha > 1.")
+    if np.any(v < 0):
+        raise ValueError("generalised_entropy: negative values are undefined here")
+
+    p = w / np.sum(w)
+    ratio = v / mean
+    if np.isclose(alpha, 0.0):
+        return float(np.sum(p * np.log(1.0 / ratio)))
+    if np.isclose(alpha, 1.0):
+        with np.errstate(divide="ignore", invalid="ignore"):
+            term = np.where(ratio > 0, ratio * np.log(ratio), 0.0)
+        return float(np.sum(p * term))
+    return float(np.sum(p * (ratio ** alpha - 1.0)) / (alpha * (alpha - 1.0)))
+
+
+def theil_t(values, weights=None) -> float:
+    """Theil's T index, GE(1). Decomposes by income share."""
+    return generalised_entropy(values, weights, alpha=1.0)
+
+
+def theil_l(values, weights=None) -> float:
+    """Mean log deviation, GE(0), sometimes Theil's L. Decomposes by population share."""
+    return generalised_entropy(values, weights, alpha=0.0)
+
+
+def atkinson(values, weights=None, *, epsilon: float = 1.0) -> float:
+    """Atkinson index at inequality aversion ``epsilon``.
+
+    Reads as a fraction of total income society would be willing to give up to
+    have the remainder distributed equally. A(0) is 0 by construction: no
+    aversion, no cost. A(2) is the usual high-aversion figure. The index is
+    weakly increasing in ``epsilon``, which `test_inequality.py` checks.
+
+    Undefined at ``epsilon`` ≥ 1 with any zero in the data: one destitute
+    household drives the equally-distributed-equivalent income to zero and the
+    index to 1, correctly but uninformatively. That case raises.
+    """
+    if epsilon < 0:
+        raise ValueError("atkinson: epsilon must be non-negative")
+    v, w = prepare(values, weights, name="atkinson")
+    if np.any(v < 0):
+        raise ValueError("atkinson: negative values are undefined here")
+    mean = np.sum(v * w) / np.sum(w)
+    if mean <= 0:
+        raise ValueError("atkinson: the weighted mean must be positive")
+    if epsilon >= 1 and np.any(v <= 0):
+        raise ValueError(
+            f"atkinson: at epsilon = {epsilon:g} a single zero forces the index to 1. "
+            "Drop or censor the zeros and say which you did.")
+
+    p = w / np.sum(w)
+    ratio = v / mean
+    if np.isclose(epsilon, 0.0):
+        return 0.0
+    if np.isclose(epsilon, 1.0):
+        return float(1.0 - np.exp(np.sum(p * np.log(ratio))))
+    ede = np.sum(p * ratio ** (1.0 - epsilon)) ** (1.0 / (1.0 - epsilon))
+    return float(1.0 - ede)
+
+
+def weighted_quantile(values, q, weights=None) -> float:
+    """Value at weighted quantile ``q``, interpolated on the midpoint rank.
+
+    Used for cut-points (a poverty line at the 20th percentile, a top-decile
+    threshold). Consistent with `fractional_rank`, so a household at the
+    returned value has rank ``q`` by this package's definition rather than by
+    numpy's, which handles ties and weights differently.
+    """
+    if not 0 <= q <= 1:
+        raise ValueError("weighted_quantile: q must be in [0, 1]")
+    v, w = prepare(values, weights, name="weighted_quantile")
+    order = np.argsort(v, kind="mergesort")
+    sv, sw = v[order], w[order]
+    cum = (np.cumsum(sw) - 0.5 * sw) / sw.sum()
+    return float(np.interp(q, cum, sv))
diff --git a/inequality/ranks.py b/inequality/ranks.py
new file mode 100644
index 0000000..8110c1c
--- /dev/null
+++ b/inequality/ranks.py
@@ -0,0 +1,128 @@
+"""Weighted fractional ranks, the shared foundation of every index here.
+
+Almost every distributional statistic in this package is a covariance between a
+variable and somebody's *position* in a distribution. The position is a weighted
+fractional rank: the share of the population poorer than you, plus half your own
+share. Get this wrong and the Gini coefficient, the concentration index and the
+Lorenz curve are all wrong together, in ways that look plausible.
+
+Two decisions are made here rather than in each index.
+
+**Ties take the midpoint of the tie group.** If four hundred households sit in
+wealth quintile 2, every one of them gets the same rank: the midpoint of the
+block of population that quintile 2 occupies. The alternative, breaking ties by
+row order, silently makes the answer depend on how the file was sorted, and
+wealth quintiles are exactly the kind of coarse variable this package is pointed
+at. Stata's `conindex` and the World Bank's ADePT do the same thing.
+
+**Weights are relative, not absolute.** Everything divides by the total weight,
+so a DHS `v005` normalised to average one and a PLFS multiplier expanding to
+1.4 billion give the same answer. That is the property that lets these functions
+be handed either without a footnote.
+
+The midpoint rank has one arithmetic property worth stating, because three
+proofs downstream lean on it: the weighted mean rank is exactly 0.5, for any
+weights and any pattern of ties. `test_inequality.py` asserts it.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+
+__all__ = ["fractional_rank", "prepare", "weighted_mean"]
+
+
+def prepare(values, weights=None, *, name: str = "values"):
+    """Coerce a value/weight pair to clean float arrays, dropping missing rows.
+
+    Returns ``(values, weights)`` with any row missing either one removed. A
+    missing value is not a zero: a household that did not report consumption is
+    absent from the distribution, and treating it as destitute would push every
+    inequality measure up. Rows are dropped pairwise and the caller can compare
+    lengths if it needs to report how many.
+
+    Raises on a negative weight, on a weight total of zero, and on an empty
+    result, because each of those produces a number rather than an error further
+    down.
+    """
+    v = np.asarray(pd.to_numeric(pd.Series(values), errors="coerce"), dtype=float)
+    if weights is None:
+        w = np.ones_like(v, dtype=float)
+    else:
+        w = np.asarray(pd.to_numeric(pd.Series(weights), errors="coerce"), dtype=float)
+        if w.shape != v.shape:
+            raise ValueError(
+                f"{name} has {v.shape[0]} rows but weights has {w.shape[0]}")
+
+    keep = np.isfinite(v) & np.isfinite(w)
+    v, w = v[keep], w[keep]
+
+    if v.size == 0:
+        raise ValueError(f"{name}: no rows left after dropping missing values")
+    if np.any(w < 0):
+        raise ValueError(f"{name}: negative weights are not a thing; check the weight column")
+    total = w.sum()
+    if total <= 0:
+        raise ValueError(f"{name}: weights sum to {total}, so nothing can be estimated")
+    return v, w
+
+
+def weighted_mean(values, weights=None) -> float:
+    """Weighted arithmetic mean, with missing rows dropped pairwise."""
+    v, w = prepare(values, weights)
+    return float(np.sum(v * w) / np.sum(w))
+
+
+def fractional_rank(values, weights=None) -> np.ndarray:
+    """Weighted fractional rank in [0, 1], ties sharing the group midpoint.
+
+    For observation *i* the rank is the population share strictly poorer than
+    *i*, plus half the share at *i*'s own value. Returned in the caller's row
+    order, so it can be assigned straight back onto a data frame column.
+
+    Missing values rank as NaN rather than sorting to one end, which is the
+    behaviour you want when the living-standards variable has gaps: those rows
+    then fall out of the covariance instead of being counted as poorest.
+
+    >>> fractional_rank([1, 2, 3, 4, 5]).round(2)
+    array([0.1, 0.3, 0.5, 0.7, 0.9])
+    >>> fractional_rank([1, 1, 2]).round(4)          # the tie shares a rank
+    array([0.3333, 0.3333, 0.8333])
+    """
+    v = np.asarray(pd.to_numeric(pd.Series(values), errors="coerce"), dtype=float)
+    if weights is None:
+        w = np.ones_like(v, dtype=float)
+    else:
+        w = np.asarray(pd.to_numeric(pd.Series(weights), errors="coerce"), dtype=float)
+    ok = np.isfinite(v) & np.isfinite(w)
+    if not ok.any():
+        raise ValueError("fractional_rank: every row is missing a value or a weight")
+    if np.any(w[ok] < 0):
+        raise ValueError("fractional_rank: negative weights are not a thing")
+
+    out = np.full(v.shape, np.nan, dtype=float)
+    vv, ww = v[ok], w[ok]
+    total = ww.sum()
+    if total <= 0:
+        raise ValueError("fractional_rank: weights sum to zero")
+
+    order = np.argsort(vv, kind="mergesort")     # stable, so ties keep input order
+    sv, sw = vv[order], ww[order]
+
+    cum_before = np.concatenate(([0.0], np.cumsum(sw)[:-1]))
+
+    # Collapse ties: every row sharing a value gets the midpoint of the whole
+    # block that value occupies, so the answer cannot depend on row order.
+    starts = np.concatenate(([True], sv[1:] != sv[:-1]))
+    group = np.cumsum(starts) - 1
+    block_start = cum_before[starts]              # weight strictly below each value
+    block_weight = np.zeros(group[-1] + 1)
+    np.add.at(block_weight, group, sw)            # total weight sharing each value
+
+    ranks_sorted = (block_start[group] + block_weight[group] / 2.0) / total
+
+    inv = np.empty_like(order)
+    inv[order] = np.arange(order.size)
+    out[ok] = ranks_sorted[inv]
+    return out
diff --git a/social_sector/composite_index.py b/social_sector/composite_index.py
new file mode 100644
index 0000000..ca10611
--- /dev/null
+++ b/social_sector/composite_index.py
@@ -0,0 +1,253 @@
+"""Composite indices: combining several indicators into one score, defensibly.
+
+Every access index, vulnerability index and readiness index is the same three
+decisions, and almost every published one leaves at least two of them implicit:
+
+1. **Direction.** Infant mortality and immunisation coverage both measure child
+   health, and they point opposite ways. Averaging them raw cancels the signal.
+2. **Normalisation.** Indicators arrive in incompatible units: a rate per
+   thousand, a percentage, a rupee figure. How you put them on a common scale
+   decides the answer, and min-max, z-score and goalpost give different
+   rankings on the same data.
+3. **Weights.** Equal weighting is a choice, not the absence of one. It says a
+   percentage point of literacy is worth a percentage point of piped water.
+
+This module makes all three arguments rather than assumptions, records what was
+chosen on the result, and refuses the combinations that produce a number
+somebody will quote and nobody can defend.
+
+The predecessor of this file computed ``df[cols].mean(axis=1)`` and called the
+result an access index. That is min-max-free, direction-free, weight-free, and
+it silently ranks a district by whichever of its indicators happens to be
+measured on the largest scale.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+
+__all__ = ["normalise", "composite_index", "rank_sensitivity"]
+
+_METHODS = ("minmax", "zscore", "goalpost", "rank")
+
+
+def normalise(series, *, method: str = "minmax", direction: str = "higher_is_better",
+              goalposts: tuple[float, float] | None = None) -> pd.Series:
+    """Put one indicator on a common 0-1 scale (or a z-score).
+
+    Parameters
+    ----------
+    method
+        ``minmax``   rescales to [0, 1] using the observed range. Simple, and
+                     the answer moves when you add a new district, so a series
+                     of annual indices built this way is not comparable
+                     year to year.
+        ``zscore``   standardises to mean 0, sd 1. Keeps relative distances,
+                     unbounded, and sensitive to outliers.
+        ``goalpost`` rescales against fixed minimum and maximum values you
+                     supply, the way the UNDP's Human Development Index does.
+                     The only one of the four that is comparable across time
+                     and across samples, because the scale does not depend on
+                     the data in hand.
+        ``rank``     replaces values with their fractional rank. Throws away
+                     magnitude, immune to outliers and to skew.
+    direction
+        ``higher_is_better`` or ``lower_is_better``. For a mortality rate or a
+        distance to the nearest facility, use the latter and the series is
+        flipped so that a high score always means a good outcome.
+
+    Returns a Series aligned to the input, NaN preserved.
+    """
+    if method not in _METHODS:
+        raise ValueError(f"normalise: method must be one of {_METHODS}")
+    if direction not in ("higher_is_better", "lower_is_better"):
+        raise ValueError("normalise: direction must be 'higher_is_better' or 'lower_is_better'")
+
+    s = pd.to_numeric(pd.Series(series), errors="coerce").astype(float)
+    if s.notna().sum() == 0:
+        raise ValueError("normalise: the indicator is entirely missing")
+
+    if method == "goalpost":
+        if goalposts is None:
+            raise ValueError(
+                "normalise: method='goalpost' needs explicit goalposts. That is the "
+                "point of it: the scale must not depend on the sample, or the index "
+                "is not comparable with the same index computed next year.")
+        lo, hi = goalposts
+        if not hi > lo:
+            raise ValueError("normalise: goalposts must satisfy max > min")
+        out = (s - lo) / (hi - lo)
+        out = out.clip(0.0, 1.0)
+    elif method == "minmax":
+        lo, hi = s.min(), s.max()
+        if hi == lo:
+            out = pd.Series(np.where(s.notna(), 0.5, np.nan), index=s.index)
+        else:
+            out = (s - lo) / (hi - lo)
+    elif method == "zscore":
+        sd = s.std(ddof=0)
+        out = (s - s.mean()) / sd if sd > 0 else pd.Series(
+            np.where(s.notna(), 0.0, np.nan), index=s.index)
+    else:  # rank
+        out = s.rank(pct=True, na_option="keep")
+
+    if direction == "lower_is_better":
+        out = -out if method == "zscore" else 1.0 - out
+    return out
+
+
+def composite_index(df: pd.DataFrame, indicators, *, method: str = "minmax",
+                    weights: dict | None = None,
+                    directions: dict | None = None,
+                    goalposts: dict | None = None,
+                    aggregation: str = "arithmetic",
+                    min_indicators: int | None = None,
+                    name: str = "index") -> pd.DataFrame:
+    """Combine indicators into one score, keeping every choice on the record.
+
+    Parameters
+    ----------
+    indicators
+        Column names to combine.
+    weights
+        ``{column: weight}``. Rescaled to sum to 1, so relative sizes are what
+        matter. Omitted columns get the residual share equally. Default: equal.
+    directions
+        ``{column: 'lower_is_better'}`` for any indicator where a high value is
+        a bad outcome. Anything unlisted is treated as higher-is-better.
+    aggregation
+        ``arithmetic`` for a weighted mean: indicators substitute freely, so a
+        district can offset no piped water with excellent schools.
+        ``geometric`` for a weighted geometric mean: they substitute poorly, so
+        a near-zero on any one component drags the whole index down. The UNDP
+        moved the HDI to geometric in 2010 for exactly that reason, and for an
+        *access* index it is usually the more honest choice, because a household
+        with no clinic within 40 km is not compensated by a nearby school.
+    min_indicators
+        Rows with fewer than this many non-missing indicators get NaN rather
+        than a score built from whatever happened to be present. Defaults to
+        all of them. This matters more than it sounds: a district reporting two
+        of six indicators, both good, otherwise scores near the top.
+
+    Returns the input frame with the score, the rank, and a per-indicator
+    normalised column appended. The choices are recorded in ``.attrs`` so a
+    saved result carries its own method statement.
+    """
+    indicators = list(indicators)
+    missing = [c for c in indicators if c not in df.columns]
+    if missing:
+        raise KeyError(f"composite_index: no column(s) {missing} in the frame")
+    if len(indicators) < 2:
+        raise ValueError("composite_index: combining fewer than two indicators is not an index")
+    if aggregation not in ("arithmetic", "geometric"):
+        raise ValueError("composite_index: aggregation must be 'arithmetic' or 'geometric'")
+
+    directions = dict(directions or {})
+    goalposts = dict(goalposts or {})
+    bad = set(directions) - set(indicators)
+    if bad:
+        raise KeyError(f"composite_index: directions names column(s) not in indicators: {sorted(bad)}")
+
+    if weights is None:
+        w = {c: 1.0 / len(indicators) for c in indicators}
+    else:
+        bad = set(weights) - set(indicators)
+        if bad:
+            raise KeyError(f"composite_index: weights names column(s) not in indicators: {sorted(bad)}")
+        if any(v < 0 for v in weights.values()):
+            raise ValueError("composite_index: negative weights")
+        named = sum(weights.values())
+        if named <= 0:
+            raise ValueError("composite_index: weights sum to zero")
+        rest = [c for c in indicators if c not in weights]
+        if rest:
+            residual = max(0.0, 1.0 - named)
+            if residual == 0:
+                raise ValueError(
+                    f"composite_index: the named weights already sum to {named:g}, leaving "
+                    f"nothing for {rest}. Name every indicator, or leave weights as None.")
+            w = {**weights, **{c: residual / len(rest) for c in rest}}
+        else:
+            w = dict(weights)
+        total = sum(w.values())
+        w = {k: v / total for k, v in w.items()}
+
+    out = df.copy()
+    norm_cols = []
+    for c in indicators:
+        col = f"{c}__norm"
+        out[col] = normalise(df[c], method=method,
+                             direction=directions.get(c, "higher_is_better"),
+                             goalposts=goalposts.get(c))
+        norm_cols.append(col)
+
+    normed = out[norm_cols]
+    present = normed.notna().sum(axis=1)
+    required = len(indicators) if min_indicators is None else min_indicators
+    if required > len(indicators):
+        raise ValueError("composite_index: min_indicators exceeds the number of indicators")
+
+    wvec = np.array([w[c] for c in indicators])
+    mask = normed.notna().to_numpy()
+    vals = normed.to_numpy(dtype=float)
+
+    # Reweight over the indicators actually present in each row, so a row with a
+    # gap is scored on the same 0-1 scale rather than being pushed down.
+    eff = np.where(mask, wvec, 0.0)
+    denom = eff.sum(axis=1)
+
+    if aggregation == "arithmetic":
+        score = np.where(denom > 0, np.nansum(np.where(mask, vals * wvec, 0.0), axis=1) / denom, np.nan)
+    else:
+        if method == "zscore":
+            raise ValueError(
+                "composite_index: a geometric mean of z-scores is undefined, because "
+                "half of them are negative. Use method='minmax', 'goalpost' or 'rank'.")
+        floor = 1e-9   # a true zero would send the whole index to zero
+        logv = np.where(mask, np.log(np.clip(vals, floor, None)), 0.0)
+        score = np.where(denom > 0, np.exp(np.nansum(logv * eff, axis=1) / denom), np.nan)
+
+    score = np.where(present.to_numpy() >= required, score, np.nan)
+    out[name] = score
+    out[f"{name}_rank"] = pd.Series(score, index=out.index).rank(ascending=False, method="min")
+    out[f"{name}_n_indicators"] = present
+
+    out.attrs["composite_index"] = {
+        "indicators": indicators, "weights": w, "method": method,
+        "aggregation": aggregation, "directions": directions,
+        "goalposts": goalposts, "min_indicators": required,
+    }
+    return out
+
+
+def rank_sensitivity(df: pd.DataFrame, indicators, *, unit: str,
+                     methods=("minmax", "zscore", "rank"),
+                     aggregations=("arithmetic", "geometric"),
+                     **kwargs) -> pd.DataFrame:
+    """How much the ranking moves when the normalisation choice moves.
+
+    Run this before publishing any index. A district whose rank swings from 3rd
+    to 24th depending on whether you min-max or z-score has not been measured,
+    it has been assigned, and the sensible thing is to report a band rather than
+    a position. Composite indices get quoted as facts and the arbitrariness in
+    them is invisible once the number is in a headline.
+
+    Returns one row per unit with the rank under each combination, plus the
+    spread. ``zscore`` is skipped for the geometric aggregation, where it is
+    undefined.
+    """
+    if unit not in df.columns:
+        raise KeyError(f"rank_sensitivity: no column {unit!r} in the frame")
+    ranks = pd.DataFrame({unit: df[unit].to_numpy()})
+    for m in methods:
+        for a in aggregations:
+            if a == "geometric" and m == "zscore":
+                continue
+            r = composite_index(df, indicators, method=m, aggregation=a, **kwargs)
+            ranks[f"{m}_{a}"] = r["index_rank"].to_numpy()
+    cols = [c for c in ranks.columns if c != unit]
+    ranks["best_rank"] = ranks[cols].min(axis=1)
+    ranks["worst_rank"] = ranks[cols].max(axis=1)
+    ranks["rank_spread"] = ranks["worst_rank"] - ranks["best_rank"]
+    return ranks.sort_values("rank_spread", ascending=False).reset_index(drop=True)
diff --git a/social_sector/public_health_access_index.py b/social_sector/public_health_access_index.py
index 2e1a4ea..5b6c0d9 100644
--- a/social_sector/public_health_access_index.py
+++ b/social_sector/public_health_access_index.py
@@ -1,6 +1,71 @@
-"""Generate simple binary access index from multiple variables."""
+"""Public health access index: a thin wrapper over the composite machinery.
+
+Kept because two tests and any existing script import `compute_access_index`
+from here. It now delegates to `social_sector.composite_index`, which makes the
+three decisions an index needs (direction, normalisation, weights) explicit
+instead of assuming them.
+
+The original was a bare row mean over the raw columns. That is fine when every
+indicator is already a 0-1 coverage share pointing the same way, which is the
+case this function still serves, and wrong the moment one of them is a mortality
+rate or a rupee figure, which is why the default now normalises.
+"""
+
+from __future__ import annotations
+
 import pandas as pd
 
-def compute_access_index(df, access_columns):
-    df['access_index'] = df[access_columns].mean(axis=1)
-    return df
\ No newline at end of file
+from .composite_index import composite_index
+
+__all__ = ["compute_access_index"]
+
+
+def compute_access_index(df, access_columns, *, normalise: bool = False,
+                         directions: dict | None = None,
+                         weights: dict | None = None,
+                         column: str = "access_index") -> pd.DataFrame:
+    """Mean access across several indicators, appended as a new column.
+
+    Parameters
+    ----------
+    normalise
+        ``False`` (the default, and the original behaviour) averages the
+        columns as they stand. Correct only where every indicator is already on
+        the same scale and pointing the same way, which for a set of coverage
+        percentages it is.
+        ``True`` min-max normalises each indicator first, which is what you
+        want the moment the units differ.
+    directions
+        ``{column: 'lower_is_better'}`` for any indicator where high is bad.
+        Requires ``normalise=True``, because flipping a raw unnormalised
+        indicator has no defined meaning.
+
+    For anything beyond this, call `composite_index` directly: it offers
+    goalpost normalisation for cross-year comparability, geometric aggregation
+    where indicators should not substitute for each other, and a minimum
+    non-missing count so a district reporting two of six indicators is not
+    scored as though it reported all six.
+
+    >>> import pandas as pd
+    >>> d = pd.DataFrame({"toilet": [0.8, 0.4], "water": [0.6, 0.2]})
+    >>> compute_access_index(d, ["toilet", "water"])["access_index"].tolist()
+    [0.7000000000000001, 0.30000000000000004]
+    """
+    df = pd.DataFrame(df)
+    missing = [c for c in access_columns if c not in df.columns]
+    if missing:
+        raise KeyError(f"compute_access_index: no column(s) {missing} in the frame")
+
+    if directions and not normalise:
+        raise ValueError(
+            "compute_access_index: directions needs normalise=True. Flipping a raw "
+            "indicator that has not been put on a 0-1 scale has no defined meaning.")
+
+    if not normalise:
+        out = df.copy()
+        out[column] = df[list(access_columns)].mean(axis=1)
+        return out
+
+    scored = composite_index(df, list(access_columns), method="minmax",
+                             directions=directions, weights=weights, name=column)
+    return scored.drop(columns=[c for c in scored.columns if c.endswith("__norm")])
diff --git a/tests/test_composite_index.py b/tests/test_composite_index.py
new file mode 100644
index 0000000..7f1b82c
--- /dev/null
+++ b/tests/test_composite_index.py
@@ -0,0 +1,211 @@
+"""Tests for social_sector.composite_index.
+
+The failure mode a composite index has is not an exception, it is a ranking
+that looks authoritative and would have come out differently under an equally
+defensible choice. So these tests check the choices bite: that direction flips
+the contribution, that weights move the answer by the amount they should, that
+a row with two of six indicators does not quietly score near the top, and that
+`rank_sensitivity` actually reports movement when movement exists.
+"""
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from social_sector.composite_index import composite_index, normalise, rank_sensitivity
+
+
+@pytest.fixture
+def districts():
+    return pd.DataFrame({
+        "district": list("ABCDE"),
+        "imm_pct": [95.0, 80.0, 60.0, 40.0, 20.0],
+        "imr": [12.0, 20.0, 35.0, 50.0, 70.0],       # lower is better
+        "piped_pct": [90.0, 70.0, 50.0, 30.0, 5.0],
+    })
+
+
+# --------------------------------------------------------------------------
+# normalise
+# --------------------------------------------------------------------------
+
+def test_minmax_puts_the_extremes_at_zero_and_one():
+    out = normalise([10.0, 20.0, 30.0])
+    assert out.tolist() == pytest.approx([0.0, 0.5, 1.0])
+
+
+def test_lower_is_better_flips_the_scale():
+    out = normalise([10.0, 20.0, 30.0], direction="lower_is_better")
+    assert out.tolist() == pytest.approx([1.0, 0.5, 0.0])
+
+
+def test_zscore_lower_is_better_negates_rather_than_subtracting_from_one():
+    """1 - z would shift the mean to 1; a z-score's neutral point is 0."""
+    out = normalise([10.0, 20.0, 30.0], method="zscore", direction="lower_is_better")
+    assert out.mean() == pytest.approx(0.0)
+    assert out.iloc[0] > 0 and out.iloc[-1] < 0
+
+
+def test_goalposts_do_not_depend_on_the_sample():
+    a = normalise([40.0, 60.0], method="goalpost", goalposts=(0, 100))
+    b = normalise([40.0, 60.0, 95.0], method="goalpost", goalposts=(0, 100))
+    assert a.tolist() == pytest.approx(b.iloc[:2].tolist())
+    assert a.tolist() == pytest.approx([0.4, 0.6])
+
+
+def test_minmax_does_depend_on_the_sample_which_is_why_goalposts_exist():
+    a = normalise([40.0, 60.0])
+    b = normalise([40.0, 60.0, 95.0])
+    assert a.tolist() != pytest.approx(b.iloc[:2].tolist())
+
+
+def test_goalpost_without_goalposts_refuses():
+    with pytest.raises(ValueError, match="needs explicit goalposts"):
+        normalise([1.0, 2.0], method="goalpost")
+
+
+def test_a_constant_indicator_becomes_the_midpoint_not_a_division_by_zero():
+    out = normalise([7.0, 7.0, 7.0])
+    assert out.tolist() == [0.5, 0.5, 0.5]
+
+
+def test_missing_stays_missing():
+    out = normalise([1.0, np.nan, 3.0])
+    assert np.isnan(out.iloc[1])
+    assert out.iloc[0] == 0.0 and out.iloc[2] == 1.0
+
+
+# --------------------------------------------------------------------------
+# composite_index
+# --------------------------------------------------------------------------
+
+def test_direction_is_applied_so_a_mortality_rate_does_not_cancel_coverage(districts):
+    """Without the direction argument, imm_pct and imr pull against each other."""
+    right = composite_index(districts, ["imm_pct", "imr"],
+                            directions={"imr": "lower_is_better"})
+    wrong = composite_index(districts, ["imm_pct", "imr"])
+
+    # Expected computed from the definition rather than pasted: min-max each
+    # column over its own observed range, flip the mortality rate, average.
+    imm = districts["imm_pct"]
+    imr = districts["imr"]
+    a = (imm - imm.min()) / (imm.max() - imm.min())
+    b = 1.0 - (imr - imr.min()) / (imr.max() - imr.min())
+    assert right["index"].tolist() == pytest.approx(((a + b) / 2).tolist())
+
+    # Left as higher-is-better the two indicators pull against each other and
+    # very nearly cancel, which is the bug the direction argument exists to stop.
+    assert wrong["index"].std() < 0.1 * right["index"].std()
+
+
+def test_weights_are_rescaled_to_sum_to_one(districts):
+    a = composite_index(districts, ["imm_pct", "piped_pct"],
+                        weights={"imm_pct": 3, "piped_pct": 1})
+    b = composite_index(districts, ["imm_pct", "piped_pct"],
+                        weights={"imm_pct": 0.75, "piped_pct": 0.25})
+    assert a["index"].tolist() == pytest.approx(b["index"].tolist())
+    assert a.attrs["composite_index"]["weights"]["imm_pct"] == pytest.approx(0.75)
+
+
+def test_unnamed_indicators_split_the_residual_weight(districts):
+    r = composite_index(districts, ["imm_pct", "imr", "piped_pct"],
+                        weights={"imm_pct": 0.5},
+                        directions={"imr": "lower_is_better"})
+    w = r.attrs["composite_index"]["weights"]
+    assert w["imm_pct"] == pytest.approx(0.5)
+    assert w["imr"] == pytest.approx(0.25)
+    assert w["piped_pct"] == pytest.approx(0.25)
+
+
+def test_weights_that_leave_no_residual_refuse(districts):
+    with pytest.raises(ValueError, match="leaving\nnothing|leaving nothing"):
+        composite_index(districts, ["imm_pct", "imr", "piped_pct"],
+                        weights={"imm_pct": 0.6, "imr": 0.4},
+                        directions={"imr": "lower_is_better"})
+
+
+def test_geometric_punishes_a_near_zero_where_arithmetic_averages_it_away():
+    df = pd.DataFrame({
+        "u": ["balanced", "lopsided"],
+        "a": [0.5, 1.0],
+        "b": [0.5, 0.0],
+    })
+    ari = composite_index(df, ["a", "b"], method="goalpost",
+                          goalposts={"a": (0, 1), "b": (0, 1)})
+    geo = composite_index(df, ["a", "b"], method="goalpost", aggregation="geometric",
+                          goalposts={"a": (0, 1), "b": (0, 1)})
+    # arithmetic scores them equally; geometric does not
+    assert ari["index"].iloc[0] == pytest.approx(ari["index"].iloc[1])
+    assert geo["index"].iloc[0] > geo["index"].iloc[1]
+
+
+def test_geometric_refuses_z_scores(districts):
+    with pytest.raises(ValueError, match="geometric mean of z-scores"):
+        composite_index(districts, ["imm_pct", "piped_pct"],
+                        method="zscore", aggregation="geometric")
+
+
+def test_a_row_scored_on_two_of_six_indicators_is_refused_by_default():
+    df = pd.DataFrame({
+        "u": ["complete", "sparse"],
+        "a": [0.5, 1.0], "b": [0.5, 1.0], "c": [0.5, np.nan],
+        "d": [0.5, np.nan], "e": [0.5, np.nan], "f": [0.5, np.nan],
+    })
+    cols = list("abcdef")
+    strict = composite_index(df, cols)
+    assert not np.isnan(strict["index"].iloc[0])
+    assert np.isnan(strict["index"].iloc[1])          # would otherwise have topped the list
+
+    lenient = composite_index(df, cols, min_indicators=2)
+    assert not np.isnan(lenient["index"].iloc[1])
+    assert lenient["index_n_indicators"].tolist() == [6, 2]
+
+
+def test_present_indicators_are_reweighted_so_a_gap_does_not_score_as_a_zero():
+    df = pd.DataFrame({"a": [1.0, 1.0], "b": [1.0, np.nan], "c": [1.0, 1.0]})
+    r = composite_index(df, ["a", "b", "c"], method="goalpost", min_indicators=2,
+                        goalposts={c: (0, 1) for c in "abc"})
+    assert r["index"].tolist() == pytest.approx([1.0, 1.0])
+
+
+def test_choices_are_recorded_on_the_result(districts):
+    r = composite_index(districts, ["imm_pct", "imr"],
+                        directions={"imr": "lower_is_better"}, method="rank")
+    meta = r.attrs["composite_index"]
+    assert meta["method"] == "rank"
+    assert meta["directions"] == {"imr": "lower_is_better"}
+    assert meta["indicators"] == ["imm_pct", "imr"]
+
+
+def test_a_typo_in_directions_raises_rather_than_being_ignored(districts):
+    with pytest.raises(KeyError, match="directions names column"):
+        composite_index(districts, ["imm_pct", "imr"], directions={"IMR": "lower_is_better"})
+
+
+def test_one_indicator_is_not_an_index(districts):
+    with pytest.raises(ValueError, match="not an index"):
+        composite_index(districts, ["imm_pct"])
+
+
+# --------------------------------------------------------------------------
+# rank_sensitivity
+# --------------------------------------------------------------------------
+
+def test_sensitivity_reports_zero_spread_when_indicators_agree(districts):
+    out = rank_sensitivity(districts, ["imm_pct", "imr", "piped_pct"], unit="district",
+                           directions={"imr": "lower_is_better"})
+    assert out["rank_spread"].max() == 0
+
+
+def test_sensitivity_finds_a_district_whose_rank_depends_on_the_method():
+    """One district with a huge outlier on one indicator and a poor showing on
+    the other. z-score rewards the outlier's magnitude; rank throws it away."""
+    df = pd.DataFrame({
+        "district": ["outlier", "steady1", "steady2", "steady3", "steady4"],
+        "a": [1000.0, 60.0, 55.0, 50.0, 45.0],
+        "b": [1.0, 60.0, 58.0, 56.0, 54.0],
+    })
+    out = rank_sensitivity(df, ["a", "b"], unit="district")
+    assert out["rank_spread"].max() > 0
+    worst = out.iloc[0]
+    assert worst["district"] == "outlier"
diff --git a/tests/test_inequality.py b/tests/test_inequality.py
new file mode 100644
index 0000000..92739d9
--- /dev/null
+++ b/tests/test_inequality.py
@@ -0,0 +1,466 @@
+"""Tests for the inequality package.
+
+Pinned to answers known in closed form, not to what this code happens to
+return. A distributional index cannot be checked by eye: a Gini of 0.31 and a
+Gini of 0.34 both look like plausible consumption inequality, and a sign error
+in a concentration index produces a number in exactly the right range that says
+the opposite of the truth. So every test here is either an identity the measure
+must satisfy, or a figure worked out by hand from the definition.
+
+Where a number is hardcoded, the arithmetic behind it is in the comment.
+"""
+
+import math
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from inequality import (
+    achievement_index, atkinson, benefit_incidence, benefit_incidence_by_level,
+    concentration_curve, concentration_index, concentration_index_by,
+    erreygers_index, fractional_rank, generalised_entropy, gini, lorenz_curve,
+    oaxaca_blinder, palma_ratio, quantile_shares, ratio_80_20, share_of_bottom,
+    share_of_top, theil_decomposition, theil_l, theil_t, wagstaff_index,
+    weighted_quantile,
+)
+
+# --------------------------------------------------------------------------
+# Fractional ranks
+# --------------------------------------------------------------------------
+
+def test_mean_fractional_rank_is_exactly_one_half():
+    """The proof of CI = 2·cov(h,r)/h̄ needs r̄ = 0.5. It holds for any weights."""
+    rng = np.random.default_rng(0)
+    for _ in range(20):
+        n = rng.integers(2, 60)
+        x = rng.normal(size=n)
+        w = rng.gamma(2.0, 1.0, size=n)
+        r = fractional_rank(x, w)
+        assert np.isclose(np.sum(w * r) / np.sum(w), 0.5)
+
+
+def test_ties_share_the_block_midpoint():
+    # values 1,1,2 with unit weights: the pair occupies [0, 2/3], midpoint 1/3;
+    # the single occupies [2/3, 1], midpoint 5/6.
+    assert np.allclose(fractional_rank([1, 1, 2]), [1 / 3, 1 / 3, 5 / 6])
+
+
+def test_rank_does_not_depend_on_row_order():
+    x = [3, 1, 2, 1, 3]
+    r1 = fractional_rank(x)
+    order = [4, 0, 3, 2, 1]
+    r2 = fractional_rank([x[i] for i in order])
+    assert np.allclose(sorted(r1), sorted(r2))
+
+
+def test_missing_values_rank_as_nan_not_as_poorest():
+    r = fractional_rank([1.0, np.nan, 3.0])
+    assert math.isnan(r[1])
+    assert np.allclose(r[[0, 2]], [0.25, 0.75])
+
+
+# --------------------------------------------------------------------------
+# Gini
+# --------------------------------------------------------------------------
+
+def test_gini_matches_hand_computed_mean_difference():
+    # y = 1..5. Sum of |yi - yj| over ordered pairs = 40, n^2 = 25, so the mean
+    # absolute difference is 1.6 and G = 1.6 / (2 * 3) = 0.2666...
+    assert gini([1, 2, 3, 4, 5]) == pytest.approx(4 / 15)
+
+
+def test_gini_equals_mean_difference_formula_on_random_data():
+    rng = np.random.default_rng(3)
+    for _ in range(15):
+        y = rng.gamma(2.0, 3.0, size=rng.integers(3, 40))
+        mad = np.abs(y[:, None] - y[None, :]).mean()
+        assert gini(y) == pytest.approx(mad / (2 * y.mean()))
+
+
+def test_gini_is_zero_at_equality_and_one_in_the_limit():
+    assert gini([7, 7, 7, 7]) == pytest.approx(0.0, abs=1e-12)
+    # one person holds everything, n large: G -> 1 - 1/n
+    n = 1000
+    y = np.zeros(n); y[-1] = 1.0
+    assert gini(y) == pytest.approx(1 - 1 / n, abs=1e-9)
+
+
+def test_gini_is_scale_invariant():
+    y = [3.0, 1.0, 4.0, 1.0, 5.0, 9.0]
+    assert gini(y) == pytest.approx(gini(np.array(y) * 137.0))
+
+
+def test_a_weight_of_two_equals_the_row_appearing_twice():
+    assert gini([1, 2], weights=[2, 1]) == pytest.approx(gini([1, 1, 2]))
+    assert theil_t([1, 2], weights=[2, 1]) == pytest.approx(theil_t([1, 1, 2]))
+    assert atkinson([1, 2], weights=[2, 1]) == pytest.approx(atkinson([1, 1, 2]))
+
+
+def test_gini_rejects_negative_values_rather_than_returning_a_number():
+    with pytest.raises(ValueError, match="negative"):
+        gini([-1.0, 2.0, 3.0])
+
+
+def test_small_sample_correction_is_the_n_over_n_minus_one_factor():
+    y = [1, 2, 3, 4, 5]
+    assert gini(y, small_sample_correction=True) == pytest.approx(gini(y) * 5 / 4)
+
+
+# --------------------------------------------------------------------------
+# Lorenz curve, shares
+# --------------------------------------------------------------------------
+
+def test_lorenz_area_reproduces_the_gini():
+    rng = np.random.default_rng(11)
+    y = rng.gamma(2.0, 5.0, size=200)
+    w = rng.gamma(3.0, 1.0, size=200)
+    curve = lorenz_curve(y, w)
+    area = np.trapezoid(curve["value_share"], curve["population_share"])
+    assert 2 * (0.5 - area) == pytest.approx(gini(y, w), abs=1e-9)
+
+
+def test_lorenz_starts_at_origin_and_ends_at_one():
+    c = lorenz_curve([4, 1, 9, 2])
+    assert c.iloc[0].tolist() == [0.0, 0.0]
+    assert c.iloc[-1]["population_share"] == pytest.approx(1.0)
+    assert c.iloc[-1]["value_share"] == pytest.approx(1.0)
+
+
+def test_quantile_shares_sum_to_one_with_awkward_weights():
+    rng = np.random.default_rng(5)
+    y = rng.gamma(2.0, 1.0, size=37)
+    w = rng.gamma(1.0, 4.0, size=37)          # nothing divides evenly by five
+    q = quantile_shares(y, w, q=5)
+    assert q["value_share"].sum() == pytest.approx(1.0)
+    assert len(q) == 5
+
+
+def test_shares_at_perfect_equality_equal_the_population_share():
+    y = [12.0] * 50
+    assert share_of_top(y, p=0.10) == pytest.approx(0.10)
+    assert share_of_bottom(y, p=0.40) == pytest.approx(0.40)
+    assert palma_ratio(y) == pytest.approx(0.25)
+    assert ratio_80_20(y) == pytest.approx(1.0)
+
+
+def test_palma_rises_when_the_top_pulls_away():
+    base = list(range(1, 101))
+    pulled = base[:-1] + [500]
+    assert palma_ratio(pulled) > palma_ratio(base)
+
+
+def test_weighted_quantile_is_consistent_with_the_rank_definition():
+    y = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
+    # midpoint ranks are 0.1, 0.3, 0.5, 0.7, 0.9
+    assert weighted_quantile(y, 0.5) == pytest.approx(3.0)
+    assert weighted_quantile(y, 0.3) == pytest.approx(2.0)
+
+
+# --------------------------------------------------------------------------
+# Generalised entropy, Theil, Atkinson
+# --------------------------------------------------------------------------
+
+def test_ge2_is_half_the_squared_coefficient_of_variation():
+    y = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
+    cv2 = (y.std() / y.mean()) ** 2
+    assert generalised_entropy(y, alpha=2.0) == pytest.approx(0.5 * cv2)
+
+
+def test_every_index_is_zero_at_perfect_equality():
+    y = [6.0] * 12
+    assert theil_t(y) == pytest.approx(0.0, abs=1e-12)
+    assert theil_l(y) == pytest.approx(0.0, abs=1e-12)
+    assert atkinson(y, epsilon=2.0) == pytest.approx(0.0, abs=1e-12)
+    assert generalised_entropy(y, alpha=2.0) == pytest.approx(0.0, abs=1e-12)
+
+
+def test_atkinson_is_zero_at_no_aversion_and_rises_with_it():
+    y = [1.0, 2.0, 3.0, 10.0]
+    assert atkinson(y, epsilon=0.0) == 0.0
+    values = [atkinson(y, epsilon=e) for e in (0.0, 0.5, 1.0, 1.5, 2.0)]
+    assert all(b >= a for a, b in zip(values, values[1:]))
+
+
+def test_atkinson_at_epsilon_one_is_one_minus_the_geometric_over_arithmetic_mean():
+    y = np.array([1.0, 4.0, 16.0])
+    expected = 1 - np.exp(np.log(y).mean()) / y.mean()
+    assert atkinson(y, epsilon=1.0) == pytest.approx(expected)
+
+
+def test_indices_are_scale_invariant():
+    rng = np.random.default_rng(13)
+    y = rng.gamma(2.0, 2.0, size=50)
+    for fn in (theil_t, theil_l, lambda v: atkinson(v, epsilon=1.5)):
+        assert fn(y) == pytest.approx(fn(y * 1000.0))
+
+
+def test_ge_below_one_refuses_zeros_instead_of_returning_inf():
+    with pytest.raises(ValueError, match="undefined"):
+        generalised_entropy([0.0, 1.0, 2.0], alpha=0.0)
+    with pytest.raises(ValueError, match="epsilon"):
+        atkinson([0.0, 1.0, 2.0], epsilon=1.0)
+
+
+# --------------------------------------------------------------------------
+# Concentration index
+# --------------------------------------------------------------------------
+
+def test_concentration_index_hand_computed():
+    # outcome on the poorest of four equal-weight people. Ranks 0.125, 0.375,
+    # 0.625, 0.875; mean outcome 0.25.
+    # CI = 2 * (1*0.125) / (4 * 0.25) - 1 = 0.25 - 1 = -0.75
+    assert concentration_index([1, 0, 0, 0], rank_by=[1, 2, 3, 4]) == pytest.approx(-0.75)
+    assert concentration_index([0, 0, 0, 1], rank_by=[1, 2, 3, 4]) == pytest.approx(0.75)
+
+
+def test_concentration_index_is_zero_when_the_outcome_is_flat():
+    assert concentration_index([1, 1, 1, 1], rank_by=[4, 3, 2, 1]) == pytest.approx(0.0, abs=1e-12)
+
+
+def test_concentration_index_depends_only_on_the_rank_order_not_the_values():
+    out = [0.2, 0.5, 0.1, 0.9]
+    a = concentration_index(out, rank_by=[1, 2, 3, 4])
+    b = concentration_index(out, rank_by=[10, 250, 3000, 99999])
+    assert a == pytest.approx(b)
+
+
+def test_concentration_index_is_scale_invariant_in_the_outcome():
+    out = np.array([0.2, 0.5, 0.1, 0.9])
+    r = [1, 2, 3, 4]
+    assert concentration_index(out, r) == pytest.approx(concentration_index(out * 77.0, r))
+
+
+def test_reversing_the_rank_flips_the_sign():
+    out = [0.9, 0.4, 0.3, 0.05]
+    assert (concentration_index(out, rank_by=[1, 2, 3, 4])
+            == pytest.approx(-concentration_index(out, rank_by=[4, 3, 2, 1])))
+
+
+def test_concentration_curve_area_reproduces_the_index():
+    rng = np.random.default_rng(21)
+    n = 300
+    rank = rng.uniform(size=n)
+    outcome = rng.binomial(1, np.clip(0.7 - 0.5 * rank, 0, 1)).astype(float)
+    w = rng.gamma(2.0, 1.0, size=n)
+    curve = concentration_curve(outcome, rank, w)
+    area = np.trapezoid(curve["outcome_share"], curve["population_share"])
+    assert 2 * (0.5 - area) == pytest.approx(concentration_index(outcome, rank, w), abs=1e-9)
+
+
+def test_erreygers_is_four_mu_ci_for_a_binary_outcome():
+    out = [1, 1, 0, 0, 0, 0, 1, 0]
+    rank = [1, 2, 3, 4, 5, 6, 7, 8]
+    mu = np.mean(out)
+    assert erreygers_index(out, rank) == pytest.approx(4 * mu * concentration_index(out, rank))
+
+
+def test_erreygers_refuses_bounds_the_data_breaks():
+    with pytest.raises(ValueError, match="outside the stated"):
+        erreygers_index([0, 1, 5], rank_by=[1, 2, 3], bounds=(0, 1))
+
+
+def test_wagstaff_normalisation_divides_by_one_minus_the_mean():
+    out = [1, 0, 0, 0]
+    rank = [1, 2, 3, 4]
+    assert wagstaff_index(out, rank) == pytest.approx(concentration_index(out, rank) / (1 - 0.25))
+
+
+def test_wagstaff_refuses_when_everybody_has_the_outcome():
+    with pytest.raises(ValueError, match="upper bound"):
+        wagstaff_index([1, 1, 1, 1], rank_by=[1, 2, 3, 4])
+
+
+def test_ranking_by_a_quintile_gives_the_same_sign_as_ranking_by_the_continuous_measure():
+    rng = np.random.default_rng(31)
+    n = 2000
+    wealth = rng.normal(size=n)
+    quintile = pd.qcut(wealth, 5, labels=False) + 1
+    stunted = rng.binomial(1, 1 / (1 + np.exp(1.2 * wealth))).astype(float)
+    fine = concentration_index(stunted, wealth)
+    coarse = concentration_index(stunted, quintile)
+    assert fine < 0 and coarse < 0
+    # the quintile throws away within-quintile ordering, so it attenuates
+    assert abs(coarse) < abs(fine)
+
+
+def test_concentration_index_by_reports_rather_than_drops_a_degenerate_group():
+    df = pd.DataFrame({
+        "stunted": [1, 0, 1, 0, 0, 0],
+        "wealth": [1, 2, 3, 4, 2, 2],
+        "state": ["A", "A", "A", "A", "B", "B"],
+    })
+    out = concentration_index_by(df, "stunted", "wealth", by="state")
+    assert set(out["state"]) == {"A", "B"}
+    b = out.loc[out["state"] == "B"].iloc[0]
+    assert math.isnan(b["concentration_index"])
+    assert b["note"]
+
+
+def test_achievement_index_discounts_a_pro_rich_distribution():
+    pro_poor = [1, 1, 0, 0]
+    pro_rich = [0, 0, 1, 1]
+    rank = [1, 2, 3, 4]
+    assert np.mean(pro_poor) == np.mean(pro_rich)          # same headline coverage
+    assert achievement_index(pro_poor, rank) > achievement_index(pro_rich, rank)
+
+
+# --------------------------------------------------------------------------
+# Decompositions
+# --------------------------------------------------------------------------
+
+@pytest.mark.parametrize("alpha", [0.0, 1.0])
+def test_theil_within_plus_between_equals_total(alpha):
+    rng = np.random.default_rng(41)
+    n = 300
+    y = rng.gamma(2.0, 3.0, size=n)
+    g = rng.integers(0, 5, size=n)
+    w = rng.gamma(2.0, 1.0, size=n)
+    d = theil_decomposition(y, g, w, alpha=alpha)
+    assert d["within"] + d["between"] == pytest.approx(d["total"])
+
+
+def test_between_is_zero_when_groups_have_the_same_distribution():
+    d = theil_decomposition([1, 2, 3, 1, 2, 3], ["a"] * 3 + ["b"] * 3)
+    assert d["between"] == pytest.approx(0.0, abs=1e-12)
+
+
+def test_within_is_zero_when_each_group_is_internally_equal():
+    d = theil_decomposition([1, 1, 1, 5, 5, 5], ["a"] * 3 + ["b"] * 3)
+    assert d["within"] == pytest.approx(0.0, abs=1e-12)
+    assert d["between_share"] == pytest.approx(1.0)
+
+
+def test_theil_decomposition_refuses_an_alpha_that_does_not_decompose():
+    with pytest.raises(ValueError, match="does not decompose"):
+        theil_decomposition([1, 2, 3], ["a", "a", "b"], alpha=2.0)
+
+
+def _oaxaca_frame(seed, beta_a, beta_b, mean_a=8.0, mean_b=5.0, n=400):
+    rng = np.random.default_rng(seed)
+    edu = np.concatenate([rng.normal(mean_a, 2, n), rng.normal(mean_b, 2, n)])
+    grp = np.array(["other"] * n + ["st"] * n)
+    coef = np.where(grp == "other", beta_a, beta_b)
+    y = 1.0 + coef * edu + rng.normal(0, 0.05, 2 * n)
+    return pd.DataFrame({"y": y, "edu": edu, "g": grp})
+
+
+def test_oaxaca_components_sum_to_the_gap():
+    df = _oaxaca_frame(51, 0.5, 0.3)
+    for ref in ("pooled", "advantaged", "disadvantaged"):
+        r = oaxaca_blinder(df, "y", "g", ["edu"], reference=ref)
+        assert r["explained"] + r["unexplained"] == pytest.approx(r["gap"])
+    r3 = oaxaca_blinder(df, "y", "g", ["edu"], reference="threefold")
+    assert (r3["endowments"] + r3["coefficients"] + r3["interaction"]
+            == pytest.approx(r3["gap"]))
+
+
+def test_oaxaca_unexplained_vanishes_when_the_returns_are_identical():
+    df = _oaxaca_frame(52, 0.5, 0.5)
+    r = oaxaca_blinder(df, "y", "g", ["edu"], reference="pooled")
+    assert abs(r["unexplained"]) < 0.02 * abs(r["gap"])
+
+
+def test_oaxaca_explained_vanishes_when_the_endowments_are_identical():
+    """Identical *realised* endowments, not merely the same drawing distribution.
+
+    Drawing both groups from N(6, 2) leaves a sampling difference in the means
+    of order 0.14 at n = 400, which the explained component correctly picks up:
+    it is a real difference in this sample's endowments. So the two groups get
+    the same education vector, and then the explained part must be zero to
+    machine precision rather than merely small.
+    """
+    rng = np.random.default_rng(53)
+    edu = rng.normal(6.0, 2.0, 400)
+    df = pd.DataFrame({
+        "edu": np.concatenate([edu, edu]),
+        "g": np.array(["other"] * 400 + ["st"] * 400),
+        "y": np.concatenate([1.0 + 0.5 * edu, 1.0 + 0.3 * edu]),
+    })
+    r = oaxaca_blinder(df, "y", "g", ["edu"], reference="pooled")
+    assert r["explained"] == pytest.approx(0.0, abs=1e-9)
+    assert r["unexplained"] == pytest.approx(r["gap"])
+
+
+def test_oaxaca_refuses_more_than_two_groups():
+    df = _oaxaca_frame(54, 0.5, 0.3)
+    df.loc[df.index[:5], "g"] = "third"
+    with pytest.raises(ValueError, match="exactly two groups"):
+        oaxaca_blinder(df, "y", "g", ["edu"])
+
+
+def test_oaxaca_refuses_an_underidentified_group():
+    df = _oaxaca_frame(55, 0.5, 0.3, n=3)
+    df = pd.concat([df[df.g == "other"], df[df.g == "st"].head(1)])
+    with pytest.raises(ValueError, match="not identified"):
+        oaxaca_blinder(df, "y", "g", ["edu"])
+
+
+# --------------------------------------------------------------------------
+# Benefit incidence
+# --------------------------------------------------------------------------
+
+def test_benefit_incidence_allocates_the_whole_budget_by_utilisation_share():
+    d = pd.DataFrame({"visits": [3, 2, 1, 0], "q": [1, 2, 3, 4]})
+    out = benefit_incidence(d, utilisation="visits", group="q", spending=600)
+    assert out["benefit"].sum() == pytest.approx(600)
+    assert out["benefit"].tolist() == pytest.approx([300, 200, 100, 0])
+    assert out.attrs["concentration_index"] < 0
+    assert "progressive" in out.attrs["verdict"]
+
+
+def test_benefit_incidence_refuses_blank_utilisation():
+    d = pd.DataFrame({"visits": [3, None, 1], "q": [1, 2, 3]})
+    with pytest.raises(ValueError, match="missing utilisation"):
+        benefit_incidence(d, utilisation="visits", group="q", spending=100)
+
+
+def test_net_benefit_subtracts_what_the_household_paid():
+    d = pd.DataFrame({"visits": [1, 1], "q": [1, 2], "fee": [10.0, 90.0]})
+    out = benefit_incidence(d, utilisation="visits", group="q", spending=100, fees="fee")
+    assert out.sort_values("q")["net_benefit"].tolist() == pytest.approx([40.0, -40.0])
+
+
+def test_level_split_reveals_what_the_aggregate_hides():
+    d = pd.DataFrame({
+        "u": [5, 4, 3, 2, 1, 0, 1, 2, 4, 8],
+        "q": [1, 2, 3, 4, 5] * 2,
+        "lvl": ["primary"] * 5 + ["tertiary"] * 5,
+    })
+    r = benefit_incidence_by_level(d, utilisation="u", group="q", level="lvl",
+                                   spending={"primary": 1000, "tertiary": 1000})
+    assert r["concentration_index_by_level"]["primary"] < 0
+    assert r["concentration_index_by_level"]["tertiary"] > 0
+    assert r["total"]["benefit"].sum() == pytest.approx(2000)
+
+
+def test_level_split_refuses_a_level_with_no_budget():
+    d = pd.DataFrame({"u": [1, 1], "q": [1, 2], "lvl": ["primary", "tertiary"]})
+    with pytest.raises(ValueError, match="no budget given"):
+        benefit_incidence_by_level(d, utilisation="u", group="q", level="lvl",
+                                   spending={"primary": 100})
+
+
+# --------------------------------------------------------------------------
+# Input handling shared across the package
+# --------------------------------------------------------------------------
+
+def test_missing_rows_are_dropped_pairwise_not_treated_as_zero():
+    with_gap = gini([1.0, 2.0, np.nan, 4.0])
+    without = gini([1.0, 2.0, 4.0])
+    assert with_gap == pytest.approx(without)
+
+
+def test_negative_weights_are_refused_everywhere():
+    for call in (lambda: gini([1, 2, 3], weights=[1, -1, 1]),
+                 lambda: concentration_index([1, 0, 1], [1, 2, 3], weights=[1, -1, 1]),
+                 lambda: theil_t([1, 2, 3], weights=[1, -1, 1])):
+        with pytest.raises(ValueError, match="negative weights"):
+            call()
+
+
+def test_mismatched_lengths_are_refused():
+    with pytest.raises(ValueError, match="rows"):
+        gini([1, 2, 3], weights=[1, 1])
+    with pytest.raises(ValueError, match="rows"):
+        concentration_index([1, 0, 1], rank_by=[1, 2])
diff --git a/workflows/ci.yml b/workflows/ci.yml
deleted file mode 100644
index 40a2838..0000000
--- a/workflows/ci.yml
+++ /dev/null
@@ -1,32 +0,0 @@
-name: CI
-
-on: [push, pull_request]
-
-jobs:
-  python-tests:
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v3
-      - uses: actions/setup-python@v4
-        with:
-          python-version: '3.10'
-      - run: pip install -r requirements.txt
-      - run: pytest tests/
-  r-tests:
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v3
-      - uses: r-lib/actions/setup-r@v2
-      - run: Rscript -e "install.packages(c('testthat','quarto'))"
-      - run: Rscript -e "testthat::test_dir('tests/R')"
-  docs:
-    needs: [python-tests, r-tests]
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v3
-      - uses: r-lib/actions/setup-r@v2
-      - run: Rscript -e "install.packages('quarto')"
-      - run: quarto render docs/index.qmd --to html
-      - uses: peaceiris/actions-gh-pages@v3
-        with:
-          publish_dir: ./docs/_site