# R/nonpar.R In evd: Functions for Extreme Value Distributions

```"abvnonpar"<-
function(x = 0.5, data, epmar = FALSE, nsloc1 = NULL, nsloc2 = NULL,
method = c("cfg","pickands","tdo","pot"), k = nrow(data)/4,
convex = FALSE, rev = FALSE, madj = 0, kmar = NULL, plot = FALSE,
add = FALSE, lty = 1, lwd = 1, col = 1, blty = 3, blwd = 1,
xlim = c(0,1), ylim = c(0.5,1), xlab = "t", ylab = "A(t)", ...)
{
if(mode(x) != "numeric" || any(x < 0, na.rm=TRUE) ||
any(x > 1, na.rm=TRUE)) stop("invalid argument for `x'")
method <- match.arg(method)

# Empirical transform to exponential
epdata <- apply(data, 2, rank, na.last = "keep")
nasm <- apply(data, 2, function(x) sum(!is.na(x)))
epdata <- epdata / rep(nasm+1, each = nrow(data))
epdata <- -log(epdata)
if(epmar) data <- epdata
# End empirical transform

if(!epmar) {

if(method == "pot") {
# Parametric pot transform to exponential
if(any(k >= nasm)) stop("k is too large")
u1 <- sort(data[,1], decreasing = TRUE)[k+1]
u2 <- sort(data[,2], decreasing = TRUE)[k+1]
d1ab <- (data[,1] > u1) & !is.na(data[,1])
d2ab <- (data[,2] > u2) & !is.na(data[,2])
if(!is.null(kmar)) {
data[d1ab,1] <- mtransform(data[d1ab,1], c(u1, kmar))
data[d2ab,2] <- mtransform(data[d2ab,2], c(u2, kmar))
}
else {
mle.m1 <- c(u1, fitted(fpot(data[d1ab,1], threshold = u1)))
mle.m2 <- c(u2, fitted(fpot(data[d2ab,2], threshold = u2)))
data[d1ab,1] <- mtransform(data[d1ab,1], mle.m1)
data[d2ab,2] <- mtransform(data[d2ab,2], mle.m2)
}
data[d1ab,1] <- -log(1 - k * data[d1ab,1] / nasm[1])
data[d2ab,2] <- -log(1 - k * data[d2ab,2] / nasm[2])
data[!d1ab, 1] <- epdata[!d1ab, 1]
data[!d2ab, 2] <- epdata[!d2ab, 2]
# End parametric pot transform
}

if(method != "pot") {
# Parametric gev transform to exponential
if(!is.null(kmar)) {
data <- mtransform(data, kmar)
}
else {
if(!is.null(nsloc1)) {
if(is.vector(nsloc1)) nsloc1 <- data.frame(trend = nsloc1)
if(nrow(nsloc1) != nrow(data))
stop("`nsloc1' and data are not compatible")
nslocmat1 <- cbind(1,as.matrix(nsloc1))
}
if(!is.null(nsloc2)) {
if(is.vector(nsloc2)) nsloc2 <- data.frame(trend = nsloc2)
if(nrow(nsloc2) != nrow(data))
stop("`nsloc2' and data are not compatible")
nslocmat2 <- cbind(1,as.matrix(nsloc2))
}
mle.m1 <- fitted(fgev(data[,1], nsloc = nsloc1, std.err = FALSE))
loc.mle.m1 <- mle.m1[grep("^loc", names(mle.m1))]
if(is.null(nsloc1)) loc.mle.m1 <- rep(loc.mle.m1, nrow(data))
else loc.mle.m1 <- nslocmat1 %*% loc.mle.m1
mle.m1 <- cbind(loc.mle.m1, mle.m1["scale"], mle.m1["shape"])
mle.m2 <- fitted(fgev(data[,2], nsloc = nsloc2, std.err = FALSE))
loc.mle.m2 <- mle.m2[grep("^loc", names(mle.m2))]
if(is.null(nsloc2)) loc.mle.m2 <- rep(loc.mle.m2, nrow(data))
else loc.mle.m2 <- nslocmat2 %*% loc.mle.m2
mle.m2 <- cbind(loc.mle.m2, mle.m2["scale"], mle.m2["shape"])
data <- mtransform(data, list(mle.m1, mle.m2))
}
# End parametric gev transform
}
}

