test_that("Constructor", {
  x <- lltd$D$x
  y <- lltd$D$y

  m <- length(unique(x))
  n <- length(y)

  w <- rep(1, n)

  max_iter <- 10000

  stats <- lltd$stats_1

  start <- c(0, 1, 1, 1, 1, 1)

  lower_bound <- c(0, -1, 0.5, 1, 0, 0.5)
  upper_bound <- c(3, 2, 2, 5, 2, 1)

  object <- loglogistic6_new(x, y, w, NULL, max_iter, NULL, NULL)

  expect_true(inherits(object, "loglogistic6"))
  expect_equal(object$x, x)
  expect_equal(object$y, y)
  expect_equal(object$w, w)
  expect_equal(object$n, n)
  expect_equal(object$m, m)
  expect_equal(object$stats, stats)
  expect_false(object$constrained)
  expect_equal(object$max_iter, max_iter)
  expect_null(object$start)
  expect_null(object$lower_bound)
  expect_null(object$upper_bound)

  object <- loglogistic6_new(x, y, w, start, max_iter, lower_bound, upper_bound)

  i <- c(1, 2)

  expect_true(inherits(object, "loglogistic6"))
  expect_equal(object$x, x)
  expect_equal(object$y, y)
  expect_equal(object$w, w)
  expect_equal(object$n, n)
  expect_equal(object$m, m)
  expect_equal(object$stats, stats)
  expect_true(object$constrained)
  expect_equal(object$max_iter, max_iter)
  expect_equal(object$start, c(start[i], log(start[-i])))
  expect_equal(object$lower_bound, c(lower_bound[i], log(lower_bound[-i])))
  expect_equal(object$upper_bound, c(upper_bound[i], log(upper_bound[-i])))

  w <- lltd$D$w
  stats <- lltd$stats_2

  object <- loglogistic6_new(x, y, w, NULL, max_iter, NULL, NULL)

  expect_true(inherits(object, "loglogistic6"))
  expect_equal(object$x, x)
  expect_equal(object$y, y)
  expect_equal(object$w, w)
  expect_equal(object$n, n)
  expect_equal(object$m, m)
  expect_equal(object$stats, stats)
  expect_false(object$constrained)
  expect_equal(object$max_iter, max_iter)
  expect_null(object$start)
  expect_null(object$lower_bound)
  expect_null(object$upper_bound)

  object <- loglogistic6_new(x, y, w, start, max_iter, lower_bound, upper_bound)

  expect_true(inherits(object, "loglogistic6"))
  expect_equal(object$x, x)
  expect_equal(object$y, y)
  expect_equal(object$w, w)
  expect_equal(object$n, n)
  expect_equal(object$m, m)
  expect_equal(object$stats, stats)
  expect_true(object$constrained)
  expect_equal(object$max_iter, max_iter)
  expect_equal(object$start, c(start[i], log(start[-i])))
  expect_equal(object$lower_bound, c(lower_bound[i], log(lower_bound[-i])))
  expect_equal(object$upper_bound, c(upper_bound[i], log(upper_bound[-i])))
})

test_that("Constructor: errors", {
  x <- lltd$D$x
  y <- lltd$D$y
  w <- lltd$D$w
  max_iter <- 10000

  expect_error(
    loglogistic6_new(x, y, w, c(0, 1, 1, 1, 1), max_iter, NULL, NULL),
    "'start' must be of length 6"
  )

  expect_error(
    loglogistic6_new(x, y, w, c(0, 1, 1, 1, 1, 1, 1), max_iter, NULL, NULL),
    "'start' must be of length 6"
  )

  expect_error(
    loglogistic6_new(x, y, w, c(0, 1, 0, 1, 1, 1), max_iter, NULL, NULL),
    "parameter 'eta' cannot be negative nor zero"
  )

  expect_error(
    loglogistic6_new(x, y, w, c(0, 1, -1, 1, 1, 1), max_iter, NULL, NULL),
    "parameter 'eta' cannot be negative nor zero"
  )

  expect_error(
    loglogistic6_new(x, y, w, c(0, 1, 1, 0, 1, 1), max_iter, NULL, NULL),
    "parameter 'phi' cannot be negative nor zero"
  )

  expect_error(
    loglogistic6_new(x, y, w, c(0, 1, 1, -1, 1, 1), max_iter, NULL, NULL),
    "parameter 'phi' cannot be negative nor zero"
  )

  expect_error(
    loglogistic6_new(x, y, w, c(0, 1, 1, 1, 0, 1), max_iter, NULL, NULL),
    "parameter 'nu' cannot be negative nor zero"
  )

  expect_error(
    loglogistic6_new(x, y, w, c(0, 1, 1, 1, -1, 1), max_iter, NULL, NULL),
    "parameter 'nu' cannot be negative nor zero"
  )

  expect_error(
    loglogistic6_new(x, y, w, c(0, 1, 1, 1, 1, 0), max_iter, NULL, NULL),
    "parameter 'xi' cannot be negative nor zero"
  )

  expect_error(
    loglogistic6_new(x, y, w, c(0, 1, 1, 1, 1, -1), max_iter, NULL, NULL),
    "parameter 'xi' cannot be negative nor zero"
  )

  expect_error(
    loglogistic6_new(x, y, w, NULL, max_iter, rep(-Inf, 5), rep(Inf, 5)),
    "'lower_bound' must be of length 6"
  )

  expect_error(
    loglogistic6_new(x, y, w, NULL, max_iter, rep(-Inf, 5), rep(Inf, 6)),
    "'lower_bound' must be of length 6"
  )

  expect_error(
    loglogistic6_new(x, y, w, NULL, max_iter, rep(-Inf, 6), rep(Inf, 5)),
    "'upper_bound' must be of length 6"
  )

  expect_error(
    loglogistic6_new(
      x, y, w, NULL, max_iter, rep(-Inf, 6), c(1, 1, 0, rep(Inf, 3))
    ),
    "'upper_bound[3]' cannot be negative nor zero",
    fixed = TRUE
  )

  expect_error(
    loglogistic6_new(
      x, y, w, NULL, max_iter, rep(-Inf, 6), c(1, 1, -1, rep(Inf, 3))
    ),
    "'upper_bound[3]' cannot be negative nor zero",
    fixed = TRUE
  )

  expect_error(
    loglogistic6_new(
      x, y, w, NULL, max_iter, rep(-Inf, 6), c(1, 1, Inf, 0, Inf, Inf)
    ),
    "'upper_bound[4]' cannot be negative nor zero",
    fixed = TRUE
  )

  expect_error(
    loglogistic6_new(
      x, y, w, NULL, max_iter, rep(-Inf, 6), c(1, 1, Inf, -1, Inf, Inf)
    ),
    "'upper_bound[4]' cannot be negative nor zero",
    fixed = TRUE
  )

  expect_error(
    loglogistic6_new(
      x, y, w, NULL, max_iter, rep(-Inf, 6), c(1, Inf, Inf, Inf, 0, Inf)
    ),
    "'upper_bound[5]' cannot be negative nor zero",
    fixed = TRUE
  )

  expect_error(
    loglogistic6_new(
      x, y, w, NULL, max_iter, rep(-Inf, 6), c(1, Inf, Inf, Inf, -1, Inf)
    ),
    "'upper_bound[5]' cannot be negative nor zero",
    fixed = TRUE
  )

  expect_error(
    loglogistic6_new(
      x, y, w, NULL, max_iter, rep(-Inf, 6), c(1, rep(Inf, 4), 0)
    ),
    "'upper_bound[6]' cannot be negative nor zero",
    fixed = TRUE
  )

  expect_error(
    loglogistic6_new(
      x, y, w, NULL, max_iter, rep(-Inf, 6), c(1, rep(Inf, 4), -1)
    ),
    "'upper_bound[6]' cannot be negative nor zero",
    fixed = TRUE
  )
})

test_that("Function value", {
  x <- lltd$stats_1[, 1]
  theta <- lltd$theta_6

  m <- length(x)

  true_value <- c(
    0.8, 0.18521184704873563, 0.092893218813452476, 0.065959193177705752,
    0.055218021012035280, 0.049970700263641732, 0.040265592016188420,
    0.040165327459277954
  )

  value <- loglogistic6_fn(x, theta)

  expect_type(value, "double")
  expect_length(value, m)
  expect_equal(value, true_value)

  object <- structure(list(stats = lltd$stats_1), class = "loglogistic6")

  value <- fn(object, object$stats[, 1], theta)

  expect_type(value, "double")
  expect_length(value, m)
  expect_equal(value, true_value)

  object <- structure(list(stats = lltd$stats_1), class = "loglogistic6_fit")

  value <- fn(object, object$stats[, 1], theta)

  expect_type(value, "double")
  expect_length(value, m)
  expect_equal(value, true_value)
})

