---
name: stats_ml
description: 统计建模与机器学习(sklearn / statsmodels / PyMC 三库合一,场景导航选库)。✅ 触发:配方-性能建模、DoE 响应面、强度预测、特征重要性、假设检验、置信区间、贝叶斯小样本估计、聚类异常实验检测。⛔ 不触发:用户只问描述性统计(均值方差),pandas 一行搞定;深度学习场景(走单独的 torch / lightning skill,本仓库未集成)。
---

# stats_ml

服务建材院典型场景:**掺合料配比 → 性能(强度 / 流动度 / 凝结时间)** 的建模与推断。三库分工:

| 你要做 | 用 | 一句话理由 |
|---|---|---|
| 要 p-value、置信区间、假设检验 | **statsmodels** | OLS / ANOVA / 假设检验是 statsmodels 主场,sklearn 不给 p-value |
| 要预测精度高,可解释性次要 | **sklearn** | RandomForest / GBDT / XGBoost 是性能预测的现代基线 |
| 样本 < 30 且要不确定度估计 | **PyMC** | 贝叶斯给出参数后验分布,小样本下比频率派的置信区间更稳 |
| 要可解释的线性关系 | statsmodels OLS / sklearn LinearRegression | 都可以,统计推断走 statsmodels,纯预测走 sklearn |
| 要 DoE 响应面拟合 | statsmodels(二次回归 + ANOVA) | 主流做法,系数显著性能直接读 |
| 要聚类找异常配方 | sklearn (KMeans / DBSCAN) | 标准工具 |
| 要降维可视化(配方空间) | sklearn (PCA / t-SNE) | 标准工具 |

## 何时用 / 何时不用

✅ 用:
- 用户给一张配方-性能表(N 行配方 + 性能列),要建模 / 预测 / 推断
- 要回答"哪个因素影响最大"(特征重要性 / 系数显著性)
- 要回答"这个新配方预测强度多少 + 置信区间多大"
- 要找异常实验点 / 配方聚类

⛔ 不用:
- 用户只要看均值方差直方图 → pandas + matplotlib 直接搞,杀鸡不用牛刀
- 用户给的"数据"只有 3-5 个点 → 任何 ML / 统计都不可靠,告诉用户先补数据
- 深度学习需求(CNN / Transformer)→ 不在本 skill 范围

## 准备

```python
# 按需 import,不要全 import
import numpy as np
import pandas as pd

# sklearn
from sklearn.model_selection import train_test_split, KFold, StratifiedKFold, GridSearchCV
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import r2_score, mean_squared_error
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor

# statsmodels
import statsmodels.api as sm
import statsmodels.formula.api as smf

# PyMC(重,按需 import)
# import pymc as pm
# import arviz as az
```

## 典型工作流

### A. 配方-性能回归(sklearn,要预测精度)

```python
# 数据:每行一个配方,X 是组分比例,y 是 28d 抗压强度
X = df[["cement", "fly_ash", "slag", "water", "sand", "stone"]].values
y = df["strength_28d"].values

# Pipeline 把预处理和模型绑死 —— 反模式见末尾
pipe = Pipeline([
    ("scaler", StandardScaler()),       # 线性 / 距离类模型必须;树模型可省
    ("model", RandomForestRegressor(n_estimators=300, random_state=42)),
])

# 数据小(< 200)用 5-fold CV,不要单次 train/test
from sklearn.model_selection import cross_val_score
scores = cross_val_score(pipe, X, y, cv=5, scoring="r2")
print(f"R2 mean={scores.mean():.3f} std={scores.std():.3f}")

# 训完看特征重要性
pipe.fit(X, y)
imp = pipe.named_steps["model"].feature_importances_
for name, val in sorted(zip(df.columns, imp), key=lambda t: -t[1]):
    print(f"{name}: {val:.3f}")
```

### B. 显著性分析(statsmodels,要 p-value)

```python
import statsmodels.formula.api as smf

# formula API 让模型 spec 看起来跟 R 一样
model = smf.ols(
    formula="strength_28d ~ cement + fly_ash + slag + water + water:cement",  # 含交互
    data=df,
).fit()
print(model.summary())  # 系数、std err、t、p>|t|、95% CI 全在这

# ANOVA 比较两个嵌套模型
m1 = smf.ols("strength_28d ~ cement + fly_ash", data=df).fit()
m2 = smf.ols("strength_28d ~ cement + fly_ash + slag", data=df).fit()
from statsmodels.stats.anova import anova_lm
print(anova_lm(m1, m2))  # slag 加进去显著吗
```

