Generating an average curve from 1,000 replicates of an ODE

Hi All,

I am currently attmepting to generate an average curve from 1,000 randomly iterated replicates of a 3rd order ODE. All of the randomly generated curves are unfiform in length, and I am only trying to assess the average for one of the three functions (y(2)). I have never worked with indexing nested cells, and I am having trouble writing a for loop that adds a single y value at each time point across all the replicates and divides by 1000 to acquire the average at each of the 4,991 points that comprise each curve. Its entirley possible there may be a better way to do this, but I have included my code below to show what I have been trying. Unfortunatley my code only seems to generate a cell array of length 4992 occupied entirley by empty []. I am very new to coding, so my for loop is probably a little wonky. Sorry for the foolishness, and thanks for taking a look!

%Total number of replicates

n=1000;

%Generating cells to store results

result=cell(n,1);

%Figure hold to plot all 1000 replicates on a single graph

figure; hold on

%Latin hypercube with 1000 values generated for 11 parameters

LHS=lhdesign(1000,11);

for k=1:n

%Initial conditions

n=16956;

y20=10;

y30=0;

%Parameters I am hoping to sample randomly using LHS instead of "rand"

p=(407-341)*rand+341;

d=(1/42-1/50)*rand+1/50;

yc=(4.923*10^-3-8.8*10^-5)*rand+8.8*10^-5;

yi=(4.923*10^-3-8.8*10^-5)*rand+8.8*10^-5;

a=(1/30-1/120)*rand+1/120;

r=(0.30-.10)*rand+0.1;

i=(1/4-1/15)*rand+1/15;

bc=(1.5*10^-3-1.5*10^-5)*rand+1.5*10^-5;

bi=(1.5*10^-3-1.5*10^-5)*rand+1.5*10^-5;

c=(50-5)*rand+5;

tr=(0.2-0.01)*rand+0.01;

%Time step and length

tspan=[(1:0.1:500)];

%Function definition

[t,y] = ode23s(@(t,y) [p*(1-yc-yi) + a*y(2) + tr*y(3) - (c*bc*y(2)*y(1))/n - (c*bi*y(3)*y(1))/n - d*y(1); p*yc + (c*bc*y(2)*y(1))/n + (c*bi*y(3)*y(1))/n - a*y(2) - r*i*y(2) - d*y(2); p*yi + r*i*y(2) + - d*y(3)-tr*y(3);], tspan, [n y20 y30]);

%Storing results in cells

result{k} = y;

%Graphing all replicates for y(2). This is a model of community transmitted MRSA in the LA county jail. Y(2) represents the total population of infected individuals

plot(t, y(:,2));

end

%Total number of replicates

n=1000;

%Total number of points that comprise each curve

j=4991;

%Generating a cell array to house the average curve

y2_mean=cell(j,1);

%Iterating across all 1000 replicates

for k=1:n

e=result{n,1};

%Iterating across all 4991 pointsa

for i=1:j

y2_mean{j} = (sum(e(j,2))/1000);

end

end

Answer 1

[Mathieu NOE]

hello

I tweaked a bit your code. Your results can be stored in cells but as your y array has always the same dimensions, I found even simpler to store the result as a numerical array , and then use mean function at the end of the loop

the mean value appears now as a thick red line

see some other minor comments in the code :

clc
clearvars

%Total number of replicates
nr=1000; %changed the name from n to nr to avoid conflict with other "n" below (n=16956;)

%Time step and length (this is fixed so no need to put it in the for loop)
tspan=[(1:0.1:500)];
samples = length(tspan);

%Generating cells to store results
% result=cell(n,1);
result=zeros(samples,nr); % simple numerical 2D array ; preallocation 

%Figure hold to plot all 1000 replicates on a single graph
figure; hold on

%Latin hypercube with 1000 values generated for 11 parameters
%LHS=lhdesign(1000,11); % commented this to avoid error (not needed anyway
%here)

%Initial conditions (this is fixed so no need to put it in the for loop)
n=16956; % 
y20=10;
y30=0;
    
for k=1:nr

    %Parameters I am hoping to sample randomly using LHS instead of "rand"
    p=(407-341)*rand+341;
    d=(1/42-1/50)*rand+1/50;
    yc=(4.923*10^-3-8.8*10^-5)*rand+8.8*10^-5;
    yi=(4.923*10^-3-8.8*10^-5)*rand+8.8*10^-5;
    a=(1/30-1/120)*rand+1/120;
    r=(0.30-.10)*rand+0.1;
    i=(1/4-1/15)*rand+1/15;
    bc=(1.5*10^-3-1.5*10^-5)*rand+1.5*10^-5;
    bi=(1.5*10^-3-1.5*10^-5)*rand+1.5*10^-5;
    c=(50-5)*rand+5;
    tr=(0.2-0.01)*rand+0.01;

    %Function definition
    [t,y] = ode23s(@(t,y) [p*(1-yc-yi) + a*y(2) + tr*y(3) - (c*bc*y(2)*y(1))/n - (c*bi*y(3)*y(1))/n - d*y(1); p*yc + (c*bc*y(2)*y(1))/n + (c*bi*y(3)*y(1))/n - a*y(2) - r*i*y(2) - d*y(2); p*yi + r*i*y(2) + - d*y(3)-tr*y(3);], tspan, [n y20 y30]);
    %Storing results in cells
%     result{k} = y(:,2);
    result(:,k) = y(:,2); % simple num array
    
    %Graphing all replicates for y(2). This is a model of community transmitted MRSA in the LA county jail. Y(2) represents the total population of infected individuals
    plot(t, y(:,2),'--','Linewidth',0.75);
end

% add mean curve
y2_mean = mean(result,2);
plot(t, y2_mean,'r','Linewidth',3);
ylim([0 3*max(y2_mean)]);
hold off

[Robert Pearhill]

This is awsome, thank you so much!

[Mathieu NOE]

My pleasure !

Answer 2

[Jan]

I cannot run your code due to a missing function lhdesign, but I guess you want:

avg = mean([result{:}(:, 2)], 2);

[Robert Pearhill]

Sorry about that, I forgot to remove the lhdesign before posting! Thank you for the response!