test_that("Gradient (1)", {
  x <- lltd$stats_1[, 1]
  theta <- lltd$theta_6

  m <- length(x)

  true_gradient <- matrix(
    c(
      # alpha
      rep(1, m),
      # delta
      0, 0.61478815295126437, 0.70710678118654752, 0.73404080682229425,
      0.74478197898796472, 0.75002929973635827, 0.75973440798381158,
      0.75983467254072205,
      # eta
      0, 0, -0.030633066983392100, -0.026013750024478273, -0.019855520341052468,
      -0.015280070764963223, -0.00039606856104051224, -6.2960660335592149e-06,
      # phi
      0, 0.087826878993037766, 0.044194173824159220, 0.023678735703944976,
      0.014322730365153168, 0.0094940417688146616, 0.00010124392430487894,
      1.0131074934809975e-06,
      # nu
      0, -0.052813436898417783, -0.050217590510744394, -0.050819955866517416,
      -0.051284379953673178, -0.051561660317557922, -0.052165854044333260,
      -0.052172731788759549,
      # xi
      0, 0.021956719748259442, 0.044194173824159220, 0.053277155333876195,
      0.057290921460612671, 0.059337761055091635, 0.063277452690549339,
      0.063319218342562343
    ),
    nrow = m,
    ncol = 6
  )

  G <- loglogistic6_gradient(x, theta)

  expect_type(G, "double")
  expect_length(G, m * 6)
  expect_equal(G, true_gradient)
})

test_that("Hessian (1)", {
  x <- lltd$stats_1[, 1]
  theta <- lltd$theta_6

  m <- length(x)

  true_hessian <- array(
    c(
      # (alpha, alpha)
      rep(0, m),
      # (alpha, delta)
      rep(0, m),
      # (alpha, eta)
      rep(0, m),
      # (alpha, phi)
      rep(0, m),
      # (alpha, nu)
      rep(0, m),
      # (alpha, xi)
      rep(0, m),
      # (delta, alpha)
      rep(0, m),
      # (delta, delta)
      rep(0, m),
      # (delta, eta)
      0, 0, 0.030633066983392100, 0.026013750024478273, 0.019855520341052468,
      0.015280070764963223, 0.00039606856104051224, 6.2960660335592149e-06,
      # (delta, phi)
      0, -0.087826878993037766, -0.044194173824159220, -0.023678735703944976,
      -0.014322730365153168, -0.0094940417688146616, -0.00010124392430487894,
      -1.0131074934809975e-06,
      # (delta, nu)
      0, 0.052813436898417783, 0.050217590510744394, 0.050819955866517416,
      0.051284379953673178, 0.051561660317557922, 0.052165854044333260,
      0.052172731788759549,
      # (delta, xi)
      0, -0.021956719748259442, -0.044194173824159220, -0.053277155333876195,
      -0.057290921460612671, -0.059337761055091635, -0.063277452690549339,
      -0.063319218342562343,
      # (eta, alpha)
      rep(0, m),
      # (eta, delta)
      0, 0, 0.030633066983392100, 0.026013750024478273, 0.019855520341052468,
      0.015280070764963223, 0.00039606856104051224, 6.2960660335592149e-06,
      # (eta, eta)
      0, 0, 0.014597841507866501, 0.023969505217162465, 0.024878903973796331,
      0.023035848915718639, 0.0015483969202501701, 0.000039127322111196110,
      # (eta, phi)
      0, 0.043913439496518883, 0.0010368533609975416, -0.0099786160395254187,
      -0.010784970510297762, -0.0095659567941658146, -0.00034518269393002339,
      -5.7894703132690180e-06,
      # (eta, nu)
      0, 0, -0.00026094464335276508, -0.00096186058064673394,
      -0.00098537935388951183, -0.00085702810352105560,
      -0.000027142584296135141, -4.3230007496653412e-07,
      # (eta, xi)
      0, 0, 0.0095728334323100312, 0.0094404737992058250, 0.0076367385927124878,
      0.0060443317899379837, 0.00016494059878086364, 2.6233468561331068e-06,
      # (phi, alpha)
      rep(0, m),
      # (phi, delta)
      0, -0.087826878993037766, -0.044194173824159220, -0.023678735703944976,
      -0.014322730365153168, -0.0094940417688146616, -0.00010124392430487894,
      -1.0131074934809975e-06,
      # (phi, eta)
      0, 0.043913439496518883, 0.0010368533609975416, -0.0099786160395254187,
      -0.010784970510297762, -0.0095659567941658146, -0.00034518269393002339,
      -5.7894703132690180e-06,
      # (phi, phi)
      0, -0.018820045498508093, 0.0082864075920298538, 0.0080202169319813628,
      0.0057841795705426254, 0.0041461321648620991, 0.000050554502181554024,
      5.0654699272656362e-07,
      # (phi, nu)
      0, -0.0050019202992314262, 0.00037646354291157338, 0.00087552323105072386,
      0.00071080095362538481, 0.00053250150061700330, 6.9382476172745710e-06,
      6.9561920578111675e-08,
      # (phi, xi)
      0, -0.015683371248756744, -0.013810679320049756, -0.0085930895699800315,
      -0.0055087424481358337, -0.0037555544971576984, -0.000042162481803404410,
      -4.2212587094577058e-07,
      # (nu, alpha)
      rep(0, m),
      # (nu, delta)
      0, 0.052813436898417783, 0.050217590510744394, 0.050819955866517416,
      0.051284379953673178, 0.051561660317557922, 0.052165854044333260,
      0.052172731788759549,
      # (nu, eta)
      0, 0, -0.00026094464335276508, -0.00096186058064673394,
      -0.00098537935388951183, -0.00085702810352105560,
      -0.000027142584296135141, -4.3230007496653412e-07,
      # (nu, phi)
      0, -0.0050019202992314262, 0.00037646354291157338, 0.00087552323105072386,
      0.00071080095362538481, 0.00053250150061700330, 6.9382476172745710e-06,
      6.9561920578111675e-08,
      # (nu, nu)
      0, 0.018733100695205814, 0.020851888481470548, 0.021700594857954726,
      0.022041978298575419, 0.022206117780312320, 0.022501045166927927,
      0.022504015626713413,
      # (nu, xi)
      0, -0.0067396600118727170, -0.010672079913128232, -0.011349361563604920,
      -0.011479526550651628, -0.011506305884916638, -0.011482958411840728,
      -0.011482184549508606,
      # (xi, alpha)
      rep(0, m),
      # (xi, delta)
      0, -0.021956719748259442, -0.044194173824159220, -0.053277155333876195,
      -0.057290921460612671, -0.059337761055091635, -0.063277452690549339,
      -0.063319218342562343,
      # (xi, eta)
      0, 0, 0.0095728334323100312, 0.0094404737992058250, 0.0076367385927124878,
      0.0060443317899379837, 0.00016494059878086364, 2.6233468561331068e-06,
      # (xi, phi)
      0, -0.015683371248756744, -0.013810679320049756, -0.0085930895699800315,
      -0.0055087424481358337, -0.0037555544971576984, -0.000042162481803404410,
      -4.2212587094577058e-07,
      # (xi, nu)
      0, -0.0067396600118727170, -0.010672079913128232, -0.011349361563604920,
      -0.011479526550651628, -0.011506305884916638, -0.011482958411840728,
      -0.011482184549508606,
      # (xi, xi)
      0, -0.0039208428121891860, -0.013810679320049756, -0.019334451532455071,
      -0.022034969792543335, -0.023472215607235615, -0.026351551127127756,
      -0.026382866934110661
    ),
    dim = c(m, 6, 6)
  )

  H <- loglogistic6_hessian(x, theta)

  expect_type(H, "double")
  expect_length(H, m * 6 * 6)
  expect_equal(H, true_hessian)
})

