How to Code Monte carlo simulation?

Question

harsh Brar 2022-4-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1707460-how-to-code-monte-carlo-simulation

编辑： Walter Roberson 2022-5-27

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Walter Roberson 2022-4-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1707460-how-to-code-monte-carlo-simulation#answer_954315

在 MATLAB Online 中打开

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 2022-4-30

在 MATLAB Online 中打开

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 2022-4-30

在 MATLAB Online 中打开

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 2022-5-5

在 MATLAB Online 中打开

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 2022-5-27

Hello,

Have you checked the units ?

Walter Roberson 2022-5-27

编辑：Walter Roberson 2022-5-27

在 MATLAB Online 中打开

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

请先登录，再进行评论。

How to Code Monte carlo simulation?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

How to Code Monte carlo simulation?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

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