Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

poispar_pdf(n, mu, alpha, bet)
```% poispar_pdf.m - evaluates a Poisson (Generalized) Pareto Probability Density.
%   See "On recursive Evaluation of Mixed Poisson Probabilities and Related Quantities",
%   Scand. Actuarial J., vol.2, p.114, 1993.
%
%  Created by Jim Huntley,  04/22/09
%

function [pdf] = poispar_pdf(n, mu, alpha, bet)

%persistent p0 p1 pnm1 pnm2
persistent pnm1 pnm2

%if(isempty(p0))
if(n < 3)
%p0 = gamma(alpha+bet) * real(KummerU(bet,1-alpha,mu)) / gamma(alpha);
p0 = exp(gammaln(alpha+bet) + log(real(KummerU(bet,1-alpha,mu))) - gammaln(alpha));
%p1 = mu * gamma(alpha+bet) * gamma(1+bet) * real(KummerU(1+bet,2-alpha,mu)) / (gamma(alpha)*gamma(bet));
p1 = mu * exp(gammaln(alpha+bet) + gammaln(1+bet) + log(real(KummerU(1+bet,2-alpha,mu))) - gammaln(alpha) - gammaln(bet));
pnm1 = p1;
pnm2 = p0;
end
%end

if(n == 0)
pdf = p0;
elseif(n == 1)
pdf = p1;
elseif(n > 1)
c1 = (n-1-alpha-mu) / (n);
c2 = mu * (n+bet-2) / ((n-1)*(n));
pdf = c1 * pnm1 + c2 * pnm2;
pnm2 = pnm1;
pnm1 = pdf;
% Check Recursion Relation.
%pdf(jn) = gamma(alpha+bet) * gamma(n(jn)+bet) * mu^(-bet) * ...
%          hypergeom([n(jn)+bet,alpha+bet],[],-1/mu) / ...
%          (gamma(alpha)*gamma(bet)*gamma(n(jn)+1));
end

return```