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