model-validation

# 模型验证 执行描述: $ARGUMENTS ## Constants - **OUTPUT_REPORT = `MODEL_VALIDATION_REPORT.md`** — 主输出文件，供 `paper-write`、`model-review` 使用。 - **FIGURE_DIR = `figures/validation/`** — 验证图表存放目录。 - **ARTIFACT_DIR = `artifacts/`** — 中间结构化文件存放目录。 - **DEFAULT_K_FOLDS = `5`** — 交叉验证默认折数。 - **DEFAULT_TEST_RATIO = `0.20`** — 留出法默认测试集比例。 - **DEFAULT_RANDOM_SEED = `42`** — 随机种子，确保可复现。 - **RESIDUAL_NORMALITY_ALPHA = `0.05`** — 残差正态性检验显著性水平。 - **HOMOSCEDASTICITY_ALPHA = `0.05`** — 残差方差齐性检验显著性水平。 - **INDEPENDENCE_ALPHA = `0.05`** — 残差独立性检验显著性水平（Durbin-Watson）。 - **MAX_BENCHMARK_SOLUTIONS = `5`** — 与已知解对比时最多使用 5 个基准。 - **REVIEWER_MODEL = `gpt-5.4`** — Codex MCP 交叉验证模型。 - **REASONING_EFFORT = `high`** — 模型验证审查使用高推理强度。 - **ACCEPTABLE_R2_THRESHOLD = `0.70`** — R-squared 低于 0.70 需要在报告中标记为"拟合不足"。 - **MAX_COMPUTE_MINUTES = `20`** — 单次验证计算的时间上限。 ## Workflow ### Phase 1: 定位模型与数据 Input: `$ARGUMENTS`、已有模型文件和数据。 Output: `artifacts/validation_scope.json`（模型信息与验证计划）。 1. 如果 `$ARGUMENTS` 是文件路径，读取模型代码或结果文件。 2. 如果无具体参数，扫描以下文件获取模型信息： - `MODEL_REPORT.md`、`SOLVE_PLAN.md`、`FINAL_PROPOSAL.md` - `scripts/`、`src/` 目录中的模型实现代码 - `results/` 目录中的求解结果 - `data/cleaned/` 目录中的清洗后数据 3. 识别模型类型并确定适用的验证方法： - **回归模型**: 留出法、K-fold CV、残差分析、R-squared、RMSE - **分类模型**: 留出法、K-fold CV、混淆矩阵、ROC/AUC - **优化模型**: 与已知最优解对比、松弛界对比、收敛性分析 - **预测模型**: 时间序列留出法、滑动窗口验证、预测区间 - **仿真模型**: 与实际数据对比、参数敏感性、极端场景测试 4. 确定因变量和自变量。 5. 确认数据量是否支持所选验证方法（如 K-fold 至少需要 5K 个样本才有意义）。 6. 制定验证计划，保存到 `artifacts/validation_scope.json`。 ```python import json from pathlib import Path def determine_validation_strategy(model_type: str, n_samples: int, n_features: int) -> dict: """Determine appropriate validation methods based on model type and data size.""" strategy = { "model_type": model_type, "n_samples": n_samples, "n_features": n_features, "methods": [], } if model_type in ["regression", "classification"]: if n_samples >= 500: strategy["methods"].append({"name": "k_fold_cv", "k": 5}) if n_samples >= 100: strategy["methods"].append({"name": "holdout", "test_ratio": 0.20}) strategy["methods"].append({"name": "residual_analysis"}) elif model_type == "optimization": strategy["methods"].append({"name": "benchmark_comparison"}) strategy["methods"].append({"name": "convergence_analysis"}) strategy["methods"].append({"name": "relaxation_bound"}) elif model_type == "prediction": strategy["methods"].append({"name": "temporal_holdout", "test_ratio": 0.20}) strategy["methods"].append({"name": "sliding_window", "window_size": "auto"}) strategy["methods"].append({"name": "prediction_interval"}) elif model_type == "simulation": strategy["methods"].append({"name": "real_data_comparison"}) strategy["methods"].append({"name": "extreme_scenario_test"}) # Universal methods strategy["methods"].append({"name": "hypothesis_tests"}) return strategy Path("artifacts").mkdir(exist_ok=True) # strategy = determine_validation_strategy("regression", n_samples=1000, n_features=10) # Path("artifacts/validation_scope.json").write_text( # json.dumps(strategy, ensure_ascii=False, indent=2), encoding="utf-8" # ) ``` ### Phase 2: 留出法与交叉验证 Input: 模型、数据、验证计划。 Output: `artifacts/cv_results.json`、性能指标汇总。 1. **留出法 (Holdout)**: - 按 `DEFAULT_TEST_RATIO` (20%) 划分训练集和测试集。 - 使用 `DEFAULT_RANDOM_SEED` 确保可复现。 - 在训练集上拟合模型，在测试集上评估。 - 记录训练集和测试集的性能差异（检测过拟合）。 2. **K-fold 交叉验证**: - 使用 `DEFAULT_K_FOLDS` (5) 折。 - 对每折记录性能指标。 - 计算均值和标准差。 - 如果标准差过大（CV > 0.2），提示模型不稳定。 3. **时间序列验证**（如适用）: - 不使用随机划分，按时间顺序分割。 - 使用扩展窗口或滑动窗口方法。 - 报告每个窗口的性能。 4. 性能指标根据模型类型选择： - 回归：R-squared、Adjusted R-squared、RMSE、MAE、MAPE - 分类：Accuracy、Precision、Recall、F1、AUC-ROC - 优化：目标函数值、约束违反度、与最优解的差距 5. 如果 R-squared < `ACCEPTABLE_R2_THRESHOLD`，在报告中标红。 6. 对比训练集和测试集性能，差距 > 10% 标记为过拟合风险。 ```python import numpy as np from sklearn.model_selection import KFold, train_test_split from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error def holdout_validation(model, X, y, test_ratio=0.20, seed=42): """Run holdout validation.""" X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_ratio, random_state=seed) model.fit(X_train, y_train) y_pred_train = model.predict(X_train) y_pred_test = model.predict(X_test) return { "train": { "r2": round(float(r2_score(y_train, y_pred_train)), 4), "rmse": round(float(np.sqrt(mean_squared_error(y_train, y_pred_train))), 4), "mae": round(float(mean_absolute_error(y_train, y_pred_train)), 4), }, "test": { "r2": round(float(r2_score(y_test, y_pred_test)), 4), "rmse": round(float(np.sqrt(mean_squared_error(y_test, y_pred_test))), 4), "mae": round(float(mean_absolute_error(y_test, y_pred_test)), 4), }, "n_train": len(X_train), "n_test": len(X_test), "overfit_risk": abs(r2_score(y_train, y_pred_train) - r2_score(y_test, y_pred_test)) > 0.10, } def kfold_validation(model, X, y, k=5, seed=42): """Run K-fold cross-validation.""" kf = KFold(n_splits=k, shuffle=True, random_state=seed) fold_results = [] for fold_idx, (train_idx, test_idx) in enumerate(kf.split(X)): X_train, X_test = X[train_idx], X[test_idx] y_train, y_test = y[train_idx], y[test_idx] model.fit(X_train, y_train) y_pred = model.predict(X_test) fold_results.append({ "fold": fold_idx + 1, "r2": round(float(r2_score(y_test, y_pred)), 4), "rmse": round(float(np.sqrt(mean_squared_error(y_test, y_pred))), 4), "mae": round(float(mean_absolute_error(y_test, y_pred)), 4), }) r2_values = [f["r2"] for f in fold_results] return { "k": k, "folds": fold_results, "r2_mean": round(float(np.mean(r2_values)), 4), "r2_std": round(float(np.std(r2_values)), 4), "r2_cv": round(float(np.std(r2_values) / np.mean(r2_values)), 4) if np.mean(r2_values) != 0 else None, "stable": float(np.std(r2_values) / max(abs(np.mean(r2_values)), 1e-10)) < 0.20, } ``` ### Phase 3: 残差分析 Input: 模型预测值、实际值。 Output: `artifacts/residual_analysis.json`、残差诊断图。 1. 计算残差：`residuals = y_actual - y_predicted`。 2. **正态性检验**: - Shapiro-Wilk 检验（n < 5000）或 Kolmogorov-Smirnov 检验（n >= 5000）。 - Q-Q 图（正态概率图）。 - 若 p-value < `RESIDUAL_NORMALITY_ALPHA`，残差不服从正态分布。 3. **方差齐性检验**: - 残差 vs 拟合值散点图，检查扇形模式。 - Breusch-Pagan 检验或 White 检验。 - 若 p-value < `HOMOSCEDASTICITY_ALPHA`，存在异方差性。 4. **独立性检验**: - 残差 vs 观测顺序图，检查自相关模式。 - Durbin-Watson 统计量（时序数据）。 - DW 值接近 2 表示无自相关。 5. **残差分布图**: 直方图+KDE。 6. **Residuals vs Fitted 图**: 检查非线性模式。 7. **Scale-Location 图**: 检查方差齐性。 8. 汇总四大诊断结果，判断模型假设是否成立。 9. 如果假设严重违反，建议模型改进方向（变换、加权、换模型）。 10. 残差分析是数模论文中评委重点关注的内容，图表质量要求高。 ```python import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from scipy import stats def residual_analysis(y_actual, y_predicted, fig_dir="figures/validation"): """Comprehensive residual analysis.""" residuals = np.array(y_actual) - np.array(y_predicted) fitted = np.array(y_predicted) n = len(residuals) # Normality test if n < 5000: stat_norm, p_norm = stats.shapiro(residuals) norm_test = "Shapiro-Wilk" else: stat_norm, p_norm = stats.kstest(residuals, "norm", args=(np.mean(residuals), np.std(residuals))) norm_test = "Kolmogorov-Smirnov" # Independence test (Durbin-Watson) diff = np.diff(residuals) dw = float(np.sum(diff ** 2) / np.sum(residuals ** 2)) # Standardized residuals std_residuals = (residuals - np.mean(residuals)) / np.std(residuals) report = { "n_observations": n, "residual_mean": round(float(np.mean(residuals)), 6), "residual_std": round(float(np.std(residuals)), 6), "normality": { "test": norm_test, "statistic": round(float(stat_norm), 6), "p_value": round(float(p_norm), 6), "is_normal": bool(p_norm >= 0.05), }, "independence": { "durbin_watson": round(dw, 4), "interpretation": "no autocorrelation" if 1.5 < dw < 2.5 else "possible autocorrelation", }, } # Generate diagnostic plots Path(fig_dir).mkdir(parents=True, exist_ok=True) fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # 1. Residuals vs Fitted axes[0, 0].scatter(fitted, residuals, alpha=0.5, s=10, color="#2196F3") axes[0, 0].axhline(y=0, color="red", linestyle="--", linewidth=1) axes[0, 0].set_xlabel("Fitted Values") axes[0, 0].set_ylabel("Residuals") axes[0, 0].set_title("Residuals vs Fitted") # 2. Q-Q Plot stats.probplot(residuals, dist="norm", plot=axes[0, 1]) axes[0, 1].set_title("Normal Q-Q Plot") # 3. Scale-Location axes[1, 0].scatter(fitted, np.sqrt(np.abs(std_residuals)), alpha=0.5, s=10, color="#4CAF50") axes[1, 0].set_xlabel("Fitted Values") axes[1, 0].set_ylabel("sqrt(|Standardized Residuals|)") axes[1, 0].set_title("Scale-Location Plot") # 4. Residual Histogram axes[1, 1].hist(residuals, bins=30, density=True, alpha=0.7, color="#FF9800", edgecolor="white") x_range = np.linspace(residuals.min(), residuals.max(), 100) axes[1, 1].plot(x_range, stats.norm.pdf(x_range, np.mean(residuals), np.std(residuals)), color="red", linewidth=2, label="Normal fit") axes[1, 1].set_title("Residual Distribution") axes[1, 1].legend()