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

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