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

set.seed(20101111)

# large sample estimation, testing Normality of estimators with and without covariates in model
myTest <- function(Mu, Phi, Xi, ..., info){
  ests <- lapply(1:nreps, function(i) evm(data=data.frame(y=x[,i],x=covariate), y=y,family=glo,...))
  par <- sapply(1:nreps,function(i)ests[[i]]$par)
  se <- sapply(1:nreps,function(i)ests[[i]]$se)
  true <- cbind(Mu,Phi,Xi)
  p.vals <- apply((t(par)-true)/t(se),2,function(x)shapiro.test(x)$p.value)

  expect_true(min(p.vals) >= alpha, info=info)
}

#*************************************************************
# Set up simulation test info

ntests <- 21 #total number of parameters estimated in tests below:
nreps <- 20
nsim <- 10000
alpha <- 0.05/ntests # bonferroni
p <- matrix(runif(6*nreps, -1, 1),ncol=6)
p[,3] <- p[,3]
p[,5] <- p[,5]/3
p[,6] <- p[,6]/4
covariate <- seq(0.1,0.9,length=nsim)

#*************************************************************
# Test stationary model

x <- sapply(1:nreps,function(i)rglo(nsim,mu=p[i,1],sigma=exp(p[i,3]), xi=p[i,5]))

myTest(Mu=p[,1], Phi=p[,3], Xi=p[,5], info="glo: no covariates, p-value")

#*************************************************************
# Test covariate in mu only

x <- sapply(1:nreps,function(i)rglo(nsim,mu=p[i,1] + p[i,2]*covariate,sigma=exp(p[i,3]), xi=p[i,5]))

myTest(Mu=p[,1:2],Phi=p[,3], Xi=p[,5], mu=~x, info="glo: covariates in mu only, p-value")

#*************************************************************
# Test covariate in sigma only

x <- sapply(1:nreps,function(i)rglo(nsim,mu=p[i,1],sigma=exp(p[i,3]+p[i,4]*covariate), xi=p[i,5]))

myTest(Mu=p[,1],Phi=p[,3:4], Xi=p[,5], phi=~x,info="glo: covariates in phi only, p-value")

#*************************************************************
# Test covariate in xi only

x <- sapply(1:nreps,function(i)rglo(nsim,mu=p[i,1],sigma=exp(p[i,3]), xi=p[i,5]+p[i,6]*covariate))

myTest(Mu=p[,1],Phi=p[,3], Xi=p[,5:6], xi=~x, info="glo: covariates in xi only, p-value")

#*************************************************************
# Test covariate in all variables

x <- sapply(1:nreps,function(i)rglo(nsim,mu=p[i,1]+ p[i,2]*covariate,sigma=exp(p[i,3]+p[i,4]*covariate), xi=p[i,5]+p[i,6]*covariate))

myTest(Mu=p[,1:2],Phi=p[,3:4], Xi=p[,5:6], mu=~x,phi=~x,xi=~x, info="glo: covariates in all parameters, p-value")


##############################################################################
## pglo
set.seed(20101111)

probabilities <- runif(10)

xi.values <- c(0, seq(-0.3, 0.3, length.out=10))

core.sanity.test <- function(xi) {
  randoms <- sort(rglo(10, 0, 1, xi))
  expect_true(all(diff(pglo(randoms, 0, 1, xi)) >= 0),
              info="pglo: ascending")
  expect_true(all(diff(pglo(randoms, 0,1,xi,lower.tail=FALSE)) <= 0),
              info="pglo: descending")

  expect_equal(log(pglo(randoms, 0, 1,xi)),
               pglo(randoms, 0, 1, xi, log.p=TRUE),
               info="pglo: log.p (1)")
  expect_equal(log(pglo(randoms, 0, 1, xi, lower.tail=FALSE)),
               pglo(randoms, 0, 1, xi, lower.tail=FALSE, log.p=TRUE),
               info="pglo: log.p (2)")

  mu <- runif(1, -5, 5)
  sigma <- rexp(1)

  expect_equal(pglo(randoms, 0, 1,xi),
               pglo(mu + sigma * randoms, mu, sigma, xi),
               info="pglo: shift and scale")
} # Close core.sanity.tests

qglo.comparison <- function(xi) {
  quantiles <- qglo(probabilities, 0, 1, xi)
  my.probs  <- pglo(quantiles, 0, 1, xi)
  expect_equal(my.probs, probabilities, info="pglo:straighttest")

  quantiles <- qglo(probabilities, 0, 1, xi, lower.tail=FALSE)
  my.probs  <- pglo(quantiles, 0, 1, xi, lower.tail=FALSE)
  expect_equal(my.probs, probabilities, info="pglo:lowertail")

  my.probs.2 <- pglo(quantiles, 0, 1, xi, lower.tail=TRUE)
  expect_equal(probabilities+my.probs.2,
               rep(1, length(probabilities)),
               info="pglo: tail flip")
} # Close qglo.comparison

lapply(xi.values, core.sanity.test)
##lapply(xi.values, qglo.comparison)

##############################################################################
## qglo

## get the probabilities that we'll use and sort them
## into ascending order for safekeeping
set.seed(20101111)
probabilities <- sort(runif(10))

