R/fmaxlo.R

Defines functions `fmaxlo1` `fmaxlo`

##=======================================================================
## AUTHOR: Yves Deville
##
## Find ML estimate of a two-parameters 'maxlo' distribution using a
## sample 'x'. The likelihood is concentrated with respect to the shape
## parameter
##
## This distribution is a special case of a beta distribution on an
## interval (0, beta)  where 'beta' is considered as an unknown parameter.
##
## It is a reparametrisation of a GPD with location mu = 0 and negative
## shape xi < 0
##
##=======================================================================

`fmaxlo` <- function(x,
                     shapeMin = 1.25,
                     info.observed = TRUE,
                     plot = FALSE,
                     scaleData = TRUE,
                     cov = TRUE) {
    
    Cvg <- TRUE
    if (any(x <= 0)) stop("all elements in 'x' must be > 0")
    if (shapeMin <= 0) stop("'shapeMin' must be positive")
    if (shapeMin < 1) warning("'shapeMin <= 1': non-identifiable model")
    
    parnames <- c("shape", "scale")
    n <- length(x)
    M1 <- mean(x)
    
    if (scaleData) {
        ## xUnscaled <- x    ## useful ???
        x <- x / M1
        CV <- sqrt(1.0 - 1.0 / n) * sd(x)
        cLogLik <- - n * log(M1)
        trans <- diag(c(1.0, 1.0 / M1))
        colnames(trans) <- rownames(trans) <- parnames
    } else {
        CV <- sqrt(1.0 - 1.0 / n) * sd(x) / M1
    }
    
    if (CV > 1.00) stop("CV > 1. Estimation impossible for \"maxlo\"")
    ## if (CV > 0.999) warning("large CV value: estimation may not converge")

    ## compute an upper bound for beta. M1 <- mean(x) was done before
    M2 <- mean(x^2)
    M3 <- mean(x^3)
    
    if (scaleData) {
        betaRoots <- polyroot(c(-12 * M3, 10 * M3 - 6 * M2, 3 * (M2 - 2.0)))
    } else {
        betaRoots <- polyroot(c(-12 * M1 * M3, 10 * M3 - 6 * M1 * M2,
                                3 * (M2 - 2 * M1^2)))
    }
    
    betaUpper <- max(Re(betaRoots))
    mind <- max(x)
    if (betaUpper < 2.0 * mind) betaUpper <- 2.0 * mind
    
    interv <- c(mind * (1.0 + 1e-6), betaUpper)

    ##=======================================================================
    ## o 'logLc' will be maximised. THis is the logLik concentrated with the
    ##    constraint shape >= shapeMin, sot its derivative is no longer
    ##    given by the formula for unconstrained maximisation!
    ## 
    ## o 'dlogLc' WAS used to find a zero and is used to check the
    ##    bounds in 'interval'
    ## 
    ## o other functions are needed to compute the covariance (if needed)
    ##=======================================================================
    
    ## dlogLc <- function (beta) {
    ##     xmod <- x / beta
    ##     R <- -mean(log(1.0 - xmod))
    ##     dR <- -mean( xmod / (1.0 - xmod) ) / beta
    ##     -n * ( (1.0 - R) * dR / R + 1.0 / beta ) 
    ## }
    
    logLc <- function (beta) {
        xmod <- x / beta
        R <- -mean(log(1.0 - xmod))
        rho <- 1.0 / R
        if (rho < shapeMin) {
            res <- n * (log(shapeMin) - log(beta) - (shapeMin - 1.0) * R) 
        } else {
            res <- -n * (log(R) + log(beta) - R + 1.0)
        }
        res
    }

    if (plot) cov <- TRUE
        
    if (cov) {
        
        logL <- function (parm) {
            rho <- parm[1L]
            beta <- parm[2L]
            xmod <- x / beta
            R <- -mean(log(1.0 - xmod))
            n * (log(rho) - log(beta) - (rho - 1.0) * R)
        }
        
        ## 'beta' is scaled, as is 'x"
        log2L <- function (alpha, beta) {
            xmod <- x / beta
            xmod1 <- 1.0 - xmod
            RL <- mean(-log(xmod1))
            R1 <- mean(1.0 / xmod1)
            R2 <- mean(1.0 / xmod1 / xmod1)
            
            alpha1 <- alpha - 1.0
            
            logL <- n * (log(alpha / beta)  - alpha1 * RL) 
            
            dlogL <- c(shape = n * (1 / alpha - RL),
                       scale = n * (-1 + alpha1 * (R1 - 1.0)) / beta)
            
            d2logL <- array(0, dim = c(2, 2), dimnames = list(parnames, parnames))
            d2logL["shape", "shape"] <- -n / alpha / alpha
            d2logL["shape", "scale"] <- n * (R1 - 1) / beta 
            d2logL["scale", "shape"] <- d2logL["shape", "scale"]
            d2logL["scale", "scale"] <-
                n * (alpha - alpha1 * R2 ) / beta / beta 
            
            if (scaleData) {
                logL <- logL + cLogLik
                dlogL <- trans %*% dlogL 
                d2logL <- trans %*% d2logL %*% trans
                ## added on 2015-12-17 Yves
                rownames(dlogL) <- parnames
                dimnames(d2logL) <- list(parnames, parnames)
            } 
            
            list(logL = logL,
                 dlogL = dlogL,
                 d2logL = d2logL)    
        }
        
    } 
 
    ## ## check the bounds
    ## checks <- unlist(sapply(interv, dlogLc))
    ##
    ## if ( (checks[1] < 0) || (checks[2] > 0) ) {
    ##     cat("fmaxlo bounds\n"); print(checks)
    ##     warning("no interval found to maximise loglik")
    ##     Cvg <- FALSE
    ## }

    res <- optimize(f = logLc, interval = interv, maximum = TRUE)
    
    ## caution beta.hatS is for salaed data, and unscaled 'x' is lost 
    beta.hatS <- res$maximum
    alpha.hat <- -1 / mean(log(1 - x / beta.hatS))
    constr <- FALSE
    if (alpha.hat <= shapeMin) {
        constr <- TRUE
        alpha.hat <- shapeMin
    }
    
    if (scaleData) beta.hat <- M1 * beta.hatS
    else beta.hat <- beta.hatS
        
    if (!cov) {
        loglik <- res$objective 
        if (scaleData) loglik <- loglik + cLogLik
        return(list(estimate = c(shape = alpha.hat, scale = beta.hat),
                    CV = CV,
                    loglik = loglik,
                    cvg = Cvg))
    }
    
    res2 <- log2L(alpha = alpha.hat, beta = beta.hatS)
    
    if (alpha.hat > 2 && !constr) {
        
        if (info.observed) {
            info <- -res2$d2logL
            vcov <- solve(info)
            ## changed 2015-12-17: the test was wrong!!!
            if (inherits(cov, "try-error")) {
                vcov <- NULL
                sds <- NULL
                warning("hessian could not be inverted")
            } else {
                sds <- sqrt(diag(vcov))
            }
        } else {
            a1 <- alpha.hat - 1
            a2 <- alpha.hat - 2
            if (alpha.hat <= shapeMin + 1e-4) {
                warning("'shape' is near the bound 'shapeMin'. ",
                        "Derivatives may be missleading")
            }
            i11 <- 1 / alpha.hat / alpha.hat
            i12 <- -1 / beta.hat / a1
            i22 <- alpha.hat / a2 / beta.hat / beta.hat
            info <- matrix(n * c(i11, i12, i12, i22), nrow = 2L, ncol = 2L)
            colnames(info) <- rownames(info) <- parnames
            c11 <- alpha.hat^2 * a1^2
            c12 <- alpha.hat * a1 * a2 * beta.hat
            c22 <- a1^2 * a2 * beta.hat^2 / alpha.hat
            vcov <- matrix(c(c11, c12, c12, c22) / n, nrow = 2L, ncol = 2L)
            colnames(vcov) <- rownames(vcov) <- parnames
            sds <- sqrt(diag(vcov))
        }
        
    } else {
        ## return matrix or vector of NA
        if (constr){
            warning("'shape' is at the bound given in 'shapeMin'. ML inference results not suitable")
         } else {
             warning("'shape' is <= 2. ML inference results not suitable")
         }
        info <- matrix(NA, nrow = 2L, ncol = 2L)
        colnames(info) <- rownames(info) <- parnames
        ## print(cbind(info, -res2$d2logL)) ## for checks
        vcov <- info
        sds <- rep(NA, 2L)
        names(sds) <- parnames
    }
    
    if (plot) {
        
        dlogLc <- function (beta) {
            xmod <- x / beta
            R <- -mean(log(1.0 - xmod))
            dR <- -mean( xmod / (1.0 - xmod) ) / beta
            -n * ( (1.0 - R) * dR / R + 1.0 / beta ) 
        }
        
        if (scaleData) beta.sol <- beta.hat / M1
        else beta.sol <- beta.hat
        
        ## prepare grid and compute the limit
        lcInf <- -n * (1 + log(mean(x)))
        betas <- seq(from = interv[1], to = interv[2], length = 200)   
        fs <- sapply(betas, logLc)
        dfs <- sapply(betas, dlogLc)
        
        ind <- 1L:length(betas)
        ind <- dfs < 20
        Stext <- ifelse(scaleData, "(scaled data)", "")
        
        ## now plot logLik derivative and logLik
        opar <- par(mfrow = c(2L, 1L))
        par(mar = c(0, 5, 5, 5))
        
        ## First plot: logLik derivative. 
        plot(betas[ind], dfs[ind],
             type = "n",
             main = sprintf("'Maxlo' concentrated log-lik CV = %4.2f %s", CV, Stext),
             xlab = " ", ylab = "dlogL UNCONSTR.",
             xaxt = "n", yaxt = "n",
             xlim = interv)
        axis(side = 4)
        abline(h = 0)
        abline(v = beta.sol, col = "orangered")
        abline(v = interv, col = "darkcyan", lwd = 2)
        abline(h = 0, col = "gray")
        lines(betas[ind], dfs[ind],
              type = "l", lty = "solid", col = "red3")
    
        par(mar = c(5, 5, 0, 5))
        
        ## Second plot: logLik function 
        plot(betas[ind], fs[ind], type = "l",
             lty = "solid", col = "red3",
             xlab = "beta (scale param.)", ylab = "logL",
             xlim = interv, ylim = range(fs[ind], lcInf))
        abline(h = lcInf, col = "orchid")
        mtext(text = "lim.", side = 4, at = lcInf,
              col = "orchid")
        abline(v = interv, col = "darkcyan", lwd = 2)
        abline(v = beta.sol, col = "orangered")
        mtext(text = "betaHat", col = "orangered",
              side = 1, at = beta.sol, line = 0.5)

        par(opar)
        
    }
    
    list(estimate = c(shape = alpha.hat, scale = beta.hat),
         sd = sds,
         CV = CV,
         loglik = res2$logL,
         dloglik = res2$dlogL,
         cov = vcov,
         info = info,
         cvg = Cvg)
    
}

