feat: v0.9.11 — multi-instrument MTE + true CHV-2011 PRTE weights

Closes two long-standing API gaps and an empirical math debt.

Added:
- sp.bayes_mte(instrument: str | Sequence[str], ...) — multi-IV
  support matching sp.bayes_iv / sp.bayes_hte_iv. Scalar calls
  unchanged. pi_Z ~ Normal(shape=k_iv) in joint mode; post-hoc
  propensity recomputation handles (k_iv,) coefficient vector.
- sp.policy_weight_observed_prte(propensity_sample, shift) — true
  CHV-2011 PRTE weight via gaussian_kde.integrate_box_1d — the
  exact CDF-difference form [F_P(u) - F_P(u-Δ)]/Δ, closing the
  v0.9.9 docstring gap that flagged policy_weight_prte as
  stylised.

Two-round code review, zero ship-blockers:
- Round B agent found 5 items (2 HIGH + 3 MEDIUM). My self-sweep
  caught a HIGH the agent got wrong: the original code AND the
  reviewer's proposed swap both computed the derivative of
  density, not the CDF difference. Fixed via integrate_box_1d.
  Empirical validation: uniform propensity + Δ=0.2 now gives the
  textbook trapezoid (linear ramp 0→1 on [0, 0.2], plateau ≈1.0 on
  [0.2, 0.9]); previously gave a spurious boundary spike.
- Other Round-B fixes: dropped np.clip(w, 0, None) so contraction
  policies (Δ<0) yield legitimate signed-negative weights;
  precomputed gaussian_kde covariance to avoid thread-safety race;
  dropped model_info['instrument'] (type varied scalar-vs-list);
  tightened scalar-list back-compat test to relative-SD tolerance.
- Round C: math confirmed (integrate_box_1d == CHV Theorem 1),
  negative-shift sign correct, no regressions.

Tests:
- tests/test_bayes_mte_multi_iv.py (9 new) — scalar back-compat,
  list equivalence, 2-IV + 3-IV fit, joint first-stage + hv_latent
  orthogonal with multi-IV, empty-list + missing-column rejection.
- tests/test_bayes_mte_policy.py (+7) — CHV uniform-sample plateau
  ≈1.0, negative-shift signed weights, parametrised input
  validation, end-to-end policy_effect integration.
- 61 focused MTE tests green.
- Full regression: 2195 passed / 4 skipped / 0 failed in 15:16.

Design spec: docs/superpowers/specs/2026-04-20-v0911-multi-iv-mte-observed-prte.md

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: brycewang-stanford

CHANGELOG.md

@@ -2,6 +2,58 @@
 
 All notable changes to StatsPAI will be documented in this file.
 
