Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

 

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

Contact us