首先看拋物線函數(shù)：

現(xiàn)在我取a=2,b=3,c=4 ,得到如下函數(shù)：

x或t都是指自變量，就不改了，一個(gè)意思。

問題是，我想知道對于此數(shù)據(jù)和模型，參數(shù)a,b,c的敏感性，也就是y的改變量與a、b、c的改變量的比值關(guān)系。

首先用1stopt分析：

1 NewCodeBlock"SA";2 Parameter a=[1,3],b=[2,4],c=[3,5];//要優(yōu)化的參數(shù)及其范圍3 Variable t,y;//變量4 Function y=a*t^2+b*t+c;5 Data;6 //t y 變量順序要和Variable后變量對應(yīng)

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

以下分別為各種參數(shù)敏感性方法(包括morris局部和全局，偏微分方法局部和全局)：

以上就是4中方法的結(jié)果。用的目標(biāo)函數(shù)是SSR,點(diǎn)的范圍我用的是上下浮動(dòng)50%，正好布滿整個(gè)給定空間：

可以看morris局部敏感度分析具體數(shù)據(jù)：

a：

b：

c：

然后將靈敏度指數(shù)平均得到morris局部方法(本質(zhì)是OAT方法，一次改變一個(gè))分析結(jié)果：

全局方法(本質(zhì)是通過拉丁抽樣實(shí)現(xiàn)同時(shí)考慮多個(gè)參數(shù)的影響)就不是這個(gè)思路了，看一下超拉丁抽樣：

總結(jié)一下1stopt中4種方法參數(shù)靈敏度結(jié)果：

morris方法局部：

名稱 靈敏度指數(shù) 靈敏度指數(shù)(%)

a3.916E1889.1113892365457

b4.125E179.38673341677096

c6.6E161.50187734668335

morris方法全局：

名稱 靈敏度指數(shù) 靈敏度指數(shù)(%)

a2.00290323137447E1881.4890482414753

b2.10289523274038E178.55572692595605

c2.44687506069835E179.95522483256862

偏微分方法局部：

名稱 靈敏度指數(shù) 靈敏度指數(shù)(%)

a3.99432000000567E1889.1113892365577

b4.20750000000213E179.38673341676365

c6.73199999998753E161.50187734667864

偏微分方法全局：

名稱 靈敏度指數(shù) 靈敏度指數(shù)(%)

a3.98876054886661E1882.084202339751

b4.20542269921853E178.65428655186941

c4.50049450186404E179.26151110837963

下面用matlab寫sobol方法進(jìn)行分析：

1 %sobol 參數(shù)敏感性分析2 %算法參考：3 % csdn : https://blog.csdn.net/xiaosebi1111/article/details/46517409

4 % wiki： https://5 %運(yùn)行環(huán)境 matlab2016b6 %作者 blzhu@buaa.edu.cn 2020年6月7日7 %%初始化8 clc;9 clear all;10 close all;11 %%設(shè)定：給定參數(shù)個(gè)數(shù)和各個(gè)參數(shù)的范圍12

13 % % 1-自定義子函數(shù)114 % D=3;%維度3,幾個(gè)參數(shù)15 % nPop=4:500:5000;%采樣點(diǎn)個(gè)數(shù),也就是參數(shù)水平數(shù) ,取大了好，比如4000，但慢16 % VarMin=[-pi -pi -pi ];%各個(gè)參數(shù)下限 SALib ：S1: [ 0.30797549 0.44776661 -0.00425452] ；ST: [0.56013728 0.4387225 0.24284474]17 % VarMax=[pi pi pi];%各個(gè)參數(shù)上限18 % myFunction=@(x) Ishigami(x);%目標(biāo)函數(shù)，也可以是個(gè)黑盒子19

20 % % 1-自定義子函數(shù)221 D=3;%維度3,幾個(gè)參數(shù)22 nPop=4:50:1000;%采樣點(diǎn)個(gè)數(shù),也就是參數(shù)水平數(shù) ,取大了好，比如4000，但慢23 VarMin=[1 2 3 ];%各個(gè)參數(shù)下限24 VarMax=[3 4 5];%各個(gè)參數(shù)上限25 myFunction=@(x) myx2(x);26