test_that("Gradient and Hessian (1)", {
  x <- lltd$stats_1[, 1]
  theta <- lltd$theta_6

  m <- length(x)

  true_gradient <- matrix(
    c(
      # alpha
      rep(1, m),
      # delta
      0, 0.61478815295126437, 0.70710678118654752, 0.73404080682229425,
      0.74478197898796472, 0.75002929973635827, 0.75973440798381158,
      0.75983467254072205,
      # eta
      0, 0, -0.030633066983392100, -0.026013750024478273, -0.019855520341052468,
      -0.015280070764963223, -0.00039606856104051224, -6.2960660335592149e-06,
      # phi
      0, 0.087826878993037766, 0.044194173824159220, 0.023678735703944976,
      0.014322730365153168, 0.0094940417688146616, 0.00010124392430487894,
      1.0131074934809975e-06,
      # nu
      0, -0.052813436898417783, -0.050217590510744394, -0.050819955866517416,
      -0.051284379953673178, -0.051561660317557922, -0.052165854044333260,
      -0.052172731788759549,
      # xi
      0, 0.021956719748259442, 0.044194173824159220, 0.053277155333876195,
      0.057290921460612671, 0.059337761055091635, 0.063277452690549339,
      0.063319218342562343
    ),
    nrow = m,
    ncol = 6
  )

  true_hessian <- array(
    c(
      # (alpha, alpha)
      rep(0, m),
      # (alpha, delta)
      rep(0, m),
      # (alpha, eta)
      rep(0, m),
      # (alpha, phi)
      rep(0, m),
      # (alpha, nu)
      rep(0, m),
      # (alpha, xi)
      rep(0, m),
      # (delta, alpha)
      rep(0, m),
      # (delta, delta)
      rep(0, m),
      # (delta, eta)
      0, 0, 0.030633066983392100, 0.026013750024478273, 0.019855520341052468,
      0.015280070764963223, 0.00039606856104051224, 6.2960660335592149e-06,
      # (delta, phi)
      0, -0.087826878993037766, -0.044194173824159220, -0.023678735703944976,
      -0.014322730365153168, -0.0094940417688146616, -0.00010124392430487894,
      -1.0131074934809975e-06,
      # (delta, nu)
      0, 0.052813436898417783, 0.050217590510744394, 0.050819955866517416,
      0.051284379953673178, 0.051561660317557922, 0.052165854044333260,
      0.052172731788759549,
      # (delta, xi)
      0, -0.021956719748259442, -0.044194173824159220, -0.053277155333876195,
      -0.057290921460612671, -0.059337761055091635, -0.063277452690549339,
      -0.063319218342562343,
      # (eta, alpha)
      rep(0, m),
      # (eta, delta)
      0, 0, 0.030633066983392100, 0.026013750024478273, 0.019855520341052468,
      0.015280070764963223, 0.00039606856104051224, 6.2960660335592149e-06,
      # (eta, eta)
      0, 0, 0.014597841507866501, 0.023969505217162465, 0.024878903973796331,
      0.023035848915718639, 0.0015483969202501701, 0.000039127322111196110,
      # (eta, phi)
      0, 0.043913439496518883, 0.0010368533609975416, -0.0099786160395254187,
      -0.010784970510297762, -0.0095659567941658146, -0.00034518269393002339,
      -5.7894703132690180e-06,
      # (eta, nu)
      0, 0, -0.00026094464335276508, -0.00096186058064673394,
      -0.00098537935388951183, -0.00085702810352105560,
      -0.000027142584296135141, -4.3230007496653412e-07,
      # (eta, xi)
      0, 0, 0.0095728334323100312, 0.0094404737992058250, 0.0076367385927124878,
      0.0060443317899379837, 0.00016494059878086364, 2.6233468561331068e-06,
      # (phi, alpha)
      rep(0, m),
      # (phi, delta)
      0, -0.087826878993037766, -0.044194173824159220, -0.023678735703944976,
      -0.014322730365153168, -0.0094940417688146616, -0.00010124392430487894,
      -1.0131074934809975e-06,
      # (phi, eta)
      0, 0.043913439496518883, 0.0010368533609975416, -0.0099786160395254187,
      -0.010784970510297762, -0.0095659567941658146, -0.00034518269393002339,
      -5.7894703132690180e-06,
      # (phi, phi)
      0, -0.018820045498508093, 0.0082864075920298538, 0.0080202169319813628,
      0.0057841795705426254, 0.0041461321648620991, 0.000050554502181554024,
      5.0654699272656362e-07,
      # (phi, nu)
      0, -0.0050019202992314262, 0.00037646354291157338, 0.00087552323105072386,
      0.00071080095362538481, 0.00053250150061700330, 6.9382476172745710e-06,
      6.9561920578111675e-08,
      # (phi, xi)
      0, -0.015683371248756744, -0.013810679320049756, -0.0085930895699800315,
      -0.0055087424481358337, -0.0037555544971576984, -0.000042162481803404410,
      -4.2212587094577058e-07,
      # (nu, alpha)
      rep(0, m),
      # (nu, delta)
      0, 0.052813436898417783, 0.050217590510744394, 0.050819955866517416,
      0.051284379953673178, 0.051561660317557922, 0.052165854044333260,
      0.052172731788759549,
      # (nu, eta)
      0, 0, -0.00026094464335276508, -0.00096186058064673394,
      -0.00098537935388951183, -0.00085702810352105560,
      -0.000027142584296135141, -4.3230007496653412e-07,
      # (nu, phi)
      0, -0.0050019202992314262, 0.00037646354291157338, 0.00087552323105072386,
      0.00071080095362538481, 0.00053250150061700330, 6.9382476172745710e-06,
      6.9561920578111675e-08,
      # (nu, nu)
      0, 0.018733100695205814, 0.020851888481470548, 0.021700594857954726,
      0.022041978298575419, 0.022206117780312320, 0.022501045166927927,
      0.022504015626713413,
      # (nu, xi)
      0, -0.0067396600118727170, -0.010672079913128232, -0.011349361563604920,
      -0.011479526550651628, -0.011506305884916638, -0.011482958411840728,
      -0.011482184549508606,
      # (xi, alpha)
      rep(0, m),
      # (xi, delta)
      0, -0.021956719748259442, -0.044194173824159220, -0.053277155333876195,
      -0.057290921460612671, -0.059337761055091635, -0.063277452690549339,
      -0.063319218342562343,
      # (xi, eta)
      0, 0, 0.0095728334323100312, 0.0094404737992058250, 0.0076367385927124878,
      0.0060443317899379837, 0.00016494059878086364, 2.6233468561331068e-06,
      # (xi, phi)
      0, -0.015683371248756744, -0.013810679320049756, -0.0085930895699800315,
      -0.0055087424481358337, -0.0037555544971576984, -0.000042162481803404410,
      -4.2212587094577058e-07,
      # (xi, nu)
      0, -0.0067396600118727170, -0.010672079913128232, -0.011349361563604920,
      -0.011479526550651628, -0.011506305884916638, -0.011482958411840728,
      -0.011482184549508606,
      # (xi, xi)
      0, -0.0039208428121891860, -0.013810679320049756, -0.019334451532455071,
      -0.022034969792543335, -0.023472215607235615, -0.026351551127127756,
      -0.026382866934110661
    ),
    dim = c(m, 6, 6)
  )

  gh <- loglogistic6_gradient_hessian(x, theta)

  expect_type(gh, "list")
  expect_type(gh$G, "double")
  expect_type(gh$H, "double")

  expect_length(gh$G, m * 6)
  expect_length(gh$H, m * 6 * 6)

  expect_equal(gh$G, true_gradient)
  expect_equal(gh$H, true_hessian)
})

test_that("Gradient (2)", {
  x <- lltd$stats_1[, 1]
  theta <- lltd$theta_6

  m <- length(x)

  true_gradient <- matrix(
    c(
      # alpha
      rep(1, m),
      # delta
      0, 0.61478815295126437, 0.70710678118654752, 0.73404080682229425,
      0.74478197898796472, 0.75002929973635827, 0.75973440798381158,
      0.75983467254072205,
      # log_eta
      0, 0, -0.061266133966784199, -0.052027500048956546, -0.039711040682104937,
      -0.030560141529926446, -0.00079213712208102448, -0.000012592132067118430,
      # log_phi
      0, 0.17565375798607553, 0.088388347648318441, 0.047357471407889951,
      0.028645460730306335, 0.018988083537629323, 0.00020248784860975788,
      2.0262149869619950e-06,
      # log_nu
      0, -0.21125374759367113, -0.20087036204297758, -0.20327982346606966,
      -0.20513751981469271, -0.20624664127023169, -0.20866341617733304,
      -0.20869092715503820,
      # log_xi
      0, 0.065870159244778325, 0.13258252147247766, 0.15983146600162859,
      0.17187276438183801, 0.17801328316527491, 0.18983235807164802,
      0.18995765502768703
    ),
    nrow = m,
    ncol = 6
  )

  G <- loglogistic6_gradient_2(x, theta)

  expect_type(G, "double")
  expect_length(G, m * 6)
  expect_equal(G, true_gradient)
})

