索博尔指数

索博尔指数（Sobol index）是以输出方差为基准的全局敏感性分析指标，用来区分输入的单独效应与相互作用。一阶指数度量只给定某个输入时，条件均值能够体现多少输出方差；总效应指数把所有涉及该输入的相互作用也包括在内。经典分解以输入相互独立、输出平方可积且方差非零为前提。^[1]

一个含相互作用的模型

设 $X, Y$ 相互独立，且都服从区间 $[- 1, 1]$ 上的均匀分布。取无量纲输出 $F = f (X, Y) = X + 2 Y + X Y .$ 这是一组为手算设置的随机输入，不是观测数据。输入值本身只是模型的一部分，分布同样重要：若改变取值范围或概率，后面的方差贡献也会改变。

三项的含义可以直接分辨。 $X$ 与 $2 Y$ 分别随单个输入变化； $X Y$ 使一个输入的效果随另一个输入改变。例如 $Y = 1$ 时， $X$ 的系数为二； $Y = - 1$ 时， $X$ 的系数为零。接下来要把这些作用变成可以相加的方差份额，而不是只看某个点上的导数。

条件均值与正交分解

先对一般的两输入模型作分解。设 $X, Y$ 相互独立，且 $𝔼 [F^{2}] < \infty$ ，定义 $f_{0} = 𝔼 [F], f_{1} (X) = 𝔼 [F ∣ X] - f_{0}, f_{2} (Y) = 𝔼 [F ∣ Y] - f_{0},$ 以及剩余项 $f_{12} (X, Y) = F - f_{0} - f_{1} (X) - f_{2} (Y) .$ 条件均值 $𝔼 [F ∣ X = x]$ 的意思是固定 $x$ ，按照 $Y$ 的概率规律平均其他输入。 $f_{1}$ 是该均值相对总体均值的偏移， $f_{12}$ 则保留不能写成两个单变量效应之和的部分。这种分解也称功能方差分析分解。

由条件期望的迭代性质， $𝔼 [f_{1}] = 𝔼 [f_{2}] = 0$ 。独立性还给出 $𝔼 [f_{2} (Y) ∣ X] = 𝔼 [f_{2} (Y)] = 0,$ 所以 $𝔼 [f_{12} ∣ X] = 𝔼 [F ∣ X] - f_{0} - f_{1} (X) - 0 = 0.$ 同理 $𝔼 [f_{12} ∣ Y] = 0$ 。这一步使各部分正交： $𝔼 [f_{1} f_{2}] = 𝔼 [f_{1}] 𝔼 [f_{2}] = 0, 𝔼 [f_{1} f_{12}] = 𝔼 [f_{1} 𝔼 [f_{12} ∣ X]] = 0,$ 而 $𝔼 [f_{2} f_{12}] = 0$ 由相同理由成立。因此将 $F - f_{0} = f_{1} + f_{2} + f_{12}$ 平方、取期望时，交叉项全部为零，得到 $V = Var (F) = V_{1} + V_{2} + V_{12}, V_{1} = 𝔼 [f_{1}^{2}], V_{2} = 𝔼 [f_{2}^{2}], V_{12} = 𝔼 [f_{12}^{2}] .$ 三项均非负。加法不是因为各项取不同名字，而是因为已证明它们的交叉内积为零。多输入情形进一步对输入子集逐级分解，包含三阶及更高阶相互作用；两输入模型已能完整展示基本机制。^[1]

逐项算出方差份额

回到主例。均匀分布的对称性给出 $𝔼 [X] = 𝔼 [Y] = 0$ ，并且 $𝔼 [X^{2}] = \int_{- 1}^{1} \frac{x^{2}}{2} d x = \frac{1}{3}, 𝔼 [Y^{2}] = \frac{1}{3} .$ 由独立性， $𝔼 [X Y] = 0$ ，所以 $f_{0} = 0$ 。固定 $X$ 平均 $Y$ ，得 $f_{1} = X$ ；固定 $Y$ 平均 $X$ ，得 $f_{2} = 2 Y$ ；剩余项为 $f_{12} = X Y$ 。于是 $V_{1} = \frac{1}{3}, V_{2} = 4 \cdot \frac{1}{3} = \frac{4}{3}, V_{12} = 𝔼 [X^{2}] 𝔼 [Y^{2}] = \frac{1}{9}, V = \frac{16}{9} .$ 最后一个乘积分解再次使用了独立性。图中把总方差分成十六等份，三份来自 $X$ ，十二份来自 $2 Y$ ，一份来自相互作用。

