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 A other th distributions di t ib ti directly di tl available il bl in i matlab? tl b? How do we accelerate the brute force approach? – Probability distributions and moments
Web links: http://www.riskglossary.com/link/monte_carlo_method.htm htt // i k l /li k/ t l th d ht http://physics.gac.edu/~huber/envision/instruct/montecar.htm
SOURCE: http://p pics.hoobly.com/full/AA A7G6VQPPN2A.jpg
Monte Carlo Simulation
• 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); ( ) yy=x.^2; %generates g 10x10 random matrix sumy=sum(y)/10 sumy =0.4017 0.5279 0.1367 0.3501 0.3072 0.3362 0.3855 0.3646 0.5033 0.2666 sum(sumy)/10 ans =0 =0.3580 3580
• What is the true mean
SOURCE: http://schools.sd68.bc.ca/ed611/akerley/question.jpg
SOURCE E: http://www w.sz-wholesale.com/uploadFiiles/041022104413s.jpg
Basic Monte Carlo
Basic Monte Carlo • Given a random variable X and a function h(X): sample X: [[x1,,x2,,…,x , n]; Calculate [[h(x ( 1), ),h(x ( 2), ),…,h(x , ( n)]; 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 sumy=sum(y)/10 sumy =0 =0.4017 4017 0 0.5279 5279 0 0.1367 1367 0 0.3501 3501 0 0.3072 3072 0.3362 0.3855 0.3646 0.5033 0.2666 sum(sumy)/10 ( y) ans =0.3580
• What is the true mean
SOURCE: http://schools.sd68.bc.ca/ed6 11/akerley/question.jpg
Obtaining distributions • Histogram: y=randn(100 y=randn(100,1); 1); hist(y) 25
20
15
10
5
0 -3
-2
-1
0
1
2
3
Cumulative density function • Cdfplot(y) p (y) Empirical CDF
1 0.9 08 0.8 0.7
F(x )
0.6 0.5 0.4 0.3 0.2 0.1 0 -3 3
-2 2
-1 1
0
1 x
2
3
4
Experimental CDF • [f,x]=ecdf(y); [f x]=ecdf(y); f’f Columns 1 through 8 0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 Columns 9 through 16 0.0800 0.0900 0.1000 0.1100 0.1200 0.1300 0.1400 0.1500
• x‘ Columns 1 through 8 -2.2529 -2.2529 -2.1763 -2.1746 -1.9539 -1.9068 -1.6354 -1.5729 Columns 9 through 16 -1.4127 1 4127 -1.2568 1 2568 -1.2203 1 2203 -1.2158 1 2158 -1.1957 1 1957 -1.1684 1 1684 -1.1645 1 1645 -1.0954 1 0954 Columns 89 through 96 0.9940 1.0437 1.0945 Columns 97 through 101 1.5908 1.7909 1.8918
1.1654
1.2020
2.4124
2.7335
1.2750
1.5081
1.5338
Histogram g of average g • x=rand(100); y=sum(x)/100; hist(y) 35
30
25
20
15
10
5
0 0.42
0.44
0.46
0.48
0.5
0.52
0.54
0.56
0.58
Histogram of average • x=rand(1000); y=sum(x)/1000; hist(y) 140
120
100
80
60
40
20
0 0.46 0 46
0 47 0.47
0 48 0.48
0 49 0.49
05 0.5
0 51 0.51
What is the law of large numbers?
0 52 0.52
0 53 0.53
Distribution of x=rand(10000,1); ( , ); x2=x.^2; hist(x2,20)
2 x
2500
2000
1500
1000
500
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Other distributions • Other distributions available in matlab • For example, Weibul distribution 600
b
b 1 x / a
f ( x) ba x e for x 0
b
500
400
r=wblrnd(1,1,1000); hist(r,20)
300
200
100
0
0
2
4
6
8
10
12
14
Correlated Variables • For normal distribution distrib tion can use se Matlab’s m mvnrnd nrnd • R = MVNRND(MU,SIGMA) returns an n-by-d matrix R of random vectors chosen from the multivariate normal distribution with mean vector MU, and covariance matrix SIGMA. MU is an n-by-d matrix, and MVNRND generates each row of R using the corresponding row of MU. SIGMA is a d-by-d symmetric positive semidefinite matrix
Example p mu = [2 3]; sigma = [1 1 1.5; 5; 1 1.5 5 3]; r = mvnrnd(mu,sigma,20); 5 4.5 plot(r(:,1),r(:,2),'+')) plot(r(:,1),r(:,2), 4 3.5 3 2.5 2 1.5 1 0.5 0
0
0.5
1
1.5
2
2.5
3
3.5
Latin hypercube yp sampling p g 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 lti i t normall di distribution, t ib ti b butt th the marginal distribution of each column is adjusted so that its sample marginal distribution is close to its theoretical normal distribution 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 3327 2 2.2184 2184
5
6
4
5 4
3
3
2 2
1 1
0
0
-1
-2
-1
0
0.5
1
1.5
2
2.5
3
3.5
4
-2 0.5
1
1.5
2
2.5
3
3.5
4
4.5
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. Pf m / N • 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 Pf (1 Pf ) i much is h llarger than h 1/Pf Pf N • For small Pf Pf Pf
1 m
Separable Monte Carlo • Usually limit state function is written in terms of response vs. capacity g C(r) R(r)>0 g=C(r)-R(r)>0 • Failure typically corresponds to structures with extremely low capacity or extremely high response but not both • Can C take t k advantage d t off that th t in i separable bl MC
Reading g assignment g SDM paper by Ben Smarslok
Source: www.library.veryhelpful.co.uk/ Page11.htm
problems 1. Derive formula for the standard deviation of estimate of Pf 2. Use MCS to calculate the probability of failure when a structure is designed with 50% safety factor and both stress and failure stress have a normal distribution with 10% COV.
Source: Smithsonian Institution Number: 2004-57325
3. Calculate the exact value of the answer to Problem 2 (that is without MCS).