#  R program to calculate maximum likelihood (ML) estimates for parameters
#  k and p in the negative binomial distribution.  The distribution is
#  defined by
#               gamma(k+x)
#  P(X=x) = ------------------ (p^k)[(1-p)^x] , for x=0,1,2,3,...
#           gamma(x+1)gamma(k)
#
#  Here 0<k  and 0<p<1.

#  Count frequencies are entered into the vector ff.  Count values are
#  entered into the vector xx.
ff=c(235,81,43,18,9,6,4,3,0,1)
xx=c(0,1,2,3,4,5,6,7,8,9)
#  Data in example are counts of Plantago major plants in 100 cm X 100
#  cm quadrats, from:
#  J.G. Skellam. 1952. Studies in statistical ecology: I. spatial
#  pattern. Biometrika 39:346-362.

#  Calculates initial values using moment estimates.
nn=sum(ff)
xbar=sum(ff*xx)/nn
ssq=sum(ff*(xx-xbar)*(xx-xbar))/nn
pp0=xbar/ssq
kk0=xbar*xbar/(ssq-xbar)

#  ML objective function "negloglike.ml" is negative of log-likelihood;
#  the Nelder-Mead optimization routine in R, "optim", is a minimization
#  routine.  The three function arguments are:  theta = vector of
#  parameters (transformed to real line), fs = vector of frequencies; and
#  xs = vector of counts.

negloglike.ml=function(theta,fs,xs)  
{
   kk=exp(theta[1])           #  Constrains kk > 0.
   pp=exp(-exp(theta[2]))     #  Constrains 0 < pp < 1.
   lnnb=lfactorial(kk+xs-1)-lfactorial(xs)-lfactorial(kk-1)+
      kk*log(pp)+xs*log(1-pp)
   ofn=-sum(fs*lnnb)
   return(ofn)
}

# The ML estimates.
NBML=optim(par=c(log(kk0),log(-log(pp0))),
   negloglike.ml,NULL,method="Nelder-Mead",fs=ff,xs=xx)
results=c(exp(NBML$par[1]),exp(-exp(NBML$par[2])),-NBML$val)
k.ml=results[1]            # These are the ML estimates.
p.ml=results[2]            #          --
loglike.ml=results[3]

expected=nn*exp(lfactorial(k.ml+xx-1)-lfactorial(xx)-lfactorial(k.ml-1)+
      k.ml*log(p.ml)+xx*log(1-p.ml))
observed=ff

#  Print the results.
k.ml
p.ml
loglike.ml
cbind(expected,observed)