Code covered by the BSD License

Generation of Random Variates

James Huntley (view profile)

generates random variates from over 870 univariate distributions

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

function[pdf] = multinom_pdf(n, sp1, nn, p)

%persistent ss gamnnp1

%if(isempty(ss))
ss = sp1 - 1;
gamnnp1 = gammaln(nn+1);
%end

pdf = p(1)^nn;
if(n > eps)
% 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(ss,size(d,2));
xs = int8(d(:,1:xlim));

% 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,:) = int8(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);
%prodpdr = p(1)^r0(jr) / gamma(r0(jr)+1);
prodpdr = (r0(jr)*log(p(1)) - gammaln(r0(jr)+1));
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;
end
sum1 = sum1 + exp(prodpdr);
end % summation loop over common solutions.
pdf = exp(log(sum1) + gamnnp1);

end % if n(jn) > eps.

return```