core.sanity.test <- function(xi) {
  base.quantiles <- qglo(probabilities, 0, 1, xi)
  ## check that the values are ascending
  expect_true(all(diff(base.quantiles)>=0), info="qglo:ascendingquantiles")
  ## and check that we're descending correctly for upper tail
  bq2 <- qglo(probabilities, 0, 1, xi, lower.tail=FALSE)
  expect_true(all(diff(bq2)<=0), info="qglo:descendingquantiles")
  ## does lower.tail work
  expect_equal(base.quantiles,
               qglo(1 - probabilities, 0, 1, xi, lower.tail=FALSE),
               info="qglo: lower.tail works correctly")
  ## does log.p work?
  expect_equal(base.quantiles,
               qglo(log(probabilities), 0, 1, xi, log.p=TRUE),
               info="qglo: log.p works")
  ## check shift and scale property
  sigma <- rexp(1)
  mu    <- runif(1, -5, 5)
  shifted <- mu + sigma * base.quantiles
  expect_equal(shifted, qglo(probabilities, mu, sigma, xi),
               info="qglo: shift and scale")
} # Close core.sanity.test

##lapply(c(0, seq(-5, 5, length.out=10)), core.sanity.test)

## known values
expect_equal(0, qglo(0.5, 0, 1, 0),
             info="qglo: median match at zero xi")
xi <- seq(-2, 2, length=10)
expect_equal(qglo(0.5, 0, 1, xi),
             rep(0,10),
             info="qglo: median match at nonzero xi")

##############################################################################
## dglo

local.dglo <- function(x,loc,scale,shape){ # function is not vectorised
  if(shape == 0){
    Ex <- exp(-(x-loc)/scale)
    out <- Ex/(scale* (1+Ex)^2)
  } else {
    Ox <- 1+shape/scale * (x-loc)
    out <- Ox^{-1/shape - 1} / (scale * (1+ Ox^(-1/shape))^2)
  }
  out
}

myTest <- function(mu, sig, xi, info){
  myd <- sapply(1:nreps, function(i) dglo(x[,i], mu[i],sig[i], xi[i]))
  ed <- sapply(1:nreps, function(i) local.dglo(x[,i], loc=mu[i], scale=sig[i], shape=xi[i]))
  expect_equal(ed, myd, info=info)
}

set.seed(20101111)

#*************************************************************
# Test dglo. Note that local.dglo is NOT vectorized.

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

x <- sapply(1:nreps,
            function(i)rglo(nsim,mu=mu[i],sigma=p[i,1], xi=p[i,2]))

myTest(mu=mu,sig=p[,1], xi=p[,2], info="dglo: random xi")

#*************************************************************
# Test dglo when some or all of xi == 0

p[sample(1:nreps,nreps/2),2] <- 0
x <- sapply(1:nreps,function(i)rglo(nsim,mu=mu[i],sigma=p[i,1],xi=p[i,2]))
myTest(mu=mu,sig=p[,1], xi=p[,2], info="dglo: some zero xi")

p[,2] <-  0
x <- sapply(1:nreps,function(i)rglo(nsim,mu=mu[i],sigma=p[i,1],xi=p[i,2]))
myTest(mu=mu,sig=p[,1], xi=p[,2], info="dglo: all zero xi")

#*************************************************************
# Test vectorization of dglo.

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

x <- rglo(nsim, mu, sig, xi)
myd <- dglo(x, mu, sig, xi)

ed <- sapply(1:nsim, function(i) local.dglo(x[i], loc=mu[i], scale=sig[i], shape=xi[i]))
expect_equal(ed, myd, info="dglo:vectorisation")

#*************************************************************
# test log.d argument

ld <- dglo(x,mu,sig,xi,log.d=TRUE)
expect_equal(myd, exp(ld), info="dglo:logdensity")

##############################################################################
## rglo
## so, how do we test an RNG...
num.simple <- 1000
num.quantile <- 1e6

xi.values  <- c(0, seq(-5, 5, length.out=10))
test.quantiles <- c(0.25, 0.5, 0.75)

core.sanity.test <- function(xi) {
  seed <- as.integer(runif(1, -1, 1)*(2**30))
  set.seed(seed)
  samples <- rglo(num.simple, 0, 1, xi)
  expect_equal(length(samples), num.simple,
               info = "rglo: output of correct length")
  if (xi > 0) {
    expect_true(all(samples>=-1/xi), info = "rglo:lowerboundcheck")
  } else if (xi < 0) {
    expect_true(all(samples<=-1/xi), info = "rglo:upperboundcheck")
  }
  ## scale and shift property
  sigma <- rexp(1)
  mu    <- runif(1, -5, 5)
  shifted <- mu + sigma * samples
  set.seed(seed)
  expect_equal(shifted, rglo(num.simple, mu, sigma, xi), info = "rglo: scale and shift")
} # Close core.sanity.test

quantile.test <- function(xi) {
  ## here are the sampled quantiles
  quantiles <- quantile(pglo(rglo(num.quantile, 0, 1, xi),
                             0, 1, xi),
                        probs=test.quantiles,
                        names=FALSE,na.rm=TRUE)
  ## this is a bit crude, but hey...
  expect_equal(test.quantiles, quantiles, tolerance=0.02, info="rglo: quantile test")
} # Close quantile.test
lapply(xi.values, core.sanity.test)
lapply(xi.values, quantile.test)