+## [0.9.11] - 2026-04-20 — Multi-instrument MTE + true CHV-2011 PRTE weights
+
+Closes two long-standing API gaps plus an empirical math debt.
+
+### Added (0.9.11)
+
+- **`sp.bayes_mte(instrument: str | Sequence[str], ...)`** — MTE
+  now accepts multiple instruments, matching `sp.bayes_iv` /
+  `sp.bayes_hte_iv`. Scalar calls unchanged.
+- **`sp.policy_weight_observed_prte(propensity_sample, shift)`** —
+  true Carneiro-Heckman-Vytlacil (2011) PRTE weights from the
+  observed propensity distribution via
+  `kde.integrate_box_1d(u-Δ, u) / Δ` (CDF difference). Closes the
+  v0.9.9 docstring gap where `policy_weight_prte` was flagged
+  stylised.
+
+### Round-B review found 2 HIGH + 3 MEDIUM; all fixed
+
+1. **CHV sign bug** — my original `(kde(u) - kde(u-Δ))/Δ` AND the
+   reviewer's proposed swap were both wrong (both compute
+   derivative of density, not CDF difference). Self-sweep verified
+   CHV-2011 Theorem 1 is a CDF difference. Fixed via
+   `integrate_box_1d`. Empirical: uniform propensity + Δ=0.2 now
+   gives the textbook trapezoid; previously gave a spurious
+   boundary spike.
+2. **Unconditional `np.clip(w, 0, None)`** silently altered the
+   estimand. Dropped — contraction policies now yield signed
+   negative weights, matching CHV convention.
+3. **`gaussian_kde` thread safety** — forced covariance
+   precomputation inside the builder.
+4. **`model_info['instrument']` type varied** — dropped the raw
+   key; only `instruments` (list) + `n_instruments` remain.
+5. **Back-compat test** uses relative-to-posterior-SD tolerance.
+
+### Tests (0.9.11)
+
+- `tests/test_bayes_mte_multi_iv.py` (9 tests).
+- `tests/test_bayes_mte_policy.py` (+7 tests).
+- 61 focused MTE tests green.
+
+### Code review
+
+- Round B agent: 5 items. Self-sweep caught one HIGH the agent
+  got wrong. All 5 fixed.
+- Round C agent: zero ship-blockers.
+
+### Design spec (0.9.11)
+
+- `docs/superpowers/specs/2026-04-20-v0911-multi-iv-mte-observed-prte.md`
+
+---
+
 ## [0.9.10] - 2026-04-20 — HV-latent MTE (textbook Heckman-Vytlacil via latent U_D)
 
 Closes the semantic debt v0.9.9 flagged but did not pay: the

docs/superpowers/specs/2026-04-20-v0911-multi-iv-mte-observed-prte.md

