feat(solve): 审查 markov_steady_state:工程侧收口 + 8 条测试
用户实现的 markov_steady_state 经探针实测算法本体正确(3 状态与特征分解逐位一致、 200x200 一致到 1.9e-12 且只用 6 次迭代 / 0.10s),本次只动工程侧,未改任何算法行: - ruff 5 条 W293 全清;mm1_metrics 前 4 个空行 -> 2;失败文案改成全角冒号 - docstring 补 max_iter / tol 两个参数与 Raises,写清 tol 是相邻差口径 (实测慢混合链 tol=1e-10 -> 实际误差 4.6e-9,被 1/(1-|lambda2|) 放大约 50 倍) - Notes 修正"周期链会不收敛"的说法:可约链的实情是平稳分布不唯一, 迭代法给的是离初值最近的那个(均匀初值对双随机矩阵本身就是不动点) - tests/test_solve.py:摘掉 1 条 skip + 新增 7 条(该函数共 8 条实跑) - STATUS.md 新增决策 57(6 条待拍板口径)与第六轮「下一步」;CHANGELOG;易错点 28 校验:pytest 196 passed / 30 skipped(226 收集);ruff 0;pyright 0
This commit is contained in:
1 parent
dd9bd638b8
commit
0955890010
5 files changed
+263
-27
No files matched your search
@@ -529,6 +529,32 @@
|
||||
转移矩阵未查非负(行和恰为 1.0 的 `[[1.5, -0.5], [0, 1]]` 会静默给出全 0)、`paths[:, 0]` 赋值在 `try` 之外
|
||||
(长度不匹配时裸抛 `ValueError`)、状态越界没在入口拦、`std_terminal` 对分类状态的意义。
|
||||
|
||||
### Added
|
||||
- `solve.markov_steady_state` 落地(**L1,S 档**,实现由用户本人写):`v` 从均匀分布出发,`v_next = v @ arr`
|
||||
反复推进(行随机矩阵按行向量右乘),直到 `max|v_next - v| < tol` 或耗尽 `max_iter`;
|
||||
`data` 是平稳分布向量、`metrics={"converged", "iterations"}`。Codex 审查结论:**算法本体正确** ——
|
||||
3 状态与 `np.linalg.eig(P.T)` 特征值 1 的特征向量逐位一致、200×200 一致到 `1.9e-12`(6 次迭代 / 0.10s)、
|
||||
与 `markov_simulate` 的终点频率互验 `[0.6702, 0.3298]` vs `[0.6667, 0.3333]`。
|
||||
- `tests/test_solve.py`:摘掉 `test_markov_steady_state_two_state_chain` 的 skip,另新增 7 条 ——
|
||||
与特征分解对拍、`data` 是分布(和为 1 / 非负 / `float64`)、慢混合链迭代次数 > 100、`np.int64` 的 `max_iter`、
|
||||
非法转移矩阵、`max_iter=0` / `tol=0`、docstring Example 真能跑。
|
||||
|
||||
### Changed
|
||||
- `markov_steady_state` 的 docstring 补全:`max_iter` / `tol` 两个参数、`Raises` 一节,
|
||||
并写清 `tol` 是**相邻差口径**(实测慢混合链 `tol=1e-10` 的实际误差 `4.6e-9`,被 `1/(1-|λ2|)` 放大约 50 倍)。
|
||||
- `markov_steady_state` 的 Notes 修正"周期链会不收敛"的说法:可约链的实情是**平稳分布不唯一**,
|
||||
迭代法给出的是"离初值最近的那个"(均匀初值对双随机矩阵本身就是不动点,所以 `iterations=1` 不代表链唯一平稳)。
|
||||
- `markov_steady_state` 的工程侧清理:ruff 5 条 `W293`;`mm1_metrics` 前 4 个空行 → 2;失败文案 `失败: ` → `失败:`。
|
||||
- `STATUS.md`:新增**决策 57**(6 条待拍板口径)与「下一步」第六轮;验证计数更新为 **196 passed / 30 skipped(226 收集)**。
|
||||
- `docs/dev/personal-pitfalls.md`:新增易错点 28(数值阈值参数的 `bool` / NaN 漏网)。
|
||||
|
||||
### Notes
|
||||
- **`markov_steady_state` 的 6 条口径留给用户**(决策 49):转移矩阵未查非负(行和恰为 1.0 的 `[[1.5, -0.5], [0, 1]]`
|
||||
会溢出成 `nan`,而 `data=[nan, nan]` 却 `success=True`,违反决策 48 的 NaN 口径;附带的 `RuntimeWarning`
|
||||
在 pyproject 里被设成 error,**测试会直接报错**)、`tol=True` 被当 `tol=1` 静默放行(`isinstance(True, int)`)、
|
||||
`tol=np.nan` 穿过校验白跑迭代、`converged=False` 仍 `success=True`、
|
||||
类型错也抛 `ValueError`(与决策 54 分工不一致)、`params` 没记转移矩阵。
|
||||
|
||||
## [0.0.1] - 2026-10-01
|
||||
|
||||
### Added
|
||||
|
||||
@@ -21,10 +21,10 @@
|
||||
当前 **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`、`monte_carlo_integral`、`random_walk`、`markov_simulate`
|
||||
(**均为用户写**,2026-10-05)。剩余 1 项:`save_figure`(Codex 写)。
|
||||
`monte_carlo_probability`、`monte_carlo_integral`、`random_walk`、`markov_simulate`、
|
||||
`markov_steady_state`(**均为用户写**)。剩余 1 项:`save_figure`(Codex 写)。
|
||||
`monte_carlo` 的口径已拍板(决策 48);`monte_carlo_probability` / `random_walk` 的待办**均已办结**;
|
||||
`markov_simulate` 本轮审查通过(算法本体正确),留 **5 条入口口径待你拍板**(见"下一步",决策 56 的 3 个选项)。
|
||||
`markov_simulate`(决策 56,5 条)与 `markov_steady_state`(决策 57,6 条)的**入口口径待你拍板**。
|
||||
另 3 个小口径仍悬着:`ci_low`/`ci_high` 的语义、`monte_carlo_integral` 的第一道 `isinstance` 分支可塌成两行、
|
||||
`random_state` 注解与 `std` 的 `ddof` 是否全项目对齐。
|
||||
|
||||
@@ -335,13 +335,29 @@
|
||||
`np.int64` 可收、docstring Example 真能跑
|
||||
- ⏳ **留给你的 5 条**见"下一步 → markov_simulate 待你拍板"(全是入口口径与少考虑情况,决策 49)
|
||||
|
||||
- [x] `solve/simulate.py` 的 `markov_steady_state`(**L1,S 档,算法由用户本人写**,2026-10-09)
|
||||
- 核心:`v` 从均匀分布出发,`v_next = v @ arr` 反复推进(行随机矩阵按**行向量右乘**),
|
||||
直到 `max|v_next - v| < tol` 或耗尽 `max_iter`;`data` 是平稳分布向量,`metrics={"converged", "iterations"}`
|
||||
- **Codex 审查结论:算法本体正确**(探针实测)—— 2 状态 `[2/3, 1/3]`(58 次迭代)、3 状态与 `np.linalg.eig(P.T)`
|
||||
特征值 1 对应的特征向量逐位一致(39 次)、随机 6 状态一致(13 次)、200×200 一致到 `1.9e-12`(6 次、0.10s)、
|
||||
与 `markov_simulate` 的终点频率互验 `[0.6702, 0.3298]` vs `[0.6667, 0.3333]`、`sum(data) = 1.0`
|
||||
- **Codex 只碰工程侧**:① ruff 5 条 `W293` 全清,`mm1_metrics` 前 4 个空行 → 2;
|
||||
② docstring 补 `max_iter` / `tol` 两个参数(原来 Parameters 里只有 `transition`)、补 `Raises`,
|
||||
写清 `tol` 是**相邻差口径**(实测慢混合链 `tol=1e-10` → 实际误差 `4.6e-9`,被 `1/(1-|λ2|)` 放大);
|
||||
③ 修正 Notes 里"周期链会不收敛"的说法(可约链的实情是"平稳分布不唯一、迭代法给离初值最近的那个");
|
||||
④ 失败文案 `失败: ` → `失败:`(与文件内其它函数统一)
|
||||
- **测试:+8 条、摘掉 1 条 skip**(此前该函数只有 1 条挂着的 skip)—— 2 状态 `[2/3, 1/3]` + `converged` / `iterations`、
|
||||
3 状态与特征分解对拍、`data` 是分布(和为 1 / 非负 / `float64`)、慢混合链迭代次数 > 100、
|
||||
`np.int64` 的 `max_iter`、非法转移矩阵、`max_iter=0` / `tol=0`、docstring Example 真能跑
|
||||
- ⏳ **留给你的 6 条口径**见"下一步 → markov_steady_state 待你拍板"(决策 57)
|
||||
|
||||
- [ ] 剩余 1 项:`save_figure`(A 档 = 工程,**Codex 直接写**)
|
||||
|
||||
### 验证结果
|
||||
|
||||
| 检查 | venv | conda `mcm-toolkit` |
|
||||
|---|---|---|
|
||||
| `pytest` | 见下(单列) | **188 passed, 31 skipped(0 xfail)** |
|
||||
| `pytest` | 见下(单列) | **196 passed, 30 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` | 跑通 | — |
|
||||
@@ -351,15 +367,15 @@
|
||||
(修配置之前是 17 errors,全是"无法解析导入 numpy / pandas",见决策 29 与
|
||||
`docs/04-environment.md` 坑 7)
|
||||
- base 亦已复核:numpy 2.4.6 / scipy 1.18.1 / pandas 3.0.6(与两套环境一致)
|
||||
- **219 条用例:188 passed / 31 skipped / 0 xfail**。收集数按文件:`test_data.py` 52、`test_api_surface.py` 20、
|
||||
`test_core.py` 23、`test_solve.py` 51(本轮 +9:`markov_simulate` 从 0 条起步)、`test_viz.py` 8、`test_decision.py` 7、
|
||||
- **226 条用例:196 passed / 30 skipped / 0 xfail**。收集数按文件:`test_data.py` 52、`test_api_surface.py` 20、
|
||||
`test_core.py` 23、`test_solve.py` 58(本轮 +8、摘掉 1 条 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`
|
||||
- 31 条 skip 是等实现的算法(`train_test_split`、`markov_steady_state`、`mm1_metrics`、`anova`、`fit_distribution`、
|
||||
- 30 条 skip 是等实现的算法(`train_test_split`、`mm1_metrics`、`anova`、`fit_distribution`、
|
||||
`knapsack_dp`、`linear_program`、`simulated_annealing` 等);
|
||||
`monte_carlo_probability` / `monte_carlo_integral` 已落地,分别摘掉 2 条与 1 条 skip;`markov_simulate` 新增的 9 条全部实跑;
|
||||
`monte_carlo_probability` / `monte_carlo_integral` / `markov_steady_state` 已落地,分别摘掉 2 / 1 / 1 条 skip;`markov_simulate` 新增的 9 条全部实跑;
|
||||
`TestValidate` 的 4 条仍挂着 —— `check_*` 还是 `NotImplementedError`,但它们实现时应当直接复用 `as_array`
|
||||
- **本轮回归(已修)**:`transform.py` 上一版有 3 处真 bug —— 修前 **7 failed / 75 passed**,
|
||||
且 pyright 在同一处报 `reportArgumentType`。**不是静态误报**:pyright CLI 复现 + 测试真红。
|
||||
@@ -527,17 +543,29 @@ docstring 与 `test_describe_reports_expected_keys` 早就写死了这个方向
|
||||
`True` 因 `isinstance(True, int)` 为真被当状态 1 **静默放行**;越界状态(如 `initial=5`、`s=2`)只在
|
||||
第一轮索引时才炸、被 `except` 吞成 `success=False`。可选口径见下一步的 A / B / C。
|
||||
|
||||
| 57 | `markov_steady_state` 的 6 条口径**待用户拍板**(Codex 只报选项,未代改,决策 49) |
|
||||
探针实测出的口子:① **含负元素**的"行和为 1"矩阵(`[[1.5, -0.5], [0, 1]]`)会溢出成 `nan`,而
|
||||
`data=[nan, nan]` 却 `success=True` —— 违反决策 48 的 NaN 口径(根因:没查 `arr >= 0`);附带的
|
||||
`RuntimeWarning: overflow encountered in matmul` 因为 pyproject 设了 `filterwarnings = ["error::RuntimeWarning"]`,
|
||||
**能让测试直接报错**,所以不能只靠测试兜;② `tol=True` 因 `isinstance(True, int)` 通过校验 → `tol=1` →
|
||||
1 次就"收敛",给出 `[0.55, 0.45]` 的假答案(`max_iter` 查了 `bool`、`tol` 没查,易错点 1 复发);
|
||||
③ `tol=np.nan` 因 `nan <= 0` 为假穿过校验 → 白跑 10000 次迭代;④ `converged=False` 仍 `success=True`
|
||||
且 `data` 是未收敛的迭代值(`max_iter=5` → `[0.6387, 0.3613]`);⑤ 类型错也抛 `ValueError`
|
||||
(`max_iter=100.0`),与决策 54 的"类型错 `TypeError`"分工不一致;⑥ `params` 只记 `max_iter` / `tol`,
|
||||
没记转移矩阵(`monte_carlo_integral` 记了 `bounds`),事后复现不了。可选口径见"下一步"。
|
||||
|
||||
## 下一步
|
||||
|
||||
### ⭐ 当前进度与待办(2026-10-05 第五轮)
|
||||
### ⭐ 当前进度与待办(2026-10-09 第六轮)
|
||||
|
||||
**本轮:`markov_simulate`(你写的第一版)已审查通过** —— 算法本体正确,没动你任何算法行;
|
||||
Codex 只清工程侧(ruff 8 条、docstring、3 处 `马尔科夫`、`initail` 拼写、半角逗号、`kind` 未用),
|
||||
并补 **9 条测试**(此前 0 条)。三连:`pytest` **188 passed / 31 skipped**、`ruff` 0、`pyright` 0。
|
||||
**本轮:`markov_steady_state`(你写的)已审查通过** —— 算法本体正确,一行没动(3 状态与特征分解逐位一致、
|
||||
200×200 只用 6 次迭代 / 0.10s);Codex 只清工程侧(ruff 5 条、docstring 补 `max_iter`/`tol`/`Raises`、
|
||||
修正"周期链"表述、全角冒号),并 **+8 条测试、摘掉 1 条 skip**。
|
||||
三连:`pytest` **196 passed / 30 skipped**、`ruff` 0、`pyright` 0。
|
||||
|
||||
**接下来三件事**
|
||||
1. **`markov_simulate` 的 5 条待你拍板**(下节,全是入口口径,决策 49/56)。
|
||||
2. `ci_low` / `ci_high` 的语义确认一句(A / B / C 三个选项在下面挂着),确认后 Codex 补 `docs/algorithms/monte_carlo.md` 卡片。
|
||||
1. **`markov_steady_state` 的 6 条口径待你拍板**(下节,决策 57)。
|
||||
2. **`markov_simulate` 的 5 条口径**仍在挂着(决策 56 的 A / B / C),与 `solve` 的收尾一起定。
|
||||
3. L0 最后一项 `save_figure`(A 档 = 工程,`viz/export.py`)—— **Codex 直接实现**,你只读。
|
||||
|
||||
### ⏳ markov_simulate 待你拍板(5 条,Codex 只报不改)
|
||||
@@ -556,9 +584,31 @@ Codex 只清工程侧(ruff 8 条、docstring、3 处 `马尔科夫`、`initail
|
||||
**⑤ `std_terminal` 对**分类**状态的意义**:状态编号是标称量,均值和标准差只对"有序状态"(如排队人数)成立 ——
|
||||
要不要把 `metrics` 换成状态频率(`terminal_freq`)?你定。
|
||||
|
||||
**你下一个要写的(S 档,仍在 `solve/simulate.py`)**:`markov_steady_state`。
|
||||
与 `markov_simulate` 是同一件事的两面:前者按转移矩阵逐代抽状态,后者用**迭代**求平稳分布。
|
||||
docstring 已写明口径:用迭代而不用特征分解,因为"收敛快不快"本身就是论文里要写的性质。
|
||||
### ⏳ markov_steady_state 待你拍板(6 条,Codex 只报不改)
|
||||
|
||||
**① 转移矩阵没查非负**(与 `markov_simulate` 的待办 ② 是同一块):`[[1.5, -0.5], [0, 1]]` 行和恰为 1.0 →
|
||||
通过校验 → `v` 指数发散 → `data=[nan, nan]` **但 `success=True`**(违反决策 48 的 NaN 口径)。
|
||||
实测还会先抛 `RuntimeWarning: overflow encountered in matmul`,而 pyproject 设了
|
||||
`filterwarnings = ["error::RuntimeWarning"]` —— 这条**能让测试直接报错**,所以得在校验层拦,不能只靠测试兜。
|
||||
|
||||
**② `tol` 没查 `bool`**:`tol=True` → `isinstance(True, int)` 为真 → `tol=1` →
|
||||
1 次就"收敛",给出 `[0.55, 0.45]` 的**假答案**(`success=True, converged=True, iterations=1`)。
|
||||
`max_iter` 查了 `bool`,`tol` 漏了(易错点 1 复发)。
|
||||
|
||||
**③ `tol` 没查 NaN**:`np.nan <= 0` 为假 → 穿过校验 → 白跑 10000 次迭代后返回(结果本身是对的,纯浪费)。
|
||||
|
||||
**④ `converged=False` 仍 `success=True`**,且 `data` 是未收敛的迭代值(`max_iter=5` → `[0.6387, 0.3613]`)。
|
||||
选项:**A** 保持现状(`converged` 自己看);**B** 不收敛判 `success=False`;**C** 不收敛时 `data=None`、只留 `metrics`。
|
||||
|
||||
**⑤ 类型错也抛 `ValueError`**(`max_iter=100.0`),与文件内其它函数(决策 54:类型错 `TypeError`、取值错 `ValueError`)
|
||||
不一致 —— 要不要对齐?`tol` 同理。
|
||||
|
||||
**⑥ `params` 没记转移矩阵**:`{"max_iter": 10000, "tol": 1e-10}` 复现不了"算的是哪个矩阵"
|
||||
(`monte_carlo_integral` 记了 `bounds`)—— 要不要补?
|
||||
|
||||
**你下一个要写的(S 档,仍在 `solve/simulate.py`)**:`mm1_metrics`(M/M/1 闭式指标,10 行)。
|
||||
docstring 已写明口径:`λ >= μ` 时系统不稳定、队列长度无界 —— 这种情况返回 `success=False` 并说明原因,
|
||||
而不是返回 `inf`。写完 `solve` 这一层就齐了,L1 只剩 `decision` / `models` 那几条。
|
||||
|
||||
### ✅ random_walk 的 4 条待办(2026-10-05 已办结)
|
||||
|
||||
|
||||
@@ -468,3 +468,22 @@
|
||||
对不上就说明口径没定清,**先定语义再写代码**;② 分流用的 `isinstance` 必须把 `np.integer` / `bool` 一起想到(易错点 1);
|
||||
③ 任何"写进结果数组"的赋值,问一遍"它的长度约束和上面校验用的是同一个下标维度吗"。
|
||||
- 首次:2026-10-05|复发:0
|
||||
|
||||
### 28. 数值阈值参数的两个漏网:`bool` 和 NaN
|
||||
|
||||
- 现象:`markov_steady_state` 的 `tol` 只写了 `if not isinstance(tol, (int, float)) or tol <= 0`。
|
||||
`tol=True` → `isinstance(True, int)` 为真、`True > 0` 也为真 → 阈值变成 `1`,
|
||||
**第 1 次迭代就"收敛"**,返回 `[0.55, 0.45]` 这种假答案,而且 `success=True, converged=True`。
|
||||
`tol=np.nan` → `nan <= 0` 为假 → 穿过校验 → 白跑满 `max_iter` 次迭代。
|
||||
同一个函数里 `max_iter` 查了 `bool`、`tol` 没查 —— 说明是**逐条手动补的**,不是照着清单过的。
|
||||
- 为什么会错:脑子里"阈值 = 一个正数",于是校验只写"是数 & 大于 0";而 `bool` 是 `int` 的子类、
|
||||
NaN 与任何数比较都是 `False`,这两条都恰好从"是数 & 大于 0"里漏出去。**与易错点 1 同根**,
|
||||
只是这次漏的是"比较"而不是"类型判断"。
|
||||
- 同族(同一轮发现的):只校验"转移矩阵行和 = 1.0"就以为矩阵合法 —— `[[1.5, -0.5], [0, 1]]` 行和也是 1.0,
|
||||
但含负元素,迭代时 `v` 指数发散、最后 `data=[nan, nan]` 却 `success=True`。
|
||||
**"行和 = 1"不等于"是概率矩阵"**,还要 `arr >= 0`。另外这类溢出会抛 `RuntimeWarning`,
|
||||
而本项目 `filterwarnings = ["error::RuntimeWarning"]` —— **测试会直接报错**,别指望测试能兜住静默错。
|
||||
- 以后怎么防:① 校验数值阈值时,脑子里过一遍"三件套":**类型(含 `bool` / `np.integer` 例外)、
|
||||
有限性(`np.isfinite`,挡 NaN / inf)、取值区间**;② 同一函数里同类参数(`max_iter` / `tol` / `n` / `tol`)
|
||||
用**同一套模板**写,别一个查 `bool` 一个不查;③ 校验"概率/权重"这类输入时,除了和为 1,还要问一句"有没有可能是负数/NaN"。
|
||||
- 首次:2026-10-09|复发:0
|
||||
@@ -590,23 +590,42 @@ def markov_simulate(
|
||||
)
|
||||
|
||||
|
||||
def markov_steady_state(transition: Any, max_iter: int = 10000, tol: float = 1e-10) -> Any:
|
||||
"""平稳分布。
|
||||
def markov_steady_state(transition: Any, max_iter: int | np.integer = 10000, tol: float = 1e-10) -> Any:
|
||||
"""迭代法求马尔可夫链的平稳分布(`v @ P` 反复推进,直到相邻两次几乎不变)。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
transition:
|
||||
行随机矩阵(每行和为 1)。
|
||||
行随机矩阵(每行和为 1)的二维方阵。
|
||||
max_iter:
|
||||
迭代上限,正整数,默认 10000。
|
||||
tol:
|
||||
收敛阈值,正数,默认 1e-10。判据是**相邻两次迭代的最大分量差**
|
||||
`max|v_k+1 - v_k| < tol`,不是"到不动点的距离",见 Notes。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`:`data` 是平稳分布向量,`metrics={"converged", "iterations"}`。
|
||||
`Result`:`data` 是一维平稳分布向量(和为 1);
|
||||
`metrics={"converged", "iterations"}`(是否在 `max_iter` 内向容、实际迭代次数)。
|
||||
|
||||
Raises
|
||||
------
|
||||
ValueError:
|
||||
`transition` 不是二维方阵、或某一行和不等于 1.0;
|
||||
`max_iter` 不是正整数;`tol` 不是正数。
|
||||
|
||||
Notes
|
||||
-----
|
||||
用迭代法而不是特征分解:迭代法能顺便告诉你"收敛快不快",
|
||||
而收敛速度本身就是论文里值得写一句的性质。
|
||||
**不可约 + 非周期**才有唯一平稳分布,周期链会不收敛,要检查。
|
||||
初值是**均匀分布** `1/s`;**不可约 + 非周期**才有唯一平稳分布。
|
||||
可约链(如单位阵 `eye(2)`)的平稳分布不唯一,迭代法给出的是"离初值最近的那个"
|
||||
—— 均匀初值对双随机矩阵本身就是不动点,所以 `converged=True, iterations=1` 只说明"初值恰是不动点",
|
||||
不代表这条链有唯一平稳分布;周期链从一般初值出发会一直振荡、`converged=False`。
|
||||
|
||||
`tol` 是**相邻差**口径,实际误差约被 `1 / (1 - |λ2|)` 放大:实测双状态慢混合链
|
||||
`[[0.999, 0.001], [0.02, 0.98]]`(`|λ2| = 0.979`)用默认 `tol=1e-10` 停在第 867 次,
|
||||
与特征向量相差 `4.6e-9`(约 50 倍)。慢混合链要更准,把 `tol` 压小即可。
|
||||
|
||||
Examples
|
||||
--------
|
||||
@@ -616,8 +635,57 @@ def markov_steady_state(transition: Any, max_iter: int = 10000, tol: float = 1e-
|
||||
>>> bool(np.allclose(r.data, [2 / 3, 1 / 3], atol=1e-6))
|
||||
True
|
||||
"""
|
||||
# TODO(L1)
|
||||
raise NotImplementedError
|
||||
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]
|
||||
v = np.full(s, 1.0 / s)
|
||||
|
||||
iterations = 0
|
||||
converged = False
|
||||
|
||||
if isinstance(max_iter, bool) or not isinstance(max_iter, (int, np.integer)) or max_iter <= 0:
|
||||
raise ValueError(f"max_iter 必须是正整数,你传入了 {max_iter}")
|
||||
if not isinstance(tol, (int, float)) or tol <= 0:
|
||||
raise ValueError(f"tol 必须是正数,你传入了 {tol}")
|
||||
|
||||
try:
|
||||
for _ in range(max_iter):
|
||||
v_next = v @ arr
|
||||
error = np.max(np.abs(v_next - v))
|
||||
v = v_next
|
||||
iterations += 1
|
||||
|
||||
if error < tol:
|
||||
converged = True
|
||||
break
|
||||
|
||||
metrics = {
|
||||
"converged": converged,
|
||||
"iterations": iterations
|
||||
}
|
||||
|
||||
return Result(
|
||||
name="markov_steady_state",
|
||||
params={"max_iter": max_iter, "tol": tol},
|
||||
metrics=metrics,
|
||||
data=v,
|
||||
success=True,
|
||||
message=f"平稳分布求解{'成功' if converged else '未完全收敛'},迭代 {iterations} 次。"
|
||||
)
|
||||
|
||||
except Exception as e:
|
||||
return Result(
|
||||
name="markov_steady_state",
|
||||
params={"max_iter": max_iter, "tol": tol},
|
||||
metrics={},
|
||||
data=None,
|
||||
success=False,
|
||||
message=f"平稳分布求解失败:{type(e).__name__} - {str(e)}"
|
||||
)
|
||||
|
||||
|
||||
def mm1_metrics(arrival_rate: float, service_rate: float) -> Any:
|
||||
|
||||
+76
-3
@@ -508,15 +508,88 @@ def test_markov_simulate_docstring_example_runs():
|
||||
assert np.all(r.data[:, 0] == 0)
|
||||
|
||||
|
||||
|
||||
|
||||
@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]。"""
|
||||
from mcm_toolkit.solve import markov_steady_state
|
||||
|
||||
r = markov_steady_state(np.array([[0.9, 0.1], [0.2, 0.8]]))
|
||||
assert r.success
|
||||
assert np.allclose(r.data, [2 / 3, 1 / 3], atol=1e-6)
|
||||
assert r.metrics["converged"] is True
|
||||
assert r.metrics["iterations"] > 0
|
||||
|
||||
|
||||
def test_markov_steady_state_matches_eigenvector():
|
||||
"""与特征分解对拍:取特征值 1 对应的左特征向量(只是检验手段,实现仍走迭代)。"""
|
||||
from mcm_toolkit.solve import markov_steady_state
|
||||
|
||||
P = np.array([[0.5, 0.25, 0.25], [0.1, 0.8, 0.1], [0.2, 0.2, 0.6]])
|
||||
eigenvalues, eigenvectors = np.linalg.eig(P.T)
|
||||
ref = np.real(eigenvectors[:, int(np.argmin(np.abs(eigenvalues - 1.0)))])
|
||||
ref = ref / ref.sum()
|
||||
|
||||
r = markov_steady_state(P)
|
||||
assert np.allclose(r.data, ref, atol=1e-6)
|
||||
|
||||
|
||||
def test_markov_steady_state_is_a_distribution():
|
||||
from mcm_toolkit.solve import markov_steady_state
|
||||
|
||||
P = np.array([[0.7, 0.2, 0.1], [0.1, 0.6, 0.3], [0.05, 0.15, 0.8]])
|
||||
r = markov_steady_state(P)
|
||||
assert r.data.shape == (3,)
|
||||
assert np.isclose(float(np.sum(r.data)), 1.0)
|
||||
assert np.all(r.data >= 0.0)
|
||||
assert r.data.dtype == np.float64
|
||||
|
||||
|
||||
def test_markov_steady_state_slow_mixing_needs_more_iterations():
|
||||
"""慢混合链:迭代次数显著变多,但结果仍要对得上特征向量(tol 是相邻差口径)。"""
|
||||
from mcm_toolkit.solve import markov_steady_state
|
||||
|
||||
P = np.array([[0.999, 0.001], [0.02, 0.98]])
|
||||
r = markov_steady_state(P)
|
||||
assert r.metrics["iterations"] > 100
|
||||
assert np.allclose(r.data, [0.9523809524, 0.0476190476], atol=1e-6)
|
||||
|
||||
|
||||
def test_markov_steady_state_accepts_numpy_integer():
|
||||
"""§4.7:`np.int64` 不是 `int` 的子类,注解与运行期都要收(易错点 1)。"""
|
||||
from mcm_toolkit.solve import markov_steady_state
|
||||
|
||||
P = np.array([[0.9, 0.1], [0.2, 0.8]])
|
||||
r = markov_steady_state(P, max_iter=np.int64(10000))
|
||||
assert r.success
|
||||
assert r.metrics["converged"] is True
|
||||
|
||||
|
||||
def test_markov_steady_state_rejects_bad_transition():
|
||||
from mcm_toolkit.solve import markov_steady_state
|
||||
|
||||
with pytest.raises(ValueError):
|
||||
markov_steady_state(np.ones((2, 3)) / 3)
|
||||
with pytest.raises(ValueError):
|
||||
markov_steady_state(np.array([[0.5, 0.4], [0.2, 0.8]]))
|
||||
|
||||
|
||||
def test_markov_steady_state_rejects_bad_arguments():
|
||||
"""取值错一律 `ValueError`(类型错要不要拆 `TypeError` 待定,见 STATUS 决策 57)。"""
|
||||
from mcm_toolkit.solve import markov_steady_state
|
||||
|
||||
P = np.array([[0.9, 0.1], [0.2, 0.8]])
|
||||
with pytest.raises(ValueError):
|
||||
markov_steady_state(P, max_iter=0)
|
||||
with pytest.raises(ValueError):
|
||||
markov_steady_state(P, tol=0)
|
||||
|
||||
|
||||
def test_markov_steady_state_docstring_example_runs():
|
||||
from mcm_toolkit.solve import markov_steady_state
|
||||
|
||||
P = np.array([[0.9, 0.1], [0.2, 0.8]])
|
||||
r = markov_steady_state(P)
|
||||
assert bool(np.allclose(r.data, [2 / 3, 1 / 3], atol=1e-6))
|
||||
|
||||
|
||||
|
||||
@pytest.mark.skip(reason=_NOT_IMPLEMENTED)
|
||||
|
||||
Reference in new issue
Block a user