# Using optim to find the MLE for the Gamm distribution;  #
#   the default is Nelder-Mead, see help(optim).          #

gnegloglike<-function(pars,data){
#Negative of log-likelihood for the Gamma #
 a<-pars[1]
 b<-pars[2]
 n<-length(data)
 t<--n*lgamma(a)-n*a*log(b)+(a-1)*sum(log(data))-sum(data)/b
 -1*t
 }

# Generate the data. #

x<-rgamma(10,shape=3,scale=2)

# Will usually estimate it pretty well if we start fairly close. #

optim(c(3,4),gnegloglike,data=x)

# But it gets way off if we start far away. #

optim(c(20,20),gnegloglike,data=x)

# So, we could try the contstrained optimization using the #
#   adaptive barrier method.  When you try running it though #
#   you will see it has problems passing the additional #
#   parameters.#

constrOptim(c(20,20),gnegloglike, ui=rbind(c(1,0),c(0,1)), ci=c(0,0),data=x)



# To get around that we can right another function that #
#   defines the data externally.  Of course now you have #
#   to be careful that you aren't storing anything else in x! #

gnegloglike2<-function(pars){
 a<-pars[1]
 b<-pars[2]
 n<-length(x)
 t<--n*lgamma(a)-n*a*log(b)+(a-1)*sum(log(x))-sum(x)/b
 -1*t
 }

constrOptim(c(20,20),gnegloglike2,grad=NULL,ui=rbind(c(1,0),c(0,1)), ci=c(0,0))

constrOptim(c(20,20),gnegloglike2,grad=NULL,ui=rbind(c(1,0),c(0,1)), ci=c(0,0),outer.eps=1e-20,control=list(reltol=1e-20,maxit=5000))

# Notice that the values are slightly different.  It isn't #
#   totally clear why this is the case.  You could try using #
#   optim with the start values equal to the final values of #
#   constrOptim. #