27 % % 1-自定義子函數(shù)328 % D=4;%維度3,幾個(gè)參數(shù)29 % nPop=4:50:1000;%采樣點(diǎn)個(gè)數(shù),也就是參數(shù)水平數(shù) ,取大了好，比如4000，但慢30 % VarMin=[1 2 3 4];%各個(gè)參數(shù)下限31 % VarMax=[3 4 5 5];%各個(gè)參數(shù)上限32 % myFunction=@(x) myx3(x);33

34

35 %%開始計(jì)算36 numnPop=numel(nPop);37 SAll=zeros(numnPop,D+1);%分別是各參數(shù)的敏感度，最后一列是各參數(shù)敏感度之和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));%把數(shù)字?jǐn)?shù)組轉(zhuǎn)化成字符串?dāng)?shù)組59 legend(whichPara,'Location','bestoutside');%加圖例60

61

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));%把數(shù)字?jǐn)?shù)組轉(zhuǎn)化成字符串?dāng)?shù)組69 legend(whichPara,'Location','bestoutside');%加圖例70

71 disp('一階影響指數(shù)(左方向收斂于1)Sobol-S：');72 disp(S);73 disp('總效應(yīng)指數(shù)(大于等于1，且僅當(dāng)myfun是純相加時(shí)取等號)Sobol-ST：');74 disp(ST);75 disp(datetime);76 disp('parameter sensitive analyse success use sobol method!');77 %%火車聲音提示已經(jīng)算完了78 load train79 sound(y,Fs)80

81

82

83 %% -------------------------子函數(shù) matlab2016之前不支持子函數(shù)寫在同一個(gè)m文檔中----------------------------

84 % 1-自定義子函數(shù)1(3個(gè)參數(shù))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用的這個(gè)87 % y=sin(x(1))+7*(sin(x(2)))^2+0.05*x(3)^4*sin(x(1));88 end89

90 % 1-自定義子函數(shù)2 (3個(gè)參數(shù))91 function y=myx2(x)92 t=-5:1:5;%與此處有t范圍和步距有關(guān)系93 % t=-5:0.1:5;%與此處有t范圍和步距有關(guān)系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作為目標(biāo)函數(shù)97 % y=sum((ylab-ytheory).^2)/numel(t);%各參數(shù)靈敏度與上式相同98 end99

100

101 % 1-自定義子函數(shù)3(4個(gè)參數(shù))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 %%產(chǎn)生所需的各水平參數(shù)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個(gè)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);%每行代表一組參數(shù)，其中每列代表每組參數(shù)的一個(gè)參數(shù)；行數(shù)就代表共有幾組參數(shù)128 B=R(:,D+1:end);129 %根據(jù)A B 產(chǎn)生矩陣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 %%求各參數(shù)解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 %%根據(jù)一階影響指數(shù)公式：148 VarX=zeros(D,1);%S的分子149 S=zeros(D,1);150

151 % 0: 估算基于給定樣本的方差(EXCEL var.p) ; 1:計(jì)算基于給定的樣本總體的方差(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的分母。 計(jì)算基于給定的樣本總體的方差(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);%蒙特卡羅估計(jì)量159 S(i)=VarX(i)/VarY;160 end161

162 %%總效應(yīng)指數(shù)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);%蒙特卡羅估計(jì)量170 ST(i)=EX(i)/VarY;171 end172 end

我分別取了不同個(gè)數(shù)的樣點(diǎn)?4:50:1000 ，結(jié)果如下，可見1000個(gè)樣點(diǎn)基本穩(wěn)定了。

各參數(shù)的靈敏度：

一階影響指數(shù)(左方向收斂于1)Sobol-S：

0.9728

0.0030

0.0001

總效應(yīng)指數(shù)(大于等于1，且僅當(dāng)myfun是純相加時(shí)取等號)Sobol-ST：

