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

generate_chernoff_pdf()
% generate_chernoff_pdf.m - generates a Chernoff Probability Density
%   and saves the values necessary for 'chernoff_pdf.m' to operate.
%   See "Computing Chernoff's Distribution", 
%   P. Groeneboom and J. A. Wellner, 3 November, 1999.
%
%  Created by Jim Huntley,  9/12/06.
%
%  Calls:  'pp1.m', 'pp2.m', 'airyzo.m'.
%
%  Inputs:  See Initializations.
%
%  Outputs: y - integration positions.
%           dy - integration step size.
%           xa - zeros of the Airy function.
%           xd - values of the derivative at locations of Airy function zeros.
%           pofy - values of p(y).
%           All outputs needed in order to estimate PDF using 'chernoff_pdf.m'!!!
%

function[] = generate_chernoff_pdf()

% Initializations.
xlim = 5;                       % Maximum value of position at which to estimate PDF.
xsteps = 2000;                  % Number of positions for PDF estimation.
dx = xlim / xsteps;             % Step size for PDF-estimation positions.             
x = -xlim:dx:xlim;              % PDF-estimation positions.
sx = size(x,2);
pdf(1:sx) = 0;                  
ylim = 50;                      % Integration maximum value (~ inf).
ysteps = 2000;                  % Number of integration steps.
dy = ylim / ysteps;             % Integration step size.
y = 0:dy:ylim-dy;               % Integration positions.
logpi = log(pi);

% Save 'y' and 'dy' values needed by 'chernoff.m'.
save yvalues y dy;
klim = 21;                      % Number of summation terms to keep (~ inf).
p(1:ysteps) = 0;

% First, Construct c's.
c(1) = 1;
for jk = 2:klim
    c(jk) = -(2*(jk-1)-3) * (2*(jk-1)+1) * c(jk-1) / ((jk-1)^2 * (2*(jk-1)-1) * 2^4);    
end
c

% Next, Construct a's and b's.
a(1) = 1;
b(1) = 2/3;
for jk = 2:klim
    suma = 0;
    sumb = 0;
    for kk = 1:jk-1
        %suma = suma + b(jk-kk)*beta(3*(jk-1)-2*(kk-1)-0.5,kk+0.5)/(pi*factorial(kk-1)*(-2)^(kk-1)); 
        suma = suma + exp(log(b(jk-kk))+log(beta(3*(jk-1)-2*(kk-1)-0.5,kk+0.5))-(logpi+gammaln(kk))) / (-2)^(kk-1);
        a(jk) = c(jk) - real(suma);
        %sumb = sumb + a(jk-kk+1)*beta(3*(jk)-2*(kk-1)-2,kk+0.5)/(factorial(kk-1)*(-2)^(kk));
        sumb = sumb + real(exp(log(a(jk-kk+1))+log(beta(3*(jk)-2*(kk-1)-2,kk+0.5))-gammaln(kk)) / (-2)^(kk));
    end
    %b(jk) = sumb + a(1)*beta(3*jk-2*(jk-1)-2,jk+0.5) / (factorial(jk-1)*(-2)^(jk));
    b(jk) = sumb + real(exp(log(a(1))+log(beta(3*jk-2*(jk-1)-2,jk+0.5)) - gammaln(jk)) / (-2)^(jk));
end
a
b

%Next, Find and Save the Zeros of the Airy Function as well as derivative values at these zeros.
xa(1:klim) = 0;
xb(1:klim) = 0;
xc(1:klim) = 0;
xd(1:klim) = 0;
for jk = 1:klim
    [nt,kf,xa,xb,xc,xd] = airyzo(klim,1,xa,xb,xc,xd);
end
xa
xd
save xaxd xa xd

% Next, Construct and save p's.
for jy = 1:ysteps
    if(y(jy) <= 1)
        sum1 = 0;
        sum2 = 0;
        for jk = 1:klim
            sum1 = sum1 + a(jk)*y(jy)^(3*(jk-1));
            sum2 = sum2 + b(jk)*y(jy)^(3*(jk-0.5));
        end
        p(jy) = -sqrt(0.5*pi)*sum1 + sum2;
    elseif(y(jy) > 1)
        sum = 0;
        for jk = 1:klim
            sum = sum + exp(2^(1/3)*xa(jk)*y(jy));
        end
        p(jy) = -y(jy)^(-1.5) + 2*sqrt(2*pi)*exp(-(y(jy)^3)/6)*sum;
    end
end
p
save pofy p

% Finally, Estimate and plot the PDF.
gx(1:sx) = 0;
for jx = 1:sx
    xj = x(jx);
    if(xj >= -1)
        z1 = pp1(y,p,xj,jx);
        z2 = pp2(y,xj,jx);
        gx(jx) = 2*xj - dy*trapz(z1)/sqrt(2*pi) + 2*sqrt(2/pi)*dy*trapz(z2);
    elseif(xj < -1)
        sum = 0;
        for jk = 1:klim
            sum = sum + exp(-2^(1/3)*xa(jk)*xj)/xd(jk);
        end
        gx(jx) = exp(2*xj^3/3)*4^(1/3)*sum;
    end
end

for jx = 1:sx
    pdf(jx) = 0.5 * gx(jx) * gx(sx+1-jx);
end

figure(1)
plot(x,pdf);

ipdf = simps(pdf) * xlim / xsteps

return

function [z1] = pp1(y,p,xj,jx)

z1 = p .* exp(-0.5.*y.*(2.*xj+y).^2);

return

function [z2] = pp2(y,xj,jx)

z2 = ((2.*xj+y.^2).*y.^2 + 0.5.*(2.*xj+y.^2).^2) .* exp(-0.5.*y.^2.*(2.*xj+y.^2).^2);

return

Contact us