214 lines
6.5 KiB
Python
214 lines
6.5 KiB
Python
"""敏感性分析与鲁棒性检验。
|
||||
|
|
|
|||
|
|
★ 建议做成招牌模块。原建议这一节方向完全正确,只是函数名偏多,
|
|||
|
|
已合并同义项、删掉 `grid_sweep`(`parameter_sweep` 的特例)。
|
|||
|
|
|
|||
|
|
论文里的"敏感性分析"一节基本就是这四个函数的组合:
|
|||
|
|
局部扰动 → OAT 单因素扫描 → 多因素参数扫描 → 全局 Sobol/Morris。
|
|||
|
|
"""
|
|||
|
|
|
|||
|
|
from __future__ import annotations
|
|||
|
|
|
|||
|
|
from collections.abc import Callable
|
|||
|
|
from typing import Any
|
|||
|
|
|
|||
|
|
__all__ = [
|
|||
|
|
"local_sensitivity",
|
|||
|
|
"one_at_a_time",
|
|||
|
|
"parameter_sweep",
|
|||
|
|
"robustness_test",
|
|||
|
|
]
|
|||
|
|
|
|||
|
|
|
|||
|
|
def local_sensitivity(
|
|||
|
|
model: Callable[..., Any],
|
|||
|
|
params: Any,
|
|||
|
|
step: float = 1e-4,
|
|||
|
|
relative: bool = True,
|
|||
|
|
**kwargs: Any,
|
|||
|
|
) -> Any:
|
|||
|
|
"""局部敏感性:各参数的偏导 / 弹性。
|
|||
|
|
|
|||
|
|
Parameters
|
|||
|
|
----------
|
|||
|
|
relative:
|
|||
|
|
`True` 时算弹性 `(∂y/∂p)·(p/y)`——**无量纲,可以跨参数比较**,
|
|||
|
|
这是局部敏感性的正确做法;`False` 时是原始偏导(有量纲,只有同量纲参数之间可比)。
|
|||
|
|
|
|||
|
|
Returns
|
|||
|
|
-------
|
|||
|
|
`Result`:`data` 是参数 → 弹性 的表,按绝对值降序。
|
|||
|
|
|
|||
|
|
Notes
|
|||
|
|
-----
|
|||
|
|
有限差分的步长是精度与截断误差的权衡:太大截断误差大,
|
|||
|
|
太小被浮点误差淹没。`1e-4` 是稳健默认值;参数本身量级差很多时用相对步长。
|
|||
|
|
|
|||
|
|
Examples
|
|||
|
|
--------
|
|||
|
|
>>> s = local_sensitivity(lambda p: p["a"] ** 2 * p["b"], {"a": 2.0, "b": 3.0}) # doctest: +SKIP
|
|||
|
|
"""
|
|||
|
|
# TODO(L1) ★ S 档
|
|||
|
|
raise NotImplementedError
|
|||
|
|
|
|||
|
|
|
|||
|
|
def one_at_a_time(
|
|||
|
|
model: Callable[..., Any],
|
|||
|
|
params: Any,
|
|||
|
|
parameter: str | None = None,
|
|||
|
|
values: Any = None,
|
|||
|
|
n_points: int = 20,
|
|||
|
|
span: float = 0.5,
|
|||
|
|
) -> Any:
|
|||
|
|
"""OAT(一次只动一个参数)。
|
|||
|
|
|
|||
|
|
Parameters
|
|||
|
|
----------
|
|||
|
|
parameter:
|
|||
|
|
`None` 时对所有参数依次扫描。
|
|||
|
|
span:
|
|||
|
|
相对基准值的扫描范围,`0.5` 表示 ±50%。
|
|||
|
|
|
|||
|
|
Returns
|
|||
|
|
-------
|
|||
|
|
`Result`:`data` 是长表 `(parameter, value, output)`,配 `viz.plot` 直接出蛛网图/敏感性曲线。
|
|||
|
|
|
|||
|
|
Notes
|
|||
|
|
-----
|
|||
|
|
OAT 的**已知缺陷**:它看不见参数之间的交互作用。
|
|||
|
|
如果两个参数是相乘关系(很常见),OAT 会严重低估整体敏感性——
|
|||
|
|
论文里必须说明这一点,并补一个全局敏感性作为交叉验证。
|
|||
|
|
"""
|
|||
|
|
# TODO(L1)
|
|||
|
|
raise NotImplementedError
|
|||
|
|
|
|||
|
|
|
|||
|
|
def parameter_sweep(
|
|||
|
|
model: Callable[..., Any],
|
|||
|
|
ranges: dict[str, Any],
|
|||
|
|
n_points: int = 20,
|
|||
|
|
method: str = "grid",
|
|||
|
|
n_samples: int = 1000,
|
|||
|
|
random_state: int | None = None,
|
|||
|
|
) -> Any:
|
|||
|
|
"""多参数扫描(合并了原建议的 `grid_sweep` 与 `random_parameter_sweep`)。
|
|||
|
|
|
|||
|
|
Parameters
|
|||
|
|
----------
|
|||
|
|
ranges:
|
|||
|
|
`{"alpha": (0.1, 0.9), "beta": (0, 5)}`。
|
|||
|
|
method:
|
|||
|
|
`"grid"`:笛卡尔积,点数 = `n_points^d`,**维度一高就爆**(d=5, n=20 是 320 万次);
|
|||
|
|
`"lhs"`:拉丁超立方,维度高时用这个;
|
|||
|
|
`"random"`:纯随机。
|
|||
|
|
|
|||
|
|
Returns
|
|||
|
|
-------
|
|||
|
|
`Result`:`data` 是 `(n_points^d, d + 1)` 的结果表,
|
|||
|
|
`params["best"]` 是让目标最优的参数组合。
|
|||
|
|
|
|||
|
|
Notes
|
|||
|
|
-----
|
|||
|
|
`d >= 4` 时默认改走 `"lhs"`,否则 `n_points` 会被指数放大到跑不完。
|
|||
|
|
|
|||
|
|
Examples
|
|||
|
|
--------
|
|||
|
|
>>> r = parameter_sweep(lambda a, b: -(a - 0.3) ** 2 - (b - 2) ** 2,
|
|||
|
|
... {"a": (0, 1), "b": (0, 4)}, n_points=11) # doctest: +SKIP
|
|||
|
|
"""
|
|||
|
|
# TODO(L0)
|
|||
|
|
raise NotImplementedError
|
|||
|
|
|
|||
|
|
|
|||
|
|
def robustness_test(
|
|||
|
|
model: Callable[..., Any],
|
|||
|
|
scenarios: Any = None,
|
|||
|
|
params: Any = None,
|
|||
|
|
perturbation: float = 0.1,
|
|||
|
|
n_samples: int = 1000,
|
|||
|
|
random_state: int | None = None,
|
|||
|
|
) -> Any:
|
|||
|
|
"""鲁棒性检验(合并了原建议的 `perturbation_test`)。
|
|||
|
|
|
|||
|
|
在参数上施加 ±`perturbation` 的随机扰动,看结论是否翻盘。
|
|||
|
|
|
|||
|
|
Returns
|
|||
|
|
-------
|
|||
|
|
`Result`:`metrics` 含 `"rank_stability"`(原结论在前 X% 的比例)、
|
|||
|
|
`"conclusion_flip_rate"`(结论反转频率)。**这两个数字才是审稿人想看的。**
|
|||
|
|
|
|||
|
|
Notes
|
|||
|
|
-----
|
|||
|
|
鲁棒性检验和敏感性分析不是一回事:
|
|||
|
|
敏感性问"参数变了输出变多少",鲁棒性问"参数变了**结论**还成立吗"。
|
|||
|
|
论文里两节都要有,别混着写。
|
|||
|
|
|
|||
|
|
Examples
|
|||
|
|
--------
|
|||
|
|
>>> r = robustness_test(model, params={"alpha": 0.5}, random_state=0) # doctest: +SKIP
|
|||
|
|
"""
|
|||
|
|
# TODO(L1)
|
|||
|
|
raise NotImplementedError
|
|||
|
|
|
|||
|
|
|
|||
|
|
def sobol_analyze(
|
|||
|
|
model: Callable[..., Any],
|
|||
|
|
bounds: Any,
|
|||
|
|
n_samples: int = 1024,
|
|||
|
|
calc_second_order: bool = False,
|
|||
|
|
random_state: int | None = None,
|
|||
|
|
) -> Any:
|
|||
|
|
"""Sobol 全局敏感性(底层 SALib,自己不要重写)。
|
|||
|
|
|
|||
|
|
Parameters
|
|||
|
|
----------
|
|||
|
|
bounds:
|
|||
|
|
`{"alpha": (0, 1), "beta": (0, 5)}` 形式的参数上下界。
|
|||
|
|
n_samples:
|
|||
|
|
SALib 的基样本数 N;**实际模型调用次数是 `N·(2d + 2)`**
|
|||
|
|
(二阶时更多),d=5、N=1024 就是 12288 次调用。先估一下模型单次耗时。
|
|||
|
|
|
|||
|
|
Returns
|
|||
|
|
-------
|
|||
|
|
`Result`:`data` 含 `S1`(一阶)、`ST`(总效应)、`S2`(二阶,可选)三列。
|
|||
|
|
|
|||
|
|
Notes
|
|||
|
|
-----
|
|||
|
|
* `ST - S1` 的差距就是**交互作用**的大小——这正是 OAT 看不见的东西。
|
|||
|
|
两者差得多,说明不能只看单因素。
|
|||
|
|
* 需要 `pip install "mcm-toolkit[sensitivity]"`;缺 SALib 时本函数给出
|
|||
|
|
可读的安装提示,而不是裸抛 `ModuleNotFoundError`。
|
|||
|
|
* 费用高:**先做 OAT 缩小参数范围,再对少数关键参数做 Sobol**,
|
|||
|
|
不要一上来就对全部参数做全局敏感性。
|
|||
|
|
|
|||
|
|
Examples
|
|||
|
|
--------
|
|||
|
|
>>> r = sobol_analyze(model, bounds={"a": (0, 1), "b": (0, 5)}) # doctest: +SKIP
|
|||
|
|
"""
|
|||
|
|
# TODO(L1) ★ S 档:配完整函数卡片
|
|||
|
|
raise NotImplementedError
|
|||
|
|
|
|||
|
|
|
|||
|
|
def morris_analyze(
|
|||
|
|
model: Callable[..., Any],
|
|||
|
|
bounds: Any,
|
|||
|
|
n_trajectories: int = 10,
|
|||
|
|
n_levels: int = 4,
|
|||
|
|
random_state: int | None = None,
|
|||
|
|
) -> Any:
|
|||
|
|
"""Morris 筛选法(底层 SALib)。
|
|||
|
|
|
|||
|
|
调用次数只与轨迹数成正比(`(d + 1) · n_trajectories`),比 Sobol 便宜一个量级,
|
|||
|
|
适合**参数多但只需要找出"哪几个重要"**的筛选阶段。
|
|||
|
|
|
|||
|
|
Returns
|
|||
|
|
-------
|
|||
|
|
`Result`:`data` 含 `mu`、`mu_star`(推荐看的)、`sigma` 三列。
|
|||
|
|
`sigma` 大表示该参数的效应是非线性的或与其他参数有交互。
|
|||
|
|
|
|||
|
|
Examples
|
|||
|
|
--------
|
|||
|
|
>>> r = morris_analyze(model, bounds={"a": (0, 1), "b": (0, 5)}) # doctest: +SKIP
|
|||
|
|
"""
|
|||
|
|
# TODO(L1)
|
|||
|
|
raise NotImplementedError
|