R/fitpareto.R

In piLaboratory/sads: Maximum Likelihood Models for Species Abundance Distributions

fitpareto <- function(x, trunc, start.value, upper = 20, ...){
  dots <- list(...)
  if (any(x <= 0)) stop ("All x must be positive")
  if (!missing(trunc)){
    if (min(x)<=trunc) stop("truncation point should be lower than the lowest data value")
  }
  if(missing(start.value)){
    alpha <- length(x)/sum(log(x)-log(min(x))) ## actually the analytical mle of shape
  } else{
    alpha <- start.value[1]
  }
  if (missing(trunc)){
    LL <- function(shape,scale) -sum(dpareto(x, shape, scale, log = TRUE))
  } else {
    LL <- function(shape, scale) -sum(dtrunc("pareto", x = x,
                                             coef = list(shape = shape, scale = scale), trunc = trunc, log = TRUE))
  }  
  ##result <- do.call("mle2", c(list(LL, start = list(shape = alpha, scale=min(x)), fixed=list(scale=min(x)),
  ##                                 data = list(x = x), method = "Brent", lower = 0, upper = upper), dots))
  ## With both parameters free, because the fit has two degrees of freedom
  result <- do.call("mle2", c(list(LL, start = list(shape = alpha, scale=min(x)),
                                   data = list(x = x), method = "L-BFGS-B", upper= c(shape = upper, scale=min(x))), dots))
  if(abs(as.numeric(result@coef[1]) - upper) < 0.001) 
    warning("mle equal to upper bound provided. \n Try value for the 'upper' argument")
  new("fitsad", result, sad="pareto", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc)) 
}