sobol灵敏度分析matlab_Python中的模型敏感度分析(使用Salib)

敏感度分析的基础概念

文本主要参考了维基百科（对其中的关键部分进行了摘选了翻译）：

https://en.wikipedia.org/wiki/Sensitivity_analysis​en.wikipedia.org

敏感度分析是指对一个模型输出中的不确定性进行研究,并进一步判断不确定性的来源,也就是研究哪个输入参数的改变造成的输出变化的程度大小. 所以灵敏度分析是进行数学建模过程中一个必不可少的常规步骤.

选择敏感度分析方法的时候需要考虑的要素：

• 每运行一次模型的计算代价
• 输入参数之间的相关性
• 模型的响应是否非线性
• 输入因素之间的相互作用
• 已有数据的输入范围

常见的敏感度分析方法（跳过了傅里叶分析相关的方法，只有最后两个方法是全局分析）：

• One at a time(OAT) 方法
  • 每次变动一个输入并检查对于输出的影响。
  • 这种方法很简单，但由于它没有考虑输入变量的同时变化，因此并未充分探索输入空间。也无法检测输入变量之间是否存在交互。
• Screening方法
  • 窗口法是一种基于采样的方法。目的是要确定哪些输入变量对高维模型中的输出不确定性有重大影响，而不是准确地量化灵敏度。
  • 它具有相对较低的计算成本，并且可以在对其余集合应用更具信息性的分析之前，用于初步分析中以清除无影响的变量。
  • 最常用的筛选方法之一是基本效应方法（moris方法）
• 散点图法
• 基于偏导数的局部分析法
  • 检查输出对于各个输入的偏导数
  • 也无法充分探索输入空间
  • Adjoint modelling and Automated Differentiation 都属于这类方法
• 回归分析
  • 在敏感性分析的背景下，回归分析包括将线性回归拟合到模型响应，并使用标准化回归系数作为敏感性的直接度量。
  • 因此，当模型响应实际上是线性时，此方法最合适。例如，如果确定系数大，则可以确认回归模型有效。
  • 回归分析的优点是简单且计算成本低。
• 基于方差的方法（sabol方法）
  • 是全局方法。可以处理非线性响应，并且可以度量非加性系统中相互作用的影响。
  • 属于概率论方法。可将输入和输出不确定性量化为概率分布，并将输出方差分解为可归因于输入变量和变量组合的部分。因此，输出对输入变量的敏感度通过该输入引起的输出变化量来度量。
  • 常常会涉及蒙特卡洛方法
• 基于响应面的方差分析（VARS）
  • 这种通过一系列多个变量、确定性的“试验”，来模拟真实极限状态曲面的方法称为响应面法。

相关软件分析工具的外部链接

• SimLab，联合研究中心用于全球敏感性分析的免费软件
• 灵敏度分析Excel加载项是一个免费的（供私人和商业使用）Excel加载项，它允许基于样本的简单灵敏度分析运行
• MUCM项目 –用于计算需求模型的不确定性和敏感性分析的大量资源。
• GEM-SA –使用高斯过程执行灵敏度分析的程序。
• Python中的SALib灵敏度分析库（Numpy）。包含Sobol，Morris，分数阶乘和FAST方法。
• (另外在matlab 和 sas， spss 等商业软件中也都有敏感度分析工具)

Python中的敏感度分析

python中的敏感度分析工具是salib,可以直接通过pip安装, 而且支持的敏感度分析方法能够满足很多需求. 下面这个例子是我完善了补充了该库官方的入门示例后的小成果, 更加通俗易懂(使用sobol方法分析了模型 y = x[0]+ x[1]* cos(2* x[2] ) 的灵敏度):

从整体步骤上说, salib先提供采样输入, 拿到采样输入后我们讲其送入我们想要分析的自定义的模型(任何python实现的模型都可以), 然后再将模型的对应输出送回 salib就可以完成分析了.

from SALib.sample import saltelli
from SALib.analyze import sobol
import numpy as np
import math

# Define the model inputs
problem = {
    'num_vars': 3,
    'names': ['x1', 'x2', 'x3'],
    'bounds': [[-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359]]
}


def evaluate(X):  # 这里是我们要进行灵敏度分析的模型,接受一个数组,每个数组元素作为模型的一个输入,模型的输出是一个float,干函数返回的时候再讲所有输出并起来
    return np.array([math.sin(x[0]) + x[1] * math.cos(2 * x[2]) for x in X])


# Generate samples
param_values = saltelli.sample(problem, 1000)

# Run model (example)
Y = evaluate(param_values)
print(param_values.shape, Y.shape)
# Perform analysis (这里运行完成后会自动对结果进行展示)
Si = sobol.analyze(problem, Y, print_to_console=True)
print()

# Print the first-order sensitivity indices  一阶灵敏度
print('S1:', Si['S1'])

# Print the second-order sensitivity indices   二阶灵敏度
print("x1-x2:", Si['S2'][0, 1])
print("x1-x3:", Si['S2'][0, 2])
print("x2-x3:", Si['S2'][1, 2])

该程序的输出如下:

(8000, 3) (8000,)

Parameter S1 S1_conf ST ST_conf
x1 0.241379 0.040549 0.233006 0.020597
x2 -0.001635 0.075150 0.757444 0.083538
x3 -0.019499 0.082548 0.746122 0.064596

Parameter_1 Parameter_2 S2 S2_conf
x1 x2 -0.002905 0.057218
x1 x3 -0.007212 0.051962
x2 x3 0.768228 0.138674

S1: [ 0.24137922 -0.00163541 -0.01949871]
x1-x2: -0.0029050009015768627
x1-x3: -0.007212159874777576
x2-x3: 0.7682277244642012

输出结果的一般形式：Si是一个字典，含有"S1"， "S2"，"ST"，"S1_conf"，"S2_conf"，和 "ST_conf"项目。

其中"S1"， "S2"，"ST" 分别表示一阶灵敏度，二阶灵敏度和总灵敏度。（一个二阶灵敏度会涉及两个变量）, 而对应他们名字加上_conf 则表示对应的置信度（95%）。

灵敏度分析结果出现负数，说明灵敏度很小，并且因为样本太小所以出现了一些统计误差。

然后我们还能使用salib自带的灵敏度结果可视化功能:

from SALib.plotting.bar import plot as barplot
import matplotlib.pyplot as plot

Si_df = Si.to_df()
barplot(Si_df[0])
plot.show()

模型参数总灵敏度可视化

需要注意的是,原salib文档中的画图说明中有错误, 需要在 Si_df 后面加入一个索引再送入画图函数,原说明文档中缺少了一个索引,所以会出错.