How to Code Monte carlo simulation?

2022 4 28
1 个回答
更新时间：2022 5 27
3 次查看（30 天）

0 个投票

In this code I am intending to get service life T(ser)
T(ser) = [ Pf(t) >= Pfmax ]
here , P is probability of an event.
Pfmax is generally between 7 to 10% and
Pf(t) = P[C >= Cr],
C = Cs[1-erf(Cd/(2*sqrt(D*T)))]
Here Cd, D, Cs, T, Cr are variables which are generated randomly, I have specified them in my code.
This is my code which is uncomplete. I am clueless on how to proceed further.
n = 20000;
%cover depth
mean_Cd = 16.1;
sigma_Cd = 0.23;
%diffusion coefficient at 28 days
mean_D = 3.87*10^(-12);
sigma_D = 0.52;
%surface chloride concn
mean_Cs = 13.1;
sigma_Cs = 0.103;
%Time exponent
mean_T = 0.2;
sigma_T = 0.2;
%critical chloride content
mean_Cr = 1.2;
sigma_Cr = 0.2;
%Now use normrnd function and lognrnd
Cd = normrnd(mean_Cd, sigma_Cd, [n,1]);
D = lognrnd(mean_D, sigma_D, [n,1]);
Cs = normrnd(mean_Cs, sigma_Cs, [n,1]);
T = normrnd(mean_T, sigma_T, [n,1]);
Cr = normrnd(mean_Cr, sigma_Cr, [n,1]);
C = zeros(n,1);
K = 0;
%Create for loop
for i=1:n
C(i,1) = Cs(i,1)[1-erf(Cd(i,1)/(2*sqrt(D(i,1)*T(i,1))))]
if
C(i,1)>=Cr(i,1)
K = K + 1;
end
end

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2022-4-30

0 个投票

C(i,1) = Cs(i,1)[1-erf(Cd(i,1)/(2*sqrt(D(i,1)*T(i,1))))]
You forgot the multiplication; MATLAB does not have implied multiplication anywhere
Also, in MATLAB, [] always has to do with building lists or arrays. Although it is often possible to use [] as a form of prioritization brackets, it is more efficient to use () instead of [] for prioritization
You should probably be considering vectorization. Generate a vector of Cs values and the corresponding C values, and then if you
Pf = mean(C > Cs)
No need to loop counting the occurances one by one.

8 个评论
显示 6更早的评论 隐藏 6更早的评论

