Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

 

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