From 4dca06cc63c13750302bd14358e389e058a28c61 Mon Sep 17 00:00:00 2001 From: phil Date: Tue, 16 Nov 2021 11:23:42 +0100 Subject: [PATCH 1/4] remove deseasonalize dependency --- DESCRIPTION | 1 - 1 file changed, 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5593cd7..2b547ba 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -18,7 +18,6 @@ Imports: grid, Rcpp, RColorBrewer, - deseasonalize, Lmoments, moments, numDeriv, From 79996d1efd9f1664d0439809332fa2e71c582a98 Mon Sep 17 00:00:00 2001 From: theGreatWhiteShark Date: Mon, 11 Jul 2022 11:17:49 +0200 Subject: [PATCH 2/4] Fix warnings in try-error checks Some objects have more than a single class / return a string list when provided to the `class()` function. All checks whether the class of an error is "try-error" is now wrapped in `any()` to ensure it both works and does not produce any warnings --- R/extremes.R | 16 ++++++++-------- R/opti.R | 28 ++++++++++++++-------------- R/utils.R | 6 +++--- 3 files changed, 25 insertions(+), 25 deletions(-) diff --git a/R/extremes.R b/R/extremes.R index e5b19e6..26eaf23 100644 --- a/R/extremes.R +++ b/R/extremes.R @@ -1371,14 +1371,14 @@ return.level.climex.fit.gev <- x$control$hessian[ 1 : 2, 1 : 2 ] ), silent = silent ) ## Augment the result again to ensure compatibility - if ( class( error.covariance ) != "try-error" ){ + if ( any( class( error.covariance ) != "try-error" ) ){ dummy.matrix <- matrix( rep( 0, 9 ), nrow = 3, ncol = 3 ) dummy.matrix[ 1 : 2, 1 : 2 ] <- error.covariance error.covariance <- dummy.matrix } } - if ( class( error.covariance ) == "try-error" ){ + if ( any( class( error.covariance ) == "try-error" ) ){ x.hessian <- numDeriv::hessian( likelihood, x = x$par, x.in = x$x, model = "gev" ) @@ -1754,14 +1754,14 @@ return.level.climex.fit.gpd <- x$control$hessian[ 1 ] ), silent = silent ) ## Augment the result again to ensure compatibility - if ( class( error.covariance ) != "try-error" ){ + if ( any( class( error.covariance ) != "try-error" ) ){ dummy.matrix <- matrix( rep( 0, 4 ), nrow = 2, ncol = 2 ) dummy.matrix[ 1 ] <- error.covariance error.covariance <- dummy.matrix } } - if ( class( error.covariance ) == "try-error" ){ + if ( any( class( error.covariance ) == "try-error" ) ){ x.hessian <- numDeriv::hessian( likelihood, x = x$par, x.in = x$x, model = "gpd" ) @@ -2938,14 +2938,14 @@ upper.limit.climex.fit.gev <- x$control$hessian[ 1 : 2, 1 : 2 ] ), silent = TRUE ) ## Augment the result again to ensure compatibility - if ( class( error.covariance ) != "try-error" ){ + if ( any( class( error.covariance ) != "try-error" ) ){ dummy.matrix <- matrix( rep( 0, 9 ), nrow = 3, ncol = 3 ) dummy.matrix[ 1 : 2, 1 : 2 ] <- error.covariance error.covariance <- dummy.matrix } } - if ( class( error.covariance ) == "try-error" ){ + if ( any( class( error.covariance ) == "try-error" ) ){ x.hessian <- numDeriv::hessian( likelihood, x = x$par, x.in = x$x, model = "gev" ) @@ -3269,14 +3269,14 @@ upper.limit.climex.fit.gpd <- x$control$hessian[ 1 ] ), silent = TRUE ) ## Augment the result again to ensure compatibility - if ( class( error.covariance ) != "try-error" ){ + if ( any( class( error.covariance ) != "try-error" ) ){ dummy.matrix <- matrix( rep( 0, 4 ), nrow = 2, ncol = 2 ) dummy.matrix[ 1 ] <- error.covariance error.covariance <- dummy.matrix } } - if ( class( error.covariance ) == "try-error" ){ + if ( any( class( error.covariance ) == "try-error" ) ){ x.hessian <- numDeriv::hessian( likelihood, x = x$par, x.in = x$x, model = "gpd" ) diff --git a/R/opti.R b/R/opti.R index 4e96478..fb5655e 100644 --- a/R/opti.R +++ b/R/opti.R @@ -822,13 +822,13 @@ fit.gev.xts <- function( x, initial = NULL, res.optim$control$hessian[ 1 : 2, 1 : 2 ] ), silent = silent ) ## Augment the result again to ensure compatibility - if ( class( error.covariance ) != "try-error" ){ + if ( any( class( error.covariance ) != "try-error" ) ) { dummy.matrix <- matrix( rep( 0, 9 ), nrow = 3, ncol = 3 ) dummy.matrix[ 1 : 2, 1 : 2 ] <- error.covariance error.covariance <- dummy.matrix } } - if ( class( error.covariance ) == "try-error" || + if ( any( class( error.covariance ) == "try-error" ) || error.estimation == "MC" || any( is.nan( res.optim$control$hessian ) ) ){ parameter.estimate <- res.optim$par @@ -895,7 +895,7 @@ fit.gev.xts <- function( x, initial = NULL, method = "Nelder-Mead" ), x.in = yy, model = "gev" )$par ) ) ) } - if ( class( samples.fit ) == "try-error" ){ + if ( any( class( samples.fit ) == "try-error" ) ){ errors <- c( NaN, NaN, NaN, rep( NaN, length( return.period ) ) ) } else { @@ -1344,13 +1344,13 @@ fit.gev.default <- function( x, initial = NULL, res.optim$control$hessian[ 1 : 2, 1 : 2 ] ), silent = silent ) ## Augment the result again to ensure compatibility - if ( class( error.covariance ) != "try-error" ){ + if ( any( class( error.covariance ) != "try-error" ) ){ dummy.matrix <- matrix( rep( 0, 9 ), nrow = 3, ncol = 3 ) dummy.matrix[ 1 : 2, 1 : 2 ] <- error.covariance error.covariance <- dummy.matrix } } - if ( class( error.covariance ) == "try-error" || + if ( any( class( error.covariance ) == "try-error" ) || error.estimation == "MC" || any( is.nan( res.optim$control$hessian ) ) ){ parameter.estimate <- res.optim$par @@ -1417,7 +1417,7 @@ fit.gev.default <- function( x, initial = NULL, method = "Nelder-Mead" ), x.in = yy, model = "gev" )$par ) ) ) } - if ( class( samples.fit ) == "try-error" ){ + if ( any( class( samples.fit ) == "try-error" ) ){ errors <- c( NaN, NaN, NaN, rep( NaN, length( return.period ) ) ) } else { @@ -2390,13 +2390,13 @@ fit.gpd.xts <- function( x, initial = NULL, threshold = NULL, res.optim$control$hessian[ 1 ] ), silent = silent ) ## Augment the result again to ensure compatibility - if ( class( error.covariance ) != "try-error" ){ + if ( any( class( error.covariance ) != "try-error" ) ){ dummy.matrix <- matrix( rep( 0, 4 ), nrow = 2, ncol = 2 ) dummy.matrix[ 1 ] <- error.covariance error.covariance <- dummy.matrix } } - if ( class( error.covariance ) == "try-error" || + if ( any( class( error.covariance ) == "try-error" ) || error.estimation == "MC" || any( is.nan( res.optim$control$hessian ) ) ){ parameter.estimate <- res.optim$par @@ -2422,7 +2422,7 @@ fit.gpd.xts <- function( x, initial = NULL, threshold = NULL, control.outer = list( trace = FALSE, method = "Nelder-Mead" ), x.in = yy, model = "gpd" )$par ) ) ) - if ( class( samples.fit ) == "try-error" ){ + if ( any( class( samples.fit ) == "try-error" ) ){ errors <- c( NaN, NaN, NaN ) } else { errors <- data.frame( @@ -2953,13 +2953,13 @@ fit.gpd.default <- function( x, initial = NULL, threshold = NULL, res.optim$control$hessian[ 1 ] ), silent = silent ) ## Augment the result again to ensure compatibility - if ( class( error.covariance ) != "try-error" ){ + if ( any( class( error.covariance ) != "try-error" ) ){ dummy.matrix <- matrix( rep( 0, 4 ), nrow = 2, ncol = 2 ) dummy.matrix[ 1 ] <- error.covariance error.covariance <- dummy.matrix } } - if ( class( error.covariance ) == "try-error" || + if ( any( class( error.covariance ) == "try-error" ) || error.estimation == "MC" || any( is.nan( res.optim$control$hessian ) ) ){ parameter.estimate <- res.optim$par @@ -2984,7 +2984,7 @@ fit.gpd.default <- function( x, initial = NULL, threshold = NULL, control.outer = list( trace = FALSE, method = "Nelder-Mead" ), x.in = yy, model = "gpd" )$par ) ) ) - if ( class( samples.fit ) == "try-error" ){ + if ( any( class( samples.fit ) == "try-error" ) ){ errors <- c( NaN, NaN, NaN ) } else { errors <- data.frame( @@ -3615,7 +3615,7 @@ likelihood.initials <- function( x, model = c( "gev", "gpd" ), ## extRemes:::initializer.lmoments function lambda <- try( Lmoments::Lmoments( x ), silent = TRUE ) - if ( class( lambda ) == "try-error" ){ + if ( any( class( lambda ) == "try-error" ) ){ initial.lmom <- c( Inf, Inf, Inf ) } else { tau3 <- lambda[ 3 ]/ lambda[ 2 ] @@ -3635,7 +3635,7 @@ likelihood.initials <- function( x, model = c( "gev", "gpd" ), ## Linux principle), I will use the interior of the ## extRemes:::initializer.lmoments function lambda <- try( Lmoments::Lmoments( x ), silent = TRUE ) - if ( class( lambda ) == "try-error" ){ + if ( any( class( lambda ) == "try-error" ) ){ initial.lmom <- c( Inf, Inf, Inf ) } else { tau2 <- lambda[ 2 ]/ lambda[ 1 ] diff --git a/R/utils.R b/R/utils.R index 1763685..2924574 100644 --- a/R/utils.R +++ b/R/utils.R @@ -29,11 +29,11 @@ print.climex.fit.gev <- function( x, ... ){ ##' @author Philipp Mueller summary.climex.fit.gev <- function( object, ... ){ cat( "\n" ) - cat( paste( length( object$x ), "block maxima fitted using then " ) ) if ( object$control$error.estimation != "none" ){ - cat( paste( " Errors using", - object$control$error.estimation, "approach." ) ) + cat( paste( "Errors using", + object$control$error.estimation, "approach.\n\n" ) ) } + cat( paste( length( object$x ), "block maxima fitted." ) ) cat( "\n\n" ) cat( "\t\tFunction evaluations:\n" ) print( data.frame( function.eval = as.numeric( object$counts[ 1 ] ), From 90463867afa31e7a395bdf407646c0723cbdbdad Mon Sep 17 00:00:00 2001 From: theGreatWhiteShark Date: Mon, 11 Jul 2022 12:09:44 +0200 Subject: [PATCH 3/4] test: fix bootstrap unit tests there seems to have been some changes in the sampling algorithms or random number generation of R invalidating the previous bootstrap error estimates. Well, since there was an update to a new major release, this is not too unexpected. Also, the errors of the same order and still reproducible. I updated the values. --- tests/testthat/test_eva.R | 21 +++++++++++++++------ tests/testthat/test_fit.R | 14 ++++++++------ 2 files changed, 23 insertions(+), 12 deletions(-) diff --git a/tests/testthat/test_eva.R b/tests/testthat/test_eva.R index f7a7fc6..d6827f7 100644 --- a/tests/testthat/test_eva.R +++ b/tests/testthat/test_eva.R @@ -311,6 +311,7 @@ test_that( "return.level get GP error estimation right for MLE", { )$error ), c( 3.71517331, 9.08015755 ) ) }) + test_that( "return.level get the error estimation right for MC", { set.seed( 123 ) expect_equal( as.numeric( @@ -347,6 +348,7 @@ test_that( "return.level get the error estimation right for MC", { error.estimation = "MC", threshold = 29, monte.carlo.sample.size = 10 ) ) }) + test_that( "return.level get the error estimation right for bootstrap", { set.seed( 123 ) expect_equal( as.numeric( @@ -354,7 +356,8 @@ test_that( "return.level get the error estimation right for bootstrap", { return.period = c( 10, 100, 332 ), error.estimation = "bootstrap", bootstrap.sample.size = 10 )$error ), - c( 0.104627306, 0.273167480, 0.358395056 ) ) + c( 0.1785051, 0.4003212, 0.5322032 ), + tolerance = 1e-6 ) set.seed( 123 ) expect_equal( as.numeric( return.level( fit.gev( temp.potsdam, blocking = TRUE, @@ -362,13 +365,15 @@ test_that( "return.level get the error estimation right for bootstrap", { return.period = c( 10, 50 ), bootstrap.sample.size = 10, error.estimation = "bootstrap" )$error ), - c( 0.304337871, 0.390096183 ) ) + c( 0.2640825, 0.3143437 ), + tolerance = 1e-6 ) set.seed( 123 ) expect_equal( as.numeric( return.level( x.thresh.fit, return.period = 42, error.estimation = "bootstrap", threshold = 29, bootstrap.sample.size = 10 )$error ), - 0.185924244 ) + 0.1108912, + tolerance = 1e-6 ) set.seed( 123 ) expect_equal( as.numeric( return.level( fit.gpd( temp.potsdam, thresholding = TRUE, @@ -377,8 +382,10 @@ test_that( "return.level get the error estimation right for bootstrap", { return.period = c( 10, 50 ), bootstrap.sample.size = 10, error.estimation = "bootstrap" )$error ), - c( 0.345277289, 0.573791290 ) ) + c( 0.4137593, 0.4658717 ), + tolerance = 1e-6 ) }) + test_that( "return.level.climex.fit.gev yield equivalent results for minima and maxima", { expect_equal( climex::return.level( @@ -585,7 +592,8 @@ test_that( "the calculation of the upper limit does work", { set.seed( 2314 ) expect_equal( upper.limit( x.block.fit, error.estimation = "bootstrap" )$error, - 1.09028119 ) + 1.184065, + tolerance = 1e-6 ) expect_equal( upper.limit( x.thresh.fit, error.estimation = "MLE" )$error, 0.358538454 ) @@ -596,7 +604,8 @@ test_that( "the calculation of the upper limit does work", { set.seed( 2314 ) expect_equal( upper.limit( x.thresh.fit, error.estimation = "bootstrap" )$error, - 0.388303764 ) } ) + 0.3307957, + tolerance = 1e-6 ) } ) test_that( "upper.limit does only work for Weibull-typed data", { set.seed( 2231 ) diff --git a/tests/testthat/test_fit.R b/tests/testthat/test_fit.R index 4d302c4..f6ebae8 100644 --- a/tests/testthat/test_fit.R +++ b/tests/testthat/test_fit.R @@ -142,8 +142,8 @@ test_that( "fit.gev's error estimation works", { error.estimation = "bootstrap", return.period = c( 100, 200 ), bootstrap.sample.size = 10 )$se ), - c( 0.1586051250, 0.0740720149, 0.0523628935, 0.2731674804, - 0.3238840494 ) ) + c( 0.13234313, 0.09495922, 0.09392357, 0.40032124, 0.47744557 ), + tolerance = 1e-6 ) expect_equal( as.numeric( climex::fit.gev( temp.potsdam, blocking = TRUE, error.estimation = "MLE", @@ -168,8 +168,8 @@ test_that( "fit.gev's error estimation works", { extreme.type = "min", silent = TRUE, return.period = c( 100, 200 ), bootstrap.sample.size = 10 )$se ), - c( 0.3263569489, 0.2585539867, 0.0450737494, 0.4656176707, - 0.5520527130 ) ) + c( 0.40302692, 0.15256461, 0.04739434, 0.41928855, 0.53768576 ), + tolerance = 1e-6 ) }) test_that( "fit.gev yield equivalent results for minima and maxima", { @@ -389,7 +389,8 @@ test_that( "fit.gpd's error estimation works", { return.period = c( 100, 200 ), silent = TRUE, total.length = length( temp.potsdam ), bootstrap.sample.size = 10 )$se ), - c( 5.88442402e-01, 6.58027664e-07, 6.11622465e-01, 6.15380493e-01 ) ) + c( 6.091499e-01, 2.295979e-06, 6.331553e-01, 6.370458e-01 ), + tolerance = 1e-6 ) set.seed( 123 ) expect_equal( as.numeric( climex::fit.gpd( temp.potsdam, thresholding = TRUE, @@ -398,7 +399,8 @@ test_that( "fit.gpd's error estimation works", { return.period = c( 100, 200 ), bootstrap.sample.size = 10, total.length = length( temp.potsdam ) )$se ), - c( 0.593621431, 0.115000002, 0.494453578, 0.719987549 ) ) + c( 0.43555351, 0.06157277, 0.45968648, 0.50308282 ), + tolerance = 1e-6 ) }) test_that( "fit.gpd's threshold argument affect the result the way it is supposed to", { From 89b90d1b4cc0d18d210719fca664cf70af1af1fe Mon Sep 17 00:00:00 2001 From: theGreatWhiteShark Date: Mon, 11 Jul 2022 13:18:16 +0200 Subject: [PATCH 4/4] update ChangeLog/release notes --- DESCRIPTION | 8 ++++---- NEWS.md | 16 ++++++++++++++++ 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 2b547ba..71852d0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,12 +1,12 @@ Package: climex Type: Package Title: Estimating Extreme Events in Climatic Time Series (using Data of the DWD) -Version: 2.0.1 -Date: 2019-02-25 +Version: 2.0.2 +Date: 2022-07-11 LazyData: true Authors@R: person("Philipp Mueller", role = c("aut","cre"), - email = "thetruephil@googlemail.com") -Maintainer: Philipp Mueller + email = "princess.trudildis@posteo.de") +Maintainer: Philipp Mueller Description: The fitting of the generalized extreme value (GEV) and generalized Pareto (GP) distribution using maximum likelihood (ML) has an inherent tendency to produce numerical artifacts. By using the well-oiled framework of constrained optimization, in this we utilize the augmented Lagrangian method, the optimization find the global optimum significantly more often then the unconstrained one. In addition the package is focused on handling climatic time series of class 'xts' and provides a couple of convenience functions to handle these objects. URL: https://gitlab.com/theGreatWhiteShark/climex Depends: diff --git a/NEWS.md b/NEWS.md index cb83e98..54bd521 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,19 @@ +# v2.0.2 +- `deseasonalize` dependency (not available anymore) to allow + installing the package again +- Fix warnings in try-error checks. Some objects have more than a + single class/return a string list when provided to the `class()` + function. All checks whether the class of an error is "try-error" is + now wrapped in `any()` to ensure it both works and does not produce + any warnings +- Fix bootstrap unit tests. There seems to have been some changes in + the sampling algorithms or random number generation of R + invalidating the previous bootstrap error estimates. Well, since + there was an update to a new major release, this is not too + unexpected. Also, the errors of the same order and still + reproducible. I updated the values. + + # v2.0.1 - Updating the vignette. Removing the part describing the removed input function `source.data()` and the customization of the no