Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

James Huntley (view profile)

 

generates random variates from over 870 univariate distributions

[nt,kf,xa,xb,xc,xd]=airyzo(nt,kf,xa,xb,xc,xd);
function [nt,kf,xa,xb,xc,xd]=airyzo(nt,kf,xa,xb,xc,xd);

%                Ai(x) and Ai'(x), a and a', and the associated
%                values of Ai(a') and Ai'(a); and the first NT
%                zeros of Airy functions Bi(x) and Bi'(x), b and
%                b', and the associated values of Bi(b') and
%                Bi'(b)
%       Input :  NT    --- Total number of zeros
%                KF    --- Function code
%                          KF=1 for Ai(x) and Ai'(x)
%                          KF=2 for Bi(x) and Bi'(x)
%       Output:  XA(m) --- a, the m-th zero of Ai(x) or
%                          b, the m-th zero of Bi(x)
%                XB(m) --- a', the m-th zero of Ai'(x) or
%                          b', the m-th zero of Bi'(x)
%                XC(m) --- Ai(a') or Bi(b')
%                XD(m) --- Ai'(a) or Bi'(b)
%                          ( m --- Serial number of zeros )
%       Routine called: AIRYB for computing Airy functions and
%                       their derivatives
%       =======================================================


x=[];ai=[];bi=[];ad=[];bd=[];
pi=3.141592653589793d0;
for  i=1:nt;
if (kf == 1) ;
u=3.0.*pi.*(4.0.*i-1)./8.0d0;
u1=1./(u.*u);
rt0=-(u.*u).^(1.0./3.0).*((((-15.5902.*u1+.929844).*u1-.138889).*u1+.10416667d0).*u1+1.0d0);
elseif (kf == 2);
if (i == 1) ;
rt0=-1.17371;
else;
u=3.0.*pi.*(4.0.*i-3.0)./8.0;
u1=1.0d0./(u.*u);
rt0=-(u.*u).^(1.0./3.0).*((((-15.5902.*u1+.929844).*u1-.138889).*u1+.10416667).*u1+1.0);
end;
end;
rt=1.0e300;
while (abs((rt-rt0)./rt) > 1.d-9);
x=rt0;
[x,ai,bi,ad,bd]=airyb(x,ai,bi,ad,bd);
if (kf == 1) rt=rt0-ai./ad; end;
if (kf == 2) rt=rt0-bi./bd; end;
if (abs((rt-rt0)./rt) > 1.d-9) ;
rt0=rt;
end;
end;
xa(i)=rt;
if (kf == 1) xd(i)=ad; end;
if (kf == 2) xd(i)=bd; end;
end;
for  i=1:nt;
if (kf == 1) ;
if (i == 1) ;
rt0=-1.01879;
else;
u=3.0.*pi.*(4.0.*i-3.0)./8.0;
u1=1./(u.*u);
rt0=-(u.*u).^(1.0./3.0).*((((15.0168.*u1-.873954).*u1+.121528).*u1-.145833d0).*u1+1.0d0);
end;
elseif (kf == 2);
if (i == 1) ;
rt0=-2.29444;
else;
u=3.0.*pi.*(4.0.*i-1.0)./8.0;
u1=1.0./(u.*u);
rt0=-(u.*u).^(1.0./3.0).*((((15.0168.*u1-.873954).*u1+.121528).*u1-.145833).*u1+1.0);
end;
end;
rt=1.0e300;
while (abs((rt-rt0)./rt) > 1.0d-9);
x=rt0;
[x,ai,bi,ad,bd]=airyb(x,ai,bi,ad,bd);
if (kf == 1) rt=rt0-ad./(ai.*x); end;
if (kf == 2) rt=rt0-bd./(bi.*x); end;
if (abs((rt-rt0)./rt) > 1.0d-9) ;
rt0=rt;
end;
end;
xb(i)=rt;
if (kf == 1) xc(i)=ai; end;
if (kf == 2) xc(i)=bi; end;
end;
return;



function [x,ai,bi,ad,bd]=airyb(x,ai,bi,ad,bd);

%       Input:   x  --- Argument of Airy function
%       Output:  AI --- Ai(x)
%                BI --- Bi(x)
%                AD --- Ai'(x)
%                BD --- Bi'(x)
%       =======================================================


