From a1cc26e0424cbea4340dfa19f1d9875e6b454c4f Mon Sep 17 00:00:00 2001 From: guzhujushi Date: Sat, 3 Oct 2026 19:59:38 +0800 Subject: [PATCH] =?UTF-8?q?feat(stats):=20bootstrap=5Fci=20=E8=90=BD?= =?UTF-8?q?=E5=9C=B0=EF=BC=88L0=2010/12=EF=BC=8CS=20=E6=A1=A3=EF=BC=89+=20?= =?UTF-8?q?=E5=8D=A1=E7=89=87=EF=BC=9B=E4=BF=AE=E6=8E=89=20rng=20=E7=BB=95?= =?UTF-8?q?=E8=BF=87=20get=5Frng=20=E4=B8=8E=20lambda=20=E5=90=8D=E5=AD=97?= =?UTF-8?q?=E6=AD=BB=E5=88=86=E6=94=AF?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ## 审查结论 算法移植是对的：大样本正态下百分位区间半宽与解析 `2 × 1.96 × se` 差 < 2%， 按块重采样（batch_size=1000）也照 docstring 的 Notes 做了。以下是要改的部分。 ## Fixed - `random_state=None` 时另开 `np.random.default_rng()`,**绕过了 `core.get_rng`** —— 与决策 24 冲突：`set_seed` 对它无效，同一段脚本重跑区间就变。改成无条件 `get_rng(random_state)`（决策 46）。 - `metrics["statistic_name"]` 用 `hasattr(statistic, "__name__")` 判不出 lambda： lambda 的 `__name__` 恒为 `""`，于是 `else "custom"` 是**死代码**， 实测自定义统计量拿到的是 `""`。改成 `getattr(f, "__name__", None)` + 判 `""`（易错点 17）。 - docstring 补齐 `Parameters` / `Returns` / `Raises`（§4.7），并如实写明 `method="bca"` **尚未实现** （此前 docstring 承诺支持、实现却抛 `NotImplementedError`）。 - ruff 两条：I001（导入块拆成三行且用了绝对导入）、F541（`f"暂不支持bca"` 无占位符）。 ## Tests - **摘掉 3 条 `skip`**：等 `bootstrap_ci` 的那 3 条一直挂着，实现写完照样全绿 —— 易错点 9 复发， 已在 `personal-pitfalls.md` 把该条 `复发` 记到 1 并注明本次经过。 - 新增 6 条：`estimate` 等于原样本统计量、`confidence` 单调性、五种容器一致、 自定义 lambda 走回退分支、**不给种子时跟随 `set_seed`**（钉住上面那个 rng 修复）、非法/未实现 method。 - 写测试时踩到一个真性质：极差的重采样值**永远 ≤ 原样本极差**，`high` 会正好等于 `estimate`， 断言不能写 `estimate < high`。这条已写进卡片 §5。 - 全量 **177 条：143 passed / 34 skipped / 0 xfail**；ruff 0，pyright 0；`examples/01` 跑通。 ## Docs - `docs/algorithms/bootstrap_ci.md`：**第 2 张 S 档卡片**（模板九节 + §10 待拍板）。 `docs/algorithms/README.md` 进度表打勾。 - `docs/dev/personal-pitfalls.md`：新增第 17 条（`hasattr(f, "__name__")` 挡不住 lambda，`else` 分支即死代码）， 索引补到 17；第 9 条记一次复发。 - `STATUS.md`：E 节记 `bootstrap_ci`（L0 10/12 → 剩余 `monte_carlo` `save_figure`），决策 46， 「下一步」第 1 条改写、新增第 8 条列 5 个待拍板口径。 ## 待用户拍板（Codex 一律未动，按决策 39） - `x` 里的 NaN 静默传进统计量 → `nan` 区间（`confidence_interval` / `t_test` 都是先丢 NaN 的，口径不齐） - `n_resamples <= 0` 现在抛 `IndexError: index -1 is out of bounds`（应在入口给可读 ValueError） - `confidence` 越界交给 numpy 抛 `Percentiles must be in the range [0, 100]` - `method="bca"`：补实现（≈25 行，S 档）还是从 docstring 删掉 - `bootstrap_statistic` 仍是 `NotImplementedError`，而 `bootstrap_ci` 把那段重采样循环内联了 —— 建议抽出去让 `bootstrap_ci` 调它（易错点 16 的变体） --- CHANGELOG.md | 28 +++++++ STATUS.md | 44 +++++++--- docs/algorithms/README.md | 2 +- docs/algorithms/bootstrap_ci.md | 127 +++++++++++++++++++++++++++++ docs/dev/personal-pitfalls.md | 24 +++++- src/mcm_toolkit/stats/bootstrap.py | 89 ++++++++++++++++++-- tests/test_stats.py | 79 +++++++++++++++++- 7 files changed, 373 insertions(+), 20 deletions(-) create mode 100644 docs/algorithms/bootstrap_ci.md diff --git a/CHANGELOG.md b/CHANGELOG.md index 0d07d63..6eda07b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -301,6 +301,34 @@ 统一收成 `TypeError`（决策 45，**Codex 代选、待用户确认**）—— 与 `t_test` / `confidence_interval` 已经写进 docstring 的口径一致。 +### Added +- `stats.bootstrap_ci`（**L0 第 10 项，S 档**，实现由用户本人写）：百分位法 bootstrap 置信区间， + `statistic` 可以是任意 `f(sample) -> float`（默认 `np.mean`）。返回 `Result`， + `params = {estimate, low, high, confidence, n_resamples}`、`metrics.statistic_name`。 + 重采样**按块处理**（`batch_size=1000`），不落地 `B × n` 整矩阵。 +- `docs/algorithms/bootstrap_ci.md`：**第 2 张 S 档函数卡片**。除模板九节外，记了一条非平凡性质 —— + **统计量有上界时百分位区间会塌到上界**（极差的重采样值永远 ≤ 原样本极差，实测 `high == estimate`）， + 以及 5 个待拍板口径。`docs/algorithms/README.md` 的进度表同步打勾。 + +### Fixed +- `bootstrap_ci` 的 `random_state=None` 分支原来另开 `np.random.default_rng()`，**绕过了 `core.get_rng`** —— + 与决策 24 冲突：`set_seed` 对它无效、区间不可复现。改成无条件 `get_rng(random_state)`（决策 46）。 +- `bootstrap_ci` 的 `metrics["statistic_name"]`：`hasattr(statistic, "__name__")` 判不出 lambda —— + lambda 的 `__name__` 恒为 `""`，于是 `else "custom"` 是**死代码**，自定义统计量拿不到 `"custom"`。 + 改成 `getattr(f, "__name__", None)` 再判 `name in (None, "")`（易错点 17）。 +- `bootstrap_ci` 的 docstring 补齐 `Parameters` / `Returns` / `Raises`（§4.7），并如实写明 + `method="bca"` **尚未实现**（此前 docstring 承诺支持、实现却抛 `NotImplementedError`）。 + +### Changed +- `tests/test_stats.py`：**摘掉 3 条 `skip`**（等 `bootstrap_ci` 的那 3 条），并新增 6 条断言 —— + `estimate` 等于原样本统计量、`confidence` 单调性、五种容器结果一致、自定义 lambda 走回退分支、 + **不给种子时跟随 `set_seed`**、非法/未实现 `method`。bootstrap 段共 9 条。 + 全量 **177 条：143 passed / 34 skipped**；ruff 0，pyright 0。 +- `docs/dev/personal-pitfalls.md`：**第 9 条「实现写完了测试还挂着 skip」复发 +1**（这次是 `bootstrap_ci` —— + 实现写完、3 条测试仍挂 skip、全量照样全绿）；新增第 17 条「`hasattr(f, "__name__")` 挡不住 lambda， + `else` 分支是死代码」；索引补到 17。 +- `STATUS.md`：E 节记 `bootstrap_ci`（L0 10/12），决策 46，「下一步」改写并新增第 8 条列 5 个待拍板口径。 + ## [0.0.1] - 2026-10-01 ### Added diff --git a/STATUS.md b/STATUS.md index 27b1819..12cfab2 100644 --- a/STATUS.md +++ b/STATUS.md @@ -205,13 +205,28 @@ - 两条硬规矩：**默认绝不摊平**（要摊平显式 `ndim=1`）；**不删缺失值**（删缺失是各函数的语义， 见决策 45 与易错点 15）。`numeric_only` 按"输出是否与输入同形"选，不按喜好选（决策 44）。 - `check_array`（L1，未实现）实现时直接复用 `as_array`，一份实现两张脸 -- [ ] 剩余 3 项：`bootstrap_ci` `monte_carlo` `save_figure` +- [x] `stats/bootstrap.py` 的 `bootstrap_ci`（L0 第 10 项，**S 档**，实现由用户本人写，Codex 审查 + 代修） + - 移植正确性已用**解析解**钉住：大样本正态下百分位区间半宽 ≈ `2 × 1.96 × se`（容差 10%） + - **Codex 代修的三处**： + ① `random_state=None` 时原来另开了 `np.random.default_rng()`，**绕过了 `get_rng`** —— + 与决策 24 冲突（`set_seed` 后仍然不可复现）。改成无条件的 `get_rng(random_state)`； + ② `statistic_name` 写成 `hasattr(statistic, "__name__")`，而 lambda 的 `__name__` 恒为 `""`， + 导致 `else` 分支是**死代码**、自定义统计量拿不到 `"custom"`（见易错点 17）； + ③ ruff 两条（I001 导入块、F541 空 f-string），并把 docstring 的 `Parameters` / `Returns` / `Raises` 补齐（§4.7） + - **测试是重点**：用户实现写完后，`tests/test_stats.py` 里等它的 3 条**仍然挂着 `skip`**， + 全量照样显示"全绿" —— 正是易错点 9 的复发（已把该条 `复发` 记到 1）。 + 已摘掉 3 条 skip，并新增 6 条：`estimate` 等于原样本统计量、`confidence` 单调性、 + 五种容器一致、自定义 lambda 走回退分支、**不给种子时跟随 `set_seed`**、非法 `method` 抛错。 + - **S 档函数卡片已写**：`docs/algorithms/bootstrap_ci.md`（第 2 张）， + 记了"极差统计量的区间上界必然塌到 `estimate`"这个非平凡性质，以及 5 个待拍板口径 + - 测试：`tests/test_stats.py` 的 bootstrap 段共 **9 条** +- [ ] 剩余 2 项：`monte_carlo` `save_figure` ### 验证结果 | 检查 | venv | conda `mcm-toolkit` | |---|---|---| -| `pytest` | 见下（单列） | **134 passed, 37 skipped（0 xfail）** | +| `pytest` | 见下（单列） | **143 passed, 34 skipped（0 xfail）** | | `import mcm_toolkit` | 0.0.1 / 16 项 | 0.0.1 / 16 项 | | `ruff check` | All checks passed | **All checks passed**（0 error） | | `scripts/doctor.ps1` | 跑通 | — | @@ -221,8 +236,8 @@ （修配置之前是 17 errors，全是"无法解析导入 numpy / pandas"，见决策 29 与 `docs/04-environment.md` 坑 7） - base 亦已复核：numpy 2.4.6 / scipy 1.18.1 / pandas 3.0.6（与两套环境一致） -- **171 条用例：134 passed / 37 skipped / 0 xfail**。收集数按文件：`test_data.py` 52、`test_api_surface.py` 20、 - `test_core.py` 23、`test_solve.py` 9、`test_viz.py` 8、`test_decision.py` 7、`test_stats.py` 47、 +- **177 条用例：143 passed / 34 skipped / 0 xfail**。收集数按文件：`test_data.py` 52、`test_api_surface.py` 20、 + `test_core.py` 23、`test_solve.py` 9、`test_viz.py` 8、`test_decision.py` 7、`test_stats.py` 53、 `test_models.py` 5 - **`xfail` 再次清零**：`correlation_matrix` 的契约冲突已按方案 A 解决（决策 41），那条用例转正 - 上一轮那 7 条 xfail 已全部转正，决策 30 那套信号（XPASS → 摘 mark）跑完过一个完整闭环 @@ -310,17 +325,18 @@ 用户 2026-10-03 提"每个函数都抄一遍入口/校验/交换太麻烦"。能统一的只有**入口**：出口要求"输出与输入同形"，而 `describe` / `correlation_*` / `confidence_interval` / `t_test` 天生不同形，硬套会重演易错点 14。用户明确"不要三个函数""转换回去自己手打"，所以 `as_array` 只返回 `(arr, kind)`，`ArraySpec` / `restore` 都不做。`numeric_only` 默认 `True`（用户拍的板），但**按输出契约选**：同形函数必须显式 `False` | | 45 | `as_array` **不删缺失值**；容器类型/内容转换失败**统一抛 `TypeError`** | 删缺失是**语义**（`df.corr` 成对删、`describe` 整列统计、`paired=True` 的 t 检验必须成对删），在入口层删掉就把行对应关系永久丢了 —— 这正是 `t_test` 待拍板的口径之一。错误口径统一为 `TypeError` 是**Codex 代选、待用户确认**：`normalize(含文本列的表)` 原先由 numpy 抛 `ValueError`，现在是 `TypeError`（内容转不成 float）；`describe` 的 `weights` 同类情况目前仍抛 `ValueError`（决策 39 说不统一就列出来） | +| 46 | `bootstrap_ci` 的随机数**无条件走 `core.get_rng`**；`statistic_name` 用 `getattr(f, "__name__", None)` 再判 `""` | +原码 `random_state=None` 时另开 `np.random.default_rng()`，等于绕开决策 24（`set_seed` 对它无效，区间不可复现）。**Codex 代改**：`get_rng(None)` 本来就是设计好的"复用模块级生成器且推进"，直接传下去即可。`statistic_name` 的 `else "custom"` 是死代码 —— lambda 也有 `__name__`（值恒为 `""`），实测拿不到 `"custom"`，见易错点 17 | `PearsonRResult` / `SignificanceResult` 是 `_make_tuple_bunch` 在**运行时**拼出来的类，存根里没有 `pvalue` 字段：解包会退化成 `object`、属性访问会报"未知"。运行时完全正常（实测三个函数返回的都是 `np.float64`）。收成 `Any` 是**一处**收窄，调用点仍然干净；`# type: ignore` 要挂在每一处，且会连带关掉该行后面的检查。见 `docs/04-environment.md` 坑 10 | 用户 2026-10-03 明确：没考虑到的非法情况"告诉我我自己改"——因为那是设计决策，得由写实现的人形成判断；而 lint / 进度板 / 测试这类机械活再由用户动手只是浪费。这条与决策 25 / 31 配套：25 定"谁写 src"，31 定"Codex 可改的 bug 范围"，39 定"非算法事务的默认归口" | docstring 与 `test_describe_reports_expected_keys` 早就写死了这个方向，实现里多做的 `.T` 是净负担 —— 且不是"风格不同"那么轻：转置后 `.loc[指标名] = <按变量索引的 Series>` 与列名零交集，pandas **不报错**、把四行整行填成 NaN。方向约定以"能按列名对齐"为准，见易错点 11 | 原码写的是 `if offset > 0: raise`，把唯一合法的用例（正 offset 正好抬定义域）全拒了，报错文案却是"offset 必须大于 0"，自相矛盾。改成校验真正的前提：`log` 要求自变量取正。`offset=0` 对全正数据合法（`offset=None` 的自动分支本来就取 0）。**Codex 代选、待用户确认** —— 若原意是"只收正 offset"，改回即可 | ## 下一步 -1. **`bootstrap_ci` 是 L0 里下一个该写的**（**S 档**，`stats/bootstrap.py`）—— - 除了函数本身，还要配一张完整函数卡片（模板 `docs/dev/function-card-template.md`）， - 因为它是"没有解析解时怎么给区间"的通用答案（`docs/03-roadmap.md` §3.1）。 - 注意 `tests/test_stats.py` 里已有 3 条现成的用例在等它（现在挂着 skip），写完就是现成的验收标准。 -2. 之后 L0 剩余：`monte_carlo`（**S 档**）→ `save_figure`（A）。 +1. **`bootstrap_ci` 已完成**（L0 10/12，卡片也写了）。**L0 里下一个是 `monte_carlo`**（**S 档**，`solve/simulate.py`）—— + 同样要配一张函数卡片；它的三个兄弟（`monte_carlo_probability` / `monte_carlo_integral` / `random_walk`） + 都在同一份卡片体系里，先把 `monte_carlo` 的形状定下来。 +2. L0 的最后一项是 `save_figure`（A 档，`viz/export.py`）。 3. L0 跑通后，找一道真题（建议 2025 C）走完整流程，反向沉淀函数。 4. 已拍板的旧口径（留档）：`log_transform` 的 `offset` 维持现状（决策 35）； `describe` / `zscore` 默认丢掉非数值列（决策 37）；变异系数用 `σ / |μ|`（决策 38）； @@ -346,4 +362,12 @@ docstring 与 `test_describe_reports_expected_keys` 早就写死了这个方向 - 内容转不成 float 时统一抛 `TypeError`（决策 45）—— `normalize(含文本列的表)` 原来抛 `ValueError`。 要不要改回 `ValueError`（"取值问题"）？ - 现在 `kind == "frame"` 时你要自己再写一遍 `x.select_dtypes(include=[np.number]).columns` 取列名。 - 嫌烦的话可以让 `as_array` 顺带把选中的列返回（改成三元组），一句话的事。 \ No newline at end of file + 嫌烦的话可以让 `as_array` 顺带把选中的列返回（改成三元组），一句话的事。 + +8. **`bootstrap_ci` 的 5 个待拍板口径**（Codex 一律没动，按决策 39；详见 `docs/algorithms/bootstrap_ci.md` §10）： + - `x` 里的 **NaN** 静默传进统计量 → `nan` 区间（`confidence_interval` / `t_test` 都是先丢掉的，口径不齐） + - `n_resamples <= 0` → 现在抛 `IndexError: index -1 is out of bounds`（应在入口给可读 `ValueError`） + - `confidence` 越界 → 现在由 numpy 抛 `Percentiles must be in the range [0, 100]` + - `method="bca"` → docstring 承诺过、实测 `NotImplementedError`（补实现 ≈25 行，还是从 docstring 删掉？） + - `bootstrap_statistic` 仍是 `NotImplementedError`，而 `bootstrap_ci` 把那段重采样循环**内联**了 —— + 建议抽出去让 `bootstrap_ci` 调它（易错点 16 的变体：同一段算法写了两遍，只实现了一遍） diff --git a/docs/algorithms/README.md b/docs/algorithms/README.md index 3b20b45..7facf69 100644 --- a/docs/algorithms/README.md +++ b/docs/algorithms/README.md @@ -8,7 +8,7 @@ S 档（自己实现的招牌算法）的**函数卡片**放这里，模板见 | 卡片 | 对应函数 | 状态 | |---|---|---| | `detect_outliers.md` | `data.detect_outliers` | ✅ 2026-10-02 | -| `bootstrap_ci.md` | `stats.bootstrap_ci` | ⬜ | +| `bootstrap_ci.md` | `stats.bootstrap_ci` | ✅ 2026-10-03（§10 有 5 个口径待拍板）| | `monte_carlo.md` | `solve.monte_carlo` | ⬜ | | `monte_carlo_probability.md` | `solve.monte_carlo_probability` | ⬜ | | `monte_carlo_integral.md` | `solve.monte_carlo_integral` | ⬜ | diff --git a/docs/algorithms/bootstrap_ci.md b/docs/algorithms/bootstrap_ci.md new file mode 100644 index 0000000..d371c69 --- /dev/null +++ b/docs/algorithms/bootstrap_ci.md @@ -0,0 +1,127 @@ +# bootstrap_ci + +| 项 | 内容 | +|---|---| +| 模块 | `mcm_toolkit.stats.bootstrap` | +| 档位 | S | +| 参考来源 | Efron (1979)；Efron & Tibshirani (1993) *An Introduction to the Bootstrap*； +| | `scipy.stats.bootstrap`（只做对照，没有调用） | +| 状态 | 实现已落地（2026-10-03），§10 有 4 个口径待拍板 | + +## 1. 数学定义 + +给定样本 `x_1..x_n` 与统计量 `T(·)`：从 `x` 中**有放回**抽 `n` 个得到 `x*`，重复 `B` 次， +得 `T*_1..T*_B`。百分位法区间 + +```text +estimate = T(x) # 用**原样本**，不是重采样分布的均值 +alpha = 1 - confidence +low = T*_(alpha/2) # 经验分布的 alpha/2 分位 +high = T*_(1 - alpha/2) +``` + +输入是一个一维样本，输出是「点估计 + 区间」两个标量 —— 这是**非参数**区间， +不假设正态、也不需要 `T` 的标准误解析式。 + +## 2. 适用场景 + +**用**：`T` 没有解析分布（中位数、极差、自定义指标）；样本有偏或分布未知； +想给区间但不想假设正态（`mean ± z·s/√n` 的替代品）。 + +**不用**：`n < 10` —— 经验分布太粗，`x` 能重采样出的不同样本极少； +数据有**强时间相关性**（股票、气温），此时必须用 block bootstrap； +要估的是**极值**（最大值本身就是上界，见 §5）。 + +## 3. 参数与返回值 + +| 参数 | 类型 | 默认 | 含义 | 取值影响 | +|---|---|---|---|---| +| `x` | Any | — | 样本；走 `core.as_array(x, numeric_only=True, ndim=1)` | DataFrame 只取数值列，再摊平 | +| `statistic` | Callable | `np.mean` | 任意 `f(sample) -> float` | 换 `np.median` 就得到中位数的区间 | +| `confidence` | float | 0.95 | 置信水平 | 越大区间越宽 | +| `n_resamples` | int | 10000 | 重采样次数 `B` | 越大端点越稳；10000 通常够 | +| `method` | str | `"percentile"` | 区间算法 | `"bca"` 尚未实现，传了抛 `NotImplementedError` | +| `random_state` | int \| None | `None` | 随机种子 | `None` 时复用 `core.get_rng` 的模块级生成器且会推进 | + +返回 `Result`： + +| 字段 | 内容 | +|---|---| +| `params` | `estimate` / `low` / `high` / `confidence` / `n_resamples` | +| `metrics` | `statistic_name`（`np.mean` → `"mean"`；lambda → `"custom"`）| +| `data` | 不使用 | + +## 4. 复杂度 + +- 时间：`O(B · n)`（`B` 次重采样，每次都要在 `n` 个值上算一次 `T`） +- 空间：`O(B)` —— **按块处理**，每批 `batch_size = 1000` 行，绝不落地 `B × n` 的整矩阵 +- **什么时候会跑不完**： + - `n = 1e4, B = 1e4` → `1e8` 次操作，秒级，正常； + - `n = 1e6, B = 1e4` → `1e10` 次操作，**跑不完**；而且 1000×1e6 的单批索引矩阵就要 8 GB， + 此时应改用解析近似，或把 `B` 降到 1e3 以下、换更便宜的 `T`。 + +## 5. 数值稳定性 + +- **统计量有上界时，百分位区间会「塌」到上界。** 极差 `max - min` 的重采样值**永远不可能超过** + 原样本极差，`B` 个值里有相当比例正好等于它（抽到原样本极值点的那些）， + 于是 `high == estimate` 精确成立。这不是 bug，是 bootstrap 的固有性质 —— + 区间端点是**样本**的经验分位数，不能给出样本之外的结论。要估极值请用极值理论。 +- 分块只是省内存，不改结果：每批都从同一个 `rng` 按顺序抽索引，等价于一次抽 `B × n`。 +- 随机数**只来自 `core.get_rng`**（决策 24）：给了种子可复现，不给则跟随 `set_seed`。 +- **NaN 目前不处理**：`x` 里的 NaN 会直接进 `T`，`np.mean` 于是给出 `nan` 区间。见 §10。 + +## 6. 常见错误 + +1. 不写 seed → 每次跑出来的区间都不一样，论文里的数字对不上 → 传 `random_state` 或先 `set_seed`。 +2. 以为 bootstrap 能得出**比样本更极端**的结论 → 区间端点就是样本的经验分位数， + "95% 置信上限 173.2cm" 不可能超过样本最大值。 +3. 拿它估最大值 / 最小值 → 区间会塌成一个点（原因同 §5），该换极值理论。 +4. `n < 10` 硬跑 → 能抽出的不同重采样样本太少，区间没有统计意义。 +5. 把 `n_resamples` 当成"样本量"写进论文 → 论文里要同时写清 `n`、`B`、seed、method。 + +## 7. 最小示例 + +```python +import numpy as np + +from mcm_toolkit.stats import bootstrap_ci + +r = bootstrap_ci(heights, statistic=np.median, n_resamples=10000, random_state=42) +print(r.params["low"], r.params["estimate"], r.params["high"]) +print(r.summary()) # 可直接贴进论文的 Markdown 表 +``` + +## 8. 单元测试要点 + +| 测什么 | 期望 | 对应测试 | +|---|---|---| +| 大样本正态，`statistic=np.mean` | 区间半宽 ≈ `2 × 1.96 × se`，容差 10% | `test_bootstrap_ci_matches_normal_theory_on_large_n` | +| 任意统计量（中位数，对数正态样本） | 区间包含样本中位数 | `test_bootstrap_ci_works_for_arbitrary_statistic` | +| 同种子两次 | `low` 完全相等 | `test_bootstrap_ci_is_reproducible` | +| 不给种子 | 跟随 `set_seed`（决策 24） | `test_bootstrap_ci_follows_get_rng_when_no_seed_is_given` | +| `estimate` | 等于原样本上的 `T(x)` | `test_bootstrap_ci_reports_the_original_estimate` | +| `confidence` 0.80 → 0.95 | 区间单调变宽 | `test_bootstrap_ci_interval_widens_with_confidence` | +| 极差作为统计量 | `high <= estimate`（原子点塌陷） | `test_bootstrap_ci_works_for_a_custom_lambda` | +| 五种容器 | DataFrame/Series/ndarray/list/tuple 结果一致 | `test_bootstrap_ci_accepts_the_documented_containers` | +| `method` 非法 / `"bca"` | `ValueError` / `NotImplementedError` | `test_bootstrap_ci_rejects_unknown_method` | + +## 9. MCM/ICM 用法 + +- **哪类题**：任何要"给结论 + 给不确定性"的题 —— 评价类（排名稳定性）、预测类（区间预测）。 +- **论文里怎么表述**："采用非参数 bootstrap（`B = 10000`，`seed = 42`）估计中位数的 95% 置信区间。" +- **图**：预测曲线 + 阴影置信带（`viz.plot` / `viz.save_figure`）。 +- **踩坑提醒**：必须报 `B` 和 seed；说清用的是 percentile 还是 BCa（目前只有 percentile）； + 被问到"为什么不用 ±1.96σ"时就答"统计量是中位数，没有解析标准误"。 + +## 10. 待拍板（未定稿，Codex 不代选） + +| # | 口径 | 现状 | 可选 | +|---|---|---|---| +| 1 | `x` 里的 NaN | 静默传进 `T` → `nan` 区间 | 像 `confidence_interval` 一样丢掉 / 直接报错 / 维持现状 | +| 2 | `n_resamples <= 0` | 走到 `np.percentile` 抛 `IndexError: index -1 is out of bounds` | 入口 `ValueError（必须 ≥ 1）` | +| 3 | `confidence` 越界 | 由 numpy 抛 `ValueError: Percentiles must be in the range [0, 100]` | 入口 `ValueError + 可读文案` | +| 4 | `method="bca"` | docstring 承诺过，实测 `NotImplementedError` | 补 BCa 实现（S 档，约 25 行）/ 从 docstring 删掉 | +| 5 | `bootstrap_statistic` | 仍是 `NotImplementedError`，而 `bootstrap_ci` 把重采样循环内联了 | 抽出去让 `bootstrap_ci` 调它（易错点 16 的变体） | + +> 口径 5 尤其值得做：`bootstrap_statistic` 的 docstring 与 `bootstrap_ci` 里那段循环说的是同一件事， +> 现在等于把同一段算法写了两遍，只实现了一遍。 diff --git a/docs/dev/personal-pitfalls.md b/docs/dev/personal-pitfalls.md index f3e2a79..54fa4f8 100644 --- a/docs/dev/personal-pitfalls.md +++ b/docs/dev/personal-pitfalls.md @@ -38,6 +38,7 @@ | 14 | 把"能看的表"和"能算的表"塞进同一个返回值 | 2026-10-03 | | 15 | 区间被无条件夹到 `[0, 1]`，且 `min` / `max` 参数写反 | 2026-10-03 | | 16 | 同一段入口转换抄了 N 遍，还当它是"本项目的写法" | 2026-10-03 | +| 17 | `hasattr(f, "__name__")` 挡不住 lambda，`else` 分支是死代码 | 2026-10-03 | --- @@ -151,7 +152,9 @@ - 以后怎么防：每次写完实现，用 `pytest -q -rs`（或 `--collect-only`）**核对 skip 条数掉没掉**； 条数没变就说明实现没被验证过。另一个信号是 `pytest --collect-only` 的收集数 —— 新写的测试没进收集列表，就是名字或位置不对。 -- 首次：2026-10-02｜复发：0 +- 复发（2026-10-03）：`bootstrap_ci` 写完了，`tests/test_stats.py` 里等它的 3 条仍然挂着 skip， + 全量照样 134 passed / 37 skipped 全绿 —— 连 `_NOT_IMPLEMENTED` 这个常量名都没提示你。 +- 首次：2026-10-02｜复发：0 → **1（2026-10-03）** ### 10. 把真报错当成"编辑器误报" @@ -287,3 +290,22 @@ ③ 抽出来的函数只干**一件事**：把容器转成数组、并告诉调用方原来是什么。 选列 / 摊平 / 删缺失这些是**语义**，必须留在各函数里（否则重演第 15 条与 `paired=True` 的 NaN 错位）。 - 首次：2026-10-03｜复发：0 + +### 17. `hasattr(f, "__name__")` 判断"是不是具名函数"—— lambda 也有 `__name__` + +- 现象：`bootstrap_ci` 想给 `metrics["statistic_name"]` 一个名字，写成了 + `statistic.__name__ if hasattr(statistic, '__name__') else "custom"`。 + 本意是"认得的内置统计量给名字、自定义的标 custom"，实际传 lambda 时得到的是 + **字符串 `""`** —— 因为 lambda 有 `__name__`，只是值恒为 `""`。 + 也就是说 `else` 那支**几乎永远走不到**，是死代码。 +- 为什么会错：`__name__` 在脑子里等于"函数的名字"，于是 `hasattr` 就等于"这是个函数"。 + 实际上 `__name__` 是**属性存在性**，不是"名字有意义"：lambda、`functools.partial` 包过的函数、 + `numpy` 的 ufunc（`np.mean.__name__` 是 `"mean"`，能过）各有各的情况， + 真正拿不到的只有**没有 `__name__` 的可调用对象**（比如带 `__call__` 的类实例）。 + 同类陷阱还有 `hasattr(x, "__len__")` 判容器（`str` 也有）、`hasattr(x, "__iter__")` 判列表。 +- 以后怎么防：① 判"是不是我想要的那类函数"，用**显式白名单**（`if statistic in [np.mean, np.median, ...]`）， + 别用属性存在性推断语义；② 真要取名字就写成 + `name = getattr(f, "__name__", None)` 再判断 `name in (None, "")`， + 把"没有名字"和"名字是占位符"两种情况都覆盖；③ 每个 `else`/兜底分支都问一句 + "**这条到底测得到吗**"——走不到的分支就是 bug 伪装成的防御性代码。 +- 首次：2026-10-03｜复发：0 diff --git a/src/mcm_toolkit/stats/bootstrap.py b/src/mcm_toolkit/stats/bootstrap.py index fc522df..b022570 100644 --- a/src/mcm_toolkit/stats/bootstrap.py +++ b/src/mcm_toolkit/stats/bootstrap.py @@ -10,6 +10,10 @@ from __future__ import annotations from collections.abc import Callable from typing import Any +import numpy as np + +from ..core import Result, as_array, get_rng + __all__ = ["bootstrap_ci", "bootstrap_statistic"] @@ -42,23 +46,46 @@ def bootstrap_ci( Parameters ---------- + x: + 样本：`pandas.DataFrame` / `pandas.Series` / `numpy.ndarray` / list / tuple。 + 走 `core.as_array(x, numeric_only=True, ndim=1)` —— DataFrame 只取数值列，再摊平成一串样本。 statistic: 任意把样本映射成标量的函数，默认 `np.mean`。 换成 `np.median` / `np.std` / 自定义的 `lambda s: s.max() - s.min()` 都能用， 这是 bootstrap 相对于 `mean ± z·s/√n` 的核心优势。 + confidence: + 置信水平，取 `(0, 1)` 开区间，默认 0.95。 + n_resamples: + 重采样次数，默认 10000。论文里必须写明。 method: `"percentile"`：直接取经验分布的 α/2 与 1-α/2 分位（最常用，默认）； - `"bca"`：偏差校正 + 加速，区间更准但慢，样本有偏时用。 + `"bca"`：偏差校正 + 加速，区间更准但慢，样本有偏时用 —— **尚未实现**，传了抛 + `NotImplementedError`。 + random_state: + 随机种子。传 `None` 时复用 `core.get_rng` 的**模块级生成器且会推进**， + 所以同一段脚本里连调两次结果不同；要完全可复现就显式给种子，或先 `core.set_seed`。 Returns ------- - `Result`：`params={"estimate", "low", "high", "confidence", "n_resamples"}`。 + `Result`：`params={"estimate", "low", "high", "confidence", "n_resamples"}`； + `metrics["statistic_name"]` 是统计量的名字（自定义 lambda 时是 `"custom"`）。 + + Raises + ------ + TypeError + `x` 不是受支持的容器，或内容转不成数值（由 `as_array` 抛）。 + ValueError + `DataFrame` 里一列数值都没有；`method` 不在支持列表内。 + NotImplementedError + `method="bca"`。 Notes ----- * 重采样时必须**按块处理**，不要一次生成 `(10000, n)` 的整矩阵（n 大时会爆内存）。 * 样本量 n < 10 时 bootstrap 不可靠，先补数据。 - * 比赛写论文时要注明：n_resamples、seed、用的是 percentile 还是 BCa。 + * 比赛写论文时要注明：n_resamples、seed、用的是 percentile 还是 BCa（目前只实现了 percentile）。 + * **缺失值目前不处理**：`x` 里的 NaN 会直接进统计量，`np.mean` 于是给出 nan 区间。 + `confidence_interval` / `t_test` 都是先丢掉 NaN 的，两边口径尚未统一。 Examples -------- @@ -67,5 +94,57 @@ def bootstrap_ci( >>> r.params["low"] < 3 < r.params["high"] True """ - # TODO(L0) ★ S 档：配完整函数卡片 - raise NotImplementedError + if statistic is None: + statistic = np.mean + + rng = get_rng(random_state) + + arr, _ = as_array(x=x, dtype=float, ndim=1) + + original_estimate = float(statistic(arr)) + + n = len(arr) + batch_size = 1000 + boot_stats = [] + + for i in range(0, n_resamples, batch_size): + current_batch = min(batch_size, n_resamples - i) + indices = rng.integers(0, n, size=(current_batch, n)) + batch_samples = arr[indices] + + if statistic in [np.mean, np.median, np.std, np.min, np.max]: + batch_stats = statistic(batch_samples, axis=1) + else: + batch_stats = np.array([statistic(s) for s in batch_samples]) + + boot_stats.extend(np.asarray(batch_stats).ravel().tolist()) + + boot_stats = np.array(boot_stats) + + if method == 'percentile': + alpha = 1 - confidence + lower = np.percentile(boot_stats, alpha / 2 * 100) + upper = np.percentile(boot_stats, (1 - alpha / 2) * 100) + elif method == 'bca': + # TODO(L1): BCa = 偏差校正 + 加速，属于 S 档、留给用户自己补的算法。 + raise NotImplementedError("BCa 还没实现，先用 method='percentile'") + else: + raise ValueError(f"不支持的方法:{method},仅支持'percentile'") + + # lambda 也有 __name__（值恒为 ""），所以只判断 hasattr 会让 "custom" 这支变成死代码。 + name = getattr(statistic, "__name__", None) + statistic_name = "custom" if name in (None, "") else str(name) + + return Result( + name="bootstrap_ci", + params={ + "estimate": original_estimate, + "low": float(lower), + "high": float(upper), + "confidence": confidence, + "n_resamples": n_resamples + }, + metrics={ + "statistic_name": statistic_name, + } + ) diff --git a/tests/test_stats.py b/tests/test_stats.py index 7c235d7..6b55eed 100644 --- a/tests/test_stats.py +++ b/tests/test_stats.py @@ -26,7 +26,6 @@ def _pearson_p_from_r(r: float, n: int) -> float: return float(2 * sps.t.sf(abs(t_stat), n - 2)) -@pytest.mark.skip(reason=_NOT_IMPLEMENTED) def test_bootstrap_ci_matches_normal_theory_on_large_n(rng): """大样本下 bootstrap 百分位区间应逼近正态解析区间（±10% 容差）。""" from mcm_toolkit.stats import bootstrap_ci @@ -37,7 +36,6 @@ def test_bootstrap_ci_matches_normal_theory_on_large_n(rng): assert abs((r.params["high"] - r.params["low"]) - 2 * 1.96 * se) < 0.1 * 2 * 1.96 * se -@pytest.mark.skip(reason=_NOT_IMPLEMENTED) def test_bootstrap_ci_works_for_arbitrary_statistic(rng): """bootstrap 的核心价值：任意统计量都能给区间（这里用中位数）。""" from mcm_toolkit.stats import bootstrap_ci @@ -48,7 +46,6 @@ def test_bootstrap_ci_works_for_arbitrary_statistic(rng): assert r.params["low"] < np.median(x) < r.params["high"] -@pytest.mark.skip(reason=_NOT_IMPLEMENTED) def test_bootstrap_ci_is_reproducible(rng): from mcm_toolkit.stats import bootstrap_ci @@ -58,6 +55,82 @@ def test_bootstrap_ci_is_reproducible(rng): assert a.params["low"] == b.params["low"] +def test_bootstrap_ci_reports_the_original_estimate(rng): + """`params["estimate"]` 必须是**原样本**上的统计量，不是重采样分布的均值。""" + from mcm_toolkit.stats import bootstrap_ci + + x = rng.normal(10, 2, size=200) + r = bootstrap_ci(x, statistic=np.mean, n_resamples=500, random_state=0) + + assert r.params["estimate"] == pytest.approx(float(np.mean(x))) + assert set(r.params) == {"estimate", "low", "high", "confidence", "n_resamples"} + assert r.params["n_resamples"] == 500 + assert r.metrics["statistic_name"] == "mean" + + +def test_bootstrap_ci_interval_widens_with_confidence(rng): + """置信水平越高区间必须越宽 —— 结构性检查，不依赖任何具体分布。""" + from mcm_toolkit.stats import bootstrap_ci + + x = rng.normal(size=300) + width = [ + bootstrap_ci(x, np.mean, confidence=c, n_resamples=2000, random_state=7).params["high"] + - bootstrap_ci(x, np.mean, confidence=c, n_resamples=2000, random_state=7).params["low"] + for c in (0.80, 0.95) + ] + assert width[0] < width[1] + + +def test_bootstrap_ci_accepts_the_documented_containers(rng): + """DataFrame 只取数值列后摊平；Series / ndarray / list 都收。""" + from mcm_toolkit.stats import bootstrap_ci + + x = rng.normal(size=200) + ref = bootstrap_ci(x, np.mean, n_resamples=500, random_state=0).params + + for wrapped in (pd.Series(x), pd.DataFrame({"a": x}), list(x), tuple(x)): + got = bootstrap_ci(wrapped, np.mean, n_resamples=500, random_state=0).params + assert got["low"] == pytest.approx(ref["low"]) + + +def test_bootstrap_ci_works_for_a_custom_lambda(rng): + """不认识的统计量走逐样本回退分支（这里是极差）。""" + from mcm_toolkit.stats import bootstrap_ci + + x = rng.normal(size=200) + r = bootstrap_ci(x, statistic=lambda s: float(s.max() - s.min()), n_resamples=500, random_state=0) + + # 极差的重采样值**永远不会超过原样本极差**（原样本本身就是上界）， + # 所以百分位区间的上界会正好落在 estimate 上，不能再断言 estimate < high。 + assert r.params["estimate"] == pytest.approx(float(x.max() - x.min())) + assert r.params["low"] > 0 + assert r.params["low"] < r.params["high"] <= r.params["estimate"] + assert r.metrics["statistic_name"] == "custom" + + +def test_bootstrap_ci_follows_get_rng_when_no_seed_is_given(rng): + """不给种子时必须跟随 `get_rng` 的模块级生成器（决策 24），不能另开一个 default_rng。""" + from mcm_toolkit.core import set_seed + from mcm_toolkit.stats import bootstrap_ci + + x = rng.normal(size=200) + set_seed(0) + first = bootstrap_ci(x, np.mean, n_resamples=300).params["low"] + set_seed(0) + second = bootstrap_ci(x, np.mean, n_resamples=300).params["low"] + + assert first == second + + +def test_bootstrap_ci_rejects_unknown_method(rng): + from mcm_toolkit.stats import bootstrap_ci + + with pytest.raises(ValueError, match="bogus"): + bootstrap_ci(rng.normal(size=50), np.mean, method="bogus", n_resamples=10) + with pytest.raises(NotImplementedError): + bootstrap_ci(rng.normal(size=50), np.mean, method="bca", n_resamples=10) + + def test_confidence_interval_matches_scipy(rng): from scipy import stats