R/misc.R

Defines functions strip.fevd.bayesian strip.fevd.lmoments strip.fevd.mle strip.fevd strip blockmaxxer.matrix blockmaxxer.data.frame blockmaxxer.vector blockmaxxer.fevd blockmaxxer revtrans.evd trans.fevd datagrabber.fevd datagrabber.declustered datagrabber.extremalindex trans

Documented in blockmaxxer blockmaxxer.data.frame blockmaxxer.fevd blockmaxxer.matrix blockmaxxer.vector datagrabber.declustered datagrabber.extremalindex datagrabber.fevd revtrans.evd strip strip.fevd strip.fevd.bayesian strip.fevd.lmoments strip.fevd.mle trans trans.fevd

trans <- function(object, ...) {
    UseMethod("trans", object)
} # end of 'trans' function.

datagrabber.extremalindex <- function(x, ...) {
    a <- attributes(x)
    # if(length(a$data.name) == 1) res <- c(get(a$data.name, ...))
    if(length(a$data.name) == 1) res <- c(a$data)
    else {
	look <- try(get(a$data.name[1], ...), silent = TRUE)
	if( is(look, "try-error" ) ) look <- a$data
	nm <- colnames(look)
	if(!is.null(nm) & is.element(a$data.name[2], nm)) res <- look[,a$data.name[2]]
	else res <- look[,as.numeric(a$data.name[2])]
    }
    res <- a$na.action(res)
    return(res)
} # end of 'datagrabber.extremalindex' function.

datagrabber.declustered <- function(x, ...) {
    a <- attributes(x)

    if(length(a$data.name) == 1) {
	res <- c(try(get(a$data.name, ...), silent=TRUE))
	if( is(res, "try-error" ) ) {
	    # res <- NULL
	    res <- a$data
	    return(res)
	}
    } else {
	look <- try(get(a$data.name[1], ...), silent=TRUE)
	if( is(look, "try-error" ) ) {
	    # res <- NULL
	    res <- a$data
	    return(res)
	}
	nm <- colnames(look)
	if(!is.null(nm)) {
	    if(!is.element(a$data.name[2], nm)) dc <- as.numeric(a$data.name[2])
	    else dc <- a$data.name[2]

	    if(!is.element(a$data.name[3], nm)) dc <- c(dc, as.numeric(a$data.name[3]))
            else dc <- c(dc, a$data.name[3])
	    
	} else dc <- as.numeric(a$data.name[2:3])
	res <- cbind(c(look[,dc[1]]), c(look[,dc[2]]))
    }
    return(res)
} # end of 'datagrabber.declustered' function.

datagrabber.fevd <- function(x, response=TRUE, cov.data=TRUE, ...) {
    # Get response and any covariate data sets.
    in.data <- x$in.data
    if(response) {
	# if(!in.data) {
         #    if(!is.null(x$data.pointer)) y <- get(x$data.pointer, ...)
          #   else y <- x$x
	# } else y <- model.response(model.frame(x$x.fun, data=get(x$data.name[2], ...)))
	y <- x$x
    } else y <- NULL

    if(cov.data) {
        # if(!is.null(x$cov.data)) cdata <- x$cov.data
        # else if(x$data.name[2] == "") cdata <- NULL
        # else cdata <- get(x$cov.pointer, ...)
	cdata <- x$cov.data
    } else cdata <- NULL

    if(response & cov.data) out <- cbind(y, cdata)
    else if(response) out <- y
    else out <- cdata
    if(!is.null(out)) out <- do.call(x$na.action, list(out))
    return(out)
} # end of 'datagrabber.fevd' function.

