Nothing
## 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)
}
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.