28 views (last 30 days)

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')

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.

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.

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

Start Hunting!
## 0 Comments

Sign in to comment.