Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

negfacmulti_pdf(n, sp1, k, N)
```% negfacmulti_pdf.m - evaluates a Negative Factorial Multinomial Distribution.
%   See "Univariate Discrete Distributions", Johnson, Kemp & Kotz,
%   J. Wiley, 2005, p.522.
%
%   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] = negfacmulti_pdf(n, sp1, k, N)

%persistent NN nnmax nn nmax zeroterm P0

%if(isempty(NN))
NN = sum(N(1:sp1));
nnmax = NN-k;
nn = 0:nnmax;
sumn = 0;
for js = 2:sp1
sumn = sumn + (js-1)*N(js);
end
nmax = sumn;
zeroterm = (N(1) - k + 1) * binomial_coef(N(1),k-1);
P0 = (N(1)-k+1)*N(1) / ((NN-k+1)*(k-1)*binomial_coef(NN,k-1));
%end

sumpdf = 0;
if(n <= nmax)
pdf = P0;
if(n > eps)
sumpdf = P0;
% 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)
% Find subset of rows of 'xs' of length 'xlim' that still satisfy
% the Diophantine Equation
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);
prodbcNr = 1;
for jc = 1:sc
bcNr = binomial_coef(N(jc+1),r(jc));
prodbcNr = prodbcNr * bcNr;
end
sum1 = sum1 + prodbcNr;
end % summation loop over common solutions.
pdf = sum1 * zeroterm / (binomial_coef(NN,k+nn(jnn)-1)*(NN-k-nn(jnn)+1));
sumpdf = sumpdf + pdf;

end % Main Loop over 'nn'.

end % if n > eps.

end % if(n <= nmax)

return```