Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

negmulti_pdf(n, sp1, k, nnmax, p)
```% negmulti_pdf.m - evaluates a Negative Multinomial Distribution.
%   See "Univariate Discrete Distributions", Johnson, Kemp & Kotz,
%   J. Wiley, 2005, p.552.
%
%   Created by: J. Huntley
%
%   NOTE: Text appears to give incomplete description of this pdf.  In
%   order to normalize this PDF, must also sum over all nn ('r' in text).
%   As a result, parameter range is quite limited for reasonable running
%   times.  In general, (sp1-1)*nnmax+1 ~< 35, even though partitions are
%   already pre-stored up to a value of 44. It also helps to have the p's
%   in descending order.
%   Also note that Pr(X=0) = p(1)^k.
%
%   Requires: 'partitionXX.mat', pre-stored partitions created by
%   'create_partition.m'. Note that max(XX) is currently 44.
%

function[sumpdf] = negmulti_pdf(n, sp1, k, nnmax, p)

% Initializations.

nn = 1:nnmax;
pdf0 = p(1)^k;
sumpdf = pdf0;

if(n > eps)
sumpdf = 0;
% Fetch pre-stored partitions of 'n'. Frequencies returned in array, "d".
pname = ['partition' num2str(n)];

% Select rows of 'd' with only min(ss,spc) columns populated and store in array, 'xs'.
% AND that sum to 'nn'.
xlim = min(sp1-1,size(d,2));
xs = int8(d(:,1:xlim));

for jnn = 1:size(nn,2)
zeroterm = exp(gammaln(k+nn(jnn)) + k*log(p(1)) - gammaln(k));
% Find subset of rows of 'xs' of length 'xlim' that still satisfy the Diophantine Equation and
sc = size(xs,2);
sd = size(d,2);
sr = size(xs,1);
ir1 = 0;
xs1 = [];
for jr = 1:sr
if(d(jr,sc+1:end) == 0  & sum(d(jr,1:sc)) == nn(jnn))
ir1 = ir1 + 1;
xs1(ir1,:) = int8(xs(jr,:));
end
if(sc >= sd  & sum(d(jr,1:sc)) == nn(jnn))
ir1 = ir1 + 1;
xs1(ir1,:) = int8(xs(jr,:));
end
end

% Calculate PDF.
% Sum over the subset of solutions to Diaphantine Eq. such that
% the sum of all the 'r's = 'nn'.
sum1 = 0;
sr = size(xs1,1);
sc = size(xs1,2);                         % Summation loop over common (intersecting) solutions.
for jr = 1:sr
r(1:sc) = xs1(jr,1:sc);
%prodpdr = 1;
prodpdr = 0;
for jc = 1:sc
%pdr = p(jc+1)^r(jc) / gamma(r(jc)+1);
pdr = (r(jc)*log(p(jc+1)) - gammaln(r(jc)+1));
%prodpdr = prodpdr * pdr;
prodpdr = prodpdr + pdr;
end
%sum1 = sum1 + prodpdr
sum1 = sum1 + exp(prodpdr);
end % summation loop over common solutions.
pdf = sum1 * zeroterm;
sumpdf = sumpdf + pdf;

end % Main Loop over 'nn'.

end % if n(jn) > eps.

return```