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

modhans_pdf(x, u2, u4, u6, u8)
% modhans_pdf.m - evaluates a modified 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.
%
%   Note: For a true Hansmann Distribution, the 4 input parameters, u2, u4, u6
%   and u8 are the 2nd, 4th, 6th and 8th moments of a distribution.  The
%   Hansman form of this distribution can then be evaluated from -b to b. 
%   However, since the distribution is not usually known a priori, these 
%   values are not available to the user. Therefore, this distribution has 
%   been modifed to take "arbitrary" positive values of u2 through u8 
%   (decreasing values work best). The modification is that now the distribution 
%   can only be evaluated from -a to a.
%
%   Vector Form of PDF !!!
%
%  Created by Jim Huntley,  12/13/07
%

function [pdf] = modhans_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 = abs(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)   % 2 Unequal Positive Roots.
    k = 1 / (2*c2*(b2-a2))
    pdf1(1:sx) = 0;
    for jx = 1:sx;
        x2 = x(jx)^2;
        pdf1(jx) = ((b2-x2)/(a2-x2))^k;
    end
    I1 = simps(pdf1) * dx;
    K1 = 1 / I1
    pdf = K1.*pdf1 
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