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

hansmann_pdf(x, u2, u4, u6, u8)
% Hansmann_pdf.m - evaluates a Hansmann Probability Density.
%   see "A Note on Hansmann's 1934 Family of Distributions", 
%   R.F. Pawula and S. O. Rice, IEEE Trans Info Theory, v.35, No.4, p.910 8/1989.
%
%   Vector Form of PDF !!!
%
%  Created by Jim Huntley,  12/13/07
%

function [pdf] = hansmann_pdf(x, u2, u4, u6, u8)

%persistent c2 b2 a2 a dx

%if(isempty(c2))
    del = 9*u8*(5*u4-9*u2^2) - 49*u6^2 + 210*u2*u4*u6 - 125*u4^3;
    c0 = 2*(14*u2*u6^2 - 9*u2*u4*u8 - 5*u4^2*u6) / del;
    c1 = (9*u8*(3*u2^2-u4) + 7*u6^2 - 50*u2*u4*u6 + 25*u4^3) / del;
    c2 = 2*(5*u2*u4^2 -u6*(6*u2^2-u4)) / del;
    c = [c2 c1 c0];
    r = roots(c);
    b2 = max(r);
    a2 = min(r);
    a = sqrt(a2);
    dx = x(2)-x(1);
%end

sx = size(x,2);
pdf(1:sx) = 0;

if(abs(b2-a2) <= eps && b2 > 0)    %Equal Positive Roots.
    for jx = 1:sx
        pdf(jx) = exp(1/(2*c2*(b2-x(jx)^2)));
    end
    I = simps(pdf) * dx;
    K = 1 / I;
    pdf = pdf .* K;
elseif(abs(b2-a2) > eps && b2 > 0)   %Unequal Positive Roots.
    k = 1 / (2*c2*(b2-a2));
    pdf1(1:sx) = 0;
    pdf2(1:sx) = 0;
    x2pdf1(1:sx) = 0;
    x2pdf2(1:sx) = 0;
    for jx = 1:sx;
        x2 = x(jx)^2;
        if(abs(x(jx)) <= a) 
            pdf1(jx) = ((b2-x2)/(a2-x2))^k;
            x2pdf1(jx) = x2 * pdf1(jx);
        elseif(abs(x(jx)) > a)
            pdf2(jx) = ((b2-x2)/(x2-a2))^k;
            x2pdf2(jx) = x2 * pdf2(jx);
        end
    end
    I1 = simps(pdf1) * dx;
    I2 = simps(pdf2) * dx;
    I3 = simps(x2pdf1) * dx;
    I4 = simps(x2pdf2) * dx;
    K1 = (I4 - I2*u2) / (I1*I4 - I2*I3);
    K2 = (1 - K1*I1) / I2;
    pdf = K1.*pdf1 + K2.*pdf2;
elseif(b2 > 0 && a2 < 0)    %Single Positive Root.
    for jx = 1:sx
        k = 1 / (2*c2*(b2-a2));
        x2 = x(jx)^2;
        pdf(jx) = ((b2-x2)/(x2-a2))^k;
    end
    I = simps(pdf) * dx;
    K = 1 / I;
    pdf = pdf .* K;
end

return

Contact us