Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

dirmultnom_pdf(n, a, b, c)
```function [pdf] = dirmultnom_pdf(n, a, b, c)

%*******************************************************************************
%
%! DIRICHLET_MULTINOMIAL_PDF evaluates a Dirichlet Multinomial PDF.
%
%
%  Formula:
%
%    PDF(X)(A,B,C) = Comb(A,B,X) * ( Gamma(C_Sum) / Gamma(C_Sum+A) )
%      Product ( 1 <= I <= B ) Gamma(C(I)+X(I)) / Gamma(C(I))
%
%    where:
%
%      Comb(A,B,X) is the multinomial coefficient C( A; X(1), X(2), ..., X(B) ),
%      C_Sum = Sum ( 1 <= I <= B ) C(I)
%
%  Reference:
%
%    Kenneth Lange,
%    Mathematical and Statistical Methods for Genetic Analysis,
%    Springer, 1997, page 45.
%
%  Modified:
%
%    17 December 1999 (pdf evaluation for single vector, x),
%    30 March 2011 (pdf evaluation for all permutation vectors of 'b' outcomes
%    in 'a' trials)
%
%  Authors:
%
%    John Burkardt (Fortran 90); Jim Huntley (Matlab)
%
%  Parameters:
%
%    Input, integer X(B); X(I) counts the number of occurrences of
%    outcome I, out of the total of A trials.
%
%    Input, integer A, the total number of trials.
%
%    Input, integer B, the number of different possible outcomes on
%    one trial.
%
%    Input, integer C(B); C(I) is the Dirichlet parameter associated
%    with outcome I.
%
%    Output, real PDF, the value of the Dirichlet multinomial PDF.
%
%   NOTE: SINCE THE ORDER OF THE PERMUTATION VECTORS IS ARBITRARY,
%         THE PDF VALUES ARE SORTED INTO DESCENDING ORDER TO PRODUCE
%         THE SMOOTHEST PDF AND CDF!!!
%
%    Calls: 'dist_nex.m'
%

[d,jmax] = dist_nex(b,a);
c_sum = sum( c(1:b) );
pdf_log0 = - gammaln( c_sum + real( a ) ) + gammaln( c_sum ) ...
+ gammaln( real( a + 1 ) );

for j = 1:jmax
pdf_log = pdf_log0;
x = d(j,:);
for i = 1: b
pdf_log = pdf_log + gammaln( c(i) + real ( x(i) ) ) ...
- gammaln( c(i) ) - gammaln( real ( x(i) + 1 ) );
end
pdf(j) = exp( pdf_log );
end
pdf = sort(pdf,'descend');

return
```