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

chernoff_pdf(x)
% chernoff_pdf.m - evaluates a Chernoff Probability Density.
%   See "Computing Chernoff's Distribution", 
%   P. Groeneboom and J. A. Wellner, 3 November, 1999.
%
%   Vector Form of PDF !!!
%
%  Created by Jim Huntley,  9/12/06.
%
%  Calls:   'pp1.m', 'pp2.m'.
%
%  Inputs:  Requires output files from 'generate_chernoff_pdf.m'.  See below.
%

function [pdf] = chernoff_pdf(x)

%persistent y dy xa xd p

%Initializations.
sx = size(x,2);
%if(isempty(y))
    % Get pre-computed values of 'y' and 'dy'.
    load yvalues y dy;
    % Get pre-computed Airy zeros and Airy derivatives at the Airy zeros.
    load xaxd xa xd;
    % Get values for pre-computed function, 'p(y)'.
    load pofy p;
%end
klim = size(xa,2);

% Estimate the function, g(x).
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

% Estimate the PDFs
pdf(1:sx) = 0;
for jx = 1:sx
    pdf(jx) = 0.5 * gx(jx) * gx(sx+1-jx);
end

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