test-egp3.R
In texmex: Statistical Modelling of Extreme Values

## texmex erroneously does not export qegp3 (2022-05-10)
qegp3 <- texmex:::qegp3

library(MASS)
rmod <- rlm(log(ALT.M) ~ log(ALT.B) + as.numeric(dose), data=liver, method="MM", c=3.44)
liver$r <- resid(rmod)

# Check that GPD and EGP3 match when log(k) is effectively constrained to 0
gpmod <- evm(liver[liver$dose == "D", "r"], qu=.6)
pp <- list(c(0, 0, 0), diag(c(10^(-10), 10^4, 10^4)))
e3mod <- evm(liver[liver$dose == "D", "r"], qu=.6, family=egp3, priorParameters=pp)

# Test the point estimates, standard errors and t-values all at once
expect_equal(coef(summary(gpmod)), coef(summary(e3mod))[-1, ], tolerance = .0001, label="egp3: matches gpd")

# Check SEs on return levels - derivatives were worked out manually (by Paul),
# also by Sage and (I think) by Wolfram
rgp <- do.call("rbind", predict(gpmod, M=seq(100, 2000, by=100), se.fit=TRUE)$obj)
re3 <- do.call("rbind", predict(e3mod,  M=seq(100, 2000, by=100), se.fit=TRUE)$obj)
expect_equal(rgp, re3, tolerance=.0001, label="egp3: return levels")

# Check that the plot for the Nidd data has generally similar shape to that in
# Papastathopoulos and Tawn: they use EGP1, not 3, so the plots won't be the same
plot(egp3RangeFit(nidd, umin=65, umax=90, nint=50), pch=16)

## regp3 behaves as it should
regp3_test <- function(n, kappa=1, sigma, xi, u=0){
  # Implement rng via direct cdf inversion and compare to version that uses qegp3
  kappa <- rep(kappa, length.out=n)
  sigma <- rep(sigma, length.out=n)
  xi <- rep(xi, length.out=n)
  u <- rep(u, length.out=n)

  nk1 <- sum(kappa != 1)

  x <- 1 - runif(nk1)^(1/kappa)

  xi.neq.0 <- function(x, xi, sigma, u) (x^(-xi) - 1) * sigma / xi + u
  xi.eq.0 <- function(x, xi, sigma, u) log(x) * sigma / xi + u

  res <- ifelse(xi == 0, xi.eq.0(x, xi, sigma, u), xi.neq.0(x, xi, sigma, u))

  if (nk1 < n){
    res2 <- rgpd(n - nk1, sigma, xi, u)
    r <- rep(0, n)
    r[kappa != 1] <- res
    r[kappa == 1] <- res2
    r
  } else res
}

par(mfrow=c(3, 3), oma=c(0, 0, 2, 0))
for (i in 1:9){
  kappa <- runif(1); sigma <- runif(1, .1, 6); xi <- runif(1); u <- runif(1, -10, 10)
  x1 <- sort(regp3_test(1000, kappa=.5, sigma=2, xi=.2, u=6))
  x2 <- sort(regp3(1000, kappa=.5, sigma=2, xi=.2, u=6))
  expect_true(cor(x1, x2) > 0.85)
  plot(x1, x2); abline(0, 1)
}
title("Comparing 2 implementations of EGP3 rng", outer=TRUE)


## qegp3 behaves as expected
# Note that regp3 calls on qegp3, so some implicit testing is performed by the
# regp3 tests. Also, qegp3 calls on qgpd, so again some implicit testing is done
# there.

# Test for error with out of range probabiliites
expect_error(qegp3(-.1, xi=.2, sigma=1))
expect_error(qegp3(1.1, xi=.2, sigma=1))

# Check boundaries
expect_equal(qegp3(0, xi=.2, sigma=1, kappa=runif(1)), 0)
expect_equal(qegp3(1, xi=.2, sigma=1, kappa=runif(1)), Inf)

# Check that qegp3 and pegp3 reverse each other
myTest <- function(sig, xi, kappa, thresh, msg){
  myq <- sapply(1:nreps,function(i) qegp3(x[,i], sig[i], xi[i], kappa[i], u=thresh[i]))
  myp <- sapply(1:nreps,function(i) pegp3(myq[,i], sig[i], xi[i], kappa[i], u=thresh[i]))
  expect_equal(x, myp, label=paste(msg,"testusingqgpd"))
}

nreps <- 100
nsim <- 1000
p <- matrix(runif(2*nreps, -1, 1), ncol=2)
p <- cbind(p, runif(nreps)*3)
p[, 1] <- p[, 1] + 1
thresh <- rep(0,nreps)
x <- matrix(runif(nreps*nsim), nrow=nsim)

myTest(sig=p[,1], xi=p[,2], kappa=p[,3], thresh=thresh, msg="qgpd: random xi")

# Check that qegp3 gives same answer on log scale
p <- qegp3(0.95, kappa=0.5, sigma=1, xi=.01)
p2 <- qegp3(log(0.95), kappa=0.5, sigma=1, xi=.01, log.p=TRUE)
expect_equal(p, p2, label="qegp3 returns same answer on log and untransformed scale, lower tail")

p <- qegp3(0.95, kappa=0.5, sigma=1, xi=.01, lower.tail=FALSE)
p2 <- qegp3(log(0.95), kappa=0.5, sigma=1, xi=.01, log.p=TRUE, lower.tail=FALSE)
expect_equal(p, p2, label="qegp3 returns same answer on log and untransformed scale, upper tail")


## pegp3 behaves as expected
# pegp3 does some /fairly/ simple stuff and calls on pgpd, so the testing of
# pgpd implicitly does some testing of pegp3. Also the test of qegp3 implicity
# tests pegp3

for (i in 1:10){
  s <- runif(100) + 1
  xi <- runif(100, -1, 1)
  k <- runif(100) * 3
  p <- runif(100)

  # Test lower.tail argument
  x <- pegp3(p, sigma=s, xi=xi, kappa=k)
  y <- 1 - pegp3(p, sigma=s, xi=xi, kappa=k, lower.tail=FALSE)
  expect_equal(x, y)

  # Test wtih log.p
  x <- exp(pegp3(p, sigma=s, xi=xi, kappa=k, log.p=TRUE))
  y <- 1 - exp(pegp3(p, sigma=s, xi=xi, kappa=k, lower.tail=FALSE, log.p=TRUE))
  expect_equal(x, y)
}

## degp3 behaves as expected
degp3_test <- function(x, kappa=1, sigma, xi){
  y <- 1 + x * xi/sigma
  res <- log(kappa/sigma) + (kappa - 1) * log(1 - y^(-1/xi)) - (1/xi + 1) * log(y)
  exp(res)
}

for (i in 1:10){
  xi <- runif(1, -0.4999, .9999)
  kappa <- runif(1) * 3
  sigma <- runif(1) + 1

  x <- regp3(100, sigma=sigma, xi=xi, kappa=kappa)
  d1 <- degp3(x, sigma=sigma, xi=xi, kappa=kappa)
  d2 <- degp3_test(x, sigma=sigma, xi=xi, kappa=kappa)

  # Clunky workaround: some combinations generate numbers absurdly small or large
  # that degp3_test evaluates to +/- Inf
  d1 <- d2[d2 > 1e-20 & d2 < 1e20]
  d2 <- d2[d2 > 1e-20 & d2 < 1e20]
  expect_equal(d1, d2, tolerance=1.0e-02)
}