Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

neg_pdf(x, lambda, gam)
```% neg_pdf.m - evaluates a Normal Exponential Gamma Probability Density.
%   See "Simultaneous Analysis of all SNPs in a Genome-Wide Association Study -
%   Supplementary material", C. J. Hoggart et al., 28 May, 2008,
%   http://ukpmc.ac.uk/ukpmc/ncbi/articles/PMC2464715/bin/pgen.1000130.s001.pdf
%
%   NOTE:  'lambda' must be >0 but cannot have values of 1/2, 1, 3/2, 2, 5/2, ...
%
%   Calls:  'pcfu.m'.
%
%  Created by Jim Huntley,  08/09/11
%

function [pdf] = neg_pdf(x, lambda, gam)

%persistent gamsq coef a

%if(isempty(coef) == 1)
gamsq = gam^2;
coef = lambda * 2^(lambda+1) * gamma(lambda+0.5) / (gam*sqrt(pi));
a = 2*lambda + 0.5;
%end

argD = abs(x) ./ gam;
D = pcfu(a,argD);
sx = size(x,2);
%Use asymptotic form from A&S to avoid numerical problems for large 'x'.
for jx = 1:sx
if(x(jx)/a >= 5)
D(jx) = exp(-0.25*x(jx)^2) * x(jx)^(-a-0.5) * (1 - 0.5*(a+0.5)*(a+1.5)/x(jx)^2 + ...
0.125*(a+0.5)*(a+1.5)*(a+2.5)*(a+3.5)/x(jx)^4);
end
end
pdf = coef .* exp(0.25.*x.^2./gamsq) .* D;

return```