sobol_indices#
- scipy.stats.sobol_indices(*, func, n, dists=None, method='saltelli_2010', rng=None)[source]#
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 是一个可调用对象，则必须指定，否则将被忽略。
- methodCallable 或 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)。这是一个高级功能，误用可能导致错误的分析。- rng
numpy.random.Generator，可选 伪随机数生成器状态。 当 rng 为 None 时，将使用来自操作系统的熵创建一个新的
numpy.random.Generator。 除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_order形状为 (s, d) 的 ndarray
一阶 Sobol’ 指标。
- total_order形状为 (s, d) 的 ndarray
总阶 Sobol’ 指标。
和方法
bootstrap(confidence_level: float, n_resamples: int) -> BootstrapSobolResult
一种在指标上提供置信区间的方法。有关更多详细信息，请参阅
scipy.stats.bootstrap。引导程序在第一阶和总阶指标上完成，它们作为属性
first_order和total_order在 BootstrapSobolResult 中可用。
注释
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.. “非线性数学模型的灵敏度分析”。 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, 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]Ishigami, T. and 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’ 指标的好处就体现出来了。