python如何做敏感度分析_1stopt、matlab和python用morris、sobol方法实现参数敏感性分析...

首先看抛物线函数：

现在我取a=2,b=3,c=4 ,得到如下函数：

x或t都是指自变量，就不改了，一个意思。

问题是，我想知道对于此数据和模型，参数a,b,c的敏感性，也就是y的改变量与a、b、c的改变量的比值关系。

首先用1stopt分析：

1 NewCodeBlock"SA";2 Parameter a=[1,3],b=[2,4],c=[3,5];//要优化的参数及其范围3 Variable t,y;//变量4 Function y=a*t^2+b*t+c;5 Data;6 //t y 变量顺序要和Variable后变量对应

7 -5 39

8 -4 24

9 -3 13

10 -2 6

11 -1 3

12 0 4

13 1 9

14 2 18

15 3 31

16 4 48

17 5 69

以下分别为各种参数敏感性方法(包括morris局部和全局，偏微分方法局部和全局)：

以上就是4中方法的结果。用的目标函数是SSR,点的范围我用的是上下浮动50%，正好布满整个给定空间：

可以看morris局部敏感度分析具体数据：

a：

b：

c：

然后将灵敏度指数平均得到morris局部方法(本质是OAT方法，一次改变一个)分析结果：

全局方法(本质是通过拉丁抽样实现同时考虑多个参数的影响)就不是这个思路了，看一下超拉丁抽样：

总结一下1stopt中4种方法参数灵敏度结果：

morris方法局部：

名称灵敏度指数灵敏度指数(%)

a3.916E1889.1113892365457

b4.125E179.38673341677096

c6.6E161.50187734668335

morris方法全局：

名称灵敏度指数灵敏度指数(%)

a2.00290323137447E1881.4890482414753

b2.10289523274038E178.55572692595605

c2.44687506069835E179.95522483256862

偏微分方法局部：

名称灵敏度指数灵敏度指数(%)

a3.99432000000567E1889.1113892365577

b4.20750000000213E179.38673341676365

c6.73199999998753E161.50187734667864

偏微分方法全局：

名称灵敏度指数灵敏度指数(%)

a3.98876054886661E1882.084202339751

b4.20542269921853E178.65428655186941

c4.50049450186404E179.26151110837963

下面用matlab写sobol方法进行分析：

1 %sobol 参数敏感性分析2 %算法参考：3 % csdn : https://blog.csdn.net/xiaosebi1111/article/details/46517409

4 % wiki： https://5 %运行环境 matlab2016b6 %作者 blzhu@buaa.edu.cn 2020年6月7日7 %%初始化8 clc;9 clear all;10 close all;11 %%设定：给定参数个数和各个参数的范围12

13 % % 1-自定义子函数114 % D=3;%维度3,几个参数15 % nPop=4:500:5000;%采样点个数,也就是参数水平数 ,取大了好，比如4000，但慢16 % VarMin=[-pi -pi -pi ];%各个参数下限 SALib ：S1: [ 0.30797549 0.44776661 -0.00425452] ；ST: [0.56013728 0.4387225 0.24284474]17 % VarMax=[pi pi pi];%各个参数上限18 % myFunction=@(x) Ishigami(x);%目标函数，也可以是个黑盒子19

20 % % 1-自定义子函数221 D=3;%维度3,几个参数22 nPop=4:50:1000;%采样点个数,也就是参数水平数 ,取大了好，比如4000，但慢23 VarMin=[1 2 3 ];%各个参数下限24 VarMax=[3 4 5];%各个参数上限25 myFunction=@(x) myx2(x);26

27 % % 1-自定义子函数328 % D=4;%维度3,几个参数29 % nPop=4:50:1000;%采样点个数,也就是参数水平数 ,取大了好，比如4000，但慢30 % VarMin=[1 2 3 4];%各个参数下限31 % VarMax=[3 4 5 5];%各个参数上限32 % myFunction=@(x) myx3(x);33

35 %%开始计算36 numnPop=numel(nPop);37 SAll=zeros(numnPop,D+1);%分别是各参数的敏感度，最后一列是各参数敏感度之和38 STAll=zeros(numnPop,D+1);39 for i=1:numnPop40 [S,ST]=sobol(D,nPop(i),VarMin,VarMax,myFunction);41 SAll(i,1:D)=S';

42 SAll(i,D+1)=sum(SAll(i,1:D));43 STAll(i,1:D)=ST';

