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)
的输出,其中:d
是 func 的输入维度(输入变量的数量),s
是 func 的输出维度(输出变量的数量),并且n
是样本的数量(见下面的 n)。
函数求值结果必须是有限的。
如果 func 是一个字典,包含来自三个不同数组的函数评估。键必须是:
f_A
、f_B
和f_AB
。f_A
和f_B
应该具有形状(s, n)
,而f_AB
应该具有形状(d, s, n)
。这是一个高级功能,误用可能导致错误的分析。- n整数
用于生成矩阵
A
和B
的样本数量。必须是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_order
和total_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()
备注
默认情况下,
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’指数的好处就显现出来了。