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

set.seed(201110101)
evd.qgpd <- texmex:::.evd.qgpd
myTest <- function(sig, xi, thresh, msg){
  myq <- sapply(1:nreps,function(i) qgpd(x[,i], sig[i], xi[i], u=thresh[i]))
  myp <- sapply(1:nreps,function(i) pgpd(myq[,i], sig[i], xi[i], u=thresh[i]))
  eq <- sapply(1:nreps, function(i) evd.qgpd(x[,i], loc=thresh[i], scale=sig[i], shape=xi[i]))
  expect_equal(eq, myq, info=paste(msg,"testusing.evd.qgpd"))
  expect_equal(x, myp, info=paste(msg,"testusingqgpd"))
}

#*************************************************************
# 6.4.0 Test exception for out of range probabilties
expect_error(qgpd(1.5, 1, 0, 2), info="qgpd:exceptionforoutofrangeprob")
expect_error(qgpd(-1, 1, 0, 2), info="qgpd:exceptionforoutofrangeprob")

#*************************************************************
# 6.4. Test qgpd. Note that .evd.qgpd is NOT vectorized.

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

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

#*************************************************************
# 6.5. Test qgpd when some or all of xi == 0. Note that .evd.qgpd is NOT vectorized.

p[sample(1:nreps,nreps/2),2] <- 0
myTest(sig=p[,1], xi = p[,2], thresh=thresh, msg="qgpd: some zero xi")
p[,2] <-  0
myTest(sig=p[,1], xi = p[,2], thresh=thresh, msg="qgpd: all zero xi")

#*************************************************************
# 6.6. Test vectorization of qgpd. Note that .evd.qgpd is NOT vectorized.

sig <- runif(nsim, 0, 2)
xi <- runif(nsim)
thresh <- rnorm(nsim)

x <- runif(nsim)

myq <- qgpd(x, sig, xi, thresh)
eq <- sapply(1:nsim, function(i) evd.qgpd(x[i], loc=thresh[i], scale=sig[i], shape=xi[i]))

expect_equal(eq, myq, info="qgpd:vectorisation")

#*************************************************************
# 6.6a Test log.p argument

lq <- qgpd(log(x), sig,xi,thresh,log.p=TRUE)

expect_equal(myq, lq, info="qgpd:log.p=TRUE")

#*************************************************************
# 6.6a Test log.p argument

LTq <- qgpd(1-x, sig,xi,thresh, lower.tail=FALSE)

expect_equal(myq, LTq, info="qgpd:lower.tail=FALSE")