![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/826010/image.png)
Generating an average curve from 1,000 replicates of an ODE
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Robert Pearhill
el 7 de Dic. de 2021
Comentada: Mathieu NOE
el 8 de Dic. de 2021
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
0 comentarios
Respuesta aceptada
Mathieu NOE
el 7 de Dic. de 2021
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/826010/image.png)
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
2 comentarios
Más respuestas (1)
Jan
el 7 de Dic. de 2021
I cannot run your code due to a missing function lhdesign, but I guess you want:
avg = mean([result{:}(:, 2)], 2);
Ver también
Categorías
Más información sobre Particle & Nuclear Physics en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!