Code covered by the BSD License

# Generation of Random Variates

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