function p = poisscdf(x,lambda)
%POISSCDF Poisson cumulative distribution function.
%   P = POISSCDF(X,LAMBDA) computes the Poisson cumulative
%   distribution function with parameter LAMBDA at the values in X.
%
%   The size of P is the common size of X and LAMBDA. A scalar input   
%   functions as a constant matrix of the same size as the other input.    
%
%   See also POISSFIT, POISSINV, POISSPDF, POISSRND, POISSTAT.

%   References:
%      [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
%      Functions", Government Printing Office, 1964, 26.1.22.

%   Copyright 1993-2007 The MathWorks, Inc. 
%   $Revision: 2.13.2.5 $  $Date: 2007/05/23 19:16:03 $
 
if nargin < 2, 
    error('stats:poisscdf:TooFewInputs','Requires two input arguments.'); 
end

scalarlambda = isscalar(lambda);

[errorcode x lambda] = distchck(2,x,lambda);

if errorcode > 0
    error('stats:poisscdf:InputSizeMismatch',...
          'Requires non-scalar arguments to match in size.');
end

% Initialize P to zero.
if isa(x,'single') || isa(lambda,'single')
   p = zeros(size(x),'single');
else
   p = zeros(size(x));
end
if ~isfloat(x)
   x = double(x);
end
xx = floor(x);

% Return NaN for invalid or indeterminate inputs
t = (lambda<0) | (xx==Inf & lambda==Inf);
if any(t)
    p(t) = NaN;
end
todo = (x>=0) & ~t & isfinite(lambda);

% Return 1 for x=Inf as long as lambda is valid
t = (x==Inf & lambda>0 & isfinite(lambda));
if any(t)
    todo(t) = false;
    p(t) = 1;
end

% Compute P when X for the remaining cases
k = find(todo);

t = (xx(k)<1e6);
if any(t)
    kt = k(t);
    p(kt) = sumfrom0(xx(kt),lambda(kt),scalarlambda);
end
if any(~t)
     kt = k(~t);
     p(kt) = sumA2B(xx(kt),lambda(kt));
end

% Make sure that round-off errors never make P greater than 1.
p(p>1) = 1;

% ----------------------
function p=sumfrom0(x,lambda,scalarlambda)
val = max(x(:));

if scalarlambda
    tmp = cumsum(poisspdf(0:val,lambda(1)));            
    p = tmp(x + 1);
else
    % Process each unique value of lambda
    p = zeros(size(x),class(x));
    u = unique(lambda);
    for j=1:length(u)
        t = (lambda == u(j));
        val = max(x(t));
        tmp = cumsum(poisspdf(0:val,u(j)));
        p(t) = tmp(x(t)+1);
    end
end


% ----------------------
function p = sumA2B(x,lambda)

p = zeros(size(lambda));

% myfun(k) is approximately the log of poisspdf(k,lambda) (here instead of gammaln,
% we use Stirling's approximation), minus log(eps(0)), explained later.
myfun =@(k,lambda)-lambda+k.*log(lambda)-gammaln(k)-log(eps(0));
% derivative of myfun(k) with respect to k is given by 
mydfun=@(k,lambda) log(lambda)-log(k)-1/2/k;

% If myfun(k) is negative, then poisspdf(k,lambda) is less than eps(0).
% Therefore we can neglect those terms.

% i0: indices when all terms before x in poisspdf(:,lambda) are going to be less than eps
i0=find(myfun(x,lambda)< 0 & x<lambda);
if ~isempty(i0);  % lambda is too large. The result is zero.  
      p(i0) = 0; 
end

% i1: indices when all terms after x in poisspdf are going to be less than eps
i1=find(myfun(x,lambda)< log(eps(1)/eps(0)) & x>lambda);
if ~isempty(i1);  % x is too large. The result is one.  
      p(i1) = 1; 
end

% Either lambda is large enough, and myfun has two zeros (which we call a
% and b, a<b), Or myfun has only one zero, which we call b. We will
% calculate the first zero a only when we are sure it exists and positive.
% If not, we'll set a=0.

% Case 1: lambda is small enough to calculate things straightforwardly:
ilam0 = (lambda < -2*log(eps(0)));    % those indices of lambda small enough to set a=0

% Case 2: lambda is too large: Use asymptotic formula for large lambda:
biglam = 1e10;
ibiglam = (lambda>biglam);
if any(ibiglam)
    p(ibiglam) = .5+0.5*erf((x(ibiglam)-lambda(ibiglam))./sqrt(2*lambda(ibiglam)));
    ismallp = find(abs(p(ibiglam))<=sqrt(eps(1)));
    p(ismallp) = poisspdf(x(ismallp),lambda(ismallp))./(1-x(ismallp)./lambda(ismallp));
end

% Remaining Case: lambda is moderate, then determine the range of terms and
% sum them up:
% ilam1: lambda indices for this (remaining) case
other = ~(ilam0 | ibiglam);
if ~any(other)
    return
end
x = x(other);
lambda = lambda(other);
u = unique(lambda);
other = find(other);

% If lambda is large, then locate the first zero of myfun(k) by Newton's
% iteration:(we will not consider terms before a because they are zero)
a = ones(size(lambda));
for j=1:length(u)
    a0 = 1;
    a1 = myfun(a(j),u(j));        % Residual
    iter = 0;
    while(abs(a1)>1)
        a0 = a0 - a1/mydfun(a0,u(j));   % Newton's iteration step
        a1 = myfun(a0,u(j));            % Residual
        iter = iter+1;
        if iter>100                     % trouble finding lower limit
            a0 = 0;                     % 0 is a conservative choice
            break
        end
    end
    a(lambda==u(j)) = a0;
end

% We will look only at the terms with indices between aa and x.
% Only those contribute, the remaining terms before aa are negligible.
aa=max(floor(a),0);  % a is generally not an integer
le = lasterr;
haderr = false;
for j=1:length(x)
    try
        % try to use summation; risk of running out of memory for large x
        xlist = aa(j):x(j);
        pj = sum(poisspdf(xlist,lambda(j)));
    catch
        % use asympotic formula if necessary
        haderr = true;
        pj = .5+0.5*erf((x(j)-lambda(j))./sqrt(2*lambda(j)));
        if abs(pj)<=sqrt(eps(1))
            pj = poisspdf(x(j),lambda(j))./(1-x(j)./lambda(j));
        end
    end
    p(other(j)) = pj;
end
if haderr
    lasterr(le);
end