44 STAll(i,D+1)=sum(STAll(i,1:D));45 end46 %%绘图47 color=[1 0 0;0 1 0;0 0 1;0.5 0.1 0;0 0.3 0.4;0.6 0.7 0.2;0.5 0.8 0.9;0 0.2 0.1;0.1 0.5 0;0.1 0 0.5;0.5 0 0.1];%12种颜色一般颜色不一样48 marker=['o','+','*','.','x','s','d','^','v','>','

52 figure(1)53 for i=1:D+1

54 plot(nPop,SAll(:,i),'Marker',marker(i),'LineStyle',char(linestyle(useL)),'Color',color(i,:),'LineWidth',1);55 hold on56 end57 title('Sobol-S');58 whichPara=sprintfc('%g',repmat(1:D+1,1,2));%把数字数组转化成字符串数组59 legend(whichPara,'Location','bestoutside');%加图例60

62 figure(2)63 for i=1:D+1

64 plot(nPop,STAll(:,i),'Marker',marker(i),'LineStyle',char(linestyle(useL)),'Color',color(i,:),'LineWidth',1);65 hold on66 end67 title('Sobol-ST');68 whichPara=sprintfc('%g',repmat(1:D+1,1,2));%把数字数组转化成字符串数组69 legend(whichPara,'Location','bestoutside');%加图例70

71 disp('一阶影响指数(左方向收敛于1)Sobol-S：');72 disp(S);73 disp('总效应指数(大于等于1，且仅当myfun是纯相加时取等号)Sobol-ST：');74 disp(ST);75 disp(datetime);76 disp('parameter sensitive analyse success use sobol method!');77 %%火车声音提示已经算完了78 load train79 sound(y,Fs)80

83 %% -------------------------子函数 matlab2016之前不支持子函数写在同一个m文档中----------------------------

84 % 1-自定义子函数1(3个参数)Ishigami https://www.sfu.ca/~ssurjano/ishigami.html

85 function y=Ishigami(x)86 y=sin(x(1))+7*(sin(x(2)))^2+0.1*x(3)^4*sin(x(1));%SALib用的这个87 % y=sin(x(1))+7*(sin(x(2)))^2+0.05*x(3)^4*sin(x(1));88 end89

90 % 1-自定义子函数2 (3个参数)91 function y=myx2(x)92 t=-5:1:5;%与此处有t范围和步距有关系93 % t=-5:0.1:5;%与此处有t范围和步距有关系94 ylab=2*t.^2+3*t+4;95 ytheory=x(1)*t.^2+x(2)*t+x(3);96 y=sum((ylab-ytheory).^2);%残差平方和SSR作为目标函数97 % y=sum((ylab-ytheory).^2)/numel(t);%各参数灵敏度与上式相同98 end99

100

101 % 1-自定义子函数3(4个参数)102 function y=myx3(x)103 t=-5:1:5;104 ylab=2*t.^3+3*t.^2+4*t+5;105 ytheory=x(1)*t.^3+x(2)*t.^2+x(3)*t+x(4);106 y=sum((ylab-ytheory).^2);107 end108

109

110

111 %% 2-求sobol敏感度112 function [S,ST]=sobol(D,nPop,VarMin,VarMax,myFunction)113 M=D*2;%

114 %%产生所需的各水平参数115 VarMin=[VarMin,VarMin];116 VarMax=[VarMax,VarMax];117 p= sobolset(M);% https://www.cnblogs.com/zhubinglong/p/12260292.html

118 % R=p(1:nPop,:);%我只用前nPop个119 R=[];120 for i=1:nPop121 r=p(i,:);122 r=VarMin+r.*(VarMax-VarMin);123 R=[R; r];124 end125 % plot(R(:,1),'b*')126 %拆分为A B127 A=R(:,1:D);%每行代表一组参数，其中每列代表每组参数的一个参数；行数就代表共有几组参数128 B=R(:,D+1:end);129 %根据A B 产生矩阵AB130 AB=zeros(nPop,D,D);131 for i=1:D132 tempA=A;133 tempA(:,i)=B(:,i);134 AB(1:nPop,1:D,i)=tempA;135 end136 %%求各参数解137 YA=zeros(nPop,1);%解138 YB=zeros(nPop,1);139 YAB=zeros(nPop,D);%分别代表YAB1,YAB2,YAB3，YAB(:,D)就代表YABD140 for i=1:nPop141 YA(i)=myFunction(A(i,:));142 YB(i)=myFunction(B(i,:));143 for j=1:D144 YAB(i,j)=myFunction(AB(i,:,j));145 end146 end147 %%根据一阶影响指数公式：148 VarX=zeros(D,1);%S的分子149 S=zeros(D,1);150

151 % 0: 估算基于给定样本的方差(EXCEL var.p) ; 1:计算基于给定的样本总体的方差(EXCEL var.p())152 % var([2.091363878 1.110366059 3.507651769 1.310950363 2.091363878 3.507651769 1.110366059 1.7066512],1);153 VarY=var([YA;YB],1,'omitnan');% S的分母。计算基于给定的样本总体的方差(EXCEL var.p())154 for i=1:D155 for j=1:nPop156 VarX(i)=VarX(i)+YB(j)*(YAB(j,i)-YA(j));157 end158 VarX(i)=1/nPop*VarX(i);%蒙特卡罗估计量159 S(i)=VarX(i)/VarY;160 end161

162 %%总效应指数163 EX=zeros(D,1);164 ST=zeros(D,1);165 for i=1:D166 for j=1:nPop167 EX(i)=EX(i)+(YA(j)-YAB(j,i))^2;168 end169 EX(i)=1/(2*nPop)* EX(i);%蒙特卡罗估计量170 ST(i)=EX(i)/VarY;171 end172 end

我分别取了不同个数的样点 4:50:1000 ，结果如下，可见1000个样点基本稳定了。

各参数的灵敏度：

一阶影响指数(左方向收敛于1)Sobol-S：

0.9728

0.0030

0.0001

总效应指数(大于等于1，且仅当myfun是纯相加时取等号)Sobol-ST：

0.9860

0.0031

0.0155

当然，也可以在matlab的fileexchange上下载各种工具箱，但这个根据csdn和wiki上写的算法相对简单些，便于魔改。

用python的SALib包分析

Method of Morris, including groups and optimal trajectories ([Morris 1991], [Campolongo et al. 2007])

Fourier Amplitude Sensitivity Test (FAST) ([Cukier et al. 1973], [Saltelli et al. 1999])

Delta Moment-Independent Measure ([Borgonovo 2007], [Plischke et al. 2013])

Derivative-based Global Sensitivity Measure (DGSM) ([Sobol and Kucherenko 2009])

Fractional Factorial Sensitivity Analysis ([Saltelli et al. 2008])

下面以sobol方法举例：

1 #https://salib.readthedocs.io/en/latest/basics.html#run-model

2 #-*- coding: utf-8 -*-

3 from SALib.sample importsaltelli4 from SALib.analyze importsobol5 from SALib.test_functions importIshigami6 importnumpy as np7 importmath8 from SALib.plotting.bar importplot as barplot9 importmatplotlib.pyplot as plot10

11 #problem = {

12 #'num_vars': 3,

13 #'names': ['x1', 'x2', 'x3'],

14 #'bounds': [[-3.14159265359, 3.14159265359],

15 #[-3.14159265359, 3.14159265359],

16 #[-3.14159265359, 3.14159265359]]

17 #}

19 problem ={20 'num_vars': 3,21 'names': ['x1', 'x2', 'x3'],22 'bounds': [[1, 3],23 [2, 4],24 [3, 5]]25 }26

28 param_values = saltelli.sample(problem, 1000)#不管用哪个方法计算y,这个要有

29 np.savetxt("param_values.txt", param_values)#将参数变化保存，其实是各参数范围内的sobol抽样

31 ## 计算Y

32 ##1-自定义-1

33 #Y = np.zeros([param_values.shape[0]])

34 #A = 7

35 #B = 0.1

36 #for i, X in enumerate(param_values):

37 #Y[i] = math.sin(X[0]) + A * math.pow(math.sin(X[1]), 2) + \

38 #B * math.pow(X[2], 4) * math.sin(X[0])

40 #1-自定义-2

41 Y =np.zeros([param_values.shape[0]])42 for i, X inenumerate(param_values):43 tarr=np.arange(-5,6,1);44 yerror=0.0;45 for t intarr:46 ylab=2*t**2+3*t+4;47 ytheory=X[0]*t**2+X[1]*t+X[2];48 yerror=yerror+(ylab-ytheory)**2;49

50 Y[i] =math.sqrt(yerror);51

52 #2-load计算好的txt

53 #Y = np.loadtxt("outputs.txt", float)

55 #3-SALib自带测试函数

56 #Y = Ishigami.evaluate(param_values)

57 ## np.savetxt("outputs.txt", Y)#将因变量变化结果保存

59 Si = sobol.analyze(problem, Y ,print_to_console=True)60 print()#自动输出S1(单个参数对因变量的影响)、ST(考虑各个变量相互影响)和S2(两两参数之间影响)，需有,print_to_console=True

62 print("all parameters first-order sensitivity indices S1:")63 print(Si['S1'])#一阶影响指数

64 print("all parameters second-order sensitivity indices S2:")65 print(Si['S2'])#二阶影响指数

66 print("all parameters total sensitivity indices ST:")67 print(Si['ST'])#总效应指数

69 #绘图 https://zhuanlan.zhihu.com/p/137953265

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

输出结果：

Parameter S1 S1_conf ST ST_conf

x1 0.969397 0.069624 0.982232 0.058139

x2 0.007222 0.009055 0.009680 0.001325

x3 0.000848 0.008468 0.011699 0.001033

Parameter_1 Parameter_2 S2 S2_conf

x1 x2 0.000330 0.070070

x1 x3 0.010014 0.069389

x2 x3 -0.000129 0.013788

all parameters first-order sensitivity indices S1:

[9.69397123e-01 7.22243327e-03 8.47690887e-04]

all parameters second-order sensitivity indices S2:

[[ nan 0.00032978 0.01001386]

[ nan nan -0.00012935]

[ nan nan nan]]

all parameters total sensitivity indices ST:

[0.9822323 0.00968001 0.01169928]

也可以用morris方法：

只需要导入

from SALib.analyze import morris

然后用

Si = morris.analyze(problem, Y ,print_to_console=True) 代替

Si = sobol.analyze(problem, Y ,print_to_console=True)

1 #https://salib.readthedocs.io/en/latest/basics.html#run-model

2 #-*- coding: utf-8 -*-

3 from SALib.sample importsaltelli4 from SALib.analyze importsobol5 from SALib.analyze importmorris6 from SALib.test_functions importIshigami7 importnumpy as np8 importmath9 from SALib.plotting.bar importplot as barplot10 importmatplotlib.pyplot as plot11

12 #problem = {

13 #'num_vars': 3,

14 #'names': ['x1', 'x2', 'x3'],

15 #'bounds': [[-3.14159265359, 3.14159265359],

16 #[-3.14159265359, 3.14159265359],

17 #[-3.14159265359, 3.14159265359]]

18 #}

20 problem ={21 'num_vars': 3,22 'names': ['x1', 'x2', 'x3'],23 'bounds': [[1, 3],24 [2, 4],25 [3, 5]]26 }27

29 param_values = saltelli.sample(problem, 1000)#不管用哪个方法计算y,这个要有

30 np.savetxt("param_values.txt", param_values)#将参数变化保存，其实是各参数范围内的sobol抽样

32 ## 计算Y

33 ##1-自定义-1

34 #Y = np.zeros([param_values.shape[0]])

35 #A = 7

36 #B = 0.1

37 #for i, X in enumerate(param_values):

38 #Y[i] = math.sin(X[0]) + A * math.pow(math.sin(X[1]), 2) + \

39 #B * math.pow(X[2], 4) * math.sin(X[0])

41 #1-自定义-2

42 Y =np.zeros([param_values.shape[0]])43 for i, X inenumerate(param_values):44 tarr=np.arange(-5,6,1);45 yerror=0.0;46 for t intarr:47 ylab=2*t**2+3*t+4;48 ytheory=X[0]*t**2+X[1]*t+X[2];49 yerror=yerror+(ylab-ytheory)**2;50

51 Y[i] =math.sqrt(yerror);52

53 #2-load计算好的txt

54 #Y = np.loadtxt("outputs.txt", float)

56 #3-SALib自带测试函数

57 #Y = Ishigami.evaluate(param_values)

58 ## np.savetxt("outputs.txt", Y)#将因变量变化结果保存

60 #Si = sobol.analyze(problem, Y ,print_to_console=True)

61 Si = morris.analyze(problem, Y ,print_to_console=True)62 print()#自动输出S1(单个参数对因变量的影响)、ST(考虑各个变量相互影响)和S2(两两参数之间影响)，需有,print_to_console=True

64 print("all parameters first-order sensitivity indices S1:")65 print(Si['S1'])#一阶影响指数

66 print("all parameters second-order sensitivity indices S2:")67 print(Si['S2'])#二阶影响指数

68 print("all parameters total sensitivity indices ST:")69 print(Si['ST'])#总效应指数

71 #绘图 https://zhuanlan.zhihu.com/p/137953265

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

python3.8.1+SALib1.3 use morris

Parameter S1 S1_conf ST ST_conf

x1 0.969397 0.072129 0.982232 0.062795

x2 0.007222 0.008033 0.009680 0.001303

x3 0.000848 0.009519 0.011699 0.001045

Parameter_1 Parameter_2 S2 S2_conf

x1 x2 0.000330 0.087924

x1 x3 0.010014 0.087743

x2 x3 -0.000129 0.013700

all parameters first-order sensitivity indices S1:

[9.69397123e-01 7.22243327e-03 8.47690887e-04]

all parameters second-order sensitivity indices S2:

[[ nan 0.00032978 0.01001386]

[ nan nan -0.00012935]

[ nan nan nan]]

all parameters total sensitivity indices ST:

[0.9822323 0.00968001 0.01169928]

小技巧：

SALib如果不便于将目标函数写为函数的形式，可以将代码：

np.savetxt("param_values.txt", param_values)# 将参数变化保存，其实是各参数范围内的sobol抽样

生成的抽样带入自己的系统，然后根据自己需要生成对应抽样的目标函数，将目标函数放入同目录下的 outputs.txt 文档中，一行一个结果，然后用这个语句代替上面求Y[i]:

Y = np.loadtxt("outputs.txt", float)

就是说提供了参数变化以及目标函数变化，用SALib就可以求参数灵敏度了。

总结：

matlab我用的sobol生成的抽样和别人的不一样，不知道为什么，这个是造成与python计算不一样的一个原因吧。但差别不大。

python如何做敏感度分析_1stopt、matlab和python用morris、sobol方法实现参数敏感性分析...相关推荐

python如何做敏感度分析_Python中的模型敏感度分析(使用Salib)
敏感度分析的基础概念文本主要参考了维基百科(对其中的关键部分进行了摘选了翻译):https://en.wikipedia.org/wiki/Sensitivity_analysisen.wikip ...
python怎么做情感分析_如何用python进行情感分析
我们在计划中遵循以下三个主要步骤:授权twitter API客户端. 向Twitter API发出GET请求以获取特定查询的推文. 解析推文.将每条推文分类为正面,负面或中立. 首先,我们创建一个Tw ...
python能做财务分析吗_您可以使用Python进行财务规划和分析吗？
python能做财务分析吗问题 (The Problem) If you work in the Financial Planning and Analysis area, chances are ...
python可以做计量分析吗_技术分享 - python数据分析（2）——数据特征分析（上）...
1 分布分析分布分析能揭示数据的分布特征和分布类型.对于定量数据,欲了解其分布形式是对称的还是非对称的,发现某些特大或特小的可疑值,可通过绘制频率分布表.绘制频率分布直方图.绘制茎叶图进行直观地分析 ...
学python可以做什么职业-不知道学了Python能干嘛？Python职业发展：7大职位供你选择！...
原标题:不知道学了Python能干嘛?Python职业发展:7大职位供你选择! 为什么那么多小伙伴都在学Python呢?Python到底有啥魔力?学了Python都能干啥?这篇文章,肉丝儿来和大家一起 ...
学python可以做什么职业好-业余学Python能做什么?对职业发展有什么帮助？
业余学Python能做什么?一般来说,Python有Web开发.数据科学和脚本三大应用.无论对于零基础小白,还是已经工作想要提升自己的在职人员,学好这些内容都会对职业发展有着重要作用.下面小编将详细为 ...
学python能做什么-非计算机专业的人学python能做什么？
Python职场必备进阶技能 python和SQL数据分析可以完全取代excel,这个在一些比较大型的咨询公司,财务公司.商业分析领域已经是一个趋势,岗位越高,python是必备技能.大数据时代,数据 ...
python能做什么工作-谁适合学Python？学了Python可以做什么工作？
Tips: 目前在很多行业中都在越来越多的应用Python,这也是很多行业学习Python的原因,Python主要的应用领域有哪些呢?今天我们就来详细看一下. 谁适合学Python? 我们首先来看一看 ...
python可以做动漫吗_用Python做一个以图搜番的应用程序，再也不用愁动漫图片的出处了！...
前言喜欢看动漫的朋友们大概都能体会到一个难受的事情,就是在论坛或者群聊里面看到一张动漫截图,很想知道它的出处,但百度搜了一圈却也没有一个可靠结果,就很郁闷.今天就来带大家用Python做一个简单的& ...

python如何做敏感度分析_1stopt、matlab和python用morris、sobol方法实现参数敏感性分析...

python如何做敏感度分析_1stopt、matlab和python用morris、sobol方法实现参数敏感性分析...相关推荐

最新文章

热门文章