sensitivity-analysis

name	sensitivity-analysis
description	灵敏度分析：参数扰动、单因素/多因素分析、Monte Carlo 模拟、龙卷风图/蛛网图。触发词: 灵敏度分析、参数敏感性、sensitivity analysis、Monte Carlo、鲁棒性测试、参数扰动。
argument-hint	["model-file-or-parameter-list"]
allowed-tools	Bash(*), Read, Write, Edit, Grep, Glob, Agent, mcp__codex__codex, mcp__codex__codex-reply

灵敏度分析

执行描述: $ARGUMENTS

Constants

• OUTPUT_REPORT = SENSITIVITY_ANALYSIS_REPORT.md — 主输出文件，供 paper-write、model-review 使用。
• FIGURE_DIR = figures/sensitivity/ — 灵敏度图表存放目录。
• ARTIFACT_DIR = artifacts/ — 中间结构化文件存放目录。
• DEFAULT_PERTURBATION_RANGE = 0.20 — 默认参数扰动幅度为基准值的 +/-20%。
• DEFAULT_PERTURBATION_STEPS = 11 — 单因素分析默认取 11 个等距点（含基准值）。
• MAX_PARAMETERS = 15 — 超过 15 个参数时先做初筛（Morris 方法或相关性排序）。
• MONTE_CARLO_SAMPLES = 10000 — Monte Carlo 模拟默认采样次数。
• MONTE_CARLO_MIN_SAMPLES = 1000 — 采样次数不得低于 1000。
• CONFIDENCE_LEVEL = 0.95 — 置信区间默认 95%。
• SIGNIFICANCE_THRESHOLD = 0.05 — 灵敏度指数超过 0.05 才视为"敏感参数"。
• REVIEWER_MODEL = gpt-5.4 — Codex MCP 交叉验证模型。
• REASONING_EFFORT = high — 灵敏度分析验证使用高推理强度。
• MAX_COMPUTE_MINUTES = 30 — 单次灵敏度分析的计算时间上限。

Workflow

Phase 1: 定位模型与参数

Input: $ARGUMENTS、已有模型文件和求解结果。 Output: artifacts/sensitivity_scope.json（参数列表与基准值）。

1. 如果 $ARGUMENTS 是文件路径，读取模型代码或参数文件。
2. 如果无具体参数，扫描以下文件获取模型信息：
   • MODEL_REPORT.md、SOLVE_PLAN.md、FINAL_PROPOSAL.md
   • scripts/、src/ 目录中的模型实现代码
   • results/ 目录中的求解结果
3. 从模型代码中提取所有可调参数：
   • 模型参数（权重、系数、阈值）
   • 超参数（步长、迭代次数、正则化参数）
   • 假设参数（折现率、增长率、概率分布参数）
   • 边界条件和约束常数
4. 为每个参数记录：名称、当前基准值、物理含义、合理取值范围、数据类型。
5. 如果参数过多（> MAX_PARAMETERS），使用 Morris 筛选法或专家判断先缩减范围。
6. 识别目标函数（要分析灵敏度的输出量）。
7. 确认模型可以在修改参数后重新求解（函数化封装）。

import json
from pathlib import Path
import re

def extract_parameters_from_code(code_path: str) -> list:
    """Extract tunable parameters from Python model code."""
    code = Path(code_path).read_text(encoding="utf-8")
    params = []
    # Match common parameter assignment patterns
    patterns = [
        r"(\w+)\s*=\s*([\d.eE+-]+)\s*#\s*(.+)",       # x = 0.5  # description
        r"(\w+)\s*=\s*([\d.eE+-]+)\s*$",                # x = 0.5
        r"params\[.(\w+).\]\s*=\s*([\d.eE+-]+)",        # params["x"] = 0.5
    ]
    for pattern in patterns:
        for match in re.finditer(pattern, code, re.MULTILINE):
            groups = match.groups()
            name = groups[0]
            value = float(groups[1])
            desc = groups[2] if len(groups) > 2 else ""
            params.append({
                "name": name,
                "baseline": value,
                "description": desc.strip(),
                "range_low": value * 0.8,
                "range_high": value * 1.2,
            })
    return params

# Save scope
Path("artifacts").mkdir(exist_ok=True)
# params = extract_parameters_from_code("scripts/model.py")
# Path("artifacts/sensitivity_scope.json").write_text(
#     json.dumps(params, ensure_ascii=False, indent=2), encoding="utf-8"
# )

Phase 2: 单因素灵敏度分析 (OAT — One-At-a-Time)

Input: 参数列表、基准值、模型求解函数。 Output: artifacts/oat_results.json、龙卷风图。

