Code to estimate GEV distribution which allows for one or more parameters to be fixed (VERY PRELIMINARY)

gevfp.fit.R

gevfp.fit <- function (xdat, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL, 
                        mulink = identity, siglink = identity, shlink = identity, 
                        muinit = NULL, siginit = NULL, shinit = NULL, show = TRUE, 
                        method = "Nelder-Mead", maxit = 10000, 
                        fixedPars = list(mu = NULL, sig = NULL, sh = NULL), ...) {
        z <- list()
        npmu <- length(mul) + 1
        npsc <- length(sigl) + 1
        npsh <- length(shl) + 1
        z$trans <- FALSE
        in2 <- sqrt(6 * var(xdat))/pi
        in1 <- mean(xdat) - 0.57722 * in2
        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(sigl)) {
          sigmat <- as.matrix(rep(1, length(xdat)))
          if (is.null(siginit)) 
            siginit <- in2
        }
        else {
          z$trans <- TRUE
          sigmat <- cbind(rep(1, length(xdat)), ydat[, sigl])
          if (is.null(siginit)) 
            siginit <- c(in2, rep(0, length(sigl)))
        }
        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)))
        }
        
        ### include in inital values only the parameters for whcih estiamtion is required 
        z$model <- list(mul, sigl, shl)
        z$link <- deparse(substitute(c(mulink, siglink, shlink)))
        if(is.null(fixedPars[["mu"]])) muinit <- muinit
        if(!is.null(fixedPars[["mu"]])) muinit <- NULL
        if(is.null(fixedPars[["sig"]])) siginit <- siginit
        if(!is.null(fixedPars[["sig"]])) siginit <-NULL
        if(is.null(fixedPars[["sh"]])) shinit <- shinit
        if(!is.null(fixedPars[["sh"]])) shinit <- NULL
        init <- c(muinit, siginit, shinit)
        gev.lik <- function(a) {
          if(is.null(fixedPars[["mu"]])) mu <- mulink(mumat %*% (a[1:npmu])) 
          if(!is.null(fixedPars[["mu"]])) mu <- fixedPars[["mu"]]
          if(is.null(fixedPars[["sig"]])) sc <- siglink(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
          if(!is.null(fixedPars[["sig"]])) sc <-fixedPars[["sig"]]
          if(is.null(fixedPars[["sh"]])) xi <- shlink(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
          if(!is.null(fixedPars[["sh"]])) xi <- fixedPars[["sh"]]
          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
        if(is.null(fixedPars[["mu"]])) mu <- mulink(mumat %*% (x$par[1:npmu])) 
        if(!is.null(fixedPars[["mu"]])) mu <- fixedPars[["mu"]]
        if(is.null(fixedPars[["sig"]])) sc <- siglink(sigmat %*% (x$par[seq(npmu + 1, length = npsc)]))
        if(!is.null(fixedPars[["sig"]])) sc <-fixedPars[["sig"]]
        if(is.null(fixedPars[["sh"]])) xi <- shlink(shmat %*% (x$par[seq(npmu + npsc + 1, length = npsh)]))
        if(!is.null(fixedPars[["sh"]])) xi <- fixedPars[["sh"]]
        # mu <- mulink(mumat %*% (x$par[1:npmu]))
        # sc <- siglink(sigmat %*% (x$par[seq(npmu + 1, length = npsc)]))
        # xi <- shlink(shmat %*% (x$par[seq(npmu + npsc + 1, length = npsh)]))
        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, sc, 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)
}