scipy.stats.

sobol_indices

scipy.stats.sobol_indices(*, func, n, dists=None, method='saltelli_2010', rng=None)

Sobol’ 的全局敏感性指数。

参数:

func可调用对象或 dict(str, array_like)

如果 func 是可调用对象，则为从中计算 Sobol’ 指数的函数。其签名必须是

func(x: ArrayLike) -> ArrayLike

其中 x 的形状为 (d, n)，输出的形状为 (s, n)，其中

• d 是 func 的输入维度（输入变量的数量），

• s 是 func 的输出维度（输出变量的数量），并且

• n 是样本数量（参见下面的 n）。

函数评估值必须是有限的。

如果 func 是字典，则包含来自三个不同数组的函数评估值。键必须是：f_A、f_B 和 f_AB。f_A 和 f_B 的形状应为 (s, n)，f_AB 的形状应为 (d, s, n)。这是一个高级功能，误用可能会导致错误的分析。

nint

用于生成矩阵 A 和 B 的样本数量。必须是 2 的幂。评估 func 的总点数将为 n*(d+2)。

distslist(distributions), 可选

每个参数分布的列表。参数的分布取决于应用程序，应仔细选择。假设参数是独立分布的，这意味着它们的值之间没有约束或关系。

分布必须是具有 ppf 方法的类的实例。

如果 func 是可调用对象，则必须指定，否则将被忽略。

method可调用对象或 str，默认值：‘saltelli_2010’

用于计算第一和总 Sobol’ 指数的方法。

如果为可调用对象，则其签名必须是

func(f_A: np.ndarray, f_B: np.ndarray, f_AB: np.ndarray)
-> Tuple[np.ndarray, np.ndarray]

其中 f_A, f_B 的形状为 (s, n)，f_AB 的形状为 (d, s, n)。这些数组包含来自三组不同样本的函数评估值。输出是第一和总指数的元组，形状为 (s, d)。这是一个高级功能，误用可能会导致错误的分析。

rngnumpy.random.Generator, 可选

伪随机数生成器状态。当 rng 为 None 时，将使用来自操作系统的熵创建一个新的 numpy.random.Generator。其他类型将传递给 numpy.random.default_rng 以实例化 Generator。

在 1.15.0 版本中更改: 作为从使用 numpy.random.RandomState 转换为 numpy.random.Generator 的 SPEC-007 过渡的一部分，此关键字已从 random_state 更改为 rng。在过渡期间，这两个关键字将继续工作，尽管一次只能指定一个。在过渡期之后，使用 random_state 关键字的函数调用将发出警告。经过弃用期后，将删除 random_state 关键字。

返回:

resSobolResult

具有以下属性的对象

first_orderndarray，形状 (s, d)

一阶 Sobol’ 指数。

total_orderndarray，形状 (s, d)

总阶 Sobol’ 指数。

和方法

bootstrap(confidence_level: float, n_resamples: int) -> BootstrapSobolResult

一种在指数上提供置信区间的方法。有关更多详细信息，请参阅 scipy.stats.bootstrap。

自举是在一阶和总阶指数上完成的，并且它们在 BootstrapSobolResult 中以属性 first_order 和 total_order 的形式提供。

注释

Sobol’ 方法 [1], [2] 是一种基于方差的敏感性分析，它获得每个参数对感兴趣的量（QoI；即 func 的输出）的方差的贡献。各个贡献可用于对参数进行排序，还可以通过计算模型的有效（或平均）维度来衡量模型的复杂性。

注意

假设参数是独立分布的。每个参数仍然可以遵循任何分布。事实上，分布非常重要，应该与参数的实际分布相匹配。

它使用函数的方差的函数分解来探索

\[\mathbb{V}(Y) = \sum_{i}^{d} \mathbb{V}_i (Y) + \sum_{i<j}^{d} \mathbb{V}_{ij}(Y) + ... + \mathbb{V}_{1,2,...,d}(Y),\]

引入条件方差

\[\mathbb{V}_i(Y) = \mathbb{\mathbb{V}}[\mathbb{E}(Y|x_i)] \qquad \mathbb{V}_{ij}(Y) = \mathbb{\mathbb{V}}[\mathbb{E}(Y|x_i x_j)] - \mathbb{V}_i(Y) - \mathbb{V}_j(Y),\]

Sobol’ 指数表示为