trans.fevd <- function(object, ..., burn.in=499, return.all = FALSE) {

    x <- object
    n <- x$n

    meth <- tolower(x$method)

    y <- c(datagrabber(x)[,1])

    type <- tolower(x$type)
    des <- setup.design(x)

    if(is.element(meth, c("mle","gmle"))) p <- x$results$par
    else if(meth == "bayesian") {

	p2 <- x$results
	np <- dim(p2)[2] - 1
	if(np == 1) p2 <- matrix(p2[,1], ncol=1)
	else p2 <- p2[,1:np]
	if(burn.in != 0) p <- colMeans(p2[-(1:burn.in),])
	else p <- colMeans(p2)

    }

    pnames <- names(p)

    # Calculate the location parameters for each variable.
    if(is.element(type, c("gev", "pp", "gumbel", "weibull", "frechet"))) {
	X.loc <- des$X.loc
	nloc <- ncol(X.loc)
	loc <- rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc)
    } else nloc <- 0

    # Scale
    X.sc <- des$X.sc
    nsc <- ncol(X.sc)
    scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc)
    if(x$par.models$log.scale) scale <- exp(scale)

    # Shape
    if(!is.element(type, c("gumbel","exponential"))) {
	X.sh <- des$X.sh
 	nsh <- ncol(X.sh)
	shape <- rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], n, nsh, byrow=TRUE)*X.sh)
    } else shape <- 0

    if(is.element(type, c("pp", "gp", "exponential", "beta", "pareto"))) {
	u <- x$threshold
	eid <- y > u
	if(!return.all) {
	    scale <- scale[eid]
	    y <- y[eid]
	    if(length(u) > 1) u <- u[eid]
	    shape <- shape[eid]
	    if(type=="pp") scale <- scale + shape*(u - loc[eid])
	} else if(type=="pp") scale <- scale + shape*(u - loc)
    }

    if(is.element(type, c("gev", "weibull", "frechet"))) z <- -log(as.vector((1 + (shape * (y - loc))/scale)^(-1/shape)))
    else if(type=="gumbel") z <- as.vector((y - loc)/scale)
    else if(is.element(type, c("gp","beta","pareto"))) z <- -log(as.vector((1 + (shape * (y - u))/scale)^(-1/shape)))
    else if(type=="pp") z <- as.vector((1 + (shape * (y - u)/scale))^(-1/shape))
    else if(type=="exponential") z <- as.vector((y - u)/scale)
    return(z)

} # end of 'trans.fevd' function.

revtrans.evd <- function(z, threshold=NULL, location=NULL, scale, shape=NULL, type=c("GEV","GP","PP","Gumbel","Weibull","Frechet","Exponential","Beta","Pareto")) {
    type <- match.arg(type)
    if(is.element(type,c("Weibull","Frechet"))) type <- "GEV"
    else if(is.element(type,c("Beta","Pareto"))) type <- "GP"
    type <- tolower(type)
    if(is.element(type, c("gev","pp","gp"))) {
	zid <- shape==0
	if(any(zid)) shape[zid] <- 1e-10
    }
    if(type=="gumbel") out <- scale * z + location
    else if(type=="exponential") out <- scale * z + threshold
    else if(type=="gev") out <- location + (scale/shape) * (exp(shape * z) - 1)
    else if(type=="gp") out <- threshold + (scale/shape) * (exp(shape * z) - 1)
    else if(type=="pp") out <- threshold - (1 + shape * log(z)) * (scale + shape * (threshold - location))/shape
    return(out)
} # end of 'revtrans.evd' function.

blockmaxxer <- function(x, ...) {
    UseMethod("blockmaxxer", x)
} # end of 'blockmaxxer' function.

blockmaxxer.fevd <- function(x, ...) {

    if(x$type != "PP") stop("blockmaxxer: method for PP models only.")

    y <- datagrabber(x)

    blocks <- rep(1:x$span, each=x$npy)
    n2 <- length(blocks)
    if(n2 < x$n) blocks <- c(blocks, rep(blocks[n2], x$n - n2))
    else if(n2 > x$n) blocks <- blocks[1:x$n]

    res <- blockmaxxer(y, blocks = blocks)

    return(res)

} # end of 'blockmaxxer.fevd' function.

blockmaxxer.vector <- function(x, ..., blocks = NULL, blen = NULL, span = NULL) {

    if(is.null(blocks) && (is.null(blen) || is.null(span))) stop("blockmaxxer: must supply one of blocks or blen and span.")

    n <- length(x)
    if(is.null(blocks)) {
        blocks <- rep(1:span, each = blen)
        nb <- length(blocks)
        if(nb < n) blocks <- c(blocks, rep(blocks[nb], n - nb))
        else if(nb > n) blocks <- blocks[1:n]
    }

    res <- c(aggregate(x, by = list(blocks), max, ...)$x)
    return(res)
} # end of 'blockmaxxer.vector' function.