1. 对每个参数，固定其他参数为基准值，只变动当前参数。
2. 在 [基准值 * (1 - PERTURBATION_RANGE), 基准值 * (1 + PERTURBATION_RANGE)] 范围内取 PERTURBATION_STEPS 个等距点。
3. 对每个点运行模型求解，记录目标函数值。
4. 计算灵敏度指数：S_i = (max(Y) - min(Y)) / Y_baseline。
5. 按灵敏度指数降序排列参数。
6. 灵敏度指数 > SIGNIFICANCE_THRESHOLD 的参数标记为"敏感参数"。
7. 对每个参数绘制"参数值 vs 目标函数值"曲线（蛛网图数据）。
8. 生成龙卷风图（Tornado Chart）：横轴为目标函数变化幅度，纵轴为各参数。
9. 龙卷风图按影响大小排序，最敏感的参数在最上方。
10. 如果某参数的影响明显非线性（非单调），在报告中特别标注。

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

def one_at_a_time(model_func, params: list, perturbation: float = 0.20, steps: int = 11):
    """Run OAT sensitivity analysis."""
    baseline_values = {p["name"]: p["baseline"] for p in params}
    y_baseline = model_func(**baseline_values)
    results = {}

    for p in params:
        name = p["name"]
        low = p["baseline"] * (1 - perturbation)
        high = p["baseline"] * (1 + perturbation)
        test_values = np.linspace(low, high, steps)
        y_values = []
        for v in test_values:
            trial = baseline_values.copy()
            trial[name] = v
            y_values.append(model_func(**trial))
        y_arr = np.array(y_values)
        sensitivity_index = (y_arr.max() - y_arr.min()) / abs(y_baseline) if y_baseline != 0 else float("inf")
        results[name] = {
            "test_values": test_values.tolist(),
            "y_values": [float(y) for y in y_arr],
            "sensitivity_index": round(float(sensitivity_index), 6),
            "y_baseline": float(y_baseline),
            "monotonic": bool(np.all(np.diff(y_arr) >= 0) or np.all(np.diff(y_arr) <= 0)),
        }
    return results

def plot_tornado(oat_results: dict, fig_path: str = "figures/sensitivity/tornado.png"):
    """Generate tornado chart from OAT results."""
    sorted_params = sorted(oat_results.items(), key=lambda x: x[1]["sensitivity_index"])
    names = [k for k, _ in sorted_params]
    low_deltas = []
    high_deltas = []
    baseline = list(oat_results.values())[0]["y_baseline"]

    for name, res in sorted_params:
        y_vals = res["y_values"]
        low_deltas.append(min(y_vals) - baseline)
        high_deltas.append(max(y_vals) - baseline)

    fig, ax = plt.subplots(figsize=(10, max(4, len(names) * 0.5)))
    y_pos = range(len(names))
    ax.barh(y_pos, high_deltas, left=0, color="#2196F3", alpha=0.8, label="Increase")
    ax.barh(y_pos, low_deltas, left=0, color="#FF5722", alpha=0.8, label="Decrease")
    ax.set_yticks(y_pos)
    ax.set_yticklabels(names)
    ax.axvline(x=0, color="black", linewidth=0.8)
    ax.set_xlabel("Change in Objective Function")
    ax.set_title("Tornado Chart — Parameter Sensitivity")
    ax.legend()
    fig.tight_layout()
    Path(fig_path).parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(fig_path, dpi=150)
    plt.close(fig)

Phase 3: 多因素灵敏度分析 (可选)

Input: Phase 2 识别出的敏感参数、模型求解函数。 Output: artifacts/multi_factor_results.json、交互效应热力图。

1. 仅对 Phase 2 中灵敏度指数排名前 3-5 的参数做多因素分析。
2. 使用拉丁超立方采样 (LHS) 或 Sobol 序列生成参数组合。
3. 计算 Sobol 全局灵敏度指数（一阶 S1、全阶 ST）。
4. S1 衡量单参数的独立贡献，ST 衡量包含交互效应的总贡献。
5. 如果 ST - S1 显著大于 0，说明该参数与其他参数存在交互效应。
6. 生成交互效应热力图，展示参数对之间的二阶灵敏度指数。
7. 如果计算量过大（预估 > MAX_COMPUTE_MINUTES），降低采样数或减少参数。
8. 对高度交互的参数对，在报告中建议后续建模关注其联合效应。
9. 如果所有参数的交互效应都很小，可以简化报告只展示一阶结果。
10. 此阶段为可选，如果时间紧迫可跳过并在报告中说明。

from scipy.stats import qmc

