% sobol 参数敏感性分析

%参考：

% csdn : https://blog.csdn.net/xiaosebi1111/article/details/46517409

% wiki： https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis

%运行环境 matlab2016b

%% 初始化

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