总方差16除以9分成三段，X单独贡献3除以9，Y单独贡献12除以9，交互贡献1除以9 — 长度表示方差，不表示函数值。三个正交部分的方差相加为总方差。

一阶指数与总效应指数

当 $V > 0$ 时，一阶指数定义为 $S_{X} = \frac{Var (𝔼 [F ∣ X])}{V} = \frac{V_{1}}{V}, S_{Y} = \frac{V_{2}}{V} .$ 两输入相互作用指数为 $S_{X Y} = V_{12} / V$ 。主例中， $S_{X} = \frac{3}{16}, S_{Y} = \frac{12}{16} = \frac{3}{4}, S_{X Y} = \frac{1}{16} .$ 一阶指数之和为 $15 / 16$ ，余下的 $1 / 16$ 是双方相互作用，不能硬塞给某一输入的单独效应。

总效应指数采用另一个条件平均次序：固定其他输入，保留当前输入的变化，再平均条件方差。两输入时定义为 $S_{T_{X}} = \frac{𝔼 [Var (F ∣ Y)]}{V} = 1 - \frac{Var (𝔼 [F ∣ Y])}{V},$ 其中等号由全方差公式得到。代入正交分解，得 $S_{T_{X}} = (V_{1} + V_{12}) / V$ ，同理 $S_{T_{Y}} = (V_{2} + V_{12}) / V$ 。因此 $S_{T_{X}} = \frac{1}{4}, S_{T_{Y}} = \frac{13}{16} .$ 总效应相加为 $17 / 16$ ，超过一，是因为相互作用在两个输入各自的总效应中都被计入。它不是可再次相加的互斥分配。理论上，每个一阶或总效应指数都在零与一之间；这由方差分解和全方差公式保证。^[2]

X的一阶与总效应分别为3除以16和4除以16，Y分别为12除以16和13除以16，两组各相差1除以16 — 每组柱子的差都是本例的相互作用份额。总效应包括一阶效应。

一阶指数为零不等于输入无作用。将模型改为 $F = X Y$ ，仍用原来的独立均匀输入，则 $𝔼 [F ∣ X] = 𝔼 [F ∣ Y] = 0$ ，故两一阶指数均为零；全部方差来自相互作用，两总效应指数均为一。只有查看单独和总效应，才能看到这个差别。

配对替换的蒙特卡洛估计

高维模型通常无法手算积分，可以通过重复模型运行估计指数。以下推导同时说明替换输入时哪些坐标必须保持相同。

令 $A = (A_{X}, A_{Y})$ 、 $B = (B_{X}, B_{Y})$ 是独立同分布的两组输入，组内也按主例独立抽样。构造 $C_{X} = (B_{X}, A_{Y})$ ：它与 B 共享 X，与 A 共享其余输入。给定 $B_{X}$ 时， $f (B)$ 与 $f (C_{X})$ 的其他坐标独立，故条件乘积的均值为 $𝔼 [F ∣ X = B_{X}]^{2}$ 。另一方面 A 与 B 独立， $𝔼 [f (A) f (B)] = f_{0}^{2}$ 。相减得 $𝔼 [f (B) (f (C_{X}) - f (A))] = V_{1} .$ 对于总效应，给定共同的 $A_{Y}$ 后， $f (A)$ 与 $f (C_{X})$ 是条件独立同分布的两个输出。同分布独立副本的差平方期望是两倍方差，因此 $\frac{1}{2} 𝔼 [(f (A) - f (C_{X}))^{2}] = 𝔼 [Var (F ∣ Y)] .$ 将这两个期望分别用样本均值替代，再除以输出方差，就得到一阶与总效应的估计。^[3]

输入向量A和B交叉组合为CX与CY，CX取B的X和A的Y，CY取A的X和B的Y — 每组样本运行 A、B、CX、CY 四个点。共享哪一个坐标，决定估计的是哪个条件方差。

下面的代码专用于主例，用已经证明的精确方差 $16 / 9$ 作分母，便于把抽样误差与分母估计误差分开。一般模型还需估计方差，并检查稳定性和置信区间；不能把下面的已知常数直接带到其他模型。

Python 3

import random

