scipy.stats.

sobol_indices#

scipy.stats.sobol_indices(*, func, n, dists=None, method='saltelli_2010', random_state=None)[源代码][源代码]#

Sobol’ 的全局敏感性指数

参数:
函数可调用对象或字典(str, array_like)

如果 func 是一个可调用的函数,用于计算 Sobol’ 指数。其签名必须是:

func(x: ArrayLike) -> ArrayLike

对于形状为 (d, n)x 和形状为 (s, n) 的输出,其中:

  • dfunc 的输入维度(输入变量的数量),

  • sfunc 的输出维度(输出变量的数量),并且

  • n 是样本的数量(见下面的 n)。

函数求值结果必须是有限的。

如果 func 是一个字典,包含来自三个不同数组的函数评估。键必须是:f_Af_Bf_ABf_Af_B 应该具有形状 (s, n),而 f_AB 应该具有形状 (d, s, n)。这是一个高级功能,误用可能导致错误的分析。

n整数

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

distslist(分布), 可选

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

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

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

方法可调用对象或字符串,默认值:’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}, 可选{None, int,}

如果 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(置信水平: float, 重采样次数: int) -> BootstrapSobolResult

提供指标置信区间的方法。更多详情请参见 scipy.stats.bootstrap

引导过程在首次和总阶敏感性指数上进行,这些结果在 BootstrapSobolResult 中作为属性 first_ordertotal_order 提供。

注释

Sobol’ 方法 [1], [2] 是一种基于方差的敏感性分析,它获取每个参数对感兴趣量(QoIs;即 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 应该会有所帮助。

进行良好分析所需的样本数量随着问题的维度增加而增加。例如,对于一个三维问题,考虑至少 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, and S. Tarantola. “全局敏感性分析。入门指南。” 2007.

[5]

Saltelli, A., P. Annoni, I. Azzini, F. Campolongo, M. Ratto, and S. Tarantola. “基于方差的模型输出敏感性分析。总敏感性指数的设计和估计器。” Computer Physics Communications, 181(2):259-270, DOI:10.1016/j.cpc.2009.09.018, 2010.

[6]

石神, T. 和 T. 本间. “计算机模型不确定性分析中的重要性量化技术.” 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])

置信区间可以通过自助法(bootstrapping)获得。

>>> 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()
../../_images/scipy-stats-sobol_indices-1_00_00.png

备注

默认情况下,scipy.stats.uniform 的支持范围是 [0, 1]。通过使用参数 locscale,可以得到在 [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()
../../_images/scipy-stats-sobol_indices-1_01_00.png

现在 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()
../../_images/scipy-stats-sobol_indices-1_02_00.png

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

这个案例在视觉上分析起来相当简单——尽管它只是一个定性分析。然而,当输入参数的数量增加时,这种分析变得不切实际,因为很难得出高阶项的结论。因此,使用Sobol’指数的好处就显现出来了。