Merge remote-tracking branch 'origin/develop'

theGreatWhiteShark · Jul 11, 2022 · 3a11d41 · 3a11d41
2 parents 3ac0ed4 + 89b90d1
commit 3a11d41
Show file tree

Hide file tree

Showing 7 changed files with 68 additions and 42 deletions.
diff --git 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 = "[email protected]")
-Maintainer: Philipp Mueller <[email protected]>
+    email = "[email protected]")
+Maintainer: Philipp Mueller <[email protected]>
 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:
@@ -18,7 +18,6 @@ Imports:
     grid,
     Rcpp,
     RColorBrewer,
-    deseasonalize,
     Lmoments,
     moments,
     numDeriv,

diff --git 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

diff --git 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
@@ -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
@@ -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 ] ),

diff --git 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,28 +348,32 @@ 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(
       return.level( x.block.fit,
                    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,
                             extreme.type = "min", silent = TRUE ),
                    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
@@ -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", {