### C. DoE 响应面(statsmodels 二次回归)

```python
# 中心复合设计后,做二阶模型
model = smf.ols(
    formula="strength ~ x1 + x2 + x3 + I(x1**2) + I(x2**2) + I(x3**2) + x1:x2 + x1:x3 + x2:x3",
    data=df,
).fit()
print(model.summary())
# 看哪些二阶项 / 交互项显著,定优化方向
```

### D. 小样本贝叶斯(PyMC,N<30 要不确定度)

```python
import pymc as pm
import arviz as az

with pm.Model() as model:
    # 系数先验:弱信息正态
    intercept = pm.Normal("intercept", mu=0, sigma=10)
    beta = pm.Normal("beta", mu=0, sigma=5, shape=X.shape[1])
    sigma = pm.HalfNormal("sigma", sigma=10)

    mu = intercept + pm.math.dot(X, beta)
    y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

    trace = pm.sample(2000, tune=1000, chains=4, target_accept=0.95)

az.summary(trace)            # 后验均值 / std / 94% HDI
az.plot_posterior(trace)     # 后验分布图
```

### E. 配方聚类找异常实验(sklearn KMeans / DBSCAN)

```python
from sklearn.cluster import DBSCAN
from sklearn.preprocessing import StandardScaler

Xs = StandardScaler().fit_transform(X)
labels = DBSCAN(eps=0.5, min_samples=5).fit_predict(Xs)
# labels == -1 的点 = 离群,可能是录入错或工艺异常,人工核实
```

## 反模式

### sklearn

1. **预处理在 pipeline 外** —— 交叉验证会数据泄漏(test fold 的统计量被 train 看到)
2. **scaler 在 test 上 fit** —— 一律 `fit(train)` / `transform(test)`,Pipeline 自动保证
3. **分类任务忘 stratified split** —— 类别不平衡时 `train_test_split(..., stratify=y)`
4. **不必要地 scale 树模型** —— RF / GBDT / XGBoost 不受单调变换影响
5. **忽略 convergence warning** —— LogisticRegression / SVM 不收敛要加 `max_iter` 或 scale 特征
6. **小数据集还 80/20 切** —— N<200 用 5-fold CV,N<50 用 LOOCV

### statsmodels

7. **直接读系数不看 condition number** —— `summary()` 末尾的 condition number > 30 = 严重共线性,系数解读无效
8. **多重比较忘修正** —— 同时检验多个系数显著性时,用 Bonferroni / FDR(`statsmodels.stats.multitest`)
9. **OLS 残差不检验** —— `model.resid` 不正态 / 异方差时,OLS 推断不可靠,改 robust SE(`fit(cov_type="HC3")`)

### PyMC

10. **不看 R-hat / ESS** —— `az.summary(trace)` 的 `r_hat` > 1.01 或 `ess_bulk` < 400 = 没收敛,加 tune / 调 target_accept,不能用
11. **强先验当数据用** —— 先验定太窄会主导后验,先用弱信息先验试
12. **PyMC 装包巨大** —— 仅小样本需要时才用,大样本(N>500)频率派结果一样可信且快 100 倍

### 通用

13. **没看数据分布就建模** —— 先 `df.describe()` + `df.hist()` 扫一遍,有 NaN / 量纲极差 / 长尾分布要先处理
14. **R2 > 0.95 就开香槟** —— 大概率是数据泄漏 / 过拟合 / 训练集等于测试集,先查
15. **预测新配方不报置信区间** —— 工程决策不能只给点估计,sklearn 用 `RandomForestRegressor` 的 tree-level prediction std,statsmodels / PyMC 直接给 CI / HDI
16. **数据点不到 10 还硬上 ML** —— 告诉用户先做 DoE 扩样本,再建模

## 依赖

`requirements.txt` 加:
```
scikit-learn>=1.4.0
statsmodels>=0.14.0
pymc>=5.10.0       # 重,装包带 pytensor,可选
arviz>=0.17.0      # PyMC 后验诊断
```

装:`.venv/Scripts/pip install scikit-learn statsmodels`
PyMC 单独装:`.venv/Scripts/pip install pymc arviz`(确认要做贝叶斯再装,首次装包 5 分钟起)