zcbot/skills/stats_ml/SKILL.md

name	description
stats_ml	统计建模与机器学习(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 范围

准备

# 按需 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,要预测精度)

# 数据:每行一个配方,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)

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 二次回归)

# 中心复合设计后,做二阶模型
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 要不确定度)

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)

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 分钟起)