blockmaxxer.data.frame <- function(x, ..., which = 1, blocks = NULL, blen = NULL, span = NULL) {

    if(is.null(blocks) && (is.null(blen) || is.null(span))) stop("blockmaxxer: must supply one of blocks or blen and span.")

    xd <- dim(x)

    if(is.null(blocks)) {
	blocks <- rep(1:span, each = blen)
	nb <- length(blocks)
	if(nb < xd[1]) blocks <- c(blocks, rep(blocks[nb], xd[1] - nb))
	else if(nb > xd[1]) blocks <- blocks[1:xd[1]]
    }

    res <- c(aggregate(x[[which]], by = list(blocks), max, ...)$x)

    if(xd[2] > 1) {

        bl <- unique(blocks)
	n <- length(bl)

	o <- 1:length(x[[which]])
	ind <- numeric(0)

	tmp <- x[[which]]

	for(i in 1:n) {
	    id <- blocks == bl[i]
	    o2 <- o[id]
	    x2 <- tmp[id]
	    id2 <- x2 == max(x2, ...)
	    id2 <- id2

	    ind <- c(ind, min(o2[id2]))

	} # end of for 'i' loop.

	res <- x[ind,]
    }

    return(res)

} # end of 'blockmaxxer.data.frame' function.

blockmaxxer.matrix <- function(x, ..., which = 1, blocks = NULL, blen = NULL, span = NULL) {

    x <- as.data.frame(x)

    res <- blockmaxxer(x, ..., which = which, blocks = blocks, blen = blen, span = span)
    res <- as.matrix(res)

    return(res)

} # end of 'blockmaxxer.matrix' function.

strip <- function( x, use.names = TRUE, ... ) {

    UseMethod( "strip", x )

} # end of 'strip' function.

strip.fevd <- function( x, use.names = TRUE, ... ) {

    if (x$method == "GMLE") newcl <- "mle"
    else newcl <- tolower(x$method)

    # class(x) <- c( paste("fevd.", newcl, sep = ""), class( x ) )
    # UseMethod("strip", x)

    get( paste( "strip.fevd.", newcl, sep = "" ) )( x = x, use.names = use.names, ... )

} # end of 'strip.fevd' function.

strip.fevd.mle <- function( x, use.names = TRUE, ... ) {

    out <- x$results$par

    if( !use.names ) names( out ) <- NULL

    return( out )

} # end of 'strip.fevd.mle' function.

strip.fevd.lmoments <- function( x, use.names = TRUE, ... ) {

    out <- x$results

    if( !use.names ) names( out ) <- NULL

    return( out )

} # end of 'strip.fevd.lmoments' function.

strip.fevd.bayesian <- function( x, use.names = TRUE, burn.in = 499, FUN = "mean", ... ) {

    f <- match.fun(FUN)
    p <- x$results
    np <- dim(p)[2] - 1
    p <- p[, -(np + 1)]

    pnames <- colnames(p)

    if (is.element("log.scale", pnames)) {

        id <- pnames == "log.scale"
        p[, "log.scale"] <- exp(p[, "log.scale"])
        pnames[id] <- "scale"
        colnames(p) <- pnames

    }

    if (burn.in != 0) {

        n <- dim(p)[1]
        if (missing(burn.in)) 
            if (burn.in <= 2 * n - 1) 
                burn.in <- floor(n/4)
            else if (burn.in <= n - 1) 
                stop("distill: number of MCMC iterations too small compared to burn.in")
        p <- p[-(1:burn.in), ]

    }

    if( FUN == "mean" ) out <- colMeans(p, na.rm = TRUE)
    else if( FUN == "postmode" || FUN == "mode" ) out <- postmode( x, burn.in = burn.in )
    else out <- apply( p, 2, f, ... )

    if( use.names && is.null( names( out ) ) ) names( out ) <- pnames
    else if( !use.names && !is.null( names( out ) ) ) names( out ) <- NULL

    return( out )

} # end of 'strip.fevd.bayesian' function.

Try the extRemes package in your browser

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

extRemes documentation built on May 29, 2024, 5:27 a.m.