% MCS_intro.m
% Monte Carlo Simulation ... an introductory example
%
% Y = g(X1,X2,X3) = sin(X1) + sqrt(X2) - exp(-X3) - 2;
%
% H.P. Gavin, Dept. Civil and Environmental Engineering, Duke Univ, Jan. 2012

ppl = [ 68, 1, 84 ]/256;    % purple to match histogram bars

%	X1		X2		X3	
%	normal		lognormal	Rayleigh
	mu1 = 6;	med2 = 2;	mu3 = 1;
	sd1 = 2;	cov2 = 0.3;	

N = 1000;		% number of random values in the sample

% (1) generate a large sample for each random variable in the problem ...

X1 = mu1 + sd1*randn(1,N);
X2 = logn_rnd(med2,cov2,1,N);
X3 = Rayleigh_rnd(mu3,1,N);

% (2) evaluate the function for each random variable to compute a new sample

Y = sin(X1) + sqrt(X2) - exp(-X3) - 2;

% suppose "probability of failure" is Prob[ g(X1,X2,X3) > 0 ] ...

Probability_of_failure = sum(Y>0) / N

% (3) plot histograms of the random variables 

sort_X1 = sort(X1);
sort_X2 = sort(X2);
sort_X3 = sort(X3);
CDF     = [1:N] / (1+N);		% empirical CDF of all quantities
% ---------------------------------------------- plots ... 
% 1-83, 42-43, 47, 52-53, 56, 60-61, 64, 69, 73-74, 78

nBins = floor(N/50);
figure(1)         
 clf
 subplot(231)
  [fx,xx] = hist( X1, nBins, nBins/(max(X1)-min(X1)) );	% histogram of X1
  hold on
   bar(xx,fx, 'FaceColor', ppl)
   plot(sort_X1,norm_pdf(sort_X1,mu1,sd1), 'LineWidth',2);
  hold off
   axis('tight')
   ylabel('P.D.F.')
 subplot(234)
  hold on
    stairs( sort_X1,CDF, 'color', ppl, 'LineWidth',2)
    plot(sort_X1,norm_cdf(sort_X1,mu1,sd1), 'LineWidth',0.5 )
  hold off
   axis('tight')
   ylabel('C.D.F.')
   xlabel('X_1 : normal')
 subplot(232)
  [fx,xx] = hist( X2, nBins, nBins/(max(X2)-min(X2)) );	% histogram of X2
  hold on
    bar(xx,fx, 'FaceColor', ppl)
    plot(sort_X2,logn_pdf(sort_X2,med2,cov2), 'LineWidth',2);
  hold off
   axis('tight')
 subplot(235)
  hold on
    stairs(sort_X2,CDF, 'color', ppl, 'LineWidth',2)
    plot(sort_X2,logn_cdf(sort_X2,[med2,cov2]), 'LineWidth',0.5)
  hold off
   axis('tight')
   xlabel('X_2 : log-normal')
 subplot(233)
  [fx,xx] = hist( X3, nBins, nBins/(max(X3)-min(X3)) );	% histogram of X2
  hold on
    bar(xx,fx, 'FaceColor', ppl)
    plot(sort_X3,Rayleigh_pdf(sort_X3,mu3), 'LineWidth',2 );
  hold off
   axis('tight')
 subplot(236)
  hold on
    stairs(sort_X3,CDF, 'color', ppl, 'LineWidth',2)
    plot(sort_X3,Rayleigh_cdf(sort_X3,mu3), 'LineWidth',0.5 )
  hold off
   axis('tight')
   xlabel('X_3 : Rayleigh')

  print('MCS_intro1.pdf', '-dpdfcrop')

% (4) plot histogram of the random function Y = g(X1,X2,X3)

nBins = floor(N/20);
figure(2)
 clf
 subplot(211)
  [fx,xx] = hist( Y, nBins, nBins/(max(Y)-min(Y)) );	% histogram of Y
  hold on
    bar(xx,fx, 'FaceColor', ppl)
    plot([0 0],[0 0.5],'-r','LineWidth',3)
   hold off
   axis('tight')
   ylabel('P.D.F.')
 subplot(212)
  hold on
    stairs(sort(Y),CDF, 'color', ppl)
    plot( [0 0],[0 1],'-r','LineWidth',3);
  hold off
   axis('tight')
   text(0.5,0.5,'Y>0')
   ylabel('C.D.F.')
   xlabel('Y = g(X_1,X_2,X_3)');

  print('MCS_intro2.pdf', '-dpdfcrop')


