multpdf(n,p,x)

function y = multpdf(n,p,x)
% MULTPDF Multinomial probability density function.
%  Y = MULTPDF(N,P,X) returns the multinomial probability density function value with parameters
%  N and P at the values in X.
%  Note that the density function is zero unless X is an integer.
%  Let {X1, X2,  . . .  , Xk}, k > 1, be a set of random variables, each of which can take 
%  the values 0, 1,  . . .  , n. Suppose there are k nonnegative numbers {p1, p2,  . . .  , pk}
%  that sum to one, such that for every set of k nonnegative integers {n1,  . . .  , nk} whose
%  sum is n, P( X1 = n1 and X2 = n1 and  . . .  and Xk = nk ) = 
%
%                          n!    
%                 ----------------------    p1^n1  p2^n2   . . .   pk^nk .  
%                n1!  n2!   . . .   nk!    
%
%  Then {X1, X1,  . . .  , Xk} have a multinomial joint distribution with parameters n and 
%  p1, p2,  . . .  , pk. The parameter n is called the number of trials; the parameters p1, 
%  p2,  . . .  , pk are called the category probabilities; k is called the number of categories. 
%
%
%  Syntax: function y = multpdf(n,p,x) 
%      
%  Inputs:
%       n - number of trials.
%       p - vector of associated probabilities. 
%       x - vector of the interested values.
%  Outputs:
%       y - multinomial probability density value.
%
%  From the SticiGui:Statistics Tools for Internet and Classroom Instruction with a Graphical
%  User Interface from the University of California at Berkeley. The Exercise 1 from Chapter 23
%  [http://www.stat.berkeley.edu/users/stark/SticiGui/Text/ch23.htm#solution_1]. The positions
%  on a roulette wheel are divided into three colors: red, black, and green. There are 18 red 
%  positions, 18 black positions, and two green positions. If the wheel is fair, the chance that
%  the ball lands in any position is equal to the chance that it lands in any other position, 
%  1/38.
%  The probability that in 7 independent spins of the wheel, the ball lands in a red slot four 
%  times, in a black slot two and in a green slot once is:
%
%  Calling on Matlab the function: 
%             y = multpdf(n,p,x)
%
%  where n = 7, x = [4,2,1]  and  p = [18/38,18/38,2/38]
%
%  Answer is:
%
%  y = 0.0624
%
%  Created by A. Trujillo-Ortiz, R. Hernandez-Walls and A. Castro-Perez
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             [e-mail]
%  Copyright (C) January 2, 2005.
%
%  To cite this file, this would be an appropriate format:
%  Trujillo-Ortiz, A., R. Hernandez-Walls and A. Castro-Perez. (2005). multpdf:
%    Multinomial probability density function. A MATLAB file. [WWW document]. URL
%    http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=6785
%
%  References:
% 
%  Weisstein, E. Mathworld:the Web's most extensive mathematics resource. A Wolfram Web Resource. 
%            http://mathworld.wolfram.com/MultinomialDistribution.html 
%

if nargin < 3, 
   error('You need to input three arguments.');
   return,
end;

if (length(n)~=1) | (fix(n) ~= n) | (n < 0),
   error('n must be a positive integer.');
   return,
end;
    
mbint(x);

k = length(x);
l = length(p);

if k ~= l, 
   error('The input arguments have unequal size.');
   return,
end;

pp = sum(p);

if sum(x) ~= n
    error('Inputs must satisfy n = x1 + x2 ... + xi.');
    return,
end;

if pp ~= 1,
    error('The sum of the input probabilities must be equal 1.')
    return
end;

cc = p.^x;
cp = prod(cc);
c = length(cc);

factor1 = sum(log(1:n));

f = [];
for i = 1:c,
    f = [f log(1:x(i))];
end;

factor2 = sum(f);

y = exp(factor1-factor2)*cp;

return,