ck=zeros(41,1);dk=zeros(41,1);

eps=1.0d-15;
pi=3.141592653589793d0;
c1=0.355028053887817d0;
c2=0.258819403792807d0;
sr3=1.732050807568877d0;
xa=abs(x);
xq=sqrt(xa);
if (x > 0.0d0) xm=5.0; end;
if (x <= 0.0d0) xm=8.0; end;
if (x == 0.0d0) ;
ai=c1;
bi=sr3.*c1;
ad=-c2;
bd=sr3.*c2;
return;
end;
if (xa <= xm) ;
fx=1.0d0;
r=1.0d0;
for  k=1:40;
r=r.*x./(3.0d0.*k).*x./(3.0d0.*k-1.0d0).*x;
fx=fx+r;
if (abs(r./fx) < eps) break; end;
end;
gx=x;
r=x;
for  k=1:40;
r=r.*x./(3.0d0.*k).*x./(3.0d0.*k+1.0d0).*x;
gx=gx+r;
if (abs(r./gx) < eps) break; end;
end;
ai=c1.*fx-c2.*gx;
bi=sr3.*(c1.*fx+c2.*gx);
df=.5d0.*x.*x;
r=df;
for  k=1:40;
r=r.*x./(3.0d0.*k).*x./(3.0d0.*k+2.0d0).*x;
df=df+r;
if (abs(r./df) < eps) break; end;
end;
dg=1.0d0;
r=1.0d0;
for  k=1:40;
r=r.*x./(3.0d0.*k).*x./(3.0d0.*k-2.0d0).*x;
dg=dg+r;
if (abs(r./dg) < eps) break; end;
end;
ad=c1.*df-c2.*dg;
bd=sr3.*(c1.*df+c2.*dg);
else;
xe=xa.*xq./1.5d0;
xr1=1.0d0./xe;
xar=1.0d0./xq;
xf=sqrt(xar);
rp=.5641895835477563d0;
r=1.0d0;
for  k=1:40;
r=r.*(6.0d0.*k-1.0d0)./216.0d0.*(6.0d0.*k-3.0d0)./k.*(6.0d0.*k-5.0d0)./(2.0d0.*k-1.0d0);
ck(k)=r;
dk(k)=-(6.0d0.*k+1.0d0)./(6.0d0.*k-1.0d0).*ck(k);
end;
km=fix(24.5-xa);
if (xa < 6.0) km=14; end;
if (xa > 15.0) km=10; end;
if (x > 0.0d0) ;
sai=1.0d0;
sad=1.0d0;
r=1.0d0;
for  k=1:km;
r=-r.*xr1;
sai=sai+ck(k).*r;
sad=sad+dk(k).*r;
end;
sbi=1.0d0;
sbd=1.0d0;
r=1.0d0;
for  k=1:km;
r=r.*xr1;
sbi=sbi+ck(k).*r;
sbd=sbd+dk(k).*r;
end;
xp1=exp(-xe);
ai=.5d0.*rp.*xf.*xp1.*sai;
bi=rp.*xf./xp1.*sbi;
ad=-.5d0.*rp./xf.*xp1.*sad;
bd=rp./xf./xp1.*sbd;
else;
xcs=cos(xe+pi./4.0d0);
xss=sin(xe+pi./4.0d0);
ssa=1.0d0;
sda=1.0d0;
r=1.0d0;
xr2=1.0d0./(xe.*xe);
for  k=1:km;
r=-r.*xr2;
ssa=ssa+ck(2.*k).*r;
sda=sda+dk(2.*k).*r;
end;
ssb=ck(1).*xr1;
sdb=dk(1).*xr1;
r=xr1;
for  k=1:km;
r=-r.*xr2;
ssb=ssb+ck(2.*k+1).*r;
sdb=sdb+dk(2.*k+1).*r;
end;
ai=rp.*xf.*(xss.*ssa-xcs.*ssb);
bi=rp.*xf.*(xcs.*ssa+xss.*ssb);
ad=-rp./xf.*(xcs.*sda+xss.*sdb);
bd=rp./xf.*(xss.*sda-xcs.*sdb);
end;
end;
return;

Contact us