\[S_i = \frac{\mathbb{V}_i(Y)}{\mathbb{V}[Y]} \qquad S_{ij} =\frac{\mathbb{V}_{ij}(Y)}{\mathbb{V}[Y]}.\]

\(S_{i}\) 对应于一阶项，它评估第 i 个参数的贡献，而 \(S_{ij}\) 对应于二阶项，它告知第 i 个和第 j 个参数之间交互作用的贡献。这些方程可以推广到计算更高阶项；但是，它们的计算成本很高，并且它们的解释很复杂。这就是为什么只提供一阶指数的原因。

总阶指数表示参数对 QoI 方差的全局贡献，并定义为

\[S_{T_i} = S_i + \sum_j S_{ij} + \sum_{j,k} S_{ijk} + ... = 1 - \frac{\mathbb{V}[\mathbb{E}(Y|x_{\sim i})]}{\mathbb{V}[Y]}.\]

一阶指数的总和最多为 1，而总阶指数的总和至少为 1。如果没有交互作用，则一阶和总阶指数相等，并且一阶和总阶指数的总和都为 1。

警告

负 Sobol’ 值是由于数值误差造成的。增加点数 n 应该会有所帮助。

进行良好分析所需的样本数量随着问题的维度而增加。例如，对于 3 维问题，请考虑至少 n >= 2**12。模型越复杂，需要的样本就越多。

即使对于纯粹的加性模型，由于数值噪声，指数的总和也可能不等于 1。

参考文献

[1]

Sobol, I. M.. “非线性数学模型的灵敏度分析。” 数学建模与计算实验, 1:407-414, 1993.

[2]

Sobol, I. M. (2001). “非线性数学模型的全局灵敏度指标及其蒙特卡洛估计。” 数学与计算机模拟, 55(1-3):271-280, DOI:10.1016/S0378-4754(00)00270-6, 2001.

[3]

Saltelli, A. “充分利用模型评估来计算灵敏度指标。” 计算机物理通讯, 145(2):280-297, DOI:10.1016/S0010-4655(02)00280-1, 2002.

[4]

Saltelli, A., M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, 和 S. Tarantola. “全局灵敏度分析：入门。” 2007.

[5]

Saltelli, A., P. Annoni, I. Azzini, F. Campolongo, M. Ratto, 和 S. Tarantola. “模型输出的基于方差的灵敏度分析。总体灵敏度指标的设计和估计器。” 计算机物理通讯, 181(2):259-270, DOI:10.1016/j.cpc.2009.09.018, 2010.

[6]

Ishigami, T. 和 T. Homma. “计算机模型不确定性分析中的重要性量化技术。” IEEE, DOI:10.1109/ISUMA.1990.151285, 1990.

示例

在您的浏览器中尝试一下！

以下是 Ishigami 函数的示例 [6]

\[Y(\mathbf{x}) = \sin x_1 + 7 \sin^2 x_2 + 0.1 x_3^4 \sin x_1,\]

其中 \(\mathbf{x} \in [-\pi, \pi]^3\)。此函数表现出强烈的非线性和非单调性。

请记住，Sobol' 指标假设样本是独立分布的。在这种情况下，我们在每个边缘上使用均匀分布。

>>> import numpy as np
>>> from scipy.stats import sobol_indices, uniform
>>> rng = np.random.default_rng()
>>> def f_ishigami(x):
...     f_eval = (
...         np.sin(x[0])
...         + 7 * np.sin(x[1])**2
...         + 0.1 * (x[2]**4) * np.sin(x[0])
...     )
...     return f_eval
>>> indices = sobol_indices(
...     func=f_ishigami, n=1024,
...     dists=[
...         uniform(loc=-np.pi, scale=2*np.pi),
...         uniform(loc=-np.pi, scale=2*np.pi),
...         uniform(loc=-np.pi, scale=2*np.pi)
...     ],
...     rng=rng
... )
>>> indices.first_order
array([0.31637954, 0.43781162, 0.00318825])
>>> indices.total_order
array([0.56122127, 0.44287857, 0.24229595])

可以使用自举法获得置信区间。

>>> boot = indices.bootstrap()

然后，可以很容易地可视化这些信息。