if(rev) data <- data[,2:1]
data <- na.omit(data)
if(plot || add) x <- seq(0, 1, length = 100)
d1 <- data[,1]; d2 <- data[,2]
sum1 <- sum(d1); slm1 <- sum(log(d1))
sum2 <- sum(d2); slm2 <- sum(log(d2))
nn <- nrow(data)
nx <- length(x)

mpmin <- function(a,b) {
a[a > b] <- b[a > b]
a
}
mpmax <- function(a,b) {
a[a < b] <- b[a < b]
a
}

if(method == "cfg") {
if(!convex) {
a <- numeric(nx)
for(i in 1:nx)
a[i] <- sum(log(mpmax((1-x[i]) * d1, x[i] * d2)))
a <- (a - (1-x) * slm1 - x * slm2)/nn
a <- pmin(1, pmax(exp(a), x, 1-x))
}
else {
x2 <- seq(0, 1, length = 250)
a <- numeric(250)
for(i in 1:250)
a[i] <- sum(log(mpmax((1-x2[i]) * d1, x2[i] * d2)))
a <- (a - (1-x2) * slm1 - x2 * slm2)/nn
a <- pmin(1, pmax(exp(a), x2, 1-x2))
inch <- chull(x2, a)
a <- a[inch] ; x2 <- x2[inch]
a <- approx(x2, a, xout = x, method="linear")\$y
}
}

if(method == "pickands") {
if(!convex) {
a <- numeric(nx)
d1 <- d1/mean(d1)
d2 <- d2/mean(d2)
}
for(i in 1:nx)
a[i] <- sum(mpmin(d1/x[i], d2/(1-x[i])))
a <- a - x * sum1 - (1-x) * sum2 + nn
a <- nn / a
a <- pmin(1, pmax(a, x, 1-x))
}
else {
x2 <- seq(0, 1, length = 250)
a <- numeric(250)
d1 <- d1/mean(d1)
d2 <- d2/mean(d2)
}
for(i in 1:250)
a[i] <- sum(mpmin(d1/x2[i], d2/(1-x2[i])))
a <- a - x2 * sum1 - (1-x2) * sum2 + nn
a <- nn / a
a <- pmin(1, pmax(a, x2, 1-x2))
inch <- chull(x2, a)
a <- a[inch] ; x2 <- x2[inch]
a <- approx(x2, a, xout = x, method="linear")\$y
}
}

# Undocumented method: Tiago de Oliveira (1997)
if(method == "tdo") {
if(!convex) {
a <- numeric(nx)
for(i in 1:nx)
a[i] <- sum(mpmin(x[i]/(1 + nn*d1), (1-x[i])/(1 + nn*d2)))
a <- 1 - a/(1 + log(nn))
a <- pmin(1, pmax(a, x, 1-x))
}
else {
x2 <- seq(0, 1, length = 250)
a <- numeric(250)
for(i in 1:250)
a[i] <- sum(mpmin(x2[i]/(1 + nn*d1), (1-x2[i])/(1 + nn*d2)))
a <- 1 - a/(1 + log(nn))
a <- pmin(1, pmax(a, x2, 1-x2))
inch <- chull(x2, a)
a <- a[inch] ; x2 <- x2[inch]
a <- approx(x2, a, xout = x, method="linear")\$y
}
}

