sensitivity-analysis

# 灵敏度分析 执行描述: $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. 确认模型可以在修改参数后重新求解（函数化封装）。 ```python 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. 如果某参数的影响明显非线性（非单调），在报告中特别标注。 ```python 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. 此阶段为可选，如果时间紧迫可跳过并在报告中说明。 ```python 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. 评估结果的鲁棒性：置信区间宽度占均值的比例。 ```python 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)