>>> import matplotlib.pyplot as plt
>>> fig, axs = plt.subplots(1, 2, figsize=(9, 4))
>>> _ = axs[0].errorbar(
...     [1, 2, 3], indices.first_order, fmt='o',
...     yerr=[
...         indices.first_order - boot.first_order.confidence_interval.low,
...         boot.first_order.confidence_interval.high - indices.first_order
...     ],
... )
>>> axs[0].set_ylabel("First order Sobol' indices")
>>> axs[0].set_xlabel('Input parameters')
>>> axs[0].set_xticks([1, 2, 3])
>>> _ = axs[1].errorbar(
...     [1, 2, 3], indices.total_order, fmt='o',
...     yerr=[
...         indices.total_order - boot.total_order.confidence_interval.low,
...         boot.total_order.confidence_interval.high - indices.total_order
...     ],
... )
>>> axs[1].set_ylabel("Total order Sobol' indices")
>>> axs[1].set_xlabel('Input parameters')
>>> axs[1].set_xticks([1, 2, 3])
>>> plt.tight_layout()
>>> plt.show()

注意

默认情况下，scipy.stats.uniform 的支持为 [0, 1]。使用参数 loc 和 scale，可以获得 [loc, loc + scale] 上的均匀分布。

这个结果特别有趣，因为一阶指标 \(S_{x_3} = 0\)，而它的总阶指标为 \(S_{T_{x_3}} = 0.244\)。这意味着与 \(x_3\) 的高阶相互作用导致了差异。QoI 上观察到的方差几乎有 25% 是由于 \(x_3\) 和 \(x_1\) 之间的相关性造成的，尽管 \(x_3\) 本身对 QoI 没有影响。

下面给出了此函数上 Sobol’ 指标的可视化解释。让我们在 \([-\pi, \pi]^3\) 中生成 1024 个样本并计算输出值。

>>> from scipy.stats import qmc
>>> n_dim = 3
>>> p_labels = ['$x_1$', '$x_2$', '$x_3$']
>>> sample = qmc.Sobol(d=n_dim, seed=rng).random(1024)
>>> sample = qmc.scale(
...     sample=sample,
...     l_bounds=[-np.pi, -np.pi, -np.pi],
...     u_bounds=[np.pi, np.pi, np.pi]
... )
>>> output = f_ishigami(sample.T)

现在我们可以绘制输出相对于每个参数的散点图。这提供了一种可视化方式来理解每个参数如何影响函数的输出。

>>> fig, ax = plt.subplots(1, n_dim, figsize=(12, 4))
>>> for i in range(n_dim):
...     xi = sample[:, i]
...     ax[i].scatter(xi, output, marker='+')
...     ax[i].set_xlabel(p_labels[i])
>>> ax[0].set_ylabel('Y')
>>> plt.tight_layout()
>>> plt.show()

现在，Sobol’ 更进一步：通过给定参数值（黑线）对输出值进行条件化，计算条件输出均值。它对应于项 \(\mathbb{E}(Y|x_i)\)。取该项的方差得到 Sobol' 指标的分子。

>>> mini = np.min(output)
>>> maxi = np.max(output)
>>> n_bins = 10
>>> bins = np.linspace(-np.pi, np.pi, num=n_bins, endpoint=False)
>>> dx = bins[1] - bins[0]
>>> fig, ax = plt.subplots(1, n_dim, figsize=(12, 4))
>>> for i in range(n_dim):
...     xi = sample[:, i]
...     ax[i].scatter(xi, output, marker='+')
...     ax[i].set_xlabel(p_labels[i])
...     for bin_ in bins:
...         idx = np.where((bin_ <= xi) & (xi <= bin_ + dx))
...         xi_ = xi[idx]
...         y_ = output[idx]
...         ave_y_ = np.mean(y_)
...         ax[i].plot([bin_ + dx/2] * 2, [mini, maxi], c='k')
...         ax[i].scatter(bin_ + dx/2, ave_y_, c='r')
>>> ax[0].set_ylabel('Y')
>>> plt.tight_layout()
>>> plt.show()

查看 \(x_3\)，均值的方差为零，导致 \(S_{x_3} = 0\)。但是我们可以进一步观察到，输出的方差在 \(x_3\) 的参数值下不是恒定的。这种异方差性是由高阶相互作用解释的。此外，在 \(x_1\) 上也注意到异方差性，导致 \(x_3\) 和 \(x_1\) 之间存在相互作用。在 \(x_2\) 上，方差似乎是恒定的，因此可以假设与该参数没有相互作用。

这种情况很容易进行可视化分析——尽管这只是一种定性分析。然而，当输入参数的数量增加时，这种分析变得不现实，因为很难得出关于高阶项的结论。因此，使用 Sobol’ 指标的好处就体现出来了。

返回在新标签页中打开