Python SciPy stats.sobol_indices用法及代码示例

本文简要介绍 python 语言中 scipy.stats.sobol_indices 的用法。

用法:

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

索博尔的全局敏感性指数。

参数 ：：

func： 可调用或 dict(str, 数组)

如果 func 是可调用的函数，则用于计算 Sobol’ 索引。其签名必须是：

func(x: ArrayLike) -> ArrayLike

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

• d是输入维度函数(输入变量的数量)，

• s是输出维度函数(输出变量的数量)，以及

• n是样本数(参见n以下)。

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

如果函数是一个字典，包含来自三个不同数组的函数计算。 key 必须是：f_A,f_B和f_AB.f_A和f_B应该有一个形状(s, n)和f_AB应该有一个形状(d, s, n)。这是一项高级函数，误用可能会导致错误的分析。

n： int

用于生成矩阵的样本数A和B。必须是 2 的幂。函数被评估将是n*(d+2).

dists： 列表(分布)，可选

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

发行版必须是具有 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]

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

random_state： {无，整数， numpy.random.Generator }，可选

如果random_state是 int 或 None，一个新的numpy.random.Generator是使用创建的np.random.default_rng(random_state).如果random_state已经是一个Generator实例，然后使用提供的实例。

返回 ：：

res： SobolResult

具有属性的对象：

first_order 形状为 (s, d) 的 ndarray

一阶 Sobol 指数。

total_order 形状为 (s, d) 的 ndarray

总订单 Sobol 指数。

及方法：

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

A method providing confidence intervals on the indices. See scipy.stats.bootstrap for more details.

The bootstrapping is done on both first and total order indices, and they are available in BootstrapSobolResult as attributes first_order and total_order.

注意：

索博尔法[1],[2]是基于方差的敏感性分析，它获得每个参数对感兴趣量方差(QoIs；即输出)的贡献函数)。各自的贡献可用于对参数进行排序，并通过计算模型的有效(或平均)维度来衡量模型的复杂性。

注意

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

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

\[\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),\]

索博尔指数表示为

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

\(S_{i}\) 对应于 first-order 项，它告知 i-th 参数的贡献，而 \(S_{ij}\) 对应于二阶项，它告知 i-th 和 j-th 参数之间相互作用的贡献。这些方程可以推广到计算高阶项；然而，它们的计算成本很高，而且它们的解释也很复杂。这就是为什么只提供一阶索引的原因。

总订单指数表示参数对 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。

警告

负索博尔值是由于数值错误造成的。增加点数 n 应该会有所帮助。

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

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

参考：

[1]

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

[2]

索博尔，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。 “模型输出的基于方差的敏感性分析。总灵敏度index的设计和估计。”计算机物理通信，181(2):259-270，DOI:10.1016/j.cpc.2009.09.018，2010。

[6]

石上，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 指数的好处。