Nothing
R/fitfuns.R
In extRemes: Extreme Value Analysis
Defines functions
fpois.list
fpois.matrix
fpois.data.frame
fpois.default
fpois
fevdProposalDefault
fevdPriorDefault
will.accept
initializer.mle
initializer.moms
initializer.lmoments
initializer
shapePriorBetaOld
shapePriorBeta
oevdgen
oevd
setup.design.fevd
setup.design.default
setup.design
check.constant
plot.grlevdTracer
grlevdTracer
grlevd
print.fevd
summary.fevd.mle
summary.fevd.bayesian
summary.fevd.lmoments
summary.fevd
distill.fevd.mle
distill.fevd.bayesian
distill.fevd.lmoments
distill.fevd
parcov.fevd
levd
fevd
Documented in
check.constant
distill.fevd
distill.fevd.bayesian
distill.fevd.lmoments
distill.fevd.mle
fevd
fevdPriorDefault
fevdProposalDefault
fpois
fpois.data.frame
fpois.default
fpois.list
fpois.matrix
grlevd
grlevdTracer
initializer
initializer.lmoments
initializer.mle
initializer.moms
levd
oevd
oevdgen
parcov.fevd
plot.grlevdTracer
print.fevd
setup.design
setup.design.default
setup.design.fevd
shapePriorBeta
shapePriorBetaOld
summary.fevd
summary.fevd.bayesian
summary.fevd.lmoments
summary.fevd.mle
will.accept
fevd <- function(x, data, threshold=NULL, threshold.fun=~1, location.fun=~1, scale.fun=~1, shape.fun=~1, use.phi=FALSE,
type=c("GEV", "GP", "PP", "Gumbel", "Exponential"),
method=c("MLE", "GMLE", "Bayesian", "Lmoments"), initial=NULL, span, units=NULL, time.units="days", period.basis="year",
na.action=na.fail, optim.args=NULL, priorFun=NULL, priorParams=NULL, proposalFun=NULL, proposalParams=NULL,
iter=9999, weights=1, blocks=NULL, verbose=FALSE) {
if(verbose) begin.tiid <- Sys.time()
out <- list()
inout <- list() # object to return initial parameter values with likelihood.
out$call <- match.call()
if(!missing(data)) {
out$data.name <- c(deparse(substitute(x)), deparse(substitute(data)))
# cov.pointer <- as.character(substitute(data))
} else {
out$data.name <- c(deparse(substitute(x)), "")
# cov.pointer <- NULL
}
# data.pointer <- as.character(substitute(x))
# if( length(data.pointer) > 1) out$x <- x
# else out$data.pointer <- data.pointer
# if(length(cov.pointer) > 1) out$cov.data <- data
# else out$cov.pointer <- cov.pointer
type <- match.arg(type)
method <- match.arg(method)
out$weights <- weights
if(!missing(data)) {
if(is.element(out$data.name[1], colnames(data))) {
out$in.data <- TRUE
wc <- out$data.name[1] == colnames(data)
x <- c(data[,wc])
x.fun <- ifelse(out$data.name[1] == "substitute(x)", deparse(x), out$data.name[1])
x.fun <- formula(paste(x.fun, "~ 1"))
out$x.fun <- x.fun
} else if(is.formula(x)) {
out$in.data <- TRUE
x.fun <- x
out$x.fun <- x.fun
x <- model.response(model.frame(x.fun, data=data))
} else out$in.data <- FALSE
if(length(x) != nrow(data)) stop("fevd: data must have same number of rows as the length of x.")
if(!identical(weights, 1) && length(x) != length(weights)) # CJP
stop("fevd: weights should be the same length as x.")
out$missing.values <- is.na( x )
tmp <- cbind(x, data)
tmp <- na.action(tmp)
x <- tmp[,1]
data <- tmp[,-1, drop=FALSE]
} else {
out$missing.values <- is.na( x )
if(is.formula(x)) stop("fevd: Must provide data argument if you supply a formula to the x argument.")
x <- na.action(x)
out$in.data <- FALSE
}
if(!out$in.data) {
data.pointer <- as.character(substitute(x))
if(length(data.pointer) > 1) out$x <- x
else out$data.pointer <- data.pointer
}
if(!is.null(blocks)){ # CJP
# test that blocks is of correct type
if(type == "PP"){
if(!is.element('nBlocks', names(blocks))){
if(is.element('data', names(blocks))) {
blocks$nBlocks <- nrow(blocks$data)
} else stop("fevd: When supplying blocks, must provide 'blocks$nBlocks' if 'blocks$data' is not provided.")
}
if(!is.element('weights', names(blocks))) blocks$weights <- 1
if(!is.element('proportionMissing', names(blocks))) blocks$proportionMissing <- 0
if(!is.element('threshold', names(blocks)) && !is.null(threshold)) {
if(length(threshold) == 1) {
blocks$threshold <- threshold
} else {
stop("fevd: No blocks$threshold specified and threshold is not constant. User must supply the threshold for each block via blocks$threshold.")
}
}
} else {
warning("fevd: Blocks are used only for type 'PP'; ignoring blocks argument.")
blocks <- NULL # code will now operate as if user did not supply 'blocks'
}
}
out$x <- x
if(!missing(data)) out$cov.data <- data
if( method == "MLE" && !is.null(priorFun) ) method <- "GMLE"
else if( method == "GMLE" && is.null(priorFun) ) {
if(shape.fun != ~1) stop("fevd: must supply a prior function for GMLE method when shape parameter varies.")
if(is.element(type, c("GEV","GP","PP"))) {
priorFun <- "shapePriorBeta"
if(is.null(priorParams)) priorParams <- list(q=6, p=9)
} else {
warning("fevd: Using method MLE. No default for specified arguments.")
method <- "MLE"
}
} # end of if else method is GMLE but no prior function input stmts.
if(method=="GMLE") {
out$priorFun <- priorFun
out$priorParams <- priorParams
}
out$method <- method
method <- tolower(method)
out$type <- type
type <- tolower(type)
out$period.basis <- period.basis
out$optim.args <- optim.args
out$units <- units
if(method=="bayesian" && missing(use.phi)) {
use.phi <- TRUE
if(verbose) cat("\n", "Setting use.phi argument to TRUE for greater stability in estimation (default with Bayesian method). Use use.phi=FALSE if you prefer that.\n")
}
# if(is.element(type, c("gp","pp","exponential","beta","pareto"))) out$npy <- npy
out$par.models <- list(threshold=threshold.fun, location=location.fun, scale=scale.fun, log.scale=use.phi, shape=shape.fun,
term.names=list(threshold=all.vars(threshold.fun), location=all.vars(location.fun),
scale=all.vars(scale.fun), shape=all.vars(shape.fun)))
pars <- list()
if(is.element(type, c("gp","pp","exponential","beta","pareto"))) {
const.thresh <- check.constant(threshold.fun) & check.constant(threshold)
out$const.thresh <- const.thresh
}
if(is.element(type, c("gev","pp","gumbel","weibull","frechet"))) {
const.loc <- check.constant(location.fun)
out$const.loc <- const.loc
}
const.scale <- check.constant(scale.fun)
out$const.scale <- const.scale
const.shape <- check.constant(shape.fun)
out$const.shape <- const.shape
if(is.element(type, c("pp", "gp", "exponential", "beta", "pareto"))) {
if(missing(span)) {
if(is.null(blocks)) {
tiden <- attributes(x)$times
n <- length(x)
if(is.null(tiden)) {
tiden <- 1:n
start <- 1
end <- n
span <- end - start
} else {
start <- tiden[1]
end <- tiden[n]
span <- as.numeric(difftime(as.POSIXlt(tiden)[n], as.POSIXlt(tiden)[1], units=time.units))
}
} else {
span <- blocks$nBlocks
}
} # end of if 'missing span' stmts.
if(time.units=="days") npy <- 365.25
else if(time.units=="months") npy <- 12
else if(time.units=="years") npy <- 1
else if(time.units == "hours") npy <- 24 * 365.25
else if(time.units == "minutes") npy <- 60 * 24 * 365.25
else if(time.units == "seconds") npy <- 60 * 60 * 24 * 365.25
else {
tmp.units <- unlist(strsplit(time.units, split="/"))
if(length(tmp.units) != 2) stop("fevd: invalid time.units argument.")
numper <- as.numeric(tmp.units[1])
if(is.na(numper)) stop("fevd: invalid time.units argument.")
pertiid <- tmp.units[2]
if(!is.element(pertiid, c("day","month","year","hour","minute","second"))) stop("fevd: invalid time.units argument.")
if(pertiid=="year") npy <- numper
else if(pertiid=="month") npy <- numper*12
else if(pertiid=="day") npy <- numper*365.25
else if(pertiid == "hour") npy <- numper * 24 * 365.25
else if(pertiid == "minute") npy <- numper * 60 * 24 * 365.25
else if(pertiid == "second") npy <- numper * 60 * 60 * 24 * 365.25
} # end of if complicated time units stmts.
if(!is.null(blocks)) span <- span*npy # CJP - so is on the same scale as span w/o blocks
out$time.units <- time.units
out$span <- span/npy
out$npy <- npy
if(verbose)
cat("\n", "Data span ", span/npy, "years", "\n")
} else npy <- NULL # end of if 'type is of the POT variety' stmts.
# if(!missing(data)) {
# if(is.element(out$data.name[1], colnames(data))) {
# out$in.data <- TRUE
# wc <- out$data.name[1] == colnames(data)
# x <- c(data[,wc])
# x.fun <- ifelse(out$data.name[1] == "substitute(x)", deparse(x), out$data.name[1])
# x.fun <- formula(paste(x.fun, "~ 1"))
# out$x.fun <- x.fun
# } else if(is.formula(x)) {
# out$in.data <- TRUE
# x.fun <- x
# out$x.fun <- x.fun
# x <- model.response(model.frame(x.fun, data=data))
# } else out$in.data <- FALSE
# if(length(x) != nrow(data)) stop("fevd: data must have same number of rows as the length of x.")
# tmp <- cbind(x, data)
# tmp <- na.action(tmp)
# x <- tmp[,1]
# data <- tmp[,-1]
# } else {
# if(is.formula(x)) stop("fevd: Must provide data argument if you supply a formula to the x argument.")
# x <- na.action(x)
# out$in.data <- FALSE
# }
#
# if(!out$in.data) {
# data.pointer <- as.character(substitute(x))
# if(length(data.pointer) > 1) out$x <- x
# else out$data.pointer <- data.pointer
# }
n <- length(x)
out$n <- n
out$na.action <- deparse(substitute(na.action))
if(!is.null(initial)) {
if(!is.list(initial)) stop("fevd: initial must be NULL or a list object.")
find.init <- FALSE
if(is.null(initial$location) && is.element(type, c("gev","pp","gumbel","weibull","frechet"))) find.init <- TRUE
if(use.phi && is.null(initial$log.scale)) find.init <- TRUE
if(!use.phi && is.null(initial$scale)) find.init <- TRUE
if(!is.element(type, c("gumbel","exponential")) && is.null(initial$shape)) find.init <- TRUE
} else {
initial <- list()
find.init <- TRUE
}
if(method != "lmoments") {
if(verbose) cat("Setting up parameter model design matrices.\n")
designs <- list()
if(!missing(data)) {
if(is.element(type, c("gp", "pp", "exponential", "beta", "pareto"))) X.u <- setup.design(x=threshold.fun, data=data, n=n, dname="threshold.fun")
if(is.element(type, c("gev", "pp", "gumbel", "weibull", "frechet"))) {
X.loc <- setup.design(x=location.fun, data=data, n=n, const=const.loc, dname="location.fun")
designs$X.loc <- X.loc
}
X.sc <- setup.design(x=scale.fun, data=data, n=n, const=const.scale, dname="scale.fun")
designs$X.sc <- X.sc
if(!is.element(type,c("gumbel","exponential"))) {
X.sh <- setup.design(x=shape.fun, data=data, n=n, const=const.shape, dname="shape.fun")
designs$X.sh <- X.sh
}
} else {
if(is.element(type, c("gp", "pp", "exponential", "beta", "pareto"))) X.u <- setup.design(x=threshold.fun, n=n, dname="threshold.fun")
if(is.element(type, c("gev", "pp", "gumbel", "weibull", "frechet"))) {
X.loc <- setup.design(x=location.fun, n=n, const=const.loc, dname="location.fun")
designs$X.loc <- X.loc
}
X.sc <- setup.design(x=scale.fun, n=n, const=const.scale, dname="scale.fun")
designs$X.sc <- X.sc
if(!is.element(type,c("gumbel","exponential"))) {
X.sh <- setup.design(x=shape.fun, n=n, const=const.shape, dname="shape.fun")
designs$X.sh <- X.sh
}
} # end of if else missing data stmts.
if(!is.null(blocks)){ # CJP
blocks$designs <- list()
if(is.element('data', names(blocks))){
blocks$X.u <- setup.design(x = threshold.fun, data=blocks$data, n=blocks$nBlocks, dname="threshold.fun")
blocks$designs$X.loc <- setup.design(x=location.fun, data=blocks$data, n=blocks$nBlocks, const=const.loc, dname="location.fun")
blocks$designs$X.sc <- setup.design(x=scale.fun, data=blocks$data, n=blocks$nBlocks, const=const.scale, dname="scale.fun")
blocks$designs$X.sh <- setup.design(x=shape.fun, data=blocks$data, n=blocks$nBlocks, const=const.shape, dname="shape.fun")
} else {
blocks$X.u <- setup.design(x=threshold.fun, n=blocks$nBlocks, dname="threshold.fun")
blocks$designs$X.loc <- setup.design(x=location.fun, n=blocks$nBlocks, const=const.loc, dname="location.fun")
blocks$designs$X.sc <- setup.design(x=scale.fun, n=blocks$nBlocks, const=const.scale, dname="scale.fun")
blocks$designs$X.sh <- setup.design(x=shape.fun, n=blocks$nBlocks, const=const.shape, dname="shape.fun")
}
}
if(verbose) cat("Parameter model design matrices set up.\n")
} # end of if method is not Lmoments to set up design matrices stmts.
if(is.element(type, c("gp", "pp", "exponential", "beta", "pareto"))) {
if(method != "lmoments") threshold <- rowSums(matrix(threshold, n, ncol(X.u), byrow=TRUE) * X.u)
if(!is.null(blocks)) # CJP
blocks$threshold <- rowSums(matrix(blocks$threshold, blocks$nBlocks, ncol(blocks$X.u), byrow=TRUE) * blocks$X.u)
excess.id <- x > threshold
if(all(threshold==threshold[1])) out$threshold <- threshold[1]
else out$threshold <- threshold
out$rate <- mean(excess.id)
} # end of if '!missing(threshold)' stmts.
out$blocks <- blocks # CJP
if(method=="lmoments" || find.init) {
if(method=="lmoments") {
if(verbose) cat("Beginning estimation procedure.\n")
is.constant <- unlist(lapply(list(u=threshold, loc=location.fun, scale=scale.fun, sh=shape.fun), check.constant))
if(!all(is.constant)) warning("fevd: For method Lmoments, this function does not handle covariates in the parameters. Fitting w/o covariates.")
if(!is.element(type, c("gev","gp"))) stop("fevd: currently, Lmoments are only handled for estimation of GEV and GP distribution parameters.")
}
xtemp <- x
class(xtemp) <- "lmoments"
ipars1 <- try(initializer(xtemp, model=type, threshold=threshold, npy=npy, blocks=blocks), silent = TRUE)
if( !is( ipars1, "try-error" ) ) { # Eric -- 8/6/13
if(ipars1["scale"] <= 0) ipars1["scale"] <- 1e-8
if(method=="lmoments") {
out$results <- ipars1
class(out) <- "fevd"
return(out)
}
} else {
ipars1 <- NULL
if(method == "lmoments") stop("fevd: Sorry, could not find L-moments estimates.")
}
} # end of if 'Lmoments method or find them anyway as potential initial estimates' stmts.
if((method != "lmoments") && find.init) {
xtemp <- x
class(xtemp) <- "moms"
ipars2 <- try(initializer(xtemp, model=type, threshold=threshold, npy=npy, blocks=blocks), silent = TRUE)
if( !is(ipars2, "try-error" ) ) {
if(ipars2["scale"] <= 0) ipars2["scale"] <- 1e-8
} else ipars2 <- NULL
if(!is.element(type, c("pp", "gp","exponential","beta","pareto","gumbel"))) {
if(!is.null(ipars1)) testLmoments <- levd(x, location=ipars1["location"], scale=ipars1["scale"], shape=ipars1["shape"], type=out$type, npy=npy)
else testLmoments <- Inf
if(!is.null(ipars2)) testMoments <- levd(x, location=ipars2["location"], scale=ipars2["scale"], shape=ipars2["shape"], type=out$type, npy=npy)
else testMoments <- Inf
} else if(type=="pp") { # CJP
if(!is.null(ipars1)) { # Eric -- 8/6/13
if(!is.null(blocks)) {
blocks$location=ipars1["location"]; blocks$scale=ipars1["scale"]; blocks$shape=ipars1["shape"]
}
testLmoments <- levd(x, threshold=threshold, location=ipars1["location"], scale=ipars1["scale"], shape=ipars1["shape"], type=out$type, npy=npy, blocks=blocks)
} else testLmoments <- Inf
# end of if 'ipars1' not NULL stmts.
if(!is.null(ipars2)) { # Eric -- 8/6/13
if(!is.null(blocks)) {
blocks$location=ipars2["location"]; blocks$scale=ipars2["scale"]; blocks$shape=ipars2["shape"]
}
testMoments <- levd(x, threshold=threshold, location=ipars2["location"], scale=ipars2["scale"], shape=ipars2["shape"], type=out$type, npy=npy, blocks=blocks)
} else testMoments <- Inf
# end of if 'ipars2' not NULL stmts.
if(!is.null(blocks)) # CJP
blocks$location <- blocks$scale <- blocks$shape <- NULL
} else if(!is.element(type, c("gumbel","exponential"))) {
if(!is.null(ipars1)) testLmoments <- levd(x, threshold=threshold, scale=ipars1["scale"], shape=ipars1["shape"], type=out$type, npy=npy)
else testLmoments <- Inf
if(!is.null(ipars2)) testMoments <- levd(x, threshold=threshold, scale=ipars2["scale"], shape=ipars2["shape"], type=out$type, npy=npy)
else testMoments <- Inf
} else if(type=="gumbel") {
if(!is.null(ipars1)) testLmoments <- levd(x, location=ipars1["location"], scale=ipars1["scale"], type=out$type, npy=npy)
else testLmoments <- Inf
if(!is.null(ipars2)) testMoments <- levd(x, location=ipars2["location"], scale=ipars2["scale"], type=out$type, npy=npy)
else testMoments <- Inf
} else if(type=="exponential") {
if(!is.null(ipars1)) testLmoments <- levd(x, threshold=threshold, scale=ipars1["scale"], shape=0, type=out$type, npy=npy)
else testLmoments <- Inf
if(!is.null(ipars2)) testMoments <- levd(x, threshold=threshold, scale=ipars2["scale"], shape=0, type=out$type, npy=npy)
else testMoments <- Inf
}
if(is.finite(testLmoments) || is.finite(testMoments)) { # Eric -- 8/6/13
if(testLmoments < testMoments) {
if(is.null(initial$location) && !is.element(type, c("gp","exponential","beta","pareto"))) initial$location <- ipars1["location"]
if(is.null(initial$log.scale) && use.phi) initial$log.scale <- log(ipars1["scale"])
else if(is.null(initial$scale)) initial$scale <- ipars1["scale"]
if(!is.element(type, c("gumbel","exponential")) && is.null(initial$shape)) initial$shape <- ipars1["shape"]
if(verbose) cat("Using Lmoments estimates as initial estimates. Initial value =", testLmoments, "\n")
} else {
if(is.null(initial$location) && !is.element(type, c("gp","exponential","beta","pareto"))) initial$location <- ipars2["location"]
if(is.null(initial$log.scale) && use.phi) initial$log.scale <- log(ipars2["scale"])
else if(is.null(initial$scale)) initial$scale <- ipars2["scale"]
if(!is.element(type, c("gumbel","exponential")) && is.null(initial$shape)) initial$shape <- ipars2["shape"]
if(verbose) cat("Initial estimates found where necessary (not from Lmoments). Initial value =", testMoments, "\n")
}
} else { # Eric -- 8/6/13; if both L-moments and Moments fails, set initial values to mu = 0, sig = 1, xi = 0.01
if(is.null(initial$location) && !is.element(type, c("gp","exponential","beta","pareto"))) initial$location <- 0
if(is.null(initial$log.scale) && use.phi) initial$log.scale <- 0
else if(is.null(initial$scale)) initial$scale <- 1
if(!is.element(type, c("gumbel","exponential")) && is.null(initial$shape)) initial$shape <- 0.01
warning("fevd: L-moments and Moment initial parameter estimates could not be calculated. Using arbitrary starting values.")
}
inout <- list(Lmoments=list(pars=ipars1, likelihood=testLmoments), MOM=list(pars=ipars2, likelihood=testMoments))
} # end of if 'method is not lmoments and need to find initial estimates' stmts.
# Set up initial estimates.
if(!is.null(initial$location)) {
if(ncol(X.loc) != length(initial$location)) {
if((length(initial$location)==1) && ncol(X.loc) > 1) initial$location <- c(initial$location, rep(0, ncol(X.loc)-1))
else stop("fevd: initial parameter estimates must have length 1 or number of parameters present. Incorrect number for location parameter.")
}
if(length(initial$location)==1) names(initial$location) <- "location"
else names(initial$location) <- paste("mu", 0:(ncol(X.loc)-1), sep="")
}
if(use.phi && (ncol(X.sc) != length(initial$log.scale))) {
if((length(initial$log.scale)==1) && ncol(X.sc) > 1) initial$log.scale <- c(initial$log.scale, rep(0, ncol(X.sc)-1))
else stop("fevd: initial parameter estimates must have length 1 or number of parameters present. Incorrect number for log(scale) parameter.")
} else if(!use.phi && (ncol(X.sc) != length(initial$scale))) {
if((length(initial$scale)==1) && ncol(X.sc) > 1) initial$scale <- c(initial$scale, rep(0, ncol(X.sc)-1))
else stop("fevd: initial parameter estimates must have length 1 or number of parameters present. Incorrect number for scale parameter.")
}
if(use.phi) {
if(length(initial$log.scale)==1) names(initial$log.scale) <- "log.scale"
else names(initial$log.scale) <- paste("phi", 0:(ncol(X.sc)-1), sep="")
} else {
if(length(initial$scale)==1) names(initial$scale) <- "scale"
else names(initial$scale) <- paste("sigma", 0:(ncol(X.sc)-1), sep="")
}
if(!is.element(type, c("gumbel","exponential"))) {
if(ncol(X.sh) != length(initial$shape)) {
if((length(initial$shape)==1) && ncol(X.sh) > 1) initial$shape <- c(initial$shape, rep(0, ncol(X.sh)-1))
else stop("fevd: initial parameter estimates must have length 1 or number of parameters present. Incorrect number for shape parameter.")
}
if(length(initial$shape)==1) names(initial$shape) <- "shape"
else names(initial$shape) <- paste("xi", 0:(ncol(X.sh)-1), sep="")
} # end of if type is not light tailed stmts.
if(is.element(method, c("mle","gmle"))) {
if(use.phi) init.pars <- c(initial$location, initial$log.scale, initial$shape)
else init.pars <- c(initial$location, initial$scale, initial$shape)
if(type == "exponential" && const.scale) {
if( method == "gmle" ) warning( "Method MLE used." )
res <- list()
excess.id <- x > threshold
mle <- mean(x[excess.id] - threshold[excess.id])
names(mle) <- "scale"
res$par <- mle
k <- sum(excess.id)
res$n <- k
res$value <- k * ( log(mle) + 1 )
} else {
if(!is.null(a <- optim.args)) {
anam <- names(a)
if(!is.element("gr",anam)) {
if(method=="mle") opt.gr <- grlevd
else opt.gr <- NULL
} else opt.gr <- a$gr
if(is.null(a$method) && use.phi) opt.method <- ifelse(is.element(type,c("gev","gp","pp", "gumbel")), "BFGS", "L-BFGS-B")
else opt.method <- a$method
if(is.element(type, c("weibull","beta","frechet","pareto"))) opt.method <- "L-BFGS-B"
if(is.element(opt.method, c("L-BFGS-B","Brent"))) {
if(is.null(a$lower)) {
if(!is.element(type,c("frechet","pareto"))) opt.lower <- -Inf
else opt.lower <- c(rep(-Inf, length(init.pars)-1), 0)
} else {
if(is.element(type,c("frechet","pareto"))) opt.lower <- c(a$lower[1:(length(init.pars)-1)], 0)
else opt.lower <- a$lower
}
if(is.null(a$upper)) {
if(!is.element(type,c("weibull","beta"))) opt.upper <- Inf
else opt.upper <- c(rep(Inf, length(init.pars)-1), 0)
} else {
if(is.element(type,c("weibull","beta"))) opt.upper <- c(a$upper[1:(length(init.pars)-1)], 0)
else opt.upper <- a$upper
}
} else {
opt.lower <- -Inf
opt.upper <- Inf
} # end of if 'optimization method allows bounds' stmts.
if(is.null(a$control)) opt.control <- list()
else opt.control <- a$control
anam <- names(a$control)
if(!is.element("trace", anam) && verbose) opt.control$trace <- 6
if(is.null(a$hessian)) opt.hessian <- TRUE
else opt.hessian <- a$hessian
} else {
if(method=="mle") opt.gr <- grlevd
else opt.gr <- NULL
if(is.element(type, c("gev","gp","pp","gumbel"))) opt.method <- "BFGS"
else opt.method <- "L-BFGS-B"
if(!is.element(type,c("frechet","pareto"))) opt.lower <- -Inf
else opt.lower <- c(rep(-Inf, length(init.pars)-1), 0)
if(!is.element(type,c("weibull","beta"))) opt.upper <- Inf
else opt.upper <- c(rep(Inf, length(init.pars)-1), 0)
if(verbose) opt.control <- list(trace=6)
else opt.control <- list()
opt.hessian <- TRUE
} # end of if else 'optim.args' is passed stmts.
parnames <- names(init.pars)
out$parnames <- parnames
if(verbose && (method != "lmoments")) {
cat("Initial estimates are:\n")
print(init.pars)
cat("Beginning estimation procedure.\n")
}
if(type=="pp" && find.init) {
# Yet a third attempt to find good starting values--this time using the MLE for the GP parameters.
if(verbose) cat("\n", "First fitting a GP-Poisson model in order to try to get a good initial estimate as PP likelihoods can be very unstable.\n")
look <- out
look$type <- "GP"
des2 <- designs
des2$X.loc <- NULL
if(!missing(data)) resGP <- optim(init.pars[-(1:ncol(X.loc))], oevd, gr=opt.gr, o=look, des=des2, x=x, data=data, u=threshold, npy=npy, phi=use.phi, method=opt.method, lower=opt.lower, upper=opt.upper, control=opt.control, hessian=opt.hessian)
else resGP <- optim(init.pars[-(1:ncol(X.loc))], oevd, gr=opt.gr, o=look, des=des2, x=x, u=threshold, npy=npy, phi=use.phi, method=opt.method, lower=opt.lower, upper=opt.upper, control=opt.control, hessian=opt.hessian)
tmpi <- resGP$par
if(is.null(blocks)) { # CJP
lrate <- npy * mean(x > threshold)
} else {
lrate <- sum(x > threshold)/(blocks$nBlocks * mean(blocks$weights))
}
xi3 <- tmpi[(ncol(X.sc)+1):length(tmpi)]
if(!use.phi) sigma3 <- exp(tmpi[1:ncol(X.sc)] + xi3*log(lrate))
else sigma3 <- tmpi[1:ncol(X.sc)] + xi3*log(lrate)
lp <- lrate^(-xi3) - 1
if(all(is.finite(lp))) mu3 <- mean(threshold) - (sigma3/xi3)*lp # CJP fixed bug - "mean(threshold)" not "threshold"
else mu3 <- mean(x) # CJP: caution here as mean(x) will differ when only exceedances supplied...
nloc <- ncol(X.loc)
if(length(mu3) < nloc) mu3 <- c(mu3, rep(0, nloc-length(mu3)))
else mu3 <- mu3[1]
if(!is.null(blocks)) { # CJP
blocks$location <- rowSums(matrix(mu3, blocks$nBlocks, nloc)*blocks$designs$X.loc)
blocks$scale=rowSums(matrix(sigma3, blocks$nBlocks, ncol(blocks$designs$X.sc))*blocks$designs$X.sc)
blocks$shape=rowSums(matrix(xi3, blocks$nBlocks, ncol(blocks$designs$X.sh))*blocks$designs$X.sh)
}
if(all(is.finite(c(mu3, sigma3, xi3)))) {
testGPmle <- try(levd(x=x, threshold=threshold, location=rowSums(matrix(mu3, n, nloc)*X.loc),
scale=rowSums(matrix(sigma3, n, ncol(X.sc))*X.sc),
shape=rowSums(matrix(xi3, n, ncol(X.sh))*X.sh), type="PP", npy=npy, blocks=blocks), silent = TRUE)
if( is(testGPmle, "try-error" ) ) testGPmle <- Inf
} else testGPmle <- Inf
if(!is.null(blocks)) # CJP
blocks$location <- blocks$scale <- blocks$shape <- NULL
if(is.finite(testLmoments) || is.finite(testMoments) || is.finite(testGPmle)) {
if((testGPmle < testLmoments) && (testGPmle < testMoments)) {
if(verbose) cat("\n", "Changing initial estimates to those based on GP MLEs. They are: \n")
if(use.phi) init.pars <- c(mu3, log(sigma3), xi3)
else init.pars <- c(mu3, sigma3, xi3)
names(init.pars) <- parnames
if(verbose) print(init.pars)
} else if(verbose) cat("\n", "Sticking with originally chosen initial estimates.\n")
} # end of if any one of the three initial estimate attempts was successfull stmts.
inout$PoissonGP <- list(pars=c(mu3, sigma3, xi3), likelihood=testGPmle)
} # end of if 'type is PP, then try to improve initial estimates' stmts.
if(method=="mle") {
if(!missing(data)) {
res <- optim(init.pars, oevd, gr=opt.gr, o=out, des=designs, x=x, data=data, u=threshold,
span=span/npy, npy=npy, phi=use.phi, blocks=blocks, method=opt.method, lower=opt.lower,
upper=opt.upper, control=opt.control, hessian=opt.hessian) # CJP
} else {
res <- optim(init.pars, oevd, gr=opt.gr, o=out, des=designs, x=x, u=threshold,
span=span/npy, npy=npy, phi=use.phi, blocks=blocks, method=opt.method, lower=opt.lower,
upper=opt.upper, control=opt.control, hessian=opt.hessian)
}
} else if(method == "gmle") {
if(!missing(data)) {
res <- optim(init.pars, oevdgen, gr = opt.gr, o = out, des = designs, x = x, data = data,
u = threshold, span = span / npy, npy = npy, phi = use.phi, blocks = blocks,
priorFun = priorFun, priorParams = priorParams, method = opt.method, lower = opt.lower,
upper = opt.upper, control = opt.control, hessian = opt.hessian)
} else {
res <- optim(init.pars, oevdgen, gr = opt.gr, o = out, des = designs, x = x, u = threshold,
span = span / npy, npy = npy, phi = use.phi, blocks = blocks, priorFun = priorFun,
priorParams = priorParams, method = opt.method, lower = opt.lower, upper = opt.upper,
control = opt.control, hessian = opt.hessian)
}
} # end of if method is GMLE stmts.
} # end of if else type is exponential.
if(is.element("shape", names(res$par))) {
if(is.element(type, c("frechet","pareto"))) {
res$par["shape"] <- abs(res$par["shape"])
if(res$par["shape"] == 0) {
warning("fevd: shape parameter estimated to be zero. Re-setting to be 1e16.")
res$par["shape"] <- 1e16
}
} else {
if(is.element(type, c("weibull","beta"))) res$par["shape"] <- -abs(res$par["shape"])
if(res$par["shape"] == 0) {
warning("fevd: shape parameter estimated to be zero. Re-setting to be -1e16.")
res$par["shape"] <- -1e16
}
}
}
res$num.pars <- list(location=ncol(designs$X.loc), scale=ncol(designs$X.sc), shape=ncol(designs$X.sh))
out$results <- res
} else if(method=="bayesian") {
if(is.element(type, c("gev","gumbel","weibull","frechet","pp"))) {
nloc <- ncol(X.loc)
loc.names <- names(initial$location)
} else {
nloc <- 0
loc.names <- NULL
}
nsc <- ncol(X.sc)
if(use.phi && is.null(initial$log.scale)) {
initial$log.scale <- log(initial$scale)
if(nsc == 1) names(initial$log.scale) <- "log.scale"
else names(initial$log.scale) <- paste("phi", 0:(nsc-1), sep="")
}
sc.names <- names(initial$log.scale)
if(!is.element(type, c("gumbel","exponential"))) {
nsh <- ncol(X.sh)
sh.names <- names(initial$shape)
} else {
nsh <- 0
sh.names <- NULL
}
np <- nloc + nsc + nsh
find.priorParams <- FALSE
if(is.null(priorFun) && is.null(priorParams)) find.priorParams <- TRUE
else if(is.null(priorFun) && (is.null(priorParams$m) || is.null(priorParams$v))) find.priorParams <- TRUE
else if(!is.null(priorFun)) {
if(priorFun == "fevdPriorDefault") {
if(is.null(priorParams)) find.priorParams <- TRUE
else if(is.null(priorParams$m) || is.null(priorParams$v)) find.priorParams <- TRUE
}
}
if(is.null(priorFun) || find.priorParams) {
if(is.null(priorFun)) priorFun <- "fevdPriorDefault"
if(find.priorParams) {
xtemp <- x
class(xtemp) <- "mle"
if(verbose) cat("\n", "Finding MLE to obtain prior means and variances.\n")
if(missing(data)) {
if(missing(span)) hold <- initializer(xtemp, u=threshold, use.phi=use.phi, type=out$type, time.units=time.units, period.basis=period.basis, blocks=blocks)
else hold <- initializer(xtemp, u=threshold, use.phi=use.phi, type=out$type, span=span, time.units=time.units, period.basis=period.basis, blocks=blocks)
} else {
if(missing(span)) hold <- initializer(xtemp, data=data, u=threshold, u.fun=threshold.fun, loc.fun=location.fun, sc.fun=scale.fun, sh.fun=shape.fun,
use.phi=use.phi, type=out$type, time.units=time.units, period.basis=period.basis, blocks=blocks)
else hold <- initializer(xtemp, data=data, u=threshold, u.fun=threshold.fun, loc.fun=location.fun, sc.fun=scale.fun, sh.fun=shape.fun,
use.phi=use.phi, type=out$type, span=span, time.units=time.units, period.basis=period.basis, blocks=blocks)
}
if(is.null(priorParams)) priorParams <- list(m=hold[1:np], v=rep(10,np))
else if(is.null(priorParams$m)) priorParams$m <- hold[1:np]
else if(is.null(priorParams$v)) priorParams$v <- rep(10, np)
} # end of if 'priorParams' is not provided by user.
} # end of if 'priorFun' is not provided by user.
out$priorFun <- priorFun
out$priorParams <- priorParams
if(is.null(proposalFun)) {
proposalFun <- "fevdProposalDefault"
if(is.null(proposalParams)) proposalParams <- list(sd=rep(0.1,np))
} # end of if 'proposalFun' not supplied by user stmts.
out$proposalFun <- proposalFun
out$proposalParams <- proposalParams
# Eric 05/30/14
chain.info <- matrix(NA, iter, np + 2)
colnames(chain.info) <- c(loc.names, sc.names, sh.names, "loglik", "prior")
chain.info[2:iter, 1:np] <- 0
res <- matrix(NA, iter, np+1)
res[,np+1] <- 0
colnames(res) <- c(loc.names, sc.names, sh.names, "new")
if(nloc > 0) res[1,1:nloc] <- initial$location
if(use.phi) res[1,(nloc+1):(nloc+nsc)] <- initial$log.scale
else res[1,(nloc+1):(nloc+nsc)] <- initial$scale
if( type != "Gumbel" ) res[1,(nloc+nsc+1):np] <- initial$shape
theta.i <- res[1,1:np]
if(verbose) {
cat("\n", "Finding log-Likelihood of initial parameter values:\n")
print(theta.i)
}
if(!missing(data)) ll.i <- -oevd(p=res[1,], o=out, des=designs, x=x, data=data, u=threshold, span=span, npy=npy, phi=use.phi, blocks=blocks) # CJP
else ll.i <- -oevd(p=res[1,], o=out, des=designs, x=x, u=threshold, span=span, npy=npy, phi=use.phi, blocks=blocks)
if(verbose) cat("\n", "Finding prior df value of initial parameter values.\n")
p.i <- do.call(priorFun, c(list(theta=theta.i), priorParams))
chain.info[1, np + 1] <- ll.i
chain.info[1, np + 2] <- p.i
if(verbose) cat("\n", "Beginning the MCMC iterations (", iter, " total iterations)\n")
for(i in 2:iter) {
if(verbose && i <= 10) cat(i, " ")
if(verbose && i %% 100 == 0) cat(i, " ")
ord <- sample(1:np, np)
theta.star <- theta.i
acc <- 0
for(j in ord) {
# par.star <- do.call(proposalFun, c(list(p=theta.i[j], ind=j), proposalParams))
# theta.star[j] <- par.star
par.star <- do.call(proposalFun, c(list(p=theta.i, ind=j), proposalParams))
theta.star[j] <- par.star[j]
if(!missing(data)) ll.star <- -oevd(p=theta.star, o=out, des=designs, x=x, data=data, u=threshold, span=span, npy=npy, phi=use.phi, blocks=blocks) # CJP
else ll.star <- -oevd(p=theta.star, o=out, des=designs, x=x, u=threshold, span=span, npy=npy, phi=use.phi, blocks=blocks)
prior.star <- do.call(priorFun, c(list(theta=theta.star), priorParams))
look <- will.accept(ll.i=ll.i, prior.i=p.i, ll.star=ll.star, prior.star=prior.star, log=TRUE)
if(look$accept) {
p.i <- prior.star
ll.i <- ll.star
theta.i <- theta.star
acc <- acc+1
chain.info[i, j] <- 1
}
} # end of inner 'j' loop.
res[i,] <- c(theta.i, acc)
chain.info[i, (np+1):(np+2)] <- c(ll.i, p.i)
} # end of for 'i' loop.
if(verbose) cat("\n", "Finished MCMC iterations.\n")
out$results <- res
out$chain.info <- chain.info
} else stop("fevd: invalid method argument.")
out$initial.results <- inout
# if( method %in% c( "MLE", "GMLE" ) && results$convergence != 0 ) warning( print( results$convergence ) )
if(verbose) print(Sys.time() - begin.tiid)
if( method == "GMLE" ) cl <- "mle"
else cl <- tolower( method )
# class(out) <- c( "fevd", paste( "fevd.", cl ) )
class( out ) <- "fevd"
return(out)
} # end of 'fevd' function.
levd <- function(x, threshold, location, scale, shape,
type=c("GEV", "GP", "PP", "Gumbel", "Weibull", "Frechet", "Exponential", "Beta", "Pareto"),
log=TRUE, negative=TRUE, span, npy=365.25, infval=Inf, weights=1, blocks=NULL) {
# CJP added blocks argument
type <- match.arg(type)
type <- tolower(type)
n <- length(x)
w <- weights
if(length(w) == 1) w <- rep(w, n)
# if(type=="pp") {
# if(!missing(span)) {
# if((length(threshold)==1) && (length(location)==1) && (length(scale)==1)) shortcut <- TRUE
# else shortcut <- FALSE
# } else shortcut <- FALSE
# }
# CJP: if we do implement 'shortcut', can't we check the const.thresh, const.location, etc. variables? I guess the issue is that we need to pass them in to levd()
shortcut <- FALSE
if(!missing(threshold)) if(length(threshold)==1) threshold <- rep(threshold, n)
if(missing(threshold) && missing(location)) stop("levd: must supply one of threshold or location.")
if(missing(location) && !missing(threshold) && is.element(type, c("gev", "weibull", "gumbel", "frechet", "pp"))) location <- threshold
if(!missing(location)) if(length(location)==1) location <- rep(location, n)
if(missing(scale)) stop("levd: must supply a scale argument.")
else if(length(scale)==1) scale <- rep(scale, n)
if(!missing(shape)) {
if(length(shape)==1) shape <- rep(shape, n)
if(any(is.na(shape))) {
if(!negative) return(-infval)
else return(infval)
}
zero.shape <- shape==0
if(all(zero.shape) && !is.element(type, c("gumbel", "exponential"))) {
if(type=="gev") type <- "gumbel"
else if(type=="gp") type <- "exponential"
}
if(any(zero.shape) && any(!zero.shape)) {
warning("levd: some shape parameters are zero, and some not. Re-setting those that are to -1e-10.")
shape[zero.shape] <- -1e-10
}
} else shape <- 0
# end of if '!missing(shape)' stmts.
if(!is.null(blocks)){ # CJP
if(any(is.na(blocks$shape))) {
if(!negative) return(-infval)
else return(infval)
}
blocks$zero.shape <- blocks$shape==0
if(any(blocks$zero.shape) && any(!blocks$zero.shape)) {
warning("levd: some shape parameters are zero, and some not. Re-setting those that are to -1e-10.")
blocks$shape[blocks$zero.shape] <- -1e-10
}
}
if(any(scale <= 0)) {
if(!negative) return(-infval)
else return(infval)
}
if(!is.null(blocks) && any(blocks$scale <= 0)) { # CJP
if(!negative) return(-infval)
else return(infval)
}
# stop("levd: invalid scale parameter. Must be positively valued.")
if(is.element(type, c("gp","exponential","beta","pareto"))) {
excess.id <- x > threshold
threshold <- threshold[excess.id]
scale <- scale[excess.id]
shape <- shape[excess.id]
y <- (x[excess.id] - threshold)/scale
m <- length(y)
}
if(type=="gumbel") {
y <- (x - location)/scale
res <- sum(colSums(w * cbind(log(scale), y, exp(-y))))
} else if(is.element(type, c("gev","weibull","frechet"))) {
y <- (x - location)/scale
y <- 1 + shape*y
if(any(y <= 0)) {
if(!negative) return(-infval)
else return(infval)
}
if(type=="frechet") shape <- abs(shape)
else if(type=="weibull") shape <- -abs(shape)
res <- sum(w * cbind(log(scale), y^(-1/shape), log(y) * (1/shape + 1))) # CJP removed colSums() here
} else if(type=="exponential") res <- sum(w[excess.id] * cbind(log(scale), y)) # CJP removed colSums() here
else if(is.element(type,c("gp","beta","pareto"))) {
y <- 1 + shape*y
if(any(y <= 0)) {
if(!negative) return(-infval)
else return(infval)
}
if(type=="pareto") shape <- abs(shape)
else if(type=="beta") shape <- -abs(shape)
res <- sum(w[excess.id] * cbind(log(scale), log(y) * (1/shape + 1))) # CJP removed colSums()
} else if(type=="pp") {
if(any(zero.shape)) {
warning("levd: some shape parameters are zero. Re-setting them to -1e-10.")
shape[zero.shape] <- -1e-10
}
if(!is.null(blocks) && any(blocks$zero.shape)) { # CJP
warning("levd: some shape parameters are zero. Re-setting them to -1e-10.")
blocks$shape[blocks$zero.shape] <- -1e-10
}
excess.id <- x > threshold
y <- (x - location)/scale
y <- 1 + shape * y
if(!shortcut) z <- 1 + ((shape*(threshold - location))/scale) # CJP made this produce a vector not a matrix
else z <- 1 + ((shape[1]*(threshold[1] - location[1]))/scale[1])
if(any(y[excess.id] <= 0) || any(z < 0)) {
if(!negative) return(-infval)
else return(infval)
}
if(!is.null(blocks)){ # CJP
blocks$z <- 1 + ((blocks$shape*(blocks$threshold - blocks$location))/blocks$scale)
if(any(blocks$z < 0)) {
if(!negative) return(-infval)
else return(infval)
}
}
th.u <- threshold[excess.id]
sc.u <- scale[excess.id]
sh.u <- shape[excess.id]
y.u <- y[excess.id]
if(is.null(blocks)) { # CJP: sum(w*vec1 + w*vec2) may be faster than cbind
if(!shortcut) res <- sum(w[excess.id] * cbind(log(sc.u), log(y.u) * (1/sh.u + 1))) + sum(w * z^(-1/shape))/npy # CJP removed colSums
else res <- sum(w[excess.id] * cbind(log(sc.u), log(y.u) * (1/sh.u + 1))) + span * z^(-1/shape[1]) # CJP removed colSums
} else {
res <- sum(w[excess.id] * cbind(log(sc.u), log(y.u) * (1/sh.u + 1))) # pp density of exceedances
res <- res + blocks$nBlocks * mean(blocks$weights * (1 - blocks$proportionMissing) * blocks$z^(-1/blocks$shape)) # this works either if stationary and each element is a scalar or if one or more vectors, though at the moment, I think $z and $shape will always be passed in as vectors
}
} # end of which type stmts.
if(!negative) res <- -res
if(!log) res <- exp(res)
if(is.na(res)) {
if(!negative) return(-infval)
else return(infval)
}
return(res)
} # end of 'levd' function.
parcov.fevd <- function(x) {
cov.theta <- NULL
theta.hat <- x$results$par
theta.names <- names(theta.hat)
if(is.element("log.scale",theta.names)) {
theta.hat[theta.names=="log.scale"] <- exp(theta.hat[theta.names=="log.scale"])
theta.names[theta.names=="log.scale"] <- "scale"
names(theta.hat) <- theta.names
phiU <- FALSE
} else phiU <- x$par.models$log.scale
y <- datagrabber(x)
if(x$data.name[2] != "") {
xdat <- y[,1]
data <- y[,-1]
} else {
xdat <- y[,1]
data <- NULL
}
designs <- setup.design(x)
if(x$method != "GMLE") {
hold <- try(suppressWarnings(optimHess(theta.hat, oevd, gr=grlevd, o=x, des=designs, x=xdat, data=data, u=x$threshold, npy=x$npy, phi=phiU, blocks=x$blocks)), silent=TRUE) # CJP
if( !is(hold, "try-error" ) && all( !is.na(hold) ) ) {
cov.theta <- try(suppressWarnings(solve(hold)), silent=TRUE)
if(any(diag(cov.theta) <= 0)) re.do <- TRUE
else re.do <- FALSE
} else re.do <- TRUE
} else re.do <- TRUE
if(re.do) {
hold <- try(optimHess(theta.hat, oevd, o=x, des=designs, x=xdat, data=data, u=x$threshold, npy=x$npy, phi=phiU, blocks=x$blocks), silent=TRUE) # CJP
if( !is(hold, "try-error" ) && all(!is.na(hold)) ) {
cov.theta <- try(solve(hold), silent=TRUE)
if( is(cov.theta, "try-error" ) || any(diag(cov.theta) <= 0) ) cov.theta <- NULL
}
} # end of if 're.do' stmts.
return(cov.theta)
} # end of 'parcov.fevd' function.
distill.fevd <- function(x, ...) {
if( x$method=="GMLE" ) newcl <- "mle"
else newcl <- tolower( x$method )
# class(x) <- c( paste("fevd.", newcl, sep=""), class( x ) )
# UseMethod("distill",x)
res <- get( paste( "distill.fevd.", tolower( x$method ), sep = "" ) )( x = x, ... )
return( res )
} # end of 'distill.fevd' function.
distill.fevd.lmoments <- function(x, ...) {
return(x$results)
} # end of 'distill.fevd.lmoments' function.
distill.fevd.bayesian <- function(x, cov=TRUE, FUN="mean", burn.in=499, ...) {
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(cov) {
cov.theta <- cov(p)
out <- c(out, c(cov.theta))
names(out) <- c(pnames, paste(pnames[rep(1:np,np)], pnames[rep(1:np,each=np)], sep="."))
}
return(out)
} # end of 'distill.fevd.bayesian' function.
distill.fevd.mle <- function(x, cov=TRUE, ...) {
theta.hat <- x$results$par
theta.names <- names(theta.hat)
np <- length(theta.hat)
if(is.element("log.scale",theta.names)) {
theta.hat[theta.names=="log.scale"] <- exp(theta.hat[theta.names=="log.scale"])
theta.names[theta.names=="log.scale"] <- "scale"
names(theta.hat) <- theta.names
}
out <- c(theta.hat, x$results$value)
if(cov) {
cov.theta <- parcov.fevd(x)
if(!is.null(cov.theta)) {
out <- c(out, cov.theta)
names(out) <- c(theta.names, "nllh", paste(theta.names[rep(1:np,np)], theta.names[rep(1:np,each=np)], sep="."))
} else names(out) <- c(theta.names, "nllh")
} else names(out) <- c(theta.names, "nllh")
return(out)
} # end of 'distill.fevd.mle' function.
summary.fevd <- function(object, ...) {
if(object$method == "GMLE") newcl <- "mle"
else newcl <- tolower( object$method )
# class(object) <- c( paste("fevd.", newcl, sep=""), class( object ) )
# UseMethod("summary", object)
get( paste( "summary.fevd.", newcl, sep = "" ) )( object = object, ... )
} # end of 'summary.fevd' function.
summary.fevd.lmoments <- function(object, ...) {
x <- object
a <- list(...)
if(!is.null(a$silent)) silent <- a$silent
else silent <- FALSE
theta.hat <- x$results
if(!silent) {
print(x$call)
if(x$data.name[2]=="") print(paste(x$type, " Fitted to ", x$data.name[1], " using L-moments estimation."))
else print(paste(x$type, " Fitted to ", all.vars(x$x.fun), " of ", x$data.name[2], "data frame, using L-moments estimation."))
print(theta.hat)
}
invisible(theta.hat)
} # end of 'summary.fevd.lmoments' function.
summary.fevd.bayesian <- function(object, FUN="mean", burn.in=499, ...) {
x <- object
a <- list(...)
if(!is.null(a$silent)) silent <- a$silent
else silent <- FALSE
np <- dim(x$results)[2] - 1
if(!is.na(burn.in) && !is.null(burn.in) && burn.in > 0) {
ll <- x$chain.info[-(1:burn.in), np + 1]
mth <- colMeans(x$results[-(1:burn.in), -(np + 1)], na.rm = TRUE)
} else {
ll <- x$chain.info[, np + 1]
mth <- colMeans(x$results[, -(np + 1)], na.rm = TRUE)
} # end of if else 'burn.in' stmts.
const <- is.fixedfevd(x)
if(!silent) {
cat("\n")
print(x$call)
cat("\n")
print(paste("Estimation Method used: ", x$method, sep=""))
cat("\n")
}
out <- list()
out$par <- ci(x, type="parameter", FUN=FUN, burn.in=burn.in)
# if(const) out$rl <- ci(x, return.period=c(2,20,100), FUN=FUN, burn.in=burn.in)
pdim <- dim(x$results)
# out$acceptance.rate <- colMeans(x$results[2:pdim[1], -pdim[2]] != x$results[1:(pdim[1] - 1), -pdim[2] ])
out$acceptance.rate <- colMeans(x$chain.info[, 2:(pdim[2] - 1), drop = FALSE ], na.rm = TRUE)
if(!silent) {
cat("\n", "Acceptance Rates:\n")
print(out$acceptance.rate)
print(out$par)
# if(const) print(out$rl)
}
cov.theta <- distill(x, FUN=FUN, burn.in=burn.in)
cov.theta <- matrix(cov.theta[(np + 1):length(cov.theta)], np, np)
colnames(cov.theta) <- rownames(cov.theta) <- colnames(x$results)[1:np]
if(!silent) {
cat("\n", "Estimated parameter covariance matrix.\n")
print(cov.theta)
}
out$cov.theta <- cov.theta
out$se.theta <- sqrt(diag(cov.theta))
designs <- setup.design(x)
y <- c(datagrabber(x, cov.data=FALSE))
if(x$data.name[2] == "") data <- NULL
else data <- datagrabber(x, response=FALSE)
if(is.element(x$type, c("GEV", "Gumbel", "Weibull", "Frechet"))) n <- x$n
else n <- sum(y > x$threshold)
nllh <- oevd(p=out$par, o=x, des=designs, x=y, data=data, u=x$threshold, span=x$span, npy=x$npy, phi=x$par.models$log.scale, blocks=x$blocks) # CJP
out$nllh <- nllh
# out$AIC <- 2 * nllh + 2 * np
if(FUN == "postmode" || FUN == "mode") out$BIC <- 2 * nllh + np * log(n)
Dev <- -2 * ll
Dbar <- mean(Dev, na.rm = TRUE)
Dmth <- -2 * oevd(p=mth, o=x, des=designs, x=y, data=data, u=x$threshold, span=x$span, npy=x$npy, phi=x$par.models$log.scale, blocks=x$blocks)
pd <- Dbar - Dmth
out$DIC <- Dmth + 2 * pd
if(!silent) {
# cat("\n", "AIC =", out$AIC, "\n")
# cat("\n", "BIC =", out$BIC, "\n")
cat("\n", "DIC = ", out$DIC, "\n")
}
invisible(out)
} # end of 'summary.fevd.bayesian' function.
summary.fevd.mle <- function(object, ...) {
x <- object
a <- list(...)
out <- list()
cov.theta <- se.theta <- NULL
if(!is.null(a$silent)) silent <- a$silent
else silent <- FALSE
if(!silent) {
cat("\n")
print(x$call)
cat("\n")
print(paste("Estimation Method used: ", x$method, sep=""))
cat("\n")
}
if(!silent) cat("\n", "Negative Log-Likelihood Value: ", x$results$value, "\n\n")
theta.hat <- x$results$par
theta.names <- names(theta.hat)
if(is.element("log.scale",theta.names)) {
theta.hat[theta.names=="log.scale"] <- exp(theta.hat[theta.names=="log.scale"])
theta.names[theta.names=="log.scale"] <- "scale"
names(theta.hat) <- theta.names
phiU <- FALSE
} else phiU <- x$par.models$log.scale
out$par <- theta.hat
np <- length(theta.hat)
designs <- setup.design(x)
if(!is.null(x$data.pointer)) xdat <- get(x$data.pointer)
else xdat <- x$x
cov.theta <- parcov.fevd(x)
out$cov.theta <- cov.theta
if(!is.null(cov.theta)) {
se.theta <- sqrt(diag(cov.theta))
names(se.theta) <- theta.names
out$se.theta <- se.theta
}
if(!silent) {
cat("\n", "Estimated parameters:\n")
print(theta.hat)
}
if(!is.null(se.theta)) {
if(!silent) {
cat("\n", "Standard Error Estimates:\n")
print(se.theta)
}
theta.hat <- rbind(theta.hat, se.theta)
if(is.matrix(theta.hat) && dim(theta.hat)[1]==2) rownames(theta.hat) <- c("Estimate", "Std. Error")
if(is.matrix(theta.hat)) colnames(theta.hat) <- theta.names
else names(theta.hat) <- theta.names
if(!silent && !is.null(cov.theta)) {
cat("\n", "Estimated parameter covariance matrix.\n")
print(cov.theta)
}
}
nllh <- x$results$value
out$nllh <- nllh
if(is.element(x$type, c("GEV","Gumbel","Weibull","Frechet"))) n <- x$n
else {
y <- c(datagrabber(x, cov.data=FALSE))
n <- sum(y > x$threshold)
}
out$AIC <- 2 * nllh + 2 * np
out$BIC <- 2 * nllh + np * log(n)
if(!silent) {
cat("\n", "AIC =", out$AIC, "\n")
cat("\n", "BIC =", out$BIC, "\n")
}
invisible(out)
} # end of 'summary.fevd.mle' function.
print.fevd <- function(x, ...) {
summary(x)
invisible()
} # end of 'print.fevd' function.
grlevd <- function(p, o, des, x, data, u=NULL, span, npy, phi=FALSE, blocks=NULL) {
##
## Function to estimate the likelihood
## parameter gradients.
##
npar <- length(p)
type <- o$type
n <- length(x)
w <- o$weights
if(is.null(w)) w <- 1
if(length(w) == 1) w <- rep(w, n)
shortcut <- FALSE
dpar <- matrix(NA, n, npar)
if(!is.null(blocks)) # CJP
blocks$dpar <- matrix(NA, blocks$nBlocks, npar)
if(is.element(type, c("PP","GEV","Gumbel","Weibull","Frechet"))) {
X.loc <- des$X.loc
nloc <- ncol(X.loc)
loc <- rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc)
dpar[,1:nloc] <- X.loc
if(!is.null(blocks)) { # CJP
blocks$loc <- rowSums(matrix(p[1:nloc], blocks$nBlocks, nloc, byrow=TRUE)*blocks$designs$X.loc)
blocks$dpar[,1:nloc] <- blocks$designs$X.loc
}
} else {
nloc <- 0
loc <- NULL
if(!is.null(blocks)) # CJP
blocks$loc <- NULL
}
X.sc <- des$X.sc
nsc <- ncol(X.sc)
if(phi) scale <- exp(rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc))
else scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc)
dpar[,(nloc+1):(nloc+nsc)] <- X.sc
if(phi) phi.gradTerm <- scale else phi.gradTerm <- 1 # CJP2
if(!is.null(blocks)) { # CJP
if(phi) {
blocks$scale <- exp(rowSums(matrix(p[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE)*blocks$designs$X.sc))
blocks$phi.gradTerm <- blocks$scale
} else {
blocks$scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE)*blocks$designs$X.sc)
blocks$phi.gradTerm <- 1
}
blocks$dpar[,(nloc+1):(nloc+nsc)] <- blocks$designs$X.sc
}
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)
dpar[,(nloc+nsc+1):(nloc+nsc+nsh)] <- X.sh
if(!is.null(blocks)) { # CJP
blocks$shape <- rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], blocks$nBlocks, nsh, byrow=TRUE)*blocks$designs$X.sh)
blocks$dpar[,(nloc+nsc+1):(nloc+nsc+nsh)] <- blocks$designs$X.sh
}
} else {
nsh <- 0
shape <- NULL
if(!is.null(blocks)) # CJP
blocks$shape <- NULL
}
if(is.element(type, c("Frechet","Pareto"))) shape <- abs(shape)
else if(is.element(type, c("Weibull","Beta"))) shape <- -abs(shape)
res <- p+NA
if(is.element(type, c("GEV","Weibull","Frechet"))) {
sh.sc <- shape/scale
sh.1 <- 1/shape + 1
x.fmu <- x - loc
z <- 1 + sh.sc*x.fmu
for(i in 1:npar) { # CJP: don't need colSums and sum(w*vec1 + w*vec2 + w*vec3) may be faster than cbind
if(i <= nloc) res[i] <- sum(colSums(w * cbind(-sh.1 * (sh.sc) * dpar[,i]/z, (dpar[,i]/scale)*z^(-sh.1))))
else if(i <= nloc + nsc) res[i] <- sum(colSums(w * cbind(phi.gradTerm*dpar[,i]/scale, -phi.gradTerm * sh.1 * (shape/(scale^2))*x.fmu*dpar[,i]/z, (z^(-sh.1)/(scale^2)) * dpar[,i] * x.fmu * phi.gradTerm))) # CJP2
else res[i] <- sum(colSums(w * cbind(sh.1 * (dpar[,i]/scale) * x.fmu/z, -dpar[,i]/(shape^2)*log(z),
z^(-1/shape) * ((dpar[,i]/(shape^2))*log(z) - (dpar[,i]/scale)*x.fmu/(z*shape)))))
} # end of for 'i' loop.
} else if(type == "Gumbel") {
x.fmu <- x - loc
z <- x.fmu/scale
for(i in 1:npar) {
if(i <= nloc) res[i] <- sum(colSums(w * cbind(-dpar[,i]/scale, (dpar[,i]/scale)*exp(-z))))
else res[i] <- sum(colSums(w * cbind(phi.gradTerm*dpar[,i]/scale, -((x.fmu/(scale^2))*dpar[,i])*phi.gradTerm, phi.gradTerm*(dpar[,i]/(scale^2))*x.fmu*exp(-z)))) # CJP2 - Eric, please check this
} # end of for 'i' loop.
} else if(is.element(type, c("GP","Beta","Pareto"))) {
excess.id <- x > u
x <- x[excess.id]
scale <- scale[excess.id]
shape <- shape[excess.id]
u <- u[excess.id]
dparu <- dpar[excess.id,]
x.fu <- x - u
z <- 1 + (shape/scale)*x.fu
sh.1 <- 1/shape + 1
if(phi) phi.gradTermExc <- scale else phi.gradTermExc <- 1 # CJP2
for(i in 1:npar) {
if(i <= nsc) res[i] <- sum(colSums(w[excess.id] * cbind(phi.gradTermExc * dparu[,i]/scale, - ((phi.gradTermExc*sh.1*dparu[,i]*shape*x.fu/(scale^2))/z)))) # CJP2 - Eric, please check this
else res[i] <- sum(colSums(w[excess.id] * cbind((-dparu[,i]/(shape^2))*log(z), sh.1*((dparu[,i]/scale)*x.fu)/z)))
} # end of for 'i' loop.
} else if(type=="Exponential") {
excess.id <- x > u
if(!is.matrix(dpar)) dpar <- cbind(dpar)
dparu <- dpar[excess.id,]
if(!is.matrix(dparu)) dparu <- cbind(dparu)
scale <- scale[excess.id]
xu <- x[excess.id] - u[excess.id]
if(phi) phi.gradTermExc <- scale else phi.gradTermExc <- 1 # CJP2
for(i in 1:npar) res[i] <- sum(colSums(w[excess.id] * cbind(phi.gradTermExc*dparu[,i]/scale, - phi.gradTermExc * xu * dparu[,i]/(scale^2)))) # CJP2 - Eric, please check this
} else if(type=="PP") {
excess.id <- x > u
mu <- loc[excess.id]
sc <- scale[excess.id]
sh <- shape[excess.id]
shu.scu <- sh/sc
uu <- u[excess.id]
xu <- x[excess.id]
dparu <- dpar[excess.id,]
x.fmu <- xu - mu
if(!shortcut) u.l <- u - loc
else u.l <- u[1] - loc[1]
zu <- 1 + shu.scu*x.fmu
if(!shortcut) Z <- matrix(1 + (shape/scale)*u.l, ncol=1)
else Z <- 1 + (shape[1]/scale[1])*u.l
shu.1 <- 1/sh + 1
if(!shortcut) sh.1 <- 1/shape + 1
else sh.1 <- 1/shape[1] + 1
zpsh1 <- Z^(-sh.1)
if(phi) phi.gradTermExc <- sc else phi.gradTermExc <- 1 # CJP2
if(!is.null(blocks)) {
blocks$u.l <- blocks$threshold - blocks$loc
blocks$Z <- 1 + (blocks$shape/blocks$scale)*blocks$u.l
blocks$sh.1 <- 1/blocks$shape + 1
blocks$zpsh1 <- blocks$Z^(-blocks$sh.1)
}
for(i in 1:npar) { # CJP: I changed sum(colSums(...)) here to sum()
if(!shortcut) { # CJP: same comment about sum(w*vec1 + w*vec2)
if(is.null(blocks)) { # CJP
if(i <= nloc) res[i] <- sum(w[excess.id] * cbind(-shu.1*dparu[,i]*shu.scu/zu)) + colSums(w * zpsh1 * dpar[,i]/scale)/npy # CJP: I think this can be sum not colSums
else if(i <= nloc + nsc) res[i] <- sum(w[excess.id] * cbind(phi.gradTermExc * dparu[,i]/sc, -(shu.1*dparu[,i]*sh*x.fmu/(zu*sc^2))*phi.gradTermExc)) + colSums(w * phi.gradTerm * dpar[,i]*u.l*zpsh1/(scale^2))/npy # CJP: I think this can be sum not colSums; CJP2 - added phi.gradTerm
else res[i] <- sum(w[excess.id] * cbind(-(dparu[,i]/(sh^2))*log(zu), (shu.1*(dparu[,i]/sc)*x.fmu)/zu)) +
colSums(w * (exp(-log(Z)/shape)*((dpar[,i]/(shape^2))*log(Z) - (dpar[,i]/(shape*scale)*u.l)/Z)))/npy
} else { # w/ blocks
if(i <= nloc) res[i] <- sum(w[excess.id] * cbind(-shu.1*dparu[,i]*shu.scu/zu)) +
blocks$nBlocks * mean(blocks$weights * (1 - blocks$proportionMissing) * blocks$zpsh1 * blocks$dpar[,i]/blocks$scale)
else if(i <= nloc + nsc) res[i] <- sum(w[excess.id] * cbind(phi.gradTermExc*dparu[,i]/sc, -(shu.1*dparu[,i]*sh*x.fmu/(zu*sc^2))*phi.gradTermExc)) +
blocks$nBlocks * mean(blocks$weights * (1 - blocks$proportionMissing) * blocks$phi.gradTerm * blocks$dpar[,i]*blocks$u.l*blocks$zpsh1/(blocks$scale^2)) # CJP2: added blocks$phi.gradTerm
else res[i] <- sum(w[excess.id] * cbind(-(dparu[,i]/(sh^2))*log(zu), (shu.1*(dparu[,i]/sc)*x.fmu)/zu)) +
blocks$nBlocks * mean(blocks$weights * (1 - blocks$proportionMissing) * (exp(-log(blocks$Z)/blocks$shape)*((blocks$dpar[,i]/(blocks$shape^2))*log(blocks$Z) - (blocks$dpar[,i]/(blocks$shape*blocks$scale)*blocks$u.l)/blocks$Z)))
}
} else {
if(i <= nloc) res[i] <- -sum(w[excess.id] * cbind(shu.1*dparu[,i]*shu.scu/zu)) + span * (zpsh1*dpar[1,i]/scale[1])
else if(i <= nloc + nsc) res[i] <- sum(w[excess.id] * cbind(phi.gradTermExc*dparu[,i]/sc, -(shu.1*dparu[,i]*sh*x.fmu/(zu*sc^2))*phi.gradTermExc)) + span * (phi.gradTerm*dpar[1,i]*u.l*zpsh1/(scale[1]^2)) # CJP2: added phi.gradTerm
else res[i] <- sum(w[excess.id] * cbind(-(dparu[,i]/(sh^2))*log(zu), shu.1*(dparu[,i]/sc)*x.fmu/zu)) +
span * (exp(-log(Z)/shape[1])*((dpar[1,i]/(shape[1]^2))*log(Z) - (dpar[1,i]/(shape[1]*scale[1])*u.l)/Z))
}
} # end of for 'i' loop.
} # end of which gradient likelihood to use stmts.
return(res)
} # end of 'grlevd' function.
grlevdTracer <- function(x, p=c(0,1,0), which.vary=1, p1.range, p2.range=NULL, threshold=NULL,
threshold.fun=~1, location.fun=~1, scale.fun=~1, shape.fun=~1, data, phi=FALSE,
type=c("GEV","GP","PP","Gumbel","Weibull","Frechet","Exponential","Beta","Pareto"),
nint=100, npy=365.25, weights = 1, blocks=NULL, na.action=na.fail, par1.name="", par2.name="", zl.lh=NULL,
plot=TRUE, ...) {
# CJP note: this function now assumes that blocks contains the blocks data and design matrices already constructed within fevd and passed out as the $blocks element of the output from fevd - I didn't want to repeat all of that code here
# CJP2 : changed 'u' to 'threshold' and 'threshold' to 'threshold.fun'
type <- match.arg(type)
id <- which.vary
data.name <- c(deparse(substitute(x)), deparse(substitute(data)))
if(!is.function(na.action)) na.action <- get(na.action)
if(!missing(data)) {
if(length(x) != nrow(data)) stop("fevd: data must have same number of rows as the length of x.")
tmp <- cbind(x, data)
tmp <- na.action(tmp)
x <- tmp[,1]
data <- tmp[,-1]
n <- length(x)
} else {
x <- na.action(x)
n <- length(x)
}
w <- weights
if(length(w) == 1) w <- rep(w, n)
if(!is.null(threshold)) {
const.thresh <- check.constant(threshold.fun)
if(!missing(data)) X.u <- setup.design(x=threshold.fun, data=data, n=n, const=const.thresh, dname="threshold")
else X.u <- setup.design(x=threshold.fun, n=n, const=const.thresh, dname="threshold")
threshold <- rowSums(threshold*X.u)
} # end of if 'threshold (u)' stmts.
p1 <- seq(p1.range[1],p1.range[2],,nint)
if(!is.null(p2.range)) p2 <- seq(p2.range[1],p2.range[2],,nint)
loc.exists <- is.element(type, c("GEV","PP","Gumbel","Weibull","Frechet"))
sh.exists <- !is.element(type, c("Gumbel","Exponential"))
designs <- list()
if(loc.exists) {
const.loc <- check.constant(location.fun)
if(!missing(data)) X.loc <- setup.design(x=location.fun, data=data, n=n, const=const.loc, dname="location.fun")
else X.loc <- setup.design(x=location.fun, n=n, const=const.loc, dname="location")
designs$X.loc <- X.loc
nloc <- ncol(X.loc)
} else {
nloc <- 0
const.loc <- NULL
}
const.scale <- check.constant(scale.fun)
if(!missing(data)) X.sc <- setup.design(x=scale.fun, data=data, n=n, const=const.scale, dname="scale.fun")
else X.sc <- setup.design(x=scale.fun, n=n, const=const.scale, dname="scale.fun")
designs$X.sc <- X.sc
nsc <- ncol(X.sc)
if(sh.exists) {
const.shape <- check.constant(shape.fun)
if(!missing(data)) X.sh <- setup.design(x=shape.fun, data=data, n=n, const=const.shape, dname="shape.fun")
else X.sh <- setup.design(x=shape.fun, n=n, const=const.shape, dname="shape")
designs$X.sh <- X.sh
nsh <- ncol(X.sh)
if(is.element(type, c("Frechet","Pareto"))) shape <- abs(shape)
else if(is.element(type, c("Weibull","Beta"))) shape <- -abs(shape)
} else {
nsh <- 0
const.shape <- NULL
shape <- 0
} # end of 'shape.range' stmts.
obj <- list(type=type, const.loc=const.loc, const.scale=const.scale, const.shape=const.shape,
par.models=list(location=location.fun, scale=scale.fun, shape=shape.fun, phi=phi, weights=w,
term.names=list(location=all.vars(location.fun), scale=all.vars(scale.fun), shape=all.vars(shape.fun))))
num.vary <- length(id)
if(num.vary==1) {
lh <- lhgr <- numeric(nint)+NA
for(i in 1:nint) {
p[id] <- p1[i]
p.tmp <- p
scale <- rowSums(matrix(p.tmp[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc)
if(phi) scale <- exp(scale)
if(!is.null(blocks)) {
blocks$location <- rowSums(matrix(p[1:nloc], blocks$nBlocks, nloc, byrow=TRUE)*blocks$designs$X.loc)
blocks$scale <- rowSums(matrix(p.tmp[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE)*blocks$designs$X.sc)
if(phi) blocks$scale <- exp(blocks$scale)
blocks$shape=rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], blocks$nBlocks, nsh, byrow=TRUE)*blocks$designs$X.sh)
}
if(loc.exists && sh.exists) lh[i] <- levd(x=x, threshold=threshold, location=rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc),
scale=scale, shape=rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], n, nsh, byrow=TRUE)*X.sh),
type=type, npy=npy, weights=w, blocks=blocks)
else if(loc.exists) lh[i] <- levd(x=x, threshold=threshold, location=rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc),
scale=scale, type=type, npy=npy, weights=w, blocks=blocks)
else if(sh.exists) lh[i] <- levd(x=x, threshold=threshold, scale=scale, shape=rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], n, nsh, byrow=TRUE)*X.sh),
type=type, npy=npy, weights=w)
else lh[i] <- levd(x=x, threshold=threshold, scale=scale, type=type, npy=npy, weights=w, blocks=blocks)
if(!missing(data)) lhgr[i] <- grlevd(p=p.tmp, o=obj, des=designs, x=x, data=data, u=threshold, npy=npy, phi=phi, blocks=blocks)[id]
else lhgr[i] <- grlevd(p=p.tmp, o=obj, des=designs, x=x, u=threshold, npy=npy, phi=phi, blocks=blocks)[id]
} # end of for 'i' loop.
res <- cbind(p1, lh, lhgr)
colnames(res) <- c("par.values","nllh","nllh.gr")
attributes(res)$data.name <- data.name
attributes(res)$par.name <- par1.name
attributes(res)$model.name <- type
} else if(num.vary==2) {
lh <- lhgr1 <- lhgr2 <- matrix(NA, nint, nint)
for(i in 1:nint) {
p[id[1]] <- p1[i]
p.tmp <- p
if(const.scale && phi) p.tmp[nloc+1] <- exp(p.tmp[nloc+1])
for(j in 1:nint) {
p[id[2]] <- p2[j]
p.tmp <- p
scale <- rowSums(matrix(p.tmp[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc)
if(phi) scale <- exp(scale)
if(!is.null(blocks)) {
blocks$location <- rowSums(matrix(p[1:nloc], blocks$nBlocks, nloc, byrow=TRUE)*blocks$designs$X.loc)
blocks$scale <- rowSums(matrix(p.tmp[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE)*blocks$designs$X.sc) # fixed by CJP2 6/11/13 (was 'blocks$X.sc')
if(phi) blocks$scale <- exp(blocks$scale)
blocks$shape=rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], blocks$nBlocks, nsh, byrow=TRUE)*blocks$designs$X.sh)
}
if(loc.exists && sh.exists) lh[i,j] <- levd(x=x, threshold=threshold, location=rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc),
scale=scale, shape=rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], n, nsh, byrow=TRUE)*X.sh),
type=type, npy=npy, weights=w, blocks=blocks)
else if(loc.exists) lh[i,j] <- levd(x=x, threshold=threshold, location=rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc),
scale=scale, type=type, npy=npy, weights=w, blocks=blocks)
else if(sh.exists) lh[i,j] <- levd(x=x, threshold=threshold, scale=scale, shape=rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], n, nsh, byrow=TRUE)*X.sh),
type=type, npy=npy, weights=w)
else lh[i,j] <- levd(x=x, threshold=threshold, scale=scale, type=type, npy=npy, weights=w, blocks=blocks)
if(!missing(data)) hold <- grlevd(p=p.tmp, o=obj, des=designs, x=x, data=data, u=threshold, npy=npy, blocks=blocks)[id]
else hold <- grlevd(p=p.tmp, o=obj, des=designs, x=x, u=threshold, npy=npy, blocks=blocks)[id]
lhgr1[i,j] <- hold[1]
lhgr2[i,j] <- hold[2]
} # end of for 'j' loop.
} # end of for 'i' loop.
par.values <- list(p1,p2)
names(par.values) <- c(paste("gr1.", par1.name, sep=""), paste("gr2.", par2.name, sep=""))
res <- list(nllh=lh, nllh.gr1=lhgr1, nllh.gr2=lhgr2, p1=p1, p2=p2, call=match.call(), data.name=data.name, par.name=c(par1.name, par2.name), model=type)
parvar <- cbind(rep(p1, nint), rep(p2, each=nint))
} # end of 1 or two varied parameters stmts.
class(res) <- "grlevdTracer"
if(plot) plot(res, type="both", which.gr=1:2, ...)
invisible(res)
} # end of 'grlevdTracer' function.
plot.grlevdTracer <- function(x, type=c("likelihood", "gradient", "both"), which.gr=1, log=TRUE, negative=TRUE, ...) {
type <- match.arg(type)
a <- list(...)
if(is.element(type, c("likelihood","both")) && is.null(a$ylab)) {
if(log && negative) nllh.ylab <- "Negative Log-Likelihood"
else if(log) nllh.ylab <- "Log-Likelihood"
else if(!negative) nllh.ylab <- "Likelihood"
else nllh.ylab <- "Negative Likelihood"
}
if(is.element(type, c("gradient", "both"))) gr.ylab <- "Negative Log-Likelihood Gradient"
if(is.matrix(x)) {
tmp <- attributes(x)
dname <- tmp$data.name[1]
pname <- tmp$par.name
mname <- tmp$model.name
if(is.element(type, c("likelihood","both"))) {
X <- x[,"nllh"]
if(!any(is.finite(X))) {
skip <- TRUE
} else {
skip <- FALSE
if(!log && !negative) X <- exp(-X)
else if(!negative) X <- -X
else if(!log) X <- -exp(-X)
}
}
if(is.null(a$ylim)) {
if(is.element(type, c("likelihood","both"))) yl <- range(X, finite=TRUE)
if(type=="both") par(mfrow=c(2,1))
if(is.element(type, c("likelihood","both"))) {
if(skip) plot(1:length(X), rep(0, length(X)), type="n", yaxt="n", main="No finite negative log-likelihood values for this range")
else {
if(is.null(a$ylab) && is.null(a$main)) {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", ylab=nllh.ylab, main=paste(mname, dname, sep="\n"), ylim=yl, ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, ylab=nllh.ylab, main=paste(mname, dname, sep="\n"), ylim=yl, ...)
else plot(x[,"par.values"], X, type="l", ylab=nllh.ylab, main=paste(mname, dname, sep="\n"), ylim=yl, ...)
} else if(is.null(a$ylab)) {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", ylab=nllh.ylab, ylim=yl, ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, ylab=nllh.ylab, ylim=yl, ...)
else plot(x[,"par.values"], X, type="l", ylab=nllh.ylab, ylim=yl, ...)
} else if(is.null(a$main)) {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", main=paste(mname, dname, sep="\n"), ylim=yl, ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, main=paste(mname, dname, sep="\n"), ylim=yl, ...)
else plot(x[,"par.values"], X, type="l", ylim=yl, main=paste(mname, dname, sep="\n"), ...)
} else {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", ylim=yl, ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, ylim=yl, ...)
else plot(x[,"par.values"], X, type="l", ylim=yl, ...)
}
} # end of if else skip stmts.
} # end of if type is likelihood or both stmts.
if(is.element(type, c("gradient", "both"))) {
if(!any(is.finite(x[,"nllh.gr"]))) {
plot(1:length(x[,"nllh.gr"]), rep(0, length(x[,"nllh.gr"])), type="n", yaxt="n", main="No finite negative log-likelihood gradient values for this range")
} else {
if(is.null(a$main) && type=="both") m1 <- ""
else if(is.null(a$main)) m1 <- paste(mname, dname, sep="\n")
if(is.null(a$ylab) && is.null(a$main)) {
# gr.ylab <- paste(gr.ylab, " wrt ", pname, sep="")
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, ylab=gr.ylab, main=m1, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", ylab=gr.ylab, main=m1, ...)
} else if(is.null(a$ylab)) {
# gr.ylab <- paste(gr.ylab, " wrt ", pname, sep="")
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, ylab=gr.ylab, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", ylab=gr.ylab, ...)
} else if(is.null(a$main)) {
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, main=m1, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", main=m1, ...)
} else {
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", ...)
} # end of if ylab and main supplied or not stmts.
abline(h=0, lty=2)
} # end of if else no finite gradient values stmts.
} # end of if type is gradient or both stmts.
} else {
if(type=="both") par(mfrow=c(2,1))
if(is.element(type, c("likelihood","both"))) {
if(is.null(a$ylab) && is.null(a$main)) {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", ylab=nllh.ylab, main=paste(mname, dname, sep="\n"), ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, ylab=nllh.ylab, main=paste(mname, dname, sep="\n"), ...)
else plot(x[,"par.values"], X, type="l", ylab=nllh.ylab, main=paste(mname, dname, sep="\n"), ...)
} else if(is.null(a$ylab)) {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", ylab=nllh.ylab, ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, ylab=nllh.ylab, ...)
else plot(x[,"par.values"], X, type="l", ylab=nllh.ylab, ...)
} else {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, ...)
else plot(x[,"par.values"], X, type="l", ...)
}
} # end of if type is likelihood or both stmts.
if(is.element(type, c("gradient", "both"))) {
if(is.null(a$main) && type=="both") m1 <- ""
else if(is.null(a$main)) m1 <- paste(mname, dname, sep="\n")
if(is.null(a$ylab) && is.null(a$main)) {
gr.ylab <- paste(gr.ylab, " wrt ", pname, sep="")
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, ylab=gr.ylab, main=m1, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", ylab=gr.ylab, main=m1, ...)
} else if(is.null(a$ylab)) {
gr.ylab <- paste(gr.ylab, " wrt ", pname, sep="")
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, ylab=gr.ylab, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", ylab=gr.ylab, ...)
} else if(is.null(a$main)) {
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, main=m1, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", main=m1, ...)
} else {
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", ...)
} # end of if ylab and main supplied or not stmts.
abline(h=0, lty=2)
} # end of if type is gradient or both stmts.
} # end of if ylim is passed or not stmts.
} else {
dname <- x$data.name[1]
p1name <- x$par.name[1]
p2name <- x$par.name[2]
p1 <- x$p1
p2 <- x$p2
nint <- length(p1)
if(is.element(type, c("likelihood","both"))) {
X <- x$nllh
if(!log && !negative) X <- exp(-X)
else if(!negative) X <- -X
else if(!log) X <- -exp(-X)
minloc <- cbind(rep(p1, nint), rep(p2, each=nint))[which.min(c(X)),]
}
if(!(which.gr[1]==1 || which.gr[1]==2 || all(which.gr==1:2))) stop("plot.grlevdTracer: invalid which.gr argument. Must be 1, 2 or 1:2.")
if(is.null(a$xlab)) xl <- p1name
else xl <- a$xlab
if(is.null(a$ylab)) yl <- p2name
else yl <- a$ylab
if(!is.null(a$zlim)) zl.lh <- zl.gr1 <- zl.gr2 <- a$zlim
else {
zl.lh <- range(c(X), finite=TRUE)
zl.gr1 <- range(c(x$nllh.gr1), finite=TRUE)
zl.gr2 <- range(c(x$nllh.gr2), finite=TRUE)
}
# if(is.null(a$col)) zcol <- tim.colors(64)
# else zcol <- a$col
if(is.null(a$horizontal)) hz <- TRUE
else hz <- a$horizontal
if(is.null(a$legend.mar)) legmar <- ifelse(hz, 3.1, 5.1)
else legmar <- a$legend.mar
if(type=="both") {
if(length(which.gr)==1) par(mfrow=c(1,2))
else par(mfrow=c(1,3))
}
if(is.null(a$main)) {
if(is.element(type, c("likelihood","both"))) m1 <- paste(nllh.ylab, "\n", x$model, ": ", x$data.name[1], sep="")
else m1 <- ""
if(length(which.gr)==1) m2 <- m3 <- paste(x$model, ": Negative Log-Likelihood\nGradient wrt ", x$par.name[which.gr], " (", x$data.name[1], ")", sep="")
else {
m2 <- paste(x$model, ": Negative Log-Likelihood\nGradient wrt ", p1name, " (", x$data.name[1], ")", sep="")
m3 <- paste(x$model, ": Negative Log-Likelihood\nGradient wrt ", p2name, " (", x$data.name[1], ")", sep="")
}
} else m1 <- m2 <- m3 <- a$main
if(is.element(type,c("likelihood","both"))) {
image(p1, p2, X, main=m1, xlab=xl, ylab=yl, zlim=zl.lh)
points(minloc[1], minloc[2], pch="x", col="darkorange")
# image.plot(p1, p2, X, zlim=zl.lh, horizontal=hz, legend.only=TRUE, legend.mar=legmar)
} # end of if likelihood or both stmts.
if(is.element(type, c("gradient","both"))) {
if(is.element(1, which.gr)) {
image(p1, p2, x$nllh.gr1, main=m2, xlab=xl, ylab=yl, zlim=zl.gr1)
# image.plot(p1, p2, x$nllh.gr1, zlim=zl.gr1, horizontal=hz, legend.only=TRUE, legend.mar=legmar)
}
if(is.element(2, which.gr)) {
image(p1, p2, x$nllh.gr2, main=m3, xlab=xl, ylab=yl, zlim=zl.gr2)
# image.plot(p1, p2, x$nllh.gr2, zlim=zl.gr2, horizontal=hz, legend.only=TRUE, legend.mar=legmar)
}
} # end of if type is gradient or both stmts.
}# end of if one or two parameters varied stmts.
invisible()
} # end of 'plot.grlevdTracer' function.
check.constant <- function(x) {
# internal function to determine whether a parameter is to be modeled
# as constant (TRUE) or not (FALSE).
if(is.formula(x)) {
if(length(as.character(x))==2 && as.character(x)[2] == "1") return(TRUE)
else return(FALSE)
} else {
if(length(x)==1) return(TRUE)
else if(all(x == x[1])) return(TRUE)
else return(FALSE)
}
} # end of 'check.constant' function.
setup.design <- function(x, ...) {
UseMethod("setup.design", x)
} # end of 'setup.design' function.
setup.design.default <- function(x, ..., data, n, const, dname) {
if(missing(data)) data <- NULL
if(missing(const)) const <- check.constant(x)
if(const || is.vector(x)) return(matrix(1, n, 1))
else if(is.formula(x)) {
if(is.null(data)) return(model.matrix(x))
else return(model.matrix(x, data))
} else stop(paste("setup.design: invalid ", dname, " parameter.", sep=""))
} # end of 'setup.design' function.
setup.design.fevd <- function(x, ...) {
out <- list()
n <- x$n
mods <- x$par.models
dname <- x$data.name[2]
data <- datagrabber(x, response=FALSE)
out$X.u <- setup.design(mods$threshold, data=data, n=n, dname="threshold")
out$X.loc <- setup.design(mods$location, data=data, n=n, dname="location")
out$X.sc <- setup.design(mods$scale, data=data, n=n, dname="scale")
out$X.sh <- setup.design(mods$shape, data=data, n=n, dname="shape")
return(out)
} # end of 'setup.design.fevd' function.
oevd <- function(p, o, des, x, data=NULL, u=NULL, span, npy, phi=FALSE, blocks=NULL) {
##
## Function to be optimized by 'optim'.
## Takes the parameter design matrices and
## calculates the resulting vector of parameters,
## then calls 'levd' to get the negative log-likelihood.
type <- o$type
n <- length(x)
w <- o$weights
if(length(w) == 1) w <- rep(w, n)
if(!is.element(type, c("GP","Beta","Pareto","Exponential"))) {
X.loc <- des$X.loc
nloc <- ncol(X.loc)
loc <- rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc)
} else {
nloc <- 0
loc <- 0
}
X.sc <- des$X.sc
nsc <- ncol(X.sc)
if(phi) scale <- exp(rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc))
else scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc)
if(!is.element(o$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 {
nsh <- 0
shape <- 0
}
if(!is.null(blocks)) { # could check for const.loc, etc. and just do sum(p[1:nloc]*blocks$designs$X.param[1,])
if(is.element(type, c("PP"))) {
blocks$location <- rowSums(matrix(p[1:nloc], blocks$nBlocks, nloc, byrow=TRUE)*blocks$designs$X.loc)
}
if(phi) blocks$scale <- exp(rowSums(matrix(p[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE)*blocks$designs$X.sc))
else blocks$scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE)*blocks$designs$X.sc)
blocks$shape <- rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], blocks$nBlocks, nsh, byrow=TRUE)*blocks$designs$X.sh)
}
if(missing(span)) res <- levd(x=x, threshold=u, location=loc, scale=scale, shape=shape, type=type, npy=npy, infval=1e10, weights=w, blocks=blocks)
else res <- levd(x=x, threshold=u, location=loc, scale=scale, shape=shape, type=type, span=span, npy=npy, infval=1e10, weights=w, blocks=blocks)
return(res)
} # end of 'oevd' function.
oevdgen <- function(p, o, des, x, data = NULL, u = NULL, span, npy, phi = FALSE, priorFun, priorParams = NULL, blocks = NULL) {
##
## Function to be optimized by 'optim'.
## Takes the parameter design matrices and
## calculates the resulting vector of parameters,
## then calls 'levd' to get the negative log-likelihood.
type <- o$type
n <- length(x)
w <- o$weights
if(length(w) == 1) w <- rep(w, n)
if(!is.element(type, c("GP","Beta","Pareto","Exponential"))) {
X.loc <- des$X.loc
nloc <- ncol(X.loc)
loc <- rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE) * X.loc)
if(nloc == 1) lname <- "location"
else lname <- paste("mu", 0:(nloc - 1), sep = "")
} else {
nloc <- 0
loc <- 0
lname <- NULL
}
X.sc <- des$X.sc
nsc <- ncol(X.sc)
if(phi) scale <- exp(rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE) * X.sc))
else scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE) * X.sc)
if(nsc == 1) scname <- "scale"
else scname <- paste("phi", 0:(nsc - 1), sep="")
if(!is.null(des$X.sh)) {
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)
if(nsh == 1) shname <- "shape"
else shname <- paste("xi", 0:(nsh - 1), sep = "")
} else {
nsh <- 0
shape <- 0
shname <- NULL
}
if(!is.null(blocks)) {
blocks$location <- rowSums(matrix(p[1:nloc], blocks$nBlocks, nloc, byrow=TRUE) * blocks$designs$X.loc)
if(phi) blocks$scale <- exp(rowSums(matrix(p[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE) * blocks$designs$X.sc))
else blocks$scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE) * blocks$designs$X.sc)
blocks$shape <- rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], blocks$nBlocks, nsh, byrow=TRUE) * blocks$designs$X.sh)
}
if(missing(span)) res <- levd(x=x, threshold=u, location=loc, scale=scale, shape=shape, type=type, npy=npy, infval=1e10, weights=w, blocks=blocks)
else res <- levd(x=x, threshold=u, location=loc, scale=scale, shape=shape, type=type, span=span, npy=npy, infval=1e10, weights=w, blocks=blocks)
# newp <- cbind(loc, scale, shape)
# colnames(newp) <- c(lname, scname, shname)
# res2 <- do.call(priorFun, c(list(x = newp), priorParams))
res2 <- do.call(priorFun, c(list(x = p), priorParams))
return(res + res2)
} # end of 'oevdgen' function.
shapePriorBeta <- function( x, p, q ) {
xi <- x[ length(x) ]
if( xi >= -0.5 & xi <= 0.5 ) {
-( 0.5 + xi )^( q-1 ) * ( 0.5 - xi )^( p-1 ) / beta(q,p)
} else {
1e16
} # end of if else 'xi between -0.5 and 0.5' stmts.
} # end of 'shapePriorBeta' function.
shapePriorBetaOld <- function(x, p, q) {
##
## Function to calculate the GMLE prior for the shape parameter.
##
# if(!is.element("shape", colnames(x))) stop("shapePriorBeta: invalid model for this prior penalty function.")
# xi <- 0.5 + x[,"shape"]
xi <- x[ length(x) ]
# res <- -sum(dbeta(x=xi, shape1=q, shape2=p, log=TRUE), na.rm=TRUE)
res <- - dbeta(x = xi, shape1 = q, shape2 = p, log = TRUE)
if(!is.finite(res)) res <- sign(res) * 1e16
return(res)
} # end of 'shapePriorBetaOld' function.
initializer <- function(x, model, threshold, ...) {
UseMethod("initializer", x)
} # end of 'initializer' function.
initializer.lmoments <- function(x, model, threshold, ..., npy, blocks=NULL) {
if(is.element(model, c("gev", "weibull", "frechet", "gumbel"))) {
lambda <- Lmoments(x)
tau3 <- lambda[3] / lambda[2]
co <- 2 / (3 + tau3) - log(2) / log(3)
kappa <- 7.8590 * co + 2.9554 * co^2
g <- gamma(1 + kappa)
sigma <- (lambda[2] * kappa)/((1 - 2^(-kappa)) * g)
mu <- lambda[1] - (sigma / kappa) * (1 - g)
xi <- -kappa
res <- c(mu, sigma, xi)
names(res) <- c("location", "scale", "shape")
} else if(is.element(model, c("gp", "exponential", "beta", "pareto", "pp"))) {
n <- length(x)
u <- threshold
if(length(u)==1) u <- rep(u, n)
id <- x > u
lambda <- Lmoments(x[id] - u[id])
tau2 <- lambda[2] / lambda[1]
sigma <- lambda[1] * (1 / tau2 - 1)
kappa <- 1 / tau2 - 2
xi <- -kappa
if(model=="pp") {
if(is.null(blocks)) {
rate <- npy * sum(id)/n
} else {
rate <- sum(id)/(blocks$nBlocks * mean(blocks$weights))
}
sigma <- exp(log(sigma) + xi*log(rate))
mu <- mean(threshold) - (sigma/xi)*(rate^(-xi) - 1)
res <- c(mu, sigma, xi, rate)
names(res) <- c("location", "scale", "shape", "rate")
} else if(model != "exponential") {
res <- c(sigma, xi)
names(res) <- c("scale", "shape")
} else {
res <- sigma
names(res) <- c("scale")
}
} # end of 'which models to estimate' stmts.
return(res)
} # end of 'initializer.lmoments' function.
initializer.moms <- function(x, model, threshold, ..., npy, blocks=NULL) {
if(is.element(model, c("gev","gumbel","frechet","weibull"))) {
m <- mean(x)
s <- sqrt(6*var(x))/pi
mu <- m - 0.57722*s
sigma <- log(s)
if(model=="weibull") xi <- -1e-8
else if(is.element(model,c("gev","frechet"))) xi <- 1e-8
else xi <- 0
res <- c(mu, sigma, xi)
names(res) <- c("location", "scale", "shape")
} else if(is.element(model, c("pp", "gp", "exponential", "beta", "pareto"))) {
n <- length(x)
u <- threshold
if(length(u)==1) u <- rep(u,n)
id <- x > u
x.u <- x[id] - u[id]
m <- mean(x.u)
sigma <- sqrt(var(x.u))
if(model=="beta") xi <- -1e-8
else if(is.element(model, c("pp","gp","beta"))) xi <- 1e-8
else xi <- 0
if(model=="pp") {
if(is.null(blocks)) {
rate <- npy * sum(id)/n
} else {
rate <- sum(id)/(blocks$nBlocks * mean(blocks$weights))
}
sigma <- exp(log(sigma) + xi*rate)
mu <- mean(u) - (sigma/xi)*(rate^(-xi) - 1)
res <- c(mu, sigma, xi)
names(res) <- c("location", "scale", "shape")
} else if(model != "exponential") {
res <- c(sigma, xi)
names(res) <- c("scale", "shape")
} else {
res <- sigma
names(res) <- "scale"
}
} else stop("initializer: invalid model argument.")
return(res)
} # end of 'initializer.moms' function.
initializer.mle <- function(x, model, threshold, ..., data=NULL, u, u.fun, loc.fun, sc.fun, sh.fun, use.phi, type, span, time.units, period.basis, blocks=NULL) {
attributes(x) <- NULL
if(is.null(data)) {
if(missing(span)) res <- fevd(x=x, threshold=u, use.phi=use.phi, type=type, time.units=time.units, period.basis=period.basis, blocks=blocks)
else res <- fevd(x=x, threshold=u, use.phi=use.phi, type=type, span=span, time.units=time.units, period.basis=period.basis, blocks=blocks)
} else {
if(missing(span)) res <- fevd(x=x, data=data, threshold=u, threshold.fun=u.fun, location.fun=loc.fun, scale.fun=sc.fun, shape.fun=sh.fun,
use.phi=use.phi, type=type, time.units=time.units, period.basis=period.basis, blocks=blocks)
else res <- fevd(x=x, data=data, threshold=u, threshold.fun=u.fun, location.fun=loc.fun, scale.fun=sc.fun, shape.fun=sh.fun,
use.phi=use.phi, type=type, span=span, time.units=time.units, period.basis=period.basis, blocks=blocks)
}
return(res$results$par)
} # end of 'initializer.bayesian' function.
will.accept <- function( ll.i, prior.i, ll.star, prior.star, log=TRUE) {
##
## Function to determine whether to accept a proposed parameter value or not.
##
## Arguments:
##
## 'll.i' numeric giving the log-likelihood value for the i-th step parameter combination.
## 'prior.i' numeric giving the (log) prior probability for the i-th step parameter combination.
## 'll.star', 'prior.star' are the log-likelihood and (log) prior values for the proposed parameter vector.
##
## Details: Calculates the ratio R=(L(theta*;data)p(theta*))/(L(theta.i;data)p(theta.i)), and determines
## whether to accept or reject the proposed parameters with probability alpha = min(1, R).
##
## Value: list with components:
## 'accept' a logical saying whether to accept the proposed parameters (if TRUE) or not (if FALSE),
## 'ratio' numeric giving the value of R (see Details above),
## 'alpha' numeric giving the value of 'alpha' (see Details above),
## 'll.star' numeric giving the log-likelihood of the proposed parameter combination,
## 'prior.star' the prior probability for the proposed parameter combination.
##
out <- list()
ratio <- (ll.star + prior.star) - (ll.i + prior.i)
if(!log) ratio <- exp( ratio)
if( is.na( ratio)) {
out$accept <- FALSE
out$ratio <- out$alpha <- NA
return( out)
}
if(!log) alpha <- min( 1, ratio)
else alpha <- min(0, ratio)
if(!log) {
if( runif( 1) < alpha) accept <- TRUE
else accept <- FALSE
} else {
if(log(runif(1)) < alpha) accept <- TRUE
else accept <- FALSE
}
out$accept <- accept
out$ratio <- ratio
out$alpha <- alpha
return( out)
} # end of 'will.accept' function.
fevdPriorDefault <- function(theta, m, v, log=TRUE) {
np <- length(theta)
dfun <- function(th) dnorm(th[1], mean=th[2], sd=th[3], log=log)
if(length(m)==1) m <- rep(m, np)
else if(length(m) != np) stop("fevdPriorDefault: m must have length 1 or same as number of parameters.")
if(length(v)==1) v <- rep(v, np)
else if(length(v) != np) stop("fevdPriorDefault: v must have length 1 or same as number of parameters.")
th <- cbind(theta, m, v)
res <- apply(th, 1, dfun)
if(log) return(sum(res))
else return(prod(res))
} # end of 'fevdPriorDefault' function.
fevdProposalDefault <- function(p, ind, ...) {
a <- list(...)
if(!is.null(a$mean)) {
#if(length(a$mean) > 1) m <- a$mean[ind]
#else m <- a$mean
m <- a$mean
} else m <- 0
if(!is.null(a$sd)) {
# if(length(a$sd) > 1) s <- a$sd[ind]
# else s <- a$sd
s <- a$sd
} else s <- 0.1
return(p + rnorm(length(p), mean=m, sd=s))
} # end of 'fevdProposalDefault' function.
fpois <- function(x, na.action = na.fail, ...) {
UseMethod("fpois", x)
} # end of 'fpois' function.
fpois.default <- function(x, na.action = na.fail, ...) {
theCall <- match.call()
dname <- deparse(substitute(x))
x <- na.action(x)
n <- length(x)
m <- mean(x)
v <- var(x)
PARAMETER <- c(m, v, n - 1)
names(PARAMETER) <- c("mean", "variance", "degrees of freedom")
CRITVAL <- (n - 1) * v / m
names(CRITVAL) <- "Chi-square(n - 1)"
PVAL <- pchisq(CRITVAL, df = n - 1, lower.tail = FALSE)
structure(list(statistic = CRITVAL, parameter = PARAMETER,
alternative = "greater", p.value = PVAL, method = "Test for Equality of (Poisson) Mean and Variance",
data.name = dname, call = theCall), class = "htest")
} # end of 'fpois.default' function.
fpois.data.frame <- function(x, na.action = na.fail, ..., which.col = 1) {
theCall <- match.call()
dname <- deparse(substitute(x))
if(length(which.col) > 1) stop("fpois: invalid which.col argument. Length must be one. To include covariates, use glm with poisson() family link.")
x <- x[, which.col]
x <- na.action(x)
n <- length(x)
m <- mean(x)
v <- var(x)
PARAMETER <- c(m, v, n - 1)
names(PARAMETER) <- c("mean", "variance", "degrees of freedom")
CRITVAL <- (n - 1) * v / m
names(CRITVAL) <- "Chi-square(n - 1)"
PVAL <- pchisq(CRITVAL, df = n - 1, lower.tail = FALSE)
structure(list(statistic = CRITVAL, parameter = PARAMETER,
alternative = "greater", p.value = PVAL, method = "Test for Equality of (Poisson) Mean and Variance",
data.name = dname, call = theCall), class = "htest")
} # end of 'fpois.data.frame' function.
fpois.matrix <- function(x, na.action = na.fail, ..., which.col = 1) {
theCall <- match.call()
dname <- deparse(substitute(x))
if(length(which.col) > 1) stop("fpois: invalid which.col argument. Length must be one. To include covariates, use glm with poisson() family link.")
x <- x[, which.col]
x <- na.action(x)
n <- length(x)
m <- mean(x)
v <- var(x)
PARAMETER <- c(m, v, n - 1)
names(PARAMETER) <- c("mean", "variance", "degrees of freedom")
CRITVAL <- (n - 1) * v / m
names(CRITVAL) <- "Chi-square(n - 1)"
PVAL <- pchisq(CRITVAL, df = n - 1, lower.tail = FALSE)
structure(list(statistic = CRITVAL, parameter = PARAMETER,
alternative = "greater", p.value = PVAL, method = "Test for Equality of (Poisson) Mean and Variance",
data.name = dname, call = theCall), class = "htest")
} # end of 'fpois.matrix' function.
fpois.list <- function(x, na.action = na.fail, ..., which.component = 1) {
theCall <- match.call()
dname <- deparse(substitute(x))
if(length(which.component) > 1) stop("fpois: invalid which.component argument. Length must be one. To include covariates, use glm with poisson() family link.")
x <- x[[which.component]]
x <- na.action(x)
n <- length(x)
m <- mean(x)
v <- var(x)
PARAMETER <- c(m, v, n - 1)
names(PARAMETER) <- c("mean", "variance", "degrees of freedom")
CRITVAL <- (n - 1) * v / m
names(CRITVAL) <- "Chi-square(n - 1)"
PVAL <- pchisq(CRITVAL, df = n - 1, lower.tail = FALSE)
structure(list(statistic = CRITVAL, parameter = PARAMETER,
alternative = "greater", p.value = PVAL, method = "Test for Equality of (Poisson) Mean and Variance",
data.name = dname, call = theCall), class = "htest")
} # end of 'fpois.list' 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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions fpois.list fpois.matrix fpois.data.frame fpois.default fpois fevdProposalDefault fevdPriorDefault will.accept initializer.mle initializer.moms initializer.lmoments initializer shapePriorBetaOld shapePriorBeta oevdgen oevd setup.design.fevd setup.design.default setup.design check.constant plot.grlevdTracer grlevdTracer grlevd print.fevd summary.fevd.mle summary.fevd.bayesian summary.fevd.lmoments summary.fevd distill.fevd.mle distill.fevd.bayesian distill.fevd.lmoments distill.fevd parcov.fevd levd fevd
Documented in check.constant distill.fevd distill.fevd.bayesian distill.fevd.lmoments distill.fevd.mle fevd fevdPriorDefault fevdProposalDefault fpois fpois.data.frame fpois.default fpois.list fpois.matrix grlevd grlevdTracer initializer initializer.lmoments initializer.mle initializer.moms levd oevd oevdgen parcov.fevd plot.grlevdTracer print.fevd setup.design setup.design.default setup.design.fevd shapePriorBeta shapePriorBetaOld summary.fevd summary.fevd.bayesian summary.fevd.lmoments summary.fevd.mle will.accept
fevd <- function(x, data, threshold=NULL, threshold.fun=~1, location.fun=~1, scale.fun=~1, shape.fun=~1, use.phi=FALSE,
type=c("GEV", "GP", "PP", "Gumbel", "Exponential"),
method=c("MLE", "GMLE", "Bayesian", "Lmoments"), initial=NULL, span, units=NULL, time.units="days", period.basis="year",
na.action=na.fail, optim.args=NULL, priorFun=NULL, priorParams=NULL, proposalFun=NULL, proposalParams=NULL,
iter=9999, weights=1, blocks=NULL, verbose=FALSE) {
if(verbose) begin.tiid <- Sys.time()
out <- list()
inout <- list() # object to return initial parameter values with likelihood.
out$call <- match.call()
if(!missing(data)) {
out$data.name <- c(deparse(substitute(x)), deparse(substitute(data)))
# cov.pointer <- as.character(substitute(data))
} else {
out$data.name <- c(deparse(substitute(x)), "")
# cov.pointer <- NULL
}
# data.pointer <- as.character(substitute(x))
# if( length(data.pointer) > 1) out$x <- x
# else out$data.pointer <- data.pointer
# if(length(cov.pointer) > 1) out$cov.data <- data
# else out$cov.pointer <- cov.pointer
type <- match.arg(type)
method <- match.arg(method)
out$weights <- weights
if(!missing(data)) {
if(is.element(out$data.name[1], colnames(data))) {
out$in.data <- TRUE
wc <- out$data.name[1] == colnames(data)
x <- c(data[,wc])
x.fun <- ifelse(out$data.name[1] == "substitute(x)", deparse(x), out$data.name[1])
x.fun <- formula(paste(x.fun, "~ 1"))
out$x.fun <- x.fun
} else if(is.formula(x)) {
out$in.data <- TRUE
x.fun <- x
out$x.fun <- x.fun
x <- model.response(model.frame(x.fun, data=data))
} else out$in.data <- FALSE
if(length(x) != nrow(data)) stop("fevd: data must have same number of rows as the length of x.")
if(!identical(weights, 1) && length(x) != length(weights)) # CJP
stop("fevd: weights should be the same length as x.")
out$missing.values <- is.na( x )
tmp <- cbind(x, data)
tmp <- na.action(tmp)
x <- tmp[,1]
data <- tmp[,-1, drop=FALSE]
} else {
out$missing.values <- is.na( x )
if(is.formula(x)) stop("fevd: Must provide data argument if you supply a formula to the x argument.")
x <- na.action(x)
out$in.data <- FALSE
}
if(!out$in.data) {
data.pointer <- as.character(substitute(x))
if(length(data.pointer) > 1) out$x <- x
else out$data.pointer <- data.pointer
}
if(!is.null(blocks)){ # CJP
# test that blocks is of correct type
if(type == "PP"){
if(!is.element('nBlocks', names(blocks))){
if(is.element('data', names(blocks))) {
blocks$nBlocks <- nrow(blocks$data)
} else stop("fevd: When supplying blocks, must provide 'blocks$nBlocks' if 'blocks$data' is not provided.")
}
if(!is.element('weights', names(blocks))) blocks$weights <- 1
if(!is.element('proportionMissing', names(blocks))) blocks$proportionMissing <- 0
if(!is.element('threshold', names(blocks)) && !is.null(threshold)) {
if(length(threshold) == 1) {
blocks$threshold <- threshold
} else {
stop("fevd: No blocks$threshold specified and threshold is not constant. User must supply the threshold for each block via blocks$threshold.")
}
}
} else {
warning("fevd: Blocks are used only for type 'PP'; ignoring blocks argument.")
blocks <- NULL # code will now operate as if user did not supply 'blocks'
}
}
out$x <- x
if(!missing(data)) out$cov.data <- data
if( method == "MLE" && !is.null(priorFun) ) method <- "GMLE"
else if( method == "GMLE" && is.null(priorFun) ) {
if(shape.fun != ~1) stop("fevd: must supply a prior function for GMLE method when shape parameter varies.")
if(is.element(type, c("GEV","GP","PP"))) {
priorFun <- "shapePriorBeta"
if(is.null(priorParams)) priorParams <- list(q=6, p=9)
} else {
warning("fevd: Using method MLE. No default for specified arguments.")
method <- "MLE"
}
} # end of if else method is GMLE but no prior function input stmts.
if(method=="GMLE") {
out$priorFun <- priorFun
out$priorParams <- priorParams
}
out$method <- method
method <- tolower(method)
out$type <- type
type <- tolower(type)
out$period.basis <- period.basis
out$optim.args <- optim.args
out$units <- units
if(method=="bayesian" && missing(use.phi)) {
use.phi <- TRUE
if(verbose) cat("\n", "Setting use.phi argument to TRUE for greater stability in estimation (default with Bayesian method). Use use.phi=FALSE if you prefer that.\n")
}
# if(is.element(type, c("gp","pp","exponential","beta","pareto"))) out$npy <- npy
out$par.models <- list(threshold=threshold.fun, location=location.fun, scale=scale.fun, log.scale=use.phi, shape=shape.fun,
term.names=list(threshold=all.vars(threshold.fun), location=all.vars(location.fun),
scale=all.vars(scale.fun), shape=all.vars(shape.fun)))
pars <- list()
if(is.element(type, c("gp","pp","exponential","beta","pareto"))) {
const.thresh <- check.constant(threshold.fun) & check.constant(threshold)
out$const.thresh <- const.thresh
}
if(is.element(type, c("gev","pp","gumbel","weibull","frechet"))) {
const.loc <- check.constant(location.fun)
out$const.loc <- const.loc
}
const.scale <- check.constant(scale.fun)
out$const.scale <- const.scale
const.shape <- check.constant(shape.fun)
out$const.shape <- const.shape
if(is.element(type, c("pp", "gp", "exponential", "beta", "pareto"))) {
if(missing(span)) {
if(is.null(blocks)) {
tiden <- attributes(x)$times
n <- length(x)
if(is.null(tiden)) {
tiden <- 1:n
start <- 1
end <- n
span <- end - start
} else {
start <- tiden[1]
end <- tiden[n]
span <- as.numeric(difftime(as.POSIXlt(tiden)[n], as.POSIXlt(tiden)[1], units=time.units))
}
} else {
span <- blocks$nBlocks
}
} # end of if 'missing span' stmts.
if(time.units=="days") npy <- 365.25
else if(time.units=="months") npy <- 12
else if(time.units=="years") npy <- 1
else if(time.units == "hours") npy <- 24 * 365.25
else if(time.units == "minutes") npy <- 60 * 24 * 365.25
else if(time.units == "seconds") npy <- 60 * 60 * 24 * 365.25
else {
tmp.units <- unlist(strsplit(time.units, split="/"))
if(length(tmp.units) != 2) stop("fevd: invalid time.units argument.")
numper <- as.numeric(tmp.units[1])
if(is.na(numper)) stop("fevd: invalid time.units argument.")
pertiid <- tmp.units[2]
if(!is.element(pertiid, c("day","month","year","hour","minute","second"))) stop("fevd: invalid time.units argument.")
if(pertiid=="year") npy <- numper
else if(pertiid=="month") npy <- numper*12
else if(pertiid=="day") npy <- numper*365.25
else if(pertiid == "hour") npy <- numper * 24 * 365.25
else if(pertiid == "minute") npy <- numper * 60 * 24 * 365.25
else if(pertiid == "second") npy <- numper * 60 * 60 * 24 * 365.25
} # end of if complicated time units stmts.
if(!is.null(blocks)) span <- span*npy # CJP - so is on the same scale as span w/o blocks
out$time.units <- time.units
out$span <- span/npy
out$npy <- npy
if(verbose)
cat("\n", "Data span ", span/npy, "years", "\n")
} else npy <- NULL # end of if 'type is of the POT variety' stmts.
# if(!missing(data)) {
# if(is.element(out$data.name[1], colnames(data))) {
# out$in.data <- TRUE
# wc <- out$data.name[1] == colnames(data)
# x <- c(data[,wc])
# x.fun <- ifelse(out$data.name[1] == "substitute(x)", deparse(x), out$data.name[1])
# x.fun <- formula(paste(x.fun, "~ 1"))
# out$x.fun <- x.fun
# } else if(is.formula(x)) {
# out$in.data <- TRUE
# x.fun <- x
# out$x.fun <- x.fun
# x <- model.response(model.frame(x.fun, data=data))
# } else out$in.data <- FALSE
# if(length(x) != nrow(data)) stop("fevd: data must have same number of rows as the length of x.")
# tmp <- cbind(x, data)
# tmp <- na.action(tmp)
# x <- tmp[,1]
# data <- tmp[,-1]
# } else {
# if(is.formula(x)) stop("fevd: Must provide data argument if you supply a formula to the x argument.")
# x <- na.action(x)
# out$in.data <- FALSE
# }
#
# if(!out$in.data) {
# data.pointer <- as.character(substitute(x))
# if(length(data.pointer) > 1) out$x <- x
# else out$data.pointer <- data.pointer
# }
n <- length(x)
out$n <- n
out$na.action <- deparse(substitute(na.action))
if(!is.null(initial)) {
if(!is.list(initial)) stop("fevd: initial must be NULL or a list object.")
find.init <- FALSE
if(is.null(initial$location) && is.element(type, c("gev","pp","gumbel","weibull","frechet"))) find.init <- TRUE
if(use.phi && is.null(initial$log.scale)) find.init <- TRUE
if(!use.phi && is.null(initial$scale)) find.init <- TRUE
if(!is.element(type, c("gumbel","exponential")) && is.null(initial$shape)) find.init <- TRUE
} else {
initial <- list()
find.init <- TRUE
}
if(method != "lmoments") {
if(verbose) cat("Setting up parameter model design matrices.\n")
designs <- list()
if(!missing(data)) {
if(is.element(type, c("gp", "pp", "exponential", "beta", "pareto"))) X.u <- setup.design(x=threshold.fun, data=data, n=n, dname="threshold.fun")
if(is.element(type, c("gev", "pp", "gumbel", "weibull", "frechet"))) {
X.loc <- setup.design(x=location.fun, data=data, n=n, const=const.loc, dname="location.fun")
designs$X.loc <- X.loc
}
X.sc <- setup.design(x=scale.fun, data=data, n=n, const=const.scale, dname="scale.fun")
designs$X.sc <- X.sc
if(!is.element(type,c("gumbel","exponential"))) {
X.sh <- setup.design(x=shape.fun, data=data, n=n, const=const.shape, dname="shape.fun")
designs$X.sh <- X.sh
}
} else {
if(is.element(type, c("gp", "pp", "exponential", "beta", "pareto"))) X.u <- setup.design(x=threshold.fun, n=n, dname="threshold.fun")
if(is.element(type, c("gev", "pp", "gumbel", "weibull", "frechet"))) {
X.loc <- setup.design(x=location.fun, n=n, const=const.loc, dname="location.fun")
designs$X.loc <- X.loc
}
X.sc <- setup.design(x=scale.fun, n=n, const=const.scale, dname="scale.fun")
designs$X.sc <- X.sc
if(!is.element(type,c("gumbel","exponential"))) {
X.sh <- setup.design(x=shape.fun, n=n, const=const.shape, dname="shape.fun")
designs$X.sh <- X.sh
}
} # end of if else missing data stmts.
if(!is.null(blocks)){ # CJP
blocks$designs <- list()
if(is.element('data', names(blocks))){
blocks$X.u <- setup.design(x = threshold.fun, data=blocks$data, n=blocks$nBlocks, dname="threshold.fun")
blocks$designs$X.loc <- setup.design(x=location.fun, data=blocks$data, n=blocks$nBlocks, const=const.loc, dname="location.fun")
blocks$designs$X.sc <- setup.design(x=scale.fun, data=blocks$data, n=blocks$nBlocks, const=const.scale, dname="scale.fun")
blocks$designs$X.sh <- setup.design(x=shape.fun, data=blocks$data, n=blocks$nBlocks, const=const.shape, dname="shape.fun")
} else {
blocks$X.u <- setup.design(x=threshold.fun, n=blocks$nBlocks, dname="threshold.fun")
blocks$designs$X.loc <- setup.design(x=location.fun, n=blocks$nBlocks, const=const.loc, dname="location.fun")
blocks$designs$X.sc <- setup.design(x=scale.fun, n=blocks$nBlocks, const=const.scale, dname="scale.fun")
blocks$designs$X.sh <- setup.design(x=shape.fun, n=blocks$nBlocks, const=const.shape, dname="shape.fun")
}
}
if(verbose) cat("Parameter model design matrices set up.\n")
} # end of if method is not Lmoments to set up design matrices stmts.
if(is.element(type, c("gp", "pp", "exponential", "beta", "pareto"))) {
if(method != "lmoments") threshold <- rowSums(matrix(threshold, n, ncol(X.u), byrow=TRUE) * X.u)
if(!is.null(blocks)) # CJP
blocks$threshold <- rowSums(matrix(blocks$threshold, blocks$nBlocks, ncol(blocks$X.u), byrow=TRUE) * blocks$X.u)
excess.id <- x > threshold
if(all(threshold==threshold[1])) out$threshold <- threshold[1]
else out$threshold <- threshold
out$rate <- mean(excess.id)
} # end of if '!missing(threshold)' stmts.
out$blocks <- blocks # CJP
if(method=="lmoments" || find.init) {
if(method=="lmoments") {
if(verbose) cat("Beginning estimation procedure.\n")
is.constant <- unlist(lapply(list(u=threshold, loc=location.fun, scale=scale.fun, sh=shape.fun), check.constant))
if(!all(is.constant)) warning("fevd: For method Lmoments, this function does not handle covariates in the parameters. Fitting w/o covariates.")
if(!is.element(type, c("gev","gp"))) stop("fevd: currently, Lmoments are only handled for estimation of GEV and GP distribution parameters.")
}
xtemp <- x
class(xtemp) <- "lmoments"
ipars1 <- try(initializer(xtemp, model=type, threshold=threshold, npy=npy, blocks=blocks), silent = TRUE)
if( !is( ipars1, "try-error" ) ) { # Eric -- 8/6/13
if(ipars1["scale"] <= 0) ipars1["scale"] <- 1e-8
if(method=="lmoments") {
out$results <- ipars1
class(out) <- "fevd"
return(out)
}
} else {
ipars1 <- NULL
if(method == "lmoments") stop("fevd: Sorry, could not find L-moments estimates.")
}
} # end of if 'Lmoments method or find them anyway as potential initial estimates' stmts.
if((method != "lmoments") && find.init) {
xtemp <- x
class(xtemp) <- "moms"
ipars2 <- try(initializer(xtemp, model=type, threshold=threshold, npy=npy, blocks=blocks), silent = TRUE)
if( !is(ipars2, "try-error" ) ) {
if(ipars2["scale"] <= 0) ipars2["scale"] <- 1e-8
} else ipars2 <- NULL
if(!is.element(type, c("pp", "gp","exponential","beta","pareto","gumbel"))) {
if(!is.null(ipars1)) testLmoments <- levd(x, location=ipars1["location"], scale=ipars1["scale"], shape=ipars1["shape"], type=out$type, npy=npy)
else testLmoments <- Inf
if(!is.null(ipars2)) testMoments <- levd(x, location=ipars2["location"], scale=ipars2["scale"], shape=ipars2["shape"], type=out$type, npy=npy)
else testMoments <- Inf
} else if(type=="pp") { # CJP
if(!is.null(ipars1)) { # Eric -- 8/6/13
if(!is.null(blocks)) {
blocks$location=ipars1["location"]; blocks$scale=ipars1["scale"]; blocks$shape=ipars1["shape"]
}
testLmoments <- levd(x, threshold=threshold, location=ipars1["location"], scale=ipars1["scale"], shape=ipars1["shape"], type=out$type, npy=npy, blocks=blocks)
} else testLmoments <- Inf
# end of if 'ipars1' not NULL stmts.
if(!is.null(ipars2)) { # Eric -- 8/6/13
if(!is.null(blocks)) {
blocks$location=ipars2["location"]; blocks$scale=ipars2["scale"]; blocks$shape=ipars2["shape"]
}
testMoments <- levd(x, threshold=threshold, location=ipars2["location"], scale=ipars2["scale"], shape=ipars2["shape"], type=out$type, npy=npy, blocks=blocks)
} else testMoments <- Inf
# end of if 'ipars2' not NULL stmts.
if(!is.null(blocks)) # CJP
blocks$location <- blocks$scale <- blocks$shape <- NULL
} else if(!is.element(type, c("gumbel","exponential"))) {
if(!is.null(ipars1)) testLmoments <- levd(x, threshold=threshold, scale=ipars1["scale"], shape=ipars1["shape"], type=out$type, npy=npy)
else testLmoments <- Inf
if(!is.null(ipars2)) testMoments <- levd(x, threshold=threshold, scale=ipars2["scale"], shape=ipars2["shape"], type=out$type, npy=npy)
else testMoments <- Inf
} else if(type=="gumbel") {
if(!is.null(ipars1)) testLmoments <- levd(x, location=ipars1["location"], scale=ipars1["scale"], type=out$type, npy=npy)
else testLmoments <- Inf
if(!is.null(ipars2)) testMoments <- levd(x, location=ipars2["location"], scale=ipars2["scale"], type=out$type, npy=npy)
else testMoments <- Inf
} else if(type=="exponential") {
if(!is.null(ipars1)) testLmoments <- levd(x, threshold=threshold, scale=ipars1["scale"], shape=0, type=out$type, npy=npy)
else testLmoments <- Inf
if(!is.null(ipars2)) testMoments <- levd(x, threshold=threshold, scale=ipars2["scale"], shape=0, type=out$type, npy=npy)
else testMoments <- Inf
}
if(is.finite(testLmoments) || is.finite(testMoments)) { # Eric -- 8/6/13
if(testLmoments < testMoments) {
if(is.null(initial$location) && !is.element(type, c("gp","exponential","beta","pareto"))) initial$location <- ipars1["location"]
if(is.null(initial$log.scale) && use.phi) initial$log.scale <- log(ipars1["scale"])
else if(is.null(initial$scale)) initial$scale <- ipars1["scale"]
if(!is.element(type, c("gumbel","exponential")) && is.null(initial$shape)) initial$shape <- ipars1["shape"]
if(verbose) cat("Using Lmoments estimates as initial estimates. Initial value =", testLmoments, "\n")
} else {
if(is.null(initial$location) && !is.element(type, c("gp","exponential","beta","pareto"))) initial$location <- ipars2["location"]
if(is.null(initial$log.scale) && use.phi) initial$log.scale <- log(ipars2["scale"])
else if(is.null(initial$scale)) initial$scale <- ipars2["scale"]
if(!is.element(type, c("gumbel","exponential")) && is.null(initial$shape)) initial$shape <- ipars2["shape"]
if(verbose) cat("Initial estimates found where necessary (not from Lmoments). Initial value =", testMoments, "\n")
}
} else { # Eric -- 8/6/13; if both L-moments and Moments fails, set initial values to mu = 0, sig = 1, xi = 0.01
if(is.null(initial$location) && !is.element(type, c("gp","exponential","beta","pareto"))) initial$location <- 0
if(is.null(initial$log.scale) && use.phi) initial$log.scale <- 0
else if(is.null(initial$scale)) initial$scale <- 1
if(!is.element(type, c("gumbel","exponential")) && is.null(initial$shape)) initial$shape <- 0.01
warning("fevd: L-moments and Moment initial parameter estimates could not be calculated. Using arbitrary starting values.")
}
inout <- list(Lmoments=list(pars=ipars1, likelihood=testLmoments), MOM=list(pars=ipars2, likelihood=testMoments))
} # end of if 'method is not lmoments and need to find initial estimates' stmts.
# Set up initial estimates.
if(!is.null(initial$location)) {
if(ncol(X.loc) != length(initial$location)) {
if((length(initial$location)==1) && ncol(X.loc) > 1) initial$location <- c(initial$location, rep(0, ncol(X.loc)-1))
else stop("fevd: initial parameter estimates must have length 1 or number of parameters present. Incorrect number for location parameter.")
}
if(length(initial$location)==1) names(initial$location) <- "location"
else names(initial$location) <- paste("mu", 0:(ncol(X.loc)-1), sep="")
}
if(use.phi && (ncol(X.sc) != length(initial$log.scale))) {
if((length(initial$log.scale)==1) && ncol(X.sc) > 1) initial$log.scale <- c(initial$log.scale, rep(0, ncol(X.sc)-1))
else stop("fevd: initial parameter estimates must have length 1 or number of parameters present. Incorrect number for log(scale) parameter.")
} else if(!use.phi && (ncol(X.sc) != length(initial$scale))) {
if((length(initial$scale)==1) && ncol(X.sc) > 1) initial$scale <- c(initial$scale, rep(0, ncol(X.sc)-1))
else stop("fevd: initial parameter estimates must have length 1 or number of parameters present. Incorrect number for scale parameter.")
}
if(use.phi) {
if(length(initial$log.scale)==1) names(initial$log.scale) <- "log.scale"
else names(initial$log.scale) <- paste("phi", 0:(ncol(X.sc)-1), sep="")
} else {
if(length(initial$scale)==1) names(initial$scale) <- "scale"
else names(initial$scale) <- paste("sigma", 0:(ncol(X.sc)-1), sep="")
}
if(!is.element(type, c("gumbel","exponential"))) {
if(ncol(X.sh) != length(initial$shape)) {
if((length(initial$shape)==1) && ncol(X.sh) > 1) initial$shape <- c(initial$shape, rep(0, ncol(X.sh)-1))
else stop("fevd: initial parameter estimates must have length 1 or number of parameters present. Incorrect number for shape parameter.")
}
if(length(initial$shape)==1) names(initial$shape) <- "shape"
else names(initial$shape) <- paste("xi", 0:(ncol(X.sh)-1), sep="")
} # end of if type is not light tailed stmts.
if(is.element(method, c("mle","gmle"))) {
if(use.phi) init.pars <- c(initial$location, initial$log.scale, initial$shape)
else init.pars <- c(initial$location, initial$scale, initial$shape)
if(type == "exponential" && const.scale) {
if( method == "gmle" ) warning( "Method MLE used." )
res <- list()
excess.id <- x > threshold
mle <- mean(x[excess.id] - threshold[excess.id])
names(mle) <- "scale"
res$par <- mle
k <- sum(excess.id)
res$n <- k
res$value <- k * ( log(mle) + 1 )
} else {
if(!is.null(a <- optim.args)) {
anam <- names(a)
if(!is.element("gr",anam)) {
if(method=="mle") opt.gr <- grlevd
else opt.gr <- NULL
} else opt.gr <- a$gr
if(is.null(a$method) && use.phi) opt.method <- ifelse(is.element(type,c("gev","gp","pp", "gumbel")), "BFGS", "L-BFGS-B")
else opt.method <- a$method
if(is.element(type, c("weibull","beta","frechet","pareto"))) opt.method <- "L-BFGS-B"
if(is.element(opt.method, c("L-BFGS-B","Brent"))) {
if(is.null(a$lower)) {
if(!is.element(type,c("frechet","pareto"))) opt.lower <- -Inf
else opt.lower <- c(rep(-Inf, length(init.pars)-1), 0)
} else {
if(is.element(type,c("frechet","pareto"))) opt.lower <- c(a$lower[1:(length(init.pars)-1)], 0)
else opt.lower <- a$lower
}
if(is.null(a$upper)) {
if(!is.element(type,c("weibull","beta"))) opt.upper <- Inf
else opt.upper <- c(rep(Inf, length(init.pars)-1), 0)
} else {
if(is.element(type,c("weibull","beta"))) opt.upper <- c(a$upper[1:(length(init.pars)-1)], 0)
else opt.upper <- a$upper
}
} else {
opt.lower <- -Inf
opt.upper <- Inf
} # end of if 'optimization method allows bounds' stmts.
if(is.null(a$control)) opt.control <- list()
else opt.control <- a$control
anam <- names(a$control)
if(!is.element("trace", anam) && verbose) opt.control$trace <- 6
if(is.null(a$hessian)) opt.hessian <- TRUE
else opt.hessian <- a$hessian
} else {
if(method=="mle") opt.gr <- grlevd
else opt.gr <- NULL
if(is.element(type, c("gev","gp","pp","gumbel"))) opt.method <- "BFGS"
else opt.method <- "L-BFGS-B"
if(!is.element(type,c("frechet","pareto"))) opt.lower <- -Inf
else opt.lower <- c(rep(-Inf, length(init.pars)-1), 0)
if(!is.element(type,c("weibull","beta"))) opt.upper <- Inf
else opt.upper <- c(rep(Inf, length(init.pars)-1), 0)
if(verbose) opt.control <- list(trace=6)
else opt.control <- list()
opt.hessian <- TRUE
} # end of if else 'optim.args' is passed stmts.
parnames <- names(init.pars)
out$parnames <- parnames
if(verbose && (method != "lmoments")) {
cat("Initial estimates are:\n")
print(init.pars)
cat("Beginning estimation procedure.\n")
}
if(type=="pp" && find.init) {
# Yet a third attempt to find good starting values--this time using the MLE for the GP parameters.
if(verbose) cat("\n", "First fitting a GP-Poisson model in order to try to get a good initial estimate as PP likelihoods can be very unstable.\n")
look <- out
look$type <- "GP"
des2 <- designs
des2$X.loc <- NULL
if(!missing(data)) resGP <- optim(init.pars[-(1:ncol(X.loc))], oevd, gr=opt.gr, o=look, des=des2, x=x, data=data, u=threshold, npy=npy, phi=use.phi, method=opt.method, lower=opt.lower, upper=opt.upper, control=opt.control, hessian=opt.hessian)
else resGP <- optim(init.pars[-(1:ncol(X.loc))], oevd, gr=opt.gr, o=look, des=des2, x=x, u=threshold, npy=npy, phi=use.phi, method=opt.method, lower=opt.lower, upper=opt.upper, control=opt.control, hessian=opt.hessian)
tmpi <- resGP$par
if(is.null(blocks)) { # CJP
lrate <- npy * mean(x > threshold)
} else {
lrate <- sum(x > threshold)/(blocks$nBlocks * mean(blocks$weights))
}
xi3 <- tmpi[(ncol(X.sc)+1):length(tmpi)]
if(!use.phi) sigma3 <- exp(tmpi[1:ncol(X.sc)] + xi3*log(lrate))
else sigma3 <- tmpi[1:ncol(X.sc)] + xi3*log(lrate)
lp <- lrate^(-xi3) - 1
if(all(is.finite(lp))) mu3 <- mean(threshold) - (sigma3/xi3)*lp # CJP fixed bug - "mean(threshold)" not "threshold"
else mu3 <- mean(x) # CJP: caution here as mean(x) will differ when only exceedances supplied...
nloc <- ncol(X.loc)
if(length(mu3) < nloc) mu3 <- c(mu3, rep(0, nloc-length(mu3)))
else mu3 <- mu3[1]
if(!is.null(blocks)) { # CJP
blocks$location <- rowSums(matrix(mu3, blocks$nBlocks, nloc)*blocks$designs$X.loc)
blocks$scale=rowSums(matrix(sigma3, blocks$nBlocks, ncol(blocks$designs$X.sc))*blocks$designs$X.sc)
blocks$shape=rowSums(matrix(xi3, blocks$nBlocks, ncol(blocks$designs$X.sh))*blocks$designs$X.sh)
}
if(all(is.finite(c(mu3, sigma3, xi3)))) {
testGPmle <- try(levd(x=x, threshold=threshold, location=rowSums(matrix(mu3, n, nloc)*X.loc),
scale=rowSums(matrix(sigma3, n, ncol(X.sc))*X.sc),
shape=rowSums(matrix(xi3, n, ncol(X.sh))*X.sh), type="PP", npy=npy, blocks=blocks), silent = TRUE)
if( is(testGPmle, "try-error" ) ) testGPmle <- Inf
} else testGPmle <- Inf
if(!is.null(blocks)) # CJP
blocks$location <- blocks$scale <- blocks$shape <- NULL
if(is.finite(testLmoments) || is.finite(testMoments) || is.finite(testGPmle)) {
if((testGPmle < testLmoments) && (testGPmle < testMoments)) {
if(verbose) cat("\n", "Changing initial estimates to those based on GP MLEs. They are: \n")
if(use.phi) init.pars <- c(mu3, log(sigma3), xi3)
else init.pars <- c(mu3, sigma3, xi3)
names(init.pars) <- parnames
if(verbose) print(init.pars)
} else if(verbose) cat("\n", "Sticking with originally chosen initial estimates.\n")
} # end of if any one of the three initial estimate attempts was successfull stmts.
inout$PoissonGP <- list(pars=c(mu3, sigma3, xi3), likelihood=testGPmle)
} # end of if 'type is PP, then try to improve initial estimates' stmts.
if(method=="mle") {
if(!missing(data)) {
res <- optim(init.pars, oevd, gr=opt.gr, o=out, des=designs, x=x, data=data, u=threshold,
span=span/npy, npy=npy, phi=use.phi, blocks=blocks, method=opt.method, lower=opt.lower,
upper=opt.upper, control=opt.control, hessian=opt.hessian) # CJP
} else {
res <- optim(init.pars, oevd, gr=opt.gr, o=out, des=designs, x=x, u=threshold,
span=span/npy, npy=npy, phi=use.phi, blocks=blocks, method=opt.method, lower=opt.lower,
upper=opt.upper, control=opt.control, hessian=opt.hessian)
}
} else if(method == "gmle") {
if(!missing(data)) {
res <- optim(init.pars, oevdgen, gr = opt.gr, o = out, des = designs, x = x, data = data,
u = threshold, span = span / npy, npy = npy, phi = use.phi, blocks = blocks,
priorFun = priorFun, priorParams = priorParams, method = opt.method, lower = opt.lower,
upper = opt.upper, control = opt.control, hessian = opt.hessian)
} else {
res <- optim(init.pars, oevdgen, gr = opt.gr, o = out, des = designs, x = x, u = threshold,
span = span / npy, npy = npy, phi = use.phi, blocks = blocks, priorFun = priorFun,
priorParams = priorParams, method = opt.method, lower = opt.lower, upper = opt.upper,
control = opt.control, hessian = opt.hessian)
}
} # end of if method is GMLE stmts.
} # end of if else type is exponential.
if(is.element("shape", names(res$par))) {
if(is.element(type, c("frechet","pareto"))) {
res$par["shape"] <- abs(res$par["shape"])
if(res$par["shape"] == 0) {
warning("fevd: shape parameter estimated to be zero. Re-setting to be 1e16.")
res$par["shape"] <- 1e16
}
} else {
if(is.element(type, c("weibull","beta"))) res$par["shape"] <- -abs(res$par["shape"])
if(res$par["shape"] == 0) {
warning("fevd: shape parameter estimated to be zero. Re-setting to be -1e16.")
res$par["shape"] <- -1e16
}
}
}
res$num.pars <- list(location=ncol(designs$X.loc), scale=ncol(designs$X.sc), shape=ncol(designs$X.sh))
out$results <- res
} else if(method=="bayesian") {
if(is.element(type, c("gev","gumbel","weibull","frechet","pp"))) {
nloc <- ncol(X.loc)
loc.names <- names(initial$location)
} else {
nloc <- 0
loc.names <- NULL
}
nsc <- ncol(X.sc)
if(use.phi && is.null(initial$log.scale)) {
initial$log.scale <- log(initial$scale)
if(nsc == 1) names(initial$log.scale) <- "log.scale"
else names(initial$log.scale) <- paste("phi", 0:(nsc-1), sep="")
}
sc.names <- names(initial$log.scale)
if(!is.element(type, c("gumbel","exponential"))) {
nsh <- ncol(X.sh)
sh.names <- names(initial$shape)
} else {
nsh <- 0
sh.names <- NULL
}
np <- nloc + nsc + nsh
find.priorParams <- FALSE
if(is.null(priorFun) && is.null(priorParams)) find.priorParams <- TRUE
else if(is.null(priorFun) && (is.null(priorParams$m) || is.null(priorParams$v))) find.priorParams <- TRUE
else if(!is.null(priorFun)) {
if(priorFun == "fevdPriorDefault") {
if(is.null(priorParams)) find.priorParams <- TRUE
else if(is.null(priorParams$m) || is.null(priorParams$v)) find.priorParams <- TRUE
}
}
if(is.null(priorFun) || find.priorParams) {
if(is.null(priorFun)) priorFun <- "fevdPriorDefault"
if(find.priorParams) {
xtemp <- x
class(xtemp) <- "mle"
if(verbose) cat("\n", "Finding MLE to obtain prior means and variances.\n")
if(missing(data)) {
if(missing(span)) hold <- initializer(xtemp, u=threshold, use.phi=use.phi, type=out$type, time.units=time.units, period.basis=period.basis, blocks=blocks)
else hold <- initializer(xtemp, u=threshold, use.phi=use.phi, type=out$type, span=span, time.units=time.units, period.basis=period.basis, blocks=blocks)
} else {
if(missing(span)) hold <- initializer(xtemp, data=data, u=threshold, u.fun=threshold.fun, loc.fun=location.fun, sc.fun=scale.fun, sh.fun=shape.fun,
use.phi=use.phi, type=out$type, time.units=time.units, period.basis=period.basis, blocks=blocks)
else hold <- initializer(xtemp, data=data, u=threshold, u.fun=threshold.fun, loc.fun=location.fun, sc.fun=scale.fun, sh.fun=shape.fun,
use.phi=use.phi, type=out$type, span=span, time.units=time.units, period.basis=period.basis, blocks=blocks)
}
if(is.null(priorParams)) priorParams <- list(m=hold[1:np], v=rep(10,np))
else if(is.null(priorParams$m)) priorParams$m <- hold[1:np]
else if(is.null(priorParams$v)) priorParams$v <- rep(10, np)
} # end of if 'priorParams' is not provided by user.
} # end of if 'priorFun' is not provided by user.
out$priorFun <- priorFun
out$priorParams <- priorParams
if(is.null(proposalFun)) {
proposalFun <- "fevdProposalDefault"
if(is.null(proposalParams)) proposalParams <- list(sd=rep(0.1,np))
} # end of if 'proposalFun' not supplied by user stmts.
out$proposalFun <- proposalFun
out$proposalParams <- proposalParams
# Eric 05/30/14
chain.info <- matrix(NA, iter, np + 2)
colnames(chain.info) <- c(loc.names, sc.names, sh.names, "loglik", "prior")
chain.info[2:iter, 1:np] <- 0
res <- matrix(NA, iter, np+1)
res[,np+1] <- 0
colnames(res) <- c(loc.names, sc.names, sh.names, "new")
if(nloc > 0) res[1,1:nloc] <- initial$location
if(use.phi) res[1,(nloc+1):(nloc+nsc)] <- initial$log.scale
else res[1,(nloc+1):(nloc+nsc)] <- initial$scale
if( type != "Gumbel" ) res[1,(nloc+nsc+1):np] <- initial$shape
theta.i <- res[1,1:np]
if(verbose) {
cat("\n", "Finding log-Likelihood of initial parameter values:\n")
print(theta.i)
}
if(!missing(data)) ll.i <- -oevd(p=res[1,], o=out, des=designs, x=x, data=data, u=threshold, span=span, npy=npy, phi=use.phi, blocks=blocks) # CJP
else ll.i <- -oevd(p=res[1,], o=out, des=designs, x=x, u=threshold, span=span, npy=npy, phi=use.phi, blocks=blocks)
if(verbose) cat("\n", "Finding prior df value of initial parameter values.\n")
p.i <- do.call(priorFun, c(list(theta=theta.i), priorParams))
chain.info[1, np + 1] <- ll.i
chain.info[1, np + 2] <- p.i
if(verbose) cat("\n", "Beginning the MCMC iterations (", iter, " total iterations)\n")
for(i in 2:iter) {
if(verbose && i <= 10) cat(i, " ")
if(verbose && i %% 100 == 0) cat(i, " ")
ord <- sample(1:np, np)
theta.star <- theta.i
acc <- 0
for(j in ord) {
# par.star <- do.call(proposalFun, c(list(p=theta.i[j], ind=j), proposalParams))
# theta.star[j] <- par.star
par.star <- do.call(proposalFun, c(list(p=theta.i, ind=j), proposalParams))
theta.star[j] <- par.star[j]
if(!missing(data)) ll.star <- -oevd(p=theta.star, o=out, des=designs, x=x, data=data, u=threshold, span=span, npy=npy, phi=use.phi, blocks=blocks) # CJP
else ll.star <- -oevd(p=theta.star, o=out, des=designs, x=x, u=threshold, span=span, npy=npy, phi=use.phi, blocks=blocks)
prior.star <- do.call(priorFun, c(list(theta=theta.star), priorParams))
look <- will.accept(ll.i=ll.i, prior.i=p.i, ll.star=ll.star, prior.star=prior.star, log=TRUE)
if(look$accept) {
p.i <- prior.star
ll.i <- ll.star
theta.i <- theta.star
acc <- acc+1
chain.info[i, j] <- 1
}
} # end of inner 'j' loop.
res[i,] <- c(theta.i, acc)
chain.info[i, (np+1):(np+2)] <- c(ll.i, p.i)
} # end of for 'i' loop.
if(verbose) cat("\n", "Finished MCMC iterations.\n")
out$results <- res
out$chain.info <- chain.info
} else stop("fevd: invalid method argument.")
out$initial.results <- inout
# if( method %in% c( "MLE", "GMLE" ) && results$convergence != 0 ) warning( print( results$convergence ) )
if(verbose) print(Sys.time() - begin.tiid)
if( method == "GMLE" ) cl <- "mle"
else cl <- tolower( method )
# class(out) <- c( "fevd", paste( "fevd.", cl ) )
class( out ) <- "fevd"
return(out)
} # end of 'fevd' function.
levd <- function(x, threshold, location, scale, shape,
type=c("GEV", "GP", "PP", "Gumbel", "Weibull", "Frechet", "Exponential", "Beta", "Pareto"),
log=TRUE, negative=TRUE, span, npy=365.25, infval=Inf, weights=1, blocks=NULL) {
# CJP added blocks argument
type <- match.arg(type)
type <- tolower(type)
n <- length(x)
w <- weights
if(length(w) == 1) w <- rep(w, n)
# if(type=="pp") {
# if(!missing(span)) {
# if((length(threshold)==1) && (length(location)==1) && (length(scale)==1)) shortcut <- TRUE
# else shortcut <- FALSE
# } else shortcut <- FALSE
# }
# CJP: if we do implement 'shortcut', can't we check the const.thresh, const.location, etc. variables? I guess the issue is that we need to pass them in to levd()
shortcut <- FALSE
if(!missing(threshold)) if(length(threshold)==1) threshold <- rep(threshold, n)
if(missing(threshold) && missing(location)) stop("levd: must supply one of threshold or location.")
if(missing(location) && !missing(threshold) && is.element(type, c("gev", "weibull", "gumbel", "frechet", "pp"))) location <- threshold
if(!missing(location)) if(length(location)==1) location <- rep(location, n)
if(missing(scale)) stop("levd: must supply a scale argument.")
else if(length(scale)==1) scale <- rep(scale, n)
if(!missing(shape)) {
if(length(shape)==1) shape <- rep(shape, n)
if(any(is.na(shape))) {
if(!negative) return(-infval)
else return(infval)
}
zero.shape <- shape==0
if(all(zero.shape) && !is.element(type, c("gumbel", "exponential"))) {
if(type=="gev") type <- "gumbel"
else if(type=="gp") type <- "exponential"
}
if(any(zero.shape) && any(!zero.shape)) {
warning("levd: some shape parameters are zero, and some not. Re-setting those that are to -1e-10.")
shape[zero.shape] <- -1e-10
}
} else shape <- 0
# end of if '!missing(shape)' stmts.
if(!is.null(blocks)){ # CJP
if(any(is.na(blocks$shape))) {
if(!negative) return(-infval)
else return(infval)
}
blocks$zero.shape <- blocks$shape==0
if(any(blocks$zero.shape) && any(!blocks$zero.shape)) {
warning("levd: some shape parameters are zero, and some not. Re-setting those that are to -1e-10.")
blocks$shape[blocks$zero.shape] <- -1e-10
}
}
if(any(scale <= 0)) {
if(!negative) return(-infval)
else return(infval)
}
if(!is.null(blocks) && any(blocks$scale <= 0)) { # CJP
if(!negative) return(-infval)
else return(infval)
}
# stop("levd: invalid scale parameter. Must be positively valued.")
if(is.element(type, c("gp","exponential","beta","pareto"))) {
excess.id <- x > threshold
threshold <- threshold[excess.id]
scale <- scale[excess.id]
shape <- shape[excess.id]
y <- (x[excess.id] - threshold)/scale
m <- length(y)
}
if(type=="gumbel") {
y <- (x - location)/scale
res <- sum(colSums(w * cbind(log(scale), y, exp(-y))))
} else if(is.element(type, c("gev","weibull","frechet"))) {
y <- (x - location)/scale
y <- 1 + shape*y
if(any(y <= 0)) {
if(!negative) return(-infval)
else return(infval)
}
if(type=="frechet") shape <- abs(shape)
else if(type=="weibull") shape <- -abs(shape)
res <- sum(w * cbind(log(scale), y^(-1/shape), log(y) * (1/shape + 1))) # CJP removed colSums() here
} else if(type=="exponential") res <- sum(w[excess.id] * cbind(log(scale), y)) # CJP removed colSums() here
else if(is.element(type,c("gp","beta","pareto"))) {
y <- 1 + shape*y
if(any(y <= 0)) {
if(!negative) return(-infval)
else return(infval)
}
if(type=="pareto") shape <- abs(shape)
else if(type=="beta") shape <- -abs(shape)
res <- sum(w[excess.id] * cbind(log(scale), log(y) * (1/shape + 1))) # CJP removed colSums()
} else if(type=="pp") {
if(any(zero.shape)) {
warning("levd: some shape parameters are zero. Re-setting them to -1e-10.")
shape[zero.shape] <- -1e-10
}
if(!is.null(blocks) && any(blocks$zero.shape)) { # CJP
warning("levd: some shape parameters are zero. Re-setting them to -1e-10.")
blocks$shape[blocks$zero.shape] <- -1e-10
}
excess.id <- x > threshold
y <- (x - location)/scale
y <- 1 + shape * y
if(!shortcut) z <- 1 + ((shape*(threshold - location))/scale) # CJP made this produce a vector not a matrix
else z <- 1 + ((shape[1]*(threshold[1] - location[1]))/scale[1])
if(any(y[excess.id] <= 0) || any(z < 0)) {
if(!negative) return(-infval)
else return(infval)
}
if(!is.null(blocks)){ # CJP
blocks$z <- 1 + ((blocks$shape*(blocks$threshold - blocks$location))/blocks$scale)
if(any(blocks$z < 0)) {
if(!negative) return(-infval)
else return(infval)
}
}
th.u <- threshold[excess.id]
sc.u <- scale[excess.id]
sh.u <- shape[excess.id]
y.u <- y[excess.id]
if(is.null(blocks)) { # CJP: sum(w*vec1 + w*vec2) may be faster than cbind
if(!shortcut) res <- sum(w[excess.id] * cbind(log(sc.u), log(y.u) * (1/sh.u + 1))) + sum(w * z^(-1/shape))/npy # CJP removed colSums
else res <- sum(w[excess.id] * cbind(log(sc.u), log(y.u) * (1/sh.u + 1))) + span * z^(-1/shape[1]) # CJP removed colSums
} else {
res <- sum(w[excess.id] * cbind(log(sc.u), log(y.u) * (1/sh.u + 1))) # pp density of exceedances
res <- res + blocks$nBlocks * mean(blocks$weights * (1 - blocks$proportionMissing) * blocks$z^(-1/blocks$shape)) # this works either if stationary and each element is a scalar or if one or more vectors, though at the moment, I think $z and $shape will always be passed in as vectors
}
} # end of which type stmts.
if(!negative) res <- -res
if(!log) res <- exp(res)
if(is.na(res)) {
if(!negative) return(-infval)
else return(infval)
}
return(res)
} # end of 'levd' function.
parcov.fevd <- function(x) {
cov.theta <- NULL
theta.hat <- x$results$par
theta.names <- names(theta.hat)
if(is.element("log.scale",theta.names)) {
theta.hat[theta.names=="log.scale"] <- exp(theta.hat[theta.names=="log.scale"])
theta.names[theta.names=="log.scale"] <- "scale"
names(theta.hat) <- theta.names
phiU <- FALSE
} else phiU <- x$par.models$log.scale
y <- datagrabber(x)
if(x$data.name[2] != "") {
xdat <- y[,1]
data <- y[,-1]
} else {
xdat <- y[,1]
data <- NULL
}
designs <- setup.design(x)
if(x$method != "GMLE") {
hold <- try(suppressWarnings(optimHess(theta.hat, oevd, gr=grlevd, o=x, des=designs, x=xdat, data=data, u=x$threshold, npy=x$npy, phi=phiU, blocks=x$blocks)), silent=TRUE) # CJP
if( !is(hold, "try-error" ) && all( !is.na(hold) ) ) {
cov.theta <- try(suppressWarnings(solve(hold)), silent=TRUE)
if(any(diag(cov.theta) <= 0)) re.do <- TRUE
else re.do <- FALSE
} else re.do <- TRUE
} else re.do <- TRUE
if(re.do) {
hold <- try(optimHess(theta.hat, oevd, o=x, des=designs, x=xdat, data=data, u=x$threshold, npy=x$npy, phi=phiU, blocks=x$blocks), silent=TRUE) # CJP
if( !is(hold, "try-error" ) && all(!is.na(hold)) ) {
cov.theta <- try(solve(hold), silent=TRUE)
if( is(cov.theta, "try-error" ) || any(diag(cov.theta) <= 0) ) cov.theta <- NULL
}
} # end of if 're.do' stmts.
return(cov.theta)
} # end of 'parcov.fevd' function.
distill.fevd <- function(x, ...) {
if( x$method=="GMLE" ) newcl <- "mle"
else newcl <- tolower( x$method )
# class(x) <- c( paste("fevd.", newcl, sep=""), class( x ) )
# UseMethod("distill",x)
res <- get( paste( "distill.fevd.", tolower( x$method ), sep = "" ) )( x = x, ... )
return( res )
} # end of 'distill.fevd' function.
distill.fevd.lmoments <- function(x, ...) {
return(x$results)
} # end of 'distill.fevd.lmoments' function.
distill.fevd.bayesian <- function(x, cov=TRUE, FUN="mean", burn.in=499, ...) {
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(cov) {
cov.theta <- cov(p)
out <- c(out, c(cov.theta))
names(out) <- c(pnames, paste(pnames[rep(1:np,np)], pnames[rep(1:np,each=np)], sep="."))
}
return(out)
} # end of 'distill.fevd.bayesian' function.
distill.fevd.mle <- function(x, cov=TRUE, ...) {
theta.hat <- x$results$par
theta.names <- names(theta.hat)
np <- length(theta.hat)
if(is.element("log.scale",theta.names)) {
theta.hat[theta.names=="log.scale"] <- exp(theta.hat[theta.names=="log.scale"])
theta.names[theta.names=="log.scale"] <- "scale"
names(theta.hat) <- theta.names
}
out <- c(theta.hat, x$results$value)
if(cov) {
cov.theta <- parcov.fevd(x)
if(!is.null(cov.theta)) {
out <- c(out, cov.theta)
names(out) <- c(theta.names, "nllh", paste(theta.names[rep(1:np,np)], theta.names[rep(1:np,each=np)], sep="."))
} else names(out) <- c(theta.names, "nllh")
} else names(out) <- c(theta.names, "nllh")
return(out)
} # end of 'distill.fevd.mle' function.
summary.fevd <- function(object, ...) {
if(object$method == "GMLE") newcl <- "mle"
else newcl <- tolower( object$method )
# class(object) <- c( paste("fevd.", newcl, sep=""), class( object ) )
# UseMethod("summary", object)
get( paste( "summary.fevd.", newcl, sep = "" ) )( object = object, ... )
} # end of 'summary.fevd' function.
summary.fevd.lmoments <- function(object, ...) {
x <- object
a <- list(...)
if(!is.null(a$silent)) silent <- a$silent
else silent <- FALSE
theta.hat <- x$results
if(!silent) {
print(x$call)
if(x$data.name[2]=="") print(paste(x$type, " Fitted to ", x$data.name[1], " using L-moments estimation."))
else print(paste(x$type, " Fitted to ", all.vars(x$x.fun), " of ", x$data.name[2], "data frame, using L-moments estimation."))
print(theta.hat)
}
invisible(theta.hat)
} # end of 'summary.fevd.lmoments' function.
summary.fevd.bayesian <- function(object, FUN="mean", burn.in=499, ...) {
x <- object
a <- list(...)
if(!is.null(a$silent)) silent <- a$silent
else silent <- FALSE
np <- dim(x$results)[2] - 1
if(!is.na(burn.in) && !is.null(burn.in) && burn.in > 0) {
ll <- x$chain.info[-(1:burn.in), np + 1]
mth <- colMeans(x$results[-(1:burn.in), -(np + 1)], na.rm = TRUE)
} else {
ll <- x$chain.info[, np + 1]
mth <- colMeans(x$results[, -(np + 1)], na.rm = TRUE)
} # end of if else 'burn.in' stmts.
const <- is.fixedfevd(x)
if(!silent) {
cat("\n")
print(x$call)
cat("\n")
print(paste("Estimation Method used: ", x$method, sep=""))
cat("\n")
}
out <- list()
out$par <- ci(x, type="parameter", FUN=FUN, burn.in=burn.in)
# if(const) out$rl <- ci(x, return.period=c(2,20,100), FUN=FUN, burn.in=burn.in)
pdim <- dim(x$results)
# out$acceptance.rate <- colMeans(x$results[2:pdim[1], -pdim[2]] != x$results[1:(pdim[1] - 1), -pdim[2] ])
out$acceptance.rate <- colMeans(x$chain.info[, 2:(pdim[2] - 1), drop = FALSE ], na.rm = TRUE)
if(!silent) {
cat("\n", "Acceptance Rates:\n")
print(out$acceptance.rate)
print(out$par)
# if(const) print(out$rl)
}
cov.theta <- distill(x, FUN=FUN, burn.in=burn.in)
cov.theta <- matrix(cov.theta[(np + 1):length(cov.theta)], np, np)
colnames(cov.theta) <- rownames(cov.theta) <- colnames(x$results)[1:np]
if(!silent) {
cat("\n", "Estimated parameter covariance matrix.\n")
print(cov.theta)
}
out$cov.theta <- cov.theta
out$se.theta <- sqrt(diag(cov.theta))
designs <- setup.design(x)
y <- c(datagrabber(x, cov.data=FALSE))
if(x$data.name[2] == "") data <- NULL
else data <- datagrabber(x, response=FALSE)
if(is.element(x$type, c("GEV", "Gumbel", "Weibull", "Frechet"))) n <- x$n
else n <- sum(y > x$threshold)
nllh <- oevd(p=out$par, o=x, des=designs, x=y, data=data, u=x$threshold, span=x$span, npy=x$npy, phi=x$par.models$log.scale, blocks=x$blocks) # CJP
out$nllh <- nllh
# out$AIC <- 2 * nllh + 2 * np
if(FUN == "postmode" || FUN == "mode") out$BIC <- 2 * nllh + np * log(n)
Dev <- -2 * ll
Dbar <- mean(Dev, na.rm = TRUE)
Dmth <- -2 * oevd(p=mth, o=x, des=designs, x=y, data=data, u=x$threshold, span=x$span, npy=x$npy, phi=x$par.models$log.scale, blocks=x$blocks)
pd <- Dbar - Dmth
out$DIC <- Dmth + 2 * pd
if(!silent) {
# cat("\n", "AIC =", out$AIC, "\n")
# cat("\n", "BIC =", out$BIC, "\n")
cat("\n", "DIC = ", out$DIC, "\n")
}
invisible(out)
} # end of 'summary.fevd.bayesian' function.
summary.fevd.mle <- function(object, ...) {
x <- object
a <- list(...)
out <- list()
cov.theta <- se.theta <- NULL
if(!is.null(a$silent)) silent <- a$silent
else silent <- FALSE
if(!silent) {
cat("\n")
print(x$call)
cat("\n")
print(paste("Estimation Method used: ", x$method, sep=""))
cat("\n")
}
if(!silent) cat("\n", "Negative Log-Likelihood Value: ", x$results$value, "\n\n")
theta.hat <- x$results$par
theta.names <- names(theta.hat)
if(is.element("log.scale",theta.names)) {
theta.hat[theta.names=="log.scale"] <- exp(theta.hat[theta.names=="log.scale"])
theta.names[theta.names=="log.scale"] <- "scale"
names(theta.hat) <- theta.names
phiU <- FALSE
} else phiU <- x$par.models$log.scale
out$par <- theta.hat
np <- length(theta.hat)
designs <- setup.design(x)
if(!is.null(x$data.pointer)) xdat <- get(x$data.pointer)
else xdat <- x$x
cov.theta <- parcov.fevd(x)
out$cov.theta <- cov.theta
if(!is.null(cov.theta)) {
se.theta <- sqrt(diag(cov.theta))
names(se.theta) <- theta.names
out$se.theta <- se.theta
}
if(!silent) {
cat("\n", "Estimated parameters:\n")
print(theta.hat)
}
if(!is.null(se.theta)) {
if(!silent) {
cat("\n", "Standard Error Estimates:\n")
print(se.theta)
}
theta.hat <- rbind(theta.hat, se.theta)
if(is.matrix(theta.hat) && dim(theta.hat)[1]==2) rownames(theta.hat) <- c("Estimate", "Std. Error")
if(is.matrix(theta.hat)) colnames(theta.hat) <- theta.names
else names(theta.hat) <- theta.names
if(!silent && !is.null(cov.theta)) {
cat("\n", "Estimated parameter covariance matrix.\n")
print(cov.theta)
}
}
nllh <- x$results$value
out$nllh <- nllh
if(is.element(x$type, c("GEV","Gumbel","Weibull","Frechet"))) n <- x$n
else {
y <- c(datagrabber(x, cov.data=FALSE))
n <- sum(y > x$threshold)
}
out$AIC <- 2 * nllh + 2 * np
out$BIC <- 2 * nllh + np * log(n)
if(!silent) {
cat("\n", "AIC =", out$AIC, "\n")
cat("\n", "BIC =", out$BIC, "\n")
}
invisible(out)
} # end of 'summary.fevd.mle' function.
print.fevd <- function(x, ...) {
summary(x)
invisible()
} # end of 'print.fevd' function.
grlevd <- function(p, o, des, x, data, u=NULL, span, npy, phi=FALSE, blocks=NULL) {
##
## Function to estimate the likelihood
## parameter gradients.
##
npar <- length(p)
type <- o$type
n <- length(x)
w <- o$weights
if(is.null(w)) w <- 1
if(length(w) == 1) w <- rep(w, n)
shortcut <- FALSE
dpar <- matrix(NA, n, npar)
if(!is.null(blocks)) # CJP
blocks$dpar <- matrix(NA, blocks$nBlocks, npar)
if(is.element(type, c("PP","GEV","Gumbel","Weibull","Frechet"))) {
X.loc <- des$X.loc
nloc <- ncol(X.loc)
loc <- rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc)
dpar[,1:nloc] <- X.loc
if(!is.null(blocks)) { # CJP
blocks$loc <- rowSums(matrix(p[1:nloc], blocks$nBlocks, nloc, byrow=TRUE)*blocks$designs$X.loc)
blocks$dpar[,1:nloc] <- blocks$designs$X.loc
}
} else {
nloc <- 0
loc <- NULL
if(!is.null(blocks)) # CJP
blocks$loc <- NULL
}
X.sc <- des$X.sc
nsc <- ncol(X.sc)
if(phi) scale <- exp(rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc))
else scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc)
dpar[,(nloc+1):(nloc+nsc)] <- X.sc
if(phi) phi.gradTerm <- scale else phi.gradTerm <- 1 # CJP2
if(!is.null(blocks)) { # CJP
if(phi) {
blocks$scale <- exp(rowSums(matrix(p[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE)*blocks$designs$X.sc))
blocks$phi.gradTerm <- blocks$scale
} else {
blocks$scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE)*blocks$designs$X.sc)
blocks$phi.gradTerm <- 1
}
blocks$dpar[,(nloc+1):(nloc+nsc)] <- blocks$designs$X.sc
}
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)
dpar[,(nloc+nsc+1):(nloc+nsc+nsh)] <- X.sh
if(!is.null(blocks)) { # CJP
blocks$shape <- rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], blocks$nBlocks, nsh, byrow=TRUE)*blocks$designs$X.sh)
blocks$dpar[,(nloc+nsc+1):(nloc+nsc+nsh)] <- blocks$designs$X.sh
}
} else {
nsh <- 0
shape <- NULL
if(!is.null(blocks)) # CJP
blocks$shape <- NULL
}
if(is.element(type, c("Frechet","Pareto"))) shape <- abs(shape)
else if(is.element(type, c("Weibull","Beta"))) shape <- -abs(shape)
res <- p+NA
if(is.element(type, c("GEV","Weibull","Frechet"))) {
sh.sc <- shape/scale
sh.1 <- 1/shape + 1
x.fmu <- x - loc
z <- 1 + sh.sc*x.fmu
for(i in 1:npar) { # CJP: don't need colSums and sum(w*vec1 + w*vec2 + w*vec3) may be faster than cbind
if(i <= nloc) res[i] <- sum(colSums(w * cbind(-sh.1 * (sh.sc) * dpar[,i]/z, (dpar[,i]/scale)*z^(-sh.1))))
else if(i <= nloc + nsc) res[i] <- sum(colSums(w * cbind(phi.gradTerm*dpar[,i]/scale, -phi.gradTerm * sh.1 * (shape/(scale^2))*x.fmu*dpar[,i]/z, (z^(-sh.1)/(scale^2)) * dpar[,i] * x.fmu * phi.gradTerm))) # CJP2
else res[i] <- sum(colSums(w * cbind(sh.1 * (dpar[,i]/scale) * x.fmu/z, -dpar[,i]/(shape^2)*log(z),
z^(-1/shape) * ((dpar[,i]/(shape^2))*log(z) - (dpar[,i]/scale)*x.fmu/(z*shape)))))
} # end of for 'i' loop.
} else if(type == "Gumbel") {
x.fmu <- x - loc
z <- x.fmu/scale
for(i in 1:npar) {
if(i <= nloc) res[i] <- sum(colSums(w * cbind(-dpar[,i]/scale, (dpar[,i]/scale)*exp(-z))))
else res[i] <- sum(colSums(w * cbind(phi.gradTerm*dpar[,i]/scale, -((x.fmu/(scale^2))*dpar[,i])*phi.gradTerm, phi.gradTerm*(dpar[,i]/(scale^2))*x.fmu*exp(-z)))) # CJP2 - Eric, please check this
} # end of for 'i' loop.
} else if(is.element(type, c("GP","Beta","Pareto"))) {
excess.id <- x > u
x <- x[excess.id]
scale <- scale[excess.id]
shape <- shape[excess.id]
u <- u[excess.id]
dparu <- dpar[excess.id,]
x.fu <- x - u
z <- 1 + (shape/scale)*x.fu
sh.1 <- 1/shape + 1
if(phi) phi.gradTermExc <- scale else phi.gradTermExc <- 1 # CJP2
for(i in 1:npar) {
if(i <= nsc) res[i] <- sum(colSums(w[excess.id] * cbind(phi.gradTermExc * dparu[,i]/scale, - ((phi.gradTermExc*sh.1*dparu[,i]*shape*x.fu/(scale^2))/z)))) # CJP2 - Eric, please check this
else res[i] <- sum(colSums(w[excess.id] * cbind((-dparu[,i]/(shape^2))*log(z), sh.1*((dparu[,i]/scale)*x.fu)/z)))
} # end of for 'i' loop.
} else if(type=="Exponential") {
excess.id <- x > u
if(!is.matrix(dpar)) dpar <- cbind(dpar)
dparu <- dpar[excess.id,]
if(!is.matrix(dparu)) dparu <- cbind(dparu)
scale <- scale[excess.id]
xu <- x[excess.id] - u[excess.id]
if(phi) phi.gradTermExc <- scale else phi.gradTermExc <- 1 # CJP2
for(i in 1:npar) res[i] <- sum(colSums(w[excess.id] * cbind(phi.gradTermExc*dparu[,i]/scale, - phi.gradTermExc * xu * dparu[,i]/(scale^2)))) # CJP2 - Eric, please check this
} else if(type=="PP") {
excess.id <- x > u
mu <- loc[excess.id]
sc <- scale[excess.id]
sh <- shape[excess.id]
shu.scu <- sh/sc
uu <- u[excess.id]
xu <- x[excess.id]
dparu <- dpar[excess.id,]
x.fmu <- xu - mu
if(!shortcut) u.l <- u - loc
else u.l <- u[1] - loc[1]
zu <- 1 + shu.scu*x.fmu
if(!shortcut) Z <- matrix(1 + (shape/scale)*u.l, ncol=1)
else Z <- 1 + (shape[1]/scale[1])*u.l
shu.1 <- 1/sh + 1
if(!shortcut) sh.1 <- 1/shape + 1
else sh.1 <- 1/shape[1] + 1
zpsh1 <- Z^(-sh.1)
if(phi) phi.gradTermExc <- sc else phi.gradTermExc <- 1 # CJP2
if(!is.null(blocks)) {
blocks$u.l <- blocks$threshold - blocks$loc
blocks$Z <- 1 + (blocks$shape/blocks$scale)*blocks$u.l
blocks$sh.1 <- 1/blocks$shape + 1
blocks$zpsh1 <- blocks$Z^(-blocks$sh.1)
}
for(i in 1:npar) { # CJP: I changed sum(colSums(...)) here to sum()
if(!shortcut) { # CJP: same comment about sum(w*vec1 + w*vec2)
if(is.null(blocks)) { # CJP
if(i <= nloc) res[i] <- sum(w[excess.id] * cbind(-shu.1*dparu[,i]*shu.scu/zu)) + colSums(w * zpsh1 * dpar[,i]/scale)/npy # CJP: I think this can be sum not colSums
else if(i <= nloc + nsc) res[i] <- sum(w[excess.id] * cbind(phi.gradTermExc * dparu[,i]/sc, -(shu.1*dparu[,i]*sh*x.fmu/(zu*sc^2))*phi.gradTermExc)) + colSums(w * phi.gradTerm * dpar[,i]*u.l*zpsh1/(scale^2))/npy # CJP: I think this can be sum not colSums; CJP2 - added phi.gradTerm
else res[i] <- sum(w[excess.id] * cbind(-(dparu[,i]/(sh^2))*log(zu), (shu.1*(dparu[,i]/sc)*x.fmu)/zu)) +
colSums(w * (exp(-log(Z)/shape)*((dpar[,i]/(shape^2))*log(Z) - (dpar[,i]/(shape*scale)*u.l)/Z)))/npy
} else { # w/ blocks
if(i <= nloc) res[i] <- sum(w[excess.id] * cbind(-shu.1*dparu[,i]*shu.scu/zu)) +
blocks$nBlocks * mean(blocks$weights * (1 - blocks$proportionMissing) * blocks$zpsh1 * blocks$dpar[,i]/blocks$scale)
else if(i <= nloc + nsc) res[i] <- sum(w[excess.id] * cbind(phi.gradTermExc*dparu[,i]/sc, -(shu.1*dparu[,i]*sh*x.fmu/(zu*sc^2))*phi.gradTermExc)) +
blocks$nBlocks * mean(blocks$weights * (1 - blocks$proportionMissing) * blocks$phi.gradTerm * blocks$dpar[,i]*blocks$u.l*blocks$zpsh1/(blocks$scale^2)) # CJP2: added blocks$phi.gradTerm
else res[i] <- sum(w[excess.id] * cbind(-(dparu[,i]/(sh^2))*log(zu), (shu.1*(dparu[,i]/sc)*x.fmu)/zu)) +
blocks$nBlocks * mean(blocks$weights * (1 - blocks$proportionMissing) * (exp(-log(blocks$Z)/blocks$shape)*((blocks$dpar[,i]/(blocks$shape^2))*log(blocks$Z) - (blocks$dpar[,i]/(blocks$shape*blocks$scale)*blocks$u.l)/blocks$Z)))
}
} else {
if(i <= nloc) res[i] <- -sum(w[excess.id] * cbind(shu.1*dparu[,i]*shu.scu/zu)) + span * (zpsh1*dpar[1,i]/scale[1])
else if(i <= nloc + nsc) res[i] <- sum(w[excess.id] * cbind(phi.gradTermExc*dparu[,i]/sc, -(shu.1*dparu[,i]*sh*x.fmu/(zu*sc^2))*phi.gradTermExc)) + span * (phi.gradTerm*dpar[1,i]*u.l*zpsh1/(scale[1]^2)) # CJP2: added phi.gradTerm
else res[i] <- sum(w[excess.id] * cbind(-(dparu[,i]/(sh^2))*log(zu), shu.1*(dparu[,i]/sc)*x.fmu/zu)) +
span * (exp(-log(Z)/shape[1])*((dpar[1,i]/(shape[1]^2))*log(Z) - (dpar[1,i]/(shape[1]*scale[1])*u.l)/Z))
}
} # end of for 'i' loop.
} # end of which gradient likelihood to use stmts.
return(res)
} # end of 'grlevd' function.
grlevdTracer <- function(x, p=c(0,1,0), which.vary=1, p1.range, p2.range=NULL, threshold=NULL,
threshold.fun=~1, location.fun=~1, scale.fun=~1, shape.fun=~1, data, phi=FALSE,
type=c("GEV","GP","PP","Gumbel","Weibull","Frechet","Exponential","Beta","Pareto"),
nint=100, npy=365.25, weights = 1, blocks=NULL, na.action=na.fail, par1.name="", par2.name="", zl.lh=NULL,
plot=TRUE, ...) {
# CJP note: this function now assumes that blocks contains the blocks data and design matrices already constructed within fevd and passed out as the $blocks element of the output from fevd - I didn't want to repeat all of that code here
# CJP2 : changed 'u' to 'threshold' and 'threshold' to 'threshold.fun'
type <- match.arg(type)
id <- which.vary
data.name <- c(deparse(substitute(x)), deparse(substitute(data)))
if(!is.function(na.action)) na.action <- get(na.action)
if(!missing(data)) {
if(length(x) != nrow(data)) stop("fevd: data must have same number of rows as the length of x.")
tmp <- cbind(x, data)
tmp <- na.action(tmp)
x <- tmp[,1]
data <- tmp[,-1]
n <- length(x)
} else {
x <- na.action(x)
n <- length(x)
}
w <- weights
if(length(w) == 1) w <- rep(w, n)
if(!is.null(threshold)) {
const.thresh <- check.constant(threshold.fun)
if(!missing(data)) X.u <- setup.design(x=threshold.fun, data=data, n=n, const=const.thresh, dname="threshold")
else X.u <- setup.design(x=threshold.fun, n=n, const=const.thresh, dname="threshold")
threshold <- rowSums(threshold*X.u)
} # end of if 'threshold (u)' stmts.
p1 <- seq(p1.range[1],p1.range[2],,nint)
if(!is.null(p2.range)) p2 <- seq(p2.range[1],p2.range[2],,nint)
loc.exists <- is.element(type, c("GEV","PP","Gumbel","Weibull","Frechet"))
sh.exists <- !is.element(type, c("Gumbel","Exponential"))
designs <- list()
if(loc.exists) {
const.loc <- check.constant(location.fun)
if(!missing(data)) X.loc <- setup.design(x=location.fun, data=data, n=n, const=const.loc, dname="location.fun")
else X.loc <- setup.design(x=location.fun, n=n, const=const.loc, dname="location")
designs$X.loc <- X.loc
nloc <- ncol(X.loc)
} else {
nloc <- 0
const.loc <- NULL
}
const.scale <- check.constant(scale.fun)
if(!missing(data)) X.sc <- setup.design(x=scale.fun, data=data, n=n, const=const.scale, dname="scale.fun")
else X.sc <- setup.design(x=scale.fun, n=n, const=const.scale, dname="scale.fun")
designs$X.sc <- X.sc
nsc <- ncol(X.sc)
if(sh.exists) {
const.shape <- check.constant(shape.fun)
if(!missing(data)) X.sh <- setup.design(x=shape.fun, data=data, n=n, const=const.shape, dname="shape.fun")
else X.sh <- setup.design(x=shape.fun, n=n, const=const.shape, dname="shape")
designs$X.sh <- X.sh
nsh <- ncol(X.sh)
if(is.element(type, c("Frechet","Pareto"))) shape <- abs(shape)
else if(is.element(type, c("Weibull","Beta"))) shape <- -abs(shape)
} else {
nsh <- 0
const.shape <- NULL
shape <- 0
} # end of 'shape.range' stmts.
obj <- list(type=type, const.loc=const.loc, const.scale=const.scale, const.shape=const.shape,
par.models=list(location=location.fun, scale=scale.fun, shape=shape.fun, phi=phi, weights=w,
term.names=list(location=all.vars(location.fun), scale=all.vars(scale.fun), shape=all.vars(shape.fun))))
num.vary <- length(id)
if(num.vary==1) {
lh <- lhgr <- numeric(nint)+NA
for(i in 1:nint) {
p[id] <- p1[i]
p.tmp <- p
scale <- rowSums(matrix(p.tmp[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc)
if(phi) scale <- exp(scale)
if(!is.null(blocks)) {
blocks$location <- rowSums(matrix(p[1:nloc], blocks$nBlocks, nloc, byrow=TRUE)*blocks$designs$X.loc)
blocks$scale <- rowSums(matrix(p.tmp[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE)*blocks$designs$X.sc)
if(phi) blocks$scale <- exp(blocks$scale)
blocks$shape=rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], blocks$nBlocks, nsh, byrow=TRUE)*blocks$designs$X.sh)
}
if(loc.exists && sh.exists) lh[i] <- levd(x=x, threshold=threshold, location=rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc),
scale=scale, shape=rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], n, nsh, byrow=TRUE)*X.sh),
type=type, npy=npy, weights=w, blocks=blocks)
else if(loc.exists) lh[i] <- levd(x=x, threshold=threshold, location=rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc),
scale=scale, type=type, npy=npy, weights=w, blocks=blocks)
else if(sh.exists) lh[i] <- levd(x=x, threshold=threshold, scale=scale, shape=rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], n, nsh, byrow=TRUE)*X.sh),
type=type, npy=npy, weights=w)
else lh[i] <- levd(x=x, threshold=threshold, scale=scale, type=type, npy=npy, weights=w, blocks=blocks)
if(!missing(data)) lhgr[i] <- grlevd(p=p.tmp, o=obj, des=designs, x=x, data=data, u=threshold, npy=npy, phi=phi, blocks=blocks)[id]
else lhgr[i] <- grlevd(p=p.tmp, o=obj, des=designs, x=x, u=threshold, npy=npy, phi=phi, blocks=blocks)[id]
} # end of for 'i' loop.
res <- cbind(p1, lh, lhgr)
colnames(res) <- c("par.values","nllh","nllh.gr")
attributes(res)$data.name <- data.name
attributes(res)$par.name <- par1.name
attributes(res)$model.name <- type
} else if(num.vary==2) {
lh <- lhgr1 <- lhgr2 <- matrix(NA, nint, nint)
for(i in 1:nint) {
p[id[1]] <- p1[i]
p.tmp <- p
if(const.scale && phi) p.tmp[nloc+1] <- exp(p.tmp[nloc+1])
for(j in 1:nint) {
p[id[2]] <- p2[j]
p.tmp <- p
scale <- rowSums(matrix(p.tmp[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc)
if(phi) scale <- exp(scale)
if(!is.null(blocks)) {
blocks$location <- rowSums(matrix(p[1:nloc], blocks$nBlocks, nloc, byrow=TRUE)*blocks$designs$X.loc)
blocks$scale <- rowSums(matrix(p.tmp[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE)*blocks$designs$X.sc) # fixed by CJP2 6/11/13 (was 'blocks$X.sc')
if(phi) blocks$scale <- exp(blocks$scale)
blocks$shape=rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], blocks$nBlocks, nsh, byrow=TRUE)*blocks$designs$X.sh)
}
if(loc.exists && sh.exists) lh[i,j] <- levd(x=x, threshold=threshold, location=rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc),
scale=scale, shape=rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], n, nsh, byrow=TRUE)*X.sh),
type=type, npy=npy, weights=w, blocks=blocks)
else if(loc.exists) lh[i,j] <- levd(x=x, threshold=threshold, location=rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc),
scale=scale, type=type, npy=npy, weights=w, blocks=blocks)
else if(sh.exists) lh[i,j] <- levd(x=x, threshold=threshold, scale=scale, shape=rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], n, nsh, byrow=TRUE)*X.sh),
type=type, npy=npy, weights=w)
else lh[i,j] <- levd(x=x, threshold=threshold, scale=scale, type=type, npy=npy, weights=w, blocks=blocks)
if(!missing(data)) hold <- grlevd(p=p.tmp, o=obj, des=designs, x=x, data=data, u=threshold, npy=npy, blocks=blocks)[id]
else hold <- grlevd(p=p.tmp, o=obj, des=designs, x=x, u=threshold, npy=npy, blocks=blocks)[id]
lhgr1[i,j] <- hold[1]
lhgr2[i,j] <- hold[2]
} # end of for 'j' loop.
} # end of for 'i' loop.
par.values <- list(p1,p2)
names(par.values) <- c(paste("gr1.", par1.name, sep=""), paste("gr2.", par2.name, sep=""))
res <- list(nllh=lh, nllh.gr1=lhgr1, nllh.gr2=lhgr2, p1=p1, p2=p2, call=match.call(), data.name=data.name, par.name=c(par1.name, par2.name), model=type)
parvar <- cbind(rep(p1, nint), rep(p2, each=nint))
} # end of 1 or two varied parameters stmts.
class(res) <- "grlevdTracer"
if(plot) plot(res, type="both", which.gr=1:2, ...)
invisible(res)
} # end of 'grlevdTracer' function.
plot.grlevdTracer <- function(x, type=c("likelihood", "gradient", "both"), which.gr=1, log=TRUE, negative=TRUE, ...) {
type <- match.arg(type)
a <- list(...)
if(is.element(type, c("likelihood","both")) && is.null(a$ylab)) {
if(log && negative) nllh.ylab <- "Negative Log-Likelihood"
else if(log) nllh.ylab <- "Log-Likelihood"
else if(!negative) nllh.ylab <- "Likelihood"
else nllh.ylab <- "Negative Likelihood"
}
if(is.element(type, c("gradient", "both"))) gr.ylab <- "Negative Log-Likelihood Gradient"
if(is.matrix(x)) {
tmp <- attributes(x)
dname <- tmp$data.name[1]
pname <- tmp$par.name
mname <- tmp$model.name
if(is.element(type, c("likelihood","both"))) {
X <- x[,"nllh"]
if(!any(is.finite(X))) {
skip <- TRUE
} else {
skip <- FALSE
if(!log && !negative) X <- exp(-X)
else if(!negative) X <- -X
else if(!log) X <- -exp(-X)
}
}
if(is.null(a$ylim)) {
if(is.element(type, c("likelihood","both"))) yl <- range(X, finite=TRUE)
if(type=="both") par(mfrow=c(2,1))
if(is.element(type, c("likelihood","both"))) {
if(skip) plot(1:length(X), rep(0, length(X)), type="n", yaxt="n", main="No finite negative log-likelihood values for this range")
else {
if(is.null(a$ylab) && is.null(a$main)) {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", ylab=nllh.ylab, main=paste(mname, dname, sep="\n"), ylim=yl, ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, ylab=nllh.ylab, main=paste(mname, dname, sep="\n"), ylim=yl, ...)
else plot(x[,"par.values"], X, type="l", ylab=nllh.ylab, main=paste(mname, dname, sep="\n"), ylim=yl, ...)
} else if(is.null(a$ylab)) {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", ylab=nllh.ylab, ylim=yl, ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, ylab=nllh.ylab, ylim=yl, ...)
else plot(x[,"par.values"], X, type="l", ylab=nllh.ylab, ylim=yl, ...)
} else if(is.null(a$main)) {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", main=paste(mname, dname, sep="\n"), ylim=yl, ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, main=paste(mname, dname, sep="\n"), ylim=yl, ...)
else plot(x[,"par.values"], X, type="l", ylim=yl, main=paste(mname, dname, sep="\n"), ...)
} else {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", ylim=yl, ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, ylim=yl, ...)
else plot(x[,"par.values"], X, type="l", ylim=yl, ...)
}
} # end of if else skip stmts.
} # end of if type is likelihood or both stmts.
if(is.element(type, c("gradient", "both"))) {
if(!any(is.finite(x[,"nllh.gr"]))) {
plot(1:length(x[,"nllh.gr"]), rep(0, length(x[,"nllh.gr"])), type="n", yaxt="n", main="No finite negative log-likelihood gradient values for this range")
} else {
if(is.null(a$main) && type=="both") m1 <- ""
else if(is.null(a$main)) m1 <- paste(mname, dname, sep="\n")
if(is.null(a$ylab) && is.null(a$main)) {
# gr.ylab <- paste(gr.ylab, " wrt ", pname, sep="")
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, ylab=gr.ylab, main=m1, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", ylab=gr.ylab, main=m1, ...)
} else if(is.null(a$ylab)) {
# gr.ylab <- paste(gr.ylab, " wrt ", pname, sep="")
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, ylab=gr.ylab, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", ylab=gr.ylab, ...)
} else if(is.null(a$main)) {
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, main=m1, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", main=m1, ...)
} else {
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", ...)
} # end of if ylab and main supplied or not stmts.
abline(h=0, lty=2)
} # end of if else no finite gradient values stmts.
} # end of if type is gradient or both stmts.
} else {
if(type=="both") par(mfrow=c(2,1))
if(is.element(type, c("likelihood","both"))) {
if(is.null(a$ylab) && is.null(a$main)) {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", ylab=nllh.ylab, main=paste(mname, dname, sep="\n"), ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, ylab=nllh.ylab, main=paste(mname, dname, sep="\n"), ...)
else plot(x[,"par.values"], X, type="l", ylab=nllh.ylab, main=paste(mname, dname, sep="\n"), ...)
} else if(is.null(a$ylab)) {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", ylab=nllh.ylab, ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, ylab=nllh.ylab, ...)
else plot(x[,"par.values"], X, type="l", ylab=nllh.ylab, ...)
} else {
if((type == "both") && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab="", ...)
else if(type == "likelihood" && is.null(a$xlab)) plot(x[,"par.values"], X, type="l", xlab=pname, ...)
else plot(x[,"par.values"], X, type="l", ...)
}
} # end of if type is likelihood or both stmts.
if(is.element(type, c("gradient", "both"))) {
if(is.null(a$main) && type=="both") m1 <- ""
else if(is.null(a$main)) m1 <- paste(mname, dname, sep="\n")
if(is.null(a$ylab) && is.null(a$main)) {
gr.ylab <- paste(gr.ylab, " wrt ", pname, sep="")
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, ylab=gr.ylab, main=m1, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", ylab=gr.ylab, main=m1, ...)
} else if(is.null(a$ylab)) {
gr.ylab <- paste(gr.ylab, " wrt ", pname, sep="")
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, ylab=gr.ylab, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", ylab=gr.ylab, ...)
} else if(is.null(a$main)) {
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, main=m1, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", main=m1, ...)
} else {
if(is.null(a$xlab)) plot(x[,"par.values"], x[,"nllh.gr"], type="l", xlab=pname, ...)
else plot(x[,"par.values"], x[,"nllh.gr"], type="l", ...)
} # end of if ylab and main supplied or not stmts.
abline(h=0, lty=2)
} # end of if type is gradient or both stmts.
} # end of if ylim is passed or not stmts.
} else {
dname <- x$data.name[1]
p1name <- x$par.name[1]
p2name <- x$par.name[2]
p1 <- x$p1
p2 <- x$p2
nint <- length(p1)
if(is.element(type, c("likelihood","both"))) {
X <- x$nllh
if(!log && !negative) X <- exp(-X)
else if(!negative) X <- -X
else if(!log) X <- -exp(-X)
minloc <- cbind(rep(p1, nint), rep(p2, each=nint))[which.min(c(X)),]
}
if(!(which.gr[1]==1 || which.gr[1]==2 || all(which.gr==1:2))) stop("plot.grlevdTracer: invalid which.gr argument. Must be 1, 2 or 1:2.")
if(is.null(a$xlab)) xl <- p1name
else xl <- a$xlab
if(is.null(a$ylab)) yl <- p2name
else yl <- a$ylab
if(!is.null(a$zlim)) zl.lh <- zl.gr1 <- zl.gr2 <- a$zlim
else {
zl.lh <- range(c(X), finite=TRUE)
zl.gr1 <- range(c(x$nllh.gr1), finite=TRUE)
zl.gr2 <- range(c(x$nllh.gr2), finite=TRUE)
}
# if(is.null(a$col)) zcol <- tim.colors(64)
# else zcol <- a$col
if(is.null(a$horizontal)) hz <- TRUE
else hz <- a$horizontal
if(is.null(a$legend.mar)) legmar <- ifelse(hz, 3.1, 5.1)
else legmar <- a$legend.mar
if(type=="both") {
if(length(which.gr)==1) par(mfrow=c(1,2))
else par(mfrow=c(1,3))
}
if(is.null(a$main)) {
if(is.element(type, c("likelihood","both"))) m1 <- paste(nllh.ylab, "\n", x$model, ": ", x$data.name[1], sep="")
else m1 <- ""
if(length(which.gr)==1) m2 <- m3 <- paste(x$model, ": Negative Log-Likelihood\nGradient wrt ", x$par.name[which.gr], " (", x$data.name[1], ")", sep="")
else {
m2 <- paste(x$model, ": Negative Log-Likelihood\nGradient wrt ", p1name, " (", x$data.name[1], ")", sep="")
m3 <- paste(x$model, ": Negative Log-Likelihood\nGradient wrt ", p2name, " (", x$data.name[1], ")", sep="")
}
} else m1 <- m2 <- m3 <- a$main
if(is.element(type,c("likelihood","both"))) {
image(p1, p2, X, main=m1, xlab=xl, ylab=yl, zlim=zl.lh)
points(minloc[1], minloc[2], pch="x", col="darkorange")
# image.plot(p1, p2, X, zlim=zl.lh, horizontal=hz, legend.only=TRUE, legend.mar=legmar)
} # end of if likelihood or both stmts.
if(is.element(type, c("gradient","both"))) {
if(is.element(1, which.gr)) {
image(p1, p2, x$nllh.gr1, main=m2, xlab=xl, ylab=yl, zlim=zl.gr1)
# image.plot(p1, p2, x$nllh.gr1, zlim=zl.gr1, horizontal=hz, legend.only=TRUE, legend.mar=legmar)
}
if(is.element(2, which.gr)) {
image(p1, p2, x$nllh.gr2, main=m3, xlab=xl, ylab=yl, zlim=zl.gr2)
# image.plot(p1, p2, x$nllh.gr2, zlim=zl.gr2, horizontal=hz, legend.only=TRUE, legend.mar=legmar)
}
} # end of if type is gradient or both stmts.
}# end of if one or two parameters varied stmts.
invisible()
} # end of 'plot.grlevdTracer' function.
check.constant <- function(x) {
# internal function to determine whether a parameter is to be modeled
# as constant (TRUE) or not (FALSE).
if(is.formula(x)) {
if(length(as.character(x))==2 && as.character(x)[2] == "1") return(TRUE)
else return(FALSE)
} else {
if(length(x)==1) return(TRUE)
else if(all(x == x[1])) return(TRUE)
else return(FALSE)
}
} # end of 'check.constant' function.
setup.design <- function(x, ...) {
UseMethod("setup.design", x)
} # end of 'setup.design' function.
setup.design.default <- function(x, ..., data, n, const, dname) {
if(missing(data)) data <- NULL
if(missing(const)) const <- check.constant(x)
if(const || is.vector(x)) return(matrix(1, n, 1))
else if(is.formula(x)) {
if(is.null(data)) return(model.matrix(x))
else return(model.matrix(x, data))
} else stop(paste("setup.design: invalid ", dname, " parameter.", sep=""))
} # end of 'setup.design' function.
setup.design.fevd <- function(x, ...) {
out <- list()
n <- x$n
mods <- x$par.models
dname <- x$data.name[2]
data <- datagrabber(x, response=FALSE)
out$X.u <- setup.design(mods$threshold, data=data, n=n, dname="threshold")
out$X.loc <- setup.design(mods$location, data=data, n=n, dname="location")
out$X.sc <- setup.design(mods$scale, data=data, n=n, dname="scale")
out$X.sh <- setup.design(mods$shape, data=data, n=n, dname="shape")
return(out)
} # end of 'setup.design.fevd' function.
oevd <- function(p, o, des, x, data=NULL, u=NULL, span, npy, phi=FALSE, blocks=NULL) {
##
## Function to be optimized by 'optim'.
## Takes the parameter design matrices and
## calculates the resulting vector of parameters,
## then calls 'levd' to get the negative log-likelihood.
type <- o$type
n <- length(x)
w <- o$weights
if(length(w) == 1) w <- rep(w, n)
if(!is.element(type, c("GP","Beta","Pareto","Exponential"))) {
X.loc <- des$X.loc
nloc <- ncol(X.loc)
loc <- rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc)
} else {
nloc <- 0
loc <- 0
}
X.sc <- des$X.sc
nsc <- ncol(X.sc)
if(phi) scale <- exp(rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc))
else scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc)
if(!is.element(o$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 {
nsh <- 0
shape <- 0
}
if(!is.null(blocks)) { # could check for const.loc, etc. and just do sum(p[1:nloc]*blocks$designs$X.param[1,])
if(is.element(type, c("PP"))) {
blocks$location <- rowSums(matrix(p[1:nloc], blocks$nBlocks, nloc, byrow=TRUE)*blocks$designs$X.loc)
}
if(phi) blocks$scale <- exp(rowSums(matrix(p[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE)*blocks$designs$X.sc))
else blocks$scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE)*blocks$designs$X.sc)
blocks$shape <- rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], blocks$nBlocks, nsh, byrow=TRUE)*blocks$designs$X.sh)
}
if(missing(span)) res <- levd(x=x, threshold=u, location=loc, scale=scale, shape=shape, type=type, npy=npy, infval=1e10, weights=w, blocks=blocks)
else res <- levd(x=x, threshold=u, location=loc, scale=scale, shape=shape, type=type, span=span, npy=npy, infval=1e10, weights=w, blocks=blocks)
return(res)
} # end of 'oevd' function.
oevdgen <- function(p, o, des, x, data = NULL, u = NULL, span, npy, phi = FALSE, priorFun, priorParams = NULL, blocks = NULL) {
##
## Function to be optimized by 'optim'.
## Takes the parameter design matrices and
## calculates the resulting vector of parameters,
## then calls 'levd' to get the negative log-likelihood.
type <- o$type
n <- length(x)
w <- o$weights
if(length(w) == 1) w <- rep(w, n)
if(!is.element(type, c("GP","Beta","Pareto","Exponential"))) {
X.loc <- des$X.loc
nloc <- ncol(X.loc)
loc <- rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE) * X.loc)
if(nloc == 1) lname <- "location"
else lname <- paste("mu", 0:(nloc - 1), sep = "")
} else {
nloc <- 0
loc <- 0
lname <- NULL
}
X.sc <- des$X.sc
nsc <- ncol(X.sc)
if(phi) scale <- exp(rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE) * X.sc))
else scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE) * X.sc)
if(nsc == 1) scname <- "scale"
else scname <- paste("phi", 0:(nsc - 1), sep="")
if(!is.null(des$X.sh)) {
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)
if(nsh == 1) shname <- "shape"
else shname <- paste("xi", 0:(nsh - 1), sep = "")
} else {
nsh <- 0
shape <- 0
shname <- NULL
}
if(!is.null(blocks)) {
blocks$location <- rowSums(matrix(p[1:nloc], blocks$nBlocks, nloc, byrow=TRUE) * blocks$designs$X.loc)
if(phi) blocks$scale <- exp(rowSums(matrix(p[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE) * blocks$designs$X.sc))
else blocks$scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], blocks$nBlocks, nsc, byrow=TRUE) * blocks$designs$X.sc)
blocks$shape <- rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], blocks$nBlocks, nsh, byrow=TRUE) * blocks$designs$X.sh)
}
if(missing(span)) res <- levd(x=x, threshold=u, location=loc, scale=scale, shape=shape, type=type, npy=npy, infval=1e10, weights=w, blocks=blocks)
else res <- levd(x=x, threshold=u, location=loc, scale=scale, shape=shape, type=type, span=span, npy=npy, infval=1e10, weights=w, blocks=blocks)
# newp <- cbind(loc, scale, shape)
# colnames(newp) <- c(lname, scname, shname)
# res2 <- do.call(priorFun, c(list(x = newp), priorParams))
res2 <- do.call(priorFun, c(list(x = p), priorParams))
return(res + res2)
} # end of 'oevdgen' function.
shapePriorBeta <- function( x, p, q ) {
xi <- x[ length(x) ]
if( xi >= -0.5 & xi <= 0.5 ) {
-( 0.5 + xi )^( q-1 ) * ( 0.5 - xi )^( p-1 ) / beta(q,p)
} else {
1e16
} # end of if else 'xi between -0.5 and 0.5' stmts.
} # end of 'shapePriorBeta' function.
shapePriorBetaOld <- function(x, p, q) {
##
## Function to calculate the GMLE prior for the shape parameter.
##
# if(!is.element("shape", colnames(x))) stop("shapePriorBeta: invalid model for this prior penalty function.")
# xi <- 0.5 + x[,"shape"]
xi <- x[ length(x) ]
# res <- -sum(dbeta(x=xi, shape1=q, shape2=p, log=TRUE), na.rm=TRUE)
res <- - dbeta(x = xi, shape1 = q, shape2 = p, log = TRUE)
if(!is.finite(res)) res <- sign(res) * 1e16
return(res)
} # end of 'shapePriorBetaOld' function.
initializer <- function(x, model, threshold, ...) {
UseMethod("initializer", x)
} # end of 'initializer' function.
initializer.lmoments <- function(x, model, threshold, ..., npy, blocks=NULL) {
if(is.element(model, c("gev", "weibull", "frechet", "gumbel"))) {
lambda <- Lmoments(x)
tau3 <- lambda[3] / lambda[2]
co <- 2 / (3 + tau3) - log(2) / log(3)
kappa <- 7.8590 * co + 2.9554 * co^2
g <- gamma(1 + kappa)
sigma <- (lambda[2] * kappa)/((1 - 2^(-kappa)) * g)
mu <- lambda[1] - (sigma / kappa) * (1 - g)
xi <- -kappa
res <- c(mu, sigma, xi)
names(res) <- c("location", "scale", "shape")
} else if(is.element(model, c("gp", "exponential", "beta", "pareto", "pp"))) {
n <- length(x)
u <- threshold
if(length(u)==1) u <- rep(u, n)
id <- x > u
lambda <- Lmoments(x[id] - u[id])
tau2 <- lambda[2] / lambda[1]
sigma <- lambda[1] * (1 / tau2 - 1)
kappa <- 1 / tau2 - 2
xi <- -kappa
if(model=="pp") {
if(is.null(blocks)) {
rate <- npy * sum(id)/n
} else {
rate <- sum(id)/(blocks$nBlocks * mean(blocks$weights))
}
sigma <- exp(log(sigma) + xi*log(rate))
mu <- mean(threshold) - (sigma/xi)*(rate^(-xi) - 1)
res <- c(mu, sigma, xi, rate)
names(res) <- c("location", "scale", "shape", "rate")
} else if(model != "exponential") {
res <- c(sigma, xi)
names(res) <- c("scale", "shape")
} else {
res <- sigma
names(res) <- c("scale")
}
} # end of 'which models to estimate' stmts.
return(res)
} # end of 'initializer.lmoments' function.
initializer.moms <- function(x, model, threshold, ..., npy, blocks=NULL) {
if(is.element(model, c("gev","gumbel","frechet","weibull"))) {
m <- mean(x)
s <- sqrt(6*var(x))/pi
mu <- m - 0.57722*s
sigma <- log(s)
if(model=="weibull") xi <- -1e-8
else if(is.element(model,c("gev","frechet"))) xi <- 1e-8
else xi <- 0
res <- c(mu, sigma, xi)
names(res) <- c("location", "scale", "shape")
} else if(is.element(model, c("pp", "gp", "exponential", "beta", "pareto"))) {
n <- length(x)
u <- threshold
if(length(u)==1) u <- rep(u,n)
id <- x > u
x.u <- x[id] - u[id]
m <- mean(x.u)
sigma <- sqrt(var(x.u))
if(model=="beta") xi <- -1e-8
else if(is.element(model, c("pp","gp","beta"))) xi <- 1e-8
else xi <- 0
if(model=="pp") {
if(is.null(blocks)) {
rate <- npy * sum(id)/n
} else {
rate <- sum(id)/(blocks$nBlocks * mean(blocks$weights))
}
sigma <- exp(log(sigma) + xi*rate)
mu <- mean(u) - (sigma/xi)*(rate^(-xi) - 1)
res <- c(mu, sigma, xi)
names(res) <- c("location", "scale", "shape")
} else if(model != "exponential") {
res <- c(sigma, xi)
names(res) <- c("scale", "shape")
} else {
res <- sigma
names(res) <- "scale"
}
} else stop("initializer: invalid model argument.")
return(res)
} # end of 'initializer.moms' function.
initializer.mle <- function(x, model, threshold, ..., data=NULL, u, u.fun, loc.fun, sc.fun, sh.fun, use.phi, type, span, time.units, period.basis, blocks=NULL) {
attributes(x) <- NULL
if(is.null(data)) {
if(missing(span)) res <- fevd(x=x, threshold=u, use.phi=use.phi, type=type, time.units=time.units, period.basis=period.basis, blocks=blocks)
else res <- fevd(x=x, threshold=u, use.phi=use.phi, type=type, span=span, time.units=time.units, period.basis=period.basis, blocks=blocks)
} else {
if(missing(span)) res <- fevd(x=x, data=data, threshold=u, threshold.fun=u.fun, location.fun=loc.fun, scale.fun=sc.fun, shape.fun=sh.fun,
use.phi=use.phi, type=type, time.units=time.units, period.basis=period.basis, blocks=blocks)
else res <- fevd(x=x, data=data, threshold=u, threshold.fun=u.fun, location.fun=loc.fun, scale.fun=sc.fun, shape.fun=sh.fun,
use.phi=use.phi, type=type, span=span, time.units=time.units, period.basis=period.basis, blocks=blocks)
}
return(res$results$par)
} # end of 'initializer.bayesian' function.
will.accept <- function( ll.i, prior.i, ll.star, prior.star, log=TRUE) {
##
## Function to determine whether to accept a proposed parameter value or not.
##
## Arguments:
##
## 'll.i' numeric giving the log-likelihood value for the i-th step parameter combination.
## 'prior.i' numeric giving the (log) prior probability for the i-th step parameter combination.
## 'll.star', 'prior.star' are the log-likelihood and (log) prior values for the proposed parameter vector.
##
## Details: Calculates the ratio R=(L(theta*;data)p(theta*))/(L(theta.i;data)p(theta.i)), and determines
## whether to accept or reject the proposed parameters with probability alpha = min(1, R).
##
## Value: list with components:
## 'accept' a logical saying whether to accept the proposed parameters (if TRUE) or not (if FALSE),
## 'ratio' numeric giving the value of R (see Details above),
## 'alpha' numeric giving the value of 'alpha' (see Details above),
## 'll.star' numeric giving the log-likelihood of the proposed parameter combination,
## 'prior.star' the prior probability for the proposed parameter combination.
##
out <- list()
ratio <- (ll.star + prior.star) - (ll.i + prior.i)
if(!log) ratio <- exp( ratio)
if( is.na( ratio)) {
out$accept <- FALSE
out$ratio <- out$alpha <- NA
return( out)
}
if(!log) alpha <- min( 1, ratio)
else alpha <- min(0, ratio)
if(!log) {
if( runif( 1) < alpha) accept <- TRUE
else accept <- FALSE
} else {
if(log(runif(1)) < alpha) accept <- TRUE
else accept <- FALSE
}
out$accept <- accept
out$ratio <- ratio
out$alpha <- alpha
return( out)
} # end of 'will.accept' function.
fevdPriorDefault <- function(theta, m, v, log=TRUE) {
np <- length(theta)
dfun <- function(th) dnorm(th[1], mean=th[2], sd=th[3], log=log)
if(length(m)==1) m <- rep(m, np)
else if(length(m) != np) stop("fevdPriorDefault: m must have length 1 or same as number of parameters.")
if(length(v)==1) v <- rep(v, np)
else if(length(v) != np) stop("fevdPriorDefault: v must have length 1 or same as number of parameters.")
th <- cbind(theta, m, v)
res <- apply(th, 1, dfun)
if(log) return(sum(res))
else return(prod(res))
} # end of 'fevdPriorDefault' function.
fevdProposalDefault <- function(p, ind, ...) {
a <- list(...)
if(!is.null(a$mean)) {
#if(length(a$mean) > 1) m <- a$mean[ind]
#else m <- a$mean
m <- a$mean
} else m <- 0
if(!is.null(a$sd)) {
# if(length(a$sd) > 1) s <- a$sd[ind]
# else s <- a$sd
s <- a$sd
} else s <- 0.1
return(p + rnorm(length(p), mean=m, sd=s))
} # end of 'fevdProposalDefault' function.
fpois <- function(x, na.action = na.fail, ...) {
UseMethod("fpois", x)
} # end of 'fpois' function.
fpois.default <- function(x, na.action = na.fail, ...) {
theCall <- match.call()
dname <- deparse(substitute(x))
x <- na.action(x)
n <- length(x)
m <- mean(x)
v <- var(x)
PARAMETER <- c(m, v, n - 1)
names(PARAMETER) <- c("mean", "variance", "degrees of freedom")
CRITVAL <- (n - 1) * v / m
names(CRITVAL) <- "Chi-square(n - 1)"
PVAL <- pchisq(CRITVAL, df = n - 1, lower.tail = FALSE)
structure(list(statistic = CRITVAL, parameter = PARAMETER,
alternative = "greater", p.value = PVAL, method = "Test for Equality of (Poisson) Mean and Variance",
data.name = dname, call = theCall), class = "htest")
} # end of 'fpois.default' function.
fpois.data.frame <- function(x, na.action = na.fail, ..., which.col = 1) {
theCall <- match.call()
dname <- deparse(substitute(x))
if(length(which.col) > 1) stop("fpois: invalid which.col argument. Length must be one. To include covariates, use glm with poisson() family link.")
x <- x[, which.col]
x <- na.action(x)
n <- length(x)
m <- mean(x)
v <- var(x)
PARAMETER <- c(m, v, n - 1)
names(PARAMETER) <- c("mean", "variance", "degrees of freedom")
CRITVAL <- (n - 1) * v / m
names(CRITVAL) <- "Chi-square(n - 1)"
PVAL <- pchisq(CRITVAL, df = n - 1, lower.tail = FALSE)
structure(list(statistic = CRITVAL, parameter = PARAMETER,
alternative = "greater", p.value = PVAL, method = "Test for Equality of (Poisson) Mean and Variance",
data.name = dname, call = theCall), class = "htest")
} # end of 'fpois.data.frame' function.
fpois.matrix <- function(x, na.action = na.fail, ..., which.col = 1) {
theCall <- match.call()
dname <- deparse(substitute(x))
if(length(which.col) > 1) stop("fpois: invalid which.col argument. Length must be one. To include covariates, use glm with poisson() family link.")
x <- x[, which.col]
x <- na.action(x)
n <- length(x)
m <- mean(x)
v <- var(x)
PARAMETER <- c(m, v, n - 1)
names(PARAMETER) <- c("mean", "variance", "degrees of freedom")
CRITVAL <- (n - 1) * v / m
names(CRITVAL) <- "Chi-square(n - 1)"
PVAL <- pchisq(CRITVAL, df = n - 1, lower.tail = FALSE)
structure(list(statistic = CRITVAL, parameter = PARAMETER,
alternative = "greater", p.value = PVAL, method = "Test for Equality of (Poisson) Mean and Variance",
data.name = dname, call = theCall), class = "htest")
} # end of 'fpois.matrix' function.
fpois.list <- function(x, na.action = na.fail, ..., which.component = 1) {
theCall <- match.call()
dname <- deparse(substitute(x))
if(length(which.component) > 1) stop("fpois: invalid which.component argument. Length must be one. To include covariates, use glm with poisson() family link.")
x <- x[[which.component]]
x <- na.action(x)
n <- length(x)
m <- mean(x)
v <- var(x)
PARAMETER <- c(m, v, n - 1)
names(PARAMETER) <- c("mean", "variance", "degrees of freedom")
CRITVAL <- (n - 1) * v / m
names(CRITVAL) <- "Chi-square(n - 1)"
PVAL <- pchisq(CRITVAL, df = n - 1, lower.tail = FALSE)
structure(list(statistic = CRITVAL, parameter = PARAMETER,
alternative = "greater", p.value = PVAL, method = "Test for Equality of (Poisson) Mean and Variance",
data.name = dname, call = theCall), class = "htest")
} # end of 'fpois.list' function.
Try the extRemes package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.