scipy.stats.

sobol_indices

scipy.stats.sobol_indices(*, func, n, dists=None, method='saltelli_2010', random_state=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), optional

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

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

如果 func 是一个可调用的对象，则必须指定，否则忽略。

methodCallable 或 str，默认值：‘saltelli_2010’

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

如果是一个可调用的对象，则其签名必须是

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

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

random_state{None, int, numpy.random.Generator}, optional

如果random_state是整数或 None，将使用np.random.default_rng(random_state)创建一个新的numpy.random.Generator。如果random_state已经是Generator实例，则使用所提供的实例。

返回 :

resSobolResult

具有以下属性的对象

first_order形状为 (s, d) 的 ndarray

一阶 Sobol’ 指数。

total_order形状为 (s, d) 的 ndarray

总阶 Sobol’ 指数。

以及方法

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

提供对指数的置信区间的函数。有关详细信息，请参阅scipy.stats.bootstrap。

自举是在一阶和总阶指数上进行的，可以在BootstrapSobolResult中以 first_order 和 total_order 形式作为属性获得。

注意

Sobol’ 方法[1]、[2]是一种基于方差的灵敏度分析，它可以获得每个参数对目标量（即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.. “非线性数学模型的敏感性分析。”Mathematical Modeling and Computational Experiment，1:407-414，1993。

[2]

Sobol, I. M. (2001)。“非线性数学模型及其蒙特卡罗估计的全局敏感性指标。”Mathematics and Computers in Simulation，55(1-3):271-280，DOI:10.1016/S0378-4754(00)00270-6，2001。

[3]

Saltelli, A. “最有效地利用模型评估计算敏感性指标。”Computer Physics Communications，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。 “基于方差的模型输出敏感性分析。总敏感性指标的设计和估计。”Computer Physics Communications，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)
...     ],
...     random_state=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 指数的优点。

返回在新标签页中打开