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

Monte Carlo Simulation SOURCE: http://pics.hoobly.com/full/AA7G6VQPPN2A.jpg 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

Basic Monte Carlo SOURCE: http://www.sz-wholesale.com/uploadFiles/041022104413s.jpg • Given 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 x=rand(10); y=x.^2; %generates 10x10 random matrix mean=sum(y)/10 x =0.4017 0.5279 0.1367 0.3501 0.3072 0.3362 0.3855 0.3646 0.5033 0.2666 mean=0.3580 • What is the true mean SOURCE: http://schools.sd68.bc.ca/ed611/akerley/question.jpg

Obtaining distributions • Histogram: y=randn(100,1); hist(y)

Cumulative density function • Cdfplot(y) • [f,x]=ecdf(y);

Histogram of average • x=rand(100); y=sum(x)/100; hist(y)

Histogram of average • x=rand(1000); y=sum(x)/1000; hist(y) What is the law of large numbers?

Distribution of x2 x=rand(10000,1); x2=x.^2; hist(x2,20)

Other distributions • Other distributions available in matlab • For example, Weibulldistribution r=wblrnd(1,1,1000); hist(r,20)

Correlated Variables • For normal distribution can use Matlab’smvnrnd • 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.

Example mu = [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 Carlo • Use 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 sampling X = 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 LHS mu = [2 2]; sigma = [1 0; 0 3]; r = lhsnorm(mu,sigma,20); sum(r)/20 ans = 1.9732 2.0259 r = mvnrnd(mu,sigma,20); sum(r)/20 ans =2.3327 2.2184

Evaluating probabilities of failure • Failure 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, N • The accuracy of the estimate is poor unless N is much larger than 1/Pf • For small Pf

problems probability of failure • Derive formula for the standard deviation of estimate of Pf • If 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 Institution Number: 2004-57325

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

Reading assignment Ravishankar, 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

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