context("Estimate Quantiles for given Return Periods") # importing the discharges at gauge "Wildungsmauer (Danube)" # truncating the data to complete hydrological years infile <- readlfdata("QMittelTag207373.txt", type="HZB", hyearstart = 4, baseflow = FALSE) wild <- subset(infile, hyear %in% 1996:2011) am <- tyears(wild, dist = "wei", plot = FALSE)$values # Zeitreihe mit Luecken wild.zeros <- wild wild.zeros$flow[c(10, 500)] <- 0 load("reference-gl.RData") test_that("94% quantile for Wildungsmauer is correct", { expect_equal(unname(Qxx(wild, 94)), 988.74, tolerance = 1e-2) }) pars <- tyears(wild, dist="wei", plot = FALSE)$parameters$wei test_that("return period for 94% event is the same as in script GL", { rp <- 1/cdfwei(Qxx(wild, 94), pars) # => 1.424714 (i.e.: Return period of 1.4 years) # numerisch nicht ganz ident mit tolerance = 1e-5 expect_equal(object = as.vector(rp), expected = 1.42471, tolerance = 1e-4) }) test_that("flows for given return periods are the same as in script GL", { expect_equal(as.vector(quawei(f=1/20, pars)), 668.4, tolerance = 1e-2) expect_equal(as.vector(quawei(f=1/5, pars)), 772.3, tolerance = 1e-2) }) # test_that("parameters estimated by evquantile() are correct", { # test_fit_weibull <- function(){ # shape = runif(n = 1, min = 0.5, max = 2) # scale = runif(n = 1, min = 10, max = 500) # am <- rweibull(500, shape = shape, scale = scale) # # y <- evfit(am, distribution = "wei", zeta = 0) # # return(c(shape.pop = shape, scale.pop = scale, # shape.est = y$parameters$wei["delta"], # scale.est = y$parameters$wei["beta"], # shape.e = unname((y$parameters$wei["delta"] - shape) / shape), # scale.e = unname((y$parameters$wei["beta"] - scale) / scale))) # } # # # tolerances computed with n=100 replications # error <- replicate(n = 10, expr = test_fit_weibull(), simplify = TRUE) # # expect_less_than(object = max(abs(error["shape.e", ])), # expected = 0.1) # # expect_less_than(object = max(abs(error["scale.e", ])), # expected = 0.1) # }) test_that("warnings are given", { # warn: "... using gevR instead." expect_warning(tyears(wild, dist = "gev", plot = FALSE)) # warn that zero flow observations are present expect_warning(tyears(wild.zeros, dist = "wei", plot = FALSE)) }) test_that("gevR and wei behave identically", { # only for time series with strictly positive values # pelwei complains about "invalid L-moments" y <- suppressWarnings(tyears(wild, plot = FALSE, dist = c("wei", "gevR"), event = c(1, 2, 7.9, 8, 8.1, 10))) rp <- y$T_Years_Event # quantiles for weibull und reversed GEV are identical expect_equal(rp[, "wei"], rp[, "gevR"], tolerance = 1e-10) # quantiles expect_equal(qua_ev(distribution = "gevR", f = seq(0, 1, 0.01), para = y$parameters$gevR), qua_ev(distribution = "wei", f = seq(0, 1, 0.01), para = y$parameters$wei), tolerance = 1e-10) # cdf expect_equal(cdf_ev(distribution = "gevR", x = seq(200, 1300, 1), para = y$parameters$gevR), cdf_ev(distribution = "wei", x = seq(200, 1300, 1), para = y$parameters$wei), tolerance = 1e-10) }) test_that("quantiles are INF for return period = 1 year", { # doesn't work for "gev", because a Frechet type GEV is fitted, which # has a finite upper bound y <- suppressWarnings(tyears(wild.zeros, plot = FALSE, zeta = 0, dist = c("wei", "gum", "gevR"), event = c(1))) expect_true(all(y$T_Years_Event == Inf)) y <- suppressWarnings(tyears(wild, plot = FALSE, zeta = 0, dist = c("wei", "gum", "gevR"), event = c(1))) expect_true(all(y$T_Years_Event == Inf)) }) test_that("quantiles for mixed distributions are plausible", { y <- suppressWarnings(tyears(wild.zeros, plot = FALSE, zeta = 0, dist = c("wei", "gum", "gevR"), event = c(1, 2, 7.9, 8, 8.1, 10, 50))) rp <- y$T_Years_Event # return period of 1 has infinite quantiles expect_true(all(rp[rownames(rp) == 1, ] == Inf)) # return periods lower than frequency of zero obersvations yield >0m³/s mask <- as.numeric(rownames(rp)) < 1 / y$freq.zeros expect_true(all(rp[mask, ] > 0)) # return periods higher than frequency of zero obersvations yield 0m³/s mask <- as.numeric(rownames(rp)) >= 1 / y$freq.zeros expect_true(all(rp[mask, ] == 0)) }) # test_that("evquantiles() gives the same results as WMO manual", { # # wmo <- list( # nicholson = list(data = c(8, 13, 2, 2.7, 9.1, 0, 0.7, 21.6, 35.1, 12.7, 1.3, # 9.1, 26.0, 22.7, 21.3, 7.9, 23, 6, 3.1, 4.1, 0, # 4.9, 2.3, 12.4, 4.4, 2.1, 3.1, 7.6, 9, 9, 20.4, # 13.6, 3.3, 10.3), # lmom = matrix(c(9.765, 10.375, 4.721, 4.670, 0.273, 0.276), # ncol = 3, nrow = 2, # dimnames = list(c("values", "censored"), # c("l_1", "l_2", "t_3"))), # pmom = c(9.765, 8.68, 1.149), # quantile = 1.2), # # timbarra = list(data = c(27.3, 38, 72.9, 53.6, 45.7, 37, 21.7, 36.6, 43.1, # 12.6, 20.7, 66.7, 76.6, 53.1, 15.7, 66.6, 78, 53.1, # 62.6, 40.7, 51.6, 23.1, 23.4, 23, 6), # lmom = matrix(c(41.977, 12.305, 0.045), # ncol = 3, nrow = 1, # dimnames = list(c("values"), # c("l_1", "l_2", "t_3"))), # pmom = c(41.977, 21.05, 0.146), # quantile = 15.4)) # # # # # # reproduce product moments (to test if input data is correct) # library(e1071) # pmom <- sapply(c(mean, sd, skewness), function(x) x(wmo$nicholson$data)) # expect_true(abs(pmom[1] - wmo$nicholson$pmom[1]) < 1e-3) # expect_true(abs(pmom[2] - wmo$nicholson$pmom[2]) < 1e-2) # expect_true(abs(pmom[3] - wmo$nicholson$pmom[3]) < 1e-3) # # pmom <- sapply(c(mean, sd, skewness), function(x) x(wmo$timbarra$data)) # expect_true(abs(pmom[1] - wmo$timbarra$pmom[1]) < 1e-3) # expect_true(abs(pmom[2] - wmo$timbarra$pmom[2]) < 1e-2) # expect_true(abs(pmom[3] - wmo$timbarra$pmom[3]) < 1e-3) # # # reproduce L-moments # nicholson <- evfit(x = wmo$nicholson$data, # distribution = "wei", zeta = 0) # # timbarra <- evfit(x = wmo$timbarra$data, # distribution = "wei", zeta = 0) # # expect_true(abs(nicholson$lmom[, c("l_1", "l_2", "t_3")] - # wmo$nicholson$lmom[, c("l_1", "l_2", "t_3")]) < 1e-3) # # expect_true(abs(timbarra$lmom["raw data", c("l_1", "l_2", "t_3")] - # wmo$timbarra$lmom[, c("l_1", "l_2", "t_3")]) < 1e-3) # # # reproduce parameters of distribution # # # # reproduce equantile of 10 year event # quant <- evquantile(nicholson, return.period = 10) # expect_equal(object = as.vector(quant$T_Years_Event), # expected = wmo$nicholson$quantile, # tolerance = 1e-2) # # quant <- evquantile(timbarra, return.period = 10) # expect_equal(object = as.vector(quant$T_Years_Event), # expected = wmo$timbarra$quantile, # tolerance = 1e-2) # }) test_that("same results as skripts from GL", { # all gl fits with pooling are incorrect, because of wrong pooling in old package # tyears_new_GL lfstat <- suppressWarnings(tyears(wild, dist = c("gevR", "wei"), plot = FALSE)) expect_equal2(object = lfstat$parameter[["wei"]], expected = reference.gl$tyears_new_GL$parameters[["wei"]]) expect_equal2(object = lfstat$T_Years_Event[1, "wei"], expected = unname(reference.gl$tyears_new_GL$T_Years_Event["wei"])) expect_equal2(object = lfstat$parameter[["gevR"]], expected = reference.gl$tyears_new_GL$parameters[["gev"]]) expect_equal2(object = lfstat$T_Years_Event[1, "gevR"], expected = unname(reference.gl$tyears_new_GL$T_Years_Event["gev"])) # tyears_MAX_D_hyear_GL-sum lfstat <- suppressWarnings(tyearsS(wild, dist = c("gev", "wei"), variable = "d", aggr = sum, plot = FALSE, threshold = function(x) quantile(x, probs = 0.06, na.rm = TRUE))) expect_equal2(object = lfstat$parameter[["gev"]], expected = reference.gl$"tyears_MAX_D_hyear_GL-sum"$parameters[["gev"]]) expect_equal2(object = lfstat$T_Years_Event[1, "gev"], expected = unname(reference.gl$"tyears_MAX_D_hyear_GL-sum"$T_Years_Event["gev"])) # tyears_MAX_D_hyear_GL-max # lfstat <- suppressWarnings(tyearsS(wild, dist = c("gev", "wei"), # pooling = pool_ic, # variable = "d", aggr = max, # plot = FALSE, # threshold = function(x) quantile(x, probs = 0.06, na.rm = TRUE))) # # expect_equal2(object = lfstat$parameter[["gev"]], # expected = reference.gl$"tyears_MAX_D_hyear_GL-max"$parameters[["gev"]]) # expect_equal2(object = lfstat$T_Years_Event[1, "gev"], # expected = unname(reference.gl$"tyears_MAX_D_hyear_GL-max"$T_Years_Event["gev"])) # # tyears_MAX_V_hyear_GL-sum lfstat <- suppressWarnings(tyearsS(wild, dist = c("gev", "wei"), variable = "v", threshold = function(x) quantile(x, probs = 0.06, na.rm = TRUE), aggr = sum, plot = FALSE)) expect_equal2(object = lfstat$parameter[["gev"]], expected = reference.gl$"tyears_MAX_V_hyear_GL-sum"$parameters[["gev"]], tolerance = 1e-4) expect_equal2(object = lfstat$T_Years_Event[1, "gev"], expected = unname(reference.gl$"tyears_MAX_V_hyear_GL-sum"$T_Years_Event["gev"]), tolerance = 1e-4) # tyears_MAX_V_hyear_GL-max # lfstat <- suppressWarnings(tyearsS(wild, dist = c("gev", "wei"), variable = "v", # threshold = function(x) quantile(x, probs = 0.06, na.rm = TRUE), # aggr = max, pooling = pool_ic, plot = FALSE)) # # expect_equal2(object = lfstat$parameter[["gev"]], # expected = reference.gl$"tyears_MAX_V_hyear_GL-max"$parameters[["gev"]], # tolerance = 1e-7) # expect_equal2(object = lfstat$T_Years_Event[1, "gev"], # expected = unname(reference.gl$"tyears_MAX_V_hyear_GL-max"$T_Years_Event["gev"]), # tolerance = 1e-6) })