tests/testthat/test_amax.R

In martindurocher/floodStat: FloodNet Regional Frequency Analysis

rm(list = ls())
######################################################
## testing FitAmax
#######################################################

uu <- (1:100 - .5) /100
xgev <- qAmax(uu, c(100,30,0), 'gev')
vname <-  c("lmom", "method", "para", "distr", "varcov", "llik", "data")

## Verify that important distribution and method works
TestObj <- function(f){
  expect_is(f, 'amax')
  expect_equal(names(f), vname)
  expect_true(is.na(f$varcov))
  expect_equal(length(f$data), 100)
  expect_equal(length(coef(f)), 3)
  expect_true(all(is.finite(coef(f))))
}

S <- function(z) signif(as.numeric(z))

TestObj(f <- FitAmax(xgev, distr = 'gev', method = 'lmom', varcov = FALSE))
expect_equal(S(f$para), c(99.8386, 30.1261,  0.0000))
expect_equal(S(f$para), c(99.8386, 30.1261,  0.0000))

TestObj(f <- FitAmax(xgev, distr = 'gno', method = 'lmom', varcov = FALSE))
expect_equal(S(f$lmom), c(117.228, 20.8818, 3.54844, 3.19116, 1.18389))

TestObj(FitAmax(xgev, distr = 'glo', method = 'lmom', varcov = FALSE))
TestObj(FitAmax(xgev, distr = 'pe3', method = 'lmom', varcov = FALSE))
TestObj(FitAmax(xgev, distr = 'gev', method = 'mle', varcov = FALSE))
TestObj(FitAmax(xgev, distr = 'gno', method = 'mle', varcov = FALSE))
TestObj(FitAmax(xgev, distr = 'glo', method = 'mle', varcov = FALSE))

TestObj(f <- FitAmax(xgev, distr = 'pe3', method = 'mle', varcov = FALSE))
expect_match(f$method , 'mle')
expect_match(f$distr , 'pe3')

## verify output quantities are fine
f <- FitAmax(xgev,'gev', method = 'lmom', varcov = TRUE, nsim = 10)
expect_true(all(is.finite(vcov(f))))
expect_equal(dim(vcov(f)), c(3,3))
expect_equal(S(AIC(f)), 1000.42)
expect_true(all(is.finite(coef(f))))


f <- FitAmax(xgev,'gev', method = 'mle', varcov = TRUE)
expect_true(all(is.finite(vcov(f))))
expect_equal(dim(vcov(f)), c(3,3))
expect_equal(S(AIC(f)), 1000.39)


## Verify the automatic selection
f <- FitAmax.auto(xgev, distr = c('gev','glo','gno','pe3'), method = 'mle')
expect_match(f$distr, 'gno')

f <- FitAmax.auto(xgev, distr = c('gev','glo','gno','pe3'),
                  method = 'mle', tol.gev = 2)
expect_match(f$distr, 'gev')

f <- FitAmax.auto(xgev, distr = c('gev','gum'), method = 'mle')
expect_match(f$distr, 'gum')

f <- FitAmax.auto(xgev, distr = c('gev','glo','gno','pe3'),
                  method = 'lmom', varcov = FALSE)
expect_match(f$distr, 'gno')

f <- FitAmax.auto(xgev, distr = c('gev','glo','gno','pe3'),
                  method = 'lmom', varcov = FALSE, tol.gev = 2)
expect_match(f$distr, 'gev')

f <- FitAmax.auto(xgev, distr = c('gev','gum'), method = 'lmom',
                  varcov = FALSE)
expect_match(f$distr, 'gum')


## Verifying return error
expect_error(FitAmax(xgev, distr = 'gev', nsim = 1))

xgev[1] <- NA
expect_error(FitAmax(xgev, distr = 'gev'))

xgev[1] <- Inf
expect_error(FitAmax(xgev, distr = 'gev'))


######################################################
## testing FitGev
#######################################################

uu <- (seq(1000)-.5)/1000
xgev <- qAmax(uu, c(100,30, -.15), 'gev')
f <- FitGev(xgev)

expect_equal(S(f$para), c(99.9884, 29.9728, -0.149058))
expect_match(f$distr, 'gev')
expect_match(f$method, 'gml')
expect_equal(f$prior, c(-0.1, 0.015))
expect_true(all(is.finite(vcov(f))))
expect_equal(length(f$data), 1000)
expect_equal(S(AIC(f)), 10134.6)


## as approximation to mle when uniform prior is used
f <- FitGev(xgev, varcov = FALSE, mu = 0, sig2 = 1/12)
f0 <- FitAmax(xgev, 'gev', method = 'mle', varcov = FALSE)
expect_equal(S(AIC(f0)), S(AIC(f)))
expect_true(all((1-coef(f)/coef(f0))^2 < 1e-6))


######################################################
## testing predict.amax function
#######################################################
x <- ExtractAmax(flow~date,flowStJohn, tol = 355)

fit <- FitAmax(x$flow,'gev', method = 'mle')

rp <- 1-1/c(10,100)
rlev <- predict(fit, q = rp)
expect_equal(signif(rlev), c(3356.64, 4250.40))
expect_is(rlev, 'numeric')

out <- predict(fit, se = TRUE, ci = 'delta')
expect_true(all(colnames(out) == c('pred','se','lower','upper')))
expect_is(out, 'data.frame')

out <- predict(fit, se = FALSE, ci = 'delta')
expect_true(all(colnames(out)== c('pred','lower','upper')))

out <- predict(fit, se = TRUE)
expect_true(all(colnames(out) == c('pred','se')))


## The bootstrap sample used for CI are returned
fit <- FitAmax(x$flow, distr = 'gev', varcov = FALSE)
boot <- predict(fit, rp, se = FALSE, ci = 'boot',
                 nsim = 5, out.matrix = TRUE)

expect_is(boot, 'list')
expect_true(all(names(boot) == c('pred','para','qua')))


######################################################
## testing fAmax functions
#######################################################

x <- rAmax(10, c(0,1,0), 'gno')
expect_equal(length(x), 10)

x <- rAmax(10, c(0,1,0), 'gev')
expect_equal(length(x), 10)

q <- qAmax(uu, c(100,30,-.5), 'gev')
p <- pAmax(q, c(100,30,-.5), 'gev')
expect_equal(uu,p)

q <- qAmax(uu, c(100,30,-.5), 'gno')
p <- pAmax(q, c(100,30,-.5), 'gno')
expect_equal(uu,p)

AIC0 <- -2*sum(dAmax(f0$data, f0$para, f0$distr, log = TRUE)) + 6
expect_equal(AIC0, AIC(f0))

f0 <- FitAmax(xgev, 'gum', method = 'lmom', varcov = FALSE)
AIC0 <- -2*sum(dAmax(f0$data, f0$para, f0$distr, log = TRUE)) + 4
expect_equal(AIC0, AIC(f0))