test_that("Hessian (2)", {
  x <- lltd$stats_1[, 1]
  theta <- lltd$theta_6

  m <- length(x)

  true_hessian <- array(
    c(
      # (alpha, alpha)
      rep(0, m),
      # (alpha, delta)
      rep(0, m),
      # (alpha, log_eta)
      rep(0, m),
      # (alpha, log_phi)
      rep(0, m),
      # (alpha, log_nu)
      rep(0, m),
      # (alpha, log_xi)
      rep(0, m),
      # (delta, alpha)
      rep(0, m),
      # (delta, delta)
      rep(0, m),
      # (delta, log_eta)
      0, 0, 0.061266133966784199, 0.052027500048956546, 0.039711040682104937,
      0.030560141529926446, 0.00079213712208102448, 0.000012592132067118430,
      # (delta, log_phi)
      0, -0.17565375798607553, -0.088388347648318441, -0.047357471407889951,
      -0.028645460730306335, -0.018988083537629323, -0.00020248784860975788,
      -2.0262149869619950e-06,
      # (delta, log_nu)
      0, 0.21125374759367113, 0.20087036204297758, 0.20327982346606966,
      0.20513751981469271, 0.20624664127023169, 0.20866341617733304,
      0.20869092715503820,
      # (delta, log_xi)
      0, -0.065870159244778325, -0.13258252147247766, -0.15983146600162859,
      -0.17187276438183801, -0.17801328316527491, -0.18983235807164802,
      -0.18995765502768703,
      # (log_eta, alpha)
      rep(0, m),
      # (log_eta, delta)
      0, 0, 0.061266133966784199, 0.052027500048956546, 0.039711040682104937,
      0.030560141529926446, 0.00079213712208102448, 0.000012592132067118430,
      # (log_eta, log_eta)
      0, 0, -0.0028747679353181964, 0.043850520819693315, 0.059804575213080389,
      0.061583254132948111, 0.0054014505589196559, 0.00014391715637766601,
      # (log_eta, log_phi)
      0, 0.17565375798607553, 0.0041474134439901663, -0.039914464158101675,
      -0.043139882041191050, -0.038263827176663258, -0.0013807307757200936,
      -0.000023157881253076072,
      # (log_eta, log_nu)
      0, 0, -0.0020875571468221206, -0.0076948846451738715,
      -0.0078830348311160946, -0.0068562248281684448, -0.00021714067436908112,
      -3.4584005997322730e-06,
      # (log_eta, log_xi)
      0, 0, 0.057437000593860187, 0.056642842795234950, 0.045820431556274927,
      0.036265990739627902, 0.00098964359268518184, 0.000015740081136798641,
      # (log_phi, alpha)
      rep(0, m),
      # (log_phi, delta)
      0, -0.17565375798607553, -0.088388347648318441, -0.047357471407889951,
      -0.028645460730306335, -0.018988083537629323, -0.00020248784860975788,
      -2.0262149869619950e-06,
      # (log_phi, log_eta)
      0, 0.17565375798607553, 0.0041474134439901663, -0.039914464158101675,
      -0.043139882041191050, -0.038263827176663258, -0.0013807307757200936,
      -0.000023157881253076072,
      # (log_phi, log_phi)
      0, 0.10037357599204316, 0.12153397801643786, 0.079438339135815403,
      0.051782179012476837, 0.035572612197077719, 0.00040470585733597398,
      4.0524029578682495e-06,
      # (log_phi, log_nu)
      0, -0.040015362393851410, 0.0030117083432925871, 0.0070041858484057909,
      0.0056864076290030784, 0.0042600120049360264, 0.000055505980938196568,
      5.5649536462489340e-07,
      # (log_phi, log_xi)
      0, -0.094100227492540464, -0.082864075920298538, -0.051558537419880189,
      -0.033052454688815002, -0.022533326982946191, -0.00025297489082042646,
      -2.5327552256746235e-06,
      # (log_nu, alpha)
      rep(0, m),
      # (log_nu, delta)
      0, 0.21125374759367113, 0.20087036204297758, 0.20327982346606966,
      0.20513751981469271, 0.20624664127023169, 0.20866341617733304,
      0.20869092715503820,
      # (log_nu, log_eta)
      0, 0, -0.0020875571468221206, -0.0076948846451738715,
      -0.0078830348311160946, -0.0068562248281684448, -0.00021714067436908112,
      -3.4584005997322730e-06,
      # (log_nu, log_phi)
      0, -0.040015362393851410, 0.0030117083432925871, 0.0070041858484057909,
      0.0056864076290030784, 0.0042600120049360264, 0.000055505980938196568,
      5.5649536462489340e-07,
      # (log_nu, log_nu)
      0, 0.088475863529621896, 0.13275985366055119, 0.14392969426120596,
      0.14753413296251400, 0.14905124321476543, 0.15135330649351379,
      0.15137332287237641,
      # (log_nu, log_xi)
      0, -0.080875920142472603, -0.12806495895753878, -0.13619233876325904,
      -0.13775431860781954, -0.13807567061899966, -0.13779550094208873,
      -0.13778621459410327,
      # (log_xi, alpha)
      rep(0, m),
      # (log_xi, delta)
      0, -0.065870159244778325, -0.13258252147247766, -0.15983146600162859,
      -0.17187276438183801, -0.17801328316527491, -0.18983235807164802,
      -0.18995765502768703,
      # (log_xi, log_eta)
      0, 0, 0.057437000593860187, 0.056642842795234950, 0.045820431556274927,
      0.036265990739627902, 0.00098964359268518184, 0.000015740081136798641,
      # (log_xi, log_phi)
      0, -0.094100227492540464, -0.082864075920298538, -0.051558537419880189,
      -0.033052454688815002, -0.022533326982946191, -0.00025297489082042646,
      -2.5327552256746235e-06,
      # (log_xi, log_nu)
      0, -0.080875920142472603, -0.12806495895753878, -0.13619233876325904,
      -0.13775431860781954, -0.13807567061899966, -0.13779550094208873,
      -0.13778621459410327,
      # (log_xi, log_xi)
      0, 0.030582573935075651, 0.0082864075920298538, -0.014178597790467052,
      -0.026441963751052002, -0.033236657299845631, -0.047331602072501791,
      -0.047488147379308920
    ),
    dim = c(m, 6, 6)
  )

  H <- loglogistic6_hessian_2(x, theta)

  expect_type(H, "double")
  expect_length(H, m * 6 * 6)
  expect_equal(H, true_hessian)
})

