Comparison of computation time

library(TLMoments)
library(lmomco)
library(Lmoments)
library(lmom)
sessionInfo()

## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] lmom_3.2          Lmoments_1.3-1    lmomco_2.5.1      TLMoments_0.7.5.3
## [5] Rcpp_1.0.14       rmarkdown_2.29   
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.37     R6_2.6.1          fastmap_1.2.0     xfun_0.51        
##  [5] maketools_1.3.2   cachem_1.1.0      knitr_1.49        htmltools_0.5.8.1
##  [9] buildtools_1.0.0  lifecycle_1.0.4   cli_3.6.4         goftest_1.2-3    
## [13] sass_0.4.9        jquerylib_0.1.4   compiler_4.4.2    sys_3.4.3        
## [17] tools_4.4.2       evaluate_1.0.3    bslib_0.9.0       yaml_2.3.10      
## [21] jsonlite_1.9.0    rlang_1.1.5       MASS_7.3-64

This document shows a comparison of computation time of TL-moments between different packages available, as well as between the different approaches built-in in this package.

This package offers the following computation methods (available via computation.method-attribute in TLMoments or TLMoment):

direct: Calculation as a weighted mean of the ordered data vector
pwm: Calculation of probabilty-weighted moments and using the conversion to TL-moments
recursive: An alternative recursive estimation of the weights of the direct approach
recurrence: Estimating the L-moments first and using the recurrence property to derive TL-moments

For a complete and thorough analysis of all these approaches and another speed comparison see Hosking & Balakrishnan (2015, A uniqueness result for L-estimators, with applications to L-moments, Statistical Methodology, 24, 69-80).

Besides our implementation, L-moments and/or TL-moments can be calculated using the packages

lmomco: L-moments and TL-moments
Lmoments: L-moments and TL(1,1)-moments
lmom: only L-moments

(all availabe at CRAN). The functions lmomco::lmoms, lmomco::TLmoms, and Lmoments::Lmoments return list objects with (T)L-moments and (T)L-moment-ratios and are therefore compared to our TLMoments. The function lmom::samlmu returns a vector of lambdas and is compared to the function TLMoment (which is a faster bare-bone function to compute TL-moments but is not suited to be transmitted to parameters or other functions of this package).

Calculation of L-moments

First we check if all calculation approaches in TLMoments give the same results (lmomco::lmoms is added as comparison):

n <- c(25, 50, 100, 200, 500, 1000, 10000, 50000)
sapply(n, function(nn) {
  x <- rgev(nn)
  check <- lmomco::lmoms(x, 4)$lambdas
  sapply(c("direct", "pwm", "recursive"), function(comp) {
    isTRUE(all.equal(TLMoment(x, order = 1:4, computation.method = comp), check, check.attributes = FALSE))
  })
})

##           [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## direct    TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## pwm       TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## recursive TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Now we compare the functions giving L-moments and L-moment-ratios simultaneously regarding computation speed:

possib <- list(
  TLMoments_direct = function(x) TLMoments(x, max.order = 4, computation.method = "direct"), 
  TLMoments_pwm = function(x) TLMoments(x, max.order = 4, computation.method = "pwm"), 
  TLMoments_recursive = function(x) TLMoments(x, max.order = 4, computation.method = "recursive"), 
  lmomco = function(x) lmomco::lmoms(x, 4), 
  Lmoments = function(x) Lmoments::Lmoments(x, returnobject = TRUE)
)

# n = 50
datalist <- replicate(200, rgev(50), simplify = FALSE)

do.call("rbind", lapply(possib, function(f) {
  system.time(lapply(datalist, f))[3]
}))

##                     elapsed
## TLMoments_direct      0.025
## TLMoments_pwm         0.022
## TLMoments_recursive   0.018
## lmomco                0.900
## Lmoments              0.009

# n = 1000
datalist <- replicate(200, evd::rgev(1000), simplify = FALSE)

do.call("rbind", lapply(possib, function(f) {
  system.time(lapply(datalist, f))[3]
}))

##                     elapsed
## TLMoments_direct      0.163
## TLMoments_pwm         0.125
## TLMoments_recursive   0.044
## lmomco               17.724
## Lmoments              0.020

Lmoments (since version 1.3-1) is the fastest implementation. Within TLMoments the recursive approach is the fastest. After this, the pwm approach is to be prefered over the direct approach. The implementation in lmomco is slow, compared to the others, especially for longer data vectors.

Comparison of functions that only return a vector of L-moments:

possib <- list(
  TLMoments_direct = function(x) TLMoment(x, order = 1:4, computation.method = "direct"), 
  TLMoments_pwm = function(x) TLMoment(x, order = 1:4, computation.method = "pwm"), 
  TLMoments_recursive = function(x) TLMoment(x, order = 1:4, computation.method = "recursive"), 
  lmom = function(x) lmom::samlmu(x, 4), 
  Lmoments = function(x) Lmoments::Lmoments(x, returnobject = FALSE)
)

# n = 50
datalist <- replicate(200, rgev(50), simplify = FALSE)

do.call("rbind", lapply(possib, function(f) {
  system.time(lapply(datalist, f))[3]
}))

##                     elapsed
## TLMoments_direct      0.012
## TLMoments_pwm         0.009
## TLMoments_recursive   0.005
## lmom                  0.004
## Lmoments              0.006

# n = 1000
datalist <- replicate(200, rgev(1000), simplify = FALSE)

do.call("rbind", lapply(possib, function(f) {
  system.time(lapply(datalist, f))[3]
}))

##                     elapsed
## TLMoments_direct      0.148
## TLMoments_pwm         0.104
## TLMoments_recursive   0.030
## lmom                  0.016
## Lmoments              0.016

