Using Latin Hyper Cube Sampling for Acquiring Multiple Replicates of an ODE
17 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Hi All,
See below code for producing 1,000 replicates of an ODE using random parameter values within a specified range. I am hoping to replace this method for acquring random parameter values with Latin Hypercube Sampling. I have been able to design a Latin Hypercube with the "lhdesign" function, but I am confused regarding how I should reference it when defining my parameters. I can see that lhdesign defaults to producing values in a unit interval from 0 to 1, but I figured I could use the same method I am currently employing with the "rand" function to specify my desired ranges for each parameter. It should also be noted that I am hoping to sample a uniform distribution. Many thanks in advance for the help!
%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
plot(t, y(:,2));
end
0 Commenti
Risposte (1)
Aditya
il 19 Feb 2024
To use Latin Hypercube Sampling (LHS) for generating your parameter values instead of the uniform random numbers generated by rand, you need to scale the LHS samples to match the range of each parameter.
Here's how you can modify your code to use LHS for your parameters:
% 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 samples generated for 11 parameters
LHS = lhdesign(n, 11);
for k = 1:n
% Initial conditions
n = 16956;
y20 = 10;
y30 = 0;
% Parameters sampled using LHS instead of "rand"
p = (407-341) * LHS(k, 1) + 341;
d = (1/42-1/50) * LHS(k, 2) + 1/50;
yc = (4.923*10^-3-8.8*10^-5) * LHS(k, 3) + 8.8*10^-5;
yi = (4.923*10^-3-8.8*10^-5) * LHS(k, 4) + 8.8*10^-5;
a = (1/30-1/120) * LHS(k, 5) + 1/120;
r = (0.30-.10) * LHS(k, 6) + 0.1;
i = (1/4-1/15) * LHS(k, 7) + 1/15;
bc = (1.5*10^-3-1.5*10^-5) * LHS(k, 8) + 1.5*10^-5;
bi = (1.5*10^-3-1.5*10^-5) * LHS(k, 9) + 1.5*10^-5;
c = (50-5) * LHS(k, 10) + 5;
tr = (0.2-0.01) * LHS(k, 11) + 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
plot(t, y(:,2));
end
In this edited version of your code, each parameter is scaled based on the corresponding column of the LHS matrix. Each row in the LHS matrix corresponds to a different replicate, ensuring that the parameter values are sampled uniformly across their ranges. This method provides a more evenly distributed sampling of the parameter space compared to random sampling and is particularly useful when the number of samples is limited.
0 Commenti
Vedere anche
Categorie
Scopri di più su Industrial Statistics in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!