LHS.zip_ARGUMENT!_LHS matlab_LHS采样_dist.m_拉丁 抽样

function S=LHS(m,dist,mu,sigma,lowb,upb) %code of Latin Hypercube Sampling% %call S=lhs(m,dist,mu,sigma,lowb,upb) % %Input argument %m: a scalar,the number of sample points %dist: A row with distribution type flags of basic random variables;the %value of the flag can be 1 (for uniform distribution, 2(for normal distribution), 3(for lognormal) %and 4(for extreme type 1). %mu: A row vector comprising the mean value of basic random variables. %sigma: A row vector with its length equaligng to mu,including the standard %deviation of basic random variables. %lowb: a row vector with its elements are the lower bound of the sampling %interval %upb:a row vector with its elements are the upper bounds of the sampling %interval %dist,mu,sigma,lowb,upb must have the same length. %Output argument %S: sampling point matrix, of which each row is a sampling point. % %last update june 13, 2009 n=length(mu); if length(dist)~=n|length(sigma)~=n|length(lowb)~=n|length(upb)~=n error('dist,mu,sigma,lowb,upb must have the same length'); end rvcom=[]; for j=1:n rv=[]; if dist(j)==1 p_low(j)=unifcdf(lowb(j),lowb(j),upb(j)); p_up(j)=unifcdf(upb(j),lowb(j),upb(j)); p_bound(j)=p_up(j)-p_low(j); p_subbound(j)=p_bound(j)./(m-1); for i=0:m-1 rv=[rv;unifinv(p_low(j)+i.*p_subbound(j),lowb(j),upb(j))]; end elseif dist(j)==2 p_low(j)=normcdf(lowb(j),mu(j),sigma(j)); p_up(j)=normcdf(upb(j),mu(j),sigma(j)); p_bound(j)=p_up(j)-p_low(j); p_subbound(j)=p_bound(j)./(m-1); for i=0:m-1 rv=[rv;norminv(p_low(j)+i.*p_subbound(j),mu(j),sigma(j))]; end elseif dist(j)==3 p_low(j)=logncdf(lowb(j),mu(j),sigma(j)); p_up(j)=logncdf(upb(j),mu(j),sigma(j)); p_bound(j)=p_up(j)-p_low(j); p_subbound(j)=p_bound(j)./(m-1); for i=0:m-1 rv=[rv;logninv(p_low(j)+i.*p_subbound(j),mu(j),sigma(j))]; end elseif dist(j)==4 p_low(j)=evcdf(lowb(j),mu(j),sigma(j)); p_up(j)=evcdf(upb(j),mu(j),sigma(j)); p_bound(j)=p_up(j)-p_low(j); p_subbound(j)=p_bound(j)./(m-1); for i=O:m-1 rv=[rv;evinv(p_low(j)+i.*p_subbound(j),mu(j),sigma(j))]; end end rvcom=[rvcom rv]; end S=[]; for i=1:n S=[S randsample(rvcom(:,i),m)]; end if n==2 figure plot(S(:,1),S(:,2),'k.'); elseif n==3 figure plot3(S(:,1),S(:,2),S(:,3),'lc.'); end