How to plot a pdf and cdf for my code

I am just scratching the surface with monte carlo and distributions and am looking for a solution to plotting a pdf and cdf for my code, aswell as a brief explanation of setting it up. My attempts used norm=normpdf(Y,averageY,sigmaY) with x=Y then figure;plot(x,norm). This was clearly inccorect as the pdf should peak around .07

clc clear

% A simple program to develop an understanding of random number generation, transformation, % Monte Carlo simulation, and analysis of the outputs.

% This code can also be used to find the overlap (probability of failure) between any two normal distributions* % *use the analytical result and/or a sufficient number of trials

clear all;

trials = 100000; fprintf('Number of trials is set to %i.\n',trials) x1 = randn(trials,1); mu1 = 500; sigma1 = 100; fprintf('The mean and standard deviation of input 1 is %f, %f, respectively\n', mu1, sigma1) x2 = randn(trials,1); mu2 = 1000; sigma2 = 100; fprintf('The mean and standard deviation of input 2 is %f, %f, respectively\n', mu2, sigma2)

y1 = []; y2 = []; Y = [];

sum = 0; for i = 1:trials y1(i) = x1(i)*sigma1 + mu1; % transform the normally distributed randon numbers for input 1 y2(i) = x2(i)*sigma2 + mu2; % transform the normally distributed randon numbers for input 2 Y(i) = (3-((4*1000000)/(30000000*2*4))*(sqrt((y2(i)./16)^2+(y1(i)./4)^2))); % The output, where "failure" is defined as y2 is less than y1 if Y(i) < 0 sum = sum + 1; end end

%y1

% Monte Carlo predicted probability fprintf('Monte Carlo predicted probability of failure for %i trials is:\n',trials); probF = sum/trials

%--------------------------------% % find the analytical "z" value analytical = ((-(mu1 - mu2))/sqrt(sigma1^2 + sigma2^2)); %equation found in Shigley, page 25 in tenth edition

% integrate the gaussian cdf to find the probability % the function, in terms of "x", to integrate: fun = @(x)(1/(sqrt(2*pi)))*exp(-((x.^2)/2));

% perform the integration from -infinity to z (used -100 for practical reasons) % this is th equivalent of looking up the z value in a table (e.g. A-10 in Shigley) % change this to the "integral" function if using Matlab area=quad(fun,-100,analytical);

% Analytical result - only valid if both distributions are normal and the integration works!!! fprintf('Analytical probability of failure is %f. Note: only valid if both input and output are normal!\n',area) %--------------------------------%

%--------------------------------% % Calculate statistics on the output, Y, and plot the results

averageY = mean(Y); sigmaY = std(Y);

fprintf('The mean and standard deviation of the output are %f, %f, respectively\n',averageY, sigmaY)

 figure
  % a line for y1 - y2 = 0, i.e. the "failure" line - Change this to define "failure" for the problem.
  plot(y1, y2, 'k.')
  hold on
  xlabel('X1')
  ylabel('X2')

Answer 1

[njj1]

If you want to compute the pdf for a normal distribution over a range of values, you must input a vector of possible values sorted from low to high, the mean of the distribution, and the standard deviation of the distribution. For example:

averageY = mean(Y);
sigmaY = std(Y);
x = linspace(min(Y),max(Y),10000); %this produces a sorted vector over which to plot the pdf
norm = normpdf(x,averageY,sigmaY);
figure;
plot(x,norm,'-k')

Hope this helps.

[njj1]

The cdf for this can be computed using the pdf.

norm_cdf = cumsum(norm)./sum(norm); %make sure your variable "sum" is cleared before doing this operation...
figure;
plot(x,norm_cdf,'-k')

[mike genzen]

My only question is if I compute the area under the curve should it not equal the probability of failure? The picture I attached was what the outcome should be.

[njj1]

If you run the code that I typed above, you will get the same results as shown in your images (with the only difference being in the resolution along the x-axis). Except the image has the pdf and cdf reversed.

Yes, if you compute area under the pdf over a specified interval, then this is equal to the probability of the event occurring in that interval (i.e., P(f|t1<t<t2) = int_{t1}^{t2} pdf(t)dt). The total area under the pdf, however, will always equal 1.

The idea that this is the probability of failure does not make sense when you use a normal distribution because a normal distribution always has a greater than zero probability at negative values. This implies that a component or system could fail in negative time, which is nonsensical. Most of the time in reliability engineering/analysis, the Weibull distribution is used (also the log-normal is used occasionally).

[mike genzen]

This is the plot that populates when I run the code. Also what do you mean by sum is cleared for the cdf?

[njj1]

The plot that you showed in Capture2.JPG looks exactly like the Figure 9.6 shown in Capture1.JPG, except that the resolution along the x-axis is higher.

You typed this code in your OP:

sum = 0; for i = 1:trials y1(i) = x1(i)*sigma1 + mu1; % transform the normally distributed randon numbers for input 1 y2(i) = x2(i)*sigma2 + mu2; % transform the normally distributed randon numbers for input 2 Y(i) = (3-((4*1000000)/(30000000*2*4))*(sqrt((y2(i)./16)^2+(y1(i)./4)^2))); % The output, where "failure" is defined as y2 is less than y1 if Y(i) < 0 sum = sum + 1; end end

In this code, you have a variable called "sum". I suggested that you use the following code to form your cdf:

norm_cdf = cumsum(norm)./sum(norm); %make sure your variable "sum" is cleared before doing this operation...

You must clear your variable "sum" before you use the Matlab command sum().

[mike genzen]

What about the Y axis? fig 9.6 shows the y axis peak close to .07 vs 1.0 from the code

[njj1]

Honestly, without seeing any of the data that went into producing that plot, I can't say for sure what the difference is. The pdf should always integrate to 1, and when I integrate the pdf I produce, the value is 1.

sum(norm.*mean(diff(x)))
ans =
  0.9999

Perhaps they used a different method of normalization.

[mike genzen]

Thank you for the help on this. I definitely learned afew things. Again thank you!