Help with running a Monte Carlo please...and another question :)

Hi there,

My code is broken up into the following segments:

Lines 7-54 : import data and preallocate for later matrices

Lines 61-148 : a bunch of variables that will be subject to a Monte Carlo simulations (normal distribution)

Lines 155 - 327 : a bunch of for-loops used to define the parameters used in the equation on Line 338.

Lines 328 and on are where I'm having difficulty... There are two issues.

1st: I'm asking the code to run the equation on L342 at sucessively later and later startpoints until the modeled serum level an the end of simulation no longer exceeds that of the measured serum. Once that is met, I'd like to save what that start time was of the last iteration (I took a stab at this on L350, but it's not working...). It should be noted that each column is a seperate person's data and each row is a 1 year time interval. So for example, in the first iteration, the model starts at 1970, if Serum(98,column) is greater than Serum_blood(column), then the simulation starts over at 1971, if it is still greater than in that iteration, run the model again starting at 1972, and so on. Ultimatly, I'm trying to, in a sense, back calculate the start of exposure to a contaminant. So, that is why I'm trying to determine/save the final start point for each participant (column). Hopefully this makes sense.

2nd: (and thank you for still giving this quandry the time of day/night!!) I want to run the monte carlo 1000 times.... so that I'm able to come up with a range and distribution of possible exposure start times per person. Do I nest the bulk of the code (L63-end) in a master for-loop? How would I index this and also, how would I save the output so I'm not overwriting previous monte carlo iteration results?

Thank you for any insights!!!!!! They are appreciated.

[darova]

This part of the code can be simpler

    if Water(t,column) == 1           
        Consumed(t,column) = Cons1; 
    else
        if Water(t,column) == 2           
        Consumed(t,column)= Cons2;
        else
            if Water(t,column) == 3
            Consumed(t,column)= Cons3;
            else
                if Water(t,column) == 4
                Consumed(t,column)= Cons4;
                else
                    if Water(t,column) == 5
                    Consumed(t,column)= Cons5;
                    else
                        if Water(t,column) == 6
                        Consumed(t,column)= Cons6;
                        else 
                            if Water(t,column) == 7
                            Consumed(t,column) = Cons7;
                            else
                                if Water(t,column) == 8
                                Consumed(t,column) = Cons8;
                                else
                                    Consumed(t,column)= Cons9;
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

Just use arrays

cons = [ 1 2 3 4 ... % some data ]
COnsumed(t,column) = cons(Water(t,column));

Answer 1

[Jan]

See darovas comment: Do not hide indices in names of variables but use arrays instead. See also: https://www.mathworks.com/matlabcentral/answers/57445-faq-how-can-i-create-variables-a1-a2-a10-in-a-loop

In your case, do not use Cons1, Cons2, ... but the vector Cons. Then your:

    if Water(t,column) == 1           
        Consumed(t,column) = Cons1; 
    else
        if Water(t,column) == 2           
        Consumed(t,column)= Cons2;
        else
           ...
    end
end

can be replaced by a single line:

Consumed(t, column) = Cons(Water(t, column));

The same applies for:

for column = (1:ncolumn)  
for t = 1:timesteps 
    
    if Zone(t,column) == 1           
    WatCon(t,column)= S1(t);        %lives in Security Zone 1
    else        

Even if there would be a reason for so many different IF conditions, it is much nicer to use ELSEIF instead of nested IF branches:

if Zone(t,column) == 2
WatCon(t,column) = S2(t);       %lives in Security Zone 2
else 
    if Zone(t,column) == 3
    WatCon(t,column) = S3(t);        %lives in Security Zone 3
    else
         if Zone(t,column) == 6
         WatCon(t,column) = W4(t);            %lives in Widefield Zone 4
         else
             ...
% Nicer:
if Zone(t,column) == 2
   WatCon(t,column) = S2(t);       %lives in Security Zone 2
elseif Zone(t,column) == 3
   WatCon(t,column) = S3(t);        %lives in Security Zone 3
elseif Zone(t,column) == 6
   WatCon(t,column) = W4(t);            %lives in Widefield Zone 4
else
   ...
       

Another option is a SWITCH block:

switch Zone(t,column)
    case 2
       WatCon(t,column) = S2(t);       %lives in Security Zone 2
    case 3
       WatCon(t,column) = S3(t);        %lives in Security Zone 3
    case 6
       WatCon(t,column) = W4(t);            %lives in Widefield Zone 4
    ...

But again an indexing is much nicer and more efficient.

Your code is easier to read if you care for a proper indentation. Press Ctrl-a Ctrl-i to fix this automatically.

Remember that statements like "it's not working" does not allow to understand, what your problem is. Post the error message instead or explain the difference between the result and your expectations.

If the nested IF branches are cleaned, the code will be much shorter. Then it is a good strategy to insert it into a function which replies the result as output. This makes it easier to call it in a loop and storing the results in an array.

Another hint:

You can replace the loop:

% I = Intake rate 
for column = (1:ncolumn)    
  for t= 1:timesteps
    I(t,column) = WatCon(t,column) * Consumed(t,column) / (W(t,column)/2.204623);   
  end 
end

by the faster and nicer:

I = WatCon .* Consumed ./ (W / 2.204623);

Matlab is such charming.

The final loop looks suspicious:

ExpoStart = zeros(1,ncolumn);   %preallocates for the start date of the last iteration of

for column = (1:ncolumn)
   for t = startpt:timesteps
      Serum(t,column) = Serum(t-1,column) + (I(t-1,column) * Ea / Vd(t-1,column)) - (Ke(t,column) * Serum(t-1,column));
   end
   
   
   while Serum(98,column)>Serum_blood(column)
      startpt=startpt+1;
      for t=startpt:timesteps
         Serum(t,column) = Serum(t-1,column) + (I(t-1,column) * Ea / Vd(t-1,column)) - (Ke(t,column) * Serum(t-1,column));
      end
   end
   
   ExpoStart(t,column) = startpt(t,column);
end

"ExpoStart = zeros(1,ncolumn)" is not a proper pre-allocation, if it is allocated with a single row only. Why not defining it with the final size?

You write to ExpoStart(t,column) only. t has the same value in all iterations. Is this really wanted?

[Wesser]

Thanks for the suggestions. I swapped the nested for-loops for swtch block per your recomendation and I'm getting the folloowing error:

SWITCH expression must be a scalar or a character vector.
Error in ToxMod_V7 (line 156)
switch Zone(t,column)

...so switch blocks aren't an optiion for matrices then? Is there a way to make a switch block work for a matrix? I like the idea of a Switch block as they are the easyest way for me to kep track of the relations being made. Also I don't quite understand how

    Consumed(t, column) = Cons(Water(t, column));

will correlate 1 to Cons1, 2 to Cons2, and so forth. Cons1, Cons2, etc are the output frrom the equations in L102-144. I'm clearly missing something here... Thank you for the thread about not creating unnessary variables...I understand the concept of that thread but still don't see how to correlate 1 to Cons1 (for example) in my code. Thank you for any feedback! *Updated code attached*

[Jan]

switch Zone(t,column)

...so switch blocks aren't an optiion for matrices then? Is there a way to make a switch block work for a matrix?

No, of course SWITCH can work with a scalar only. in the v7 version you have removed the loops over t and column. Instead of scalars as in v6 these two variables are vectors now. This cannot work. Insert the loops over these two variables.

If you do not create a list of variables Cons1, Cons2, ... but an array:

ConsList = {Cons1, Cons2, ...}

This code:

switch Water(t,column)
    case 1
       Consumed(t,column) = Cons1;               
    case 2
       Consumed(t,column) = Cons2;
    case 3
       Consumed(t,column) = Cons3;   
    case 4
       Consumed(t,column) = Cons4;
    case 5
       Consumed(t,column) = Cons5;
    case 6
       Consumed(t,column) = Cons6;   
    case 7
       Consumed(t,column) = Cons7;
    case 8
       Consumed(t,column) = Cons8;   
    case 9
       Consumed(t,column) = Cons9;     
end

can be simplified to:

Consumed(t, column) = ConsList{Water(t,column)};

The readability of your code will profit from a proper indentation: Ctrl-a Ctrl-i . This needs less than a second.

[Wesser]

Thank you for your feedback!