function PR=MonteCarlo(N,q,X,T)
%   N is the total number of ossuaries in the tomb
%   q is a vector of probabilities of K<= names.
%   q must add to less than or equal to one.
%       else PR=666 is returned. 
%   If there are no successes, PR=777 is returned
%   X is a vector of length K giving the
%       threshold of each name. There is a default.
%   T is the number of Bernoulli trials
%   
%   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;T=10^7;D=2509;q=[103*231/(D^2) 2/317 26/317 9/D];q=q/2;p=MonteCarlo(N,q,0,T)
%                   p=1.9000e-006                            1.9154e-006
%   
%   CASE 2
%   clear all;N=20;T=10^7;D=2509;q=[103*231/(D^2) 2/317 26/317 9/D];q=q/2;p=MonteCarlo(N,q,0,T)
%                   3.5100e-005         vs.                 3.5229e-005

%   CASE 3
%   clear all;N=10;T=10^7; D=2509;q=[103*231/(D^2) 80/317 34/D ];q=q/2;X=[1 2 1 0];p=MonteCarlo(N,q,X,T)
%                   2.9350e-004       vs.                   2.9960e-004
%
%   CASE 4
%   clear all;N=20;T=10^7;D=2509;q=[103*231/(D^2) 80/317 34/D ];q=q/2;X=[1 2 1 0];p=MonteCarlo(N,q,X,T)
%                   3.0963e-006   vs                    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)];
P=cumsum(p);
num=0;
den=0;
for t=1:T
    tomb=zeros(K,1);
    for n=1:N
        r=rand;
        for k=1:K
            if r<P(k);
                tomb(k)=tomb(k)+1;
                break
            end
        end
    end   
    if tomb>=X'
        num=num+1;
    end
    den=den+1;
end
if den>0
    PR=num/den;
else
    PR=777;
end