bisect<-function(x,fn,a,b,tol=1e-10){
##
#This function performs the univariate bisection #
#  method to find the root of an equation. It begins#
#  with the boundaries a<b.   The function fn #
#  must be of the form f(x,a) where x is the data #
#  and a is the single parameter. #
#It outputs all of the iterations as a matrix. #
#                           BH 03/03/04 #
##
  i<-0
  cur<-(a+b)/2
  res<-c(i,cur,fn(x,cur))
  if ((fn(x,a)*fn(x,b)>0)|(a>b)) {
    return('Improper start values') }
  else
    while (abs(fn(x,cur))>tol) {
      i<-i+1
      if (fn(x,a)*fn(x,cur)>0) {
        a<-cur
        cur<-(a+b)/2 
      }
      else {
        b<-cur
        cur<-(a+b)/2 
      }
      res<-rbind(res,c(i,cur,fn(x,cur)))  
  }
return(res)
}


gammamle<-function(x,a){
##
#The equation whose root is the MLE of the #
#  a parameter. #
##
  log(a)-digamma(a)-log(mean(x))+mean(log(x))}


x<-rgamma(10,shape=3,scale=2)
bisect(x,gammamle,0.00001,10.00001,tol=10e-16)



newton<-function(x,fn,fprime,start,tol=1e-10){
##
#This function performs the univariate Newton's #
#  method to find the root of an equation. It begins#
#  with the initial value start.   The function fn #
#  must be of the form f(x,a) where x is the data #
#  and a is the single parameter. #
#It outputs all of the iterations as a matrix. #
#                           BH 03/03/04 #
##
  i<-0
  cur<-start
  res<-c(i,cur,fn(x,cur))
  while (abs(fn(x,cur))>tol) {
    i<-i+1
    cur<-cur-fn(x,cur)/fprime(x,cur)
    res<-rbind(res,c(i,cur,fn(x,cur))) 
    }  
  return(res)
}



gammamlep<-function(x,a){
##
#The derivative of the equation whose root  #
#  is the MLE of a gamma's a parameter. #
##
  1/a-trigamma(a)}

newton(x,gammamle,gammamlep,5,tol=10e-16) 



#Starting values matter!#

x<-rgamma(10,shape=3,scale=2)
moma<-mean(x)^2/var(x)
newton(x,gammamle,gammamlep,moma,tol=10e-16)
newton(x,gammamle,gammamlep,5,tol=10e-16)