test_that("Gradient and Hessian (2)", {
  x <- lltd$stats_1[, 1]
  theta <- lltd$theta_6

  m <- length(x)

  true_gradient <- matrix(
    c(
      # alpha
      rep(1, m),
      # delta
      0, 0.61478815295126437, 0.70710678118654752, 0.73404080682229425,
      0.74478197898796472, 0.75002929973635827, 0.75973440798381158,
      0.75983467254072205,
      # log_eta
      0, 0, -0.061266133966784199, -0.052027500048956546, -0.039711040682104937,
      -0.030560141529926446, -0.00079213712208102448, -0.000012592132067118430,
      # log_phi
      0, 0.17565375798607553, 0.088388347648318441, 0.047357471407889951,
      0.028645460730306335, 0.018988083537629323, 0.00020248784860975788,
      2.0262149869619950e-06,
      # log_nu
      0, -0.21125374759367113, -0.20087036204297758, -0.20327982346606966,
      -0.20513751981469271, -0.20624664127023169, -0.20866341617733304,
      -0.20869092715503820,
      # log_xi
      0, 0.065870159244778325, 0.13258252147247766, 0.15983146600162859,
      0.17187276438183801, 0.17801328316527491, 0.18983235807164802,
      0.18995765502768703
    ),
    nrow = m,
    ncol = 6
  )

  true_hessian <- array(
    c(
      # (alpha, alpha)
      rep(0, m),
      # (alpha, delta)
      rep(0, m),
      # (alpha, log_eta)
      rep(0, m),
      # (alpha, log_phi)
      rep(0, m),
      # (alpha, log_nu)
      rep(0, m),
      # (alpha, log_xi)
      rep(0, m),
      # (delta, alpha)
      rep(0, m),
      # (delta, delta)
      rep(0, m),
      # (delta, log_eta)
      0, 0, 0.061266133966784199, 0.052027500048956546, 0.039711040682104937,
      0.030560141529926446, 0.00079213712208102448, 0.000012592132067118430,
      # (delta, log_phi)
      0, -0.17565375798607553, -0.088388347648318441, -0.047357471407889951,
      -0.028645460730306335, -0.018988083537629323, -0.00020248784860975788,
      -2.0262149869619950e-06,
      # (delta, log_nu)
      0, 0.21125374759367113, 0.20087036204297758, 0.20327982346606966,
      0.20513751981469271, 0.20624664127023169, 0.20866341617733304,
      0.20869092715503820,
      # (delta, log_xi)
      0, -0.065870159244778325, -0.13258252147247766, -0.15983146600162859,
      -0.17187276438183801, -0.17801328316527491, -0.18983235807164802,
      -0.18995765502768703,
      # (log_eta, alpha)
      rep(0, m),
      # (log_eta, delta)
      0, 0, 0.061266133966784199, 0.052027500048956546, 0.039711040682104937,
      0.030560141529926446, 0.00079213712208102448, 0.000012592132067118430,
      # (log_eta, log_eta)
      0, 0, -0.0028747679353181964, 0.043850520819693315, 0.059804575213080389,
      0.061583254132948111, 0.0054014505589196559, 0.00014391715637766601,
      # (log_eta, log_phi)
      0, 0.17565375798607553, 0.0041474134439901663, -0.039914464158101675,
      -0.043139882041191050, -0.038263827176663258, -0.0013807307757200936,
      -0.000023157881253076072,
      # (log_eta, log_nu)
      0, 0, -0.0020875571468221206, -0.0076948846451738715,
      -0.0078830348311160946, -0.0068562248281684448, -0.00021714067436908112,
      -3.4584005997322730e-06,
      # (log_eta, log_xi)
      0, 0, 0.057437000593860187, 0.056642842795234950, 0.045820431556274927,
      0.036265990739627902, 0.00098964359268518184, 0.000015740081136798641,
      # (log_phi, alpha)
      rep(0, m),
      # (log_phi, delta)
      0, -0.17565375798607553, -0.088388347648318441, -0.047357471407889951,
      -0.028645460730306335, -0.018988083537629323, -0.00020248784860975788,
      -2.0262149869619950e-06,
      # (log_phi, log_eta)
      0, 0.17565375798607553, 0.0041474134439901663, -0.039914464158101675,
      -0.043139882041191050, -0.038263827176663258, -0.0013807307757200936,
      -0.000023157881253076072,
      # (log_phi, log_phi)
      0, 0.10037357599204316, 0.12153397801643786, 0.079438339135815403,
      0.051782179012476837, 0.035572612197077719, 0.00040470585733597398,
      4.0524029578682495e-06,
      # (log_phi, log_nu)
      0, -0.040015362393851410, 0.0030117083432925871, 0.0070041858484057909,
      0.0056864076290030784, 0.0042600120049360264, 0.000055505980938196568,
      5.5649536462489340e-07,
      # (log_phi, log_xi)
      0, -0.094100227492540464, -0.082864075920298538, -0.051558537419880189,
      -0.033052454688815002, -0.022533326982946191, -0.00025297489082042646,
      -2.5327552256746235e-06,
      # (log_nu, alpha)
      rep(0, m),
      # (log_nu, delta)
      0, 0.21125374759367113, 0.20087036204297758, 0.20327982346606966,
      0.20513751981469271, 0.20624664127023169, 0.20866341617733304,
      0.20869092715503820,
      # (log_nu, log_eta)
      0, 0, -0.0020875571468221206, -0.0076948846451738715,
      -0.0078830348311160946, -0.0068562248281684448, -0.00021714067436908112,
      -3.4584005997322730e-06,
      # (log_nu, log_phi)
      0, -0.040015362393851410, 0.0030117083432925871, 0.0070041858484057909,
      0.0056864076290030784, 0.0042600120049360264, 0.000055505980938196568,
      5.5649536462489340e-07,
      # (log_nu, log_nu)
      0, 0.088475863529621896, 0.13275985366055119, 0.14392969426120596,
      0.14753413296251400, 0.14905124321476543, 0.15135330649351379,
      0.15137332287237641,
      # (log_nu, log_xi)
      0, -0.080875920142472603, -0.12806495895753878, -0.13619233876325904,
      -0.13775431860781954, -0.13807567061899966, -0.13779550094208873,
      -0.13778621459410327,
      # (log_xi, alpha)
      rep(0, m),
      # (log_xi, delta)
      0, -0.065870159244778325, -0.13258252147247766, -0.15983146600162859,
      -0.17187276438183801, -0.17801328316527491, -0.18983235807164802,
      -0.18995765502768703,
      # (log_xi, log_eta)
      0, 0, 0.057437000593860187, 0.056642842795234950, 0.045820431556274927,
      0.036265990739627902, 0.00098964359268518184, 0.000015740081136798641,
      # (log_xi, log_phi)
      0, -0.094100227492540464, -0.082864075920298538, -0.051558537419880189,
      -0.033052454688815002, -0.022533326982946191, -0.00025297489082042646,
      -2.5327552256746235e-06,
      # (log_xi, log_nu)
      0, -0.080875920142472603, -0.12806495895753878, -0.13619233876325904,
      -0.13775431860781954, -0.13807567061899966, -0.13779550094208873,
      -0.13778621459410327,
      # (log_xi, log_xi)
      0, 0.030582573935075651, 0.0082864075920298538, -0.014178597790467052,
      -0.026441963751052002, -0.033236657299845631, -0.047331602072501791,
      -0.047488147379308920
    ),
    dim = c(m, 6, 6)
  )

  gh <- loglogistic6_gradient_hessian_2(x, theta)

  expect_type(gh, "list")
  expect_type(gh$G, "double")
  expect_type(gh$H, "double")

  expect_length(gh$G, m * 6)
  expect_length(gh$H, m * 6 * 6)

  expect_equal(gh$G, true_gradient)
  expect_equal(gh$H, true_hessian)

  object <- structure(list(stats = lltd$stats_1), class = "loglogistic6")

  gh <- gradient_hessian(object, theta)

  expect_type(gh, "list")
  expect_type(gh$G, "double")
  expect_type(gh$H, "double")

  expect_length(gh$G, m * 6)
  expect_length(gh$H, m * 6 * 6)

  expect_equal(gh$G, true_gradient)
  expect_equal(gh$H, true_hessian)
})

test_that("Value of the RSS", {
  theta <- lltd$theta_6
  theta[3:6] <- log(theta[3:6])

  true_value <- 2.8517831854811524

  object <- structure(
    list(stats = lltd$stats_1, m = nrow(lltd$stats_1)),
    class = "loglogistic6"
  )

  rss_fn <- rss(object)

  expect_type(rss_fn, "closure")

  value <- rss_fn(theta)

  expect_type(value, "double")
  expect_length(value, 1)
  expect_equal(value, true_value)

  known_param <- c(theta[1], NA, NA, theta[4], theta[5], NA)
  rss_fn <- rss_fixed(object, known_param)

  expect_type(rss_fn, "closure")

  value <- rss_fn(theta[c(2, 3, 6)])

  expect_type(value, "double")
  expect_length(value, 1)
  expect_equal(value, true_value)
})

test_that("Gradient and Hessian of the RSS", {
  theta <- lltd$theta_6
  theta[3:6] <- log(theta[3:6])

  true_gradient <- c(
    -6.0009199747246628, -4.0522246811859619, 0.19326331591384029,
    -0.54349897146732128, 1.2055648493511092, -0.74130668456282984
  )

  true_hessian <- matrix(
    c(
      # alpha
      19, 11.458419973724663, -0.48193239808186501, 0.99849932252717809,
      -3.2962982209552613, 2.3653553321675766,
      # delta
      11.458419973724663, 8.2475430534317714, -0.54553034548108759,
      1.2100724872004481, -3.5647974500569246, 2.4699056900542093,
      # log_eta
      -0.48193239808186501, -0.54553034548108759, -0.10640462804110470,
      -0.25893162252936833, 0.12101531495410192, -0.27866626083924837,
      # log_phi
      0.99849932252717809, 1.2100724872004481, -0.25893162252936833,
      -0.39071424549993905, -0.15025685364370334, 0.49787524455008154,
      # log_nu
      -3.2962982209552613, -3.5647974500569246, 0.12101531495410192,
      -0.15025685364370334, -0.048738529350026712, 0.19698523460411099,
      # log_xi
      2.3653553321675766, 2.4699056900542093, -0.27866626083924837,
      0.49787524455008154, 0.19698523460411099, 0.36853880023851156
    ),
    nrow = 6,
    ncol = 6
  )

  object <- structure(
    list(stats = lltd$stats_1, m = nrow(lltd$stats_1)),
    class = "loglogistic6"
  )

  rss_gh <- rss_gradient_hessian(object)

  expect_type(rss_gh, "closure")

  gh <- rss_gh(theta)

  expect_type(gh$G, "double")
  expect_type(gh$H, "double")

  expect_length(gh$G, 6)
  expect_length(gh$H, 6 * 6)

  expect_equal(gh$G, true_gradient)
  expect_equal(gh$H, true_hessian)

  known_param <- c(theta[1], NA, NA, theta[4], theta[5], NA)
  rss_gh <- rss_gradient_hessian_fixed(object, known_param)

  expect_type(rss_gh, "closure")

  gh <- rss_gh(theta[c(2, 3, 6)])

  expect_type(gh$G, "double")
  expect_type(gh$H, "double")

  expect_length(gh$G, 3)
  expect_length(gh$H, 3 * 3)

  expect_equal(gh$G, true_gradient[c(2, 3, 6)])
  expect_equal(gh$H, true_hessian[c(2, 3, 6), c(2, 3, 6)])
})

