MATLAB Answers

How to plot a pdf and cdf for my code

45 views (last 30 days)
mike genzen
mike genzen on 24 Apr 2018
Commented: mike genzen on 24 Apr 2018
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')

  0 Comments

Sign in to comment.

Answers (1)

njj1
njj1 on 24 Apr 2018
Edited: njj1 on 24 Apr 2018
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.

  8 Comments

Show 5 older comments
mike genzen
mike genzen on 24 Apr 2018

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

njj1
njj1 on 24 Apr 2018

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
mike genzen on 24 Apr 2018
Thank you for the help on this. I definitely learned afew things. Again thank you!

Sign in to comment.

Sign in to answer this question.