0.9860

0.0031

0.0155

當(dāng)然，也可以在matlab的fileexchange上下載各種工具箱，但這個(gè)根據(jù)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 #}

18

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

27

28 param_values = saltelli.sample(problem, 1000)#不管用哪個(gè)方法計(jì)算y,這個(gè)要有

29 np.savetxt("param_values.txt", param_values)#將參數(shù)變化保存，其實(shí)是各參數(shù)范圍內(nèi)的sobol抽樣

30

31 ## 計(jì)算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])

39

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計(jì)算好的txt

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

54

55 #3-SALib自帶測試函數(shù)

56 #Y = Ishigami.evaluate(param_values)

57 ## np.savetxt("outputs.txt", Y)#將因變量變化結(jié)果保存

58

59 Si = sobol.analyze(problem, Y ,print_to_console=True)60 print()#自動(dòng)輸出S1(單個(gè)參數(shù)對因變量的影響)、ST(考慮各個(gè)變量相互影響)和S2(兩兩參數(shù)之間影響)，需有,print_to_console=True

61

62 print("all parameters first-order sensitivity indices S1:")63 print(Si['S1'])#一階影響指數(shù)

64 print("all parameters second-order sensitivity indices S2:")65 print(Si['S2'])#二階影響指數(shù)

66 print("all parameters total sensitivity indices ST:")67 print(Si['ST'])#總效應(yīng)指數(shù)

68

69 #繪圖 https://zhuanlan.zhihu.com/p/137953265

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

輸出結(jié)果：

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方法：

只需要導(dǎo)入

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 #}

19

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

28

29 param_values = saltelli.sample(problem, 1000)#不管用哪個(gè)方法計(jì)算y,這個(gè)要有

30 np.savetxt("param_values.txt", param_values)#將參數(shù)變化保存，其實(shí)是各參數(shù)范圍內(nèi)的sobol抽樣

31

32 ## 計(jì)算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])

40

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計(jì)算好的txt

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

55

56 #3-SALib自帶測試函數(shù)

57 #Y = Ishigami.evaluate(param_values)

58 ## np.savetxt("outputs.txt", Y)#將因變量變化結(jié)果保存

59

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

61 Si = morris.analyze(problem, Y ,print_to_console=True)62 print()#自動(dòng)輸出S1(單個(gè)參數(shù)對因變量的影響)、ST(考慮各個(gè)變量相互影響)和S2(兩兩參數(shù)之間影響)，需有,print_to_console=True

63

64 print("all parameters first-order sensitivity indices S1:")65 print(Si['S1'])#一階影響指數(shù)

66 print("all parameters second-order sensitivity indices S2:")67 print(Si['S2'])#二階影響指數(shù)

68 print("all parameters total sensitivity indices ST:")69 print(Si['ST'])#總效應(yīng)指數(shù)

70

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如果不便于將目標(biāo)函數(shù)寫為函數(shù)的形式，可以將代碼：

np.savetxt("param_values.txt",?param_values)#?將參數(shù)變化保存，其實(shí)是各參數(shù)范圍內(nèi)的sobol抽樣

生成的抽樣帶入自己的系統(tǒng)，然后根據(jù)自己需要生成對應(yīng)抽樣的目標(biāo)函數(shù)，將目標(biāo)函數(shù)放入同目錄下的 outputs.txt 文檔中，一行一個(gè)結(jié)果，然后用這個(gè)語句代替上面求Y[i]:

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

就是說提供了參數(shù)變化以及目標(biāo)函數(shù)變化，用SALib就可以求參數(shù)靈敏度了。

總結(jié)：

matlab我用的sobol生成的抽樣和別人的不一樣，不知道為什么，這個(gè)是造成與python計(jì)算不一樣的一個(gè)原因吧。但差別不大。

總結(jié)

以上是生活随笔為你收集整理的python如何做敏感度分析_1stopt、matlab和python用morris、sobol方法实现参数敏感性分析...的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。