set.seed(121, "Mersenne-Twister") ## Set the rng seed
## We will need stats4 library
library(stats4)

## Define function newFit
newFit <- function(isi, ## The data set
                   isi.density = "gamma", ## The name of the density used as model
                   initial.para = list( shape = (mean(isi)/sd(isi))^2,
                     scale = sd(isi)^2 / mean(isi) ), ## Inital parameters passed to optim
                   optim.method = "BFGS", ## optim method
                   optim.lower = numeric(length(initial.para)) + 0.00001,
                   optim.upper = numeric(length(initial.para)) + Inf,
                   ...) {

  require(stats4)

  para.names <- paste(names(initial.para)[1], " = ", names(initial.para)[1], sep = "")
  
  if (length(names(initial.para)) > 1)
    for ( this.name in names(initial.para)[-1] ) {
      para.names <- paste(para.names,
                          ", ",
                          this.name,
                          " = ",
                          this.name,
                          sep = ""
                          )
    }
    
  ## Create a string with the log likelihood definition
  minusLogLikelihood.txt <- paste("function( ",
                                  paste(names(initial.para), collapse = ", "),
                                  " )",
                                  "-sum(",
                                  paste("d", isi.density,
                                        "(",
                                        deparse(substitute(isi)),
                                        ",",
                                        para.names,
                                        ", log = TRUE) )",
                                        sep = ""
                                        )
                                  )

  ## Define logLikelihood function
  minusLogLikelihood <- eval( parse(text = minusLogLikelihood.txt) )
  environment(minusLogLikelihood) <- .GlobalEnv

 
  if ( all( is.infinite( c(optim.lower,optim.upper) ) ) ) {
      getFit <- mle(minusLogLikelihood,
                    start = initial.para,
                    method = optim.method,
                    ...
                    )
  } else {
    getFit <- mle(minusLogLikelihood,
                  start = initial.para,
                  method = "L-BFGS-B",
                  lower = optim.lower,
                  upper = optim.upper,
                  ...
                  )
  }  ## End of conditional on all(is.infinite(c(optim.lower,optim.upper)))
 
  getFit
 
}
## End of function newFit definition

isi1 <- rgamma(100, shape = 2, scale = 1) ## Generate 100 rn with a Gamma distribution
## Define function llgama which return the oposite of the log-likelihood of isi1 assuming
## a Gamma density with parameters shape = shape and scale = scale
llgamma <- function(shape, scale) -sum(dgamma(isi1, shape = shape, scale = scale, log = TRUE))
## Call mle to perform the fit using the method of moments for the initial guess
fitA <- mle(llgamma, start = list( shape = (mean(isi1)/sd(isi1))^2, scale = sd(isi1)^2 / mean(isi1) ),lower = c(0.00001,0.00001), method = "L-BFGS-B")
## Do in principle the same thing by calling newFit
fit1 <- newFit(isi = isi1, isi.density = "gamma", initial.para = list( shape = (mean(isi1)/sd(isi1))^2, scale = sd(isi1)^2 / mean(isi1) ), optim.method = "L-BFGS-B", optim.lower = c(0.00001,0.00001))