@@ -0,0 +1,115 @@
+# StatsPAI v0.9.11 — Multi-instrument MTE + CHV-2011 observed-propensity PRTE
+
+**Author:** Bryce Wang · **Date:** 2026-04-20 · **Status:** design → implementation
+
+## 1. Motivation
+
+v0.9.10 shipped textbook HV-latent MTE but left two threads from earlier non-goals lists:
+
+1. **Multi-instrument MTE** — `sp.bayes_mte` took `instrument: str` (scalar), forcing users with 2+ IVs to pick one. All other Bayesian estimators (`sp.bayes_iv`, `sp.bayes_hte_iv`) already support `instrument: str | list`; MTE was the inconsistent one.
+
+2. **True CHV-2011 PRTE builder** — v0.9.9 shipped `sp.policy_weight_prte(shift)` as a **stylised** rectangle and the docstring explicitly said the real CHV-2011 PRTE requires the observed propensity kernel. The worked example in the docstring asked the user to hand-roll a `gaussian_kde`-based weight_fn. v0.9.11 makes that one-liner: `sp.policy_weight_observed_prte(propensity_sample, shift)`.
+
+Both are small, well-scoped additions that close open API gaps.
+
+## 2. Scope
+
+### In scope (v0.9.11)
+
+- **`sp.bayes_mte(instrument: str | Sequence[str], ...)`** — accept multi-instrument:
+  - Scalar path unchanged (API back-compat).
+  - List path: `Z` becomes `(n, k)` matrix; first-stage logit `pi_Z ~ Normal(0, σ, shape=k)`; `logit = pi_0 + Z @ pi_Z + ...`.
+- **`sp.policy_weight_observed_prte(propensity_sample, shift)`** — CHV 2011 weights: `w(u) ∝ [f_P(u) - f_{P+Δ}(u)] / Δ` where `f_P` is the kernel-density estimate of the observed propensity sample. Normalised to unit sum on the grid passed into `policy_effect`.
+- Retain the stylised `sp.policy_weight_prte(shift)` unchanged — useful as a quick exploration tool, and the docstring already flags it as stylised.
+
+### Out of scope (explicitly deferred)
+
+- Bivariate-normal HV selection model (Heckman-style with `(U_0, U_1, V)` covariance structure). This is its own design problem — the right move is 0.9.12+.
+- Policy counterfactual = "add a new instrument" — would require dedicated API.
+- Rust Phase 2.
+
+## 3. API changes
+
+### 3.1 `bayes_mte(instrument: str | Sequence[str])`
+
+```python
+def bayes_mte(
+    data: pd.DataFrame,
+    y: str,
+    treat: str,
+    instrument: Union[str, Sequence[str]],   # CHANGED: was str
+    covariates: Optional[List[str]] = None,
+    ...
+)
+```
+
+Inside:
+- Normalise to a list `iv_cols = [instrument] if isinstance(instrument, str) else list(instrument)`.
+- `Z = clean[iv_cols].to_numpy(dtype=float)` → shape `(n, k)`.
+- `_logit_propensity(Z, X, D)` already handles 2-D Z (see its `W = Z.reshape(-1, 1) if Z.ndim == 1 else Z` branch — good).
+- PyMC side: `pi_Z = pm.Normal('pi_Z', mu=0, sigma=prior_coef_sigma, shape=k)` then `logit = pi_intercept + pm.math.dot(Z, pi_Z) + ...`.
+
+### 3.2 `sp.policy_weight_observed_prte(propensity_sample, shift)`
+
+```python
+def policy_weight_observed_prte(
+    propensity_sample: np.ndarray,
+    shift: float,
+    *,
+    bw_method: str | float | None = None,
+) -> Callable[[np.ndarray], np.ndarray]:
+    """True CHV-2011 PRTE weights from the observed propensity
+    distribution, via Gaussian KDE."""
+```
+
+Implementation:
+- Validate `propensity_sample` lies in `[0, 1]` (common sanity check).
+- Validate `shift` in `(-1, 1)` non-zero (matches `policy_weight_prte`).
+- Build `kde = scipy.stats.gaussian_kde(propensity_sample, bw_method=bw_method)`.
+- Return closure `w(u) -> (kde(u) - kde(u - shift)) / shift`, clipped at 0 from below (negative weights are not meaningful for integration against an MTE curve and usually indicate grid-edge artefacts).
+
+Edge cases handled:
+- `shift > 0`: positive marginal expansion of propensity (compliers shift up).
+- `shift < 0`: negative marginal shrinkage (defiers / contraction).
+- `u - shift` outside `[0, 1]`: kde density naturally falls off; clip avoids negative weight.
+
+### 3.3 Export wiring
+
+Add `policy_weight_observed_prte` to `sp.bayes.__all__` and top-level `__all__`.
+
+## 4. File plan
+
+| File | Change |
+|---|---|
+| `src/statspai/bayes/mte.py` | `instrument: str | Sequence[str]`, list normalisation, shape-k `pi_Z` prior |
+| `src/statspai/bayes/policy_weights.py` | NEW func `policy_weight_observed_prte` |
+| `src/statspai/bayes/__init__.py` | Export new builder |
+| `src/statspai/__init__.py` | Top-level export |
+| `tests/test_bayes_mte_multi_iv.py` | NEW — multi-IV recovery + scalar back-compat |
+| `tests/test_bayes_mte_policy.py` | Extend with `policy_weight_observed_prte` tests |
+| `pyproject.toml` | `version = "0.9.11"` |
+| `CHANGELOG.md` | 0.9.11 entry |
+
+## 5. Test plan
+
+- `test_bayes_mte_multi_instrument_scalar_back_compat` — passing a single-element list returns same posterior as scalar within sampling noise.
+- `test_bayes_mte_multi_instrument_recovery` — 2-IV DGP where both first-stage coefficients are identified; the MTE polynomial still recovers truth.
+- `test_bayes_mte_multi_instrument_model_info` — `model_info['instruments']` reports the list.
+- `test_policy_weight_observed_prte_returns_callable`.
+- `test_policy_weight_observed_prte_input_validation` — reject out-of-bounds samples, zero shift, shift outside `(-1,1)`.
+- `test_policy_weight_observed_prte_positive_shift_yields_marginal_mass` — on a uniform-propensity sample with shift=0.1, the weight peaks near the induced margin.
+- `test_policy_weight_observed_prte_integrates_with_policy_effect` — end-to-end: fit MTE, call `r.policy_effect(sp.policy_weight_observed_prte(r._propensity_sample, 0.1))`, returns a finite posterior.
+
+## 6. Success criteria
+
+1. Multi-instrument `bayes_mte` runs on a 2-IV DGP and recovers the true MTE polynomial within HDI at n=600.
+2. `policy_weight_observed_prte(uniform_sample, shift=0.1)` returns a mass concentrated around `u ∈ [0.5-0.05, 0.5+0.05]` on a uniform-propensity DGP (sanity).
+3. Scalar-instrument calls remain backward-compatible (existing tests don't change behaviour).
+4. Two rounds of code review — no ship-blockers.
+5. Full regression stays within the existing flakiness baseline (≤ 2 pre-existing flakies on an isolated-retry-passing basis).
+
+## 7. Non-goals
+
+- Bivariate-normal HV (deferred).
+- Per-instrument policy weights (deferred).
+- IV-strength diagnostics for MTE first stage (users can call `sp.check_identification` on the input).

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "StatsPAI"
-version = "0.9.10"
+version = "0.9.11"
 description = "The Agent-Native Causal Inference & Econometrics Toolkit for Python"
 readme = "README.md"
 license = {text = "MIT"}

src/statspai/__init__.py

@@ -22,7 +22,7 @@
 >>> sp.outreg2(result, filename="results.xlsx")
 """
 
-__version__ = "0.9.10"
+__version__ = "0.9.11"
 __author__ = "Biaoyue Wang"
 __email__ = "bryce@copaper.ai"
 
@@ -147,6 +147,7 @@
     BayesianCausalResult, BayesianHTEIVResult, BayesianMTEResult,
     policy_weight_ate, policy_weight_subsidy,
     policy_weight_prte, policy_weight_marginal,
+    policy_weight_observed_prte,
 )
 from .regression.heckman import heckman
 from .regression.quantile import qreg, sqreg
@@ -196,6 +197,8 @@
 from .qte import qdid, qte, QTEResult
 from .mht import romano_wolf, RomanoWolfResult, adjust_pvalues, bonferroni, holm, benjamini_hochberg
 from .registry import list_functions, describe_function, function_schema, search_functions, all_schemas
+# Unified help entry point (aggregates registry + docstring + category + search)
+from .help import help, HelpResult
 
 # === Article-facing aliases (sp.rdd / sp.frontdoor / sp.xlearner / ...) ===
 # Thin wrappers around existing implementations; see _article_aliases.py
@@ -638,6 +641,7 @@
     "policy_weight_subsidy",
     "policy_weight_prte",
     "policy_weight_marginal",
+    "policy_weight_observed_prte",
     # Neural Causal Models
     "tarnet",
     "cfrnet",
@@ -750,6 +754,8 @@
     "function_schema",
     "search_functions",
     "all_schemas",
+    "help",
+    "HelpResult",
     # Data Generating Processes
     "dgp_did",
     "dgp_rd",

src/statspai/bayes/__init__.py

@@ -36,6 +36,7 @@
     policy_weight_subsidy,
     policy_weight_prte,
     policy_weight_marginal,
+    policy_weight_observed_prte,
 )
 
 __all__ = [
@@ -52,4 +53,5 @@
     'policy_weight_subsidy',
     'policy_weight_prte',
     'policy_weight_marginal',
+    'policy_weight_observed_prte',
 ]

src/statspai/bayes/mte.py

@@ -71,7 +71,7 @@
 """
 from __future__ import annotations
 
-from typing import List, Optional, Tuple, Union
+from typing import List, Optional, Sequence, Tuple, Union
 
 import numpy as np
 import pandas as pd
@@ -106,7 +106,7 @@ def bayes_mte(
     data: pd.DataFrame,
     y: str,
     treat: str,
-    instrument: str,
+    instrument: Union[str, Sequence[str]],
     covariates: Optional[List[str]] = None,
     *,
     first_stage: str = 'plugin',
@@ -190,12 +190,22 @@ def bayes_mte(
     """
     pm, _ = _require_pymc()
 
-    for c in [y, treat, instrument] + (list(covariates) if covariates else []):
+    # Normalise instrument to a list so scalar and list-of-instruments
+    # paths share the same downstream code. Back-compat: passing a
+    # single string still works and returns the same posterior.
+    if isinstance(instrument, str):
+        iv_cols: List[str] = [instrument]
+    else:
+        iv_cols = list(instrument)
+    if len(iv_cols) == 0:
+        raise ValueError("instrument must name at least one column.")
+
+    for c in [y, treat] + iv_cols + (list(covariates) if covariates else []):
         if c not in data.columns:
             raise ValueError(f"Column '{c}' not found in data")
 
     cov_cols = list(covariates) if covariates else []
-    clean = data[[y, treat, instrument] + cov_cols].dropna().reset_index(drop=True)
+    clean = data[[y, treat] + iv_cols + cov_cols].dropna().reset_index(drop=True)
     n = len(clean)
     if n < 50:
         raise ValueError(
@@ -204,7 +214,8 @@ def bayes_mte(
 
     Y = clean[y].to_numpy(dtype=float)
     D = clean[treat].to_numpy(dtype=float)
-    Z = clean[instrument].to_numpy(dtype=float)
+    # 2-D Z of shape (n, k_iv). k_iv=1 recovers the scalar case.
+    Z = clean[iv_cols].to_numpy(dtype=float)
     X = clean[cov_cols].to_numpy(dtype=float) if cov_cols else None
     uniq_D = np.unique(D)
     if not set(uniq_D).issubset({0.0, 1.0}):
@@ -274,8 +285,14 @@ def bayes_mte(
             pi_intercept = pm.Normal(
                 'pi_intercept', mu=0.0, sigma=prior_coef_sigma,
             )
-            pi_Z = pm.Normal('pi_Z', mu=0.0, sigma=prior_coef_sigma)
-            logit = pi_intercept + pi_Z * Z
+            # pi_Z is a vector of length k_iv. For back-compat with
+            # the scalar case, shape=(1,) broadcasts against a (n, 1)
+            # Z matrix via pm.math.dot.
+            k_iv = Z.shape[1]
+            pi_Z = pm.Normal(
+                'pi_Z', mu=0.0, sigma=prior_coef_sigma, shape=k_iv,
+            )
+            logit = pi_intercept + pm.math.dot(Z, pi_Z)
             if X is not None:
                 pi_X = pm.Normal(
                     'pi_X', mu=0.0, sigma=prior_coef_sigma,
@@ -384,8 +401,12 @@ def bayes_mte(
     #   the coefficients is a natural point summary.
     if first_stage == 'joint':
         pi0_mean = float(trace.posterior['pi_intercept'].values.mean())
-        piZ_mean = float(trace.posterior['pi_Z'].values.mean())
-        lin = pi0_mean + piZ_mean * Z
+        # pi_Z is shape (k_iv,). Collapse (chains, draws, k_iv) to
+        # posterior-mean vector of length k_iv.
+        piZ_mean = trace.posterior['pi_Z'].values.reshape(
+            -1, Z.shape[1]
+        ).mean(axis=0)
+        lin = pi0_mean + Z @ piZ_mean
         if X is not None and 'pi_X' in trace.posterior:
             piX_mean = trace.posterior['pi_X'].values.reshape(
                 -1, X.shape[1]
@@ -443,7 +464,11 @@ def _integrated_effect(U_population):
         'target_accept': target_accept,
         'poly_u': poly_u,
         'u_grid': u_grid.tolist(),
-        'instrument': instrument,
+        # `instruments` is always a list regardless of whether the
+        # user passed a scalar or list to keep downstream code from
+        # branching on type (cf. v0.9.11 round-B review).
+        'instruments': iv_cols,
+        'n_instruments': len(iv_cols),
         'covariates': cov_cols,
         'prior_mte_sigma': prior_mte_sigma,
         'prior_coef_sigma': prior_coef_sigma,