test_that("mle_asy", {
  x <- lltd$D$x
  y <- lltd$D$y
  w <- rep(1, length(y))

  max_iter <- 10000

  theta <- c(
    0, 1, 3.3501562135542870, 2.3530975491734142, 3.0041670187744469,
    3.4109688853855420
  )

  true_value <- c(
    0.86567490772243801, -0.91541486214812162, 3.3501562135542870,
    2.3530975491734142, 3.0041670187744469, 3.4109688853855420
  )

  object <- loglogistic6_new(x, y, w, NULL, max_iter, NULL, NULL)

  result <- mle_asy(object, theta)

  expect_type(result, "double")
  expect_length(result, 6)
  expect_equal(result, true_value)
})

test_that("fit", {
  x <- lltd$D$x
  y <- lltd$D$y

  n <- length(y)
  w <- rep(1, n)

  k <- as.numeric(table(x))

  max_iter <- 10000

  # loglogistic6 model is basically unidentifiable: many parameters are
  # associated with the same residual sum of squares
  # there is no point in testing the values of `result$coefficients`
  estimated <- c(
    alpha = TRUE, delta = TRUE, eta = TRUE, phi = TRUE, nu = TRUE, xi = TRUE
  )

  rss_value <- 0.058276476180507351

  fitted_values <- rep(
    c(
      0.86567490772243801, 0.7901639477115967, 0.6645425346486514,
      0.508922320838699, 0.329951365685670, 0.142289339645596,
      0.092688548125015, 0.092688548125015
    ),
    k
  )

  residuals <- c(
    -0.01297490772243801, -0.10857490772243801, 0.07132509227756199,
    0.0373360522884033, -0.0123639477115967, 0.0944360522884033,
    -0.1084425346486514, 0.0402574653513486, 0.036377679161301,
    -0.022922320838699, 0.026177679161301, -0.069022320838699,
    0.025348634314330, -0.016751365685670, 0.019948634314330,
    -0.000189339645596, -0.075888548125015, 0.045911451874985, 0.030011451874985
  )

  object <- loglogistic6_new(x, y, w, NULL, max_iter, NULL, NULL)

  result <- fit(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_false(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 6)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  object <- loglogistic6_new(x, y, w, c(0, 1, 1, 1, 1, 1), max_iter, NULL, NULL)

  result <- fit(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_false(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 6)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)
})

test_that("fit_constrained: inequalities", {
  x <- lltd$D$x
  y <- lltd$D$y

  n <- length(y)
  w <- rep(1, n)

  k <- as.numeric(table(x))

  max_iter <- 10000

  estimated <- c(
    alpha = TRUE, delta = TRUE, eta = TRUE, phi = TRUE, nu = TRUE, xi = TRUE
  )

  rss_value <- 0.058276476180507351

  fitted_values <- rep(
    c(
      0.86567490772243801, 0.7901639477115967, 0.6645425346486514,
      0.508922320838699, 0.329951365685670, 0.142289339645596,
      0.092688548125015, 0.092688548125015
    ),
    k
  )

  residuals <- c(
    -0.01297490772243801, -0.10857490772243801, 0.07132509227756199,
    0.0373360522884033, -0.0123639477115967, 0.0944360522884033,
    -0.1084425346486514, 0.0402574653513486, 0.036377679161301,
    -0.022922320838699, 0.026177679161301, -0.069022320838699,
    0.025348634314330, -0.016751365685670, 0.019948634314330,
    -0.000189339645596, -0.075888548125015, 0.045911451874985, 0.030011451874985
  )

  object <- loglogistic6_new(
    x, y, w, NULL, max_iter,
    c(0.5, -1, 25, 8, 15, 20), c(1, -0.5, 30, 12, 30, 40)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_true(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 6)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  # initial values within the boundaries
  object <- loglogistic6_new(
    x, y, w, c(0.7, -0.6, 29, 11, 16, 38), max_iter,
    c(0.5, -1, 25, 8, 15, 20), c(1, -0.5, 30, 12, 30, 40)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_true(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 6)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  # initial values outside the boundaries
  object <- loglogistic6_new(
    x, y, w, c(-2, 2, 1, 1, 1, 1), max_iter,
    c(0.5, -1, 25, 8, 15, 20), c(1, -0.5, 30, 12, 30, 40)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_true(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 6)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)
})

test_that("fit_constrained: equalities", {
  x <- lltd$D$x
  y <- lltd$D$y

  n <- length(y)
  w <- rep(1, n)

  k <- as.numeric(table(x))

  max_iter <- 10000

  estimated <- c(
    alpha = FALSE, delta = FALSE, eta = TRUE, phi = TRUE, nu = TRUE, xi = TRUE
  )

  rss_value <- 0.074747066165150929

  fitted_values <- rep(
    c(
      0.8, 0.77225113247950056, 0.6741844832087593, 0.5057443762603047,
      0.317532596086360, 0.190181348756377, 0.088632229933508, 0.088631934629668
    ),
    k
  )

  residuals <- c(
    0.0527, -0.0429, 0.1370, 0.05524886752049944, 0.00554886752049944,
    0.11234886752049944, -0.1180844832087593, 0.0306155167912407,
    0.0395556237396953, -0.0197443762603047, 0.0293556237396953,
    -0.0658443762603047, 0.037767403913640, -0.004332596086360,
    0.032367403913640, -0.048081348756377, -0.071832229933508,
    0.049967770066492, 0.034068065370332
  )

  object <- loglogistic6_new(
    x, y, w, NULL, max_iter,
    c(0.8, -0.9, rep(-Inf, 4)), c(0.8, -0.9, rep(Inf, 4))
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_false(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 4)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  # initial values with same equalities
  object <- loglogistic6_new(
    x, y, w, c(0.8, -0.9, 1, 1, 1, 1), max_iter,
    c(0.8, -0.9, rep(-Inf, 4)), c(0.8, -0.9, rep(Inf, 4))
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_false(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 4)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  # initial values with different equalities
  object <- loglogistic6_new(
    x, y, w, c(0, 1, 1, 1, 1, 1), max_iter,
    c(0.8, -0.9, rep(-Inf, 4)), c(0.8, -0.9, rep(Inf, 4))
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_false(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 4)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)
})

test_that("fit_constrained: equalities and inequalities", {
  x <- lltd$D$x
  y <- lltd$D$y

  n <- length(y)
  w <- rep(1, n)

  k <- as.numeric(table(x))

  max_iter <- 10000

  estimated <- c(
    alpha = FALSE, delta = FALSE, eta = TRUE, phi = TRUE, nu = TRUE, xi = TRUE
  )

  rss_value <- 0.074747066165150929

  fitted_values <- rep(
    c(
      0.8, 0.77225113247950056, 0.6741844832087593, 0.5057443762603047,
      0.317532596086360, 0.190181348756377, 0.088632229933508, 0.088631934629668
    ),
    k
  )

  residuals <- c(
    0.0527, -0.0429, 0.1370, 0.05524886752049944, 0.00554886752049944,
    0.11234886752049944, -0.1180844832087593, 0.0306155167912407,
    0.0395556237396953, -0.0197443762603047, 0.0293556237396953,
    -0.0658443762603047, 0.037767403913640, -0.004332596086360,
    0.032367403913640, -0.048081348756377, -0.071832229933508,
    0.049967770066492, 0.034068065370332
  )

  object <- loglogistic6_new(
    x, y, w, NULL, max_iter,
    c(0.8, -0.9, 4, 5, 2, 1), c(0.8, -0.9, 8, 10, 3, 3)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_true(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 4)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  # initial values within the boundaries
  object <- loglogistic6_new(
    x, y, w, c(0.8, -0.9, 7, 9, 2.1, 1.2), max_iter,
    c(0.8, -0.9, 4, 5, 2, 1), c(0.8, -0.9, 8, 10, 3, 3)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_true(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 4)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  # initial values outside the boundaries
  object <- loglogistic6_new(
    x, y, w, c(0, 1, 0.5, 0.5, 10, 5), max_iter,
    c(0.8, -0.9, 4, 5, 2, 1), c(0.8, -0.9, 8, 10, 3, 3)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_true(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 4)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)
})

test_that("fit (weighted)", {
  x <- lltd$D$x
  y <- lltd$D$y
  w <- lltd$D$w

  n <- length(y)

  k <- as.numeric(table(x))

  max_iter <- 10000

  # loglogistic6 model is basically unidentifiable: many parameters are
  # associated with the same residual sum of squares
  # there is no point in testing the values of `result$coefficients`
  estimated <- c(
    alpha = TRUE, delta = TRUE, eta = TRUE, phi = TRUE, nu = TRUE, xi = TRUE
  )

  rss_value <- 0.022085036915158753

  fitted_values <- rep(
    c(
      0.90681779366690136, 0.82534227701819370, 0.68768283546185108,
      0.51593867110192256, 0.31929217045238740, 0.14655690461413620,
      0.093029886494351149, 0.093029886494351068
    ),
    k
  )

  residuals <- c(
    -0.054117793666901360, -0.14971779366690136, 0.030182206333098640,
    0.0021577229818062970, -0.047542277018193703, 0.059257722981806297,
    -0.13158283546185108, 0.017117164538148918, 0.029361328898077444,
    -0.029938671101922556, 0.019161328898077444, -0.076038671101922556,
    0.036007829547612596, -0.0060921704523874042, 0.030607829547612596,
    -0.0044569046141361993, -0.076229886494351149, 0.045570113505648851,
    0.029670113505648932
  )

  object <- loglogistic6_new(x, y, w, NULL, max_iter, NULL, NULL)

  result <- fit(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_false(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 6)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  object <- loglogistic6_new(
    x, y, w, c(1, -1, 1, 1, 1, 1), max_iter, NULL, NULL
  )

  result <- fit(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_false(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 6)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)
})

test_that("fit_constrained (weighted): inequalities", {
  x <- lltd$D$x
  y <- lltd$D$y
  w <- lltd$D$w

  n <- length(y)

  k <- as.numeric(table(x))

  max_iter <- 10000

  estimated <- c(
    alpha = TRUE, delta = TRUE, eta = TRUE, phi = TRUE, nu = TRUE, xi = TRUE
  )

  rss_value <- 0.022085036915158753

  fitted_values <- rep(
    c(
      0.90681779366690136, 0.82534227701819370, 0.68768283546185108,
      0.51593867110192256, 0.31929217045238740, 0.14655690461413620,
      0.093029886494351149, 0.093029886494351068
    ),
    k
  )

  residuals <- c(
    -0.054117793666901360, -0.14971779366690136, 0.030182206333098640,
    0.0021577229818062970, -0.047542277018193703, 0.059257722981806297,
    -0.13158283546185108, 0.017117164538148918, 0.029361328898077444,
    -0.029938671101922556, 0.019161328898077444, -0.076038671101922556,
    0.036007829547612596, -0.0060921704523874042, 0.030607829547612596,
    -0.0044569046141361993, -0.076229886494351149, 0.045570113505648851,
    0.029670113505648932
  )

  object <- loglogistic6_new(
    x, y, w, NULL, max_iter,
    c(0.5, -1, 10, 8, 5, 0), c(1, -0.5, 20, 10, 15, 2)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_true(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 6)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  # initial values within the boundaries
  object <- loglogistic6_new(
    x, y, w, c(0.7, -0.6, 18, 8.5, 7, 1.8), max_iter,
    c(0.5, -1, 10, 8, 5, 0), c(1, -0.5, 20, 10, 15, 2)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_true(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 6)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  # initial values outside the boundaries
  object <- loglogistic6_new(
    x, y, w, c(-2, -5, 0.5, 20, 0.1, 3), 10000,
    c(0.5, -1, 10, 8, 5, 0), c(1, -0.5, 20, 10, 15, 2)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_true(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 6)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)
})

test_that("fit_constrained (weighted): equalities", {
  x <- lltd$D$x
  y <- lltd$D$y
  w <- lltd$D$w

  n <- length(y)

  k <- as.numeric(table(x))

  max_iter <- 10000

  estimated <- c(
    alpha = FALSE, delta = FALSE, eta = TRUE, phi = TRUE, nu = TRUE, xi = TRUE
  )

  rss_value <- 0.022165226781174076

  fitted_values <- rep(
    c(
      0.9, 0.82394807262993030, 0.68895716760665773, 0.51663843729538372,
      0.31768368751798345, 0.14888199437111868, 0.092568060779228249,
      0.092568060779225243
    ),
    k
  )

  residuals <- c(
    -0.0473, -0.1429, 0.037, 0.0035519273700696982, -0.046148072629930302,
    0.060651927370069698, -0.13285716760665773, 0.015842832393342270,
    0.028661562704616279, -0.030638437295383721, 0.018461562704616279,
    -0.076738437295383721, 0.037616312482016546, -0.0044836875179834541,
    0.032216312482016546, -0.0067819943711186764, -0.075768060779228249,
    0.046031939220771751, 0.030131939220774757
  )

  object <- loglogistic6_new(
    x, y, w, NULL, max_iter,
    c(0.9, -0.9, rep(-Inf, 4)), c(0.9, -0.9, rep(Inf, 4))
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_false(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 4)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  # initial values with same equalities
  object <- loglogistic6_new(
    x, y, w, c(0.9, -0.9, 1, 1, 1, 1), max_iter,
    c(0.9, -0.9, rep(-Inf, 4)), c(0.9, -0.9, rep(Inf, 4))
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_false(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 4)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  # initial values with different equalities
  object <- loglogistic6_new(
    x, y, w, c(0, 1, 1, 1, 1, 1), max_iter,
    c(0.9, -0.9, rep(-Inf, 4)), c(0.9, -0.9, rep(Inf, 4))
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_false(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 4)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)
})

test_that("fit_constrained (weighted): equalities and inequalities", {
  x <- lltd$D$x
  y <- lltd$D$y
  w <- lltd$D$w

  n <- length(y)

  k <- as.numeric(table(x))

  max_iter <- 10000

  estimated <- c(
    alpha = FALSE, delta = FALSE, eta = TRUE, phi = TRUE, nu = TRUE, xi = TRUE
  )

  rss_value <- 0.022165226781174076

  fitted_values <- rep(
    c(
      0.9, 0.82394807262993030, 0.68895716760665773, 0.51663843729538372,
      0.31768368751798345, 0.14888199437111868, 0.092568060779228249,
      0.092568060779225243
    ),
    k
  )

  residuals <- c(
    -0.0473, -0.1429, 0.037, 0.0035519273700696982, -0.046148072629930302,
    0.060651927370069698, -0.13285716760665773, 0.015842832393342270,
    0.028661562704616279, -0.030638437295383721, 0.018461562704616279,
    -0.076738437295383721, 0.037616312482016546, -0.0044836875179834541,
    0.032216312482016546, -0.0067819943711186764, -0.075768060779228249,
    0.046031939220771751, 0.030131939220774757
  )

  object <- loglogistic6_new(
    x, y, w, NULL, max_iter,
    c(0.9, -0.9, 5, 5, 7, 0), c(0.9, -0.9, 15, 10, 12, 3)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_true(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 4)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  # initial values within the boundaries
  object <- loglogistic6_new(
    x, y, w, c(0.9, -0.9, 6, 7, 11, 0.5), max_iter,
    c(0.9, -0.9, 5, 5, 7, 0), c(0.9, -0.9, 15, 10, 12, 3)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_true(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 4)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)

  # initial values outside the boundaries
  object <- loglogistic6_new(
    x, y, w, c(0, 1, 0.5, 0.5, 5, 6), max_iter,
    c(0.9, -0.9, 5, 5, 7, 0), c(0.9, -0.9, 15, 10, 12, 3)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "loglogistic6_fit"))
  expect_true(inherits(result, "loglogistic"))
  expect_true(result$converged)
  expect_true(result$constrained)
  expect_equal(result$estimated, estimated)
  expect_equal(result$rss, rss_value)
  expect_equal(result$df.residual, object$n - 4)
  expect_equal(result$fitted.values, fitted_values, tolerance = 1.0e-6)
  expect_equal(result$residuals, residuals, tolerance = 1.0e-6)
  expect_equal(result$weights, w)
})

test_that("fisher_info", {
  x <- lltd$D$x
  y <- lltd$D$y
  w <- lltd$D$w

  max_iter <- 10000

  theta <- lltd$theta_6
  names(theta) <- c("alpha", "delta", "eta", "phi", "nu", "xi")

  sigma <- lltd$sigma

  true_value <- matrix(c(
      # alpha
      6206.96, 4033.6259327887426, -78.699146035632538, 198.77058917560704,
      -293.58875257662276, 269.87863134052620, 89181.461311549703,
      # delta
      4033.6259327887426, 2878.4597528145919, -88.705348162354274,
      244.00408263047195, -320.80516117638226, 282.55798267636782,
      59530.094792948388,
      # eta
      -78.699146035632538, -88.705348162354274, -24.139034413674205,
      -33.757408298046454, 4.9250747626882528, -15.028191756557660,
      -1242.7147074919567,
      # phi
      198.77058917560704, 244.00408263047195, -33.757408298046454,
      21.182660748584033, -6.5834237992925952, 32.161335478721606,
      4533.2517221704692,
      # nu
      -293.58875257662276, -320.80516117638226, 4.9250747626882528,
      -6.5834237992925952, -29.341061059068340, 6.4659688730596288,
      -4510.5256202602449,
      # xi
      269.87863134052620, 282.55798267636782, -15.028191756557660,
      32.161335478721606, 6.4659688730596288, 41.574031108086305,
      3449.7573253555426,
      # sigma
      89181.461311549703, 59530.094792948388, -1242.7147074919567,
      4533.2517221704692, -4510.5256202602449, 3449.7573253555426,
      1.3580467956656330e+06
    ),
    nrow = 7,
    ncol = 7
  )

  rownames(true_value) <- colnames(true_value) <- c(
    "alpha", "delta", "eta", "phi", "nu", "xi", "sigma"
  )

  object <- loglogistic6_new(x, y, w, NULL, max_iter, NULL, NULL)

  fim <- fisher_info(object, theta, sigma)

  expect_type(fim, "double")
  expect_length(fim, 7 * 7)
  expect_equal(fim, true_value)
})

test_that("drda: 'lower_bound' argument errors", {
  x <- lltd$D$x
  y <- lltd$D$y

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      lower_bound = c("a", "b", "c", "d", "e", "f")
    ),
    "'lower_bound' must be a numeric vector"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      lower_bound = matrix(-Inf, nrow = 6, ncol = 2),
      upper_bound = rep(Inf, 6)
    ),
    "'lower_bound' must be a numeric vector"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      lower_bound = rep(-Inf, 7),
      upper_bound = rep(Inf, 6)
    ),
    "'lower_bound' and 'upper_bound' must have the same length"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      lower_bound = c( 0, -Inf, -Inf, -Inf, -Inf, -Inf),
      upper_bound = c(-1, Inf, Inf, Inf, Inf, Inf)
    ),
    "'lower_bound' cannot be larger than 'upper_bound'"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      lower_bound = c(Inf, -Inf, -Inf, -Inf, -Inf, -Inf),
      upper_bound = c(Inf, Inf, Inf, Inf, Inf, Inf)
    ),
    "'lower_bound' cannot be equal to infinity"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      lower_bound = rep(-Inf, 7),
      upper_bound = rep(Inf, 7)
    ),
    "'lower_bound' must be of length 6"
  )
})

test_that("drda: 'upper_bound' argument errors", {
  x <- lltd$D$x
  y <- lltd$D$y

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      upper_bound = c("a", "b", "c", "d", "e", "f")
    ),
    "'upper_bound' must be a numeric vector"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      lower_bound = rep(-Inf, 6),
      upper_bound = matrix(Inf, nrow = 6, ncol = 2)
    ),
    "'upper_bound' must be a numeric vector"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      lower_bound = c(-Inf, -Inf, -Inf, -Inf, -Inf, -Inf),
      upper_bound = c(-Inf, Inf, Inf, Inf, Inf, Inf)
    ),
    "'upper_bound' cannot be equal to -infinity"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      lower_bound = rep(-Inf, 7),
      upper_bound = rep(Inf, 7)
    ),
    "'lower_bound' must be of length 6"
  )
})

