180 likes | 426 Views
Latin Hypercube Sampling Example. Jake Blanchard Spring 2010. Example. Z=XY X and Y follow exponential distributions x =1 y =1/2. Step 1. Divide cdfs into even intervals (vertical axis). Step 2. Now sample a number in each section
E N D
Latin Hypercube Sampling Example Jake Blanchard Spring 2010 Uncertainty Analysis for Engineers
Example • Z=XY • X and Y follow exponential distributions • x=1 • y=1/2 Uncertainty Analysis for Engineers
Step 1 • Divide cdfs into even intervals (vertical axis) Uncertainty Analysis for Engineers
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
Step 3 • Sort the values (shuffle) • nux=xx(randperm(n)); Uncertainty Analysis for Engineers
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
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
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
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
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
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
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
Results Uncertainty Analysis for Engineers