How to Code Monte carlo simulation?

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

Answer 1

[Walter Roberson]

  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.

[harsh Brar]

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]

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]

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]

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]

i got the answer. thank you for your efforts Mr. roberson.

[harsh Brar]

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]

Hello,

Have you checked the units ?

[Walter Roberson]

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