model-validation

name	model-validation
description	模型验证：交叉验证、留出法、残差分析、与已知解对比、假设检验。触发词: 模型验证、交叉验证、残差分析、model validation、留出法、误差分析、假设检验。
argument-hint	["model-file-or-results-path"]
allowed-tools	Bash(*), Read, Write, Edit, Grep, Glob, Agent, mcp__codex__codex, mcp__codex__codex-reply

模型验证

执行描述: $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。

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% 标记为过拟合风险。

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. 残差分析是数模论文中评委重点关注的内容，图表质量要求高。

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()