function y = cbinom(n,p,k)
% CBINOM y = cbinom(n,p,k) Cumulative binomial probabilities
% Version of 10/5/93
% n is a positive integer; p is a probability
% k is a matrix of integers between 0 and n
% y = P(X>=k) (a matrix of probabilities)
if p > 0.5
  a  = [1 ((1-p)/p)*ones(1,n)];
  b  = [1 n:-1:1];
  c  = [1 1:n];
  br = (p^n)*cumprod(a.*b./c);
  bcr = cumsum(br);
  bc = fliplr(bcr);
else
  a  = [1 (p/(1-p))*ones(1,n)];
  b  = [1 n:-1:1];
  c  = [1 1:n];
  bi = ((1-p)^n)*cumprod(a.*b./c);
  b  = cumsum(bi);
  bc = [1 1-b(1:n)];
end
y = bc(k+1);