Nothing
lr.test <- function(x, y, alpha=0.05, df=1, ...) {
cx <- class( x )
if( "fevd" %in% cx ) {
if(!is.element(x$method, c("MLE", "GMLE"))) stop("lr.test: fit method must be MLE or GMLE")
l1 <- x$results$value
df1 <- length(x$results$par)
dname <- x$data.name[1]
} else if(is.numeric(x) && length(x)==1) l1 <- x
else if(is.numeric(x) && length(x)==2) l1 <- x[1]
else stop("lr.test: invalid x argument. Must be a single number, length two numeric vector, or an fevd object.")
if( !("fevd" %in% cx ) ) dname <- deparse(substitute(x))
cy <- class( y )
if( "fevd" %in% cy ) {
if(!is.element(y$method, c("MLE", "GMLE"))) stop("lr.test: fit method must be MLE or GMLE")
l2 <- y$results$value
df2 <- length(y$results$par)
dname <- c(dname, y$data.name[1])
} else if(is.numeric(y) && length(y) == 1) l2 <- y
else stop("lr.test: invalid y argument. Must be a single number or an fevd object.")
if(("fevd" %in% cx) && ( "fevd" %in% cy ) ) {
if(df2 < df1) return(lr.test(x=y, y=x, alpha=alpha, df=df, ...))
else df <- df2 - df1
}
if(!("fevd" %in% cy ) ) dname <- c(dname, deparse(substitute(y)))
names(dname) <- c("x", "y")
STATISTIC <- -2*(l2 - l1)
names(STATISTIC) <- "Likelihood-ratio"
CRITVAL <- qchisq(alpha, df=df, lower.tail=FALSE)
PVAL <- pchisq(STATISTIC, df=df, lower.tail=FALSE)
PARAMETER <- c(CRITVAL, alpha, df)
names(PARAMETER) <- c("chi-square critical value", "alpha", "Degrees of Freedom")
structure(list(statistic=STATISTIC, parameter=PARAMETER, alternative="greater", p.value=PVAL, method="Likelihood-ratio Test", data.name=dname), class="htest")
} # end of 'lr.test' function.
plot.fevd <- function(x, type = c("primary", "probprob", "qq", "qq2", "Zplot", "hist", "density", "rl", "trace"),
rperiods = c(2, 5, 10, 20, 50, 80, 100, 120, 200, 250, 300, 500, 800), a=0, hist.args=NULL, density.args=NULL, d = NULL, ...) {
# x2 <- x
# if(is.element(x$method, c("MLE", "GMLE"))) class(x2) <- c( "fevd.mle", class( x2 ) )
# else class( x2 ) <- c( paste( "fevd.", tolower( x$method ), sep="" ), class( x2 ) )
# UseMethod("plot", x2)
# UseMethod( "plot", x )
newcl <- if(is.element(x$method, c("MLE", "GMLE")))
"fevd.mle"
else
paste( "fevd.", tolower( x$method ), sep="" )
get(paste0("plot.", newcl))(x = x, type = type, rperiods = rperiods,
a = a, hist.args = hist.args, density.args = density.args,
d = d, ...)
} # end of 'plot.fevd' function.
plot.fevd.lmoments <- function(x, type=c("primary", "probprob", "qq", "qq2", "Zplot", "hist", "density", "rl", "trace"),
rperiods=c(2, 5, 10, 20, 50, 80, 100, 120, 200, 250, 300, 500, 800), a=0, hist.args=NULL, density.args=NULL, d = NULL, ...) {
# type <- tolower(type)
type <- match.arg(type)
type <- tolower(type)
if( type == "trace" ) {
cat("\n", "No trace plot available for L-moments estimation.\n")
invisible()
}
args <- list(...)
model <- x$type
p <- x$results
pnames <- names(p)
if(is.element("location",pnames)) loc <- p["location"]
else loc <- 0
scale <- p["scale"]
if(is.element("shape",pnames)) shape <- p["shape"]
else shape <- 0
if(!is.null(x$threshold)) u <- x$threshold
else u <- 0
y <- datagrabber(x, cov.data=FALSE)
if(is.element(model, c("PP","GP","Exponential","Beta","Pareto"))) eid <- y > u
else eid <- !logical(x$n)
m <- sum(eid)
if(is.element(type, c("primary","probprob","qq","rl"))) xp <- ppoints(m, a = a)
if(type=="probprob") out <- probprob.plot.evd(xp=xp, y=y, model=model, loc=loc, scale=scale, shape=shape, u=u, tform=FALSE, eid=eid, obj=x, npy=x$npy, ...)
if( type == "primary" ) {
op <- par()
if(model != "PP") par(mfrow=c(2,2), oma=c(0,0,2,0))
else par(mfrow = c(3, 2), oma = c(0,0,2,0))
}
if(is.element(type,c("primary","qq"))) out <- quantquant.plot.evd(x=x, xp=xp, y=y, u=u, loc=loc, scale=scale, shape=shape, tform=FALSE, eid=eid, model=model, type=type, ...)
if(is.element(type,c("primary","qq2"))) out <- quantquant2.plot.evd(x=x, y=y, eid=eid, model=model, type=type)
if(type=="hist") out <- histplot.evd(x=x, y=y, u=u, loc=loc, scale=scale, shape=shape, eid=eid, model=model, hist.args=hist.args, ...)
if(is.element(type, c("primary","density"))) out <- densplot.evd(x=x, y=y, u=u, loc=loc, scale=scale, shape=shape, eid=eid,
model=model, tform=FALSE, density.args=density.args, type=type, ...)
if(is.element(type, c("primary", "Zplot")) && model == "PP") {
if(type == "primary") {
eeplot(x = x, type = "Zplot", main = "Zplot", d = d, ...)
} else out <- eeplot(x = x, type = "Zplot", d = d, ...)
}
if(is.element(type, c("primary","rl"))) out <- rlplot.evd(x=x, xp=xp, y=y, u=u, eid=eid, rperiods=rperiods, tform=FALSE, model=model, type=type, a=a, ...)
if( type == "primary" ) {
mtext(deparse(x$call), line=0.5, outer=TRUE)
par( mfrow = op$mfrow, oma = op$oma )
} # end of if 'type' is primary stmt.
if( type != "primary" ) invisible( out )
else invisible()
} # end of 'plot.fevd.lmoments' function.
plot.fevd.bayesian <- function(x, type=c("primary", "probprob", "qq", "qq2", "Zplot", "hist", "density", "rl", "trace"),
rperiods=c(2, 5, 10, 20, 50, 80, 100, 120, 200, 250, 300, 500, 800), a=0, hist.args=NULL, density.args=NULL,
burn.in=499, d = NULL, ...) {
type <- match.arg(type)
if(type != "Zplot") type <- tolower(type)
model <- x$type
p <- p2 <- x$results
np <- ncol(p) - 1
iters <- nrow(p)
p <- p2 <- p[,1:np]
pnames <- colnames(p)
if(is.element("log.scale",pnames)) {
id <- pnames == "log.scale"
p[,id] <- exp(p[,id])
pnames[id] <- "scale"
colnames(p) <- pnames
}
p <- apply(p[-(1:burn.in), ], 2, mean, na.rm=TRUE)
tform <- !is.fixedfevd(x)
if(tform && is.element(model, c("PP", "GP", "Beta", "Exponential", "Pareto")) && is.element(type, c("hist", "density")))
stop("plot.fevd: invalid type argument for this model.")
if(!tform) {
if(is.element(model, c("PP","GEV","Gumbel","Weibull","Frechet"))) {
loc <- p["location"]
nloc <- 1
} else loc <- nloc <- 0
scale <- p["scale"]
nsc <- 1
if(!is.element(model, c("Gumbel","Exponential"))) {
shape <- p["shape"]
nsh <- 1
} else nsh <- shape <- 0
ytrans <- NULL
} else {
ytrans <- trans(x)
if(model == "PP") ytransPP <- trans(x, return.all = TRUE)
designs <- setup.design(x)
if(is.element(model, c("PP","GEV","Gumbel","Weibull","Frechet"))) {
X.loc <- designs$X.loc
nloc <- ncol(X.loc)
loc <- rowSums(matrix(p[1:nloc], x$n, nloc, byrow=TRUE) * X.loc)
} else {
loc <- 0
nloc <- 0
}
X.sc <- designs$X.sc
nsc <- ncol(X.sc)
scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], x$n, nsc, byrow=TRUE) * X.sc)
if(x$par.models$log.scale) scale <- exp(scale)
if(!is.element(model, c("Gumbel","Exponential"))) {
X.sh <- designs$X.sh
nsh <- ncol(X.sh)
shape <- rowSums(matrix(p[(nloc+nsc+1):np], x$n, nsh, byrow=TRUE) * X.sh)
} else {
shape <- 0
nsh <- 0
}
if(is.element(type, c("primary","return.level"))) if(missing(rperiods)) rperiods <- c(2, 20, 100)
} # end of if '!tform' stmts.
if(!is.null(x$threshold)) u <- x$threshold
else u <- 0
y <- datagrabber(x, cov.data=FALSE)
if(model == "PP" && is.element(type, c("primary", "hist", "density"))) {
blocks <- rep(1:round(x$span), each=round(x$npy))
n2 <- length(blocks)
if(n2 < x$n) blocks <- c(blocks, rep(blocks[n2], x$n - n2))
else if(n2 > x$n) blocks <- blocks[1:x$n]
ybm <- c(aggregate(y, by=list(blocks), max)$x)
}
if(is.element(model, c("PP","GP","Exponential","Beta","Pareto"))) eid <- y > u
else eid <- !logical(x$n)
m <- sum(eid)
if(is.element(type, c("primary","probprob","qq","rl"))) xp <- ppoints(m, a = a)
if( type == "probprob" ) out <- probprob.plot.evd(xp=xp, y=y, model=model, loc=loc, scale=scale, shape=shape, u=u, tform=tform, eid=eid, obj=x, ytrans=ytrans, npy=x$npy, ...)
if(type=="primary") {
op <- par()
if(model != "PP") par(mfrow=c(2,2), oma=c(0,0,2,0))
else if(!tform) par(mfrow = c(3,2), oma = c(0,0,2,0))
else par(mfrow = c(2,2), oma = c(0,0,2,0))
} else if(type=="trace") {
op <- par()
par(mfcol=c(2,np), oma=c(0,0,2,0))
}
if(is.element(type,c("primary","qq"))) out <- quantquant.plot.evd(x=x, xp=xp, y=y, u=u, loc=loc, scale=scale, shape=shape,
tform=tform, ytrans=ytrans, eid=eid, model=model, type=type, ...)
if(is.element(type,c("primary","qq2"))) out <- quantquant2.plot.evd(x=x, y=y, eid=eid, model=model, type=type)
if(type=="hist") {
if(model != "PP") out <- histplot.evd(x=x, y=y, u=u, loc=loc, scale=scale, shape=shape, ytrans=ytrans,
tform=tform, eid=eid, model=model, hist.args=hist.args, ...)
else if(!tform) out <- histplot.evd(x=x, y=ybm, u=u, loc=loc, scale=scale, shape=shape, ytrans=ytrans,
tform=tform, eid=eid, model=model, hist.args=hist.args, ...)
}
if(is.element(type, c("primary", "Zplot")) && model == "PP") {
if(type == "primary") {
eeplot(x = x, type = "Zplot", main = "Z plot", d = d, ...)
} else out <- eeplot(x = x, type = type, d = d, ...)
}
if(type == "primary" && tform && is.element(model, c("PP", "GP", "Beta", "Exponential", "Pareto"))) {
## Eric -- At some point, should change ths to be a little cleaner, but ok for now.
## Just do not want to plot the density plot if it is a non-stationary POT model.
} else if(is.element(type, c("primary","density"))) {
if(model != "PP") out <- densplot.evd(x=x, y=y, u=u, loc=loc, scale=scale, shape=shape, tform=tform, eid=eid, ytrans=ytrans,
model=model, density.args=density.args, type=type, ...)
else if(!tform) out <- densplot.evd(x=x, y=y, u=u, loc=loc, scale=scale, shape=shape, tform=tform, eid=eid, ytrans=ybm,
model=model, density.args=density.args, type=type, ...)
}
if(is.element(type, c("primary","rl"))) {
# Eric 8/27/13 -- make sure for POT models that if the threshold varies, so do parameters, or else don't do rl plot.
do.rlplot <- TRUE
if(is.element(model, c("PP", "GP", "Exponential", "Beta", "Pareto"))) {
const.thresh <- check.constant(x$par.models$threshold)
const.loc <- check.constant(x$par.models$location)
const.scale <- check.constant(x$par.models$scale)
const.shape <- check.constant(x$par.models$shape)
if(!const.thresh && all(c(const.loc, const.scale, const.shape))) do.rlplot <- FALSE
if(!do.rlplot && type == "rl")
stop("plot.fevd: invalid type argument for POT models with varying thresholds but constant parameters (are you sure you about this model choice?).")
}
if(do.rlplot) out <- rlplot.evd(x=x, xp=xp, y=y, u=u, eid=eid, rperiods=rperiods, tform=tform, model=model, type=type, a=a, ...)
}
if(type=="trace") {
pnames2 <- colnames(p2)
id <- pnames2 == "log.scale"
if(any(id)) p2[,id] <- exp(p2[,id])
msg <- paste("Posterior Density\n", pnames, sep="")
for(i in 1:np) {
out1 <- list()
if(is.null(density.args)) {
if(i==1) plot( out1[[ i ]] <- density(p2[-(1:burn.in),i]), main=msg[i], ... )
else plot( out1[[ i ]] <- density(p2[-(1:burn.in),i]), ylab="", main=msg[i], ... )
} else {
yd <- do.call("density", c(list(x=p2[-(1:burn.in),i]), density.args))
if(i==1) plot(yd, main=msg[i], ...)
else plot(yd, ylab="", main=msg[i], ...)
out1[[ i ]] <- yd
}
xt <- 1:iters
if(i==1) plot(xt, p2[,i], type="l", xlab=pnames[i], ylab="MCMC trace", ...)
else plot(xt, p2[,i], type="l", xlab=pnames[i], ylab="", ...)
abline(h=p[i], v=burn.in, lty=2, col="gray")
# xt2 <- xt[(burn.in+1):iters]
# mt <- (cumsum(p2[(burn.in+1):iters,i])/xt2) - p[i]
# par(usr=c(range(xt), 1.5 * range(mt,finite=TRUE)))
# lines(xt2, mt, col="darkorange", lty=3)
# axis(4, at=pretty(mt), labels=pretty(mt), col="darkorange", col.ticks="darkorange")
# legend("topleft", legend=c("trace","posterior mean","cumsum diff"), col=c("black","gray","darkorange"), lty=1:3, bty="n")
} # end of for 'i' loop.
out2 <- data.frame( xt, p2 )
colnames( out2 ) <- c( "iteration", pnames )
out <- list( posterior.density = out1, MCMC.trace = out2 )
} # end of if 'trace' stmts.
if(is.element(type, c("primary","trace"))) {
mtext(deparse(x$call), line=0.5, outer=TRUE)
par( mfrow = op$mfrow, mfcol = op$mfcol, oma = op$oma )
}
if( type != "primary" ) invisible( out )
else invisible()
} # end of 'plot.fevd.bayesian' function.
plot.fevd.mle <- function(x, type=c("primary", "probprob", "qq", "qq2", "Zplot", "hist", "density", "rl", "trace"),
rperiods=c(2, 5, 10, 20, 50, 80, 100, 120, 200, 250, 300, 500, 800), a=0, hist.args=NULL, density.args=NULL, period="year",
prange=NULL, d = NULL, ...) {
type <- match.arg(type)
model <- x$type
pars <- x$results$par
args <- list(...)
# Get response and any covariate data sets.
ytmp <- datagrabber(x)
if(x$data.name[2] != "") data <- ytmp[,-1]
else data <- NULL
y <- c(ytmp[,1])
const.thresh <- check.constant(x$par.models$threshold)
const.loc <- check.constant(x$par.models$location)
const.scale <- check.constant(x$par.models$scale)
const.shape <- check.constant(x$par.models$shape)
# Eric -- 8/27/13 -- Don't allow rl plot if model is POT, threshold varies, but parameters are all constant.
if(is.element(model, c("PP", "GP", "Exponential", "Beta", "Pareto")) && !const.thresh && all(c(const.loc, const.scale, const.shape)) && type == "rl")
stop("plot.fevd: invalid type argument for POT models with varying thresholds but constant parameters (are you sure about this model choice?).")
tform <- !is.fixedfevd(x)
if(tform) {
ytrans <- trans(x)
if(model == "PP" && is.element(type, c("primary", "hist", "density"))) ytransPP <- trans(x, return.all = TRUE)
}
if(!tform) {
if(!is.element(model, c("GP","Beta","Pareto","Exponential"))) loc <- pars["location"]
else loc <- NULL
if(is.element("scale", names(pars))) scale <- pars["scale"]
else scale <- exp(pars["log.scale"])
if(!is.element(model, c("Gumbel","Exponential"))) shape <- pars["shape"]
else shape <- 0
} else if(is.element(type, c("primary", "rl"))) if(missing(rperiods)) rperiods <- c(2, 20, 100)
# end of if else '!tform' stmts.
npy <- x$npy
if(is.element(model, c("PP","GP","Beta","Pareto","Exponential"))) {
u <- x$threshold
eid <- y > u
lam <- mean(eid)
} else u <- lam <- NULL
if(is.element(type,c("primary","probprob","qq", "rl"))) {
if(is.element(model, c("GP","PP","Beta","Pareto","Exponential"))) {
if(!tform) n <- sum(y > x$threshold)
else n <- sum(!is.na(ytrans) & !is.nan(ytrans))
} else n <- x$n
xp <- ppoints(n = n, a = a)
}
if(type=="primary") {
op <- par()
# if(is.element(model, c("GP", "Exponential", "Beta", "Pareto")) && tform) par(mfrow = c(1, 2), oma = c(0,0,2,0))
if(model != "PP") par(mfrow=c(2,2), oma=c(0,0,2,0))
else if(!tform) par(mfrow = c(3, 2), oma = c(0,0,2,0))
else par(mfrow = c(2,2), oma = c(0,0,2,0))
}
if(type=="probprob") {
if(is.null(args$main)) m1 <- deparse(x$call)
if(!tform) {
if(is.element(model, c("PP","GP","Beta","Pareto","Exponential"))) {
yp <- pevd(sort(y[eid]), loc=loc, scale=scale, shape=shape, threshold=u, npy=npy, type=model)
} else yp <- pevd(sort(y), loc=loc, scale=scale, shape=shape, threshold=u, npy=npy, type=model)
if(is.null(args$main)) plot(xp, yp, main=m1, xlab="Empirical Probabilities", ylab="Model Probabilities", ...)
else plot(xp, yp, xlab="Empirical Probabilities", ylab="Model Probabilities", ...)
abline(0,1)
} else {
# yp <- pevd(sort(ytrans), loc=0, scale=1, shape=0, threshold=0, npy=npy, type=model)
if(is.element(model, c("GEV","Gumbel","Weibull","Frechet"))) yp <- exp(-exp(-sort(ytrans)))
else if(model=="PP") yp <- sort(ytrans)
else yp <- 1 - exp(-sort(ytrans))
if(is.null(args$main)) {
plot(xp, yp, main=m1, xlab="Residual Empirical Probabilities", ylab="Residual Model Probabilities", ...)
} else plot(xp, yp, xlab="Residual Empirical Probabilities", ylab="Residual Model Probabilities", ...)
abline(0,1)
}
out <- data.frame( empirical = xp, model = yp )
} # end of 'probprob' stmts.
if(is.element(type,c("primary","qq"))) {
if(!tform) {
if(is.element(model, c("Weibull","Frechet"))) mod2 <- "GEV"
else if(is.element(model, c("Exponential","Beta","Pareto"))) mod2 <- "GP"
else mod2 <- model
if(!is.element(model,c("PP","GP","Beta","Pareto","Exponential"))) yq <- qevd(xp, loc=loc, scale=scale, shape=shape, type=mod2)
else yq <- qevd(xp, threshold=u, loc=loc, scale=scale, shape=shape, type=mod2)
if(is.null(args$main)) {
if(type=="primary") m2 <- ""
else m2 <- deparse(x$call)
if(is.element(model, c("GEV","Weibull","Frechet","Gumbel"))) {
plot(yq, sort(y), xlab="Model Quantiles", ylab="Empirical Quantiles", main=m2)
if( type == "qq" ) out <- data.frame( empirical = sort( y ), model = yq )
} else {
plot(yq, sort(y[eid]), xlab="Model Quantiles", ylab="Empirical Quantiles", main=m2, ...)
if( type == "qq" ) out <- data.frame( empirical = sort( y[ eid ] ), model = yq )
}
} else {
if(is.element(model, c("GEV","Weibull","Frechet","Gumbel"))) {
plot(yq, sort(y), xlab="Model Quantiles", ylab="Empirical Quantiles", ...)
if( type == "qq" ) out <- data.frame( empirical = sort( y ), model = yq )
} else {
plot(yq, sort(y[eid]), xlab="Model Quantiles", ylab="Empirical Quantiles", ...)
if( type == "qq" ) out <- data.frame( empirical = sort( y[ eid ] ), model = yq )
}
}
} else {
if(is.null(args$main)) {
if(type=="primary") m2 <- ""
else m2 <- deparse(x$call)
if(is.element(model, c("GEV","Weibull","Frechet"))) m2 <- paste(m2, "(Gumbel Scale)", sep="\n")
else if(is.element(model, c("PP", "GP", "Beta", "Pareto"))) m2 <- paste(m2, "Exponential Scale", sep="\n")
if(is.element(model, c("GEV","Weibull","Gumbel","Frechet"))) {
plot(-log(-log(sort(xp))), sort(ytrans), main=m2, xlab="(Standardized) Model Quantiles",
ylab="Empirical Residual Quantiles", ...)
if( type == "qq" ) out <- data.frame( empirical = sort( ytrans ), model = -log(-log(sort(xp))) )
} else if(is.element(model, c("GP","Beta","Exponential","Pareto"))) {
plot(-log(1 - xp), sort(ytrans), main=m2, xlab="(Standardized) Residual Quantiles",
ylab="Empirical Residual Quantiles", ...)
if( type == "qq" ) out <- data.frame( empirical = sort( ytrans ), model = -log(1 - xp) )
} else if(model=="PP") {
plot(-log(1 - xp), sort(-log(ytrans)), main=m2, xlab="(Standardized) Residual Quantiles",
ylab="Empirical Residual Quantiles", ...)
if( type == "qq" ) out <- data.frame( empirical = sort(-log(ytrans)), model = -log(1 - xp) )
}
} else {
if(is.element(model, c("GEV","Weibull","Gumbel","Frechet"))) {
plot(-log(-log(sort(xp))), sort(ytrans), xlab="Model", ylab="Empirical", ...)
if( type == "qq" ) out <- data.frame( empirical = sort( ytrans ), model = -log(-log(sort(xp))) )
} else if(is.element(model, c("GP","Beta","Exponential","Pareto"))) {
plot(-log(1 - xp), sort(ytrans), xlab="Model", ylab="Empirical", ...)
if( type == "qq" ) out <- data.frame( empirical = sort( ytrans ), model = -log(1 - xp) )
} else if(model=="PP") {
plot(-log(1 - xp), sort(-log(ytrans)), xlab="Model", ylab="Empirical", ...)
if( type == "qq" ) out <- data.frame( empirical = sort( -log( ytrans ) ), model = -log(1 - xp) )
}
}
} # end of if else 'tform' stmts.
abline(0,1)
} # end of if 'qq' stmts.
if(is.element(type, c("primary","qq2"))) {
if(is.fixedfevd(x) && !is.null(x$blocks)) {
z <- rextRemes(x, round(x$blocks$nBlocks * x$npy))
# CJP: for stationary, it generates based on number of obs, so need to tell it how many obs there would be
} else {
if(!is.element(model, c("PP","GP","Beta","Exponential","Pareto"))) z <- rextRemes(x)
else z <- rextRemes( x, n = sum( eid, na.rm = TRUE ) )
# for non-stationary, it generates based on number of exceedances, so this is fine
}
if(!is.element(model, c("PP","GP","Beta","Exponential","Pareto"))) {
yQQ <- y
if(is.null(args$xlab)) xl <- paste(x$data.name[1], " Empirical Quantiles", sep="")
else xl <- args$xlab
} else {
yQQ <- y[eid]
if(is.null(args$xlab)) {
if(length(x$threshold)==1) xl <- paste(x$data.name[1], "( > ", x$threshold, ") Empirical Quantiles", sep="")
else xl <- paste(x$data.name[1], "( > threshold) Empirical Quantiles", sep="")
} else xl <- args$xlab
}
if(!is.null(args$main)) {
if(type=="primary") mQQ <- ""
else mQQ <- deparse(x$call)
if( type != "qq2" ) qqplot(yQQ, z, main=mQQ, xlab=xl, ylab="Quantiles from Model Simulated Data")
else out <- qqplot(yQQ, z, main=mQQ, xlab=xl, ylab="Quantiles from Model Simulated Data")
} else {
if( type != "qq2" ) qqplot(yQQ, z, xlab=xl, ylab="Quantiles from Model Simulated Data")
else out <- qqplot(yQQ, z, xlab=xl, ylab="Quantiles from Model Simulated Data")
}
} # end of if do second qq-plot using simulated data from the model stmts.
if(!(type == "primary" && model == "PP" && tform)) {
if(is.element(type, c("primary","density")) && is.null(x$blocks)) {
# CJP -- this needs the full data vector (or perhaps more extensive rewriting to work with blocks)
if(!tform) {
if(!is.element(model, c("PP","GP","Beta","Exponential","Pareto"))) yd <- y
else if(model != "PP") yd <- y[eid] - x$threshold
else {
if(x$span %% 1 != 0 || x$npy %% 1 != 0)
warning("plot.fevd.mle: span or npy not integers; determination of max in each block may be substantially in error when there are many blocks.")
blocks <- rep(1:round(x$span), each=round(x$npy))
# CJP2: hard to deal exactly with inhomogeneity in number of values in each block;
# hopefully rounding here is reasonable in most cases but with many blocks,
# the indexing could be substantially off.
n2 <- length(blocks)
if(n2 < x$n) blocks <- c(blocks, rep(blocks[n2], x$n - n2))
else if(n2 > x$n) blocks <- blocks[1:x$n]
yd <- c(aggregate(y, by=list(blocks), max)$x)
}
if(is.null(density.args)) yd <- density(yd)
else yd <- do.call("density", c(list(yd), density.args))
if(is.element(type, c("primary","density"))) xd <- seq(min(yd$x, na.rm=TRUE), max(yd$x, na.rm=TRUE),,100)
else xd <- seq(min(yh, na.rm=TRUE), max(yh, na.rm=TRUE),,100)
if(is.element(model, c("PP", "Gumbel", "Weibull", "Frechet"))) mod2 <- "GEV"
else if(is.element(model, c("Pareto", "Frechet", "Beta", "Exponential"))) mod2 <- "GP"
else mod2 <- model
if(tform) {
mu <- 0
sig <- 1
xi <- 0
u2 <- 0
} else {
mu <- loc
sig <- scale
xi <- shape
u2 <- u[1]
# if(model=="PP") sig <- sig + xi * (u2 - loc)
}
yd2 <- devd(xd, loc=mu, scale=sig, shape=xi, threshold=u2, type=mod2)
# lines(xd, yd2, lty=2, col="blue", lwd=1.5)
if(is.null(args$ylim)) {
yld <- range(c(yd$y, yd2), finite=TRUE)
yld[1] <- min(yld[1], 0)
# yld[2] <- max(yld[2]+0.5, 1)
}
if(is.null(args$main)) {
if(type=="primary") m3 <- ""
else m3 <- deparse(x$call)
if(is.null(args$ylim)) plot(yd, main=m3, ylim=yld, ...)
else plot(yd, main=m3, ...)
} else {
if(is.null(args$ylim)) plot(yd, ylim=yld, ...)
else plot(yd, ...)
}
} else {
if(type == "density" && is.element(model, c("PP", "GP", "Exponential", "Pareto", "Beta")))
stop("plot.fevd: invalid type argument for this model.")
if(!is.element(model, c("PP", "GP", "Exponential", "Pareto", "Beta"))) {
# if(model == "PP") { # Eric 8/14/13
# if(x$span %% 1 != 0 || x$npy %% 1 != 0)
# warning("plot.fevd.mle: span or npy not integers; determination of max in each block may be substantially in error when there are many blocks.")
# blocks <- rep(1:round(x$span), each=round(x$npy))
#
# # CJP2: hard to deal exactly with inhomogeneity in number of values in each block;
# # hopefully rounding here is reasonable in most cases but with many blocks,
# # the indexing could be substantially off
#
#
# n2 <- length(blocks)
# if(n2 < x$n) blocks <- c(blocks, rep(blocks[n2], x$n - n2))
# else if(n2 > x$n) blocks <- blocks[1:x$n]
#
# ytmp <- c(aggregate(y, by = list(blocks), max)$x)
#
# mbid <- logical(x$n)
#
# for(i in 1:(length(unique(blocks)))) {
#
# tmpind <- blocks == blocks[i]
# tmpind2 <- y == ytmp[i]
# tmpind2[is.na(tmpind2)] <- FALSE
#
# finind <- tmpind & tmpind2
#
# if(sum(finind) > 1) {
#
# numsind <- (1:x$n)[finind]
# numsind <- numsind[1]
# finind[-numsind] <- FALSE
#
# }
#
# mbid[finind] <- TRUE
#
# } # end of for 'i' loop.
#
# ytrans2 <- ytransPP[mbid]
# # ytrans2 <- c(aggregate(ytransPP, by=list(blocks), max)$x)
#
# } else ytrans2 <- ytrans
if(is.null(density.args)) yd <- density(ytrans)
else yd <- do.call("density", c(list(ytrans), density.args))
if(is.null(args$ylim)) {
yld <- range(yd$y, finite=TRUE)
yld[1] <- min(yld[1], 0)
# yld[2] <- max(yld[2]+0.5, 1)
}
if(is.null(args$main)) {
if(type=="primary") m3 <- "Transformed Data"
else m3 <- paste(deparse(x$call), "Transformed Data", sep="\n")
if(is.null(args$ylim)) plot(yd, main=m3, ylim=yld, ...)
else plot(yd, main=m3, ...)
} else {
if(is.null(args$ylim)) plot(yd, ylim=yld, ...)
else plot(yd, ...)
}
} # end of if model is non-stationary GEV stmts.
} # end of if else '!tform' stmts.
if( type == "density" ) out <- yd
} # end of if 'primary plots or density' stmts.
}
if(type=="hist" && is.null(x$blocks)) { # CJP2: I previously forgot to exclude this when there are blocks
if(!tform) {
if(is.null(args$main)) {
if(model != "PP") m4 <- paste(deparse(x$call), "Histogram", sep="\n")
else m4 <- paste(deparse(x$call), "\n", paste("Histogram (", x$period.basis, " maxima)", sep=""))
}
if(is.element(model, c("GP", "Beta", "Exponential", "Pareto"))) yh <- y[eid] - x$threshold
else if(model != "PP") yh <- y
else {
blocks <- rep(1:x$span, each=x$npy)
n2 <- length(blocks)
if(n2 < x$n) blocks <- c(blocks, rep(blocks[n2], x$n - n2))
else if(n2 > x$n) blocks <- blocks[1:x$n]
yh <- c(aggregate(y, by=list(blocks), max)$x)
}
} else {
if(is.element(model, c("PP", "GP", "Exponential", "Beta", "Pareto"))) stop("plot.fevd: invalid type argument for this model.")
if(is.null(args$main)) m4 <- paste(deparse(x$call), "Histogram of Transformed Data", sep="\n")
# if(model != "PP") yh <- ytrans
# else {
# blocks <- rep(1:x$span, each=x$npy)
# n2 <- length(blocks)
# if(n2 < x$n) blocks <- c(blocks, rep(blocks[n2], x$n - n2))
# else if(n2 > x$n) blocks <- blocks[1:x$n]
# yh <- c(aggregate(ytransPP, by=list(blocks), max)$x)
# } # end of if 'model != PP' stmts.
}
if(is.null(hist.args)) {
if(!is.null(args$ylim)) {
if(is.null(args$main)) {
if(is.null(args$col)) {
out <- hist(yh, col="darkblue", freq=FALSE, breaks="FD", xlab=x$data.name[1], main=m4, ...)
} else out <- hist(yh, freq=FALSE, breaks="FD", main=m4, xlab=x$data.name[1], ...)
} else {
if(is.null(args$col)) out <- hist(yh, col="darkblue", freq=FALSE, breaks="FD", xlab=x$data.name[1], ...)
else out <- hist(yh, freq=FALSE, breaks="FD", xlab=x$data.name[1], ...)
}
} else {
if(is.null(args$main)) {
if(is.null(args$col)) out <- hist(yh, col="darkblue", freq=FALSE, breaks="FD", xlab=x$data.name[1], main=m4, ylim=c(0,1.5), ...)
else out <- hist(yh, freq=FALSE, breaks="FD", main=m4, xlab=x$data.name[1], ylim=c(0,1.5), ...)
} else {
if(is.null(args$col)) out <- hist(yh, col="darkblue", freq=FALSE, breaks="FD", xlab=x$data.name[1], ylim=c(0,1.5), ...)
else out <- hist(yh, freq=FALSE, breaks="FD", xlab=x$data.name[1], ylim=c(0,1.5), ...)
}
}
} else out <- do.call("hist", c(list(yh), hist.args))
} # end of if 'type = hist' stmts.
if(!(tform && is.element(model, c("PP", "GP", "Beta", "Exponential", "Pareto")))) {
if(is.element(type, c("primary", "density", "hist")) && is.null(x$blocks)) {
# CJP: this needs the full data vector (or perhaps more extensive rewriting to work with blocks)
if(is.element(type, c("primary","density"))) xd <- seq(min(yd$x, na.rm=TRUE), max(yd$x, na.rm=TRUE),,100)
else xd <- seq(min(yh, na.rm=TRUE), max(yh, na.rm=TRUE),,100)
if(is.element(model, c("PP", "Gumbel", "Weibull", "Frechet"))) mod2 <- "GEV"
else if(is.element(model, c("Pareto", "Frechet", "Beta", "Exponential"))) mod2 <- "GP"
else mod2 <- model
if(tform) {
mu <- 0
sig <- 1
xi <- 0
u2 <- 0
} else {
mu <- loc
sig <- scale
xi <- shape
u2 <- u[1]
# if(model=="PP") sig <- sig + xi * (u2 - loc)
}
yd2 <- devd(xd, loc=mu, scale=sig, shape=xi, threshold=u2, type=mod2)
lines(xd, yd2, lty=2, col="blue", lwd=1.5)
if(model != "PP") {
if(type=="hist" || !is.null(density.args)) legend("topright", legend="Modeled Density", col="blue", lty=2, lwd=1.5, bty="n")
else if(is.null(density.args)) legend("topright", legend=c("Empirical", "Modeled"), col=c("black","blue"), lty=c(1,2), lwd=c(1, 1.5), bty="n")
} else {
if(type=="hist" || !is.null(density.args)) legend("topright", legend="Modeled Density", col="blue", lty=2, lwd=1.5, bty="n")
else if(is.null(density.args)) legend("topright", legend=c(paste("Empirical (", x$period.basis, " maxima)", sep=""), "Modeled"),
col=c("black","blue"), lty=c(1,2), lwd=c(1, 1.5), bty="n")
}
} # end of if add density lines stmts.
} # end of making sure that it is not a non-stationary POT model stmts.
# Z and W plots for all PP models.
if(is.element(type, c("primary", "Zplot")) && model == "PP") {
if(type == "primary") {
eeplot(x = x, type = "Zplot", main = "Z plot", d = d, ...)
} else out <- eeplot(x = x, type = type, d = d, ...)
}
if(type == "Zplot" && model != "PP") stop("plot.fevd.mle: invalid type argument for this model.")
if(is.element(type, c("primary", "rl")) && is.null(x$blocks)) {
# Eric 8/27/13 -- do not plot return levels if model is POT, threshold varies and parameters do not vary.
if(!(is.element(model, c("PP", "GP", "Exponential", "Beta", "Pareto")) && !const.thresh && all(c(const.loc, const.scale, const.shape)))) {
# CJP: this needs the full data vector
# (or perhaps more extensive rewriting to work with blocks,
# but would be hard to deal with blocks with no exceedances)
if(is.null(args$main)) {
if(type=="primary") m5 <- ""
else m5 <- deparse(x$call)
}
if(model=="PP") {
# if(!tform) mod2 <- "GEV"
# else mod2 <- "GP"
mod2 <- "GEV"
if(is.null(args$main)) {
# if(!tform) {
if(type=="rl") m5 <- paste(m5, "Return Levels based on approx. equivalent GEV", sep="\n")
else m5 <- paste("Return Levels based on approx.", "equivalent GEV", sep="\n")
# } else {
# if(type=="rl") m5 <- paste(m5, "Return Levels based on approx. equivalent GP df", sep="\n")
# else m5 <- paste("Return Levels based on approx.", "equivalent GP df", sep="\n")
# }
}
blocks <- rep(1:x$span, each=x$npy)
n2 <- length(blocks)
if(n2 < x$n) blocks <- c(blocks, rep(blocks[n2], x$n - n2))
else if(n2 > x$n) blocks <- blocks[1:x$n]
yEmp <- c(aggregate(y, by=list(blocks), max)$x)
} else mod2 <- model
if(is.null(x$units)) ylb <- "Return Level"
else ylb <- paste("Return Level (", x$units, ")", sep="")
if(!tform) {
yrl <- rlevd(rperiods, loc=loc, scale=scale, shape=shape, threshold=u, type=mod2, npy=npy, rate=lam)
bds <- ci(x, return.period=rperiods)
if(is.null(args$ylim)) yl <- range(c(bds), finite=TRUE)
if(is.element(model, c("PP", "GEV", "Gumbel", "Weibull", "Frechet"))) xrl <- -1/(log(1 - 1/rperiods))
else xrl <- rperiods
# if(is.null(x$units)) ylb <- "Return Level"
# else ylb <- paste("Return Level (", x$units, ")", sep="")
xlb <- paste("Return Period (", x$period.basis, "s)", sep="")
if(is.null(args$main)) {
if(!is.null(args$ylim)) plot(xrl, yrl, type="l", log="x", xlab=xlb, ylab=ylb, main=m5, ...)
else plot(xrl, yrl, type="l", log="x", ylim=yl, xlab=xlb, ylab=ylb, main=m5, ...)
} else {
if(!is.null(args$ylim)) plot(xrl, yrl, type="l", log="x", xlab=xlb, ylab=ylb, ...)
else plot(xrl, yrl, type="l", log="x", ylim=yl, xlab=xlb, ylab=ylb, ...)
}
lines(xrl, bds[,1], col="gray", lty=2, lwd=2)
lines(xrl, bds[,3], col="gray", lty=2, lwd=2)
if(is.element(model, c("GEV", "Gumbel", "Weibull", "Frechet"))) {
points(-1/log(xp), sort(y))
if( type == "rl" ) out <- list( model = bds,
empirical = data.frame( transformed.period = -1/log(xp), sorted.level = sort( y ) ) )
} else if(is.element(model, c("GP", "Beta", "Pareto", "Exponential"))) {
n2 <- x$n
if(is.null(a)) xp2 <- ppoints(n2)
else xp2 <- ppoints(n2, a=a)
sdat <- sort(y)
points(-1/log(xp2)[sdat > u]/npy, sdat[sdat > u])
if( type == "rl" ) {
out <- list( model = bds,
empirical = data.frame( transformed.period = -1/log(xp2)[sdat > u]/npy,
sorted.level = sdat[sdat > u] ) )
}
} else if(model == "PP") {
if(is.null(a)) xp2 <- ppoints(length(yEmp))
else xp2 <- ppoints(length(yEmp), a=a)
points(-1/log(xp2), sort(yEmp))
if( type == "rl" ) out <- list( model = bds,
empirical = data.frame( transformed.period = -1/log(xp2),
sorted.level = sort( yEmp) ) )
}
} else {
np <- length(rperiods)
effrl <- matrix(NA, length(y), np)
for(i in 1:np) effrl[,i] <- erlevd(x = x, period = rperiods[i])
if(is.null(args$ylim)) {
yl <- range(c(c(y), c(effrl)), finite=TRUE)
yl[2] <- yl[2] + sign(yl[2]) * log2(abs(yl[2]))
if(is.null(args$main)) {
# if(!is.element(model, c("GP", "Beta", "Exponential", "Pareto")))
plot(y, type="l", xlab="index", ylab=ylb, main=m5, ylim=yl, ...)
# else plot(y[eid], type="l", xlab="index", ylab=ylb, main=m5, ylim=yl, ...)
} else {
# if(!is.element(model, c("GP", "Beta","Exponential","Pareto")))
plot(y, type="l", xlab="index", ylab=ylb, ...)
# else plot(y[eid], type="l", xlab="index", ylab=ylb, ylim=yl, ...)
}
if( type == "rl" ) out <- list( series = y, level = effrl )
} else {
if(is.null(args$main)) {
# if(!is.element(model, c("GP", "Beta","Exponential","Pareto")))
plot(y, type="l", xlab="index", ylab=ylb, main=m5, ...)
# else plot(y[eid], type="l", xlab="index", ylab=ylb, main=m5, ...)
} else {
# if(!is.element(model, c("GP", "Beta", "Exponential", "Pareto")))
plot(y, type="l", xlab="index", ylab=ylb, ...)
# else plot(y[eid], type="l", xlab="index", ylab=ylb, ...)
}
if( type == "rl" ) out <- list( series = y[ eid ], level = effrl )
} # end of if else 'ylim passed via '...' stmts.
for(i in 1:np) lines(effrl[,i], lty=i, col=i+1)
if(is.element(model, c("PP", "GP", "Beta", "Exponential", "Pareto"))) {
if(length(x$threshold) == 1) abline(h=x$threshold, col="darkorange", lwd=2)
else lines(x$threshold, col="darkorange", lwd=2)
}
if(!is.element(model, c("GEV", "Gumbel", "Weibull", "Frechet"))) {
legend("topleft", legend=c(paste(rperiods, "-", period, " level", sep=""), "threshold"), lty=c(1:np,1), col=c(2:(np+1),"darkorange"), bg="white")
} else legend("topleft", legend=paste(rperiods, "-", period, " level", sep=""), lty=1:np, col=2:(np + 1), bg="white")
} # end of if else '!tform' stmts.
} # end of if model is not POT and threshold varies with non-varying parameters stmts.
} # end of if 'make return level plots' stmts.
if(type=="trace") {
op <- par()
ntheta <- length(pars)
if(is.null(prange)) {
theta.hat <- pars
if(check.constant(x$par.models$scale)) {
phiU <- FALSE
if(x$par.models$log.scale) theta.hat["log.scale"] <- exp(theta.hat["log.scale"])
names(theta.hat)[names(theta.hat)=="log.scale"] <- "scale"
} else phiU <- x$par.models$log.scale
prange <- matrix(NA, 2, ntheta)
tmp <- summary(x, silent=TRUE)
tmp <- rbind(tmp$par, tmp$se.theta)
if(is.matrix(tmp)) for(j in 1:ntheta) prange[,j] <- c(tmp[1,j] - 2*tmp[2,j], tmp[1,j] + 2*tmp[2,j])
else {
tmp <- c(tmp)
for(j in 1:ntheta) prange[,j] <- c(tmp[j] - 2*log2(abs(tmp[j])), tmp[j] + 2*log2(abs(tmp[j])))
}
if(any(tmp.id <- names(tmp)=="scale")) {
if(!phiU) {
prange[1,tmp.id] <- max(prange[1,tmp.id], 1e-8)
if(prange[1,tmp.id] > prange[2, tmp.id]) prange[2, tmp.id] <- 2*prange[1, tmp.id]
}
}
} # end of if 'prange' is null or not stmts.
hold <- list()
par(mfrow=c(2,ntheta), oma=c(0,0,2,0))
for(i in 1:ntheta) {
if(!is.null(data)) hold[[i]] <- grlevdTracer(x=y, p=theta.hat, which.vary=i, p1.range=c(prange[,i]), threshold=u, threshold.fun=x$par.models$threshold,
location.fun=x$par.models$location, scale.fun=x$par.models$scale, shape.fun=x$par.models$shape, data=data, phi=phiU, blocks=x$blocks, # CJP2 - changed u/threshold args
type=model, npy=x$npy, na.action=x$na.action, par1.name=names(theta.hat)[i], plot=FALSE)
else hold[[i]] <- grlevdTracer(x=y, p=theta.hat, which.vary=i, p1.range=c(prange[,i]), threshold=u, threshold.fun=x$par.models$threshold, phi=phiU, blocks=x$blocks, # CJP2 - changed u/threshold arg calls
type=model, npy=x$npy, na.action=x$na.action, par1.name=names(theta.hat)[i], plot=FALSE)
if(i==1) plot(hold[[i]], type="likelihood", main="", xlab="")
else plot(hold[[i]], type="likelihood", main="", xlab="", ylab="")
abline(v=theta.hat[i], lty=2)
} # end of for 'i' loop.
for(i in 1:ntheta) {
if(i==1) plot(hold[[i]], type="gradient", main="")
else plot(hold[[i]], type="gradient", main="", ylab="")
} # end of second for 'i' loop.
out <- hold
} # end of if 'type == trace' stmts.
if(is.element(type, c("primary","trace"))) {
mtext(deparse(x$call), line=0.5, outer=TRUE)
par( mfrow = op$mfrow, oma = op$oma )
}
if( type != "primary" ) invisible( out )
else invisible()
} # end of 'plot.fevd.mle' function.
erlevd <- function(x, period=100) {
# CJP: this function now returns results per block rather than per observation when blocks is part of 'x'
type <- x$type
y <- datagrabber(x)
if(x$data.name[2] != "") {
data <- y[,-1]
y <- y[,1]
} else data <- NULL
pars <- findpars(x, use.blocks = TRUE) # CJP2: when x$blocks not NULL, findpars will now return pars with one row per block not per obs.
if(is.null(pars$location)) pars$location <- rep(0, x$n)
if(is.null(x$threshold)) u <- rep(0, x$n)
else if(length(x$threshold)==1) u <- rep(x$threshold, x$n)
else u <- x$threshold
if(!is.null(x$threshold)) lam <- mean(y > u)
else lam <- rep(1, x$n) # Eric -- changed 'len' to 'x$n'
if(!is.null(x$blocks)) { # CJP2
u <- x$blocks$threshold
if(is.null(pars$location)) pars$location <- rep(0, x$nBlocks) # this should never happen since when have blocks is only for PP model...
}
# CJP2: commented this out
# if(type=="PP") {
# scale <- pars$scale
# scale <- scale + pars$shape * (u - pars$location)
#} else scale <- pars$scale
theta <- cbind(pars$location, pars$scale, pars$shape, u) # CJP2: now 'pars$scale' not 'scale'
# if(is.element(type,c("PP","GP","Beta","Exponential","Pareto"))) theta <- theta[y > u,]
if(type=="PP") type2 <- "GEV" # CJP2
else type2 <- type
ifun <- function(theta, pd, npy, rate) rlevd(period=pd, loc=theta[1], scale=theta[2], shape=theta[3], threshold=theta[4], type=type2, npy=npy, rate=rate)
out <- apply(theta, 1, ifun, pd=period, npy=x$npy, rate=lam)
return(out)
} # end of 'erlevd' function.
findpars <- function(x, ...) {
UseMethod("findpars", x)
} # end of 'findpars' function.
findpars.fevd <- function(x, ...) {
meth <- tolower(x$method)
if( meth == "gmle" ) meth <- "mle"
# newcl <- paste("fevd.", meth, sep="")
# class(x) <- c( newcl, class( x ) )
# UseMethod( "findpars", x )
get( paste( "findpars.fevd.", meth, sep = "" ) )( x = x, ... )
} # end of 'findpars.fevd' function.
findpars.fevd.lmoments <- function(x, ...) {
return(x$results)
} # end of 'findpars.fevd.lmoments' function.
findpars.fevd.bayesian <- function(x, burn.in=499, FUN="mean", use.blocks=FALSE, ..., qcov = NULL) {
model <- x$type
tform <- !is.fixedfevd(x)
f <- match.fun(FUN)
p <- p2 <- x$results
np <- ncol(p) - 1
if(burn.in > 0) p <- p[-(1:burn.in),1:np]
pnames <- colnames(p)
if(is.element("log.scale", pnames)) {
id <- pnames == "log.scale"
p[,id] <- exp(p[,id])
pnames[id] <- "scale"
colnames(p) <- pnames
}
if(FUN=="mean") p <- colMeans(p)
else p <- apply(p, 2, f)
if(!tform) {
if(is.element(model, c("PP","GEV","Gumbel","Weibull","Frechet"))) {
loc <- p["location"]
nloc <- 1
} else {
loc <- NULL
nloc <- 0
}
scale <- p["scale"]
if(!is.element(model, c("Gumbel","Exponential"))) {
shape <- p["shape"]
} else shape <- NULL
ytrans <- NULL
} else {
# Eric -- 8/8/13 -- I don't think this was used.
# ytrans <- trans(x)
if(is.null(qcov)) { # Eric 8/8/13 -- Attempting to allow for 'qcov' argument.
if(is.null(x$blocks) || !use.blocks) { # CJP2
designs <- setup.design(x)
if(is.element(model, c("PP","GEV","Gumbel","Weibull","Frechet"))) {
X.loc <- designs$X.loc
nloc <- ncol(X.loc)
loc <- rowSums(matrix(p[1:nloc], x$n, nloc, byrow=TRUE) * X.loc)
} else {
loc <- NULL
nloc <- 0
}
X.sc <- designs$X.sc
nsc <- ncol(X.sc)
scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], x$n, nsc, byrow=TRUE) * X.sc)
if(x$par.models$log.scale && nsc > 1) scale <- exp(scale)
if(!is.element(model, c("Gumbel","Exponential"))) {
X.sh <- designs$X.sh
nsh <- ncol(X.sh)
shape <- rowSums(matrix(p[(nloc+nsc+1):np], x$n, nsh, byrow=TRUE) * X.sh)
} else shape <- NULL
} else {
# blocks$designs should already exist
if(is.element(model, c("PP","GEV","Gumbel","Weibull","Frechet"))) {
X.loc <- x$blocks$designs$X.loc
nloc <- ncol(X.loc)
loc <- rowSums(matrix(p[1:nloc], x$blocks$nBlocks, nloc, byrow=TRUE) * X.loc)
} else {
loc <- NULL
nloc <- 0
}
X.sc <- x$blocks$designs$X.sc
nsc <- ncol(X.sc)
scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], x$blocks$nBlocks, nsc, byrow=TRUE) * X.sc)
if(x$par.models$log.scale) scale <- exp(scale)
if(!is.element(model, c("Gumbel","Exponential"))) {
X.sh <- x$blocks$designs$X.sh
nsh <- ncol(X.sh)
shape <- rowSums(matrix(p[(nloc+nsc+1):np], x$blocks$nBlocks, nsh, byrow=TRUE) * X.sh)
} else shape <- NULL
}
} else {
nr <- nrow(qcov)
if(is.element(model, c("GP", "Beta", "Exponential", "Pareto"))) nloc <- 0
else if(is.element("location", pnames)) nloc <- 1
else nloc <- sum(substring(pnames, 1, 2) == "mu")
if(x$par.models$log.scale) {
if(!any(substring(pnames, 1, 3) == "phi")) nsc <- 1
else nsc <- sum(substring(pnames, 1, 3) == "phi")
} else {
if(is.element("scale", pnames)) nsc <- 1
else nsc <- sum(substring(pnames, 1, 3) == "sig")
}
if(is.element("shape", pnames)) nsh <- 1
else nsh <- sum(substring(pnames, 1, 2) == "xi")
if(!is.element(model, c("GP", "Beta", "Pareto"))) loc <- rowSums(matrix(rep(p[1:nloc], nr), nloc, nr, byrow = TRUE) * qcov[,1:nloc])
else loc <- NULL
scale <- rowSums(matrix(rep(p[(nloc + 1):(nloc + nsc)], nr), nsc, nr, byrow=TRUE) * qcov[,(nloc + 1):(nloc + nsc)])
if(x$par.models$log.scale && nsc > 1) scale <- exp(scale)
if(!is.element(model, c("Gumbel", "Exponential"))) shape <- rowSums(matrix(rep(p[(nloc + nsc + 1):(nloc + nsc + nsh)], nr), nsh, nr, byrow = TRUE) *
qcov[,(nloc + nsc + 1):(nloc + nsc + nsh)])
else shape <- NULL
} # end of if else 'qcov' is null stmts.
} # end of if '!tform' stmts.
return(list(location=loc, scale=scale, shape=shape))
} # end of 'findpars.fevd.bayesian' function.
findpars.fevd.mle <- function(x, use.blocks=FALSE, ..., qcov = NULL) { # Eric 8/8/13 -- added 'qcov' argument.
type <- tolower(x$type)
p <- x$results$par
np <- length(p)
if(is.null(qcov)) { # Eric 8/8/13
if(is.null(x$blocks) || !use.blocks) { # CJP2
n <- x$n
designs <- setup.design(x)
} else {
n <- x$blocks$nBlocks
designs <- x$blocks$designs
}
if(is.element(type, c("gev", "weibull", "gumbel", "frechet","pp"))) {
X.loc <- designs$X.loc
nloc <- ncol(X.loc)
loc <- rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE)*X.loc)
} else {
nloc <- 0
loc <- NULL
}
X.sc <- designs$X.sc
nsc <- ncol(X.sc)
scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc)
if(x$par.models$log.scale) scale <- exp(scale)
if(!is.element(type, c("gumbel","exponential"))) {
X.sh <- designs$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
}
} else {
pnames <- names(p)
model <- x$type
nr <- nrow(qcov)
if(is.element(model, c("GP", "Beta", "Exponential", "Pareto"))) nloc <- 0
else if(is.element("location", pnames)) nloc <- 1
else nloc <- sum(substring(pnames, 1, 2) == "mu")
if(x$par.models$log.scale) {
if(!any(substring(pnames, 1, 3) == "phi")) nsc <- 1
else nsc <- sum(substring(pnames, 1, 3) == "phi")
} else {
if(is.element("scale", pnames)) nsc <- 1
else nsc <- sum(substring(pnames, 1, 3) == "sig")
}
if(is.element("shape", pnames)) nsh <- 1
else nsh <- sum(substring(pnames, 1, 2) == "xi")
if(!is.element(model, c("GP", "Beta", "Pareto"))) loc <- rowSums(matrix(rep(p[1:nloc], nr), nloc, nr, byrow = TRUE) * qcov[,1:nloc])
else loc <- NULL
scale <- rowSums(matrix(rep(p[(nloc + 1):(nloc + nsc)], nr), nsc, nr, byrow=TRUE) * qcov[,(nloc + 1):(nloc + nsc)])
if(x$par.models$log.scale) scale <- exp(scale)
if(!is.element(model, c("Gumbel", "Exponential"))) shape <- rowSums(matrix(rep(p[(nloc + nsc + 1):(nloc + nsc + nsh)], nr), nsh, nr, byrow = TRUE) *
qcov[,(nloc + nsc + 1):(nloc + nsc + nsh)])
else shape <- NULL
} # end of if else 'qcov' is null stmts.
return(list(location=loc, scale=scale, shape=shape))
} # end of 'findpars.fevd.mle' function.
profliker <- function(object, type=c("return.level", "parameter"), xrange=NULL, return.period=100, which.par=1, nint=20, plot=TRUE,
gr=NULL, method="BFGS", lower=-Inf, upper=Inf, control=list(), ...) {
type <- tolower(type)
type <- match.arg(type)
designs <- setup.design(object)
tmp <- datagrabber(object)
if(object$data.name[2] == "") {
y <- tmp
data <- NULL
} else {
y <- tmp[,1]
data <- tmp[,-1]
}
thresh <- object$threshold
pars <- object$results$par
pnames <- names(pars)
args <- list(...)
if(missing(nint) && object$type=="PP") nint <- 20
if(object$type=="PP" && type=="return.level") {
object$type <- "GEV"
plot.msg <- "Return levels based on equivalent GEV df"
} else if(plot) plot.msg <- NULL
if(type=="return.level") {
if(!is.fixedfevd(object)) stop("profliker: Sorry, currently no support for profile likelihood with effective return levels.")
if(is.null(xrange)) {
rl.hat <- rl.fevd(object, period=return.period)
if(rl.hat != 0) xrange <- c(rl.hat - 2 * abs(log2(abs(rl.hat))), rl.hat + 2 * abs(log2(abs(rl.hat))))
else xrange <- range(y, finite=TRUE)
}
} else {
fv <- pars[which.par]
if(is.null(xrange)) {
tmp <- summary(object, silent=TRUE)
tmp <- rbind(tmp$par, tmp$se.theta)
if(is.matrix(tmp)) {
look <- c(tmp[,which.par])
xrange <- c(look[1] - 4 * look[2], look[1] + 4 * look[2])
th.names <- colnames(tmp)
} else {
look <- tmp[which.par]
xrange <- c(look - 4 * abs(log2(abs(look))), look + 4 * abs(log2(abs(look))))
th.names <- names(tmp)
}
if((th.names[which.par] == "scale") && !object$par.models$log.scale) {
xrange[1] <- max(xrange[1], 1e-10)
if(xrange[1] > xrange[2]) xrange[2] <- 2 * xrange[1]
}
} # end of if find 'xrange' stmts.
} # end of if else 'return.levels' stmts.
fixedvals <- seq(xrange[1], xrange[2],, nint)
res <- numeric(nint) + NA
if(type=="parameter") ipars <- pars[-which.par]
else {
if(is.element(object$type, c("GEV", "Gumbel", "Weibull", "Frechet"))) {
if(is.element("location", pnames)) ipars <- pars[pnames != "location"]
else {
id <- substring(pnames, 1, 2) == "mu"
ipars <- pars[!id]
}
} else {
if(is.element("scale", pnames)) ipars <- pars[pnames != "scale"]
else if(is.element("log.scale", pnames)) ipars <- pars[pnames != "log.scale"]
else {
istr <- substring(pnames, 1, 3)
id <- (istr == "sig") | (istr == "phi")
ipars <- pars[!id]
}
}
}
for(i in 1:nint) {
hold <- optim(par = ipars, fn = oevd.profpar, o=object, des=designs, x=y, data=data, u=thresh,
fixed.index=which.par, fixed.value=fixedvals[i], which.type=type, rperiod=return.period,
gr=gr, method=method, lower=lower, upper=upper, control=control)
ipars <- hold$par
res[ i ] <- hold$value
} # end of for 'i' loop.
res <- -res
if(plot) {
if(type == "parameter") xl <- names(pars)[which.par]
else xl <- paste(return.period, "-", object$period.basis, " return level", sep="")
xl <- paste(xl, plot.msg, sep="\n")
if(is.null(args$main)) {
msg <- deparse(object$call)
if(is.null(args$xlab)) plot(fixedvals, res, type="l", xlab=xl, ylab="Profile Log-Likelihood", main=msg, ...)
else plot(fixedvals, res, type="l", ylab="Profile Log-Likelihood", main=msg, ...)
} else {
if(is.null(args$xlab)) plot(fixedvals, res, type="l", xlab=xl, ylab="Profile Log-Likelihood", ...)
else plot(fixedvals, res, type="l", ylab="Profile Log-Likelihood", ...)
}
if(type=="parameter") abline(v=pars[which.par], lty=2, lwd=1.5)
else abline(v=rl.fevd(object, period=return.period), lty=2, lwd=1.5)
} # end of if 'plot' stmts.
invisible(res)
} # end of 'profliker' function.
oevd.profpar <- function(p, o, des, x, data=NULL, u=NULL, fixed.index, fixed.value, which.type=c("return.level","parameter"), rperiod=100) {
##
## 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.
## The only difference between this function and that of
## 'oevd' is that one of the parameters is held fixed in
## order to minimize a profile negative log-likelihood.
## Called by 'profliker' when a parameter's profile likelihood
## is of interest.
which.type <- match.arg(which.type)
type <- o$type
span <- o$span
npy <- o$npy
n <- length(x)
phi <- o$par.models$log.scale
np <- length(p)+1
if(which.type == "parameter") {
p.tmp <- numeric(np) + NA
p.tmp[-fixed.index] <- p
p.tmp[fixed.index] <- fixed.value
} else p.tmp <- p
if(which.type=="parameter") {
if(!is.element(type, c("GP","Beta","Pareto","Exponential"))) {
if(which.type=="parameter") {
X.loc <- des$X.loc
nloc <- ncol(X.loc)
loc <- rowSums(matrix(p.tmp[1:nloc], n, nloc, byrow=TRUE)*X.loc)
} else nloc <- 1
} else {
nloc <- 0
loc <- 0
}
} else {
nloc <- 0
loc <- 0
} # end of if 'which.type' parameter stmts.
if(which.type=="parameter" || is.element(type, c("GEV", "Gumbel", "Weibull", "Frechet", "PP"))) {
X.sc <- des$X.sc
nsc <- ncol(X.sc)
if(phi) scale <- exp(rowSums(matrix(p.tmp[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc))
else scale <- rowSums(matrix(p.tmp[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE)*X.sc)
} else nsc <- 0
if(!is.element(type, c("Gumbel","Exponential"))) {
X.sh <- des$X.sh
nsh <- ncol(X.sh)
shape <- rowSums(matrix(p.tmp[(nloc+nsc+1):(nloc+nsc+nsh)], n, nsh, byrow=TRUE)*X.sh)
} else {
nsh <- 0
shape <- 0
}
if(!is.element(type, c("GP","Beta","Pareto","Exponential")) && (which.type == "return.level")) {
if(type != "Gumbel") loc <- fixed.value + (scale/shape) * (1 - (-log(1 - 1/rperiod))^(-shape))
else loc <- fixed.value + scale * log(-log(1 - 1/rperiod))
} else if(which.type == "return.level") {
lam <- mean(x > u)
m <- rperiod * npy
if(type != "Exponential") {
zid <- shape == 0
if(any(zid)) shape[zid] <- 1e-10
scale <- ((fixed.value - u) * shape)/((m * lam)^shape - 1)
} else scale <- (fixed.value - u)/log(m * lam)
}
if(is.null(span)) res <- levd(x=x, threshold=u, location=loc, scale=scale, shape=shape, type=type, npy=npy, infval=1e10)
else res <- levd(x=x, threshold=u, location=loc, scale=scale, shape=shape, type=type, span=span, npy=npy, infval=1e10)
# if(class(res) == "try-error") res <- 1e10
if( o$method == "GMLE" ) res <- res + do.call( o$priorFun, c( list( x = p.tmp ), o$priorParams ) )
return(res)
} # end of 'oevd.profpar' function.
rl.fevd <- function(x, period=100) {
pmods <- x$par.models
if(is.fixedfevd(x)) {
p <- x$results$par
if(any(names(p) == "log.scale")) {
id <- names(p) == "log.scale"
p[id] <- exp(p[id])
names(p)[id] <- "scale"
}
if(!is.null(x$threshold)) lam <- mean(datagrabber(x)[,1] > x$threshold)
else lam <- NULL
if( "shape" %in% names( p ) ) {
return(rlevd(period=period, loc=p["location"], scale=p["scale"], shape=p["shape"], threshold=x$threshold, type=x$type, npy=x$npy, rate=lam))
} else return(rlevd(period=period, loc=p["location"], scale=p["scale"], shape=0, threshold=x$threshold, type=x$type, npy=x$npy, rate=lam))
} else return(erlevd(x, period=period))
} # end of 'rl.fevd' function.
# TO DO: Add possibility to make random PP draws using the rate from a Poisson
# instead of forcing the number of excesses to be n.
rextRemes <- function(x, n, ...) {
UseMethod("rextRemes", x)
} # end of 'rextRemes' function.
rextRemes.fevd <- function(x, n, ...) {
if( x$method == "GMLE" ) newcl <- "mle"
else newcl <- tolower( x$method )
# class(x) <- c( paste( "fevd.", newcl, sep="" ), class( x ) )
# UseMethod("rextRemes", x)
get( paste( "rextRemes.fevd.", newcl, sep = "" ) )( x = x, n = n, ... )
} # end of 'rextRemes.fevd' function.
rextRemes.fevd.lmoments <- function(x, n, ...) {
if( missing( n ) ) n <- x$n
p <- x$results
pnames <- names(p)
if(any(pnames == "location")) loc <- p["location"]
else loc <- 0
scale <- p["scale"]
if(any(pnames == "shape")) shape <- p["shape"]
else shape <- 0
type <- x$type
if(is.element(type, c("PP","GP","Exponential","Beta","Pareto"))) u <- x$threshold
else u <- 0
if(type=="PP") {
scale <- scale + shape * (u - loc)
type <- "GP"
}
if(is.element(type, c("Gumbel","Weibull","Frechet"))) type <- "GEV"
else if(is.element(type, c("Beta","Exponential","Pareto"))) type <- "GP"
out <- revd(n=n, loc=loc, scale=scale, shape=shape, threshold=u, type=type)
return(out)
} # end of 'rextRemes.fevd.lmoments' function.
rextRemes.fevd.bayesian <- function(x, n, ..., burn.in=499, FUN="mean", qcov = NULL) {
type <- x$type
f <- match.fun(FUN)
if(is.fixedfevd(x)) {
p <- x$results
np <- ncol(p) - 1
p <- p[,1:np]
pnames <- colnames(p)
if(burn.in > 0) p <- p[-(1:burn.in),]
if(is.element("log.scale",pnames)) {
id <- pnames == "log.scale"
p[,id] <- exp(p[,id])
pnames[id] <- "scale"
colnames(p) <- pnames
}
if(FUN == "mean") p <- colMeans(p)
else p <- apply(p, 2, f)
if(is.element(type, c("PP","GEV","Gumbel","Weibull","Frechet"))) {
loc <- p["location"]
nloc <- 1
} else loc <- 0
scale <- p["scale"]
if(!is.element(type, c("Gumbel","Exponential"))) shape <- p["shape"]
else shape <- 0
if(type == "PP") {
scale <- scale + shape * (x$threshold - loc)
type <- "GP"
}
return(revd(n=n, loc=loc, scale=scale, shape=shape, threshold=x$threshold, type=type))
} else {
# ytrans <- trans(x)
if(missing(n)) n <- 1
p <- findpars(x, burn.in=burn.in, FUN=FUN, qcov = qcov)
if(is.null(qcov)) K <- x$n
else K <- dim(qcov)[1]
if(!is.null(qcov)) u <- qcov[,"threshold"]
else if(is.null(x$threshold)) u <- rep(0, K)
else u <- x$threshold
if(is.null(p$location)) loc <- rep(0, K)
else loc <- p$location
if(is.null(p$shape)) {
if(type=="Gumbel") shape <- rep(0, K)
else shape <- rep(1e-10, K)
} else shape <- p$shape
if(type != "PP") theta <- cbind(u, loc, p$scale, shape)
else theta <- cbind(u, loc, p$scale + shape * (u - loc), shape)
if(is.null(qcov)) y <- c(datagrabber(x, cov.data=FALSE))
if(is.element(type,c("GEV","Gumbel","Weibull","Frechet"))) {
out <- matrix(NA, K, n)
if(is.null(qcov)) o <- order(y)
for(i in 1:n) {
z <- revd(K, type="GEV")
if(is.null(qcov)) z <- sort(z)[o]
out[,i] <- revtrans.evd(z=z, threshold=theta[,1], location=theta[,2], scale=theta[,3], shape=theta[,4], type=type)
} # end of for 'i' loop.
} else {
if(is.null(qcov)) {
eid <- y > u
m <- sum(eid)
theta <- theta[eid,]
# y <- y[eid]
# o <- order(y) # Eric 8/8/13 -- looks like this was not used.
} else {
m <- K
}
out <- matrix(NA, m, n)
for(i in 1:n) {
z <- revd(m, type="GP")
out[,i] <- revtrans.evd(z=z, threshold=theta[,1], location=theta[,2], scale=theta[,3], shape=theta[,4], type="GP")
} # end of for 'i' loop.
} # end of if else type of EVD stmts.
return(out)
} # end of if 'is.fixedfevd' stmts.
} # end of 'rextRemes.fevd.bayesian' function.
rextRemes.fevd.mle <- function(x, n, ..., qcov = NULL) {
type <- x$type
if(is.fixedfevd(x)) {
if(missing(n)) n <- x$n
p <- x$results$p
pnames <- names(p)
if(is.element("log.scale", pnames)) {
id <- pnames == "log.scale"
p[id] <- exp(p[id])
pnames[id] <- "scale"
names(p) <- pnames
}
scale <- p["scale"]
if(x$par.models$log.scale) scale <- exp(scale)
if(type=="GEV") {
loc <- p["location"]
shape <- p["shape"]
u <- 0
} else if(type=="GP") {
loc <- 0
u <- x$threshold
shape <- p["shape"]
} else if(type=="PP") {
loc <- p["location"]
u <- x$threshold
shape <- p["shape"]
scale <- scale + shape * (u - loc)
type <- "GP"
} else if(type=="Gumbel") {
loc <- p["location"]
shape <- 0
type <- "GEV"
u <- 0
} else if(is.element(type, c("Weibull","Frechet"))) {
loc <- p["location"]
shape <- p["shape"]
u <- 0
type <- "GEV"
} else if(type=="Exponential") {
loc <- 0
shape <- 0
u <- x$threshold
type <- "GP"
} else if(is.element(type, c("Beta","Pareto"))) {
loc <- 0
shape <- p["shape"]
u <- x$threshold
type <- "GP"
}
return(revd(n=n, loc=loc, scale=scale, shape=shape, threshold=u, type=type))
} else {
if(missing(n)) n <- 1
p <- findpars(x, qcov = qcov)
if(is.null(qcov)) K <- x$n
else K <- dim(qcov)[1]
if(!is.null(qcov)) u <- qcov[,"threshold"]
else if(is.null(x$threshold)) u <- rep(0, K)
else u <- x$threshold
if(is.null(p$location)) loc <- rep(0, K)
else loc <- p$location
if(is.null(p$shape)) {
if(type=="Gumbel") shape <- rep(0, K)
else shape <- rep(1e-10, K)
} else shape <- p$shape
scale <- p$scale
if(type=="PP") scale <- scale + shape * (u - loc)
theta <- cbind(u, loc, scale, shape)
if(is.null(qcov)) y <- c(datagrabber(x, cov.data=FALSE))
if(is.element(type,c("GEV","Gumbel","Weibull","Frechet"))) {
out <- matrix(NA, K, n)
for(i in 1:n) {
z <- revd(K, type="GEV")
out[,i] <- revtrans.evd(z=z, threshold=theta[,1], location=theta[,2], scale=theta[,3], shape=theta[,4], type=type)
} # end of for 'i' loop.
} else {
if(is.null(qcov)) {
eid <- y > u
m <- sum(eid)
theta <- theta[eid,]
} else {
m <- K
}
out <- matrix(NA, m, n)
if(type=="PP") type2 <- "GP" # Eric -- 8/8/13 -- should we change this to GEV?
else type2 <- type
for(i in 1:n) {
z <- revd(m, type="GP")
out[,i] <- revtrans.evd(z=z, threshold=theta[,1], location=theta[,2], scale=theta[,3], shape=theta[,4], type=type2)
} # end of for 'i' loop.
}
return(out)
}
} # end of 'rextRemes.fevd' function.
pextRemes <- function(x, q, lower.tail=TRUE, ...) {
if(missing(q)) stop("pextRemes: Must specify q argument.")
UseMethod("pextRemes", x)
} # end of 'pextRemes' function.
pextRemes.fevd <- function(x, q, lower.tail=TRUE, ..., qcov=NULL) {
if( x$method == "GMLE" ) newcl <- "mle"
else newcl <- tolower(x$method)
# class(x) <- c( paste("fevd.", newcl, sep=""), class( x ) )
# UseMethod("pextRemes", x)
get( paste( "pextRemes.fevd.", newcl, sep = "" ) )( x = x, q = q, lower.tail = lower.tail, ..., qcov = qcov )
} # end of 'pextRemes.fevd' function.
pextRemes.fevd.lmoments <- function(x, q, lower.tail=TRUE, ...) {
type <- x$type
p <- x$results
pnames <- names(p)
if(is.element("location", pnames)) loc <- p["location"]
else loc <- 0
scale <- p["scale"]
if(is.element("shape", pnames)) shape <- p["shape"]
else shape <- 0
if(!is.null(x$threshold)) u <- x$threshold
else u <- 0
if(type=="Gumbel") type <- "GEV"
else if(type=="Exponential") type <- "GP"
else if(type=="PP") {
scale <- scale + shape * (u - loc)
type <- "GP"
} else if(is.element(type, c("Weibull","Frechet"))) type <- "GEV"
else if(is.element(type, c("Beta","Pareto"))) type <- "GP"
return(pevd(q=q, loc=loc, scale=scale, shape=shape, threshold=u, lambda=x$rate, npy=x$npy, type=type, lower.tail=lower.tail))
} # end of 'pextRemes.fevd.lmoments' function.
pextRemes.fevd.bayesian <- function(x, q, lower.tail=TRUE, ..., qcov=NULL, burn.in=499, FUN="mean") {
type <- x$type
p <- x$results
np <- ncol(p) - 1
p <- p[,1:np]
pnames <- colnames(p)
if(burn.in > 0) p <- p[-(1:burn.in),]
f <- match.fun(FUN)
p <- apply(p, 2, f)
phi <- x$par.models$log.scale
if(is.fixedfevd(x)) {
if(is.element("location", pnames)) loc <- p["location"]
else loc <- 0
if(phi) scale <- exp(p["log.scale"])
else scale <- p["scale"]
if(is.element("shape", pnames)) shape <- p["shape"]
else shape <- 0
if(!is.null(x$threshold)) u <- x$threshold
else u <- 0
if(type=="Gumbel") type <- "GEV"
else if(type=="Exponential") type <- "GP"
else if(type=="PP") {
scale <- scale + shape * (u - loc)
type <- "GP"
} else if(is.element(type, c("Weibull","Frechet"))) type <- "GEV"
else if(is.element(type, c("Beta","Pareto"))) type <- "GP"
return(pevd(q=q, loc=loc, scale=scale, shape=shape, threshold=u, lambda=x$rate, npy=x$npy, type=type, lower.tail=lower.tail))
} else {
n <- length(q)
if(is.element("mu", substring(pnames, 1, 2))) nloc <- sum(substring(pnames, 1, 2) == "mu")
else if(is.element("location", pnames)) nloc <- 1
else nloc <- 0
if(is.element("phi", substring(pnames, 1, 3))) nsc <- sum(substring(pnames, 1, 3) == "phi")
else if(is.element("sig", substring(pnames, 1, 3))) nsc <- sum(substring(pnames, 1, 3) == "sig")
else nsc <- 1
if(is.element("shape", pnames)) nsh <- 1
else if(is.element("xi", substring(pnames, 1, 2))) nsh <- sum(substring(pnames, 1, 2) == "xi")
else nsh <- 0
if(is.null(qcov)) {
warning("pextRemes: qcov argument not supplied to non-stationary model. Using intercept terms only and/or if non-constant threshold, using first value.")
qcov <- matrix(0, n, np+1)
qcov[,1] <- 1
if(nloc > 0) qcov[,nloc+1] <- 1
qcov[,nloc+nsc+1] <- 1
if(!is.null(x$threshold)) qcov[,np+1] <- x$threshold[1]
else qcov[,np + 1] <- 0
colnames(qcov) <- c(pnames, "threshold")
}
if(!is.matrix(qcov)) stop("pextRemes: invalid type for qcov argument. Must be a matrix.")
if(dim(qcov)[2] == np+1) {
if(is.null(colnames(qcov))) colnames(qcov) <- c(pnames, "threshold")
else if(!is.element("threshold", colnames(qcov))) {
warning("pextRemes: qcov columns are named and dimension suggests threshold, but no threshold name. Assuming last column is threshold.")
colnames(qcov) <- c(pnames, "threshold")
}
} # end of if 'dim(qcov)[2] == np + 1' stmts.
u <- qcov[,"threshold"]
if(nloc==0) loc <- rep(0, n)
else loc <- rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE) * qcov[,1:nloc])
if(nsh == 0) {
if(type=="Gumbel") shape <- rep(0, n)
else shape <- rep(1e-10, n)
} else shape <- rowSums(matrix(p[(nloc+nsc+1):np], n, nsh, byrow=TRUE) * qcov[,(nloc+nsc+1):np])
scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE) * qcov[,(nloc+1):(nloc+nsc)])
if(phi) scale <- exp(scale)
if(type=="PP") {
scale <- scale + shape * (u - loc)
type <- "GP"
} else if(is.element(type, c("Gumbel","Weibull","Frechet"))) type <- "GEV"
else if(is.element(type, c("Exponential","Beta","Pareto"))) type <- "GP"
theta <- cbind(q, u, loc, scale, shape)
pfun <- function(th, type, rate, npy, lower.tail){
return(pevd(q=th[1], loc=th[3], scale=th[4], shape=th[5], threshold=th[2], lambda=rate, npy=npy, lower.tail=lower.tail))
} # end of internal 'pfun' function.
out <- apply(theta, 1, pfun, type=type, rate=x$rate, npy=x$npy, lower.tail=lower.tail)
return(out)
} # end of if else 'stationary model' stmts.
} # end of 'pextRemes.fevd.bayesian' function.
pextRemes.fevd.mle <- function(x, q, lower.tail=TRUE, ..., qcov=NULL) {
type <- x$type
p <- x$results$par
pnames <- names(p)
np <- length(p)
phi <- x$par.models$log.scale
if(is.fixedfevd(x)) {
if(phi) scale <- exp(p["log.scale"])
else scale <- p["scale"]
if(type=="GEV") {
loc <- p["location"]
shape <- p["shape"]
u <- 0
} else if(type=="GP") {
loc <- 0
u <- x$threshold
shape <- p["shape"]
} else if(type=="PP") {
loc <- p["location"]
u <- x$threshold
shape <- p["shape"]
scale <- scale + shape * (u - loc)
type <- "GP"
} else if(type=="Gumbel") {
loc <- p["location"]
shape <- 0
type <- "GEV"
u <- 0
} else if(is.element(type, c("Weibull","Frechet"))) {
loc <- p["location"]
shape <- p["shape"]
u <- 0
type <- "GEV"
} else if(type=="Exponential") {
loc <- 0
shape <- 0
u <- x$threshold
type <- "GP"
} else if(is.element(type, c("Beta","Pareto"))) {
loc <- 0
shape <- p["shape"]
u <- x$threshold
type <- "GP"
}
return(pevd(q=q, loc=loc, scale=scale, shape=shape, threshold=u, lambda=x$rate, npy=x$npy, type=type, lower.tail=lower.tail))
} else {
n <- length(q)
if(is.element("mu", substring(pnames, 1, 2))) nloc <- sum(substring(pnames, 1, 2) == "mu")
else if(is.element("location", pnames)) nloc <- 1
else nloc <- 0
if(is.element("phi", substring(pnames, 1, 3))) nsc <- sum(substring(pnames, 1, 3) == "phi")
else if(is.element("sig", substring(pnames, 1, 3))) nsc <- sum(substring(pnames, 1, 3) == "sig")
else nsc <- 1
if(is.element("shape", pnames)) nsh <- 1
else if(is.element("xi", substring(pnames, 1, 2))) nsh <- sum(substring(pnames, 1, 2) == "xi")
else nsh <- 0
if(is.null(qcov)) {
warning("pextRemes: qcov argument not supplied to non-stationary model. Using intercept terms only and/or if non-constant threshold, using first value.")
qcov <- matrix(0, n, np+1)
qcov[,1] <- 1
if(nloc > 0) qcov[,nloc+1] <- 1
qcov[,nloc+nsc+1] <- 1
if(!is.null(x$threshold)) qcov[,np+1] <- x$threshold[1]
else qcov[,np + 1] <- 0
colnames(qcov) <- c(pnames, "threshold")
}
if(!is.matrix(qcov)) stop("pextRemes: invalid type for qcov argument. Must be a matrix.")
if(dim(qcov)[2] == np+1) {
if(is.null(colnames(qcov))) colnames(qcov) <- c(pnames, "threshold")
else if(!is.element("threshold", colnames(qcov))) {
warning("pextRemes: qcov columns are named and dimension suggests threshold, but no threshold name. Assuming last column is threshold.")
colnames(qcov) <- c(pnames, "threshold")
}
}
if(!is.element("threshold", colnames(qcov))) {
if(is.null(x$threshold)) qcov <- cbind(qcov, 0)
else if(length(x$threshold)==1) qcov <- cbind(qcov, x$threshold)
else {
warning("pextRemes: non-constant threshold, but threshold values not supplied to qcov argument. Using first value only")
qcov <- cbind(qcov, x$threshold[1])
}
colnames(qcov) <- c(pnames, "threshold")
}
if((dim(qcov)[1] != n) || (!is.element(dim(qcov)[2], c(np,np+1)))) {
stop("pextRemes: invalid qcov argument. Must be a matrix with rows equal to the length of q and columns equal to the number of parameters, np, or np + 1.")
}
u <- qcov[,"threshold"]
if(nloc==0) loc <- rep(0, n)
else loc <- rowSums(matrix(p[1:nloc], n, nloc, byrow=TRUE) * qcov[,1:nloc])
if(nsh == 0) {
if(type=="Gumbel") shape <- rep(0, n)
else shape <- rep(1e-10, n)
} else shape <- rowSums(matrix(p[(nloc+nsc+1):np], n, nsh, byrow=TRUE) * qcov[,(nloc+nsc+1):np])
scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], n, nsc, byrow=TRUE) * qcov[,(nloc+1):(nloc+nsc)])
if(phi) scale <- exp(scale)
if(type=="PP") {
scale <- scale + shape * (u - loc)
type <- "GP"
} else if(is.element(type, c("Gumbel","Weibull","Frechet"))) type <- "GEV"
else if(is.element(type, c("Exponential","Beta","Pareto"))) type <- "GP"
theta <- cbind(q, u, loc, scale, shape)
pfun <- function(th, type, rate, npy, lower.tail) {
return(pevd(q=th[1], loc=th[3], scale=th[4], shape=th[5], threshold=th[2], lambda=rate, npy=npy, type=type, lower.tail=lower.tail))
} # end of internal 'pfun' function.
out <- apply(theta, 1, pfun, type=type, rate=x$rate, npy=x$npy, lower.tail=lower.tail)
return(out)
} # end of if else stationary model stmts.
} # end of 'pextRemes.fevd.mle' function
is.fixedfevd <- function(x) {
p <- x$par.models
return(all(c(check.constant(x$threshold), check.constant(p$threshold), check.constant(p$location), check.constant(p$scale), check.constant(p$shape))))
} # end of 'is.fixedfevd' function.
ci.fevd <- function(x, alpha=0.05, type=c("return.level", "parameter"), return.period=100, which.par, R=502, ...) {
if(x$method == "GMLE") newcl <- "fevd.mle"
else newcl <- paste("fevd.", tolower(x$method), sep="")
# class(x) <- c( newcl, class( x ) )
# UseMethod("ci",x)
get( paste( "ci.", newcl, sep = "" ) )( x = x, alpha = alpha, type = type, return.period = return.period, which.par = which.par, R = R, ... )
} # end of 'ci.fevd' function.
ci.fevd.bayesian <- function(x, alpha=0.05, type=c("return.level", "parameter"), return.period=100, which.par=1, FUN="mean", burn.in=499, tscale=FALSE, ...) {
type <- tolower(type)
type <- match.arg(type)
f <- match.fun(FUN)
const <- is.fixedfevd(x)
if(type=="return.level" && !const) return(ci.rl.ns.fevd.bayesian(x = x, alpha = alpha, return.period = return.period, FUN = FUN, burn.in = burn.in, ...))
pars <- x$results
np <- dim(pars)[2] - 1
pars <- pars[,1:np]
if(burn.in > nrow(pars)) stop("ci: burn.in > number of MCMC iterations (iter). Must be < iter.")
if(burn.in != 0) pars <- pars[-(1:burn.in),]
pnames <- colnames(pars)
iters <- dim(pars)[1]
if(type=="parameter" && missing(which.par)) which.par <- 1:np
if(is.element("log.scale",pnames)) {
id <- pnames=="log.scale"
pars[,id] <- exp(pars[,id])
pnames[id] <- "scale"
colnames(pars) <- pnames
}
conf.level <- paste(round((1 - alpha)*100, digits=2), "%", sep="")
if(type == "return.level") {
if(is.element("location",pnames)) loc <- pars[,"location"]
else loc <- rep(0, iters)
scale <- pars[,"scale"]
if(is.element("shape",pnames)) shape <- pars[,"shape"]
else shape <- rep(0, iters)
if(!is.null(x$threshold)) u <- x$threshold
else u <- 0
theta.sam <- cbind(loc, scale, shape)
rlfun <- function(th, pd, u, type, npy, rate) rlevd(period=pd, loc=th[1], scale=th[2], shape=th[3], threshold=u, type=type, npy=npy, rate=rate)
theta.sam <- apply(theta.sam, 1, rlfun, pd=return.period, u=u, type=x$type, npy=x$npy, rate=x$rate)
if(is.matrix(theta.sam)) {
if(FUN=="mean") theta.hat <- rowMeans(theta.sam)
else theta.hat <- apply(theta.sam, 1, f, ...)
bds <- apply(theta.sam, 1, quantile, probs=c(alpha/2,1-alpha/2))
est.names <- paste(names(theta.hat), "-", x$period.basis, " level", sep="")
bds.names <- rownames(bds)
out <- cbind(bds[1,], theta.hat, bds[2,])
rownames(out) <- est.names
colnames(out) <- c(bds.names[1], "Posterior Mean", bds.names[2])
} else {
if(FUN=="mean") theta.hat <- mean(theta.sam)
else theta.hat <- f(theta.sam)
bds <- quantile(theta.sam, probs=c(alpha/2,1-alpha/2))
out <- c(bds[1], theta.hat, bds[2])
names(out) <- c(paste(conf.level, " lower CI", sep=""),
paste("Posterior Mean ", return.period, "-", x$period.basis, " level", sep=""),
paste(conf.level, " upper CI", sep=""))
}
} else if(type=="parameter") {
theta.hat <- colMeans(pars)
est.names <- names(theta.hat)
if(tscale) {
if(!const && !is.element("scale",est.names) && !is.element("shape",est.names)) stop("ci: invalid argument configurations.")
if(!is.element(x$type, c("GP","Beta","Pareto"))) stop("ci: invalid argument configurations.")
pars[,"scale"] <- pars[,"scale"] - pars[,"shape"] * x$threshold
est.names[est.names == "scale"] <- "tscale"
colnames(pars) <- est.names
theta.hat["scale"] <- theta.hat["scale"] - theta.hat["shape"] * x$threshold
names(theta.hat) <- est.names
}
theta.hat <- theta.hat[which.par]
est.names <- est.names[which.par]
bds <- apply(pars, 2, quantile, probs=c(alpha/2,1-alpha/2))
if(is.matrix(bds)) bds <- bds[,which.par]
else bds <- bds[which.par]
if(is.matrix(bds)) {
bds.names <- rownames(bds)
out <- cbind(bds[1,], theta.hat, bds[2,])
rownames(out) <- est.names
colnames(out) <- c(bds.names[1], "Posterior Mean", bds.names[2])
} else {
out <- c(bds[1], theta.hat, bds[2])
names(out) <- c(paste(conf.level, " lower CI", sep=""),
paste("Posterior Mean ", est.names, " parameter", sep=""),
paste(conf.level, " upper CI", sep=""))
}
} else stop("ci: invalid type parameter.")
attr(out, "data.name") <- x$call
attr(out, "method") <- "Quantiles of MCMC Sample from Posterior Distribution"
attr(out, "conf.level") <- (1 - alpha) * 100
class(out) <- "ci"
return(out)
} # end of 'ci.fevd.bayesian' function.
ci.fevd.lmoments <- function(x, alpha=0.05, type=c("return.level", "parameter"), return.period=100, which.par, R=502, tscale=FALSE, return.samples=FALSE, ...) {
type <- tolower(type)
type <- match.arg(type)
p <- x$results
pnames <- names(p)
np <- length(p)
if(type=="parameter" && missing(which.par)) which.par <- 1:np
if(is.element(x$type, c("GEV","Gumbel","Frechet","Weibull"))) n <- x$n
else {
look <- datagrabber(x, cov.data=FALSE)
eid <- look > x$threshold
n <- sum(eid)
}
z <- rextRemes(x, n=R * n)
z <- matrix(z, n, R)
if(!is.element(x$type, c("GEV","Gumbel","Frechet","Weibull"))) {
z2 <- matrix(x$threshold, x$n, R)
z2[eid,] <- z
z <- z2
}
bfun <- function(y, ...) distill(fevd(y, ...))
if(x$type=="GEV") sam <- apply(z, 2, bfun, method="Lmoments")
else sam <- apply(z, 2, bfun, threshold=x$threshold, type=x$type, method="Lmoments", span=x$span)
if(tscale) {
if(!is.element(x$type, c("GP","Beta","Pareto"))) stop("ci: invalid argument configurations.")
sam["scale",] <- sam["scale",] - sam["shape",] * x$threshold
pnames[pnames == "scale"] <- "tscale"
rownames(sam) <- pnames
}
if(return.samples) {
if(is.matrix(sam)) out <- t(sam)
else out <- cbind(sam)
if(type=="parameter") return(out)
}
if(type=="return.level") {
if(is.element("location",pnames)) {
loc <- sam["location",]
loc0 <- p["location"]
} else loc <- loc0 <- 0
scale <- sam["scale",]
scale0 <- p["scale"]
if(is.element("shape",pnames)) {
shape <- sam["shape",]
shape0 <- p["shape"]
} else shape <- shape0 <- 0
if(!is.null(x$threshold)) u <- x$threshold
else u <- 0
theta <- rbind(loc, scale, shape)
rlfun <- function(th, ...) rlevd(period=return.period, loc=th[1], scale=th[2], shape=th[3], ...)
sam <- apply(theta, 2, rlfun, threshold=u, type=x$type, npy=x$npy, rate=x$rate)
if(is.matrix(sam)) rownames(sam) <- paste(return.period, "-", x$period.basis, sep="")
else sammy.name <- paste(return.period, "-", x$period.basis, " level", sep="")
if(return.samples) {
if(is.matrix(sam)) out <- cbind(out, t(sam))
else {
out <- cbind(out, c(sam))
onames <- colnames(out)
colnames(out) <- c(onames[-length(onames)], sammy.name)
}
return(out)
}
theta.hat <- rlevd(period=return.period, loc=loc0, scale=scale0, shape=shape0, threshold=u, type=x$type, npy=x$npy, rate=x$rate)
} else {
theta.hat <- p
if(tscale) {
if(!is.element(x$type, c("GP","Beta","Pareto"))) stop("ci: scale parameter is not a function of threshold for this df.")
sam["scale",] <- sam["scale",] - sam["shape",] * x$threshold
pnames[pnames == "scale"] <- "tscale"
rownames(sam) <- pnames
theta.hat["scale"] <- theta.hat["scale"] - theta.hat["shape"] * x$threshold
names(theta.hat) <- pnames
}
sam <- sam[which.par,]
if(!is.matrix(sam)) names(sam) <- pnames[which.par]
}
if(is.matrix(sam)) {
out <- apply(sam, 1, quantile, probs=c(alpha/2, 1 - alpha/2))
out.names <- rownames(out)
out <- rbind(out[1,], theta.hat, out[2,])
rownames(out) <- c(out.names[1], "Estimate", out.names[2])
colnames(out) <- rownames(sam)
out <- t(out)
} else {
out <- quantile(sam, probs=c(alpha/2, 1 - alpha/2))
onames <- names(out)
out <- c(out[1], theta.hat, out[2])
if(type=="parameter") names(out) <- c(onames[1], pnames[which.par], onames[2])
else names(out) <- c(onames[1], paste(return.period, "-", x$period.basis, " level", sep=""), onames[2])
}
attr(out, "method") <- "Parametric Bootstrap"
attr(out, "data.name") <- x$call
attr(out, "conf.level") <- (1 - alpha) * 100
attr(out, "R") <- R
class(out) <- "ci"
return(out)
} # end of 'ci.fevd.lmoments' function.
ci.fevd.mle <- function(x, alpha=0.05, type=c("return.level", "parameter"), return.period=100, which.par=1, R=502,
method=c("normal","boot","proflik"), xrange=NULL, nint=20, verbose=FALSE, tscale=FALSE,
return.samples=FALSE, ...) {
if(missing(method)) miss.meth <- TRUE
else miss.meth <- FALSE
method <- tolower(method)
method <- match.arg(method)
type <- tolower(type)
type <- match.arg(type)
theta.hat <- x$results$par
theta.names <- names(theta.hat)
np <- length(theta.hat)
if(type=="parameter" && missing(which.par)) which.par <- 1:np
if(any(theta.names=="log.scale")) {
id <- theta.names=="log.scale"
theta.hat[id] <- exp(theta.hat[id])
theta.names[id] <- "scale"
names(theta.hat) <- theta.names
}
const <- is.fixedfevd(x)
if(type=="return.level") par.name <- paste(return.period, "-", x$period.basis, " return level", sep="")
else if(type=="parameter") par.name <- theta.names[which.par]
if(type=="return.level" && !const) {
return(ci.rl.ns.fevd.mle(x = x, alpha = alpha, return.period = return.period, method = method, verbose = verbose, return.samples = return.samples, ...))
# stop("ci: Sorry, no current functionality for finding CI of effective return levels.")
}
if(type=="parameter") p <- theta.hat[which.par]
else {
if(is.element(x$type, c("PP","GP","Beta","Pareto","Exponential"))) lam <- mean(c(datagrabber(x)[,1]) > x$threshold)
else lam <- 1
if(is.element(x$type, c("PP","GEV","Gumbel","Weibull","Frechet"))) loc <- theta.hat["location"]
else loc <- 0
scale <- theta.hat["scale"]
if(!is.element(x$type, c("Gumbel","Exponential"))) shape <- theta.hat["shape"]
else shape <- 0
if(x$type == "PP") mod <- "GEV"
else mod <- x$type
p <- rlevd(period=return.period, loc=loc, scale=scale, shape=shape, threshold=x$threshold, type=mod, npy=x$npy, rate=lam)
}
if(verbose) {
cat("\n", "Preparing to calculate ", (1 - alpha)*100, "% CI for ",
ifelse(type=="return.level", paste(return.period, "-", x$period.basis, " return level", sep=""),
paste(par.name, " parameter", sep="")), "\n")
cat("\n", "Model is ", ifelse(const, " fixed", " non-stationary."), "\n")
if(method=="normal") cat("\n", "Using Normal Approximation Method.\n")
else if(method=="boot") cat("\n", "Using Bootstrap Method.\n")
else if(method=="proflik") cat("\n", "Using Profile Likelihood Method.\n")
}
if(method=="normal") {
method.name <- "Normal Approx."
z.alpha <- qnorm(alpha/2, lower.tail=FALSE)
cov.theta <- parcov.fevd(x)
if(is.null(cov.theta)) stop("ci: Sorry, unable to calculate the parameter covariance matrix. Maybe try a different method.")
var.theta <- diag(cov.theta)
if(any(var.theta < 0)) stop("ci: negative Std. Err. estimates obtained. Not trusting any of them.")
if(type=="parameter") {
se.theta <- sqrt(var.theta)
if(tscale) {
if(!const && !is.element("scale",theta.names) && !is.element("shape",theta.names) && !all(x$threshold == x$threshold[1])) {
stop("ci: invalid argument configurations.")
}
if(!is.element(x$type, c("GP","Beta","Pareto"))) stop("ci: invalid argument configurations.")
theta.hat["scale"] <- theta.hat["scale"] - theta.hat["shape"] * x$threshold
theta.names[theta.names == "scale"] <- "tscale"
if(!any(theta.names[which.par] == "tscale")) stop("ci: invalid argument configurations.")
names(theta.hat) <- theta.names
p <- theta.hat[which.par]
d <- rbind(1, -x$threshold)
names(se.theta) <- theta.names
se.theta["tscale"] <- sqrt(t(d) %*% cov.theta %*% d)
} else se.theta <- sqrt(var.theta)[which.par]
se.theta <- se.theta[which.par]
par.name <- theta.names[which.par]
} else if(type=="return.level") {
grads <- rlgrad.fevd(x, period=return.period)
grads <- t(grads)
if(is.element(x$type,c("GP","Beta","Pareto", "Exponential"))) {
if(x$type=="Exponential") cov.theta <- diag(c(lam * (1 - lam)/x$n, var.theta))
else cov.theta <- rbind(c(lam * (1 - lam)/x$n, 0, 0), cbind(0, cov.theta))
} else lam <- 1
var.theta <- t(grads) %*% cov.theta %*% grads # CJP2: took out sqrt as taking sqrt of off-diags
} else stop("ci: invalid type argument. Must be return.level or parameter.")
if(length(p) > 1) {
if(type=="return.level") se.theta <- sqrt(diag(var.theta))
out <- cbind(p - z.alpha * se.theta, p, p + z.alpha * se.theta)
rownames(out) <- par.name
conf.level <- paste(round((1 - alpha)*100, digits=2), "%", sep="")
colnames(out) <- c(paste(conf.level, " lower CI", sep=""), "Estimate", paste(conf.level, " upper CI", sep=""))
attr(out, "data.name") <- x$call
attr(out, "method") <- method.name
attr(out, "conf.level") <- (1 - alpha) * 100
class(out) <- "ci"
return(out)
} else out <- c(p - z.alpha * sqrt(var.theta[ which.par ]), p, p + z.alpha * sqrt(var.theta[ which.par ])) # CJP2 - using var.theta
} else if(method=="boot") {
method.name <- "Parametric Bootstrap"
if(verbose) cat("\n", "Simulating data from fitted model. Size = ", R, "\n")
if(const) {
if(is.null(x$blocks)) { # CJP
Z <- rextRemes(x, n=R * x$n)
Z <- matrix(Z, x$n, R)
} else {
Z <- rextRemes(x, n = round(R*x$npy*x$blocks$nBlocks))
Z <- matrix(Z, round(x$npy*x$blocks$nBlocks), R)
}
} else Z <- rextRemes(x, n=R)
if(verbose) cat("\n", "Simulated data found.\n")
y <- datagrabber(x, response=FALSE)
if(is.element(x$type, c("PP", "GP", "Exponential", "Beta", "Pareto"))) {
x2 <- datagrabber(x, cov.data=FALSE)
eid <- x2 > x$threshold
Z2 <- matrix(x$threshold, x$n, R)
Z2[eid,] <- Z[eid,]
Z <- Z2
lam <- mean(eid)
} else {
eid <- !logical(x$n)
lam <- 1
}
ipar <- list()
if(any(is.element(c("location","mu0"), theta.names))) {
if(is.element("location", theta.names)) ipar$location <- theta.hat["location"]
else {
id <- substring(theta.names,1,2) == "mu"
ipar$location <- theta.hat[id]
}
}
if(is.element("scale",theta.names)) ipar$scale <- theta.hat["scale"]
else {
if(!x$par.models$log.scale) id <- substring(theta.names,1,3) == "sig"
else id <- substring(theta.names, 1, 3) == "phi"
ipar$scale <- theta.hat[id]
}
if(!is.element(x$type, c("Gumbel","Exponential"))) {
if(is.element("shape",theta.names)) ipar$shape <- theta.hat["shape"]
else {
id <- substring(theta.names, 1, 2) == "xi"
ipar$shape <- theta.hat[id]
}
}
bfun <- function(z, x, y, p, ipar, eid, rate) {
pm <- x$par.models
if(is.null(y)) fit <- fevd(x=z, threshold=x$threshold, location.fun=pm$location, scale.fun=pm$scale, shape.fun=pm$shape,
use.phi=pm$log.scale, type=x$type, method=x$method, initial=ipar, span=x$span,
time.units=x$time.units, period.basis=x$period.basis, optim.args=x$optim.args,
priorFun=x$priorFun, priorParams=x$priorParams, verbose=FALSE)
else fit <- fevd(x=z, data=y, threshold=x$threshold, location.fun=pm$location, scale.fun=pm$scale, shape.fun=pm$shape,
use.phi=pm$log.scale, type=x$type, method=x$method, initial=ipar, span=x$span,
time.units=x$time.units, period.basis=x$period.basis, optim.args=x$optim.args,
priorFun=x$priorFun, priorParams=x$priorParams, verbose=FALSE)
fit$cov.data <- y
res <- distill(fit, cov=FALSE)
return(res)
} # end of internal 'bfun' function.
if(verbose) cat("\n", "Fitting model to simulated data sets (this may take a while!).")
if(type=="parameter") {
sam <- apply(Z, 2, bfun, x=x, y=y, ipar=ipar)
if(tscale) {
if(!const && !is.element("scale",theta.names) && !is.element("shape",theta.names)) stop("ci: invalid argument configurations.")
if(!is.element(x$type, c("GP","Beta","Pareto"))) stop("ci: invalid argument configurations.")
sam["scale",] <- sam["scale",] - sam["shape",] * x$threshold
theta.hat["scale"] <- theta.hat["scale"] - theta.hat["shape"] * x$threshold
theta.names[theta.names == "scale"] <- "tscale"
rownames(sam) <- theta.names
names(theta.hat) <- theta.names
} # end of if 'tscale' stmts.
sam <- sam[which.par,]
if(return.samples) return(t(sam))
} else if(type=="return.level") {
pars <- apply(Z, 2, bfun, x=x, y=y, ipar=ipar)[1:np,]
th.est <- numeric(3)
if(is.element(x$type, c("PP","GEV","Gumbel","Weibull","Frechet"))) {
loc <- pars["location",]
th.est[1] <- theta.hat["location"]
} else loc <- rep(0, R)
scale <- pars["scale",]
th.est[2] <- theta.hat["scale"]
if(!is.element(x$type, c("Gumbel","Exponential"))) {
shape <- pars["shape",]
th.est[3] <- theta.hat["shape"]
} else {
shape <- rep(1e-10, R)
th.est[3] <- 1e-8
}
if(return.samples) out <- t(pars)
th <- rbind(loc, scale, shape)
rlfun <- function(theta, p, u, type, npy, rate) rlevd(period=p, loc=theta[1], scale=theta[2], shape=theta[3], threshold=u, type=type, npy=npy, rate=rate)
if(x$type=="PP") mod <- "GEV"
else mod <- x$type
sam <- apply(th, 2, rlfun, p=return.period, u=x$threshold, type=mod, npy=x$npy, rate=lam)
if(is.matrix(sam)) rownames(sam) <- paste(rownames(sam), "-", x$period.basis, sep="")
else sammy.name <- paste(return.period, "-", x$period.basis, sep="")
if(return.samples) {
if(is.matrix(sam)) out <- cbind(pars, t(sam))
else {
onames <- colnames(out)
out <- cbind(out, sam)
colnames(out) <- c(onames, sammy.name)
}
return(out)
}
theta.hat <- rlevd(period=return.period, loc=th.est[1], scale=th.est[2], shape=th.est[3], threshold=x$threshold, type=x$type, npy=x$npy, rate=lam)
} else stop("ci: invalid type argument. Must be return.level or parameter.")
if(is.matrix(sam)) {
out <- apply(sam, 1, quantile, probs=c(alpha/2, 1 - alpha/2))
out.names <- rownames(out)
out <- rbind(out[1,], theta.hat, out[2,])
rownames(out) <- c(out.names[1], "Estimate", out.names[2])
colnames(out) <- rownames(sam)
out <- t(out)
attr(out, "data.name") <- x$call
attr(out, "method") <- method.name
attr(out, "conf.level") <- (1 - alpha) * 100
attr(out, "R") <- R
class(out) <- "ci"
return(out)
} else {
out <- quantile(sam, probs=c(alpha/2, 1 - alpha/2))
out <- c(out[1], mean(sam), out[2])
attr(out, "R") <- R
}
if(verbose) cat("\n", "Finished fitting model to simulated data.\n")
} else if( method == "proflik" ) {
if(x$type == "PP" && !is.null(x$blocks)) stop("ci: cannot do profile likelihood with blocks.") # CJP
if(tscale) stop("ci: invalid argument configurations.")
if( type == "parameter" && length( which.par ) > 1 ) stop("ci: can only do one parameter at a time with profile likelihood method.")
else if( type == "return.level" && length( return.period ) > 1 ) stop("ci: can only do one return level at a time with profile likelihood method.")
method.name <- "Profile Likelihood"
if(verbose) {
if(x$type != "PP") cat("\n", "Calculating profile likelihood. This may take a few moments.\n")
else cat("\n", "Calculating profile likelihood. This may take several moments.\n")
}
if(is.null(xrange)) {
hold2 <- c(ci(x, alpha=alpha, method="normal", type=type, return.period=return.period, which.par=which.par))[c(1,3)]
# if(!any(is.na(hold2))) xrange <- c(p - hold2[1], p + hold2[2])
if(!any(is.na(hold2))) xrange <- range(c(hold2, log2(hold2), 4 * hold2, hold2 - 4 * hold2, hold2 + 4 * hold2), finite = TRUE)
else if(!is.na(hold2[2])) xrange <- range(c(p - 2 * abs(log2(abs(p))), hold2[2], 4 * hold2[2], -4 * hold2[2], log2(p)), finite = TRUE)
else if(!is.na(hold2[1])) xrange <- range(c(p - 2 * abs(log2(abs(p))), hold2[1], 4 * hold2[1], -4 * hold2[1], log2(p)), finite = TRUE)
else if(all(is.na(hold2))) xrange <- c(p - 2 * abs(log2(abs(p))), p + 2 * abs(log2(abs(p))))
# else if(!is.na(hold2[2])) xrange <- c(p - 2 * abs(log2(abs(p))), p + hold2[2])
# else if(!is.na(hold2[1])) xrange <- c(p - hold2[1], p + 2 * abs(log2(abs(p))))
# else if(all(is.na(hold2))) xrange <- c(p - 2 * abs(log2(abs(p))), p + 2 * abs(log2(abs(p))))
if(verbose) cat("\n", "Using a range of ", xrange[1], " to ", xrange[2], "\n")
}
if(is.null(x$blocks)) {
if(!is.null(xrange)) hold <- profliker(x, type=type, xrange=xrange, return.period=return.period, which.par=which.par, nint=nint, plot=verbose, ...)
else hold <- profliker(x, type=type, return.period=return.period, which.par=which.par, nint=nint, plot=verbose, ...)
} else stop("Sorry: profile likelihood with blocks is not supported.")
# CJP: I haven't figured out if we can implement the profile likelihood with the blocks approach
ma <- -x$results$value
crit <- ma - 0.5 * qchisq(1 - alpha, 1)
if(verbose) {
cat("\n", "Profile likelihood has been calculated. Now, trying to find where it crosses the critical value = ", crit, "\n")
abline(h=crit, col="blue")
}
crit2 <- ma - 0.5 * qchisq((1 - alpha)+abs(log2(1-alpha))/2, 1)
id <- hold > crit2
z <- seq(xrange[1], xrange[2], length=length(hold))
z <- z[id]
parlik <- hold[id]
smth <- spline(z, parlik, n=200)
ind <- smth$y > crit
out <- range(smth$x[ind])
if(verbose) abline(v=out, lty=2, col="darkblue", lwd=2)
out <- c(out[1], p, out[2])
} else stop("ci: invalid method argument.")
conf.level <- paste(round((1 - alpha)*100, digits=2), "%", sep="")
names(out) <- c(paste(conf.level, " lower CI", sep=""), par.name, paste(conf.level, " upper CI", sep=""))
attr(out, "data.name") <- x$call
attr(out, "method") <- method.name
attr(out, "conf.level") <- (1 - alpha) * 100
class(out) <- "ci"
return(out)
} # end of 'ci.fevd.mle' function.
return.level.ns.fevd.bayesian <- function(x, return.period = 100, ..., burn.in = 499, FUN = "mean", do.ci = FALSE, verbose = FALSE, qcov = NULL, qcov.base = NULL) {
if(do.ci) res <- ci.rl.ns.fevd.bayesian(x = x, return.period = return.period, burn.in = burn.in, FUN = FUN, qcov = qcov, qcov.base = qcov.base, verbose = verbose, ...)
else res <- erlevd(x, period = return.period)
attr(res, "return.period") <- return.period
attr(res, "data.name") <- x$data.name
attr(res, "fit.call") <- x$call
attr(res, "call") <- match.call()
attr(res, "fit.type") <- x$type
attr(res, "data.assumption") <- "non-stationary"
attr(res, "period") <- x$period.basis
attr(res, "units") <- x$units
attr(res, "class") <- "return.level"
if(!do.ci) {
attr(res, "conf.level") <- NULL
class(res) <- "return.level"
} else class(res) <- c("return.level", "ci")
return(res)
} # end of 'return.level.ns.fevd.bayesian'
ci.rl.ns.fevd.bayesian <- function(x, alpha = 0.05, return.period = 100, FUN = "mean", burn.in = 499, ..., qcov = NULL, qcov.base = NULL, verbose = FALSE) {
if(verbose) begin.tiid <- Sys.time()
if(length(return.period) > 1) stop("ci.rl.ns.fevd.bayesian: return.period must have length 1.")
conf.level <- paste(round((1 - alpha) * 100, digits = 2), "%", sep = "")
np <- dim(x$results)[2] - 1
if(burn.in > 0) p <- x$results[-(1:burn.in), 1:np]
else p <- x$results[, 1:np]
pnames <- colnames(p)
ni <- dim(p)[1]
if(is.null(qcov) && !is.null(qcov.base)) {
warning("ci.rl.ns.fevd.bayesian: qcov must be supplied if qcov.base is supplied. Continuing as if qcov.base were qcov, and qcov.base were NULL.")
qcov <- qcov.base
qcov.base <- NULL
}
rlfun <- function(th, pd, type, npy, rate) {
return(rlevd(period = pd, loc = th[1], scale = th[2], shape = th[3], threshold = th[4], type = type, npy = npy, rate = rate))
} # end of internal 'rlfun' function.
if(!is.null(qcov)) {
if(verbose) cat("\n", "Calculating ", return.period, "-year effective return levels based on qcov.\n")
nq <- dim(qcov)[1]
## Setting up parameter matrices so that they have the same number of rows as the MCMC (less burn.in)
## and columns equal to the number of rows in the 'qcov' matrix.
##
# Set up location matrix
if(is.element("location", pnames)) {
loc <- matrix(p[, "location"], ni, nq)
nloc <- 1
} else if(is.element("mu", substring(pnames, 1, 2))) {
nloc <- sum(substring(pnames, 1, 2) == "mu")
loc <- matrix(NA, ni, nq)
for(i in 1:nq) loc[,i] <- rowSums(p[, 1:nloc, drop = FALSE] * matrix(qcov[i, 1:nloc], ni, nloc, byrow = TRUE))
} else {
loc <- matrix(0, ni, nq)
nloc <- 0
} # end of setting up location matrix stmts.
# Set up scale matrix.
if(is.element("log.scale", pnames)) {
nsc <- 1
scale <- matrix(exp(p[, "log.scale"]), ni, nq)
} else if(is.element("scale", pnames)) {
nsc <- 1
scale <- matrix(p[, "scale"], ni, nq)
} else if(is.element("sig", substring(pnames, 1, 3))) {
nsc <- sum(substring(pnames, 1, 3) == "sig")
scale <- matrix(NA, ni, nq)
for(i in 1:nq) scale[,i] <- rowSums(p[, (nloc + 1):(nloc + nsc), drop = FALSE] * matrix(qcov[i, (nloc + 1):(nloc + nsc)], ni, nsc, byrow = TRUE))
} else if(is.element("phi", substring(pnames, 1, 3))) {
nsc <- sum(substring(pnames, 1, 3) == "phi")
scale <- matrix(NA, ni, nq)
for(i in 1:nq) scale[,i] <- rowSums(p[, (nloc + 1):(nloc + nsc), drop = FALSE] * matrix(qcov[i, (nloc + 1):(nloc + nsc)], ni, nsc, byrow = TRUE))
scale <- exp(scale)
}
# Set up shape matrix
if(is.element("shape", pnames)) {
nsh <- 1
shape <- matrix(p[, "shape"], ni, nq)
} else if(is.element("xi", substring(pnames, 1, 2))) {
nsh <- sum(substring(pnames, 1, 2) == "xi")
shape <- matrix(NA, ni, nq)
for(i in 1:nq) rowSums(p[, (nloc + nsc + 1):(nloc + nsc + nsh), drop = FALSE] * matrix(qcov[i, (nloc + nsc + 1):(nloc + nsc + nsh)], ni, nsh, byrow = TRUE))
} else {
nsh <- 0
shape <- matrix(0, ni, nq)
}
u <- matrix(qcov[, "threshold"], ni, nq, byrow = TRUE)
rlmat <- matrix(NA, ni, nq)
for(i in 1:nq) rlmat[,i] <- apply(cbind(loc[,i], scale[,i], shape[,i], u[,i]), 1, rlfun, pd = return.period, type = x$type, npy = x$npy, rate = x$rate)
if(!is.null(qcov.base)) {
if(dim(qcov.base)[1] != nq) stop("ci.rl.ns.fevd.bayesian: qcov.base must have same number of rows as qcov.")
# Set up location matrix
if(is.element("location", pnames)) {
loc <- matrix(p[, "location"], ni, nq)
nloc <- 1
} else if(is.element("mu", substring(pnames, 1, 2))) {
nloc <- sum(substring(pnames, 1, 2) == "mu")
loc <- matrix(NA, ni, nq)
for(i in 1:nq) loc[,i] <- rowSums(p[, 1:nloc, drop = FALSE] * matrix(qcov.base[i, 1:nloc], ni, nloc, byrow = TRUE))
} else {
loc <- matrix(0, ni, nq)
nloc <- 0
} # end of setting up location matrix stmts.
# Set up scale matrix.
if(is.element("log.scale", pnames)) {
nsc <- 1
scale <- matrix(exp(p[, "log.scale"]), ni, nq)
} else if(is.element("scale", pnames)) {
nsc <- 1
scale <- matrix(p[, "scale"], ni, nq)
} else if(is.element("sig", substring(pnames, 1, 3))) {
nsc <- sum(substring(pnames, 1, 3) == "sig")
scale <- matrix(NA, ni, nq)
for(i in 1:nq) scale[,i] <- rowSums(p[, (nloc + 1):(nloc + nsc), drop = FALSE] * matrix(qcov.base[i, (nloc + 1):(nloc + nsc)], ni, nsc, byrow = TRUE))
} else if(is.element("phi", substring(pnames, 1, 3))) {
nsc <- sum(substring(pnames, 1, 3) == "phi")
scale <- matrix(NA, ni, nq)
for(i in 1:nq) scale[,i] <- rowSums(p[, (nloc + 1):(nloc + nsc), drop = FALSE] * matrix(qcov.base[i, (nloc + 1):(nloc + nsc)], ni, nsc, byrow = TRUE))
scale <- exp(scale)
}
# Set up shape matrix
if(is.element("shape", pnames)) {
nsh <- 1
shape <- matrix(p[, "shape"], ni, nq)
} else if(is.element("xi", substring(pnames, 1, 2))) {
nsh <- sum(substring(pnames, 1, 2) == "xi")
shape <- matrix(NA, ni, nq)
for(i in 1:nq) rowSums(p[, (nloc + nsc + 1):(nloc + nsc + nsh), drop = FALSE] * matrix(qcov.base[i, (nloc + nsc + 1):(nloc + nsc + nsh)], ni, nsh, byrow = TRUE))
} else {
nsh <- 0
shape <- matrix(0, ni, nq)
}
u <- matrix(qcov[, "threshold"], ni, nq, byrow = TRUE)
rlmat2 <- matrix(NA, ni, nq)
for(i in 1:nq) rlmat2[,i] <- apply(cbind(loc[,i], scale[,i], shape[,i], u[,i]), 1, rlfun, pd = return.period, type = x$type, npy = x$npy, rate = x$rate)
res <- rlmat - rlmat2
} else {
res <- rlmat
} # end of if 'qcov.base' stmts.
} else {
if(verbose) cat("\n", "Calculating ", return.period, "-year effective return levels based on data covariates.\n")
# Cannot use 'findpars' because we want parameters for every value of covariate at every iteration of MCMC chain.
designs <- setup.design(x)
X.loc <- designs$X.loc
if(!is.null(X.loc)) nloc <- ncol(X.loc)
else nloc <- 0
X.sc <- designs$X.sc
nsc <- ncol(X.sc)
X.sh <- designs$X.sh
if(!is.null(X.sh)) nsh <- ncol(X.sh)
else nsh <- 0
# loc <- scale <- shape <- res <- matrix(NA, x$n, ni)
if(is.null(X.loc)) loc <- matrix(0, x$n, ni)
if(is.null(X.sh)) shape <- matrix(0, x$n, ni)
if(is.null(x$threshold)) u <- matrix(0, x$n, ni)
else u <- matrix(x$threshold, x$n, ni)
parfinder <- function(z, y) {
return(rowSums(t(z * t(y)), na.rm = TRUE))
} # end of internal 'parfinder' function.
if(verbose) cat("\n", "Finding effective parameters for each iteration of MCMC sample.\n")
if(is.null(X.loc)) loc <- matrix(0, ni, x$n)
else loc <- apply(p[, 1:nloc, drop = FALSE], 1, parfinder, y = X.loc)
scale <- apply(p[, (nloc+1):(nloc+nsc), drop = FALSE], 1, parfinder, y = X.sc)
if(is.element("log.scale", pnames) || is.element("phi", substring(pnames, 1, 3))) scale <- exp(scale)
if(is.null(X.sh)) shape <- matrix(0, ni, x$n)
else shape <- apply(p[, (nloc+nsc+1):np, drop = FALSE], 1, parfinder, y = X.sh)
do.rl <- function(id, z, u, pd, type, npy, rate, verbose) {
if(verbose && id <= 5) cat(id, " ")
else if(verbose && (id %% 100 == 0)) cat(id, " ")
theta <- cbind(z$loc[,id], z$scale[,id], z$shape[,id], z$threshold[,id])
id <- theta[,3] == 0
res <- numeric(dim(theta)[1]) + NA
if(is.element(type, c("GEV","PP", "Gumbel"))) {
p <- 1 - 1/pd
if(any(id)) res[id] <- theta[id,1] - theta[id,2] * log(-log(p))
if(any(!id)) res[!id] <- theta[!id,1] + theta[!id,2] * ((-log(p))^(-theta[!id,3]) - 1)/theta[!id,3]
} else if(is.element(type, c("GP", "Exponential"))) {
m <- pd * npy * rate
if(any(id)) res[id] <- theta[id,4] + theta[id,2] * log(m)
if(any(!id)) res[!id] <- theta[!id,4] + (theta[!id,2]/theta[!id,3]) * (m^(theta[!id,3]) - 1)
}
return(res)
} # end of internal 'do.rl' function.
hold <- list(loc = loc, scale = scale, shape = shape, threshold = u)
if(verbose) cat("\n", "Calculating the return levels for ", ni, " samples after burn in period.\n")
res <- t(apply(matrix(1:ni, ncol = 1), 1, do.rl, z = hold, pd = 100, type = x$type, npy = x$npy, rate = x$rate, verbose = verbose))
} # end of if else 'qcov' stmts.
out <- cbind(c(apply(res, 2, quantile, probs = alpha / 2)), c(apply(res, 2, FUN, ...)), c(apply(res, 2, quantile, probs = 1 - alpha / 2)))
colnames(out) <- c(paste(conf.level, " lower CI", sep = ""), "Estimate", paste(conf.level, " upper CI", sep = ""))
attr(out, "conf.level") <- (1 - alpha) * 100
class(out) <- "ci"
attr(out, "return.period") <- return.period
attr(out, "data.name") <- x$data.name
attr(out, "fit.call") <- x$call
attr(out, "call") <- match.call()
attr(out, "fit.type") <- x$type
attr(out, "data.assumption") <- "non-stationary"
attr(out, "period") <- x$period.basis
attr(out, "units") <- x$units
attr(out, "class") <- "ci"
if(verbose) {
cat("\n", "Return levels and CIs estimated.\n")
print(Sys.time() - begin.tiid)
}
return(out)
} # end of 'ci.rl.ns.fevd.bayesian' function.
ci.rl.ns.fevd.mle <- function(x, alpha = 0.05, return.period = 100, method = c("normal"), verbose = FALSE, qcov = NULL, qcov.base = NULL, ...) {
method <- tolower(method)
method <- match.arg(method)
# if(method != "normal") stop("ci.rl.ns.fevd.mle: Currently only normal approximation method available for non-stationary MLE models.")
par.name <- paste(return.period, "-", x$period.basis, " return level", sep = "")
if (verbose) {
cat("\n", "Preparing to calculate ", (1 - alpha) *
100, "% CI for ", paste(return.period, "-", x$period.basis,
" return level", sep = ""), "\n")
cat("\n", "Model is non-stationary.\n")
# if(method == "normal")
cat("\n", "Using Normal Approximation Method.\n")
# else if(method == "boot") cat("\n", "Using Parametric Boot strap method with ", R, " replicate samples.\n")
}
if(method == "normal") method.name <- "Normal Approx."
# else if(method == "boot") method.name <- "Parametric Bootstrap"
if(method == "normal") {
res <- return.level.ns.fevd.mle(x = x, return.period = return.period, ..., do.ci = FALSE, verbose = verbose, qcov = qcov, qcov.base = qcov.base)
z.alpha <- qnorm(alpha/2, lower.tail = FALSE)
cov.theta <- parcov.fevd(x)
if (is.null(cov.theta))
stop("ci: Sorry, unable to calculate the parameter covariance matrix. Maybe try a different method.")
var.theta <- diag(cov.theta)
if (any(var.theta < 0))
stop("ci: negative Std. Err. estimates obtained. Not trusting any of them.")
grads <- t(rlgrad.fevd(x, period = return.period, qcov = qcov,
qcov.base = qcov.base))
se.theta <- sqrt(diag(t(grads) %*% cov.theta %*% grads))
out <- cbind(c(res) - z.alpha * se.theta, c(res), c(res) + z.alpha * se.theta, se.theta)
if (length(return.period) > 1)
rownames(out) <- par.name
else rownames(out) <- NULL
conf.level <- paste(round((1 - alpha) * 100, digits = 2),
"%", sep = "")
colnames(out) <- c(paste(conf.level, " lower CI", sep = ""),
"Estimate", paste(conf.level, " upper CI", sep = ""),
"Standard Error")
} else if(method == "boot") {
stop("ci.rl.ns.fevd.mle: Sorry, this functionality has not yet been added.")
}
attr(out, "data.name") <- x$call
attr(out, "method") <- method.name
attr(out, "conf.level") <- (1 - alpha) * 100
class(out) <- "ci"
return(out)
} # end of 'ci.rl.ns.fevd.mle' function.
rlgrad.fevd <- function(x, period=100, qcov=NULL, qcov.base=NULL) { # CJP2; many changes here for nonstationary models
# qcov.base is for when one wants to work with the difference in return levels for different covariate sets in a nonstationary model
type <- tolower(x$type)
if(!is.element(x$method, c("MLE","GMLE"))) stop("rlgrad.fevd: Estimation method must be MLE/GMLE.")
p <- x$results$par
if(is.element("log.scale",names(p))) {
id <- names(p) == "log.scale"
p[id] <- exp(p[id])
names(p)[id] <- "scale"
}
if(is.element("shape",names(p))) {
if(p["shape"] == 0) {
if(is.element(type, c("gev","pp","gumbel"))) type <- "gumbel"
else if(is.element(type, c("gp","exponential"))) type <- "exponential"
else stop("rlgrad.fevd: invalid type for the shape parameter.")
}
}
if(!is.fixedfevd(x)) { # CJP2: this whole block of code
if(is.null(qcov)) stop("rlgrad.fevd: qcov required for nonstationary models.")
if(!is.matrix(qcov)) qcov <- matrix(qcov, nrow = 1)
if(!is.null(qcov.base)) {
if(!is.matrix(qcov.base)) qcov.base <- matrix(qcov.base, nrow = 1)
if(nrow(qcov) != nrow(qcov.base) ||
ncol(qcov) != ncol(qcov.base))
stop("rlgrad.fevd: qcov and qcov.base must have the same number of covariates and values.")
}
if(length(period) > 1 && nrow(qcov) > 1)
stop("rlgrad.fevd: Cannot compute gradient for multiple return periods and multiple covariate values simultaneously.")
loc.id <- 1:x$results$num.pars$location
sc.id <- (1 + x$results$num.pars$location) : (x$results$num.pars$location + x$results$num.pars$scale)
sh.id <- (1 + x$results$num.pars$location + x$results$num.pars$scale) : (x$results$num.pars$location + x$results$num.pars$scale + x$results$num.pars$shape)
# loc <- qcov[ , loc.id, drop=FALSE] %*% p[loc.id] # not needed
sc <- c(qcov[ , sc.id, drop=FALSE] %*% p[sc.id])
sh <- c(qcov[ , sh.id, drop=FALSE] %*% p[sh.id])
if(!is.null(qcov.base)) {
# loc.base <- qcov.base[ , loc.id, drop=FALSE] %*% p[loc.id] # not needed
sc.base <- c(qcov.base[ , sc.id, drop=FALSE] %*% p[sc.id])
sh.base <- c(qcov.base[ , sh.id, drop=FALSE] %*% p[sh.id])
}
# CJP2 - this next block of code deals with the only case that params are on log scale when considering return levels, namely when use.phi=TRUE for nonstationarity in the scale
pnames <- names(p)
if(is.element("phi", substring(pnames, 1, 3))) {
sc <- exp(sc)
phi.gradTerm <- sc
if(!is.null(qcov.base)) {
sc.base <- exp(sc.base)
phi.gradTerm.base <- sc.base
}
} else{
phi.gradTerm <- 1
if(!is.null(qcov.base)) phi.gradTerm.base <- 1
}
if(is.element(type, c("pp","gev","weibull","frechet"))) {
yp <- -log(1 - 1/period)
if(length(period) == 1) {
res <- cbind(qcov[ , loc.id, drop = FALSE], -phi.gradTerm/sh * (1-yp^(-sh)) * qcov[ , sc.id, drop = FALSE],
(sc * sh^(-2) * (1 - yp^(-sh)) - sc/sh * yp^(-sh) * log(yp)) * qcov[ , sh.id, drop = FALSE])
if(!is.null(qcov.base))
res <- res - cbind(qcov.base[ , loc.id, drop = FALSE], -phi.gradTerm.base/sh.base * (1-yp^(-sh.base)) * qcov.base[ , sc.id, drop = FALSE], (sc.base * sh.base^(-2) * (1 - yp^(-sh.base)) - sc.base/sh.base * yp^(-sh.base) * log(yp)) * qcov.base[ , sh.id, drop = FALSE])
} else {
res <- cbind(outer(rep(1, length(period)), qcov[1, loc.id]), outer(-phi.gradTerm/sh * (1-yp^(-sh)), qcov[1, sc.id]), outer(sc * sh^(-2) * (1 - yp^(-sh)) - sc/sh * yp^(-sh) * log(yp), qcov[1, sh.id]))
if(!is.null(qcov.base))
res <- res - cbind(outer(rep(1, length(period)), qcov.base[1, loc.id]), outer(-phi.gradTerm.base/sh.base * (1-yp^(-sh.base)), qcov.base[1, sc.id]), outer(sc.base * sh.base^(-2) * (1 - yp^(-sh.base)) - sc.base/sh.base * yp^(-sh.base) * log(yp), qcov.base[1, sh.id]))
}
} else if(type=="gumbel") {
yp <- -log(1 - 1/period)
if(length(period) == 1) {
res <- cbind(qcov[ , loc.id, drop = FALSE], -log(yp)*phi.gradTerm * qcov[ , sc.id, drop = FALSE])
if(!is.null(qcov.base))
res <- res - cbind(qcov.base[ , loc.id, drop = FALSE], -log(yp)*phi.gradTerm.base * qcov.base[ , sc.id, drop = FALSE])
} else {
res <- cbind(outer(rep(1, length(period)), qcov[1, loc.id]), outer(-log(yp)*phi.gradTerm, qcov[1, loc.id]))
if(!is.null(qcov.base))
res <- res - cbind(outer(rep(1, length(period)), qcov.base[1, loc.id]), outer(-log(yp)*phi.gradTerm.base, qcov.base[1, loc.id]))
}
} else stop("rlgrad.fevd: not implemented for nonstationary models for GP, beta, pareto, exponential models.")
} else { # fixed.fevd
if(is.element(type, c("pp","gev","weibull","frechet"))) {
yp <- -log(1 - 1/period)
res <- cbind(1,
(-1/p["shape"]) * (1 - yp^(-p["shape"])),
p["scale"] * (p["shape"])^(-2) * (1 - yp^(-p["shape"])) - (p["scale"]/p["shape"]) * yp^(-p["shape"]) * log(yp))
} else if(type=="gumbel") {
yp <- -log(1 - 1/period)
res <- cbind(1, -log(yp))
} else if(is.element(type, c("gp","beta","pareto"))) {
lam <- mean(c(datagrabber(x)[,1]) > x$threshold)
m <- period * x$npy
mlam <- m * lam
res <- cbind(p["scale"] * m^(-p["shape"]) * lam^(-p["shape"] - 1),
(p["shape"])^(-1) * ((mlam)^(p["shape"]) - 1),
-p["scale"] * (p["shape"])^(-2)*((mlam)^(p["shape"]) - 1) + (p["scale"]/p["shape"]) * (mlam)^(p["shape"]) * log(mlam))
} else if(type=="exponential") {
lam <- mean(c(datagrabber(x)[,1]) > x$threshold)
m <- period * x$npy
mlam <- m * lam
res <- cbind(p["scale"]/mlam, log(mlam))
} else stop("rlgrad.fevd: Hmmm. You should not be getting this error message. Something is horribly wrong.")
} # end of if else use actual gradients or finite differences stmts.
return(res)
} # end of 'rlgrad.fevd' function.
return.level <- function(x, return.period=c(2, 20, 100), ...) {
UseMethod("return.level", x)
} # end of 'return.level' function.
return.level.fevd <- function(x, return.period=c(2, 20, 100), ...) {
newcl <- tolower(x$method)
if(newcl == "gmle") newcl <- "mle"
# class(x) <- c( paste("fevd.", newcl, sep=""), class( x ) )
# UseMethod("return.level", x)
get( paste( "return.level.fevd.", newcl, sep = "" ) )( x = x, return.period = return.period, ... )
} # end of 'return.level.fevd' function.
return.level.fevd.lmoments <- function(x, return.period=c(2, 20, 100), ..., do.ci = FALSE) {
model <- x$type
p <- x$results
pnames <- names(p)
if(model=="PP") mod2 <- "GEV"
else mod2 <- model
if(!do.ci) {
if(all(is.element(c("location","shape"),pnames))) {
res <- rlevd(return.period, loc=p["location"], scale=p["scale"], shape=p["shape"], threshold=x$threshold, type=mod2, npy=x$npy, rate=x$rate)
} else if(is.element("shape", pnames)) {
res <- rlevd(return.period, scale=p["scale"], shape=p["shape"], threshold=x$threshold, type=mod2, npy=x$npy, rate=x$rate)
} else if(is.element("location", pnames)) {
res <- rlevd(return.period, loc=p["location"], scale=p["scale"], threshold=x$threshold, type=mod2, npy=x$npy, rate=x$rate)
} else res <- rlevd(return.period, scale=p["scale"], threshold=x$threshold, type=mod2, npy=x$npy, rate=x$rate)
attr(res, "return.period") <- return.period
attr(res, "data.name") <- x$data.name
attr(res, "fit.call") <- x$call
attr(res, "call") <- match.call()
attr(res, "fit.type") <- x$type
attr(res, "data.assumption") <- "stationary"
attr(res, "period") <- x$period.basis
attr(res, "units") <- x$units
attr(res, "class") <- "return.level"
} else if(do.ci) res <- ci(x, return.period=return.period, ...)
return(res)
} # end of 'return.level.fevd.lmoments' function.
return.level.fevd.bayesian <- function(x, return.period = c(2, 20, 100), ..., do.ci = FALSE, burn.in = 499, FUN = "mean", qcov = NULL, qcov.base = NULL) {
model <- x$type
if(model=="PP") mod2 <- "GEV"
else mod2 <- model
# tform <- !is.fixedfevd(x) # CJP2
if(model=="PP") mod2 <- "GEV"
else mod2 <- model
f <- match.fun(FUN)
p <- x$results
np <- dim(p)[2] - 1
p <- p[,1:np]
pnames <- colnames(p)
if(FUN=="mean") p <- colMeans(p)
else p <- apply(p, 2, f)
if(is.fixedfevd(x)) { # CJP2
if(is.element("log.scale",pnames)) {
p["log.scale"] <- exp(p["log.scale"])
pnames[pnames == "log.scale"] <- "scale"
names(p) <- pnames
}
if(!do.ci) {
if(all(is.element(c("location","shape"),pnames))) {
res <- rlevd(return.period, loc=p["location"], scale=p["scale"], shape=p["shape"], threshold=x$threshold, type=mod2, npy=x$npy, rate=x$rate)
} else if(is.element("shape", pnames)) {
res <- rlevd(return.period, scale=p["scale"], shape=p["shape"], threshold=x$threshold, type=mod2, npy=x$npy, rate=x$rate)
} else if(is.element("location", pnames)) {
res <- rlevd(return.period, loc=p["location"], scale=p["scale"], threshold=x$threshold, type=mod2, npy=x$npy, rate=x$rate)
} else res <- rlevd(return.period, scale=p["scale"], threshold=x$threshold, type=mod2, npy=x$npy, rate=x$rate)
attr(res, "return.period") <- return.period
attr(res, "data.name") <- x$data.name
attr(res, "fit.call") <- x$call
attr(res, "call") <- match.call()
attr(res, "fit.type") <- x$type
attr(res, "data.assumption") <- "stationary"
attr(res, "period") <- x$period.basis
attr(res, "units") <- x$units
attr(res, "class") <- "return.level"
} else if(do.ci) res <- ci(x, return.period=return.period, ...)
} else {
if(missing(return.period)) return.period <- 100
res <- return.level.ns.fevd.bayesian(x = x, return.period = return.period, ..., do.ci = do.ci, qcov = qcov, qcov.base = qcov.base)
if(do.ci) return(res)
if(length(return.period)==1) res <- matrix(res, ncol=1)
colnames(res) <- paste(return.period, "-", x$period.basis, " level", sep="")
attr(res, "return.period") <- return.period
attr(res, "data.name") <- x$data.name
attr(res, "fit.call") <- x$call
attr(res, "call") <- match.call()
attr(res, "fit.type") <- x$type
attr(res, "data.assumption") <- "non-stationary"
attr(res, "period") <- x$period.basis
attr(res, "units") <- x$units
if(is.null(qcov)) attr(res, "qcov") <- x$data.name[2]
else attr(res, "qcov") <- deparse(substitute(qcov))
attr(res, "class") <- "return.level"
}
return(res)
} # end of 'return.level.fevd.bayesian' function.
return.level.ns.fevd.mle <-
function (x, return.period = c(2, 20, 100), ..., alpha = 0.05,
method = c("normal"), do.ci = FALSE, verbose = FALSE, qcov = NULL,
qcov.base = NULL)
{
if(missing(return.period)) return.period <- 100
if (do.ci && method != "normal")
stop("return.level.ns.fevd.mle: only normal approximation CI calculations currently available for nonstationary return levels.")
model <- x$type
if(!(tolower(model) %in% c("pp", "gev", "weibull", "frechet", "gumbel")))
stop("return.level.ns.fevd.mle: not implemented for GP, beta, pareto, exponential models.")
if (do.ci && length(return.period) > 1 && nrow(qcov) > 1)
stop("return.level.ns.fevd.mle:: Cannot calculate confidence intervals for multiple return periods and multiple covariate values simultaneously.")
if(!do.ci) {
if (model == "PP") {
mod2 <- "GEV"
} else mod2 <- model
p <- x$results$par
pnames <- names(p)
if(is.fixedfevd(x)) stop("return.level.ns.fevd.mle: this function is for nonstationary models.")
if(is.null(qcov)) stop("return.level.ns.fevd.mle: qcov required for this function.")
if(!is.matrix(qcov)) qcov <- matrix(qcov, nrow = 1)
if(!is.qcov(qcov)) qcov <- make.qcov(x = x, vals = qcov, nr = nrow(qcov))
if(!is.null(qcov.base)) {
if(!is.matrix(qcov.base)) qcov.base <- matrix(qcov.base, nrow = 1)
if(nrow(qcov) != nrow(qcov.base) || ncol(qcov) != ncol(qcov.base))
stop("return.level.ns.fevd.mle: qcov and qcov.base must have the same number of covariates and values.")
if(!is.qcov(qcov.base)) qcov.base <- make.qcov(x = x, vals = qcov.base, nr = nrow(qcov.base))
}
loc.id <- 1:x$results$num.pars$location
sc.id <- (1 + x$results$num.pars$location):(x$results$num.pars$location + x$results$num.pars$scale)
sh.id <- (1 + x$results$num.pars$location + x$results$num.pars$scale):(x$results$num.pars$location +
x$results$num.pars$scale + x$results$num.pars$shape)
loc <- qcov[, loc.id, drop = FALSE] %*% p[loc.id]
sc <- qcov[, sc.id, drop = FALSE] %*% p[sc.id]
sh <- qcov[, sh.id, drop = FALSE] %*% p[sh.id]
if(!is.null(qcov.base)) {
loc.base <- qcov.base[, loc.id, drop = FALSE] %*% p[loc.id]
sc.base <- qcov.base[, sc.id, drop = FALSE] %*% p[sc.id]
sh.base <- qcov.base[, sh.id, drop = FALSE] %*% p[sh.id]
}
if(x$par.models$log.scale) {
sc <- exp(sc)
if(!is.null(qcov.base)) sc.base <- exp(sc.base)
}
theta <- cbind(qcov[, "threshold"], loc, sc, sh)
rlfun2 <- function(th, pd, type, npy, rate) rlevd(pd, loc = th[2],
scale = th[3], shape = th[4], threshold = th[1], type = type,
npy = npy)
res <- apply(theta, 1, rlfun2, pd = return.period, type = mod2, npy = x$npy)
if (!is.null(qcov.base)) {
theta.base <- cbind(qcov.base[, "threshold"], loc.base, sc.base, sh.base)
res <- res - apply(theta.base, 1, rlfun2, pd = return.period, type = mod2, npy = x$npy)
}
res <- t(matrix(res, nrow = length(return.period)))
colnames(res) <- paste(return.period, "-", x$period.basis,
" level", sep = "")
attr(res, "return.period") <- return.period
attr(res, "data.name") <- x$data.name
attr(res, "fit.call") <- x$call
attr(res, "call") <- match.call()
attr(res, "fit.type") <- x$type
attr(res, "data.assumption") <- "non-stationary"
attr(res, "period") <- x$period.basis
attr(res, "units") <- x$units
attr(res, "qcov") <- deparse(substitute(qcov))
attr(res, "class") <- "return.level"
if (!is.null(qcov.base)) {
attr(res, "qcov.base") <- deparse(substitute(qcov.base))
attr(res, "class") <- "return.level.diff"
}
return(res)
} else {
out <- ci.rl.ns.fevd.mle(x = x, alpha = alpha, return.period = return.period, qcov = qcov, qcov.base = qcov.base, ...)
return(out)
}
} # end of 'return.level.ns.fevd.mle' function.
return.level.fevd.mle <- function(x, return.period = c(2, 20, 100), ..., do.ci = FALSE, qcov = NULL, qcov.base = NULL) {
model <- x$type
# tform <- !is.fixedfevd(x) # removed by CJP2
if(model=="PP") mod2 <- "GEV"
else mod2 <- model
p <- x$results$par
pnames <- names(p)
if(is.fixedfevd(x)) { # CJP2
if(is.element("log.scale",pnames)) {
p["log.scale"] <- exp(p["log.scale"])
pnames[pnames == "log.scale"] <- "scale"
names(p) <- pnames
}
if(!do.ci) {
if(all(is.element(c("location","shape"),pnames))) {
res <- rlevd(return.period, loc=p["location"], scale=p["scale"], shape=p["shape"], threshold=x$threshold, type=mod2, npy=x$npy, rate=x$rate)
} else if(is.element("shape", pnames)) {
res <- rlevd(return.period, scale=p["scale"], shape=p["shape"], threshold=x$threshold, type=mod2, npy=x$npy, rate=x$rate)
} else if(is.element("location", pnames)) {
res <- rlevd(return.period, loc=p["location"], scale=p["scale"], threshold=x$threshold, type=mod2, npy=x$npy, rate=x$rate)
} else res <- rlevd(return.period, scale=p["scale"], threshold=x$threshold, type=mod2, npy=x$npy, rate=x$rate)
attr(res, "return.period") <- return.period
attr(res, "data.name") <- x$data.name
attr(res, "fit.call") <- x$call
attr(res, "call") <- match.call()
attr(res, "fit.type") <- x$type
attr(res, "data.assumption") <- "stationary"
attr(res, "period") <- x$period.basis
attr(res, "units") <- x$units
attr(res, "class") <- "return.level"
} else if(do.ci) res <- ci(x, return.period=return.period, ...)
} else {
if(missing(return.period)) return.period <- 100
if(do.ci) {
res <- ci.rl.ns.fevd.mle(x = x, return.period = return.period, ..., qcov = qcov, qcov.base = qcov.base)
} else {
if(is.null(qcov) && !is.null(qcov.base)) {
qcov <- qcov.base
qcov.base <- NULL
warning("return.level.fevd.mle: attempt to set qcov to null but not qcov.base. Setting qcov to qcov.base and qcov.base to NULL.")
}
if(is.null(qcov)) {
class(x) <- "fevd"
rlfun <- function(p, x) return(erlevd(x=x, period=p))
res <- apply(matrix(return.period, ncol=1), 1, rlfun, x=x)
} else {
res <- return(return.level.ns.fevd.mle(x = x, return.period = return.period, ..., do.ci = FALSE, qcov = qcov, qcov.base = qcov.base))
# return(return.level.ns.fevd.mle(x = x, return.period = return.period, ..., do.ci = do.ci, qcov = qcov, qcov.base = qcov.base))
#
# if(!is.matrix(qcov)) qcov <- matrix(qcov, nrow=1)
# if(!is.qcov(qcov)) qcov <- make.qcov(x=x, vals=qcov, nr=nrow(qcov))
#
# nr <- nrow(qcov)
#
# if(is.element("location", pnames)) {
#
# nloc <- 1
# loc <- p["location"]
#
# } else if(is.element("mu", substring(pnames,1,2))) {
#
# id <- substring(pnames,1,2) == "mu"
# nloc <- sum(id)
# loc <- rowSums(matrix(p[1:nloc], nr, nloc, byrow=TRUE) * qcov[,1:nloc])
#
# } else loc <- nloc <- 0
#
# if(is.element("scale", pnames)) {
#
# nsc <- 1
# scale <- p["scale"]
#
# } else if(is.element("phi", substring(pnames, 1, 3))) {
#
# id <- substring(pnames,1,3) == "phi"
# nsc <- sum(id)
# scale <- exp(rowSums(matrix(p[(nloc+1):(nloc+nsc)], nr, nsc, byrow=TRUE) * qcov[,(nloc+1):(nloc+nsc)]))
#
# } else if(is.element("sig", substring(pnames, 1, 3))) {
#
# id <- substring(pnames,1,3) == "sig"
# nsc <- sum(id)
# scale <- rowSums(matrix(p[(nloc+1):(nloc+nsc)], nr, nsc, byrow=TRUE) * qcov[,(nloc+1):(nloc+nsc)])
#
# }
#
# if(is.element("shape", pnames)) {
#
# nsh <- 1
# shape <- p["shape"]
#
# } else if(is.element("xi", substring(pnames,1,2))) {
#
# id <- substring(pnames,1,2) == "xi"
# nsh <- sum(id)
# shape <- rowSums(matrix(p[(nloc+nsc+1):(nloc+nsc+nsh)], nr, nsh, byrow=TRUE) * qcov[,(nloc+nsc+1):(nloc+nsc+nsh)])
#
# } else nsh <- shape <- 0
#
# theta <- cbind(qcov[,"threshold"], loc, scale, shape)
# rlfun2 <- function(th, pd, type, npy, rate) rlevd(return.period, loc=th[2], scale=th[3], shape=th[4], threshold=th[1], type=type, npy=npy, rate=rate)
#
# res <- t(apply(theta, 1, rlfun2, pd=return.period, type=mod2, npy=x$npy, rate=x$rate))
#
} # end of if else 'qcov' argument is NULL stmts.
if(length(return.period)==1) res <- matrix(res, ncol=1)
colnames(res) <- paste(return.period, "-", x$period.basis, " level", sep="")
attr(res, "return.period") <- return.period
attr(res, "data.name") <- x$data.name
attr(res, "fit.call") <- x$call
attr(res, "call") <- match.call()
attr(res, "fit.type") <- x$type
attr(res, "data.assumption") <- "non-stationary"
attr(res, "period") <- x$period.basis
attr(res, "units") <- x$units
if(is.null(qcov)) attr(res, "qcov") <- x$data.name[2]
else attr(res, "qcov") <- deparse(substitute(qcov))
attr(res, "class") <- "return.level"
} # end of if else 'do.ci' stmts.
} # end of if else fixed model stmts.
return(res)
} # end of 'return.level.fevd.mle' function.
print.return.level <- function(x, ...) {
tmp <- attributes(x)
print(tmp$fit.call)
print(tmp$call)
if(!is.null(tmp$units)) print(paste(tmp$fit.type, " model fitted to ", tmp$data.name[1], " (", tmp$units, ")", sep=""))
else cat("\n", tmp$fit.type, "model fitted to ", tmp$data.name, "\n")
if(tmp$fit.type=="PP") print(paste("Return levels based on GEV equivalency (i.e., return levels are for block maxima, where the blocks are ", tmp$period, "s)", sep=""))
cat("Data are assumed to be ", tmp$data.assumption, "\n")
if(tmp$data.name[2] != "") {
print(paste("Covariate data = ", tmp$data.name[2], sep=""))
if(!is.null(tmp$qcov)) if(tmp$qcov != tmp$data.name[2]) print(paste("Covariate data used for effective return levels here = ", tmp$qcov, sep=""))
}
print(paste("Return Levels for period units in ", tmp$period, "s", sep=""))
# cat("Return Period(s) = \n", tmp$return.period, "\n")
if(is.null(tmp$dim)) {
y <- c(x)
names(y) <- paste(tmp$return.period, "-", tmp$period, " level", sep="")
} else {
y <- matrix(c(x), tmp$dim[1], tmp$dim[2])
colnames(y) <- paste(tmp$return.period, "-", tmp$period, " level", sep="")
}
print(y)
invisible()
} # end of 'print.return.level' function.
make.qcov <- function(x, vals, nr=1, ...) {
if(x$method != "Bayesian") p <- x$results$par
else p <- apply(x$results[, -dim(x$results)[2] ], 2, mean, na.rm=TRUE)
np <- length(p)
pnames <- names(p)
if(is.element("mu", substring(pnames, 1, 2))) nloc <- sum(substring(pnames, 1, 2) == "mu")
else if(is.element("location", pnames)) nloc <- 1
else nloc <- 0
if(is.element("phi", substring(pnames, 1, 3))) nsc <- sum(substring(pnames, 1, 3) == "phi")
else if(is.element("sig", substring(pnames, 1, 3))) nsc <- sum(substring(pnames, 1, 3) == "sig")
else nsc <- 1
if(is.element("shape", pnames)) nsh <- 1
else if(is.element("xi", substring(pnames, 1, 2))) nsh <- sum(substring(pnames, 1, 2) == "xi")
else nsh <- 0
if(missing(vals)) {
out <- matrix(0, nr, np+1)
out[,1] <- 1
if(nloc > 0) out[,nloc+1] <- 1
out[,nloc+nsc+1] <- 1
if(!is.null(x$threshold)) {
if(length(x$threshold) >= nr) out[,np+1] <- x$threshold[1:nr]
else out[,np+1] <- x$threshold[1]
}
colnames(out) <- c(pnames, "threshold")
return(out)
}
if(is.list(vals)) {
vnames <- names(vals)
if(is.null(vnames)) stop("make.qcov: vals must be a named list.")
if(!all(is.element(vnames, c(pnames, "threshold")))) stop("make.qcov: names of vals list must match parameter names.")
if(missing(nr)) {
nv <- lapply(vals, length)
if(length(unique(nv)) != 1) stop("make.qcov: Sorry, this function is limited. Length of each component must be the same or nr must be specified.")
nr <- length(vals[[1]])
}
out <- matrix(NA, nrow=nr, ncol=np + 1)
for(i in 1:np) {
if(is.element(pnames[i], vnames)) {
id <- (1:length(vnames))[ vnames == pnames[i] ]
out[,i] <- vals[[ id ]]
} else out[,i] <- 1
} # end of for 'i' loop.
if(is.element("threshold", vnames)) {
id <- (1:length(vnames))[ vnames == "threshold" ]
out[,np + 1] <- vals[[ id ]]
} else if(!is.element(x$type, c("GEV", "Gumbel", "Frechet"))) out[,np + 1] <- x$threshold[1]
} else if(is.matrix(vals)) out <- vals
else if(is.numeric(vals)) out <- matrix(vals, nrow=nr, ...)
else stop("make.qcov: invalid vals argument.")
if(dim(out)[2]==np) {
if(!is.null(x$threshold)) {
if(length(x$threshold) >= nr) out <- cbind(out, x$threshold[1:nr])
else out <- cbind(out, x$threshold[1])
} else out <- cbind(out, 0)
colnames(out) <- c(pnames, "threshold")
} else if(dim(out)[2]==np+1) {
if(is.null(colnames(out))) colnames(out) <- c(pnames, "threshold")
} else stop("make.qcov: length of (or number of columns of) vals must be equal to number of model parameters, np, or np + 1.")
# Some final checks.
if(is.element("location", pnames)) {
if(!all(out[,"location"]==1)) {
warning("make.qcov: invalid qcov values for fixed location parameter. Re-setting to one.")
out[,"location"] <- 1
}
} else if(is.element("mu0", pnames)) {
if(!all(out[,"mu0"] == 1)) {
warning("make.qcov: invalid qcov values for fixed mu0 parameter. Re-setting to one.")
out[,"mu0"] <- 1
}
}
if(any(is.element(c("log.scale","scale"), pnames))) {
id <- (pnames == "scale") | (pnames == "log.scale")
id <- c(id, FALSE)
if(!all(out[,id] == 1)) {
warning("make.qcov: invalid qcov values for fixed scale parameter. Re-setting to one.")
out[,id] <- 1
}
} else if(any(is.element(c("sig0","phi0"), pnames))) {
id <- (pnames == "sig0") | (pnames == "phi0")
id <- c(id, FALSE)
if(!all(out[,id] == 1)) {
warning("make.qcov: invalid qcov values for fixed scale intercept term. Re-setting to one.")
out[,id] <- 1
}
}
if(is.element("shape", pnames)) {
if(!all(out[,"shape"] == 1)) {
warning("make.qcov: invalid qcov values for fixed shape parameter. Re-setting to one.")
out[,"shape"] <- 1
}
} else if(is.element("xi0", pnames)) {
if(!all(out[,"xi0"] == 1)) {
warning("make.qcov: invalid qcov values for fixed xi0 parameter. Re-setting to one.")
out[,"xi0"] <- 1
}
}
return(out)
} # end of 'make.qcov' function.
is.qcov <- function(x) {
if(!is.matrix(x)) return(FALSE)
d <- dim(x)
if(is.null(d)) return(FALSE)
dnames <- colnames(x)
if(is.null(dnames)) return(FALSE)
if(dnames[d[2]] != "threshold") return(FALSE)
if(d[2] > 1) return(TRUE)
else return(FALSE)
} # end of 'is.qcov' function.
probprob.plot.evd <- function(xp, y, model, loc, scale, shape, u, tform=FALSE, eid, obj, ytrans=NULL, npy=NULL, ...) {
args <- list(...)
if(is.null(args$main)) m1 <- deparse(obj$call)
if(!tform) {
if(is.element(model, c("PP","GP","Beta","Pareto","Exponential"))) {
yp <- pevd(sort(y[eid]), loc=loc, scale=scale, shape=shape, threshold=u, npy=npy, type=model)
} else yp <- pevd(sort(y), loc=loc, scale=scale, shape=shape, threshold=u, npy=npy, type=model)
if(is.null(args$main)) plot(xp, yp, main=m1, xlab="Empirical Probabilities", ylab="Model Probabilities", ...)
else plot(xp, yp, xlab="Empirical Probabilities", ylab="Model Probabilities", ...)
abline(0,1)
} else {
# yp <- pevd(sort(ytrans), loc=0, scale=1, shape=0, threshold=0, npy=npy, type=model)
if(is.element(model, c("GEV","Gumbel","Weibull","Frechet"))) yp <- exp(-exp(-sort(ytrans)))
else if(model=="PP") yp <- sort(ytrans)
else yp <- 1 - exp(-sort(ytrans))
if(is.null(args$main)) {
plot(xp, yp, main=m1, xlab="Residual Empirical Probabilities", ylab="Residual Model Probabilities", ...)
} else plot(xp, yp, xlab="Residual Empirical Probabilities", ylab="Residual Model Probabilities", ...)
abline(0,1)
}
invisible( data.frame( empirical = xp, model = yp ) )
} # end of 'probprob.plot.evd' stmts.
quantquant.plot.evd <- function(x, xp, y, u, loc, scale, shape, tform=FALSE, eid, ytrans=NULL,
model= c("GEV", "GP", "PP", "Gumbel", "Frechet", "Weibull", "Exponential", "Beta", "Pareto"),
type=c("primary","qq"), ...) {
args <- list(...)
type <- match.arg(type)
if(!tform) {
if(is.element(model, c("Weibull","Frechet"))) mod2 <- "GEV"
else mod2 <- model
if(!is.element(model,c("PP","GP","Beta","Pareto"))) yq <- qevd(xp, loc=loc, scale=scale, shape=shape, type=mod2)
else yq <- qevd(xp, threshold=u, loc=loc, scale=scale, shape=shape, type=mod2)
if(is.null(args$main)) {
if(type=="primary") m2 <- ""
else m2 <- deparse(x$call)
if(is.element(model, c("GEV","Weibull","Frechet","Gumbel"))) {
plot(yq, sort(y), xlab="Model Quantiles", ylab="Empirical Quantiles", main=m2)
if( type == "qq" ) out <- data.frame( empirical = sort( y ), model = yq )
} else {
plot(yq, sort(y[eid]), xlab="Model Quantiles", ylab="Empirical Quantiles", main=m2, ...)
if( type == "qq" ) out <- data.frame( empirical = sort( y[ eid ] ), model = yq )
}
} else {
if(is.element(model, c("GEV","Weibull","Frechet","Gumbel"))) {
plot(yq, sort(y), xlab="Model Quantiles", ylab="Empirical Quantiles", ...)
if( type == "qq" ) out <- data.frame( empirical = sort( y ), model = yq )
} else {
plot(yq, sort(y[eid]), xlab="Model Quantiles", ylab="Empirical Quantiles", ...)
if( type == "qq" ) out <- data.frame( empirical = sort( y[ eid ] ), model = yq )
}
}
} else {
if(is.null(args$main)) {
if(type=="primary") m2 <- ""
else m2 <- deparse(x$call)
if(is.element(model, c("GEV","Weibull","Frechet"))) m2 <- paste(m2, "(Gumbel Scale)", sep="\n")
else if(is.element(model, c("PP", "GP", "Beta", "Pareto"))) m2 <- paste(m2, "Exponential Scale", sep="\n")
if( is.element(model, c("GEV","Weibull","Gumbel","Frechet")) ) {
plot(-log(-log(sort(xp))), sort(ytrans), main=m2, xlab="(Standardized) Model Quantiles",
ylab="Empirical Residual Quantiles", ...)
if( type == "qq" ) out <- data.frame( empirical = sort( ytrans ), model = -log(-log(sort(xp))) )
} else if(is.element(model, c("GP","Beta","Exponential","Pareto"))) {
plot(-log(1 - xp), sort(ytrans), main=m2, xlab="(Standardized) Residual Quantiles",
ylab="Empirical Residual Quantiles", ...)
if( type == "qq" ) out <- data.frame( empirical = sort( ytrans ), model = -log(1 - xp) )
} else if(model=="PP") {
plot(-log(1 - xp), sort(-log(ytrans)), main=m2, xlab="(Standardized) Residual Quantiles",
ylab="Empirical Residual Quantiles", ...)
if( type == "qq" ) out <- data.frame( empirical = sort(-log(ytrans)), model = -log(1 - xp) )
}
} else {
if(is.element(model, c("GEV","Weibull","Gumbel","Frechet"))) {
plot(-log(-log(sort(xp))), sort(ytrans), xlab="Model", ylab="Empirical", ...)
if( type == "qq" ) out <- data.frame( empirical = sort(ytrans), model = -log(-log(sort(xp))) )
} else if(is.element(model, c("GP","Beta","Exponential","Pareto"))) {
plot(-log(1 - xp), sort(ytrans), xlab="Model", ylab="Empirical", ...)
if( type == "qq" ) out <- data.frame( empirical = sort( ytrans ), model = -log(1 - xp) )
} else if(model=="PP") {
plot(-log(1 - xp), sort(-log(ytrans)), xlab="Model", ylab="Empirical", ...)
if( type == "qq" ) out <- data.frame( empirical = sort(-log(ytrans)), model = -log(1 - xp) )
}
}
}
abline(0,1)
if( type == "qq" ) invisible( out )
else invisible()
} # end of 'quantquant.plot.evd' funciton.
quantquant2.plot.evd <- function(x, y, eid, model=c("GEV", "GP", "PP", "Gumbel", "Frechet", "Weibull", "Exponential", "Beta", "Pareto"), type=c("primary","qq2"), ...) {
args <- list(...)
type <- match.arg(type)
if(is.fixedfevd(x)) z <- rextRemes(x, n=x$n)
else {
if(!is.element(model, c("PP","GP","Beta","Exponential","Pareto"))) z <- rextRemes(x)
else z <- rextRemes( x, n = sum( eid, na.rm = TRUE ) )
}
if(!is.element(model, c("PP","GP","Beta","Exponential","Pareto"))) {
yQQ <- y
if(is.null(args$xlab)) xl <- paste(x$data.name[1], " Empirical Quantiles", sep="")
else xl <- args$xlab
} else {
yQQ <- y[eid]
if(is.null(args$xlab)) {
if(length(x$threshold)==1) xl <- paste(x$data.name[1], "( > ", x$threshold, ") Empirical Quantiles", sep="")
else xl <- paste(x$data.name[1], "( > threshold) Empirical Quantiles", sep="")
} else xl <- args$xlab
}
if(!is.null(args$main)) {
if(type=="primary") mQQ <- ""
else mQQ <- deparse(x$call)
if( type != "qq2" ) qqplot(yQQ, z, main=mQQ, xlab=xl, ylab="Quantiles from Model Simulated Data" )
else out <- qqplot(yQQ, z, main=mQQ, xlab=xl, ylab="Quantiles from Model Simulated Data" )
} else {
if( type != "qq2" ) qqplot(yQQ, z, xlab=xl, ylab="Quantiles from Model Simulated Data")
else out <- qqplot(yQQ, z, xlab=xl, ylab="Quantiles from Model Simulated Data")
}
if( type == "qq2" ) invisible( out )
else invisible()
} # end of 'quantquant2.plot.evd' function.
histplot.evd <- function(x, y, u=u, loc=loc, scale=scale, shape=shape, ytrans=NULL, tform=FALSE, eid,
model= c("GEV", "GP", "PP", "Gumbel", "Frechet", "Weibull", "Exponential", "Beta", "Pareto"), hist.args=NULL, ...) {
args <- list(...)
if(!tform) {
if(is.null(args$main)) {
if(model != "PP") m4 <- paste(deparse(x$call), "Histogram", sep="\n")
else m4 <- paste(deparse(x$call), "\n", paste("Histogram (", x$period.basis, " maxima)", sep=""))
}
if(is.element(model, c("GP", "Beta", "Exponential", "Pareto"))) yh <- y[eid]
else if(model != "PP") yh <- y
else {
blocks <- rep(1:x$span, each=x$npy)
n2 <- length(blocks)
if(n2 < x$n) blocks <- c(blocks, rep(blocks[n2], x$n - n2))
else if(n2 > x$n) blocks <- blocks[1:x$n]
yh <- c(aggregate(y, by=list(blocks), max)$x)
}
} else { # Eric 8/14/13
if(is.element(model, c("PP", "GP", "Exponential", "Beta", "Pareto"))) stop("plot.fevd: invalid type argument for this model.")
if(is.null(args$main)) m4 <- paste(deparse(x$call), "Histogram of Transformed Data", sep="\n")
# if(model != "PP") yh <- ytrans
# else {
# blocks <- rep(1:x$span, each=x$npy)
# n2 <- length(blocks)
# if(n2 < x$n) blocks <- c(blocks, rep(blocks[n2], x$n - n2))
# else if(n2 > x$n) blocks <- blocks[1:x$n]
# yh <- c(aggregate(ytrans, by=list(blocks), max)$x)
# } # end of if else 'model != PP' stmts.
}
if(is.null(hist.args)) {
if(!is.null(args$ylim)) {
if(is.null(args$main)) {
if(is.null(args$col)) out <- hist(yh, col="darkblue", freq=FALSE, breaks="FD", xlab=x$data.name[1], main=m4, ...)
else out <- hist(yh, freq=FALSE, breaks="FD", main=m4, xlab=x$data.name[1], ...)
} else {
if(is.null(args$col)) out <- hist(yh, col="darkblue", freq=FALSE, breaks="FD", xlab=x$data.name[1], ...)
else out <- hist(yh, freq=FALSE, breaks="FD", xlab=x$data.name[1], ...)
}
} else {
if(is.null(args$main)) {
if(is.null(args$col)) out <- hist(yh, col="darkblue", freq=FALSE, breaks="FD", xlab=x$data.name[1], main=m4, ylim=c(0,1.5), ...)
else out <- hist(yh, freq=FALSE, breaks="FD", main=m4, xlab=x$data.name[1], ylim=c(0,1.5), ...)
} else {
if(is.null(args$col)) out <- hist(yh, col="darkblue", freq=FALSE, breaks="FD", xlab=x$data.name[1], ylim=c(0,1.5), ...)
else out <- hist(yh, freq=FALSE, breaks="FD", xlab=x$data.name[1], ylim=c(0,1.5), ...)
}
}
} else out <- do.call("hist", c(list(yh), hist.args))
xh <- seq(min(yh, na.rm=TRUE), max(yh, na.rm=TRUE),,100)
if(is.element(model, c("Gumbel", "Weibull", "Frechet"))) mod2 <- "GEV"
else if(is.element(model, c("PP", "Pareto", "Frechet", "Beta", "Exponential"))) mod2 <- "GP"
else mod2 <- model
if(mod2 != "GP") ymod <- devd(xh, loc=loc, scale=scale, shape=shape, threshold=u, type=mod2)
else ymod <- devd(xh - u, loc=loc, scale=scale, shape=shape, threshold=u, type=mod2)
lines(xh, ymod, lty=2, col="blue", lwd=1.5)
invisible( out )
} # end of 'hitsplot.evd' function.
densplot.evd <- function(x, y, u, loc, scale, shape, tform=FALSE, eid, ytrans=NULL,
model= c("GEV", "GP", "PP", "Gumbel", "Frechet", "Weibull", "Exponential", "Beta", "Pareto"), density.args=NULL,
type=c("primary","density"), leg=TRUE, ...) {
args <- list(...)
type <- match.arg(type)
if(!tform) {
if(!is.element(model, c("PP","GP","Beta","Exponential","Pareto"))) yd <- y
else if(model != "PP") yd <- y[eid]
else {
blocks <- rep(1:x$span, each=x$npy)
n2 <- length(blocks)
if(n2 < x$n) blocks <- c(blocks, rep(blocks[n2], x$n - n2))
else if(n2 > x$n) blocks <- blocks[1:x$n]
ytmp <- c(datagrabber(x, cov.data = FALSE))
ytmp2 <- c(aggregate(ytmp, by = list(blocks), max)$x)
mbid <- logical(x$n)
for(i in 1:(length(unique(blocks)))) {
tmpind <- blocks == blocks[i]
tmpind2 <- ytmp == ytmp2[i]
tmpind2[is.na(tmpind2)] <- FALSE
finind <- tmpind & tmpind2
if(sum(finind) > 1) {
numsind <- (1:x$n)[finind]
numsind <- numsind[1]
finind[-numsind] <- FALSE
}
mbid[finind] <- TRUE
} # end of for 'i' loop.
yd <- y[mbid]
# yd <- c(aggregate(y, by=list(blocks), max)$x)
}
if(is.null(density.args)) yd <- density(yd)
else yd <- do.call("density", c(list(yd), density.args))
if(is.null(args$ylim)) {
yld <- range(yd$y, finite=TRUE)
yld[1] <- min(yld[1], 0)
# yld[2] <- max(yld[2]+0.5, 1)
}
if(is.null(args$main)) {
if(type=="primary") m3 <- ""
else m3 <- deparse(x$call)
if(is.null(args$ylim)) plot(yd, main=m3, ylim=yld, ...)
else plot(yd, main=m3, ...)
} else {
if(is.null(args$ylim)) {
plot(yd, ylim=yld, ...)
} else plot(yd, ...)
}
if(is.element(model, c("PP", "Gumbel", "Weibull", "Frechet"))) mod2 <- "GEV"
else if(is.element(model, c("Pareto", "Frechet", "Beta", "Exponential"))) mod2 <- "GP"
else mod2 <- model
xd <- seq(min(yd$x, na.rm=TRUE), max(yd$x, na.rm=TRUE),,100)
if(mod2 != "GP") ymod <- devd(xd, loc=loc, scale=scale, shape=shape, threshold=u, type=mod2)
else ymod <- devd(xd - u, loc=loc, scale=scale, shape=shape, threshold=u, type=mod2)
lines(xd, ymod, lty=2, col="blue", lwd=1.5)
if(leg) legend("topright", legend=c("Data","Model"), col=c("black","blue"), lty=c(1,2), lwd=c(1,1.5), bty="n")
} else {
# if(model != "PP") ytrans2 <- ytrans
# else {
# blocks <- rep(1:x$span, each=x$npy)
# n2 <- length(blocks)
# if(n2 < x$n) blocks <- c(blocks, rep(blocks[n2], x$n - n2))
# else if(n2 > x$n) blocks <- blocks[1:x$n]
# ytrans2 <- c(aggregate(ytrans, by=list(blocks), max)$x)
# }
if(is.null(density.args)) yd <- density(ytrans)
else yd <- do.call("density", c(list(ytrans), density.args))
if(is.null(args$ylim)) {
yld <- range(yd$y, finite=TRUE)
yld[1] <- min(yld[1], 0)
# yld[2] <- max(yld[2]+0.5, 1)
}
if(is.null(args$main)) {
if(type=="primary") m3 <- ""
else m3 <- paste(deparse(x$call), "Transformed", sep="\n")
if(is.null(args$ylim)) plot(yd, main=m3, ylim=yld, ...)
else plot(yd, main=m3, ...)
} else {
if(is.null(args$ylim)) plot(yd, ylim=yld, ...)
else plot(yd, ...)
}
if(is.element(model, c("Gumbel", "Weibull", "Frechet"))) mod2 <- "GEV"
else if(is.element(model, c("PP", "Pareto", "Frechet", "Beta", "Exponential"))) mod2 <- "GP"
else mod2 <- model
xd <- seq(min(yd$x, na.rm=TRUE), max(yd$x, na.rm=TRUE),,100)
ymod <- devd(xd, loc=0, scale=1, shape=0, threshold=0, type=mod2)
lines(xd, ymod, lty=2, col="blue", lwd=1.5)
if(leg) legend("topright", legend=c("Transformed Data","Standardized Model"), col=c("black","blue"), lty=c(1,2), lwd=c(1,1.5), bty="n")
} # end of if else '!tform' stmts.
# if(!is.element(model, c("pp","gp","beta","exponential","pareto"))) points(y, rep(0, length(y)), pch="|", cex=0.5)
# else points(y[eid], rep(0, length(y[eid])), pch="|", cex=0.5)
invisible( yd )
} # end of 'densplot.evd' function.
rlplot.evd <- function(x, xp, y, u, eid, rperiods, tform=FALSE, model= c("GEV", "GP", "PP", "Gumbel", "Frechet", "Weibull", "Exponential", "Beta", "Pareto"),
type=c("primary","rl"), leg=TRUE, a=0, ...) {
args <- list(...)
type <- match.arg(type)
if(x$method != "Lmoments") { # Eric -- 8/27/13
const.thresh <- check.constant(x$par.models$threshold)
const.loc <- check.constant(x$par.models$location)
const.scale <- check.constant(x$par.models$scale)
const.shape <- check.constant(x$par.models$shape)
if(is.element(model, c("PP", "GP", "Exponential", "Beta", "Pareto")) && !const.thresh && all(c(const.loc, const.scale, const.shape)) && type == "rl")
stop("rlplot.evd: invalid type argument for POT models with varying thresholds but constant parameters (are you sure you about this model choice?).")
}
if(is.null(args$main)) {
if(type=="primary") m5 <- ""
else m5 <- deparse(x$call)
}
if(model=="PP") {
if(!tform) mod2 <- "GEV"
else mod2 <- "GEV"
if(is.null(args$main)) {
if(!tform) {
if(type=="rl") m5 <- paste(m5, "Return Levels based on approx. equivalent GEV", sep="\n")
else m5 <- paste("Return Levels based on approx.", "equivalent GEV", sep="\n")
} else {
if(type=="rl") m5 <- paste(m5, "Return Levels based on approx. equivalent GEV df", sep="\n")
else m5 <- paste("Return Levels based on approx.", "equivalent GEV df", sep="\n")
}
}
blocks <- rep(1:x$span, each=x$npy)
n2 <- length(blocks)
if(n2 < x$n) blocks <- c(blocks, rep(blocks[n2], x$n - n2))
else if(n2 > x$n) blocks <- blocks[1:x$n]
yEmp <- c(aggregate(y, by=list(blocks), max)$x)
} else mod2 <- model
if(!tform) {
bds <- ci(x, return.period=rperiods)
yrl <- bds[,2]
if(is.null(args$ylim)) yl <- range(c(bds), finite=TRUE)
if(is.element(model, c("PP", "GEV", "Gumbel", "Weibull", "Frechet"))) xrl <- -1/(log(1 - 1/rperiods))
else xrl <- rperiods
if(is.null(x$units)) ylb <- "Return Level"
else ylb <- paste("Return Level (", x$units, ")", sep="")
xlb <- paste("Return Period (", x$period.basis, "s)", sep="")
if(is.null(args$main)) {
if(!is.null(args$ylim)) plot(xrl, yrl, type="l", log="x", xlab=xlb, ylab=ylb, main=m5, ...)
else plot(xrl, yrl, type="l", log="x", ylim=yl, xlab=xlb, ylab=ylb, main=m5, ...)
} else {
if(!is.null(args$ylim)) plot(xrl, yrl, type="l", log="x", xlab=xlb, ylab=ylb, ...)
else plot(xrl, yrl, type="l", log="x", ylim=yl, xlab=xlb, ylab=ylb, ...)
}
lines(xrl, bds[,1], col="gray", lty=2, lwd=2)
lines(xrl, bds[,3], col="gray", lty=2, lwd=2)
if(is.element(model, c("GEV", "Gumbel", "Weibull", "Frechet"))) points(-1/log(xp), sort(y))
else if(is.element(model, c("GP", "Beta", "Pareto", "Exponential"))) {
n2 <- x$n
if(is.null(a)) xp2 <- ppoints(n2)
else xp2 <- ppoints(n2, a=a)
sdat <- sort(y)
points(-1/log(xp2)[sdat > u]/x$npy, sdat[sdat > u])
out <- list( model = bds, empirical = data.frame( transformed.points = -1/log(xp2)[sdat > u]/x$npy,
sorted.level = sdat[sdat > u] ) )
} else if(model == "PP") {
if(is.null(a)) xp2 <- ppoints(length(yEmp))
else xp2 <- ppoints(length(yEmp), a=a)
points(-1/log(xp2), sort(yEmp))
out <- list( model = bds, empirical = data.frame( transformed.points = -1/log(xp2),
sorted.level = sort( yEmp ) ) )
}
} else {
np <- length(rperiods)
n <- length(y)
effrl <- matrix(NA, n, np)
for(i in 1:np) effrl[,i] <- return.level(x, return.period=rperiods[i])
if(is.null(args$ylim)) {
yl <- range(c(y, c(effrl)), finite=TRUE)
yl[2] <- yl[2] + sign(yl[2]) * log2(abs(yl[2]))
if(is.null(args$main)) {
if(!is.element(model, c("GP", "Beta", "Exponential", "Pareto"))) plot(y, type="l", xlab="index", main=m5, ylim=yl, ...)
else plot(y[eid], type = "l", xlab = "index", ylim = yl, main = m5, ...)
} else {
if(!is.element(model, c("GP", "Beta", "Exponential", "Pareto"))) plot(y, type="l", xlab="index", ...)
else plot(y[eid], type="l", xlab="index", ylim=yl, ...)
}
} else {
if(is.null(args$main)) {
if(!is.element(model, c("GP", "Beta", "Exponential", "Pareto"))) plot(y, type="l", xlab="index", main=m5, ...)
else plot(y[eid], type = "l", xlab = "index", main = m5, ...)
} else {
if(!is.element(model, c("GP", "Beta","Exponential","Pareto"))) plot(y, type="l", xlab="index", ...)
else plot(y[eid], type="l", xlab="index", ...)
}
} # end of if else 'ylim passed via '...' stmts.
for(i in 1:np) lines(effrl[,i], lty=i, col=i+1)
out <- list( series = y[ eid ], level = effrl )
if(is.element(model, c("PP", "GP", "Beta", "Exponential", "Pareto"))) {
if(length(u) > 1 && all(u == u[1])) u <- u[1]
if(length(u) == 1) abline(h=u, col="darkorange")
else lines(u, col="darkorange")
if(leg) legend("topleft", legend=c(paste(rperiods, "-", x$period.basis, " level", sep=""), "threshold"), lty=c(1:np,1), col=c(2:(np+1),"darkorange"), bg="white")
} else if(leg) legend("topleft", legend=paste(rperiods, "-", x$period.basis, " level", sep=""), lty=1:np, col=2:(np+1), bg="white")
} # end of if else '!tform' stmts.
if( type == "rl" ) invisible( out )
else invisible()
} # end of 'rlplot.evd' function.
eeplot <- function(x, type = c("both", "Zplot", "Wplot"), set.pw = FALSE, d = NULL, ...) {
type <- match.arg(type)
theCall <- match.call()
if(x$type != "PP") stop("eeplot: Sorry, this funciton is for PP fits only.")
y <- datagrabber(x, cov.data = FALSE)
u <- x$threshold
p <- findpars(x)
Tk <- (1:x$n)[y > u]
out <- Tk
if(set.pw) {
op <- par()
if(type == "both") par(mfrow = c(1, 2), oma = c(0, 0, 2, 0))
else par(mfrow = c(1, 1), oma = c(0, 0, 2, 0))
}
if(is.element(type, c("both", "Zplot"))) {
lambda <- (1 + p$shape * (u - p$location)/p$scale)^(-1/p$shape)
if(is.null(d)) {
if(x$time.units == "days") {
if(x$period.basis == "year") d <- 365.25
else stop("eeplot: invalid period basis, perhaps set the d argument.")
} else if(x$time.units == "months") {
if(x$period.basis == "year") d <- 12
else stop("eeplot: invalid period basis, perhaps set the d argument.")
} else if(x$time.units == "years") {
if(x$period.basis == "year") d <- 1
else stop("eeplot: invalid period basis, perhaps set the d argument.")
} else if(x$time.units == "hours") {
if(x$period.basis == "year") d <- 8766 # 24 * 365.25
else stop("eeplot: invalid period basis, perhaps set the d argument.")
} else if(x$time.units == "minutes") {
if(x$period.basis == "year") d <- 525960 # 60 * 24 * 365.25
else stop("eeplot: invalid period basis, perhaps set the d argument.")
} else if(x$time.units == "seconds") {
if(x$period.basis == "year") d <- 31557600 # 60 * 60 * 24 * 365.25
else stop("eeplot: invalid period basis, perhaps set the d argument.")
} else {
tmp.units <- unlist(strsplit(x$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") d <- numper
else if(pertiid=="month") d <- numper*12
else if(pertiid=="day") d <- numper*365.25
else if(pertiid == "hour") d <- numper * 24 * 365.25
else if(pertiid == "minute") d <- numper * 60 * 24 * 365.25
else if(pertiid == "second") d <- numper * 60 * 60 * 24 * 365.25
}
}
lambda <- lambda/d
if(length(lambda) == 1) lambda <- rep(lambda, x$n)
Zk <- diff(c(0, cumsum(lambda)[Tk]))
out <- cbind(out, Zk)
colnames(out) <- c("Exceed Time", "Zk")
# qqPlot(Zk, distribution = "exp", xlab = "Expected Values\nUnder exponential(1)", ylab = "Observed Z_k Values",
# col.lines = "grey", grid = FALSE, ...)
qqplot( qexp( (1:length(Zk) - 0.5)/length(Zk) ), Zk, xlab = "Expected Values\nUnder exponential(1)", ylab = "Observed Z_k Values" )
abline(0, 1, col = "darkorange", lty = 2)
legend("topleft", legend = c("1-1 line", "regression line", "95% confidence bands"), col = c("darkorange", "grey", "grey"), lty = c(2, 1, 2))
}
if(is.element(type, c("both", "Wplot"))) {
if(length(u) == 1) u <- rep(u, x$n)
Wk <- (1/p$shape[Tk]) * log(1 + p$shape[Tk] * (y[Tk] - u[Tk])/(p$scale[Tk] + p$shape[Tk] * (u[Tk] - p$location[Tk])))
# qqPlot(Wk, distribution = "exp", xlab = "Expected Values\nUnder exponential(1)", ylab = "Observed W_k Values",
# col.lines = "grey", grid = FALSE, ...)
qqplot( qexp( (1:length(Wk) - 0.5)/length(Wk) ), Wk, xlab = "Expected Values\nUnder exponential(1)", ylab = "Observed W_k Values" )
abline(0, 1, col = "darkorange", lty = 2)
legend("topleft", legend = c("1-1 line", "regression line", "95% confidence bands"), col = c("darkorange", "grey", "grey"), lty = c(2, 3, 1))
out <- cbind(out, Wk)
if(is.matrix(out)) colnames(out) <- c("Exceed Time", "Zk", "Wk")
else colnames(out) <- c("Exceed Time", "Wk")
}
if(set.pw) {
mtext(deparse(x$call), line=0.5, outer=TRUE)
par( mfrow = op$mfrow, oma = op$oma )
}
attr(out, "call") <- theCall
class(out) <- "ee"
invisible(out)
} # end of 'eeplot' function.
plot.ee <- function(x, pch2 = "+", col2 = "gray", xlab = "Time of Exceeding Threshold", ylab = "", ...) {
args <- list(...)
if(is.null(args$pch)) pch1 <- "o"
else pch1 <- args$pch
if(is.null(args$col)) col1 <- 1
else col1 <- args$col
if(is.null(args$ylim)) {
yl <- range(x[,-1], finite = TRUE)
plot(x[,1], x[,2], ylim = yl, xlab = xlab, ylab = ylab, ...)
} else plot(x[,1], x[,2], xlab = xlab, ylab = ylab, ...)
if(dim(x)[2] == 3) {
points(x[,1], x[,3], pch = pch2, col = col2)
legend("topleft", legend = c("Zk", "Wk"), pch = c(pch1, pch2), col = c(col1, col2), bty = "n")
}
if(is.null(args$main)) title(attributes(x)$call)
invisible( x )
} # end of 'plot.ee' function.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.