Code covered by the BSD License

Generation of Random Variates

James Huntley (view profile)

generates random variates from over 870 univariate distributions

ipnegbin_pdf(n, ppi, rho, r)
% ipnegbin_pdf.m - evaluates an Inflated Parameter Negative Binomial Probability denisity.
%   See "Inflated-Parameter Family of Generalized Power Series Distributions
%   and their Application in Analysis of Overdispersed Insurance Data",
%   N. Kolev et al. ARCH 00V211.PDF
%
%  Created by Jim Huntley,  10/27/09
%
%   Loads files 'partition1-N.dat'
%

function [pdf] = ipnegbin_pdf(n, ppi, rho, r)

%persistent log1mpr logrho pn

%Initializations.
%if(isempty(pn))
pn = ppi^r;
log1mpr = log((1-ppi)*(1-rho));
logrho = log(rho);
%end

% Evaluate PDF.
if(n == 0)
pdf = pn;

elseif(n > 0)
% Fetch pre-stored partitions of 'n'. Frequencies returned in array, "d".
pname = ['partition' num2str(n)];
spc = size(d,2);
d = double(d);

% Calculate PDF.
% Sum over all partition frequencies to get solution to Diophantine Equation.
sum1 = 0;
d1 = size(d,1);
for jr = 1:d1
for jc = 1:spc
x(jc) = d(jr,jc);
dx(jc) = (jc-1) * x(jc);
%fx(jc) = factorial(x(jc));
fx(jc) = gammaln(x(jc)+1);
end
sumx = sum(x);
%sumxprm1 = sumx + r - 1;
sumxpr = sumx + r;
sumdx = sum(dx);
%pfx = prod(fx)*factorial(r-1);
pfx = sum(fx)+ gammaln(r);
%sum1 = sum1 + factorial(sumxprm1)*((1-ppi)*(1-rho))^sumx*rho^sumdx / pfx;
sum1 = sum1 + exp(gammaln(sumxpr) + sumx*log1mpr + sumdx*logrho - pfx);
end
pdf = pn * sum1;

end % conditional for n > 0.

return