MATLAB 实现sobol参数敏感性分析

% sobol 参数敏感性分析
%参考：
% csdn : https://blog.csdn.net/xiaosebi1111/article/details/46517409
% wiki： https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis
%运行环境 matlab2016b
%作者 blzhu@buaa.edu.cn 2020年6月7日
%% 初始化
clc;
clear all;
close all;
%% 设定：给定参数个数和各个参数的范围
D=3;% 维度3,几个参数
M=D*2;% 
nPop=4;% 采样点个数,也就是参数水平数 ,取大了好，比如4000，但慢
VarMin=[0 0 0 ];%各个参数下限
VarMax=[1 1 1];%各个参数上限
%% 产生所需的各水平参数
VarMin=[VarMin,VarMin];
VarMax=[VarMax,VarMax];
p= sobolset(M);% https://www.cnblogs.com/zhubinglong/p/12260292.html
% R=p(1:nPop,:);% 我只用前nPop个
R=[];
for i=1:nPop
    r=p(i,:);
    r=VarMin+r.*(VarMax-VarMin);
    R=[R; r];
end
% plot(R(:,1),'b*')
% 拆分为A B
A=R(:,1:D);% 每行代表一组参数，其中每列代表每组参数的一个参数；行数就代表共有几组参数
B=R(:,D+1:end);
% 根据A B 产生矩阵AB
AB=zeros(nPop,D,D);
for i=1:D
    tempA=A;
    tempA(:,i)=B(:,i);
    AB(1:nPop,1:D,i)=tempA;
end
%% 求各参数解
YA=zeros(nPop,1);% 解
YB=zeros(nPop,1);
YAB=zeros(nPop,D);%分别代表YAB1,YAB2,YAB3，YAB(:,D)就代表YABD
for i=1:nPop
    YA(i)=myfun(A(i,:));
    YB(i)=myfun(B(i,:));
    for j=1:D
        YAB(i,j)=myfun(AB(i,:,j));
    end
end
%%  根据一阶影响指数公式：
VarX=zeros(D,1);% S的分子
S=zeros(D,1);

% 0: 估算基于给定样本的方差(EXCEL var.p) ;   1:计算基于给定的样本总体的方差(EXCEL var.p())
% var([2.091363878    1.110366059    3.507651769    1.310950363    2.091363878    3.507651769    1.110366059    1.7066512],1);
VarY=var([YA;YB],1);% S的分母。 计算基于给定的样本总体的方差(EXCEL var.p())
for i=1:D
    for j=1:nPop
         VarX(i)=VarX(i)+YB(j)*(YAB(j,i)-YA(j));
    end
     VarX(i)=1/nPop*VarX(i);
     S(i)=VarX(i)/VarY;
end

%% 总效应指数
EX=zeros(D,1);
ST=zeros(D,1);
for i=1:D
    for j=1:nPop
         EX(i)=EX(i)+(YA(j)-YAB(j,i))^2;
    end
      EX(i)=1/(2*nPop)* EX(i);
     ST(i)=EX(i)/VarY;
end
disp('一阶影响指数：S');
disp(S);
disp('总效应指数：ST');
disp(ST);
disp('success!');


%% 子函数 matlab2016之前不支持子函数写在同一个m文档中
function y=myfun(x)
y=sin(x(1))+7*(sin(x(2)))^2+0.1*x(3)^4*sin(x(1));
end