quasirandn_.m

产生正常分布的sobol sequences

function x=quasirandn
%%% purpose:generating of normally distributed 
%%% quasi uniform distributed: sobol sequences
%%%  mean: 0
%%% variance: 1
%%% Moro method employed 
%%% Box-Muller是经典的算法，实际应用中偏向Moro的算法，因为在Quasi Monte
%%% Carlo下，Box-Muller会扭曲结果,本算法采用的是Moro算法

%% initi of the coff  
a=[2.50662823884,-18.61500062529,41.39119773534,-25.44106049637];
b=[1.00,-8.4735109309,23.08336743743,-21.06224101826,3.13082909833];
c=[7.7108870705487895,2.7772013533685169,0.3614964129261002,0.0373418233434554,0.0028297143036967,0.0001625716917922,0.0000080173304740,0.0000003840919865,0.0000000129707170];
k=[0.4179886424926431,4.2454686881376569];

y=halton1(1);
u=y-0.5;

if abs(u)<=0.42
   x1=0;%initi
   x2=0;
   for n=0:3
       x1=x1+a(n+1)*u^(2*n);
   end
   for n=0:4
       x2=x2+b(n+1)*u^(2*n);
   end
   x=u*x1/x2; 
else
   z=k(1)*(2*log(-log(0.5-abs(u)))-k(2));
   d=zeros(1,10);
   d(8)=2*z*d(9)-d(10)+c(8+1);
   d(7)=2*z*d(8)-d(9)+c(7+1);
   d(6)=2*z*d(7)-d(8)+c(6+1);
   d(5)=2*z*d(6)-d(7)+c(5+1);
   d(4)=2*z*d(5)-d(6)+c(4+1);
   d(3)=2*z*d(4)-d(5)+c(3+1);
   d(2)=2*z*d(3)-d(4)+c(2+1);
   d(1)=2*z*d(2)-d(3)+c(1+1);
   fz=z*d(1)-d(2)+c(1)/2;
   if 0.42<u<0.5
      x=fz;
   else
      x=-fz;
   end
end