前言

最近实习，敏度部长给了我一个灵敏度分析算法的分析法工作。灵敏度分析分为局部灵敏度分析和全局灵敏度分析；局部灵敏度分析包括：直接求导法、敏度有限差分法、分析法格林函数法。敏度全局灵敏度分析算法有筛选法、分析法蒙特卡洛方法、敏度基于方差的分析法方法。筛选法主要有IFFD、敏度MOAT、分析法COTTER法等，敏度一般用于分析包含大量输入变量的分析法系统模型、计算量相对较小。敏度蒙塔卡罗方法是分析法一种主要基于统计学理论的方法，包括散点图法、敏度相关系数法、回归分析法等，对于线性单调模型分析具有较强的适用性。基于方差的方法是近年来研究较多且应用最广的一类方法，主要有重要性估计法、傅里叶幅值灵敏度测试法（FAST）、以及拓展傅里叶幅值灵敏度检验方法（EFAST）等。SOBOL全局灵敏度分析法同时具备了蒙塔卡罗法和方差法的特点。

这里我查阅资料，有一些博主所写的只是在sobol上进行修改来实现morri，这里我试了反正效果不尽如人意，这里我最后还是上到SALib官网上取查对应的模块进行修正，得到我想要的预期结果。这里SALib的官网为：SALib网址，感谢博主忘荃小学和博主猪冰龙

MOAT筛选法（morris筛选法）

通过对可行域内随机取值，大大降低了对初始值的依赖程度，但在随机采样过程中容易出现误差，因而需要多次取平均结果来体现全局性。MOAT方法把灵敏度结果表示为EE（elementary effect），这里在下面的python代码中可以有着很清楚的体现。

#先引入我们需要的包from SALib.sample import saltellifrom SALib.analyze import sobolfrom SALib.analyze import morrisfrom SALib.test_functions import Ishigamiimport numpy as np import math from SALib.plotting.bar import plot as barplotimport matplotlib.pyplot as pltimport pandas as pd

定义模型的输入

我们这里验证的函数是： y = a x 2 + b x + c y=ax^2+bx+c y=ax2+bx+c这里对函数灵敏度有影响的参数为 a , b , c a,b,c a,b,c事先定义他们的范围。

#定义一个字典problem = {     'num_vars':3,    'names':['x1', 'x2', 'x3'],    'bounds':[[1, 3],[2, 4],[3, 5]]}

生成样本（samples）

我们要生成一样本(samples),由于我们要进行morris灵敏度分析我们需要用Saltelli samplers生成很多samples。

param_values = saltelli.sample(problem,1000)np.savetxt("param_values.txt", param_values)#将参数变化保存，其实是个参数范围内的sobol抽样print(param_values.shape)

(8000, 3)

param_values 是一个Numpy的矩阵，我们将打印出param_values的shape为（8000，3），Saltelli生成 N ∗ ( 2 D + 2 ) N*(2D+2) N∗(2D+2)个samples，N是我们在函数中所给定的，D是模型的三输入个数。

我们需要经过循环求出每个samples的模型输出值

#1-自定义-2Y = np.zeros([param_values.shape[0]])for i,X in enumerate(param_values):    tarr = np.arange(-5,6,1);    yerror = 0.0;    for t in tarr:        ylab = 2*t**2+3*t+4;        ytheory = X[0]*t**2+X[1]*t+X[2];        yerror = yerror+(ylab-ytheory)**2;    Y[i] = math.sqrt(yerror);

这里我们要计算敏感指数，在这个例子中我们采用的是 m o r r i s . a n a l y z e morris.analyze morris.analyze将会计算出平均灵敏、平均灵敏度的绝对值和灵敏度的标准差。

Si = morris.analyze(problem,param_values,Y,print_to_console=True)print("the mean elementary effect mu:")print(Si['mu'])#the mean elementary effectprint("the absolute of the mean elementary effect mu_star:")print(Si['mu_star'])#the absolute of the mean elementary effectprint("the standard deviation of the elementary effect sigma sigma:")print(Si['sigma'])#the standard deviation of the elementary effect

mu    mu_star      sigma  mu_star_confx1 -0.239573  31.922459  41.718545      1.009824x2 -0.035074  35.178599  39.714890      0.781839x3  0.051262  18.772476  27.984910      0.873359the mean elementary effect mu:[-0.23957343 -0.03507359  0.05126165]the absolute of the mean elementary effect mu_star:[31.922458653757158 35.17859916653295 18.772476282420456]the standard deviation of the elementary effect sigma:[41.71854476 39.71488989 27.98491008]

我们利用SALib提供的可视化将上述的灵敏度展示出来：

from SALib.plotting.bar import plot as barplotimport matplotlib.pyplot as plotSi_df = Si.to_df()print(Si_df)barplot(Si_df)

mu    mu_star      sigma  mu_star_confx1 -0.239573  31.922459  41.718545      1.310684x2 -0.035074  35.178599  39.714890      0.812093x3  0.051262  18.772476  27.984910      0.863634