function pb=mcbinom(M,p,n)
% MCBINOM - cumulative binomial distribution, in a form for MCSSA.
% Syntax:  pb=mcbinom(M,p); pb=mcbinom(M,p,n);
% 
% MCBINOM(M,p,n), where n is an integer, returns the probability that n
% *or more* trials will be successful, given the total number of trials M
% and the probability p that an individual trial is *unsuccessful*.
%
% MCBINOM(M,p) returns a vector pb where the n-th element is the 
% probability of n *or more* trials being successful, given the
% total number of trials M and the probability p that an individual
% trial is *unsuccessful*.
%
% Input: M - number of trials.
%        p - probability of *no* success for an individual trial.
%        n - (optional) number of successful trials.
% 
% In the context of MCSSA, n is the number of significant eigenvalues, M is
% the embedding dimension, and p is the confidence level. The n-th element
% of the output pb is the probability that n or more eigenvalues will be
% found significant at level p in M trials.
%
% For example:  pb=mcbinom(100, .99, 2) yields the probability that 2 or 
% more significant eigenvalues will be found in Monte Carlo testing at the
% 99% level, for an SSA with M=100. These values yield pb = 0.2642, so
% there is a 26% chance of getting two significant eigenvalues with one
% of the surrogates.
%
% This routine is to be used in conjunction with CHI2CONF. The output
% from CHI2CONF and MCBINOM approximates the confidence limits (c and c2)
% which are output from MCSSA1 and MCSSA2. Note that the M tests in
% MCSSA are not strictly independent, so the binomial distribution is 
% not strictly valid - although still useful.
%
% Written by Eric Breitenberger.      Version 3/10/96
% Please send comments and suggestions to eric@gi.alaska.edu       
%

if nargin==2
  n=0:M; 
else
  n=n:M;
end 

p=1-p; % put p in the 'usual' form

% Compute all the binomial coefficients using:
%   bc(M,n)=gamma(M+1)/(gamma(n+1)*gamma(M-n+1))
% (use gammaln rather than gamma to avoid overflow)

bc=round(exp(ones(1,length(n))*gammaln(M+1)-gammaln(n+1)-gammaln(M-n+1))); 

pb=bc.*p.^n.*(1-p).^(M-n); 

if nargin==2
  pb=cumsum(pb(M:-1:1)); 
  pb=pb(M-1:-1:1); 
else
  pb=sum(pb);
end