def estimate_example(n=100000, seed=27):
    if isinstance(n, bool) or not isinstance(n, int) or n < 1:
        raise ValueError("n must be a positive integer")
    rng = random.Random(seed)
    first_x = first_y = total_x = total_y = 0.0

    def model(x, y):
        return x + 2.0 * y + x * y

    for _ in range(n):
        ax, ay = rng.uniform(-1, 1), rng.uniform(-1, 1)
        bx, by = rng.uniform(-1, 1), rng.uniform(-1, 1)
        fa, fb = model(ax, ay), model(bx, by)
        cx, cy = model(bx, ay), model(ax, by)
        first_x += fb * (cx - fa)
        first_y += fb * (cy - fa)
        total_x += 0.5 * (fa - cx) ** 2
        total_y += 0.5 * (fa - cy) ** 2

    denominator = n * (16.0 / 9.0)
    return {
        "S_X": first_x / denominator,
        "S_Y": first_y / denominator,
        "ST_X": total_x / denominator,
        "ST_Y": total_y / denominator,
    }

print(estimate_example())
print({"S_X": 3 / 16, "S_Y": 12 / 16,
       "ST_X": 4 / 16, "ST_Y": 13 / 16})

输入 n 为正整数，seed 控制重复实验；返回四个估计，不改写外部输入。每组四次模型计算，总时间为 $O (n)$ ，额外空间为 $O (1)$ 。四个累加器分别对应刚证明的期望恒等式，因此随着独立样本组数增加，由大数定律趋向相应指数。有限样本的一阶估计可能为负，或估计值超过一；这不改变理论界限，也不应通过直接裁剪掩盖估计误差。两次运行看起来接近，仍不能替代误差评估。

独立性与解释边界

独立性用于正交分解，也用于构造共享部分坐标的独立副本。若 $X = Y$ 、 $F = X$ 且方差非零，那么给定 X 或给定 Y 都完全确定 F，两个条件均值的方差比都为一。它们不能被解释成互斥的两份方差贡献。依赖输入需要另行规定共享信息的归属或采用适合相关性的分析，不能照用独立输入的加法结论。

若输出恒为常数，总方差为零，所有上述比值均未定义；不应把它们报告为零。若输入不确定性很大但输出近乎常数，估计时也可能受数值误差影响。即使定义成立，索博尔指数描述的仍是指定模型和输入概率下的方差贡献，不表示控制某输入的实际干预效果，也不保证对尾部事件最有解释力。

作为可复算的变式，将主例去掉乘积项，得到 $F = X + 2 Y$ 。此时总方差为 $5 / 3$ ，两一阶指数为 $1 / 5$ 、 $4 / 5$ ，各自总效应与一阶相等。比较它与主例，能够看出加入相互作用后，分子和总方差分母都会改变。

参考文献

↑ ^1.0 ^1.1 I. Azzini、G. Listorti、T. A. Mara、R. Rosati，2020，Uncertainty and Sensitivity Analysis for policy decision making: An introductory guide，Joint Research Centre，EUR 30432 EN，关于方差分解和全局敏感度的章节。
↑ SALib，Basics of sensitivity analysis，对一阶、二阶及总效应的解释；实现文档见 SALib.analyze.sobol.analyze。
↑ A. Saltelli、P. Annoni、I. Azzini、F. Campolongo、M. Ratto、S. Tarantola，2010，Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index，Computer Physics Communications，181(2):259—270，doi:10.1016/j.cpc.2009.09.018。公式实现及该出处可在 SALib 官方 Sobol 源码的 first_order、total_order 中核对；本文通过条件期望直接推导所用恒等式。

[jrcguide-1] 1.0 ^1.1 I. Azzini、G. Listorti、T. A. Mara、R. Rosati，2020，Uncertainty and Sensitivity Analysis for policy decision making: An introductory guide，Joint Research Centre，EUR 30432 EN，关于方差分解和全局敏感度的章节。

[2] SALib，Basics of sensitivity analysis，对一阶、二阶及总效应的解释；实现文档见 SALib.analyze.sobol.analyze。

[3] A. Saltelli、P. Annoni、I. Azzini、F. Campolongo、M. Ratto、S. Tarantola，2010，Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index，Computer Physics Communications，181(2):259—270，doi:10.1016/j.cpc.2009.09.018。公式实现及该出处可在 SALib 官方 Sobol 源码的 first_order、total_order 中核对；本文通过条件期望直接推导所用恒等式。

[1]

[2]

[3]