# Comparison of computation time

#### 2019-12-04

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

    ## R version 3.6.1 (2019-07-05)
    ## Platform: x86_64-pc-linux-gnu (64-bit)
    ## Running under: Ubuntu 18.04.3 LTS
    ##
    ## Matrix products: default
    ## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
    ## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
    ##
    ## locale:
    ##  [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C
    ##  [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=C
    ##  [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=de_DE.UTF-8
    ##  [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C
    ##  [9] LC_ADDRESS=C               LC_TELEPHONE=C
    ## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
    ##
    ## attached base packages:
    ## [1] stats     graphics  grDevices utils     datasets  methods   base
    ##
    ## other attached packages:
    ## [1] lmom_2.8        Lmoments_1.3-1  lmomco_2.3.2    TLMoments_0.7.5
    ## [5] Rcpp_1.0.3
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] digest_0.6.23   MASS_7.3-51.4   magrittr_1.5    evaluate_0.14
    ##  [5] rlang_0.4.2     stringi_1.4.3   goftest_1.2-2   rmarkdown_1.18
    ##  [9] tools_3.6.1     stringr_1.4.0   xfun_0.11       yaml_2.2.0
    ## [13] compiler_3.6.1  htmltools_0.4.0 knitr_1.26

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.037
    ## TLMoments_pwm         0.032
    ## TLMoments_recursive   0.026
    ## lmomco                0.509
    ## Lmoments              0.018

    # 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.242
    ## TLMoments_pwm         0.189
    ## TLMoments_recursive   0.054
    ## lmomco                9.276
    ## Lmoments              0.027

`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.018
    ## TLMoments_pwm         0.014
    ## TLMoments_recursive   0.007
    ## lmom                  0.008
    ## Lmoments              0.010

    # 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.217
    ## TLMoments_pwm         0.163
    ## TLMoments_recursive   0.034
    ## lmom                  0.019
    ## Lmoments              0.020

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.101
    ## TLMoments_pwm          0.033
    ## TLMoments_recurrence   0.025
    ## lmomco                 0.495

    # 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       1.730
    ## TLMoments_pwm          0.247
    ## TLMoments_recurrence   0.058
    ## lmomco                 9.040

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.099
    ## TLMoments_pwm          0.040
    ## TLMoments_recurrence   0.025
    ## lmomco                 0.449

    # 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       1.784
    ## TLMoments_pwm          0.459
    ## TLMoments_recurrence   0.076
    ## lmomco                 8.974

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.