Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

factmult_pdf(n, sp1, nn, N)
```% factmult_pdf.m - tests a Factorial Multinomial Distribution
%   See "Univariate Discrete Distributions", Johnson, Kemp & Kotz,
%   J. Wiley, 2005, p.523.
%
%   Created by: J. Huntley, 11/20/08
%
%

function[pdf] = factmult_pdf(n, sp1, nn, N)

%persistent ss NN denom

%if(isempty(ss))
ss = sp1 - 1;
NN = sum(N(1:ss+1));
denom = 1 / binomial_coef(NN,nn);
%end

pdf = 0;
if(n > eps)
% Fetch pre-stored partitions of 'n'. Frequencies returned in array, "d".
pname = ['partition' num2str(n)];
spc = size(d,2)
spr = size(d,1);
xlim = min(ss,spc);
d = double(d)

% Select rows of 'd' with only min(nn,spc) columns populated and store in array, 'xs'.
% AND that sum to 'nn' or less.
xs = d(:,1:xlim);
%xsn = xs

% 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)
ir1 = ir1 + 1;
xs1(ir1,:) = xs(jr,:);
end
if(sc >= sd  & sum(d(jr,1:sc)) <= nn)
ir1 = ir1 + 1
xs1(ir1,:) = xs(jr,:);
end
end

% Estimate 'r0'.
r0 = [];
for jr = 1:ir1
sumd = sum(xs1(jr,1:xlim));
r0(jr) = max(0,nn-sumd);
end

%xs1n = xs1
%r0

% Calculate PDF.
% Sum over the subset of solutions to Diaphantine Eq. such that
% the sum of all the 'r's (including 'r0') = '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 = binomial_coef(N(1),r0(jr));
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 * denom;

end % if n(jn) > eps.

return```