tests/t-mledist-cvg-NelderMead.R

In aursiber/fitdistrplus: Help to Fit of a Parametric Distribution to Non-Censored or Censored Data

if(FALSE)
{
  
  
  x <- c(37.37, 40.97, 30, 23.1, 22.97, 24.95, 17.73, 18.63, 18.22, 32.65, 56.62, 69.74, 
         59.18, 54.47, 37.91, 24.48, 20.26, 20.51, 16.54, 18.94, 19.47, 20.96, 34.02, 
         42.58, 32.74, 33.62, 24.58, 19.55, 17.94, 18.45, 18.8, 18.26, 19.1, 24.01, 30.32, 
         62.84, 87.75, 54.34, 36.42, 25.07, 19.62, 19.32, 19.53, 16, 15.38, 21.89, 27.31, 
         26.73, 38.58, 33.54, 24.4, 21.59, 19.23, 23.07, 20.59, 20.29, 24, 31.66, 52.58, 
         47.07, 72.04, 45.1, 32.38, 28.74, 24.2, 22.25, 22.47, 22.66, 21.09, 28.31, 45.42, 
         47.87, 63.96, 55.67, 35.57, 27.73, 23.42, 21.03, 19.67, 19.53, 18.46, 19.99, 29.49, 
         48.55, 46.29, 43.74, 31.97, 19.71, 17.39, 16.12, 15.62, 14.78, 13.17, 12.01, 19.8, 
         22.41, 32.57, 34.15, 24.59, 14.41, 15.71, 16.57, 15.65, 16.19, 17.66, 19.62, 30.68, 
         64.23, 45.28, 40.55, 28.2, 20.05, 17.5, 15.58, 14.22, 15.29, 15.28, 16.22, 23.38, 
         29.99, 24.76, 28.1, 22.36, 16.23, 12.2, 11.05, 14.67, 14.64, 17.99, 28.87, 45.8, 
         50.56, 36.44, 22.78, 16.27, 24.21, 24.73, 18.35, 23.87, 25.44, 29.87, 48.47, 98.02, 
         113.65, 83.34, 51.41, 32.58, 30.93, 25.84, 26.48, 23.95, 23.05, 24.64, 43.15, 57.83, 
         59.68, 48.1, 32.01, 24.76, 22.56, 20.28, 22.45, 21.74, 24.45, 31.74, 46.12, 54.64, 
         77.71, 46.13, 28.79, 24.26, 22.97, 24.38, 20.95, 21.22, 25.5, 29.55, 37.22, 64.53, 
         61.69, 43.41, 31.6, 26.02, 20.05, 21.81, 19.43, 19.24, 19.19, 20.51, 33.84, 41.1, 
         43.41, 29.16, 28.54, 19, 20.85, 18.18, 17.13, 18.9, 23.91, 38.15, 49.57, 76.37, 
         62.73, 46.32, 30.45, 21.25, 20.78, 20.78, 25.02, 21.85, 24.08, 33.47, 44.53, 89.55, 
         97.05, 59.12, 37.81, 27.99, 25.57, 21.53, 20.15, 23.23, 23.42, 27.61, 40.73, 60.87, 
         94.03, 68.09, 44.47, 40.25, 35.11, 31.67, 32.65, 31.08, 31.13, 35.64, 47.05, 90.5, 
         73.74, 66.07, 46.82, 26.75, 33.9, 25.7, 24.02, 24.19, 23.06, 27.73, 41.85, 51.53, 
         58.44, 45.39, 37.15, 28.05, 27.39, 26.82, 34.61, 31.05, 49.03, 51.79, 79.63, 130.68, 
         156.68, 119.39, 69, 45.43, 36.11, 33.1, 33.9, 32.4, 32.25, 47.55, 68.08, 98.85, 
         78.32, 62.22, 44.85, 31.58, 29.32, 27.15, 26.65, 25.24, 30.45, 43.39, 50.17, 88.16, 
         98, 68.32, 55.06, 35.28, 34, 32.43, 29.03, 23.68, 23.73, 24.98, 40.28, 60.1, 55.24, 
         53.04, 34.5, 27.27, 27, 34.08, 25.83, 26.73, 28.63, 37.93, 49.98, 95.06, 92.85, 
         63.44, 54.62, 32.62, 28.91, 26.12, 28.41, 25.79, 27.21, 38.84, 72.92, 92.98, 87.7, 
         66.72, 48.38, 36.31, 29.72, 25.01, 25.15, 24.8, 22.56, 27.19, 45.14, 49.81, 61.59,
         56.25, 35.97, 23.58, 22.83, 20.82, 17.16, 20.44, 22.51, 28.58, 47.17, 48.19, 49.5, 
         48.32, 32.43, 26.21, 22.61, 21.06, 19.64, 19.21, 20.05, 19.98, 33.69, 36.74, 41.86, 
         39.67, 31.89, 24.22, 24.77, 23.32, 23.38, 21.03, 25.59, 27.83, 42.69, 44.23, 57.09, 
         44.12, 37.79, 23.94, 26.67, 24.56, 20.49, 21.41, 25.14, 32.87, 56.61, 63.28, 69.32, 
         54.11, 36.99, 29.29, 29.34, 26.49, 26.07, 24.35, 27.12, 34.09, 44.23, 57.17, 73.7, 
         49.37, 40.08, 26.57, 25.89, 23.89, 23.66, 23.35, 26.51, 33.06, 54.13, 75.73, 57.07, 
         43.91, 33.23, 25.83, 24.14, 26.8, 23.08, 21.83, 25.34, 29.17, 47.65, 68.89, 49.86, 
         39.97, 36.04, 27.08, 22.77, 22.17, 20.39, 18.92, 20.51, 19.76, 24.82, 30.41, 39.41, 
         34.22, 29.98, 21.02, 19.76, 21.29, 22.8, 24.28, 30.11, 29.84, 49.04, 69.08, 83.77, 
         59.48, 38.73, 33.33, 25.29, 23.36, 21.16, 17.57, 18.87, 18.99, 23.14, 33.25, 40.72, 
         36.29, 25.86, 19.22, 18.47, 18.35, 18.01, 18.64, 26.27, 38.11, 46.6, 41.36, 32.16, 
         26.98, 22.78, 20.95, 22.06, 22.56, 21.81, 24.27, 36.45, 50.89, 85.68, 78.65, 64.6, 
         40.89, 32.57, 33.52, 31.16, 27.65, 30.43, 31.38, 42.25, 51.06, 95.54, 73.87, 59.43, 
         44.12, 29.82, 27.56, 27.08, 25.01, 36.9, 33.42, 44.05, 65.06, 91.67, 101.25, 71.29, 
         54.64, 33.9, 31.19, 31.35, 26.54, 26.96, 27.75, 30.91, 42.86, 49.9, 57.43, 52.39, 
         39.16, 37.1, 29.37, 27.77, 25.66, 24.82, 29.68, 32.42, 34.04, 43.42, 49.7, 39.72, 
         31.13, 22.77, 26.46, 29.01, 24.12, 34.24, 30.43, 41.73, 66.21, 96.4, 119.75, 93.03, 
         51.91, 40.59, 34.24, 37.53, 37.58, 33.54, 33.02, 41.69, 61.18, 85.78, 103.03, 59.1, 
         46.59, 37.95, 30.45, 29.46, 27.18, 25.26, 26.52, 29.74, 39.26, 45.58, 51.21, 48.02, 
         31.16, 26.21, 30.6, 26.5, 25.78, 25.68, 27.75, 35.1, 56.41, 65.78, 55.85, 42.79, 
         37.08, 34.66, 29.68, 19.74, 19.75, 19.88, 22.75, 24.83, 34.12, 50.5, 54.59, 39.74, 
         32.85, 27.21, 26.08, 24.82, 23.67, 23.22, 23.65, 25.46, 26.66, 29.21, 33.91, 31.65, 
         27.98, 21.47, 18.1, 17.32, 16.78, 15.75, 15.86, 19.8, 29.82, 38.21, 37.43, 30.52, 
         25.19, 17.53, 15.88, 17.73, 19.83, 35.86, 38.04, 50.61, 38.37, 25.58, 20.42, 18.43, 
         18.93, 16.46, 17.46, 18.23, 19.67, 31.11, 37.68, 35.82, 32.17, 20.06, 16.57, 18.91, 
         16.74, 16.16, 17.86, 17.15, 21.52, 31.73, 33.6, 37.36, 32.27, 26.99, 18.9, 16.54, 
         18.92, 16.89, 17.78, 18.59, 19.28, 30.47, 57.84, 61.37, 37.52, 24.32, 27.28, 26.77, 
         24.26, 22.94, 21.89, 23.9, 28.22, 39.27, 49.72, 50.01, 41.14, 30, 24.07, 23.92, 
         20.54, 18.61, 17.74, 19.82, 23.77, 32.17, 40.81, 35.09, 31.9, 22.15, 16.79, 15, 
         14.54, 15.52, 14.25, 13.96, 16.4, 28.13, 40.19, 36.19, 33.16, 21.8, 17.33, 17.48, 
         14.56, 13.81, 13.68, 13.7, 14.32, 20.33, 21.14, 26.52, 23.45, 19.32, 13.81, 13.08, 
         12.45, 13.72, 13.72, 14.97, 21.83, 35.08, 37.23, 35.64, 26.3, 23.39, 18.37, 17.01, 
         16.38, 12.27, 11.87, 13.01, 18.88, 27.35, 27.84, 25.99, 21.4, 15.68, 13.55, 13.86, 
         12.82, 11.37, 12.2, 12.82, 15.74, 31.26, 31.95, 27.58, 23.02, 19.35, 14.85, 13.99, 
         16.58, 11.32, 13.1, 19.09, 25.72, 40.09, 73.78, 75.82, 55.02, 24.44, 19.5, 21.79, 
         18.7, 17.59, 21.11, 20.9, 30.82, 45.72, 54.72, 52.81, 39.04)
  
  #Taken from require(evd); evd::dgev
  dgev <- function (x, loc = 0, scale = 1, shape = 0, log = FALSE) 
  {
    if (min(scale) <= 0) 
      stop("invalid scale")
    if (length(shape) != 1) 
      stop("invalid shape")
    x <- (x - loc)/scale
    if (shape == 0) 
      d <- log(1/scale) - x - exp(-x)
    else {
      nn <- length(x)
      xx <- 1 + shape * x
      xxpos <- xx[xx > 0 | is.na(xx)]
      scale <- rep(scale, length.out = nn)[xx > 0 | is.na(xx)]
      d <- numeric(nn)
      d[xx > 0 | is.na(xx)] <- log(1/scale) - xxpos^(-1/shape) - 
        (1/shape + 1) * log(xxpos)
      d[xx <= 0 & !is.na(xx)] <- -Inf
    }
    if (!log) 
      d <- exp(d)
    d
  }
  
  #m1 <- fgev(x) ; as.list(m1)
  mystart <- list(loc=24.6, scale=9.93, shape=0.3)
  data.size <- length(x)
  
  #check MLE no weight complete data with log argument
  mledist(x, "gev", optim.method="Nelder-Mead", start=mystart)[c("estimate", "value", "counts")]
  
  #check MLE with weight complete data with log argument
  mledist(x, "gev", optim.method="Nelder-Mead", start=mystart, weights=rep(1, data.size))[c("estimate", "value", "counts")]
  
  #check MLE no weight censored data
  mledist(data.frame(left=x, right=x), "gev", optim.method="Nelder-Mead", start=mystart)[c("estimate", "value", "counts")]
  
  #check MLE no weight censored data
  mledist(data.frame(left=x, right=x), "gev", optim.method="Nelder-Mead", 
          start=mystart, weights=rep(1, data.size))[c("estimate", "value", "counts")]
  
  #check manually with optim()
  fnobj_old <- function(par, obs, ddistnam)
  {
    T1 <- do.call(ddistnam, c(list(obs), as.list(par), log=TRUE) )
    -mean(T1)
  }
  
  optim(mystart, fnobj_old, method="Nelder-Mead", obs=x, ddistnam="dgev", control=list(trace=0))
  
  
  #check MLE no weight complete data without log argument
  dgev2 <- function(x, loc = 0, scale = 1, shape = 0)
    dgev(x)
  fnobj_old <- function(par, obs, ddistnam)
  {
    T1 <- do.call(ddistnam, c(list(obs), as.list(par)) )
    -mean(log(T1))
  }
  
  mledist(x, "gev2", optim.method="Nelder-Mead", start=mystart)[c("estimate", "value", "counts")]
  optim(mystart, fnobj_old, method="Nelder-Mead", obs=x, ddistnam="dgev2", control=list(trace=0))
}