test_that("drda: 'start' argument errors", {
  x <- lltd$D$x
  y <- lltd$D$y

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      start = c("a", "b", "c", "d", "e", "f")
    ),
    "'start' must be a numeric vector"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      start = c(0, Inf, 1, 1, 1, 1)
    ),
    "'start' must be finite"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      start = c(-Inf, 1, 1, 1, 1, 1)
    ),
    "'start' must be finite"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      start = c(1, 1, 1, 1, 1, 1, 1)
    ),
    "'start' must be of length 6"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      start = c(0, 1, -1, 1, 1, 1)
    ),
    "parameter 'eta' cannot be negative nor zero"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      start = c(0, 1, 0, 1, 1, 1)
    ),
    "parameter 'eta' cannot be negative nor zero"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      start = c(0, 1, 1, -1, 1, 1)
    ),
    "parameter 'phi' cannot be negative nor zero"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      start = c(0, 1, 1, 0, 1, 1)
    ),
    "parameter 'phi' cannot be negative nor zero"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      start = c(0, 1, 1, 1, -1, 1)
    ),
    "parameter 'nu' cannot be negative nor zero"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      start = c(0, 1, 1, 1, 0, 1)
    ),
    "parameter 'nu' cannot be negative nor zero"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      start = c(0, 1, 1, 1, 1, -1)
    ),
    "parameter 'xi' cannot be negative nor zero"
  )

  expect_error(
    drda(
      y ~ x, mean_function = "loglogistic6",
      start = c(0, 1, 1, 1, 1, 0)
    ),
    "parameter 'xi' cannot be negative nor zero"
  )
})

