Code to estimate GEV distribution based on coefficient of variation tau. See also the gevcvd.fit in https://github.com/ilapros/ilaprosUtils.

gevcv.fit.R

gevcv.fit <- function (xdat, ydat = NULL, mul = NULL, taul = NULL, shl = NULL, 
          mulink = identity, taulink = identity, shlink = identity, 
          muinit = NULL, tauinit = NULL, shinit = NULL, show = TRUE, 
          method = "Nelder-Mead", maxit = 10000, ...) {
  z <- list()
  npmu <- length(mul) + 1
  npcv <- length(taul) + 1
  npsh <- length(shl) + 1
  z$trans <- FALSE
  in2 <- sqrt(6 * var(xdat))/pi
  in1 <- mean(xdat) - 0.57722 * in2
  ins <- in2/in1
  if (is.null(mul)) {
    mumat <- as.matrix(rep(1, length(xdat)))
    if (is.null(muinit)) 
      muinit <- in1
  }
  else {
    z$trans <- TRUE
    mumat <- cbind(rep(1, length(xdat)), ydat[, mul])
    if (is.null(muinit)) 
      muinit <- c(in1, rep(0, length(mul)))
  }
  if (is.null(taul)) {
    taumat <- as.matrix(rep(1, length(xdat)))
    if (is.null(tauinit)) 
      tauinit <- in2
  }
  else {
    z$trans <- TRUE
    taumat <- cbind(rep(1, length(xdat)), ydat[, taul])
    if (is.null(tauinit)) 
      tauinit <- c(in2, rep(0, length(taul)))
  }
  if (is.null(shl)) {
    shmat <- as.matrix(rep(1, length(xdat)))
    if (is.null(shinit)) 
      shinit <- 0.1
  }
  else {
    z$trans <- TRUE
    shmat <- cbind(rep(1, length(xdat)), ydat[, shl])
    if (is.null(shinit)) 
      shinit <- c(0.1, rep(0, length(shl)))
  }
  z$model <- list(mul, taul, shl)
  z$link <- deparse(substitute(c(mulink, taulink, shlink)))
  init <- c(muinit, tauinit, shinit)
  gev.lik <- function(a) {
    mu <- mulink(mumat %*% (a[1:npmu]))
    cv <- taulink(taumat %*% (a[seq(npmu + 1, length = npcv)]))
    xi <- shlink(shmat %*% (a[seq(npmu + npcv + 1, length = npsh)]))
    sc <- cv*mu
    y <- (xdat - mu)/sc
    y <- 1 + xi * y
    if (any(y <= 0) || any(sc <= 0)) 
      return(10^6)
    sum(log(sc)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 
                                                    1))
  }
  x <- optim(init, gev.lik, hessian = TRUE, method = method, 
             control = list(maxit = maxit, ...))
  z$conv <- x$convergence
  mu <- mulink(mumat %*% (x$par[1:npmu]))
  cv <- taulink(taumat %*% (x$par[seq(npmu + 1, length = npcv)]))
  xi <- shlink(shmat %*% (x$par[seq(npmu + npcv + 1, length = npsh)]))
  sc <- mu * cv
  z$nllh <- x$value
  z$data <- xdat
  if (z$trans) {
    z$data <- -log(as.vector((1 + (xi * (xdat - mu))/sc)^(-1/xi)))
  }
  z$mle <- x$par
  z$cov <- solve(x$hessian)
  z$se <- sqrt(diag(z$cov))
  z$vals <- cbind(mu, cv, xi)
  if (show) {
    if (z$trans) 
      print(z[c(2, 3, 4)])
    else print(z[4])
    if (!z$conv) 
      print(z[c(5, 7, 9)])
  }
  class(z) <- "gev.fit"
  invisible(z)
}