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

airy_pdf(x)
% airy_pdf.m - evaluates an Airy Probability Density.
%   See "Airy Distribution Function: From the Area under a Brownian
%   Excursion to the Maximal Height of Fluctuating Interfaces", S.N. Majumdar
%   & A.Comtet, http://arxiv.org/PS_cache/cond-mat/pdf/0409/0409566v2.pdf, 1994.
%
%   Calls:  'airyzo.m'
%
%  Created by Jim Huntley,  06/07/11
%

function [pdf] = airy_pdf(x)

%Initializations.
kmax = 5;                       % Converges after very few terms
coef = 2*sqrt(6);
ai = [];
bi = [];
ad = [];
bd = [];
kf = 1;
a = -5/6;
b = 4/3;
gam1 = gamma(1-b) / gamma(a-b+1);
gam2 = gamma(b-1) / gamma(a);   
[kmax,kf,ai,bi,ad,bd] = airyzo(kmax,kf,ai,bi,ad,bd);
ai = abs(ai)
bk = 2 .* ai.^3 ./ 27
sx = size(x,2);

sumk(1:sx) = 0;
for jk = 1:kmax
    argz = bk(jk) ./ x.^2;
    M1 = KummerComplex([a],[b],argz);
    M2 = KummerComplex([a-b+1],[2-b],argz);
    U = real(M1 .* gam1 + argz.^(1-b) .* M2 .* gam2);
    sumk = sumk + exp(-bk(jk)./x.^2) .* bk(jk)^(2/3) .* U;
end      
pdf = coef .* sumk ./ x.^(10/3);
for jx = 1:sx
    if(~isfinite(pdf(jx)))
        pdf(jx) = exp(-bk(1)/x(jx)^2) / x(jx)^5;
    end
    if(pdf(jx) < 0)
        pdf(jx) = 0;
    end
end

return

Contact us