Problem with calculating Sobol Indices, some sort of background error is making answers that should be zero non-zero.

Hi all,

I'm currently trying to calculate Sobol indices for a set of parameters in Matlab. It involves manipulation of a couple of matrices and then using the results to calculate the indices. We expect the indices to be between 0 and 1, where 1 indicates a very sensitive parameter.

However, when I run the program I seem to get some background noise. For example the latest output was

s =

      0.2761    0.3196    0.2730    0.6570    0.2771    0.2780    0.2666    0.4497    0.2790    0.2786

From other analysis we have done, we expect the 4th and 8th parameter to be most sensitive. As you can see here, they're more sensitive than the others but I would expect the rest to be around 0, however there is some error where each parameter is about 0.26 higher than it should be. Sobol indices should add up to be less than or equal to 1 when finished so I know this answer to be incorrect.

Basically, if all these answers had 0.2666 subtracted from them I feel they would be correct, but something in the code that I've been unable to spot is not working correctly.

Any help with this would be excellent, perhaps it's something simple I've missed because I've been staring at it too long! Here's the code anyway

function Sobol
clear all
tic;

%Declare Parameters
parameters=[ 2  1.374  1.6e4  2.99  2.26 2.65  2  0.315  0.0315 110];

%Specify max parameter level
parameter_max=3*parameters;

%Specify min parameter level
parameter_min=parameters*0.333;

%Calculate no. of parameters
pn=length(parameter_max);

%no. of tests
n = 500;

%Generate hypercube sampled matrix
X = lhsdesign(n, pn);

%Generate random parameters in specified parameter space
parameter_sets=exp(log(ones(n, 1)*parameter_min)+((log(ones(n, 1)*parameter_max-log(ones(n,1)*parameter_min)))).*X);

%Generate A and B
A = parameter_sets(1:(n/2), :);
B = parameter_sets((n/2 + 1):n,:);

%Initial Conditions
ic1 = [4*1.32e1 2.65/2 0 1.375/2 0];

%Time steps
tg = [0 5];

%Options
options = odeset('RelTol',1e-6,'AbsTol',1e-6);

%Find f(A)
for i=1:(n/2)

%Read parameters from first row
a=num2cell(A(i, :));
[dgsh bgsh kgsh kg ks bs ds k450 kn kpsh]=deal(a{:});

%Run Model
[t y]=ode15s(@model,tg,ic1,options);

%Save fA as (N x 1) matrix
fA(i) = y(end,5);
end

%Find f(B)
%for i=1:N

%Read parameters from first line
%b=num2cell(B(i, :));
%[k1 k2 k3 k4]=deal(b{:});

%Run Model
%[t y]=ode15s(@model,tg,ic1,options);

%Save fB as N x pn matrix
%fB(i) = max(y(:,2));

%end

%Create Ci

for i=1:pn
C = B;
C(:,i) = A(:,i);

for j=1:(n/2)

      %Read parameters from first line
      c=num2cell(C(j, :));
      [dgsh bgsh kgsh kg ks bs ds k450 kn kpsh]=deal(c{:});

      %Run Model
      [t y]=ode15s(@model,tg,ic1,options);

      fC(i,j) = y(end,5);  

end

end

%Calculate Si

%Calculate term 2 of Si equation      
    b = ((1/n)*sum(fA))^2

%Calculate term 3 of Si equation  
    c = (1/n)*sum(fA.^2)

for i = 1:pn

      %Calculate term 1 of Si equation 
      a = (1/n)*sum(fA.*fC(i,:))

      %Calculate Si    
      s(i) = (a - b)/(c - b)
  end  

sum(s)
toc;

%Model
function yp=model(t,y)
P=y(1);
S=y(2);
N=y(3);
G=y(4);
T=y(5);

yp1 = -ks*S*P - kg*P - k450*P + kn*N;
yp2 = -ks*S*P + bs - ds*S;
yp3 = k450*P - kn*N - kgsh*N*G - kpsh*N;
yp4 = -kgsh*N*G + bgsh - dgsh*G;
yp5 = kpsh*N;

yp=[yp1; yp2; yp3; yp4; yp5];
end

end

Any help would be much appreciated,

Thank you.