Latin Hypercube Sampling Example Jake Blanchard Spring 2010
Latin Hypercube Sampling Example Jake Blanchard Spring 2010 Uncertainty Analysis for Engineers 1
Example �Z=XY �X and Y follow exponential distributions � x=1 � y=1/2 Uncertainty Analysis for Engineers 2
Step 1 �Divide cdfs into even intervals (vertical axis) Uncertainty Analysis for Engineers 3
Step 2 �Now sample a number in each section �For example, pick a random number between 0. 8 and 1. 0 and use it to get a random value for x �xx=expinv(rx, mux); Uncertainty Analysis for Engineers 4
Step 3 �Sort the values (shuffle) �nux=xx(randperm(n)); Uncertainty Analysis for Engineers 5
An Alternative �Instead of using expinv, we can generate the inverse ourselves �Just take the CDF and solve for x �xx=-mux*log(1 -rx); Uncertainty Analysis for Engineers 6
First Script n=10000000; mux=1; muy=2; x=exprnd(mux, n, 1); y=exprnd(muy, n, 1); mz=mean(x. *y); error=abs(mz-mux*muy)/mux/muy d=linspace(0, 1, n+1); rx=unifrnd(d, d+1/n); ry=unifrnd(d, d+1/n); xx=expinv(rx, mux); yy=expinv(ry, muy); nux=xx(randperm(n)); nuy=yy(randperm(n)); mz=mean(nux. *nuy); error=abs(mz-mux*muy)/mux/muy Uncertainty Analysis for Engineers 7
Alternative n=10000000; mux=1; muy=2; x=exprnd(mux, n, 1); y=exprnd(muy, n, 1); mz=mean(x. *y); error=abs(mz-mux*muy)/mux/muy d=linspace(0, 1, n+1); rx=unifrnd(d, d+1/n); ry=unifrnd(d, d+1/n); xx=-mux*log(1 -rx); yy=-muy*log(1 -ry); nux=xx(randperm(n)); nuy=yy(randperm(n)); mz=mean(nux. *nuy); error=abs(mz-mux*muy)/mux/muy Uncertainty Analysis for Engineers 8
Test �Z=XY �X and Y are beta with mean of 1 and 2, respectively �Use simple Monte Carlo �Then use LHS without sorting �Then use LHS with sorting �N=100, 000 �For each case, find mean 100 times and then take standard deviation of results Uncertainty Analysis for Engineers 9
Case 1 n=100000; ntrials=100; mz=zeros(ntrials, 1); for i=1: ntrials x=exprnd(mux, n, 1); y=exprnd(muy, n, 1); mz(i)=mean(x. *y); end std(mz) Uncertainty Analysis for Engineers 10
Case 2 d=linspace(0, 1, n+1); for i=1: ntrials rx=unifrnd(d, d+1/n); rx=rx(1: end-1); ry=unifrnd(d, d+1/n); ry=ry(1: end-1); x=expinv(rx, mux); y=expinv(ry, muy); mz(i)=mean(x. *y); end std(mz) Uncertainty Analysis for Engineers 11
Case 3 d=linspace(0, 1, n+1); for i=1: ntrials rx=unifrnd(d, d+1/n); rx=rx(1: end-1); ry=unifrnd(d, d+1/n); ry=ry(1: end-1); x=expinv(rx, mux); y=expinv(ry, muy); nux=x(randperm(n)); nuy=y(randperm(n)); mz(i)=mean(nux. *nuy); end std(mz) Uncertainty Analysis for Engineers 12
Results Mean(mz) Std(mz) Case 1 1. 9999 0. 0098 Case 2 4. 0000 0. 00029 Case 3 1. 9998 0. 0064 Uncertainty Analysis for Engineers 13
- Slides: 13