##=======================================================================
## Author: Yves Deville
##
## Find ML estimate of a two parameter 'maxlo' distribution WITH KNOWN
## SHAPE.
##
## NB. When the SCALE 'beta' is known, -log(1 - x / beta) is exponential
## with rate 'alpha'.
##
## TODO: add the 'scaleData' argument, find bounds for the optim...
##
##=======================================================================

`fmaxlo1` <- function(x,
                      shape = 1.5,
                      info.observed = TRUE,
                      scaleData = TRUE,
                      cov = TRUE,
                      plot = FALSE) {

    Cvg <- TRUE
    
    if (any(x <= 0)) stop("all elements in 'x' must be >0")  
    if (shape < 1.0) stop("'shape <= 1': non-identifiable model")
    
    n <- length(x)
    parnames <- c("shape", "scale")
    M1 <- mean(x)

    if (scaleData) {
        x <- x / M1
        CV <- sqrt(1.0 - 1.0 / n) * sd(x)
        cLogLik <- - n * log(M1)
        trans <- 1.0 / M1
    } else {
        CV <- sqrt(1.0 - 1.0 / n) * sd(x) / M1
    }
    
    alpha <- unname(shape)
    alpha1 <- alpha - 1.0
    
    ## could be used to find a zero
    dlogL1 <- function (beta) {
        xmod <- x / beta
        R <- -mean(log(1.0 - xmod))
        S1 <- mean(xmod / (1.0 - xmod))
        n * (-1.0  + alpha1 * S1) / beta 
    }
    
    logL1 <- function (beta) {
        xmod <- x / beta
        R <- -mean(log(1.0 - xmod))
        n * (log(alpha / beta) - alpha1 * R)
    }

    if (plot) cov <- TRUE

    if (cov) {
        log2L1 <- function (beta) {
            
            xmod <- x / beta
            R <- -mean(log(1.0 - xmod))
            S1 <- mean(xmod / (1.0 - xmod))
            S2 <- mean(xmod / (1.0 - xmod) / (1.0 - xmod))
            
            logL <- n * (log(alpha) - log(beta)  - alpha1 * R)
            dlogL <- n * (-1.0 + alpha1 * S1) / beta
            d2logL <- n * (1.0 - alpha1 * S1 - alpha1 * S2) / beta / beta
            
            if (scaleData) {
                logL <- logL + cLogLik
                dlogL <- trans * dlogL 
                d2logL <- trans * d2logL * trans
            }
            
            list(logL = logL,
                 dlogL = dlogL,
                 d2logL = d2logL)
        }
    }
        
    ## find interval
    mind <- max(x)
    betaL <- mind * (1.0 + alpha1 / n)
    if (scaleData) betaU <- 2 * pmax(mind, alpha1)
    else betaU <- 2 * pmax(mind, alpha1 * M1)
    
    interv <- c(betaL, betaU)
    checks <- unlist(sapply(interv, dlogL1))
    if ( (checks[1L] < 0) ||  (checks[2L] > 0) ) {
        warning("no interval found to maximise loglik")
        Cvg <- FALSE
    }
    
    res <- optimize(f = logL1, interval = interv, maximum = TRUE,
                    tol = .Machine$double.eps^0.3)
    beta.hatS <- res$maximum
    alpha.hat <- alpha
    
    if (scaleData) beta.hat <- M1 * beta.hatS
    else beta.hat <- beta.hatS

    if (!cov) {
        loglik <- res$objective
        if (scaleData) loglik <- loglik + cLogLik
        return(list(estimate = c(scale = beta.hat),
                    CV = CV,
                    loglik = loglik,
                    cvg = Cvg))
    }

    res2 <- log2L1(beta = beta.hatS)
    loglik <- res2$logL
    
    if (alpha >= 2.0) {
        if (info.observed) {
            info <- -res2$d2logL
        } else {
            info <-  n * alpha / (alpha - 2.0) / beta.hat / beta.hat
        }
        cov <- solve(info)
        sds <- sqrt(diag(cov))
        
    } else {
        warning("'shape' is < 2 ML inference results not suitable")
        info <- NA
        cov <- NA
        sds <- NA
    }
    
    if (plot) {
        Stext <- ifelse(scaleData, "(scaled data)", "")
        betas <- seq(from = betaL, to = betaU, length = 200)
        fs <- sapply(betas, logL1)
        dfs <- c(NA, diff(fs) / diff(betas))
        ind <- 1L:length(betas)
        ind <- dfs < 20
        plot(betas[ind], dfs[ind], type = "l",
             lty = "dotted",
             main = sprintf("'Lomax' log-lik derivative %s", Stext),
             xlab = "beta (scale)", ylab = "dlogL")
        lines(betas[ind], sapply(betas[ind], dlogL1), col = "orangered")
        abline(v = mind, col = "purple")
        abline(h = 0, v = beta.hatS, col = "SpringGreen3")
    }
    
    list(estimate = c(scale = beta.hat),
         sd = sds,
         CV = CV,
         loglik = loglik,
         dloglik = res2$dlogL,
         cov = cov,
         info = info,
         cvg = Cvg)
    
}

Try the Renext package in your browser

Any scripts or data that you put into this service are public.

Renext documentation built on Aug. 30, 2023, 1:06 a.m.