# Undocumented pot method: Beirlant et al (2004)
if(method == "pot") {
a <- numeric(nx)
rr <- rowSums(1/data)
rrk <- sort(rr, decreasing = TRUE)[k+1]
for(i in 1:nx)
a[i] <- sum(mpmax(x[i]/(d1 * rr), (1 - x[i])/(d2 * rr))[rr > rrk])
a <- 2/k * a
a0 <- 2/k * sum(1/(d2 * rr)[rr > rrk])
a1 <- 2/k * sum(1/(d1 * rr)[rr > rrk])
a <- a + 1 - (1-x) * a0 - x * a1
a <- pmin(1, pmax(a, x, 1-x))
}

plot(x, a, type="n", xlab = xlab, ylab = ylab, xlim = xlim,
ylim = ylim, ...)
polygon(c(0, 0.5, 1), c(1, 0.5, 1), lty = blty, lwd = blwd)
}
lines(x, a, lty = lty, lwd = lwd, col = col)
return(invisible(list(x = x, y = a)))
}
a
}

"qcbvnonpar"<-
function(p = seq(0.75, 0.95, 0.05), data, epmar = FALSE, nsloc1 = NULL,
nsloc2 = NULL, mint = 1, method = c("cfg","pickands","tdo"),
convex = FALSE, madj = 0, kmar = NULL, plot = FALSE, add = FALSE,
lty = 1, lwd = 1, col = 1, xlim = range(data[,1], na.rm = TRUE),
ylim = range(data[,2], na.rm = TRUE), xlab = colnames(data)[1],
ylab = colnames(data)[2], ...)
{
if(mode(p) != "numeric" || p <= 0 || p >= 1)
stop("`p' must be a vector of probabilities")
method <- match.arg(method)
nxv <- 100
x <- seq(0, 1, length = nxv)
ax <- abvnonpar(x = x, data = data, epmar = epmar, nsloc1 = nsloc1,
nsloc2 = nsloc2, method = method, convex = convex, madj = madj,
kmar = kmar, plot = FALSE)
np <- length(p)
qct <- list()
p <- p^mint
xlim <- par("usr")[1:2]
ylim <- par("usr")[1:2]
if(par("xlog")) xlim <- 10^xlim
if(par("ylog")) ylim <- 10^ylim
}
for(i in 1:np) {
qct[[i]] <- -cbind(x/ax * log(p[i]), (1-x)/ax * log(p[i]))
if(epmar) {
qct[[i]] <- cbind(quantile(data[,1], probs = exp(-qct[[i]][,1]),
na.rm = TRUE), quantile(data[,2], probs = exp(-qct[[i]][,2]),
na.rm = TRUE))
}
else {
if(is.null(kmar)) {
# Transform from exponential margins
mle.m1 <- fitted(fgev(data[,1], nsloc = nsloc1, std.err = FALSE))
mle.m2 <- fitted(fgev(data[,2], nsloc = nsloc2, std.err = FALSE))
mle.m1 <- mle.m1[c("loc","scale","shape")]
mle.m2 <- mle.m2[c("loc","scale","shape")]
qct[[i]] <- mtransform(qct[[i]], list(mle.m1, mle.m2), inv = TRUE)
}
else {
if(!is.null(nsloc1) || !is.null(nsloc2))
warning("ignoring `nsloc1' and `nsloc2' arguments")
qct[[i]] <- mtransform(qct[[i]], kmar, inv = TRUE)
}
}
qct[[i]][1,1] <- 1.5 * xlim[2]
qct[[i]][nxv,2] <- 1.5 * ylim[2]
}
if((!is.null(nsloc1) || !is.null(nsloc2)) && !epmar && is.null(kmar)) {
data <- fbvevd(data, model = "log", dep = 1, nsloc1 = nsloc1,
nsloc2 = nsloc2, std.err = FALSE)\$tdata
}
if(plot) {
plot(data, xlab = xlab, ylab = ylab, xlim = xlim, ylim = ylim, ...)
for(i in 1:np) lines(qct[[i]], lty = lty, lwd = lwd, col = col)
return(invisible(qct))
}
for(i in 1:np) lines(qct[[i]], lty = lty, lwd = lwd, col = col)
return(invisible(qct))
}
qct
}

