敏感性分析 sobol matlab 代码

敏感性分析 sobol matlab 代码

start之前是写出 sobol 序列，后面直接运行就行， 如有错误请指正 ，谢谢。
下面展示一些 内联代码片。

// A code block
% 样本数为4(N = 4)，自变量数目为3(D = 3)。 
% 1、 生成N * 2D（即4行6列）的样本矩阵。这个就是我们Sobol sequence做的事情。这边我们生成的矩阵为：
clc;
clear all;
close all;
M=2*3;% 维度，几个参数
nPop=5;% 样本数为4
VarMin=[0,  0 ,  0.0, 0,   0 ,   0.0];%各个参数下限
VarMax=[1.0, 1.0, 1.0,1.0,   1.0, 1.0];%各个参数上限
p = sobolset(M);
% R=p(1:nPop,:);% 我只用前nPop个
R=[];
for i=1:nPop
    r=p(i,:);
    r=VarMin+r.*(VarMax-VarMin);
    R=[R; r];
end
plot(R(:,2),'b*');

R2=R(2:5,:);
%===============================start ===========
R22=[0.5 0.5 0.5 0.5 0.5 0.5;
     0.75 0.25 0.25 0.25 0.75 0.75;
     0.25 0.75 0.75 0.75 0.25 0.25
     0.375 0.375 0.625 0.875 0.375 0.125];
A=R22( : ,1:3);
B=R22( : ,end-2:end);
%---------------------------
A1=A; % 记录A B 初始值的中间变量
B1=B;
A(: ,1 )=B(: ,1 );
AB_1=A ;
%---------------------------
A=A1;
B=B1;
A(: ,2 )=B(: ,2 );
AB_2=A ;
%---------------------------
A=A1;
B=B1;
A(: ,3 )=B(: ,3 );
AB_3=A ;


% (D + 2) * N (即20)组x1、x2、x3输入数据，因此我们将有20组Y值。
% 样本数为4(N = 4)，自变量数目为3(D = 3)。 
% 构造了A、B、AB1、AB2、AB3这五个矩阵
syms x1 x2 x3 Y
% Y=sin(x1)+7*sin(x2)^2+0.1*x3^4*sin(x1);
for i=1:4
x11=A1(i,1);x12=A1(i,2);x13=A1(i,3);
YA(i)=sin(x11)+7*sin(x12)^2+0.1*x13^4*sin(x11);
end
for i=1:4
x11=B1(i,1);x12=B1(i,2);x13=B1(i,3);
YB(i)=sin(x11)+7*sin(x12)^2+0.1*x13^4*sin(x11);
end
for i=1:4
x11=AB_1(i,1);x12=AB_1(i,2);x13=AB_1(i,3);
YAB1(i)=sin(x11)+7*sin(x12)^2+0.1*x13^4*sin(x11);
end
for i=1:4
x11=AB_2(i,1);x12=AB_2(i,2);x13=AB_2(i,3);
YAB2(i)=sin(x11)+7*sin(x12)^2+0.1*x13^4*sin(x11);
end
for i=1:4
x11=AB_3(i,1);x12=AB_3(i,2);x13=AB_3(i,3);
YAB3(i)=sin(x11)+7*sin(x12)^2+0.1*x13^4*sin(x11);
end
syms Mean_Y
Y=zeros(8,1);
Y(1:4,:)=YA';
Y(5:8,:)=YB';
Mean_Y=mean(Y);
Var_Y=var(Y);
%---------------------
% 4、 根据一介影响指数公式：==================
YAj=YA;
YBj=YB;
YAB1j=YAB1;
YAB2j=YAB2;
YAB3j=YAB3;
Varx1=0;
for i=1:4
Varx1=Varx1+YBj(i)*(YAB1j(i)-YAj(i));
Varx1
end

Varx1=Varx1/4;
s1=Varx1/Var_Y;

 % ————————————————
 % 原文链接：https://blog.csdn.net/xiaosebi1111/article/details/46517409

Varx2=0;
for i=1:4
Varx2=Varx2+YBj(i)*(YAB2j(i)-YAj(i));
Varx2
end

Varx2=Varx2/4;
s2=Varx2/Var_Y;
%--------------------
Varx3=0;
for i=1:4
Varx3=Varx3+YBj(i)*(YAB3j(i)-YAj(i));
Varx3
end

Varx3=Varx3/4;
s3=Varx3/Var_Y;
%--------------
%--------------------the following 全局影响指数==================
Ex_1=0;
for i=1:4
Ex_1=Ex_1+(YAB1j(i)-YAj(i))^2;
Ex_1
end

Ex_1=Ex_1/(2*4);
st1=Ex_1/Var_Y;
%--------------
%========================
Ex_2=0;
for i=1:4
Ex_2=Ex_2+(YAB2j(i)-YAj(i))^2;
Ex_2
end

Ex_2=Ex_2/(2*4);
st2=Ex_2/Var_Y;
%--------------
%========================

Ex_3=0;
for i=1:4
Ex_3=Ex_3+(YAB3j(i)-YAj(i))^2;
Ex_3
end

Ex_3=Ex_3/(2*4);
st3=Ex_3/Var_Y;

参考链接: link