feat(solve): 审查 markov_simulate:工程侧收口 + 9 条测试

用户实现的第一版 markov_simulate 算法本体经探针实测正确(终点频率 [0.657, 0.343]
对照平稳分布 [2/3, 1/3],同种子可复现),本次只动工程侧,未改任何算法行:

- ruff 8 条全清(I001 导入排序 + 5 处 W293 / 2 处 W291)
- docstring 从一行补全为 Parameters / Returns / Raises / Notes / Examples
- 文案:马尔科夫→马尔可夫(3 处)、initail→initial、半角逗号→全角、paths[:,0]→paths[:, 0]
- arr, kind = as_array(...) 的 kind 全程未用 → arr, _
- tests/test_solve.py 新增 9 条 markov_simulate 测试(此前 0 条)
- STATUS.md 新增决策 56(initial 语义待用户拍板,A/B/C 三选项)与 5 条待办
- CHANGELOG 与 docs/dev/personal-pitfalls.md 易错点 27(入口分流后两条分支口径矛盾)

校验:pytest 188 passed / 31 skipped(219 收集);ruff 0;pyright 0
This commit is contained in:
guzhujushi committed 2026-10-05 23:11:04 +08:00
1 parent 9368dd91cb
commit dd9bd638b8
5 files changed
+323 -23

No files matched your search

+124 -6
View File
@@ -13,7 +13,7 @@ from typing import Any
import numpy as np
from scipy.stats import norm
from mcm_toolkit.core import Result, get_rng
from mcm_toolkit.core import Result, as_array, get_rng
__all__ = [
"markov_simulate",
@@ -463,13 +463,131 @@ def random_walk(
def markov_simulate(
transition: Any,
initial: Any,
steps: int = 100,
n_paths: int = 1,
steps: int | np.integer = 100,
n_paths: int | np.integer = 1,
random_state: int | None = None,
) -> Any:
"""马尔可夫链模拟,返回状态序列。"""
# TODO(L1)
raise NotImplementedError
"""马尔可夫链模拟:从 `initial` 出发,按 `transition` 逐行走 `steps` 步。
Parameters
----------
transition:
行随机矩阵(每行和为 1)的二维方阵,状态数就是它的阶数 `s`。
initial:
初始状态。**标量整数**时所有路径都从同一状态出发;
否则按一维数组处理,长度必须等于状态数 `s`,再写进每条路径的第 0 列
—— 只有 `n_paths == s` 时两者才对得上,见 Notes。
steps:
每条路径的步数,正整数。`data` 的列数是 `steps + 1`(多出的第 0 列是初始状态)。
n_paths:
路径条数,正整数,默认 1 条。
random_state:
随机种子。`None` 时跟随 `core.set_seed` 的全局流,见 `core.random.get_rng`。
Returns
-------
`Result`:`data` 是形状 `(n_paths, steps + 1)` 的整数状态序列,第 0 列是初始状态;
`metrics={"mean_terminal", "std_terminal"}`(终点状态的均值与标准差,`std` 用 `ddof=0`)。
Raises
------
TypeError:
`transition` 不是二维数组;`initial` 数组长度不等于状态数;
`steps` / `n_paths` 不是整数。
ValueError:
`transition` 不是方阵,或某一行和不等于 1.0;`steps` / `n_paths` 不是正整数。
Notes
-----
抽样用逆变换法:`cumsum` 得到累积分布,再拿 `[0, 1)` 的均匀数去"落格"。
一次抽完 `n_paths` 个状态,路径之间相互独立;步与步之间只能串行,因为下一步依赖当前状态。
`cum_probs[:, -1] = 1.0` 是为了压掉浮点累加误差,避免随机数落在 `[1 - 1e-16, 1)` 时索引越界。
终点分布随步数收敛到平稳分布:这既是本函数的自检手段,也是论文里值得写一句的性质。
遗留(决策 56 待定):`initial` 现在靠 `isinstance(initial, int)` 区分"标量状态 / 状态数组",
`np.int64`、`bool` 会掉进数组分支,"按状态数校验长度"与"按路径数写入"也互相矛盾。
口径定下来后,上面 `initial` 这一节与本段随之改写。
Examples
--------
>>> import numpy as np
>>> P = np.array([[0.9, 0.1], [0.2, 0.8]])
>>> r = markov_simulate(P, initial=0, steps=50, n_paths=4, random_state=0)
>>> r.data.shape
(4, 51)
>>> bool(np.all(r.data[:, 0] == 0))
True
"""
arr, _ = as_array(transition)
if arr.ndim != 2 or arr.shape[0] != arr.shape[1]:
raise ValueError(f"transition 必须是二维方阵,你传入的形状是 {arr.shape}")
if not np.allclose(arr.sum(axis=1), 1.0):
raise ValueError("transition 每一行的和必须是 1.0")
s = arr.shape[0]
if isinstance(n_paths, bool) or not isinstance(n_paths, (int, np.integer)):
raise TypeError(f"n_paths 必须是整数,你传入了 {type(n_paths).__name__}")
if n_paths <= 0:
raise ValueError(f"n_paths 必须是正整数,你传入了 {n_paths}")
if isinstance(steps, bool) or not isinstance(steps, (int, np.integer)):
raise TypeError(f"steps 必须是整数,你传入了 {type(steps).__name__}")
if steps <= 0:
raise ValueError(f"steps 必须是正整数,你传入了 {steps}")
rng = get_rng(random_state)
if isinstance(initial, int):
current_states = np.full(n_paths, initial)
else:
current_states, _ = as_array(initial, ndim=1)
if len(current_states) != s:
raise TypeError(f"initial 的长度必须是 {s},而你传入的是 {len(current_states)}")
paths = np.zeros((n_paths, steps + 1), dtype=int)
paths[:, 0] = current_states
call_params = {
"n_steps": steps,
"n_paths": n_paths,
"initial": initial,
"random_state": random_state,
}
try:
for t in range(steps):
probs = arr[current_states]
cum_probs = np.cumsum(probs, axis=1)
cum_probs[:, -1] = 1.0
r = rng.random(n_paths)
next_states = (cum_probs < r[:, None]).sum(axis=1)
paths[:, t + 1] = next_states
current_states = next_states
terminal_positions = paths[:, -1]
metrics = {
"mean_terminal": float(np.mean(terminal_positions)),
"std_terminal": float(np.std(terminal_positions)),
}
return Result(
name="markov_simulate",
params=call_params,
metrics=metrics,
data=paths,
success=True,
message=f"马尔可夫链模拟成功,生成了 {n_paths} 条路径,每条路径 {steps} 步",
)
except Exception as e:
return Result(
name="markov_simulate",
params=call_params,
metrics={},
data=None,
success=False,
message=f"马尔可夫链模拟失败:{type(e).__name__} - {str(e)}",
)
def markov_steady_state(transition: Any, max_iter: int = 10000, tol: float = 1e-10) -> Any: