tests/testthat/test-subbofit.R

In Rsubbotools: Fast Estimation of Subbottin and AEP Distributions (Generalized Error Distribution)

##############################################################################

# Subbofit
# checks against original implementation of subbotools
# Depends on using the same sample to fit as the subbotools package
# since R changes the rng in each version, to work this test must be manually
# updated each time or use a constant sample

skip()
paste0("Subbofit")

test_that("Subbofit Method of moments:", {
  set.seed(1)
  x <- rsubbo(1e6, 0, 1, 1.5)
  para_mm <- sd(x) / (sum(abs(x - mean(x))) / length(x))

  mm_objf <- function(x, l) {
    # lgamma(x) == log(abs(gamma(x)))
    dtmp1 <- exp(x)
    y <- .5 * lgamma(3. * dtmp1) + .5 * lgamma(dtmp1) -
      lgamma(2. * dtmp1) - (l)

    return(y)
  }

  mm_objdf <- function(x) {
    dtmp1 <- exp(x)
    dy <- dtmp1 * (1.5 * digamma(3. * dtmp1) + .5 * digamma(dtmp1) -
      2. * digamma(2. * dtmp1))

    return(dy)
  }


  mm_fl <- function(l) {
    function(x) {
      mm_objf(x, l)
    }
  }

  mm_f <- mm_fl(log(para_mm))

  test_steffensen <- steffensen(0, mm_f, mm_objdf, max_iter = 3)

  expect_equal(mm_subbotin(para_mm), exp(-test_steffensen$x[3]))
})

##############################################################################

paste0("Test fit against subfamilies of the AEP\n")

paste0(
  " This test check if the following fits work as expected:
- sepfit
- subboafish
- subboafit
- subbofish
- subbofit
- subbolafit
- laplafit
- alaplafit

The tests are conducted in the following way.
First, we generate 5 datasets contemplating each of the special cases of the
rpower generating function, which are related to the b parameter:
- b < 1 or b >4
- b == 1 (Laplacian special case)
- 1 < b < 2
- b == 2 (Gaussian special case)
- 2 < b < 4

We then run each of the original subbotools functions on these datasets with
default parameters (we hard code the results on the test) and compare the
answers against the Rsubbotools implementation.

"
)

test_that("SubLaplace:", {
  # subbofit -V 1 < sublaplace.txt
  #
  #
  #--- FINAL RESULT ----------------------------------
  #
  #     value     std.err     |  b       a       m
  # b =  0.4993    0.0008212  |    --   0.0000  0.0000
  # a =  3.995     0.006409   |  0.4709    --   0.0000
  # m =  0.0006382 -nan       | -nan -nan    --
  #
  #           Upper triangle: covariances
  #           Lower triangle: correlation coefficients
  #---------------------------------------------------
  #
  # b            a            m            log-like
  # 4.9932e-01   3.9950e+00   6.3821e-04   3.3857e+00

  orig_value <-
    generate_orig_dt(
      coef = c(4.9932e-01, 3.9950e+00, 6.3821e-04),
      log_likelihood = 3.3857e+00,
      std_error = c(0.0008212, 0.006409, NaN)
      # we pass the transposed matrix and the code corrects it
      , matrix =
        c(
          NA, 0.0000, 0.0000,
          0.4709, NA, 0.0000,
          NaN, NaN, NA
        ),
      distribution = "subbofit"
    )

  check_fits(orig_value, .5, subbofit)
})

test_that("Laplace:", {
  # subbofit -V 1 < laplace.txt
  #
  #--- FINAL RESULT ----------------------------------
  #
  #     value     std.err     |  b       a       m
  # b =  1.002     0.001861   |    --   0.0000 -0.0000
  # a =  0.9998    0.001271   |  0.6180    --  -0.0000
  # m =  0.0007125 0.001001   | -0.0000 -0.0000    --
  #
  #           Upper triangle: covariances
  #           Lower triangle: correlation coefficients
  #---------------------------------------------------
  #
  # b            a            m            log-like
  # 1.0017e+00   9.9983e-01   7.1252e-04   1.6923e+00

  orig_value <-
    generate_orig_dt(
      coef = c(1.0017e+00, 9.9983e-01, 7.1252e-04),
      log_likelihood = 1.6923e+00,
      std_error = c(0.001861, 0.001271, 0.001001),
      matrix =
        c(
          NA, 0.0000, -0.0000,
          0.6180, NA, -0.0000,
          -0.0000, -0.0000, NA
        ),
      distribution = "subbofit"
    )
  check_fits(orig_value, 1, subbofit)
})

