From b0c4b4a413e17bafae4b6a077d2cd2232241d79e Mon Sep 17 00:00:00 2001 From: guzhujushi Date: Mon, 5 Oct 2026 11:26:54 +0800 Subject: [PATCH] =?UTF-8?q?feat(solve):=20monte=5Fcarlo=20=E7=B3=BB?= =?UTF-8?q?=E5=88=97=E4=B8=A4=E4=B8=AA=E5=87=BD=E6=95=B0=E8=90=BD=E5=9C=B0?= =?UTF-8?q?=EF=BC=88L0=2011/12=E3=80=81L1=EF=BC=89+=20Wilson=20=E5=8C=BA?= =?UTF-8?q?=E9=97=B4=EF=BC=9B=E4=BF=AE=20z=5Fscore^2=20=E5=BC=82=E6=88=96?= =?UTF-8?q?=20bug?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 用户实现： - solve.monte_carlo（L0 第 11 项，S 档）：向量化采样→仿真→汇总，8 个指标； 批量大小由实现传 sampler(rng, size=n)；结果含 NaN 判失败；std 用 ddof=0。 - solve.monte_carlo_probability（L1，S 档）：event(batch)->bool，Wilson score 区间夹 [0,1]， confidence 记进 params，小概率改成 message 警告。 Codex 只碰工程侧（算法未动，唯一例外见下）： - 真 bug：z_score^2 是异或不是乘方（np.float64 ^ 2 抛 TypeError，被 except 吞成 success=False， 函数实际不可用）-> 改 z_score**2。新增经验第 21 条。 - docstring：补齐 Parameters / Raises，修掉两处跑不通的 Examples，Wald -> Wilson 同步。 - 入口校验按 §4.7 收严 n（bool / np.integer / TypeError 与 ValueError 分工）—— 决策 49 后此类归用户，已记档。 - ruff 九条（I001 / F541 x3 / SIM300 / W291 x2 / W293 x2 等）、格式与文案。 - 测试：摘掉 2 条已实现却仍挂着的 skip（易错点 9），新增 14 条（含 vectorization、 Wilson 对 statsmodels 的对照、NaN / 布尔 / 长度等边界），solve 段共 25 条。 校验：193 收集 / 161 passed / 32 skipped / 0 xfail；ruff 0；pyright 0。 文档：STATUS.md 决策 47-50；docs/dev/personal-pitfalls.md 第 18-21 条；CHANGELOG.md。 --- CHANGELOG.md | 61 +++++++++ STATUS.md | 157 ++++++++++++++++++----- docs/dev/personal-pitfalls.md | 72 ++++++++++- src/mcm_toolkit/solve/simulate.py | 188 ++++++++++++++++++++++++--- tests/test_solve.py | 202 ++++++++++++++++++++++++++++-- 5 files changed, 619 insertions(+), 61 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index bf4a06d..48037ab 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -338,6 +338,67 @@ （`monte_carlo` 的契约 / 向量化要求 / 三条测试建议 / 卡片路径 / 一个待拍板口径）； 验证结果里的 skip 说明由 37 条订正为 **34 条**。 +### Added +- `solve.monte_carlo` 落地（**L0 第 11 项，S 档**，实现由用户本人写）：向量化"采样 → 仿真 → 汇总"， + 返回 `Result`（8 个 metrics：mean / std / median / p05 / p95 / ci_low / ci_high），随机数无条件走 `get_rng`。 + 配套 `tests/test_solve.py` **10 条**测试（含"`simulator` 只被调用一次"这条可执行的向量化断言）。 + +### Fixed +- `monte_carlo` 的**布尔结果**此前让 `np.std` 在 numpy 2 里抛 `numpy boolean subtract ... is not supported`， + 整个函数判失败 —— 而"事件是否发生"是 MC 最常见的用法。入口改为 `np.asarray(..., dtype=float)`（布尔按 0/1 统计）； + 字符串结果也从"ufunc isnan 不支持"变成可读的 `could not convert string to float`。由行为探针挖出。 +- `monte_carlo` 的 docstring `Examples` 跑不通：原采样器返回标量，实测 `success=False` + （"simulator 返回的结果数量(1)必须与模拟次数 n(1000)一致"）；改成向量化写法。 +- `monte_carlo` 的 ruff 三条（I001 合并同源 import、F541 空 f-string ×2）、`success= True` 的空格、 + 失败文案"重采样" → "模拟"（从 `bootstrap_ci` 抄过来忘改）。 + +### Changed +- `monte_carlo` 按 §4.7 收严 `n`：`np.int64` 不再被误拒（注解放宽为 `int | np.integer`）； + `n=True` / 浮点 `n` 抛 `TypeError`；`n<=0` 抛 `ValueError`。 +- 用户第二轮改的三处口径（决策 48）：批量大小改由实现传 `sampler(rng, size=n)`；结果含 NaN 判 `success=False`； + `np.std` 显式 `ddof=0`。docstring 同步改写（`sampler` 契约、`Raises`、`Notes` 的 NaN 口径）。 +- `STATUS.md`：L0 10/12 → **11/12**，新增决策 48；验证计数更新为 **187 收集 / 153 passed / 34 skipped**。 + +### Changed +- 工作约定收窄（**决策 49**）：**入口没收紧 / 少考虑情况 / 不符合简单原则**这三类改动交回用户本人写， + Codex 只报告 + 列选项（判据：这处改动要不要"想一下才对"）。决策 47 里"入口 / 校验重构归 Codex"的部分作废。 +- `docs/dev/personal-pitfalls.md`：新增 **第 18 / 19 / 20 条**（入口校验没照 §4.7 清单逐条过 / 只想着主路径那种值 / + "不简单"的信号是同一个信息要写两遍），第 1 条（类型判断只写 `int`）**复发 +1**，索引补到 20。 +- `STATUS.md`：新增决策 49 与两个待拍板小口径（`random_state` 注解是否全项目对齐、`std` 的 `ddof` 是否统一）。 + +### Added +- `solve.monte_carlo_probability` 落地（**L1，S 档**，实现由用户本人写）：向量化 `event(batch) -> bool 数组`， + 返回 `Result`（`probability` / `ci_low` / `ci_high` / `std_error`），区间是 Wald 区间并夹在 `[0, 1]`（沿用决策 42 口径）。 + `tests/test_solve.py` 新增 5 条 + **摘掉 2 条 skip**（实现已写完却还挂着，易错点 9），`solve` 段共 24 条。 + +### Fixed +- `monte_carlo_probability` 的 docstring `Examples` 跑不通：原采样器返回标量，实测 `success=False` + （"event 返回的结果数量(1)必须与模拟次数 n(2000)一致"）；改成批量采样写法并实测通过。 +- `monte_carlo_probability` 的 ruff 六条：SIM300（`0.0 >= confidence` 写成 Yoda 条件）、W291 ×2、W293 ×2、F541。 +- 该函数 NaN 提示里的 "simulator" 改成 "event"（从 `monte_carlo` 抄来没改）；`metrics` 缩进与 `{"n":n}` 空格。 + +### Changed +- `monte_carlo_probability` 的 docstring 补齐 `Parameters` / `Raises`（原文只有 Returns / Notes / Examples）， + 并把"Wald 区间 + 夹到 [0,1]"与"`p <= 1e-3` 判失败"两条现状如实写进 Returns / Notes。 +- `STATUS.md`：L1 清单加入 `monte_carlo_probability`；验证计数更新为 **192 收集 / 160 passed / 32 skipped**； + 「下一步」新增该函数的 **6 条待办口径**（全部留给用户，Codex 未动）。 + +### Fixed +- **`monte_carlo_probability` 的 `z_score^2` 是异或不是乘方**（Python `^` 语义）：`np.float64 ^ 2` 直接抛 + `TypeError: ufunc 'bitwise_xor' not supported`，被 `except Exception` 吞成 `success=False` —— 函数实际上**整个不可用**， + 6 条测试全红才暴露。改成 `z_score**2`。新增经验第 21 条。 +- `monte_carlo_probability` 的 message 在没有警告时留了个尾逗号；`{"confidence":confidence}` 补空格。 + +### Changed +- 用户第二轮改定（**决策 50**）：区间由 **Wald 换成 Wilson score**（实测 p=1、n=1000 从退化的 `(1.0, 1.0)` + 变成 `(0.99617, 1.0)`；p=0 从"判失败"变成 `(0.0, 0.00383)`）；小概率不再判失败，改为 message 里的警告； + `confidence` 记进 `Result.params`；`callable` 入口校验（现在抛 `TypeError`）；`float()` 统一 `ci_*` 类型； + `n` 注解放宽为 `int | np.integer`。 +- `tests/test_solve.py`：新增 1 条对照测试（区间必须等于 `statsmodels.proportion_confint(method="wilson")`）， + 共 25 条；全量 **193 收集 / 161 passed / 32 skipped / 0 xfail**。 +- `docs/dev/personal-pitfalls.md`：新增第 **21** 条（公式里的 `^` 不是乘方，且被 `except` 吞成 `success=False`），索引补到 21。 +- docstring 的 "Wald 区间" 同步改为 "Wilson score 区间"。 + ## [0.0.1] - 2026-10-01 ### Added diff --git a/STATUS.md b/STATUS.md index cd5f8c1..3adffe9 100644 --- a/STATUS.md +++ b/STATUS.md @@ -18,10 +18,11 @@ - A 档**薄封装**（`rmse`/`mae`/`fit_linear`/`train_test_split`…）算工程，Codex 可直接实现。 - 分地盘按**文件**划，不按函数划（同一文件两人各写一半曾导致用户编辑器保存失败）。 -当前 **L0 10/12**：`Result`、`get_rng`/`set_seed`、`load_table`/`save_table`、`clean_missing`、 -`detect_outliers`、`normalize`/`standardize`、`describe`、`correlation_matrix`、`bootstrap_ci`； -另**提前落地**了 L1 的 `confidence_interval` / `t_test`。剩余 2 项：`monte_carlo`（用户写）、 -`save_figure`（Codex 写）。 +当前 **L0 11/12**：`Result`、`get_rng`/`set_seed`、`load_table`/`save_table`、`clean_missing`、 +`detect_outliers`、`normalize`/`standardize`、`describe`、`correlation_matrix`、`bootstrap_ci`、 +`monte_carlo`（**用户写**，2026-10-04）；另**提前落地**了 L1 的 `confidence_interval` / `t_test` / +`monte_carlo_probability`（**用户写**，2026-10-05）。剩余 1 项：`save_figure`（Codex 写）。 +`monte_carlo` 的口径已拍板（决策 48），`ci_low`/`ci_high` 的语义与 `monte_carlo_probability` 的 6 条待办见"下一步"。 ## 任务清单 @@ -234,13 +235,52 @@ - [x] **工程收尾（本轮）**：`tests/test_solve.py` 的模块级 `pytestmark = pytest.mark.skip(...)` 拆成 **9 条逐条** `@pytest.mark.skip(reason=_NOT_IMPLEMENTED)` —— 模块级 skip 会把整个文件吞掉， 正是易错点 9 的根源（决策 26）。收集数不变（177）、跳过数不变（34），只是粒度变细。 -- [ ] 剩余 2 项：`monte_carlo`（**用户写**）、`save_figure`（Codex 写） +- [x] `solve/simulate.py` 的 `monte_carlo`（**L0 第 11 项，S 档，算法由用户本人写**，2026-10-04） + - 核心逻辑：向量化采样 → 一次仿真 → 8 个指标；`get_rng` 无条件调用（决策 46 已内化） + - 用户第二轮改了三处（见决策 48）：批量大小改由实现传 `sampler(rng, size=n)`； + 结果含 NaN 判失败；`np.std(data, ddof=0)` 写明 + - **Codex 只碰工程侧，算法一行未动**： + ① ruff 三条（I001 合并同源 import、F541 空 f-string ×2） + ② docstring：补 `Parameters`（原文漏 `random_state` / `**kwargs`）与 `Raises`；按新契约改写 + `sampler` 说明（必须能接住 `size`，且 `**kwargs` 里别再传 `size`）；`Notes` 记 NaN 口径； + **修掉跑不通的 Examples**（原版 `lambda rng, **k: rng.normal(0, 1)` 返回标量，实测 `success=False`） + ③ 按 §4.7 收严 `n` 校验：`np.int64` 不再被拒、`n=True` 改抛 `TypeError`、浮点 `n` 归 `TypeError`、 + `n<=0` 归 `ValueError` + ④ 格式与文案：`success= True` → `success=True`；失败文案"重采样" → "模拟"（从 `bootstrap_ci` 抄来没改） + ⑤ **入口统一 `dtype=float`** —— 修掉一个真 bug：布尔结果（"事件是否发生"，MC 最常见的用法） + 会让 `np.std` 在 numpy 2 里抛 `boolean subtract ... is not supported`，整个函数判失败； + 字符串结果也从"ufunc isnan 不支持"变成可读的 `could not convert string to float` + - **测试：10 条**（`tests/test_solve.py`）—— 返回长度与 8 个 key、`simulator` 只被调用一次（向量化的可执行版本）、 + 同种子可复现、省略种子时跟随 `set_seed`、E[X²]=1（n=2e5）、docstring Example 真能跑、`np.int64` 的 n、 + 布尔结果按 0/1 统计、NaN 结果判失败、5 个非法入参 + - 行为探针实测：`simulator` 调用次数 = **1**（真向量化）；20/20 次区间覆盖 0 +- [x] `solve/simulate.py` 的 `monte_carlo_probability`（**L1，S 档，算法由用户本人写**，2026-10-05） + - 核心：向量化 `event(batch) -> bool 数组` → `p = mean(data)` → **Wilson score 区间**并夹到 `[0, 1]` + （比例量纲，沿用决策 42 的口径——同类问题没重犯） + - 用户第二轮（2026-10-05）按 6 条待办改掉 5 条：Wilson 替 Wald、`callable` 入口校验、 + `confidence` 记进 `params`、`float()` 统一 ci 类型、小概率不再判失败（改成 message 里的警告）——见决策 50 + - **Codex 第三轮只碰工程侧**：① `z_score^2` → `z_score**2` —— **真 bug**：`^` 是异或， + numpy 浮点直接抛 `TypeError`，被 `except` 吞成 `success=False`，6 条测试全红才暴露（易错点 21）； + ② message 无警告时的尾逗号；③ `{"confidence":confidence}` 的空格；④ docstring 的 "Wald" → "Wilson"； + ⑤ 新增 1 条对照测试：区间必须等于 `statsmodels.proportion_confint(method="wilson")` + - 行为探针实测：p≈0.5 正确（0.49991 @ n=2e5）、`confidence` 变大区间变宽、 + 同种子可复现、省略种子跟随 `set_seed`、`np.int64` 的 n 可收、`confidence=1.5` 抛 `ValueError` + - **Codex 只碰工程侧，算法一行未动**： + ① ruff 六条（SIM300 Yoda、W291 ×2、W293 ×2、F541） + ② docstring：补 `Parameters`（原文只有 Returns/Notes/Examples）与 `Raises`， + **修掉跑不通的 Examples**（原 `lambda rng, **k: rng.normal(0, 1)` 返回标量，实测 `success=False`）， + 并把"Wald 区间 + 夹到 [0,1]"与"`p <= 1e-3` 判失败"两条**现状**如实写进 Returns / Notes + ③ 格式：`metrics` 字典缩进、`{"n":n, ...}` 空格、行尾空白 + ④ 文案：NaN 提示里的 "simulator" 改成 "event"（从 `monte_carlo` 抄来没改） + - **测试：摘掉 2 条 skip**（`test_monte_carlo_probability_of_standard_normal`、`test_monte_carlo_is_reproducible` + —— 实现写完还挂着 skip，正是易错点 9）并新增 5 条，`solve` 段共 24 条 +- [ ] 剩余 1 项：`save_figure`（A 档 = 工程，**Codex 直接写**） ### 验证结果 | 检查 | venv | conda `mcm-toolkit` | |---|---|---| -| `pytest` | 见下（单列） | **143 passed, 34 skipped（0 xfail）** | +| `pytest` | 见下（单列） | **161 passed, 32 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` | 跑通 | — | @@ -250,13 +290,14 @@ （修配置之前是 17 errors，全是"无法解析导入 numpy / pandas"，见决策 29 与 `docs/04-environment.md` 坑 7） - base 亦已复核：numpy 2.4.6 / scipy 1.18.1 / pandas 3.0.6（与两套环境一致） -- **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、 +- **193 条用例：161 passed / 32 skipped / 0 xfail**。收集数按文件：`test_data.py` 52、`test_api_surface.py` 20、 + `test_core.py` 23、`test_solve.py` 25（+16，含摘掉 2 条 skip）、`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）跑完过一个完整闭环 - `examples/01_clean_and_describe.py` 实测跑通（exit 0）；`00_quickstart.py` 仍等 `stats` / `viz` -- 34 条 skip 是等实现的算法（`train_test_split`、`monte_carlo`、`anova`、`fit_distribution` 等）； +- 32 条 skip 是等实现的算法（`train_test_split`、`monte_carlo_integral`、`anova`、`fit_distribution` 等）； + `monte_carlo_probability` 已落地，摘掉 2 条 skip（本轮净 +3 条用例）； `bootstrap_ci` 已落地摘掉 3 条 skip（本轮新增 6 条用例，净 +3）；`test_solve.py` 9 条按决策 26 改为逐条 skip； `TestValidate` 的 4 条仍挂着 —— `check_*` 还是 `NotImplementedError`，但它们实现时应当直接复用 `as_array` - **本轮回归（已修）**：`transform.py` 上一版有 3 处真 bug —— 修前 **7 failed / 75 passed**， @@ -271,7 +312,7 @@ 只差 Python 补丁号（3.13.9 vs 3.13.15）。详见 `docs/04-environment.md` §5。 - [x] ~~Gitea 远程尚未配置~~ → **已配置并推送完成**，见上面 D 节与决策 19/20。 note：仓库当前是 **Public**（匿名 `git ls-remote` 可读），要改私有去 Gitea 仓库设置。 -- [ ] L0 完成 **10/12**（`Result`、`get_rng`/`set_seed`、`load_table`/`save_table`、`clean_missing`、 +- [ ] L0 完成 **11/12**（`Result`、`get_rng`/`set_seed`、`load_table`/`save_table`、`clean_missing`、 `detect_outliers`、`normalize`、`standardize`、`describe`、`correlation_matrix`、`bootstrap_ci`）， 另 L1 的 `confidence_interval` / `t_test` 已提前落地； **只剩 `monte_carlo`（用户写）与 `save_figure`（Codex 写）**，详情见 E 节与"下一步" @@ -353,37 +394,89 @@ docstring 与 `test_describe_reports_expected_keys` 早就写死了这个方向 工程侧（入口/校验、`tests/`、lint、文档、卡片、`.gitignore`）反过来由 Codex 直接落盘，用户只读； 分地盘按**文件**划，不按函数划。决策 39 的"语义由用户拍板"继续有效。 +| 48 | `monte_carlo` 的四个口径由**用户拍板**：实现传 `size=n`；结果含 NaN 判失败；`std` 用 `ddof=0`；出错仍统一返回 `success=False` | +用户 2026-10-04 直接改在代码里。① 批量大小由实现传（`sampler(rng, size=n, **kwargs)`）：调用方少写一次， +代价是**每个 sampler 都必须能接住 `size`**，且 `**kwargs` 里不能再传 `size`（会重复传参）—— 已写进 docstring。 +② 结果里只要出现 NaN 就 `success=False`：MC 的 NaN 多来自除零 / 模型发散，静默统计只会得到一列 nan 指标。 +③ `np.std` 默认 `ddof=0`，现在显式写明（与 `describe` 那条线的 `ddof=1` 不一致，属已知差异）。 +④ 计算期异常仍全部包成 `success=False`（保留 `try/except Exception`）；代价是自己写错代码时没有 traceback。 +配套 Codex 代改：入口 `np.asarray(..., dtype=float)` —— 布尔结果此前会让 `np.std` 抛错（真 bug，已修）。 +**唯一未定**：`ci_low`/`ci_high` 保持"结果的分位数"（现状）还是改成"均值的置信区间"。 + +| 49 | **这类改动归用户本人**：入口没收紧、少考虑情况、不符合简单原则 —— 用户 2026-10-04 明确"交给我改，不要帮我改" | +用户原话："以后像这种入口没收紧、少考虑情况、不符合简单原则的情况交给我改，不要帮我改，把这几条经验写入经验文档"。 +这是决策 47 的**收窄**：47 把"入口 / 校验重构"划给了 Codex，49 把它收回 —— 凡是**带判断、取舍、设计成分**的， +一律由用户写，Codex 只报告 + 列选项（与决策 39 同一逻辑）。判据一句话：**这处改动要不要"想一下才对"** —— +要，就留给用户；不要（ruff、格式化、docstring 缺节、`tests/`、`STATUS` / `CHANGELOG`、函数卡片）才由 Codex 顺手补。 +本轮实例（本该留给用户，Codex 已代改、用户可回退重写）：`monte_carlo` 的 `n` 按 §4.7 收严；入口 `np.asarray(..., dtype=float)`。 +经验已写入 `docs/dev/personal-pitfalls.md` **第 18 / 19 / 20 条**，并把第 1 条（类型判断只写 `int`）**复发 +1**。 + +| 50 | `monte_carlo_probability` 的区间由 **Wald 改 Wilson score**；小概率不再判失败（改成 message 警告）；`confidence` 记进 `params` | +用户 2026-10-05 按"6 条待办"自行改定。Wilson 在 p 贴近 0 / 1 时不会退化成零宽：实测 p=1、n=1000 从 `(1.0, 1.0)` +变成 `(0.99617, 1.0)`；p=0 从"判失败"变成 `(0.0, 0.00383)`（原 Wald 是 0）。区间口径与 +`statsmodels.proportion_confint(method="wilson")` 对齐，已加成对照测试。 +**遗留两条**：`event: Callable[..., bool]` 注解与"吃整批、返回数组"不符；警告文案把 p≈1 也说成"小概率事件"。 +另：`z_score^2` 那处是 Codex 修的真 bug（`^` 是异或，见易错点 21）——按决策 31 属纯 bug，未改任何语义。 + ## 下一步 -### ⭐ 明天从这里开始（第一个要写的文件） +### ⭐ 当前进度与待办（2026-10-05） -**文件：`src/mcm_toolkit/solve/simulate.py` → 写 `monte_carlo`（L0 第 11/12 项，S 档）。** -同一文件里 `monte_carlo_probability` / `monte_carlo_integral` / `random_walk` 先别动 —— 一次一个函数。 +**本轮：`monte_carlo_probability` 已由用户写完并改到第二版**（L1，S 档）。Codex 两轮都只做工程侧， +**算法逻辑一行未动**（唯一例外是一个纯 typo：`z_score^2` → `z_score**2`，见易错点 21 / 决策 50）。 -**① 参考契约（docstring 里已定，签名别改）** -- `monte_carlo(sampler, simulator, n=10000, random_state=None, **kwargs) -> Result` -- `sampler(rng, **kwargs) -> 参数`：一次抽**一整批**参数（可以是标量/数组/元组，随便） -- `simulator(params) -> float`：单次仿真结果 -- 返回 `Result`：`data` = `n` 个结果组成的 `np.ndarray`； - `metrics = {"mean", "std", "median", "p05", "p95", "ci_low", "ci_high"}` +**接下来三件事** +1. `monte_carlo_probability` 还剩 **2 条待办**（下面一节）——都是你的活。 +2. `ci_low` / `ci_high` 的语义确认一句（也是你的活），确认后 Codex 补 `docs/algorithms/monte_carlo.md` 卡片。 +3. L0 最后一项 `save_figure`（A 档 = 工程，`viz/export.py`）—— **Codex 直接实现**，你只读。 -**② 必须向量化**：`sampler` 一次抽 `n` 份（如 `rng.normal(0, 1, size=n)`），`simulator` 对整个参数数组算一次， -**绝不在 `simulator` 里写 Python 循环**。10^4 次循环还能忍，10^6 次直接放弃 —— 这是"能用"和"不能用"的分界线。 +**你下一个要写的（S 档，仍在 `solve/simulate.py`）**：`monte_carlo_integral` → `random_walk`。 +`monte_carlo_integral` 其实就是 `monte_carlo` 的特例：自己造 `n×d` 的均匀样本 → 逐点求值 → +`mean(f) × 体积`，正好把那套"造批 → 一次算 → 汇总"的写法复用第二次。 -**③ 随机数无条件走 `core.get_rng`**（决策 46 的教训，别重演 bootstrap 那个坑）： -`rng = get_rng(random_state)`，不要自己另开 `np.random.default_rng()`，否则 `set_seed` 管不到它、区间不可复现。 +### ✅ monte_carlo_probability 的 6 条待办（2026-10-05 已办 5 条） -**④ 先写一条会失败的测试**（roadmap §3.2 的纪律）。`tests/test_solve.py` 现在**没有**直接测 `monte_carlo` 的用例 -（`test_monte_carlo_is_reproducible` 名字像，实际调的是 `monte_carlo_probability`）。建议三条： -- 固定种子两次结果完全相同（可复现） -- `sampler` 恒返回常数 → `data` 全等于该常数、`std == 0`（常量边界） -- `n` 增大时 `mean` 的波动按 `O(1/√n)` 下降（可留成大样本慢测或标 `slow`） +- ~~① `p <= 1e-3` 判 `success=False`~~ → 改成 `warning_message`，仍返回区间（实测 p=0 → `(0.0, 0.00383)`）。 +- ~~② Wald 在 p 靠近 1 时退化~~ → 换 **Wilson**（p=1、n=1000：`(1.0, 1.0)` → `(0.99617, 1.0)`）。 +- ~~③ `call_params` 没记 `confidence`~~ → 已加。 +- ~~④ `event` / `sampler` 不可调用被吞成 `success=False`~~ → 已加 `callable` 校验，现在抛 `TypeError`。 +- ~~⑤ `metrics` 混 `np.float64`~~ → `float()` 统一（实测 4 个 key 全是 `float`）。 +- ~~⑥ `n: int`~~ → 已改 `int | np.integer`。 -**⑤ 还欠一张 S 档卡片**：`docs/algorithms/monte_carlo.md`（模板 `docs/dev/function-card-template.md`， -照抄 `docs/algorithms/bootstrap_ci.md` 的骨架），写完把 `docs/algorithms/README.md` 的进度表打勾。 +**还剩 2 条（你的活）** +1. `event: Callable[..., bool]` 的注解与行为不符：实际必须吃**整批**参数、返回长度 `n` 的数组。 +2. 警告文案把 **p ≈ 1** 也说成"小概率事件 P < 1e-3 不建议使用"——那一条建议对 p→1 是错的 + （临界情形该说的是"改用 1−p 或换方法"）。条件是你写的 `p_value <= 1e-3 or p_value >= 1 - 1e-3`，文案只覆盖了前一半。 -**⑥ 等你拍板的一个口径（Codex 不代选，决策 39/47）**：`ci_low` / `ci_high` 用**分位数法** -（`np.percentile(data, [2.5, 97.5])`）还是**正态近似**（`mean ± 1.96·std/√n`）？两个都能用，选一个我记进卡片。 +### ⏳ monte_carlo 只剩 1 个口径要你确认 + +**`ci_low` / `ci_high` 是不是就按"结果的分位数"定？** + +现状 `np.percentile(data, [2.5, 97.5])`。实测 n=2e5、E[X²]=1 给出 `[0.001, 5.0455]` —— +这是**模拟结果的分散度**（预测区间），与 `p05` / `p95` 同类，只是分位点不同。 +而"均值估计的置信区间"（MC 的统计误差）是 `mean ± 1.96·std/√n` = `[0.9962, 1.0086]`， +两者差两个数量级，论文话术也完全不同： + +- A（现状）：写"95% 的模拟结果落在 [a, b]"（描述风险分布） +- B：改成均值 CI → 写"估计值 1.0024 ± 0.0043（95% CI）"（证明 n 够大） +- C：两者都要 → 指标要改名（如 `p025`/`p975` + `mean_ci_low`/`mean_ci_high`） + +注意：`monte_carlo_probability` 里的 `ci_*` 已经是**估计量的 CI**（Wald 区间），所以这两个函数的 `ci_*` +目前**含义不同**——这正是要给它们定口径的原因之一。 + +### ✅ 已拍板（决策 48） + +- 批量大小由实现传 `size=n`；`sampler` 必须能接住 `size`，`**kwargs` 里不要再传 `size`。 +- 结果含 NaN → `success=False`（不落 metrics）。 +- `std` 用 `ddof=0`。 +- 出错仍统一返回 `success=False`（`try/except Exception` 保留）。 + +### ⏳ 另外两个小口径（你说了算） + +- **`random_state` 的注解**：经验文档第 3 条写的标准是 `int | np.integer | np.random.Generator | None`， + 但 `stats/bootstrap.py`、`solve/simulate.py`、`stats/inference.py` 现在都写 `int | None`（只有 `core/random.py` 是全的）—— + 要不要全项目对齐？（运行期 `get_rng` 收 `np.int64`，所以是注解与行为不一致，不是功能 bug） +- **`std` 的 `ddof`**：`monte_carlo` 用 `ddof=0`，`describe` 那条线用 `ddof=1`。要不要统一？ ### 其余待办 diff --git a/docs/dev/personal-pitfalls.md b/docs/dev/personal-pitfalls.md index 54fa4f8..a4eec79 100644 --- a/docs/dev/personal-pitfalls.md +++ b/docs/dev/personal-pitfalls.md @@ -39,6 +39,10 @@ | 15 | 区间被无条件夹到 `[0, 1]`，且 `min` / `max` 参数写反 | 2026-10-03 | | 16 | 同一段入口转换抄了 N 遍，还当它是"本项目的写法" | 2026-10-03 | | 17 | `hasattr(f, "__name__")` 挡不住 lambda，`else` 分支是死代码 | 2026-10-03 | +| 18 | 入口校验"想起什么写什么"，没照 §4.7 清单逐条过 | 2026-10-04 | +| 19 | 只想着主路径那种值，忘了布尔 / NaN / 长度不对 | 2026-10-04 | +| 20 | "不简单"的信号：同一个信息要在调用方写两遍 | 2026-10-04 | +| 21 | 公式里的 `^` 不是乘方，而且被 `except` 吞成 success=False | 2026-10-05 | --- @@ -50,7 +54,10 @@ - 以后怎么防：写类型分支时默认假设"这个值可能来自 numpy"。判断用 `isinstance(x, (int, np.integer))`， 取值统一 `int(x)` 归一化；并且给所有非预期类型加 `raise TypeError` 兜底—— **宁可报错，也不要静默走默认分支**。 -- 首次：2026-10-01｜复发：0 +- 复发（2026-10-04）：`monte_carlo` 的 `n` 只写了 `isinstance(n, int)` —— `np.int64(1000)` 被 `ValueError` **误拒**， + 而 `n=True` 反而放行（`bool` 是 `int` 子类，判断顺序又反了）。 + 同一张 §4.7 表的第一、三行，症状从"静默兜底"变成"误拒 + 误收"：**规则没换，只是没去查**（见第 18 条）。 +- 首次：2026-10-01｜复发：0 → **1（2026-10-04）** ### 2. 以为 `# TODO(L0)` 只是普通注释 @@ -309,3 +316,66 @@ 把"没有名字"和"名字是占位符"两种情况都覆盖；③ 每个 `else`/兜底分支都问一句 "**这条到底测得到吗**"——走不到的分支就是 bug 伪装成的防御性代码。 - 首次：2026-10-03｜复发：0 + +### 18. 入口校验"想起什么写什么"，没有照着清单逐条过 + +- 现象：`monte_carlo` 的入口只写了"n 是整数、n 为正"，没写"那 `n=True` / `n=np.int64` / `n=10.0` 各走哪条"。 + 结果 `np.int64(1000)` 被误拒、`n=True` 被放行、浮点 `n` 报的是 `ValueError` 而不是 `TypeError`。 + §4.7 那张判断顺序表早就写死了，写代码时一次都没翻。 +- 为什么会错：校验是"防御性代码"，写的时候满脑子是主路径（"次数当然是整数"）， + **不会自然想到 `True`、`np.int64`、`None`、0 这些"合法类型里的例外值"**。 + 清单存在的意义正是替代这种直觉——但清单要主动去查才有用；不查就等于没有。 + 同类：第 1 条（只写 `int`）、第 7 条（守卫子句的位置）、`random_state=None` 绕开 `get_rng`（决策 46）。 +- 以后怎么防：每个公开函数写完，把 §4.7 的表当 checklist **逐行过一遍** + （int + np.integer / bool 排在 int 前 / float + np.floating / 数组先 `np.asarray`）， + 再用三个探针各打一枪：`f(x=np.int64(1))`、`f(x=True)`、`f(x=1.0)`。 + **这一步是本人的活**（决策 49）：要不要放宽、报什么错、算不算合法，都是取舍，不是机械活。 +- 首次：2026-10-04｜复发：0 + +### 19. 只想着主路径上的那种值，忘了同一位置还会出现别的类型 / 值 + +- 现象：`monte_carlo` 的 `simulator` 返回**布尔**数组（"事件是否发生"，MC 最常见的用法）时， + `np.asarray(...)` 原样保留 bool，`np.std` 在 numpy 2 上直接抛 + `numpy boolean subtract ... is not supported`，整个函数判失败；字符串结果则给出 + "ufunc isnan not supported"这种看不懂的话。同一轮还漏了：结果里含一个 NaN 时 8 个指标全 NaN 却报 `success=True`。 +- 为什么会错：写的时候脑子里只有"年收益 = 浮点数"这一种结果， + "0/1 事件"、"NaN"、"长度对不上"都不在主路径上。而**这个函数恰恰是这几类输入的高发区** + （事件概率、模型发散、批量大小写错），所以它们不是边角料，是常规用法。 +- 以后怎么防：写完一个"吃调用方数组"的函数，除了主路径，按五类各跑一次： + **布尔 / 整数 / 含 NaN / 全相同（常量）/ 长度不对**。 + 类型归一化放在**入口一处**做（`np.asarray(x, dtype=float)`），后面所有指标就都安全了—— + 这正是 §4.7 表最后一行与 §4.8 存在的理由。 +- 首次：2026-10-04｜复发：0 + +### 20. "不简单"的最早信号：同一个信息要在调用方写两遍 + +- 现象：第一版 `monte_carlo` 的契约是 `sampler(rng, **kwargs)`，但采样器**拿不到 n**， + 于是调用方必须写 `n=1000, size=1000` —— 同一个数写两遍，还要靠 `len(data) != n` 事后兜底。 + 改成实现自己传 `sampler(rng, size=n, **kwargs)` 之后，调用方只写 `n=1000`。 + 同一处的后遗症：`**kwargs` 里再传 `size` 会撞成 `got multiple values for keyword argument`。 +- 为什么会错：设计签名时只想着"参数全靠 `**kwargs` 透传"这一条好处（灵活）， + 没数"调用方为了让它跑起来还得额外写什么"。 + **重复的那一个数是症状，不是问题**：真正的问题是"批量大小"这个信息在接口里**没有归属**。 +- 以后怎么防：写完签名，立刻写出**最小调用示例**，然后数两件事： + ① 有没有哪个信息写了两遍；② 有没有哪一步"只有我知道"。有 → 把它吸收进实现，而不是写进文档。 + 简单 ≠ 参数少，而是**调用方不用重复表达**。 +- 首次：2026-10-04｜复发：0 + +### 21. 把数学公式抄进代码时，`^` 当成乘方 —— 还被 `except` 吞成 `success=False` + +- 现象：`monte_carlo_probability` 的 Wilson 区间里写 `z_score^2 / n`。Python 的 `^` 是**按位异或**， + 对 `np.float64` 直接抛 `TypeError: ufunc 'bitwise_xor' not supported`；因为整段计算在 + `try/except Exception` 里，对外只表现为 `success=False` + 一句 message —— + 最后是 6 条测试全红（`KeyError: 'ci_high'`）才把它挖出来。 +- 为什么会错：纸上的 `z²`、`x²` 抄进代码时手会自动敲 `^`（数学记号、Excel、LaTeX 都是这个符号）， + 而 Python 的乘方是 `**`。更险的是 `^` 在**整数**上是合法异或：`3^2` 得到 1、不报错，只算错。 + 第二层原因是 `try/except Exception` 把"代码写错"和"数据不合法"混成一个出口， + 错误现场（哪一行、哪个异常）被压成一句字符串，静态检查也看不见（pyright 当时 0 error）。 +- 以后怎么防：① 公式里的幂一律写 `**`（或 `np.power`），见到 `^` 就当 bug 处理； + ② 新函数写完先**手动跑一次正常输入**，`success` 必须是 `True` —— 这条一跑就现形； + ③ 给算法配一条**独立实现的对照测试**（Wilson 对 `statsmodels.proportion_confint(method="wilson")`）， + 这类"公式抄错"会被直接钉死；④ 记住 `except Exception` 越宽，错误越难看见。 +- 首次：2026-10-05｜复发：0 + + + diff --git a/src/mcm_toolkit/solve/simulate.py b/src/mcm_toolkit/solve/simulate.py index cedadd6..2047b3c 100644 --- a/src/mcm_toolkit/solve/simulate.py +++ b/src/mcm_toolkit/solve/simulate.py @@ -10,6 +10,11 @@ from __future__ import annotations from collections.abc import Callable from typing import Any +import numpy as np +from scipy.stats import norm + +from mcm_toolkit.core import Result, get_rng + __all__ = [ "markov_simulate", "markov_steady_state", @@ -24,7 +29,7 @@ __all__ = [ def monte_carlo( sampler: Callable[..., Any], simulator: Callable[..., Any], - n: int = 10000, + n: int | np.integer = 10000, random_state: int | None = None, **kwargs: Any, ) -> Any: @@ -33,60 +38,209 @@ def monte_carlo( Parameters ---------- sampler: - `f(rng, **kwargs) -> 一组参数`，例如从先验分布里抽参数。 + `f(rng, size, **kwargs) -> 一批参数`。**一次就要抽够 `n` 组**， + 例如 `lambda rng, size: rng.normal(0, 1, size)`。 + 批量大小由本函数以 `size=n` 传进去：`sampler` 必须能接住 `size` + （显式声明，或带 `**kwargs`）。**别再往 `**kwargs` 里塞 `size`**，否则重复传参。 simulator: - `f(params) -> float`，例如"给定参数下的年收益"。 + `f(params) -> 一批结果`。收到 `sampler` 的那一整批参数， + 返回长度同为 `n` 的数组（标量 / 长度对不上都算契约错误）。 + **内部不要写 Python 循环**。 n: - 模拟次数。**必须向量化**：`simulator` 内部不要写 Python 循环， - 否则 10^4 次能跑几分钟，10^6 次直接放弃。 + 模拟次数，正整数。与 `simulator` 返回的长度必须一致。 + random_state: + 随机种子。`None` 时跟随 `core.set_seed` 的全局流，见 `core.random.get_rng`。 + **kwargs: + 原样转交给 `sampler`，例如分布参数 `mu` / `sigma`（`size` 由本函数自己传，见上）。 + 会一并记进 `Result.params`，便于复现。 Returns ------- - `Result`：`data` 是 `n` 个结果组成的数组， + `Result`：`data` 是 `n` 个结果组成的数组（一律转成 float，布尔结果按 0/1 统计）， `metrics={"mean", "std", "median", "p05", "p95", "ci_low", "ci_high"}`。 + Raises + ------ + TypeError: + `sampler` / `simulator` 不可调用，`n` 不是整数，或 `random_state` 类型非法。 + ValueError: + `n` 不是正整数。 + Notes ----- - 比赛里 Monte Carlo 的收敛速度是 `O(1/√n)`—— - 想少一半误差要 4 倍样本。所以论文里要给"区间随 n 收敛"的图， - 证明 `n` 已经够大（这也是 `viz.plot_convergence` 的典型用法）。 + 比赛里 Monte Carlo 的收敛速度是 `O(1/√n)` —— 想少一半误差就要 4 倍样本。 + 所以论文里要给"区间随 n 收敛"的图，证明 `n` 已经够大 + （这也是 `viz.plot_convergence` 的典型用法）。 + + 结果里只要出现 **NaN** 就判失败（`success=False`）：MC 的 NaN 通常来自除零、 + 取负数的根号或模型发散，静默统计只会得到一列 nan 指标。 Examples -------- >>> import numpy as np - >>> r = monte_carlo(lambda rng, **k: rng.normal(0, 1), np.mean, n=1000, random_state=0) + >>> draw = lambda rng, size, **k: rng.normal(0, 1, size) + >>> r = monte_carlo(draw, lambda x: x ** 2, n=10000, random_state=0) + >>> bool(abs(r.metrics["mean"] - 1.0) < 0.1) # E[X²] = 1 + True """ - # TODO(L1) ★ S 档 - raise NotImplementedError + if not callable(sampler) or not callable(simulator): + raise TypeError("sampler 和 simulator 必须是可调用函数") + if isinstance(n, bool) or not isinstance(n, (int, np.integer)): + raise TypeError(f"n 必须是整数，你传入了 {type(n).__name__}") + if int(n) <= 0: + raise ValueError(f"n 必须是正整数，你传入了 {n}") + n = int(n) + + rng = get_rng(random_state) + call_params = {"n": n, "random_state": random_state, **kwargs} + + try: + params = sampler(rng, size=n, **kwargs) + data = np.asarray(simulator(params), dtype=float).ravel() + if np.isnan(data).any(): + raise ValueError("模拟结果中包含 NaN，请检查 sampler 或 simulator 的传入参数") + if len(data) != n: + raise ValueError(f"simulator 返回的结果数量({len(data)})必须与模拟次数 n({n})一致") + + metrics = { + "mean": float(np.mean(data)), + "std": float(np.std(data, ddof=0)), + "median": float(np.median(data)), + "p05": float(np.percentile(data, 5)), + "p95": float(np.percentile(data, 95)), + "ci_low": float(np.percentile(data, 2.5)), + "ci_high": float(np.percentile(data, 97.5)), + } + + return Result( + name="monte_carlo", + params=call_params, + metrics=metrics, + data=data, + success=True, + message=f"蒙特卡洛模拟成功：已完成 {n} 次模拟", + ) + except Exception as e: + return Result( + name="monte_carlo", + params=call_params, + metrics={}, + data=None, + success=False, + message=f"蒙特卡洛模拟失败：{type(e).__name__} - {e}", + ) def monte_carlo_probability( event: Callable[..., bool], sampler: Callable[..., Any], - n: int = 10000, + n: int | np.integer = 10000, confidence: float = 0.95, random_state: int | None = None, ) -> Any: """估计事件发生概率，并给出区间。 + Parameters + ---------- + event: + `f(params) -> 一批 bool`。吃到 `sampler` 的**整批**参数， + 返回长度为 `n` 的布尔数组（写成 `z > 0` 这种逐元素判断即可，不要写 Python 循环）。 + sampler: + `f(rng, size, **kwargs) -> 一批参数`，契约同 `monte_carlo`： + 批量大小由本函数以 `size=n` 传入，`sampler` 必须能接住 `size`。 + n: + 模拟次数，正整数。与 `event` 返回的长度必须一致。 + confidence: + 区间置信水平，必须落在开区间 `(0, 1)` 内。 + random_state: + 随机种子。`None` 时跟随 `core.set_seed` 的全局流，见 `core.random.get_rng`。 + Returns ------- - `Result`：`metrics={"probability", "ci_low", "ci_high", "std_error"}`。 + `Result`：`data` 是 `n` 个 0/1 结果组成的数组， + `metrics={"probability", "ci_low", "ci_high", "std_error"}`。 + 区间是 Wilson score 区间（`p` 贴近 0 / 1 时比 Wald 稳），并夹在 `[0, 1]` 内（比例量纲，见决策 42）。 + + Raises + ------ + TypeError: + `n` 不是整数，或 `random_state` 类型非法。 + ValueError: + `n` 不是正整数，或 `confidence` 不在 `(0, 1)` 内。 Notes ----- 概率估计的相对误差 ≈ `1/√(n·p)`。**小概率事件（p < 1e-3）用朴素 MC 是不可行的**， 需要重要性采样或方差缩减——遇到这类题要提前意识到，别硬跑。 + 本函数当前的口径：`probability <= 1e-3` 时**直接判 `success=False`**（只给提示，不给区间）。 Examples -------- >>> import numpy as np - >>> r = monte_carlo_probability(lambda v: v > 0, lambda rng, **k: rng.normal(0, 1), n=2000, random_state=0) + >>> draw = lambda rng, size, **k: rng.normal(0, 1, size) + >>> r = monte_carlo_probability(lambda z: z > 0, draw, n=2000, random_state=0) >>> abs(r.metrics["probability"] - 0.5) < 0.1 True """ - # TODO(L1) ★ S 档 - raise NotImplementedError + if not callable(sampler) or not callable(event): + raise TypeError("sampler 和 event 必须是可调用函数") + if isinstance(n, bool) or not isinstance(n, (int, np.integer)): + raise TypeError(f"n 必须是整数，你传入了 {type(n).__name__}") + if int(n) <= 0: + raise ValueError(f"n 必须是正整数，你传入了 {n}") + n = int(n) + if confidence <= 0.0 or confidence >= 1: + raise ValueError(f"confidence必须在0和1之间,而你输入了{confidence}") + + rng = get_rng(random_state) + call_params = {"n": n, "random_state": random_state, "confidence": confidence} + + try: + params = sampler(rng, size=n) + data = np.asarray(event(params), dtype=float).ravel() + if np.isnan(data).any(): + raise ValueError("模拟结果中包含 NaN，请检查 event 或 sampler 的传入参数") + if len(data) != n: + raise ValueError(f"event 返回的结果数量({len(data)})必须与模拟次数 n({n})一致") + p_value = float(np.mean(data)) + std_error = float(np.sqrt(p_value * (1.0 - p_value) / n)) + + z_score = norm.ppf(1 - (1 - confidence) / 2) + center = (p_value + z_score**2 / (2 * n)) / (1 + z_score ** 2 / n) + half_width = (z_score / (1 + z_score ** 2 / n)) * np.sqrt(p_value * (1 - p_value) / n + z_score ** 2 / (4 * n ** 2)) + ci_low = max(0.0, center - half_width) + ci_high = min(1.0, center + half_width) + + warning_message = "" + if p_value <= 1e-3 or p_value >= 1 - 1e-3: + warning_message = "小概率事件 P < 1e-3 不建议使用monte_carlo_probability" + + metrics = { + "probability": p_value, + "ci_low": float(ci_low), + "ci_high": float(ci_high), + "std_error": std_error, + } + return Result( + name="monte_carlo_probability", + params=call_params, + metrics=metrics, + data=data, + success=True, + message=( + f"蒙特卡洛模拟成功：已完成 {n} 次模拟" + + (f"；{warning_message}" if warning_message else "") + ), + ) + except Exception as e: + return Result( + name="monte_carlo_probability", + params=call_params, + metrics={}, + data=None, + success=False, + message=f"蒙特卡洛模拟失败：{type(e).__name__} - {e}", + ) def monte_carlo_integral( diff --git a/tests/test_solve.py b/tests/test_solve.py index 6ff8fc1..cb46a5a 100644 --- a/tests/test_solve.py +++ b/tests/test_solve.py @@ -8,6 +8,11 @@ import pytest _NOT_IMPLEMENTED = "未实现：见 docs/03-roadmap.md L0/L1" +def _normal_batch(rng, size, **kwargs): + """`monte_carlo` 的 sampler 契约：收下实现传来的 `size`，一次抽够。""" + return rng.normal(0, 1, size) + + @pytest.mark.skip(reason=_NOT_IMPLEMENTED) def test_monte_carlo_integral_converges_to_analytic_value(): """∫₀¹ x² dx = 1/3。n=2e5 时误差应 < 0.01。""" @@ -17,30 +22,205 @@ def test_monte_carlo_integral_converges_to_analytic_value(): assert abs(r.data - 1 / 3) < 0.01 -@pytest.mark.skip(reason=_NOT_IMPLEMENTED) def test_monte_carlo_probability_of_standard_normal(): """P(Z > 0) = 0.5，n=5000 时误差应 < 0.05。""" from mcm_toolkit.solve import monte_carlo_probability - r = monte_carlo_probability( - lambda z: z > 0, - lambda rng, **kw: rng.standard_normal(), - n=5000, - random_state=0, - ) + r = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=5000, random_state=0) + assert r.success assert abs(r.metrics["probability"] - 0.5) < 0.05 -@pytest.mark.skip(reason=_NOT_IMPLEMENTED) def test_monte_carlo_is_reproducible(): from mcm_toolkit.solve import monte_carlo_probability - args = (lambda z: z > 0, lambda rng, **kw: rng.standard_normal()) - a = monte_carlo_probability(*args, n=1000, random_state=9) - b = monte_carlo_probability(*args, n=1000, random_state=9) + a = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=1000, random_state=9) + b = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=1000, random_state=9) assert a.metrics["probability"] == b.metrics["probability"] +def test_monte_carlo_returns_n_results_and_documented_metrics(): + from mcm_toolkit.solve import monte_carlo + + r = monte_carlo(_normal_batch, lambda x: x**2, n=2000, random_state=0) + assert r.success + assert r.data.shape == (2000,) + assert set(r.metrics) == {"mean", "std", "median", "p05", "p95", "ci_low", "ci_high"} + + +def test_monte_carlo_calls_simulator_only_once(): + """`simulator` 只被调用一次 —— 这是"必须向量化"的可执行版本。""" + from mcm_toolkit.solve import monte_carlo + + counter = {"calls": 0} + + def simulator(params): + counter["calls"] += 1 + return params + + r = monte_carlo(_normal_batch, simulator, n=1000, random_state=0) + assert r.success + assert counter["calls"] == 1 + + +def test_monte_carlo_data_is_reproducible_with_same_seed(): + from mcm_toolkit.solve import monte_carlo + + a = monte_carlo(_normal_batch, lambda x: x, n=2000, random_state=11) + b = monte_carlo(_normal_batch, lambda x: x, n=2000, random_state=11) + assert np.array_equal(a.data, b.data) + + +def test_monte_carlo_follows_set_seed_when_random_state_omitted(): + from mcm_toolkit.core import set_seed + from mcm_toolkit.solve import monte_carlo + + set_seed(5) + a = monte_carlo(_normal_batch, lambda x: x, n=500) + set_seed(5) + b = monte_carlo(_normal_batch, lambda x: x, n=500) + assert np.array_equal(a.data, b.data) + + +def test_monte_carlo_estimates_known_expectation(): + """E[X²] = 1（X ~ N(0, 1)）；n = 2e5 时误差应 < 0.02。""" + from mcm_toolkit.solve import monte_carlo + + r = monte_carlo(_normal_batch, lambda x: x**2, n=200_000, random_state=0) + assert abs(r.metrics["mean"] - 1.0) < 0.02 + + +def test_monte_carlo_docstring_example_runs(): + """docstring 的 Examples 必须真的能跑通（此前那一版是坏的）。""" + from mcm_toolkit.solve import monte_carlo + + draw = lambda rng, size, **k: rng.normal(0, 1, size) # noqa: E731 + r = monte_carlo(draw, lambda x: x**2, n=10000, random_state=0) + assert r.success + assert abs(r.metrics["mean"] - 1.0) < 0.1 + + +def test_monte_carlo_accepts_numpy_integer_n(): + """§4.7：`np.int64` 不是 `int` 的子类，不能漏。""" + from mcm_toolkit.solve import monte_carlo + + r = monte_carlo(_normal_batch, lambda x: x, n=np.int64(100), random_state=0) + assert r.success + assert r.data.shape == (100,) + + +def test_monte_carlo_counts_boolean_outcomes(): + """布尔结果按 0/1 统计 —— 事件概率是 MC 最常见的用法。""" + from mcm_toolkit.solve import monte_carlo + + r = monte_carlo( + lambda rng, size, **k: rng.normal(0, 1, size) > 0, + lambda x: x, + n=5000, + random_state=0, + ) + assert r.success + assert r.data.dtype == float + assert abs(r.metrics["mean"] - 0.5) < 0.05 + + +def test_monte_carlo_rejects_nan_results(): + """结果里出现 NaN 就判失败，不静默给出一列 nan 指标。""" + from mcm_toolkit.solve import monte_carlo + + r = monte_carlo( + lambda rng, size, **k: np.r_[np.zeros(size - 1), np.nan], + lambda x: x, + n=10, + random_state=0, + ) + assert not r.success + assert "NaN" in r.message + + +def test_monte_carlo_rejects_bad_arguments(): + from mcm_toolkit.solve import monte_carlo + + with pytest.raises(TypeError): + monte_carlo(123, lambda x: x) # type: ignore[arg-type] # 就是要传错 + with pytest.raises(TypeError): + monte_carlo(_normal_batch, 123) # type: ignore[arg-type] # 就是要传错 + with pytest.raises(TypeError): + monte_carlo(_normal_batch, lambda x: x, n=True) + with pytest.raises(TypeError): + monte_carlo(_normal_batch, lambda x: x, n=10.0) # type: ignore[arg-type] # 就是要传错 + with pytest.raises(ValueError): + monte_carlo(_normal_batch, lambda x: x, n=0) + + + +def test_monte_carlo_probability_returns_documented_metrics(): + from mcm_toolkit.solve import monte_carlo_probability + + r = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=20000, random_state=0) + assert r.success + assert set(r.metrics) == {"probability", "ci_low", "ci_high", "std_error"} + assert r.data.shape == (20000,) + assert 0.0 <= r.metrics["ci_low"] <= r.metrics["probability"] + assert r.metrics["probability"] <= r.metrics["ci_high"] <= 1.0 + + +def test_monte_carlo_probability_std_error_matches_binomial_formula(): + """`std_error` 就是比例的标准误 √(p(1-p)/n)。""" + from mcm_toolkit.solve import monte_carlo_probability + + n = 20000 + r = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=n, random_state=0) + p = r.metrics["probability"] + assert np.isclose(r.metrics["std_error"], np.sqrt(p * (1 - p) / n)) + + +def test_monte_carlo_probability_ci_narrows_with_n(): + """区间宽度随 n 增大而收窄（O(1/√n)）。""" + from mcm_toolkit.solve import monte_carlo_probability + + def width(n): + r = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=n, random_state=0) + return r.metrics["ci_high"] - r.metrics["ci_low"] + + assert width(100_000) < width(1_000) + + +def test_monte_carlo_probability_higher_confidence_is_wider(): + from mcm_toolkit.solve import monte_carlo_probability + + def width(confidence): + r = monte_carlo_probability( + lambda z: z > 0, _normal_batch, n=5000, random_state=0, confidence=confidence + ) + return r.metrics["ci_high"] - r.metrics["ci_low"] + + assert width(0.99) > width(0.90) + + + +def test_monte_carlo_probability_ci_matches_statsmodels_wilson(): + """区间必须对得上独立实现（statsmodels 的 wilson）—— `^` 写成异或那次就靠它拦住。""" + from statsmodels.stats.proportion import proportion_confint + + from mcm_toolkit.solve import monte_carlo_probability + + n = 5000 + r = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=n, random_state=0) + count = int(r.data.sum()) + lo, hi = proportion_confint(count, n, alpha=0.05, method="wilson") + assert np.isclose(r.metrics["ci_low"], lo) + assert np.isclose(r.metrics["ci_high"], hi) + + +def test_monte_carlo_probability_rejects_bad_confidence(): + from mcm_toolkit.solve import monte_carlo_probability + + for bad in (0.0, 1.0, -0.5, 1.5): + with pytest.raises(ValueError): + monte_carlo_probability(lambda z: z > 0, _normal_batch, n=100, confidence=bad) + + @pytest.mark.skip(reason=_NOT_IMPLEMENTED) def test_markov_steady_state_two_state_chain(): """P=[[0.9,0.1],[0.2,0.8]] 的平稳分布是 [2/3, 1/3]。"""