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

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