function PR=Multinomial(N,q,X)
%   N is the total number of ossuaries in the tomb
%   p is a vector of probabilities of K<= names.
%   p must add to less than or equal to one.
%       else PR=666 is returned.   
%   See end for some simple sanity checks
%   X is a vector of length K giving the
%       threshold of each name. There is a default.  
%  
%   
%   The following CASES are from "The Jesus Tomb Math" by William A. Dembski and Robert J. Marks II, a
%   Chapter in Buried Hopes or Risen Saviour (B&H Publishing Group, 2007).
% 
%
%   CASE 1
%   clear all;N=10; D=2509;q=[103*231/(D^2) 2/317 26/317 9/D];q=q/2;p=Multinomial(N,q,0)
%                   p=1.9154e-006
%   CASE 2
%   clear all;N=20; D=2509;q=[103*231/(D^2) 2/317 26/317 9/D];q=q/2;PR=Multinomial(N,q,0)
%                   PR=3.5229e-005
%   CASE 3
%   clear all;N=10; D=2509;q=[103*231/(D^2) 80/317 34/D ];q=q/2;X=[1 2 1 0];PR=Multinomial(N,q,X)
%                   PR=2.9960e-004
%
%   CASE 4
%   clear all;N=20; D=2509;q=[103*231/(D^2) 80/317 34/D ];q=q/2;X=[1 2 1 0];PR=Multinomial(N,q,X)
%                   PR=3.0838e-003
%
K=max(size(q))+1;           % number of named people -1 
if X==0
    X=ones(K,1);
    X(K)==0;
    X=X';
end
if or(sum(q)>1,K~=max(size(X)))
    PR=666;
    return;
end
p=[q 1-sum(q)];
fN=factorial(N);
k=zeros(K,1);
PR=0;
for kk=0:(N^(K))-1 
    for mm=1:K
        k(mm)=k(mm)+1;
        if k(mm) <= N
            break
        else
            k(mm)=0;
        end
    end
    if and(sum(k)==N,k >= X')
        fac=prod(factorial(k));
        PR=PR+(fN/fac)*prod(p'.^k);
    end
end
% %
% %SANITY CHECK #1
% clear all; N=2; q=[1/4 1/2]; PR=tomb(N,q)
% % RETURNS: PR = 0.2500
% % ANALYSIS
% % We have 1,1,0 as the ONLY solution.  1/8 x two combinations = 0.25.
% %
% %SANITY CHECK #2
% clear all; N=3; q=[1/3 1/2]; PR=tomb(N,q)
% % RETURNS: PR = 0.5833
% % Ordering solutions are
% %      2     1     0     
% %      1     2     0
% %      1     1     1
% % With probabilities 
%     PPR=(6/2)*((1/3)^2)*(1/2)+(6/2)*(1/3)*((1/2)^2)+6*(1/2)*(1/3)*(1/6)
% % RETURNS: PPR = 0.5833
% % so PR = PRR
% %
% %SANITY CHECK #3
% clear all; N=5; q=[1/3 0 1/2]; PR=tomb(N,q)
% % RETURNS: PR = 0







