Monte Carlo Simulation

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

• Given a random variable

X

and a function

h(X):

sample

X: [x 1 ,x 2 ,…,x n ];

Calculate

[h(x 1 ),h(x 2 ),…,h(x n )];

statistics of h.

use to approximate • Example:

X

is

U[0,1].

Use MCS to find mean of

X 2 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) 25 20 15 10 5 0 -3 -2 -1 0 1 2 3

Cumulative density function

• Cdfplot(y) 1 0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0 -3 • [f,x]=ecdf(y); -2 -1 Empirical CDF 0 x 1 2 3 4

Histogram of average

• 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.47

0.48

0.49

0.5

What is the law of large numbers?

0.51

0.52

0.53

Distribution of

x 2

x=rand(10000,1); x2=x.^2; hist(x2,20) 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, Weibull distribution 600 

for x

 0

b

 1

e

  

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

   

x

 

Y

 

y

 ] 

xy

  

x y

• For normal distribution can use Matlab’s mvnrnd • R = MVNRND(MU,SIGMA,N) SIGMA.

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

Example

mu = [2 3]; sigma = [1 1.5; 1.5 3]; plot(r(:,1),r(:,2),'+') 4.5

4 What is the correlation 3.5

coefficient?

3 2.5

2 1.5

1 0.5

0 0 0.5

1 1.5

2 2.5

3 3.5

Problems Monte Carlo

• Use Monte Carlo simulation to estimate the mean and standard deviation of

x 2

, 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

x 2

, 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

3 2 1 0 -1 -2 0 mu = [2 2]; sigma = [1 0; 0 3]; r = lhsnorm(mu,sigma,20); sum(r)/20 ans = 1.9732 2.0259

5 4 0.5

1 1.5

2 2.5

3 3.5

4 4 3 6 5 2 1 0 -1 -2 0.5

r = mvnrnd(mu,sigma,20); sum(r)/20 ans =2.3327 2.2184

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.

P f

 • 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/P f

 (1 

P f

) • For small

P f

  

P f

1

N P f m

problems probability of failure

1. Derive formula for the standard deviation of estimate of P f 2. 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