function pdf = dirichlet_multinomial_pdf ( x, a, b, c )

%*****************************************************************************80
%
%% DIRICHLET_MULTINOMIAL_PDF evaluates a Dirichlet Multinomial PDF.
%
%  Discussion:
%
%    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)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Kenneth Lange,
%    Mathematical and Statistical Methods for Genetic Analysis,
%    Springer, 1997, page 45.
%
%  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, real C(B); C(I) is the Dirichlet parameter associated
%    with outcome I.
%
%    Output, real PDF, the value of the Dirichlet multinomial PDF.
%
  c_sum = sum ( c(1:b) );

  pdf_log = -gammaln ( c_sum + a ) + gammaln ( c_sum ) + gammaln ( a + 1 );

  for i = 1 : b
    pdf_log = pdf_log + gammaln ( c(i) + x(i) ) - gammaln ( c(i) ) ...
      - gammaln ( x(i) + 1 );
  end

  pdf = exp ( pdf_log );

  return
end