Calculate confidence interval of data sets.

I have n=10 and calculated random values with two variables (SYSTEM.m_BW and SYSTEM.w_L). 2 ODEs and a loop is mentioned to run for every 10 iterations (for i=1:10) and plotted the graph for the same. Now I want to calculate mean of all the lines from the plot of simulated data (x(:,2)+x(:,1)) mentioned in the loop, then save x & t (time span) in a matrix to calculate and plot confidence interval from the simulated data given in the code (x(:,2)+x(:,1)) for every x= 1 to 10 lines. How do I calculate the mean of all the lines in the plot and how do I calculate confidence interval based on the mean and standard deviation of the data points? I want the graph something as given in attached graph with confidence interval of 95% calculated from mean shown with the upper and lower boundaries.

Could someone please help me with this, it would be really helpful and appreciated. Following is the given code, i have worked on so far:

function toy_example

clear all;

close all;

clc;

% Body weight and organ weight fraction [Min Max] (random numbers for 10 mice)
%Min=a
%Max=b
%r=a+(b-a)*r(n,1);
n=10;
SYSTEM.m_BW= 24+(28-24).*rand(10,1);                     % Body weight [g]
SYSTEM.w_L= 0.0417+(0.0681-0.0417).*rand(10,1);          % Liver [g/g]
SYSTEM.w_S= 0.019+(0.0051-0.0019).*rand(10,1);           % Spleen [g/g]
SYSTEM.w_K= 0.015+(0.0184-0.015).*rand(10,1);            % Kidney [g/g]
SYSTEM.w_Lu= 0.0065+(0.0081-0.0065).*rand(10,1);         % Lung [g/g]     
SYSTEM.w_BR= 0.0139+(0.0191-0.0139).*rand(10,1);         % Brain [g/g]
SYSTEM.w_Blood= 0.049+(0.049-0.049).*rand(10,1);         % Blood [g/g]
SYSTEM.w_Plasma= 0.029+(0.029-0.029).*rand(10,1);        % Plasma [g/g]
EXP.t=[0 6 12 18 24];                                    %h (time for blood pharmacokinetics)
DRUG.w_Au_iv = [15E-03 183E-03 693E-03 1058E-03];        %[mg/g]
DRUG.m_Au_iv = DRUG.w_Au_iv.*SYSTEM.m_BW(1);             %[mg]
EXP.Blood_15=(([16.32 3.8 2.87 2.33 2.06]/100).*SYSTEM.m_BW.*SYSTEM.w_Plasma.*DRUG.w_Au_iv(1).*SYSTEM.m_BW)./(SYSTEM.m_BW.*SYSTEM.w_Plasma);       
EXP.Blood_183=(([31.3 1.48 0.98 0.84 0.8]/100).*SYSTEM.m_BW.*SYSTEM.w_Plasma.*DRUG.w_Au_iv(2).*SYSTEM.m_BW)./(SYSTEM.m_BW.*SYSTEM.w_Plasma);
EXP.Blood_1058=(([57.82 0.42 0.35 0.35 0.35]/100).*SYSTEM.m_BW.*SYSTEM.w_Plasma.*DRUG.w_Au_iv(4).*SYSTEM.m_BW)./(SYSTEM.m_BW.*SYSTEM.w_Plasma);

figure
hold on
for i=1:10
       tspan=[0:24];
       MAT=[SYSTEM.m_BW(i) SYSTEM.w_L(i)]';
       x0=zeros(2,1);
       tic;
       [t,x]=ode15s(@ode_toy, tspan, x0, [ ], MAT);
       toc;
       plot(t,(x(:,2)+x(:,1)),'-.','Color', 'blue', 'LineWidth', 2)
end
end

function dxdt = ode_toy(t,x,SYSTEM)
%% Split State Vector
m_BW=SYSTEM(1); 
w_L=SYSTEM(2); 
m_Au_A=x(1);
m_Au_V=x(2);

%% ODE system
% Dosing
dm_Au_A_dt=m_BW - w_L*m_Au_A;
dx(1) = dm_Au_A_dt;
dm_Au_V_dt= w_L - m_BW*m_Au_V;
dx(2) = dm_Au_V_dt; 
dxdt = dx(:);
end

[dpb]

... 
[t,x]=ode15s(@ode_toy, tspan, x0, [ ], MAT);
...

overwrites t,x each pass through the loop.

Two choices -- either save them as cell arrays or vectors in a 2D array (which may not work without extra effort if ODE15S returns different sample times/lengths for the results) or retrieve the XData,YData arrays from the figure when it's done.

The mean/std dev of the mean are relatively easy to calculate from either excepting for the possibility that the time points/lengths don't match up. If that is the case, you'll need to fit an interpolating spline to the results to evaluate at a set of specific times and use that for mean().

Sounds more complicated than it is...give it a go! :)

[Rajvi Amle]

Thank you for your response @dpb. However I am new on Matlab, could you please show me an example of saving x,t results in a cell array or a vector? I used a code as C=cell(10,1); before the for loop as shown below.

C=cell(10,1);
figure
hold on
for i=1:10
       tspan=[0:24];
       MAT=[SYSTEM.m_BW(i) SYSTEM.w_L(i)]';
       x0=zeros(2,1);
       tic;
       [t,x]=ode15s(@ode_toy, tspan, x0,[], MAT);
       toc;
       plot(t,(x(:,2)+x(:,1)),'-.','Color', 'blue', 'LineWidth', 2)
end

Also, I didn't understand what you explained to derive confidence interval. If possible, could you please show that too with an example?

[dpb]

I'm really pressed for time, so the following is untested "air code" typed in here...hopefully it'll be enough to give you the idea...

t=cell(1,10); y=t;  % preallocate the output cell arrays as row to hold column vectors
figure
hold on
tspan=[0:24];       % invariant, don't need in loop
for i=1:10
  MAT=[SYSTEM.m_BW(i) SYSTEM.w_L(i)]';
  x0=zeros(2,1);
  [t{i},x]=ode15s(@ode_toy, tspan, x0,[], MAT); % assign output to cell arrays
  y{i}=(x(:,2)+x(:,1);                          % compute, save y also 
  plot(t{i},y{i},'-.','Color', 'blue', 'LineWidth', 2)
end

IFF the returned t vector is the same for every iteration, then you can convert it to a single vector as would all be redundant and turn y into a normal 2D array.

Then the mean, sd of mean would be

mny=mean(cell2mat(y),2);        % convert cell array to matrix, mean by row
sdmny=mny/size(y,2);            % std dev of mean(y) is mean(y)/nObs