Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

ncchisqpr_pdf(x, lambda1, lambda2, nu1, nu2)
```% ncchisqpr_pdf.m - evaluates a Non-Central Chi-Square Product Probability Density.
%   See "Continuous Univariate Distributions", v.2, Johnaon, Kotz, & Balakrishnan,
%   J. Wiley, 1995, p.472.
%
%  Vector form of PDF!!!
%
%  Created by Jim Huntley,  8/03/07
%

function [pdf] = ncchisqpr_pdf(x, lambda1, lambda2, nu1, nu2)

%persistent coef ohl1 ohl2 ohn1 ohn2 ofn1pn2 ohn1mn2 i1max i2max

%if(isempty(coef))
coef = exp(-0.5*(lambda1+lambda2));
ohl1 = 0.5 * lambda1;
ohl2 = 0.5 * lambda2;
ohn1 = 0.5 * nu1;
ohn2 = 0.5 * nu2;
ofn1pn2 = 0.25*(nu1+nu2) - 1;
ohn1mn2 = 0.5*(nu1-nu2);
i1max = fix((ohn1+ohn2)*5 + eps);
i2max = fix((ohn1+ohn2)*5 + eps);
%end

for jx = 1:size(x,2)
sum1 = 0;
rt2x = sqrt(2*x(jx));
ohx = 0.5*x(jx);
for i1 = 1:i1max
sum2 = 0;
i1m1 = i1-1;
for i2 = 1:i2max
i2m1 = i2-1;
%pr = ohl1^i1m1 * ohl2^i2m1 / (factorial(i1m1) * factorial(i2m1));
pr = exp(i1m1*log(ohl1) + i2m1*log(ohl2) - (gammaln(i1) + gammaln(i2)));
num = ohx^(ofn1pn2+0.5*(i1m1+i2m1)) * besselk(ohn1mn2+i1-i2,rt2x);
denom = gammaln(ohn1+i1m1) + gammaln(ohn2+i2m1);
sum2 = sum2 + exp(log(pr) + log(num) - denom);
end
sum1 = sum1 + sum2;
end
pdf(jx) = coef * sum1;
end

return
```