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

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