function pdf = multinomial_pdf ( x, a, b, c ) %*****************************************************************************80 % %% MULTINOMIAL_PDF computes a Multinomial PDF. % % Discussion: % % PDF(X)(A,B,C) = Comb(A,B,X) * Product ( 1 <= I <= B ) C(I)^X(I) % % where Comb(A,B,X) is the multinomial coefficient % C( A; X(1), X(2), ..., X(B) ) % % PDF(X)(A,B,C) is the probability that in A trials there % will be exactly X(I) occurrences of event I, whose probability % on one trial is C(I), for I from 1 to B. % % As soon as A or B gets large, the number of possible X's explodes, % and the probability of any particular X can become extremely small. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 19 October 2004 % % Author: % % John Burkardt % % 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 probability of outcome I on % any one trial. % % Output, real PDF, the value of the multinomial PDF. % % % To try to avoid overflow, do the calculation in terms of logarithms. % Note that Gamma(A+1) = A factorial. % pdf_log = gammaln ( a + 1 ); for i = 1 : b pdf_log = pdf_log + x(i) * log ( c(i) ) - gammaln ( x(i) + 1 ); end pdf = exp ( pdf_log ); return end