"amvnonpar"<-
function(x = rep(1/d,d), data, d = 3, epmar = FALSE, nsloc = NULL,
madj = 0, kmar = NULL, plot = FALSE, col = heat.colors(12),
blty = 0, grid = if(blty) 150 else 50, lower = 1/3, ord = 1:3,
lab = as.character(1:3), lcex = 1)
{
if(!plot) {
if(is.vector(x)) x <- as.matrix(t(x))
if(!is.matrix(x) || ncol(x) != d)
stop("`x' must be a vector/matrix with `d' elements/columns")
if(any(x < 0, na.rm = TRUE))
stop("`x' must be non-negative")
rs <- rowSums(x)
if(any(rs <= 0, na.rm = TRUE))
stop("row(s) of `x' must have a positive sum")
if(max(abs(rs[!is.na(rs)] - 1)) > 1e-6)
warning("row(s) of `x' will be rescaled")
x <- x/rs
}
if(missing(data) || ncol(data) != d)
stop("data must have `d' columns")
if(plot) {
if(d == 2) stop("use abvnonpar for bivariate plots")
if(d >= 4) stop("cannot plot in high dimensions")
}
if(epmar) {
# Empirical transform to exponential
data <- apply(data, 2, rank, na.last = "keep")
nasm <- apply(data, 2, function(x) sum(!is.na(x)))
data <- data / rep(nasm+1, each = nrow(data))
data <- -log(data)
# End empirical transform
}
if(!epmar) {
# Parametric gev transform to exponential
if(!is.null(kmar)) {
data <- mtransform(data, kmar)
}
else {
if(!is.null(nsloc)) {
nslocmat <- list()
if(!is.list(nsloc)) nsloc <- rep(list(nsloc), d)
if(length(nsloc) != d)
stop("`nsloc' should have `d' elements")
for(k in 1:d) {
if(!is.null(nsloc[[k]])) {
if(is.vector(nsloc[[k]]))
nsloc[[k]] <- data.frame(trend = nsloc[[k]])
if(nrow(nsloc[[k]]) != nrow(data))
stop("`nsloc' and data are not compatible")
nslocmat[[k]] <- cbind(1, as.matrix(nsloc[[k]]))
}
}
}
mles <- list()
for(k in 1:d) {
mle.m <- fitted(fgev(data[,k], nsloc = nsloc[[k]], std.err = FALSE))
loc.mle.m <- mle.m[grep("^loc", names(mle.m))]
if(is.null(nsloc[[k]])) loc.mle.m <- rep(loc.mle.m, nrow(data))
else loc.mle.m <- nslocmat[[k]] %*% loc.mle.m
mles[[k]] <- cbind(loc.mle.m, mle.m["scale"], mle.m["shape"])
}
data <- mtransform(data, mles)
}
# End parametric gev transform
}

data <- na.omit(data)
{
# quicker apply(am, 2, min)
mpmin <- function(am,nr) {
a <- am[1,]
for(i in 2:nr)
a[a > am[i,]] <- am[i, a > am[i,]]
a
}
nn <- nrow(data)
nx <- nrow(x)
csum <- colSums(data)
a <- numeric(nx)
data <- nn * sweep(data, 2, csum, "/")
for(i in 1:nx) {
a[i] <- sum(mpmin(t(data)/x[i,], d))
}
a <- a - colSums(t(x) * csum) + nn
a <- nn / a
xrmax <- apply(x, 1, max)
pmin(1, pmax(a, xrmax))
}

if(plot) {
mz <- tvdepfn(depfn = depfn, col = col, blty = blty, grid = grid,
lower = lower, ord = ord, lab = lab, lcex = lcex, data = data,