def latin_hypercube_sample(params: list, n_samples: int = 1000) -> np.ndarray:
    """Generate LHS samples for multi-factor analysis."""
    n_params = len(params)
    sampler = qmc.LatinHypercube(d=n_params)
    samples_unit = sampler.random(n=n_samples)
    lower = np.array([p["range_low"] for p in params])
    upper = np.array([p["range_high"] for p in params])
    samples = qmc.scale(samples_unit, lower, upper)
    return samples

def compute_sobol_indices(model_func, params: list, n_samples: int = 1024):
    """Compute first-order and total Sobol indices using Saltelli scheme."""
    from SALib.sample import saltelli
    from SALib.analyze import sobol as sobol_analyze

    problem = {
        "num_vars": len(params),
        "names": [p["name"] for p in params],
        "bounds": [[p["range_low"], p["range_high"]] for p in params],
    }
    X = saltelli.sample(problem, n_samples)
    Y = np.array([model_func(**dict(zip(problem["names"], x))) for x in X])
    Si = sobol_analyze.analyze(problem, Y)
    return {
        "S1": dict(zip(problem["names"], Si["S1"].tolist())),
        "ST": dict(zip(problem["names"], Si["ST"].tolist())),
        "S2": Si.get("S2", None),
    }

Phase 4: Monte Carlo 模拟

Input: 参数分布假设、模型求解函数。 Output: artifacts/monte_carlo_results.json、目标函数分布图。

1. 为每个敏感参数假设概率分布（正态、均匀、三角、对数正态等）。
2. 分布参数基于：数据拟合结果、物理约束、专家经验、文献参考。
3. 使用 MONTE_CARLO_SAMPLES 次随机采样。
4. 对每次采样运行模型求解，收集目标函数值序列。
5. 统计输出分布：均值、标准差、中位数、分位数、置信区间。
6. 计算 CONFIDENCE_LEVEL 置信区间的上下界。
7. 绘制目标函数的直方图+KDE 曲线。
8. 绘制累积分布函数 (CDF) 曲线。
9. 如果有多个目标函数，分别做 Monte Carlo 分析。
10. 评估结果的鲁棒性：置信区间宽度占均值的比例。

def monte_carlo_simulation(
    model_func,
    param_distributions: dict,
    n_samples: int = 10000,
    confidence: float = 0.95,
) -> dict:
    """
    Run Monte Carlo simulation.
    param_distributions: {"param_name": {"dist": "normal", "mean": 5.0, "std": 0.5}, ...}
    """
    samples = {}
    for name, dist_info in param_distributions.items():
        dist_type = dist_info["dist"]
        if dist_type == "normal":
            samples[name] = np.random.normal(dist_info["mean"], dist_info["std"], n_samples)
        elif dist_type == "uniform":
            samples[name] = np.random.uniform(dist_info["low"], dist_info["high"], n_samples)
        elif dist_type == "triangular":
            samples[name] = np.random.triangular(dist_info["left"], dist_info["mode"], dist_info["right"], n_samples)
        elif dist_type == "lognormal":
            samples[name] = np.random.lognormal(dist_info["mean"], dist_info["sigma"], n_samples)

    results = []
    for i in range(n_samples):
        kwargs = {name: float(samples[name][i]) for name in samples}
        results.append(model_func(**kwargs))

    results = np.array(results)
    alpha = 1 - confidence
    ci_low = np.percentile(results, 100 * alpha / 2)
    ci_high = np.percentile(results, 100 * (1 - alpha / 2))

    return {
        "n_samples": n_samples,
        "mean": float(np.mean(results)),
        "std": float(np.std(results)),
        "median": float(np.median(results)),
        "ci_low": float(ci_low),
        "ci_high": float(ci_high),
        "ci_level": confidence,
        "min": float(np.min(results)),
        "max": float(np.max(results)),
        "cv": float(np.std(results) / np.mean(results)) if np.mean(results) != 0 else float("inf"),
    }

def plot_mc_distribution(results: np.ndarray, fig_path: str = "figures/sensitivity/mc_distribution.png"):
    """Plot Monte Carlo result distribution."""
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Histogram + KDE
    axes[0].hist(results, bins=50, density=True, alpha=0.7, color="#2196F3", edgecolor="white")
    from scipy.stats import gaussian_kde
    kde = gaussian_kde(results)
    x_range = np.linspace(results.min(), results.max(), 200)
    axes[0].plot(x_range, kde(x_range), color="#FF5722", linewidth=2)
    axes[0].set_title("Monte Carlo Distribution (PDF)")
    axes[0].set_xlabel("Objective Function Value")
    axes[0].set_ylabel("Density")

    # CDF
    sorted_results = np.sort(results)
    cdf = np.arange(1, len(sorted_results) + 1) / len(sorted_results)