feat(solve): monte_carlo 系列两个函数落地（L0 11/12、L1）+ Wilson 区间；修 z_score^2 异或 bug
用户实现： - 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。
This commit is contained in:
guzhujushi committed
1 parent
58f2236703
commit
b0c4b4a413
5 files changed
+619
-61
No files matched your search
@@ -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(
|
||||
|
||||
Reference in new issue
Block a user