For smaller data vectors our recursive-implementation is as fast as lmom::samlmu, but for longer data vectors lmom::samlmu and Lmoments::Lmoments are faster.

Calculation of TL-moments

Again, first we check if all approaches give the same results (lmomco::Tlmoms is added as comparison):

n <- c(25, 50, 100, 150, 200, 500, 1000, 10000)
names(n) <- paste("n", n, sep = "=")
sapply(n, function(nn) {
  x <- rgev(nn)
  check <- lmomco::TLmoms(x, 4, leftrim = 0, rightrim = 1)$lambdas
  sapply(c("direct", "pwm", "recursive", "recurrence"), function(comp) {
    tlm <- suppressWarnings(TLMoments(x, rightrim = 1, computation.method = comp)$lambdas)
    isTRUE(all.equal(tlm, check, check.attributes = FALSE))
  })
})

##            n=25 n=50 n=100 n=150 n=200 n=500 n=1000 n=10000
## direct     TRUE TRUE  TRUE  TRUE  TRUE  TRUE   TRUE    TRUE
## pwm        TRUE TRUE  TRUE  TRUE  TRUE  TRUE   TRUE    TRUE
## recursive  TRUE TRUE  TRUE  TRUE FALSE FALSE  FALSE   FALSE
## recurrence TRUE TRUE  TRUE  TRUE  TRUE  TRUE   TRUE    TRUE

sapply(n, function(nn) {
  x <- rgev(nn)
  check <- lmomco::TLmoms(x, 4, leftrim = 2, rightrim = 4)$lambdas
  sapply(c("direct", "pwm", "recursive", "recurrence"), function(comp) {
    tlm <- suppressWarnings(TLMoments(x, leftrim = 2, rightrim = 4, computation.method = comp)$lambdas)
    isTRUE(all.equal(tlm, check, check.attributes = FALSE))
  })
})

##            n=25 n=50 n=100 n=150 n=200 n=500 n=1000 n=10000
## direct     TRUE TRUE  TRUE  TRUE  TRUE  TRUE   TRUE    TRUE
## pwm        TRUE TRUE  TRUE  TRUE  TRUE  TRUE   TRUE    TRUE
## recursive  TRUE TRUE  TRUE  TRUE FALSE FALSE  FALSE   FALSE
## recurrence TRUE TRUE  TRUE  TRUE  TRUE  TRUE   TRUE    TRUE

The recursive approach fails when n exceeds 150. All other implementations give the same results.

Speed comparison for TL(0,1)-moments:

possib <- list(
  TLMoments_direct = function(x) TLMoments(x, leftrim = 0, rightrim = 1, max.order = 4, computation.method = "direct"), 
  TLMoments_pwm = function(x) TLMoments(x, leftrim = 0, rightrim = 1, max.order = 4, computation.method = "pwm"), 
  TLMoments_recurrence = function(x) TLMoments(x, leftrim = 0, rightrim = 1, max.order = 4, computation.method = "recurrence"), 
  lmomco = function(x) lmomco::TLmoms(x, 4, leftrim = 0, rightrim = 1)
)

# n = 50
datalist <- replicate(200, rgev(50), simplify = FALSE)

do.call("rbind", lapply(possib, function(f) {
  system.time(lapply(datalist, f))[3]
}))

##                      elapsed
## TLMoments_direct       0.067
## TLMoments_pwm          0.024
## TLMoments_recurrence   0.018
## lmomco                 0.896

# n = 1000
datalist <- replicate(200, rgev(1000), simplify = FALSE)

do.call("rbind", lapply(possib, function(f) {
  system.time(lapply(datalist, f))[3]
}))

##                      elapsed
## TLMoments_direct       0.947
## TLMoments_pwm          0.159
## TLMoments_recurrence   0.048
## lmomco                17.844

Speed comparison for TL(2,4)-moments:

possib <- list(
  TLMoments_direct = function(x) TLMoments(x, leftrim = 2, rightrim = 4, max.order = 4, computation.method = "direct"), 
  TLMoments_pwm = function(x) TLMoments(x, leftrim = 2, rightrim = 4, max.order = 4, computation.method = "pwm"), 
  TLMoments_recurrence = function(x) TLMoments(x, leftrim = 2, rightrim = 4, max.order = 4, computation.method = "recurrence"), 
  lmomco = function(x) lmomco::TLmoms(x, 4, leftrim = 2, rightrim = 4)
)

# n = 50
datalist <- replicate(200, evd::rgev(50), simplify = FALSE)

do.call("rbind", lapply(possib, function(f) {
  system.time(lapply(datalist, f))[3]
}))

##                      elapsed
## TLMoments_direct       0.066
## TLMoments_pwm          0.029
## TLMoments_recurrence   0.019
## lmomco                 0.821

# n = 1000
datalist <- replicate(200, evd::rgev(1000), simplify = FALSE)

do.call("rbind", lapply(possib, function(f) {
  system.time(lapply(datalist, f))[3]
}))

##                      elapsed
## TLMoments_direct       0.971
## TLMoments_pwm          0.273
## TLMoments_recurrence   0.063
## lmomco                17.967

In this calculations the recurrence approach clearly outperforms the other implementations. Calculation using probabilty-weighted moments is relatively fast, but using the direct calculation should be avoided, regarding calculation speed. This package’s implementation is clearly faster than those in lmomco.

Conclusion

This results encourage to use the recursive approach for L-moments and the recurrence approach when calculating TL-moments. Therefore these are the defaults in this package, but the other computation methods (direct and pwm) are still available (by using the argument computation.method).

In comparison to other packages Lmoments is faster but only supports L-moments and TL(1,1)-moments.