Using Latin Hyper Cube Sampling for Acquiring Multiple Replicates of an ODE

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

Answer 1

[Aditya]

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.