How to generate and store 3 sets of values that change every time the loop is ran?

My code is supposed to run inside two loops. The first loop generates two random variables and from them it calculates the values over time for n1, n2 and n3 from time ranging from 1 to 20, with step tau_w (which also changes according to the random numbers generated). Once I am able to find this initial three sets of values for n1,n2 and n3 over time, I have to repeat the entire process 763 times. That would give me 763 sets for n1, n2 and n3 over time, from which I will obtain the mean. The issues I am having are that I can't figure out how to store the values obtained every repetition, it seems to only stores the last ones, and also the values are not changing, it's like they are just repeating 763 times (I'm attaching a figure of the stored values of N2 that I am currently getting below). How can I do that?

Thanks in advance!

clc

clear all

tif = 0.00218;%day^-1

tif_2 = 0.44036; %day^-1

pop = 763; % Total number of repetitions

% Initial values of

n1(1)=762; %

n2(1)=1; %

n3(1)=0; %

rep = 1; % number of repetitions

tf= 20; %days

t2(1) = 0;

i=1;

N1 = [];

N2 = [];

N3 = [];

while rep <= pop

while t2(end)< tf

% Step 3: determine length of t_w. Generate random number r1, from uniform

% distribution in [0,1], then calculate t_w for a system's transition at

% state (n1,n2,n3).

r1(i)=rand(1);

t_w(i) = (1/( -( n1(i)*n2(i)*tif + n2(i)*tif_2 ) ))* log(1-r1(i));

% Step 4: Update computer clock:

t2(i+1) = t2(i) + t_w(i);

% Step 5: Calculate transition probabilities that system will transfer from state n

%to other states Qi:

Q1 = n1(i)*n2(i)*tif/(n1(i)*n2(i)*tif + n2(i)*tif_2);

Q2 = n2(i)*tif_2/(n1(i)*n2(i)*tif + n2(i)*tif_2);

%Generate another random number r2 from uniform distribution in[0,1]:

r2(i) = rand(1);

%Determine transition type based on:

if r2(i) <=Q1

n1(i+1) = n1(i) -1;

n2(i+1) = n2(i) +1;

n3(i+1) = n3(i);

else

n1(i+1) = n1(i);

n2(i+1) = n2(i)-1;

n3(i+1) = n3(i)+1;

end

i = i +1;

%Step 6: Repeat steps 3–5 until the total time exceeds Tf;

%this terminates one replication of simulation.

%Step 7: Repeat steps 2-6 763 times and store the values to obtain the average

end

rep = rep+1;

N1 = [N1; n1]; % Store values of n1 at each repetition

N2 = [N2; n2]; % Store values of n2 at each repetition

N3 = [N3; n3]; % Store values of n3 at each repetition

end

What I am getting for N2 (similar to N1 and N3) every repetition:

[Image Analyst]

I'm looking for a for loop.

for rep = 1 : pop

Why are you using a while loop?

[Rebeca Hannah Oliveira]

With the while loop I actually could get a 763 x (time) table. With for, I only get a 763 x 3. So I guessed a while loop would be the way to proceed.

(with while loop) (with for loop)

Answer 1

[Image Analyst]

You need to reinitialize some iterators at the beginning of each iteration. Also, you need to expand the width of N1, N2, and N3 if the next n1, n2, or n3 are wider than the old N1, N2, or N3 matrices. Try this and see if it does what you want:

clc;    % Clear the command window.
fprintf('Beginning to run %s.m ...\n', mfilename);
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 18;

tif = 0.00218;		%day^-1
tif_2 = 0.44036;	%day^-1
pop = 763; % Total number of repetitions
% Initial values of
n1(1)=762; %
n2(1)=1; %
n3(1)=0; %
rep = 1; % number of repetitions
tf = 20; %days
t2 = 0;
N1 = [];
N2 = [];
N3 = [];
for rep = 1 : pop
    fprintf('Starting iteration #%d of %d...\n', rep, pop)
    % Re-initialize i and t2.
    i = 1;
    t2 =  0;
    while t2(end)< tf
        % Step 3: determine length of t_w.
        % First generate random number r1, from uniform distribution in [0,1].
        r1(i)=rand(1);
        % Then calculate t_w for a system's transition at state (n1,n2,n3).
        t_w(i) = (1/( -( n1(i)*n2(i)*tif + n2(i)*tif_2 ) ))* log(1-r1(i));
        % Step 4: Update computer clock:
        t2(i+1) = t2(i) + t_w(i);
        % Step 5: Calculate transition probabilities that system will transfer from state n
        %to other states Qi:
        Q1 = n1(i)*n2(i)*tif/(n1(i)*n2(i)*tif + n2(i)*tif_2);
        Q2 = n2(i)*tif_2/(n1(i)*n2(i)*tif + n2(i)*tif_2);
        % Generate another random number r2 from uniform distribution in[0,1]:
        r2(i) = rand(1);
        % Determine transition type based on:
        if r2(i) <= Q1
            n1(i+1) = n1(i) -1;
            n2(i+1) = n2(i) +1;
            n3(i+1) = n3(i);
        else
            n1(i+1) = n1(i);
            n2(i+1) = n2(i)-1;
            n3(i+1) = n3(i)+1;
        end
        i = i +1;
        % Step 6: Repeat steps 3–5 until the total time exceeds Tf;
        % this terminates one replication of simulation.
        % Step 7: Repeat steps 2-6 763 times and store the values to obtain the average
    end
    fprintf('  n1 is %d elements long.  n2 is %d elements long.  n3 is %d elements long.\n', length(n1), length(n2), length(n3))
    % Expand matrices if new n1, n2, n3 are wider than N1, N2, or N3.
    if rep > 1
        if size(N1, 2) < length(n1)
            N1(end, length(n1)) = 0;
        end
        if size(N2, 2) < length(n2)
            N2(end, length(n2)) = 0;
        end
        if size(N3, 2) < length(n3)
            N3(end, length(n3)) = 0;
        end
    end
    % Tack on n1, n2, and n3 as rows at the bottom of existing N1, N2, and N3.
    N1 = [N1; n1]; % Store values of n1 at each repetition
    N2 = [N2; n2]; % Store values of n2 at each repetition
    N3 = [N3; n3]; % Store values of n3 at each repetition
    % Plot
    yyaxis left
    cla;
    plot(t_w, 'r-')
    ylabel('t_w', 'Interpreter', 'none')
    hold on;
    yyaxis right
    cla
    plot(t2, 'b-')
    grid on;
    ylabel('t2')
    drawnow;
end

[Rebeca Hannah Oliveira]

Thank you! It is exactly what I wanted!

Thanks.