harsh Brar
harsh Brar 2022-4-30
Thank You Mr. Roberson! I tried what you said but
When i applied this
for i=1:n
if (D(i,1).*T(i,1) > 0 ) % Check to see if this number is real and positive before calculating sqrt
con = D(i,1).*T(i,1);
else
con = 1;
%break;
end
C(i,1) = Cs(i,1).*[1-erf(Cd(i,1)./(2.*sqrt(con)))]
Pf = mean(C >= Cr);
end
this is my output
Arrays have incompatible sizes for this operation.
Error in mcs (line 45)
Pf = mean(C > Cs);
so, i made if loop that i was applying earlier.
K = 0;
%Create for loop
for i=1:n
if (D(i,1).*T(i,1) > 0 ) % Check to see if this number is real and positive before calculating sqrt
con = D(i,1).*T(i,1);
else
con = 1;
%break;
end
C(i,1) = Cs(i,1).*[1-erf(Cd(i,1)./(2.*sqrt(con)))]
if C(i,1)>=Cr(i,1)
K = K + 1;
end
end
P_fail = K./n;
display(P_fail)
Pfmax = [0.07 0.010];
Tser = [ P_fail >= Pfmax ]
display(Tser)
But my output is all zeroes. Can you please look into it?
C =
0
C =
0
0
C =
0
0
0
C =
1.0e-09 *
0
0
0
0.7618
C =
1.0e-09 *
0
0
0
0.7618
0.0000
P_fail =
0
Tser =
1×2 logical array
0 0
Tser =
1×2 logical array
0 0
Walter Roberson
Walter Roberson 2022-4-30
I said, "Generate a vector of Cs values and the corresponding C values," -- as in generate the entire Cs and C vector before calling mean() in the way I showed.
harsh Brar
harsh Brar 2022-4-30
Thank you for replying again!
I have made the vectors as u have said. And There's one thing. I intend to calculate
Pf = P[C >= Cr], %not Cs
These are the changes
for i=1:n
if (D(i,1).*T(i,1) > 0 ) % Check to see if this number is real and positive before calculating sqrt
con = D(i,1).*T(i,1);
else
con = 1;
%break;
end
rng(0,'twister'); %generating a vector of 1000 random values for Cr values
a = 0.2;
b = 1.2;
Cr = a.*randn(100,1) + b;
display(Cr)
C(i,1) = Cs(i,1).*(1-erf(Cd(i,1)./(2.*sqrt(con))));
a = i;
b = 1;
C = a.*randn(100,1) + b;
display(C)
Pf = mean(C > Cr);
display(Pf)
end
Now, I am facing the difficulty in calculating the final output(Tser), which will be calculated on the basis of
T(ser) = [ Pf >= Pfmax ] and here Pfmax is in the range of 7 to 10%( or we can take it only 10%). I have tried making something but it's not working.
for
Pfmax = 0.07:0.010 ;
Tser = mean(Pf >= Pfmax);
display(Tser)
end
Walter Roberson
Walter Roberson 2022-4-30
here , P is probability of an event.
Pf(t) = P[C >= Cr],
Okay, P[something] means probability
T(ser) = [ Pf(t) >= Pfmax ]
That uses the [] part of the notation, but not the P part of the notation. Is it supposed to be a probability? Does it make sense to ask about the probability of a probability ? Or is it just intended to be a logical vector same length as Pfmax ?
You are asked to compute T(ser) but the right hand side of that equation does not mention ser so it is not clear what is being calculated.
for
Pfmax = 0.07:0.010 ;
Tser = mean(Pf >= Pfmax);
display(Tser)
end
Pf is a scalar at that point, and you are looping over Pfmax values so each of them will be a scalar. mean() of a scalar is the same as the scalar (except possibly converted from logical to numeric.)
Also, 0.07:0.010 is the same as 0.07 : 1 : 0.010 so the next potential term after the starting point of 0.07 would be 1.07 which is greater than the end point of 0.010, so there would be at most one entry in the vector. Are there any entries in the vector? No, there are not: 0.010 is less than 0.07 so the loop is not going to execute at all.
harsh Brar
harsh Brar 2022-5-1
i got the answer. thank you for your efforts Mr. roberson.
harsh Brar
harsh Brar 2022-5-5
hello Mr. walter, I would like to ask for help again. My code is not giving me the desired result. It always shows 1. For your ease let me explain the base of code again
The failure occur when "C >= Cr". so I intend to calculate the probabilty of that failure, lets say its P_fail. Here, C will be Calculated using
C(i,1) = Cs(i,1).*(1-erf(Cd(i,1)./(2.*sqrt(con))));
and Cr = 1.2
and ultimately final output is calculated when P_fail becomes greater than P_fail_Max (which equals to 0.10).
clear all
clc
n = 200;
%cover depth
mean_Cd = 16.1;
sigma_Cd = 0.23;
%diffusion coefficient at 28 days
mean_D = 3.87*10^(-12);
sigma_D = 0.52;
%surface chloride concn
mean_Cs = 13.1;
sigma_Cs = 0.103;
%Time exponent
mean_T = 0.2;
sigma_T = 0.2;
%critical chloride content
mean_Cr = 1.2;
sigma_Cr = 0.2;
%Now use normrnd function and lognrnd to generate values randomly
Cd = normrnd(mean_Cd, sigma_Cd, [n,1]);
D = lognrnd(mean_D, sigma_D, [n,1]);
Cs = normrnd(mean_Cs, sigma_Cs, [n,1]);
T = normrnd(mean_T, sigma_T, [n,1]);
Cr = normrnd(mean_Cr, sigma_Cr, [n,1]);
%Create for loop
for i=1:n
rng(0,'twister'); %generating a vector of 200 random values for Cr values
a = 0.52;
b = 3.87*10^(-12);
D = a.*randn(200,1) + b;
%display(D)
rng(0,'twister'); %generating a vector of 200 random values for Cr values
a = 0.2;
b = 0.2;
T = a.*randn(200,1) + b;
% display(T)
if (D(i,1).*T(i,1) > 0 ) % Check to see if this number is real and positive before calculating sqrt
con = D(i,1).*T(i,1);
else
con = 1;
%break;
end
rng(0,'twister'); %generating a vector of 200 random values for Cr values
a = 0.2;
b = 1.2;
Cr = a.*randn(200,1) + b;
%display(Cr)
rng(0,'twister'); %generating a vector of 200 random values for Cr values
a = 0.23;
b = 16.1;
Cd = a.*randn(200,1) + b;
% display(Cd)
rng(0,'twister'); %generating a vector of 1000 random values for Cr values
a = 0.103;
b = 13.1;
Cs = a.*randn(200,1) + b;
%display(Cs)
% MAIN EQUATION TO CALCULATE C
C(i,1) = Cs(i,1).*(1-erf(Cd(i,1)./(2.*sqrt(con))));
a = 1;
b = i;
C = a.*randn(20,1) + b;
display(C)
Pf = mean(C > Cr);
display(Pf)
end
Pfmax = 0.10;
Tser = mean(Pf >= Pfmax);
display(Tser) % final output
the output should be greater than 20. I hope you will look into it.
abdellah Chougrani
abdellah Chougrani 2022-5-27
Hello,
Have you checked the units ?
Walter Roberson
Walter Roberson 2022-5-27
编辑：Walter Roberson 2022-5-27
Cr = a.*randn(200,1) + b;
A column vector of length 200
C(i,1) = Cs(i,1).*(1-erf(Cd(i,1)./(2.*sqrt(con))));
C has not been initialized. In that context you would be creating a column vector.
C = a.*randn(20,1) + b;
That overwrites all of C with a column vector of length 20
Pf = mean(C > Cr);
column vectors of different lengths. Error.
If your actual code makes C the right length then the mean would be a scalar
Tser = mean(Pf >= Pfmax);
After the for loop, Pf is still a scalar. mean of a logical scalar is 0 or 1
mean of logical is a probability, never greater than 1

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Elementary Math 的更多信息

产品

版本

R2018a

标签

Community Treasure Hunt

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

Start Hunting!

Translated by