function [prob] = MultinomialDistribution(P, N);
%[prob] = MultinomialDistribution(P, N)
%   Evaluate the probability mass function of the multinomial
%   distribution with parameter P (vector) at N (integer vector of
%   same size as P) and return the scalar result. 
%   ($RCSfile: MultinomialDistribution.m,v $ $Revision: 1.1.1.1 $)

%%% Check arguments:
P=P(:); N=N(:);
if ( sum([P ; N] <0) )
  error('MultinomialDistribution: arguments must be nonnegative.')
end
if (sum(P) ~= 1)
  error('MultinomialDistribution: Sum of P must be equal to 1.')
end
if ( length(P)~=length(N) )
  error('MultinomialDistribution: P and N must have same length.')
end
if ( sum(round(N)~=N) )
  error('MultinomialDistribution: N must be an integer vector.')
end
%%% Do computation:
NN=sum(N); num=[1:NN]; den=zeros(1,NN); fak2=den; index=1;
for i=1:length(N)
  if (N(i)>0)
    den(index:index+N(i)-1)  = [1:N(i)];
    fak2(index:index+N(i)-1) = P(i)*ones(1,N(i));
    index=index+N(i);
  end
end
prob=prod((num./den).*fak2);
return;
