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

3 views (last 30 days)
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.
  1 Comment
darova
darova on 13 Mar 2021
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));

Sign in to comment.

Accepted Answer

Jan
Jan on 13 Mar 2021
Edited: Jan on 13 Mar 2021
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?
  3 Comments

Sign in to comment.

More Answers (0)

Products


Release

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!