首先看拋物線函數：
現在我取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方法局部：
名稱 靈敏度指數 靈敏度指數(%)
a 3.916E18 89.1113892365457
b 4.125E17 9.38673341677096
c 6.6E16 1.50187734668335
morris方法全局：
名稱 靈敏度指數 靈敏度指數(%)
a 2.00290323137447E18 81.4890482414753
b 2.10289523274038E17 8.55572692595605
c 2.44687506069835E17 9.95522483256862
偏微分方法局部：
名稱 靈敏度指數 靈敏度指數(%)
a 3.99432000000567E18 89.1113892365577
b 4.20750000000213E17 9.38673341676365
c 6.73199999998753E16 1.50187734667864
偏微分方法全局：
名稱 靈敏度指數 靈敏度指數(%)
a 3.98876054886661E18 82.084202339751
b 4.20542269921853E17 8.65428655186941
c 4.50049450186404E17 9.26151110837963
- 下面用matlab寫sobol方法進行分析：
1 % sobol 參數敏感性分析 2 %算法參考： 3 % csdn : https://blog.csdn.net/xiaosebi1111/article/details/46517409 4 % wiki： https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis 具體參考文獻這里也有 5 %運行環境 matlab2016b 6 %作者 blzhu@buaa.edu.cn 2020年6月7日 7 %% 初始化 8 clc; 9 clear all; 10 close all; 11 %% 設定：給定參數個數和各個參數的范圍 12 13 % % 1-自定義子函數1 14 % 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-自定義子函數2 21 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-自定義子函數3 28 % 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 34 35 %% 開始計算 36 numnPop=numel(nPop); 37 SAll=zeros(numnPop,D+1);%分別是各參數的敏感度，最后一列是各參數敏感度之和 38 STAll=zeros(numnPop,D+1); 39 for i=1:numnPop 40 [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 end 46 %% 繪圖 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','>','<','p','h'];%13種標記 一般標記也不一樣; 字符數組，每個字符占1個位置 49 linestyle=[string('-'),string('--'),string(':'),string('-.'),string('none')]; 50 useL=1; 51 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 on 56 end 57 title('Sobol-S'); 58 whichPara=sprintfc('%g',repmat(1:D+1,1,2));%把數字數組轉化成字符串數組 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 on 66 end 67 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 train 79 sound(y,Fs) 80 81 82 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 end 89 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 end 99 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 end 108 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:nPop 121 r=p(i,:); 122 r=VarMin+r.*(VarMax-VarMin); 123 R=[R; r]; 124 end 125 % plot(R(:,1),'b*') 126 % 拆分為A B 127 A=R(:,1:D);% 每行代表一組參數，其中每列代表每組參數的一個參數；行數就代表共有幾組參數 128 B=R(:,D+1:end); 129 % 根據A B 產生矩陣AB 130 AB=zeros(nPop,D,D); 131 for i=1:D 132 tempA=A; 133 tempA(:,i)=B(:,i); 134 AB(1:nPop,1:D,i)=tempA; 135 end 136 %% 求各參數解 137 YA=zeros(nPop,1);% 解 138 YB=zeros(nPop,1); 139 YAB=zeros(nPop,D);%分別代表YAB1,YAB2,YAB3，YAB(:,D)就代表YABD 140 for i=1:nPop 141 YA(i)=myFunction(A(i,:)); 142 YB(i)=myFunction(B(i,:)); 143 for j=1:D 144 YAB(i,j)=myFunction(AB(i,:,j)); 145 end 146 end 147 %% 根據一階影響指數公式： 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:D 155 for j=1:nPop 156 VarX(i)=VarX(i)+YB(j)*(YAB(j,i)-YA(j)); 157 end 158 VarX(i)=1/nPop*VarX(i);% 蒙特卡羅估計量 159 S(i)=VarX(i)/VarY; 160 end 161 162 %% 總效應指數 163 EX=zeros(D,1); 164 ST=zeros(D,1); 165 for i=1:D 166 for j=1:nPop 167 EX(i)=EX(i)+(YA(j)-YAB(j,i))^2; 168 end 169 EX(i)=1/(2*nPop)* EX(i);% 蒙特卡羅估計量 170 ST(i)=EX(i)/VarY; 171 end 172 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包分析
python中也有參數靈敏度分析的包SALib（https://salib.readthedocs.io/en/latest/index.html），這個包支持以下算法：
- Sobol Sensitivity Analysis ([Sobol 2001], [Saltelli 2002], [Saltelli et al. 2010])
- 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])
- Random Balance Designs - Fourier Amplitude Sensitivity Test (RBD-FAST) ([`Tarantola et al. 2006 <https://hal.archives-ouvertes.fr/hal-01065897/file/Tarantola06RESS_HAL.pdf`_], [Elmar Plischke 2010], [Tissot et al. 2012])
- 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 import saltelli 4 from SALib.analyze import sobol 5 from SALib.test_functions import Ishigami 6 import numpy as np 7 import math 8 from SALib.plotting.bar import plot as barplot 9 import matplotlib.pyplot as plot 10 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)# 不管用哪個方法計算y,這個要有 29 np.savetxt("param_values.txt", param_values)# 將參數變化保存，其實是各參數范圍內的sobol抽樣 30 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]) 39 40 #1-自定義-2 41 Y = np.zeros([param_values.shape[0]]) 42 for i, X in enumerate(param_values): 43 tarr=np.arange(-5,6,1); 44 yerror=0.0; 45 for t in tarr: 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) 54 55 # 3-SALib自帶測試函數 56 # Y = Ishigami.evaluate(param_values) 57 ## np.savetxt("outputs.txt", Y)#將因變量變化結果保存 58 59 Si = sobol.analyze(problem, Y ,print_to_console=True) 60 print()#自動輸出S1（單個參數對因變量的影響）、ST（考慮各個變量相互影響）和S2(兩兩參數之間影響)，需有,print_to_console=True 61 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'])# 總效應指數 68 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) 代替
python3.8.1+SALib1.3 use morris
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 import saltelli 4 from SALib.analyze import sobol 5 from SALib.analyze import morris 6 from SALib.test_functions import Ishigami 7 import numpy as np 8 import math 9 from SALib.plotting.bar import plot as barplot 10 import matplotlib.pyplot as plot 11 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)# 不管用哪個方法計算y,這個要有 30 np.savetxt("param_values.txt", param_values)# 將參數變化保存，其實是各參數范圍內的sobol抽樣 31 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]) 40 41 #1-自定義-2 42 Y = np.zeros([param_values.shape[0]]) 43 for i, X in enumerate(param_values): 44 tarr=np.arange(-5,6,1); 45 yerror=0.0; 46 for t in tarr: 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) 55 56 # 3-SALib自帶測試函數 57 # Y = Ishigami.evaluate(param_values) 58 ## np.savetxt("outputs.txt", Y)#將因變量變化結果保存 59 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 63 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'])# 總效應指數 70 71 #繪圖 https://zhuanlan.zhihu.com/p/137953265 72 Si_df = Si.to_df() 73 barplot(Si_df[0]) 74 plot.show()
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計算不一樣的一個原因吧。但差別不大。
靈敏度分析的概念在此處沒有詳細講，可以參考：https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis