function y=binomial(n,N,p)
%BINOMIAL Generate binomial random numbers
%according to the inversion method.
% n is the number of samples desired.
% N is the number of possible outcomes
% p is the probability the event occurs.
x=rand(1,n);
% Calculate binomial CDF
q = 1 - p;
CDF(1) = q^n; % prob(k) when k = 0
for k = 2:N
   prob(k) = nchoosek(n,k-1) * p^(k-1) * q^(n-k+1);
   CDF(k) = prob(k)+CDF(k-1);
end
% Generate Poisson-distributed random numbers
% using the inversion method.
for i = 1:n
   temp=find(CDF >= x(i));
   if size(temp,2) == 0
      y(i) = 0;
   else
      y(i) = temp(1)-1;
   end
end