test_that("nauc: decreasing", {
  x <- lltd$D$x
  y <- lltd$D$y
  w <- lltd$D$w

  result <- drda(y ~ x, weights = w, mean_function = "loglogistic6")

  expect_equal(nauc(result), 0.097909681210381725)
  expect_equal(nauc(result, xlim = c(0, 2)), 0.87325260264665631)
  expect_equal(nauc(result, ylim = c(0.3, 0.7)), 0.0061191350934437442)
  expect_equal(nauc(result, xlim = c(0, 2), ylim = c(0.3, 0.7)), 1.0)
  expect_equal(
    nauc(result, xlim = c(5, 8), ylim = c(0.3, 0.7)), 0.41629510891428020
  )
  expect_equal(nauc(result, xlim = c(10, 15), ylim = c(0.3, 0.7)), 0.0)
})

test_that("naac: decreasing", {
  x <- lltd$D$x
  y <- lltd$D$y
  w <- lltd$D$w

  result <- drda(y ~ x, weights = w, mean_function = "loglogistic6")

  expect_equal(naac(result), 1 - 0.097909681210381725)
  expect_equal(naac(result, xlim = c(0, 2)), 1 - 0.87325260264665631)
  expect_equal(naac(result, ylim = c(0.3, 0.7)), 1 - 0.0061191350934437442)
  expect_equal(naac(result, xlim = c(0, 2), ylim = c(0.3, 0.7)), 0.0)
  expect_equal(
    naac(result, xlim = c(5, 8), ylim = c(0.3, 0.7)), 1 - 0.41629510891428020
  )
  expect_equal(naac(result, xlim = c(10, 15), ylim = c(0.3, 0.7)), 1.0)
})

test_that("nauc: increasing", {
  x <- lltd$D$x
  y <- rev(lltd$D$y)
  w <- lltd$D$w

  result <- drda(y ~ x, weights = w, mean_function = "loglogistic6")

  expect_equal(nauc(result), 0.84987748063128600)
  expect_equal(nauc(result, xlim = c(0, 2)), 0.16335338880478015)
  expect_equal(nauc(result, ylim = c(0.3, 0.7)), 0.99497853656837985)
  expect_equal(nauc(result, xlim = c(0, 2), ylim = c(0.3, 0.7)), 0.0)
  expect_equal(
    nauc(result, xlim = c(5, 8), ylim = c(0.3, 0.7)), 0.79822066734104781
  )
  expect_equal(nauc(result, xlim = c(9, 12), ylim = c(0.3, 0.7)), 1.0)
})

test_that("naac: increasing", {
  x <- lltd$D$x
  y <- rev(lltd$D$y)
  w <- lltd$D$w

  result <- drda(y ~ x, weights = w, mean_function = "loglogistic6")

  expect_equal(naac(result), 1 - 0.84987748063128600)
  expect_equal(naac(result, xlim = c(0, 2)), 1 - 0.16335338880478015)
  expect_equal(naac(result, ylim = c(0.3, 0.7)), 1 - 0.99497853656837985)
  expect_equal(naac(result, xlim = c(0, 2), ylim = c(0.3, 0.7)), 1.0)
  expect_equal(
    naac(result, xlim = c(5, 8), ylim = c(0.3, 0.7)), 1 - 0.79822066734104781
  )
  expect_equal(naac(result, xlim = c(9, 12), ylim = c(0.3, 0.7)), 0.0)
})