Code covered by the BSD License

# Generation of Random Variates

### 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```