diff --git a/CHANGELOG.md b/CHANGELOG.md index 80fb8a6..58bb5c4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -46,6 +46,13 @@ 单列 DataFrame、均值为 0 抛 `ValueError`、空数组与二维输入。 - `examples/01_clean_and_describe.py`：填掉末尾的 `TODO(L0)`，加一段 `describe` 输出（实测 exit 0）。 +### Added +- `stats.correlation_pvalues` / `stats.correlation_table`：`correlation_matrix` 的伴随函数 + （**方案 A**，用户 2026-10-03 拍板，见决策 41）—— + `correlation_matrix` 只返回**数值**矩阵；`correlation_pvalues` 给数值 p 值矩阵、 + `correlation_table` 给"系数 + 显著性星标"的论文表格。 + 三者共用 `_check_corr_args` 入参守卫与**成对删缺失**口径，`correlation_table` 在最外层调另外两个。 + ### Changed - `core.get_rng` / `core.set_seed`：实现（L0 第 2 项）。`get_rng` 收 `int` / `np.integer` / `Generator` / `None`，其余类型抛 `TypeError`；无参调用复用同一个会推进的共享生成器。 @@ -190,6 +197,21 @@ 另新增 2 条（多列数值表能跑通、非法输入报错）。 - 测试计数：**103 passed / 39 skipped → 105 passed / 38 skipped / 1 xfail**。 +- `stats.correlation_matrix`：按**方案 A** 拆开 —— 只返回**数值**矩阵，不再返回带星标的字符串表。 + 字符串表做不了后续计算（挑高相关对、喂聚类 / PCA、画热力图），而反过来随时能格式化； + 一旦返回字符串，数值信息不可逆地丢失（见易错点 14）。 +- `stats.correlation_pvalues`：**对角线由 `0.0` 改成 `NaN`** —— 写成 0 会被下游读成"极显著"， + 而"自己和自己相关"根本没有 p 值可言。 +- `stats.correlation_pvalues`：p 值与 r 现在用**同一批观测** —— 先成对 `dropna()` 再算， + 并同样尊重 `min_periods`（不足给 NaN）。此前 `corr` 是成对删缺失、`pearsonr` 不删， + 两者基于不同的样本量，r 和 p 对不上。 +- `stats.correlation_pvalues`：只算上三角再镜像到对称位置（成对变量的 p 天然对称）， + 省掉一半重复计算。 +- `tests/test_stats.py`：`test_correlation_matrix_diagonal_is_one` 摘掉 `xfail` 转正； + 新增 8 条（数值性 + 对称性、非数值列跳过、p 值对照**解析式** `t = r·sqrt((n-2)/(1-r²))`、 + 对角线 NaN、成对删缺失口径、`min_periods`、对角线无星标、表格与数值输出一致）。 +- 测试计数：**105 passed / 38 skipped / 1 xfail → 114 passed / 38 skipped / 0 xfail**。 + ### Notes - `detect_outliers` 原有的三个 bug 加一处设计空白**已全部收口**：`threshold` 范围校验挪进 `clean_missing`、`is_series` / `is_ld_array` 提到类型转换**之前**、统计量与掩码都限制在数值列、 @@ -202,9 +224,10 @@ - 新增易错点第 10 条：**把真报错当成"编辑器误报"** —— 红线先跑 CLI 与 `pytest` 再判断（`docs/dev/personal-pitfalls.md`）。 - `stats` 这条线开工了：`describe` 落地，填掉 `examples/01_clean_and_describe.py` 的 `TODO(L0)`； `correlation_matrix` 功能上也能跑了（硬 bug 已修）。 -- `correlation_matrix` 的**返回契约待拍板**：带星标的字符串表（好贴论文）vs 数值矩阵 - （能做后续计算）。另有两处一并待定：**对角线不该带 `***`**；`df_num` 含 NaN 时 - `corr` 成对删失而 `pearsonr` 不删，r 与 p 会对不上。A/B/C 三个选项见 `STATUS.md`「下一步」第 1 条。 +- `correlation_matrix` 的返回契约已按**方案 A** 定案（决策 41）：数值归 `correlation_matrix`、 + p 值归 `correlation_pvalues`、星标表归 `correlation_table`。对角线不带星、p 与 r 同批观测两条也一并修掉。 +- `correlation_matrix` 仍只收 `DataFrame`（`describe` / `zscore` / `transform` 都收 ndarray）—— + 要不要放开由用户决定，记在 `STATUS.md`「下一步」第 5 条。 - 三处口径已拍板（2026-10-03）：`log_transform` 的 `offset` 维持现状（决策 35）； `describe` / `zscore` **默认丢掉**非数值列、与 `detect_outliers` 的"保留列"刻意不统一（决策 37）； 变异系数改用 `σ/|μ|`（决策 38）。另立决策 39 明确分工：语义口径用户拍板、工程事务 Codex 顺手补。 diff --git a/STATUS.md b/STATUS.md index bd0fe48..533045b 100644 --- a/STATUS.md +++ b/STATUS.md @@ -10,7 +10,7 @@ **第二阶段（进行中）**：按 `docs/03-roadmap.md` 的 L0 → L1 顺序**逐个实现函数**， 用户本人写 `src/`、Codex 负责审查 + `tests/`（决策 25）。 -当前 8/12：`Result`、`get_rng`/`set_seed`、`load_table`/`save_table`、`clean_missing`、`detect_outliers`、`normalize`/`standardize`、`describe`。 +当前 9/12：`Result`、`get_rng`/`set_seed`、`load_table`/`save_table`、`clean_missing`、`detect_outliers`、`normalize`/`standardize`、`describe`、`correlation_matrix`。 ## 任务清单 @@ -62,7 +62,7 @@ - [x] 推送通道：公网 22 被防火墙过滤 → 改用 Tailscale（`100.115.86.48:22`）， `~/.ssh/config` 加了一条域名别名，remote 里仍写域名 -**E. L0 实现进度（8/12）** +**E. L0 实现进度（9/12）** - [x] `core/result.py` 的 `Result` - `__repr__`：单行紧凑；`data` 只给类型与规模，**绝不**整数组打印 （十万级数组的 repr 实测仍 < 200 字符）；`success` 恒在 @@ -154,7 +154,7 @@ 决策 26 说的就是这个，实现写完了却零覆盖），改成逐条 skip 未实现的 5 条， 并新写 20 条 describe / zscore / 变异系数用例 - [x] `examples/01_clean_and_describe.py`：填掉末尾的 `TODO(L0)`，加一段 `describe` 输出 -- [ ] `stats/inference.py` 的 `correlation_matrix`（L0 第 9 项，**进行中**；**由用户本人实现**，Codex 审查 + 代修） +- [x] `stats/inference.py` 的 `correlation_matrix`（L0 第 9 项，**由用户本人实现**，Codex 审查 + 代修） - **硬 bug**：`if numeric_cols == 0:` 拿 `Index` 和标量比，得到的是布尔**数组** —— 0 列时静默不报错、1 列时碰巧能跑、**≥2 列直接 `ValueError` 崩**； 而"≥2 个数值列"正是最常见的情况，等于整个函数不可用。改成 `len(numeric_cols) == 0`（易错点 13） @@ -163,18 +163,24 @@ `_make_tuple_bunch` 运行时拼的类、存根里看不到 `pvalue`（见 `docs/04-environment.md` 坑 10）。 加了 `_pvalue(result: Any)` 小助手统一取 p 值，**一处都没挂 `# type: ignore`** - `method` 校验提到最前、清掉 3 条 F541 / 1 条 RET505 / I001，ruff 归零 - - **L0 契约待拍板**：返回值现在是**带星标的字符串表**（`1.000***`）， - 而 docstring 契约与 L0 用例都要求**数值矩阵**（对角 = 1.0）； - 该用例已按决策 30 从 `skip` 改成 `xfail(strict=False)`，见「下一步」第 1 条 - - 另两处待确认：**对角线也被标了 `***`**（自己和自己相关不该有显著性星标）； - `df_num` 含 NaN 时 `corr` 走成对删失、而 `pearsonr` 不删 → r 与 p 对不上 + - **现在返回**：`correlation_matrix` → 数值系数矩阵（对角线 1.0）； + `correlation_pvalues` → 数值 p 值矩阵（**对角线 NaN**）； + `correlation_table` → 星标字符串表（**对角线 `1.000` 不带星**）。 + - **返回契约冲突（已解决）**：上一版返回的是**带星标的字符串表**（`1.000***`）， + 而 docstring 契约与 L0 用例都要求**数值矩阵**（对角 = 1.0）；按方案 A 拆成三个函数（见下） +- [x] `correlation_matrix` 的**返回契约**按**方案 A** 落地（2026-10-03，用户拍板）： + `correlation_matrix` 只返回**数值**矩阵；星标表拆成 `correlation_table`、p 值拆成 + `correlation_pvalues`（三者共用 `_check_corr_args` 守卫与**成对删缺失**口径，见决策 41）。 + 顺带修掉两处：**对角线不再带 `***`**（p 值矩阵的对角线是 NaN、"自己和自己"没有显著性）； + **p 值与 r 用同一批观测**（都成对删缺失、都尊重 `min_periods`），不再各算各的样本量。 + p 值只算上三角再镜像，省一半重复计算。测试 **+8 条**（含用解析式 `t = r·sqrt((n-2)/(1-r²))` 对照） - [ ] 剩余 3 项：`bootstrap_ci` `monte_carlo` `save_figure` ### 验证结果 | 检查 | venv | conda `mcm-toolkit` | |---|---|---| -| `pytest` | 见下（单列） | **105 passed, 38 skipped（1 xfail）** | +| `pytest` | 见下（单列） | **114 passed, 38 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` | 跑通 | — | @@ -184,10 +190,10 @@ （修配置之前是 17 errors，全是"无法解析导入 numpy / pandas"，见决策 29 与 `docs/04-environment.md` 坑 7） - base 亦已复核：numpy 2.4.6 / scipy 1.18.1 / pandas 3.0.6（与两套环境一致） -- **144 条用例：105 passed / 38 skipped / 1 xfail**。收集数按文件：`test_data.py` 52、`test_api_surface.py` 20、 - `test_core.py` 15、`test_solve.py` 9、`test_viz.py` 8、`test_decision.py` 7、`test_stats.py` 28、 +- **152 条用例：114 passed / 38 skipped / 0 xfail**。收集数按文件：`test_data.py` 52、`test_api_surface.py` 20、 + `test_core.py` 15、`test_solve.py` 9、`test_viz.py` 8、`test_decision.py` 7、`test_stats.py` 36、 `test_models.py` 5 -- 目前 **1 条 xfail**：`test_correlation_matrix_diagonal_is_one` —— 契约（数值矩阵）与实现（带星标字符串表）冲突，口径待拍板 +- **`xfail` 再次清零**：`correlation_matrix` 的契约冲突已按方案 A 解决（决策 41），那条用例转正 - 上一轮那 7 条 xfail 已全部转正，决策 30 那套信号（XPASS → 摘 mark）跑完过一个完整闭环 - `examples/01_clean_and_describe.py` 实测跑通（exit 0）；`00_quickstart.py` 仍等 `stats` / `viz` - 38 条 skip 是等实现的算法（`train_test_split`、`describe`、`bootstrap_ci`、`monte_carlo` 等） @@ -262,25 +268,22 @@ 均值为负时 `σ/μ` 会给出负值，容易被读成"波动更小"，其实只是符号跟着均值走。取绝对值后 CV 恒非负。可比性仍只在均值远离 0 的比率尺度数据上成立（已写进 docstring 的 Notes） | | 39 | 工作约定：**语义 / 非法输入口径由用户拍板，Codex 只列选项与影响、不代选**；**ruff、`STATUS.md`、`CHANGELOG.md`、`tests/`、docstring 这类与算法无关的工程事务，Codex 顺手补掉** | | 40 | scipy 的检验结果统一走 `_pvalue(result: Any) -> float` 助手取 p 值，不逐处挂 `# type: ignore` | +| 41 | `correlation_matrix` **只返回数值矩阵**；星标表拆 `correlation_table`、p 值拆 `correlation_pvalues` | +用户 2026-10-03 选方案 A。数值矩阵能做的事（挑高相关对、喂聚类 / PCA、画热力图）字符串表都做不了，而反过来随时能格式化 —— 展示需求属于**输出端**，不该侵占返回值；一旦返回字符串，数值信息不可逆地丢失。三个函数共用 `_check_corr_args` 守卫与成对删缺失口径，`correlation_table` 在**最外层**调另外两个（计算型当主、展示型当辅）。见易错点 14 | `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. **`correlation_matrix` 的返回契约待你拍板**（这是 L0 第 9 项唯一还没定的地方）： - 现在返回的是**带星标的字符串表**（`1.000***`），方便直接贴论文； - 但 docstring 的契约与 L0 用例要的是**数值矩阵**（对角 = 1.0），两者只能留一个。 - 三个选项： - - **A（推荐）**：返回数值矩阵，另给一个 `correlation_table(corr, p)` 出星标字符串表 - —— 数值矩阵能做的事（挑高相关对、喂给聚类/PCA、画热力图）字符串表都做不了， - 而反过来随时能格式化。代价：多一个函数。 - - **B**：就返回字符串表，把 docstring 契约改掉、L0 用例改成断言 `1.000`。代价：丢数值能力。 - - **C**：返回数值矩阵，星标表挂到 `df.attrs["display"]` 之类的地方。代价：藏得深，不显式。 - 顺带两处一并定了：**对角线不该带 `***`**；`df_num` 含 NaN 时 `corr` 成对删失但 `pearsonr` 不删， - r 与 p 会对不上（要不要先 `dropna` 成对再算 p）。 -2. 之后 L0 剩余：`bootstrap_ci`（**S 档，要写函数卡片**）→ `monte_carlo`（**S 档**）→ `save_figure`（A）。 +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）。 3. L0 跑通后，找一道真题（建议 2025 C）走完整流程，反向沉淀函数。 4. 已拍板的旧口径（留档）：`log_transform` 的 `offset` 维持现状（决策 35）； `describe` / `zscore` 默认丢掉非数值列（决策 37）；变异系数用 `σ / |μ|`（决策 38）； - 分工见决策 39。 + 分工见决策 39；`correlation_matrix` 走方案 A（决策 41）。 +5. **仍未答的一个小问题**：`correlation_matrix` 现在只收 `DataFrame`， + 而 `describe` / `zscore` / `transform` 都收 ndarray —— 要不要也放开？（你说了算） diff --git a/docs/03-roadmap.md b/docs/03-roadmap.md index d8497cd..ac3e950 100644 --- a/docs/03-roadmap.md +++ b/docs/03-roadmap.md @@ -26,6 +26,9 @@ | `monte_carlo` | `solve/simulate.py` | **S** | | `save_figure` | `viz/export.py` | A | +> `correlation_matrix` 有两个**伴随函数**（决策 41）：`correlation_pvalues` 给 p 值矩阵、 +> `correlation_table` 给"系数 + 星标"的论文表格。三者共用同一套入参守卫与**成对删缺失**口径。 + ### L1 —— v0.1（共约 40 个） 在 L0 基础上补齐： diff --git a/docs/dev/personal-pitfalls.md b/docs/dev/personal-pitfalls.md index 78fd561..532665b 100644 --- a/docs/dev/personal-pitfalls.md +++ b/docs/dev/personal-pitfalls.md @@ -217,3 +217,20 @@ `if x:` 对 `DataFrame` / `Index` / 多元素 `ndarray` 会抛同一个歧义错误； ③ 新写的分支，**至少拿 0 / 1 / 多三种形状各跑一遍**，单列能过不代表多列能过。 - 首次：2026-10-03｜复发：0 + +### 14. 把"能看的表"和"能算的表"塞进同一个返回值 + +- 现象：`correlation_matrix` 返回的是一张**带显著性星标的字符串表**（`1.000***` 这种）， + 直接贴论文很爽。但同一个函数在 docstring 里承诺的是"值 = 相关系数"， + L0 用例也断言 `np.allclose(np.diag(...), 1.0)` —— 字符串表两样都做不到： + 拿到结果的人**没法**挑高相关对、喂给聚类 / PCA、画热力图。 + 于是这条用例只能先挂 `xfail`，实现和契约顶在那里。 +- 为什么会错：写的时候想的是"我最后要看到什么"（论文里的那张表）， + 没想"调用方还要拿它**算什么**"。展示需求是**输出端**的事，不该侵占返回值。 + 一旦返回字符串，数值信息就**不可逆地丢了**（你可以随时把数字格式化成字符串， + 但没法把 `"0.570***"` 还原成 `0.570` 和小数点后丢失的精度）。 +- 以后怎么防：① 一个函数只返回**一种东西** —— 要么是数据，要么是给人看的展示； + ② 判据是"这个返回值还能不能参与下一步计算"：不能，就该拆成两个函数； + ③ 拆的时候让**计算型**当主函数、**展示型**去调它（`correlation_table` 调 + `correlation_matrix` + `correlation_pvalues`），而不是反过来。 +- 首次：2026-10-03｜复发：0 diff --git a/src/mcm_toolkit/stats/__init__.py b/src/mcm_toolkit/stats/__init__.py index 4d6850f..a1bb5bb 100644 --- a/src/mcm_toolkit/stats/__init__.py +++ b/src/mcm_toolkit/stats/__init__.py @@ -11,6 +11,8 @@ from .inference import ( chi_square_test, confidence_interval, correlation_matrix, + correlation_pvalues, + correlation_table, fit_distribution, goodness_of_fit, t_test, @@ -21,6 +23,8 @@ __all__ = [ "zscore", "coefficient_of_variation", "correlation_matrix", + "correlation_pvalues", + "correlation_table", "confidence_interval", "t_test", "anova", diff --git a/src/mcm_toolkit/stats/inference.py b/src/mcm_toolkit/stats/inference.py index c4c1576..1c5abdc 100644 --- a/src/mcm_toolkit/stats/inference.py +++ b/src/mcm_toolkit/stats/inference.py @@ -18,6 +18,8 @@ __all__ = [ "chi_square_test", "confidence_interval", "correlation_matrix", + "correlation_pvalues", + "correlation_table", "fit_distribution", "goodness_of_fit", "t_test", @@ -34,12 +36,28 @@ def _pvalue(result: Any) -> float: return float(result.pvalue) +_VALID_CORR_METHODS = ("pearson", "kendall", "spearman") + + +def _check_corr_args(data: Any, method: str) -> pd.DataFrame: + """三个 `correlation_*` 共用的入参守卫，返回收窄后的**数值子表**。""" + if not isinstance(data, pd.DataFrame): + raise TypeError(f"data必须是 pd.DataFrame,而你传的是{type(data)}") + if method not in _VALID_CORR_METHODS: + raise ValueError(f"不支持的方法:{method},仅支持{list(_VALID_CORR_METHODS)}") + numeric_cols = data.select_dtypes(include=[np.number]).columns + if len(numeric_cols) == 0: + raise ValueError("没有找到数值列") + # data[Index] 在 pandas 存根里是 DataFrame | Series 重载，按运行时真值收窄（坑 8） + return cast(pd.DataFrame, data[numeric_cols]) + + def correlation_matrix( data: Any, method: Literal["pearson", "kendall", "spearman"] = "pearson", min_periods: int = 1, -) -> Any: - """相关矩阵 + 显著性星标。 +) -> pd.DataFrame: + """相关系数矩阵（**数值**，不带星标）。 Parameters ---------- @@ -49,12 +67,12 @@ def correlation_matrix( `"pearson"` / `"kendall"` / `"spearman"`。 min_periods: 每对变量至少要有这么多个**共同非缺失**的观测才给相关系数，透传给 `DataFrame.corr`。 + （pandas 对每对变量都是**成对删缺失**，不是整行删。） Returns ------- - DataFrame，值 = 相关系数；`min_periods` 控制缺失值容差。 - ⚠️ 当前实现返回的是**带星标的字符串表**（`1.000***` 这种），与上面这句契约不一致 —— - 口径待定，见 `STATUS.md`。 + pandas.DataFrame + 对称矩阵，对角线恒为 `1.0`；观测不足的格子是 NaN。 Raises ------ @@ -65,58 +83,121 @@ def correlation_matrix( Notes ----- - 星标读法：`***` p < 0.001，`**` p < 0.01，`*` p < 0.05，其余不带星。 + 要"系数 + 显著性星标"的论文表格用 `correlation_table`；要 p 值矩阵用 `correlation_pvalues`。 + 之所以拆开：这张表是**数值**，能挑高相关对、喂聚类 / PCA、画热力图；星标字符串表这些都做不了。 Examples -------- >>> correlation_matrix(df).round(3) # doctest: +SKIP """ - if not isinstance(data, pd.DataFrame): - raise TypeError(f"data必须是 pd.DataFrame,而你传的是{type(data)}") + df_num = _check_corr_args(data, method) + return df_num.corr(method=method, min_periods=min_periods) - valid_method = ["pearson", "kendall", "spearman"] - if method not in valid_method: - raise ValueError(f"不支持的方法:{method},仅支持{valid_method}") - numeric_cols = data.select_dtypes(include=[np.number]).columns - if len(numeric_cols) == 0: - raise ValueError("没有找到数值列") - # data[Index] 在 pandas 存根里是 DataFrame | Series 重载，按运行时真值收窄（坑 8） - df_num = cast(pd.DataFrame, data[numeric_cols]) +def correlation_pvalues( + data: Any, + method: Literal["pearson", "kendall", "spearman"] = "pearson", + min_periods: int = 1, +) -> pd.DataFrame: + """`correlation_matrix` 对应的 **p 值矩阵**（数值，不含星标）。 - corr_matrix = df_num.corr(method=method, min_periods=min_periods) + Parameters + ---------- + data / method / min_periods: + 含义与 `correlation_matrix` 完全一致，两者一定同形。 - p_values = pd.DataFrame(0.0, index=corr_matrix.index, columns=corr_matrix.columns, dtype=float) - for i in range(len(corr_matrix.columns)): - for j in range(len(corr_matrix.columns)): - if i == j: + Returns + ------- + pandas.DataFrame + 对称矩阵。**对角线是 NaN** —— "自己和自己相关"没有 p 值可言， + 写成 `0.0` 会被下游读成"极显著"。 + + Raises + ------ + TypeError / ValueError: + 同 `correlation_matrix`。 + + Notes + ----- + 每对变量都先**成对删缺失**再算，与 `DataFrame.corr` 的口径一致 —— + 否则 r 和 p 会基于不同的样本量，对不上；成对有效观测少于 `min_periods` 时给 NaN。 + 只算上三角再镜像，成对变量的 p 天然对称，省掉一半重复计算。 + """ + df_num = _check_corr_args(data, method) + corr = df_num.corr(method=method, min_periods=min_periods) + + p_values = pd.DataFrame(np.nan, index=corr.index, columns=corr.columns, dtype=float) + for i in range(len(corr.columns)): + for j in range(i + 1, len(corr.columns)): + pair = df_num.iloc[:, [i, j]].dropna() + if len(pair) < min_periods: continue + x = pair.iloc[:, 0] + y = pair.iloc[:, 1] if method == "pearson": - p = _pvalue(stats.pearsonr(df_num.iloc[:, i], df_num.iloc[:, j])) + p = _pvalue(stats.pearsonr(x, y)) elif method == "kendall": - p = _pvalue(stats.kendalltau(df_num.iloc[:, i], df_num.iloc[:, j])) + p = _pvalue(stats.kendalltau(x, y)) else: - p = _pvalue(stats.spearmanr(df_num.iloc[:, i], df_num.iloc[:, j])) + p = _pvalue(stats.spearmanr(x, y)) p_values.iloc[i, j] = p + p_values.iloc[j, i] = p + return p_values - def add_stars(val: float, p: float) -> str: - # 保留三位小数 - val_str = f"{val:.3f}" - if p < 0.001: - return val_str + "***" - if p < 0.01: - return val_str + "**" - if p < 0.05: - return val_str + "*" - return val_str - formatted_corr = corr_matrix.copy().astype(str) - for i in range(len(corr_matrix)): - for j in range(len(corr_matrix)): - formatted_corr.iloc[i, j] = add_stars(corr_matrix.iloc[i, j], p_values.iloc[i, j]) +def _star(value: float, p: float) -> str: + """把系数与 p 值拼成 `"0.812**"`。p 为 NaN（对角线 / 观测不足）时不加星标。""" + text = f"{value:.3f}" + if p < 0.001: + return text + "***" + if p < 0.01: + return text + "**" + if p < 0.05: + return text + "*" + return text - return formatted_corr +def correlation_table( + data: Any, + method: Literal["pearson", "kendall", "spearman"] = "pearson", + min_periods: int = 1, +) -> pd.DataFrame: + """可直接贴论文的相关性表：**系数 + 显著性星标**。 + + Parameters + ---------- + data / method / min_periods: + 含义与 `correlation_matrix` 完全一致。 + + Returns + ------- + pandas.DataFrame + 每个格子是字符串，形如 `"0.812**"`。 + 星标：`***` p < 0.001、`**` p < 0.01、`*` p < 0.05，其余不带星。 + **对角线只显示 `1.000` 不带星标** —— 自己和自己没有显著性可言。 + + Raises + ------ + TypeError / ValueError: + 同 `correlation_matrix`。 + + Notes + ----- + 这张表**只能看**：格子里是字符串，做不了后续计算。 + 要数值请用 `correlation_matrix`（系数）与 `correlation_pvalues`（p 值）。 + + Examples + -------- + >>> print(correlation_table(df)) # doctest: +SKIP + """ + corr = correlation_matrix(data, method=method, min_periods=min_periods) + p_values = correlation_pvalues(data, method=method, min_periods=min_periods) + + table = corr.astype(str) + for i in range(len(corr)): + for j in range(len(corr.columns)): + table.iloc[i, j] = _star(corr.iloc[i, j], p_values.iloc[i, j]) + return table def confidence_interval( x: Any, confidence: float = 0.95, diff --git a/tests/test_stats.py b/tests/test_stats.py index 1df10df..ff94ad3 100644 --- a/tests/test_stats.py +++ b/tests/test_stats.py @@ -15,6 +15,17 @@ import pytest _NOT_IMPLEMENTED = "未实现：见 docs/03-roadmap.md L1" +def _pearson_p_from_r(r: float, n: int) -> float: + """Pearson 相关系数的 p 值**解析式**（测解析解，不测实现细节）： + + ``t = r * sqrt((n - 2) / (1 - r^2))``，``p = 2 * P(T > |t|)``，``T ~ t(n - 2)``。 + """ + from scipy import stats as sps + + t_stat = r * np.sqrt((n - 2) / (1 - r**2)) + 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% 容差）。""" @@ -59,18 +70,27 @@ def test_confidence_interval_matches_scipy(rng): assert np.allclose([lo, hi], ref, rtol=1e-6) -@pytest.mark.xfail( - strict=False, - reason="契约要求返回**数值**矩阵（对角 = 1.0），实现返回带星标的字符串表；口径待定，见 STATUS.md", -) def test_correlation_matrix_diagonal_is_one(rng): from mcm_toolkit.stats import correlation_matrix x = pd.DataFrame(rng.normal(size=(200, 4)), columns=["a", "b", "c", "d"]) c = correlation_matrix(x) + assert np.allclose(np.diag(np.asarray(c)), 1.0) +def test_correlation_matrix_is_numeric_and_symmetric(rng): + """方案 A：返回**数值**矩阵（不是带星标的字符串表），能直接参与后续计算。""" + from mcm_toolkit.stats import correlation_matrix + + x = pd.DataFrame(rng.normal(size=(200, 3)), columns=["a", "b", "c"]) + c = correlation_matrix(x) + + assert all(np.issubdtype(dtype, np.floating) for dtype in c.dtypes) + assert np.allclose(np.asarray(c), np.asarray(x.corr(method="pearson", min_periods=1))) + assert np.allclose(np.asarray(c), np.asarray(c.T)) + + def test_correlation_matrix_runs_on_multicolumn_frame(rng): """回归：`if numeric_cols == 0` 拿 Index 跟 0 比，多列数值表当场 ValueError 崩掉。""" from mcm_toolkit.stats import correlation_matrix @@ -83,6 +103,13 @@ def test_correlation_matrix_runs_on_multicolumn_frame(rng): assert list(out.index) == ["a", "b", "c"] +def test_correlation_matrix_skips_non_numeric_columns(rng): + from mcm_toolkit.stats import correlation_matrix + + x = pd.DataFrame({"a": rng.normal(size=30), "b": rng.normal(size=30), "city": ["X"] * 30}) + assert list(correlation_matrix(x).columns) == ["a", "b"] + + def test_correlation_matrix_rejects_bad_input(): from mcm_toolkit.stats import correlation_matrix @@ -94,6 +121,85 @@ def test_correlation_matrix_rejects_bad_input(): correlation_matrix(pd.DataFrame({"a": [1.0, 2.0]}), method="bogus") # type: ignore[arg-type] +def test_correlation_pvalues_matches_analytic_formula(rng): + from mcm_toolkit.stats import correlation_matrix, correlation_pvalues + + x = pd.DataFrame(rng.normal(size=(80, 3)), columns=["a", "b", "c"]) + corr = correlation_matrix(x) + p = correlation_pvalues(x) + + for i in ["a", "b", "c"]: + for j in ["a", "b", "c"]: + if i == j: + continue + assert p.loc[i, j] == pytest.approx(_pearson_p_from_r(corr.loc[i, j], 80)) + + +def test_correlation_pvalues_diagonal_is_nan(rng): + """对角线写 0 会被下游读成"极显著"；"自己和自己相关"没有 p 值，必须是 NaN。""" + from mcm_toolkit.stats import correlation_pvalues + + p = correlation_pvalues(pd.DataFrame(rng.normal(size=(50, 2)), columns=["a", "b"])) + + assert np.isnan(p.loc["a", "a"]) + assert np.isnan(p.loc["b", "b"]) + + +def test_correlation_pvalues_uses_the_same_observations_as_r(rng): + """含缺失时 p 要和 r 用**同一批**观测：成对删缺失，而不是整行删、也不是干脆不删。""" + from mcm_toolkit.stats import correlation_pvalues + + a = rng.normal(size=40) + b = rng.normal(size=40) + a[[3, 7, 11]] = np.nan + x = pd.DataFrame({"a": a, "b": b}) + + pair = x.dropna() + r = float(np.corrcoef(pair["a"], pair["b"])[0, 1]) + + assert correlation_pvalues(x).loc["a", "b"] == pytest.approx(_pearson_p_from_r(r, len(pair))) + + +def test_correlation_pvalues_respects_min_periods(rng): + from mcm_toolkit.stats import correlation_pvalues + + a = rng.normal(size=10) + a[:5] = np.nan + x = pd.DataFrame({"a": a, "b": rng.normal(size=10)}) + + assert np.isnan(correlation_pvalues(x, min_periods=6).loc["a", "b"]) + assert np.isfinite(correlation_pvalues(x, min_periods=5).loc["a", "b"]) + + +def test_correlation_table_diagonal_has_no_stars(): + """对角线是"自己和自己"，不该有显著性星标；完全负相关则必然是 ***。""" + from mcm_toolkit.stats import correlation_table + + x = pd.DataFrame({"a": np.arange(30.0), "b": np.arange(30.0) * -2.0}) + table = correlation_table(x) + + assert table.loc["a", "a"] == "1.000" + assert table.loc["a", "b"] == "-1.000***" + + +def test_correlation_table_matches_the_numeric_outputs(rng): + from mcm_toolkit.stats import correlation_matrix, correlation_pvalues, correlation_table + + x = pd.DataFrame(rng.normal(size=(60, 3)), columns=["a", "b", "c"]) + corr = correlation_matrix(x) + p = correlation_pvalues(x) + table = correlation_table(x) + + assert list(table.columns) == list(corr.columns) + assert table.loc["a", "a"] == "1.000" + assert np.isnan(p.loc["a", "a"]) + + pv = p.loc["a", "b"] + stars = "***" if pv < 0.001 else "**" if pv < 0.01 else "*" if pv < 0.05 else "" + value = corr.loc["a", "b"] + assert table.loc["a", "b"] == f"{value:.3f}" + stars + + # ------------------------------------------------------ stats: describe