How do we generate the statistics of a function of a random variable? - PowerPoint Presentation

Slide1

How do we generate the statistics of a function of a random variable?Why is the method called “Monte Carlo?” How do we use the uniform random number generator to generate other distributions?Are other distributions directly available in matlab? How do we accelerate the brute force approach?Probability distributions and moments Web links: http://www.riskglossary.com/link/monte_carlo_method.htm http://physics.gac.edu/~huber/envision/instruct/montecar.htm

Monte Carlo Simulation

SOURCE: http://pics.hoobly.com/full/AA7G6VQPPN2A.jpg

Basic Monte CarloGiven a random variable X and a function h(X): sample X: [x1,x2,…,xn]; Calculate [h(x1),h(x2),…,h(xn)]; use to approximate statistics of h.Example: X is U[0,1]. Use MCS to find mean of X2

rng(9); %set seed so that script will produce repeatable results x=rand(10);y=x.^2; meancol=mean(y)meancol =0.0947 0.4567 0.5221 0.4693 0.3655 0.1874 0.2758 0.2041 0.2619 0.2485meanall=mean(meancol)meanall =0.3086What is the true mean

SOURCE: http://schools.sd68.bc.ca/ed611/akerley/question.jpg

SOURCE:

http://www.sz-wholesale.com/uploadFiles/041022104413s.jpg

Visualizing distributions x=randn(10000,1); hist(x,20) vs. x=randn(100,1); hist(x,10)

Cumulative density function (CDF)CDF shows the probability F of X<xCdfplot(x)With 100samples[cdf,xaxis]=ecdf(x); See notes page

With 10,000 samples it looks betterUsing 10,000 samples we get

Histogram of averagex=rand(1000);y=mean(x);hist(y)We start with uniform distribution, but averaging moves it towards normal narrowing toward the meanstd(y)=0.0092std(x)=0.2887

Distribution of x2x=rand(1000);z=reshape(x,[1000000,1]);z2=z.^2; cdfplot(z2)hist(z2,500)As expected the CDF(0.81)=0.9

Other distributionsOther distributions available in matlabFor example, Weibull distribution

y=wblrnd(100,1.5,1,10000000);>> hist(y,500)

Correlated VariablesFor normal distribution can use Matlab’s mvnrnd R = MVNRND(MU,SIGMA,N) returns a N-by-D matrix R of random vectors chosen from the multivariate normal distribution with 1-by-D mean vector MU, and D-by-D covariance matrix SIGMA.

Examplemu = [2 3];sigma = [1 1.5; 1.5 3];r = mvnrnd(mu,sigma,20);plot(r(:,1),r(:,2),'+')What is the correlation coefficient?

Problems Monte CarloUse Monte Carlo simulation to estimate the mean and standard deviation of x2, when X follows a Weibull distribution with a=b=1.Calculate by Monte Carlo simulation and check by integration the correlation coefficient between x and x2, when x is uniformly distributed in [0,1]

Latin hypercube samplingX = lhsnorm(mu,SIGMA,n) generates a latin hypercube sample X of size n from the multivariate normal distribution with mean vector mu and covariance matrix SIGMA. X is similar to a random sample from the multivariate normal distribution, but the marginal distribution of each column is adjusted so that its sample marginal distribution is close to its theoretical normal distribution.

Comparing MCS to LHSmu = [2 2];sigma = [1 0.2; 0.2 3];r = lhsnorm(mu,sigma,20);mean(r)= 2.0389 2.0316cov(r) 1.0173 0.8496 0.8496 3.3286

rng(9); r = mvnrnd(mu,sigma,20);plot(r(:,1),r(:,2),'+')mean(r)=1.8526 2.9458cov(r) 1.1215 1.1709 1.1709 3.5390

Evaluating probabilities of failureFailure is defined in terms of a limit state function that must satisfy g(r)>0, where r is a vector of random variables.Probability of failure is estimated as the ratio of number of negative g’s, m, to total MC sample size, NThe accuracy of the estimate is poor unless N is much larger than 1/PfFor small Pf

problems probability of failureDerive formula for the standard deviation of estimate of PfIf x is uniformly distributed in [0,1], use MCS to estimate the probability that x2>0.95 and estimate the accuracy of your estimate from the formula.3. Calculate the exact value of the answer to Problem 2 (that is without MCS).

Source: Smithsonian InstitutionNumber: 2004-57325

Separable Monte CarloUsually limit state function is written in terms of response vs. capacity g=C(r)-R(r)>0Failure typically corresponds to structures with extremely low capacity or extremely high response but not bothCan take advantage of that in separable MC

Reading assignmentRavishankar, Bharani, Smarslok B.P., Haftka R.T., Sankar B.V. (2010)“Error Estimation and Error Reduction in Separable Monte Carlo Method ” AIAA Journal ,Vol 48(11), 2225–2230 .

Source: www.library.veryhelpful.co.uk/ Page11.htm