test_that("Subnormal:", {
  # subbofit -V 1 < subnormal.txt
  #
  #--- FINAL RESULT ----------------------------------
  #
  #     value     std.err     |  b       a       m
  # b =  1.504     0.003088   |    --   0.0000 -0.0000
  # a =  0.7637    0.0008741  |  0.7018    --  -0.0000
  # m = -0.0002804 0.0008213  | -0.0000 -0.0000    --
  #
  #           Upper triangle: covariances
  #           Lower triangle: correlation coefficients
  #---------------------------------------------------
  #
  # b            a            m            log-like
  # 1.5041e+00   7.6375e-01  -2.8042e-04   1.2572e+00

  orig_value <-
    generate_orig_dt(
      coef = c(1.5041e+00, 7.6375e-01, -2.8042e-04),
      log_likelihood = 1.2572e+00,
      std_error = c(0.003088, 0.0008741, 0.0008213),
      matrix =
        c(
          NA, 0.0000, -0.0000,
          0.7018, NA, -0.0000,
          -0.0000, -0.0000, NA
        ),
      distribution = "subbofit"
    )
  check_fits(orig_value, 1.5, subbofit)
})

test_that("Normal:", {
  # subbofit -V 1 < normal.txt
  #
  #--- FINAL RESULT ----------------------------------
  #
  #     value     std.err     |  b       a       m
  # b =  2.004     0.004473   |    --   0.0000 -0.0000
  # a =  0.7078    0.0007626  |  0.7551    --  -0.0000
  # m =  3.197e-05 0.0007072  | -0.0000 -0.0000    --
  #
  #           Upper triangle: covariances
  #           Lower triangle: correlation coefficients
  #---------------------------------------------------
  #
  # b            a            m            log-like
  #  2.0044e+00   7.0781e-01   3.1969e-05   1.0725e+00

  orig_value <-
    generate_orig_dt(
      coef = c(2.0044e+00, 7.0781e-01, 3.1969e-05),
      log_likelihood = 1.0725e+00,
      std_error = c(0.004473, 0.0007626, 0.0007072),
      matrix =
        c(
          NA, 0.0000, -0.0000,
          0.7551, NA, -0.0000,
          -0.0000, -0.0000, NA
        ),
      distribution = "subbofit"
    )
  check_fits(orig_value, 2, subbofit)
})

test_that("SuperNormal:", {
  # subbofit -V 1 < supernormal.txt
  #
  #--- FINAL RESULT ----------------------------------
  #
  #     value     std.err     |  b       a       m
  # b =  2.499     0.005986   |    --   0.0000 -0.0000
  # a =  0.6929    0.0007172  |  0.7915    --  -0.0000
  # m = -0.0003205 0.0006301  | -0.0000 -0.0000    --
  #
  #           Upper triangle: covariances
  #           Lower triangle: correlation coefficients
  #---------------------------------------------------
  #
  # b            a            m            log-like
  #  2.4988e+00   6.9286e-01  -3.2053e-04   9.7328e-01

  orig_value <-
    generate_orig_dt(
      coef = c(2.4988e+00, 6.9286e-01, -3.2053e-04),
      log_likelihood = 9.7328e-01,
      std_error = c(0.005986, 0.0007172, 0.0006301),
      matrix =
        c(
          NA, 0.0000, -0.0000,
          0.7915, NA, -0.0000,
          -0.0000, -0.0000, NA
        ),
      distribution = "subbofit"
    )
  check_fits(orig_value, 2.5, subbofit)
})

##############################################################################