test_logistic6.R
In drda: Dose-Response Data Analysis

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

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

  w <- rep(1, n)

  max_iter <- 10000

  stats <- ltd$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)

  i <- c(1, 2, 4)

  s <- start
  s[-i] <- log(s[-i])

  lb <- lower_bound
  lb[-i] <- log(lb[-i])

  ub <- upper_bound
  ub[-i] <- log(ub[-i])

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

  expect_true(inherits(object, "logistic6"))
  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 <- logistic6_new(x, y, w, start, max_iter, lower_bound, upper_bound)

  expect_true(inherits(object, "logistic6"))
  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, s)
  expect_equal(object$lower_bound, lb)
  expect_equal(object$upper_bound, ub)

  w <- ltd$D$w
  stats <- ltd$stats_2

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

  expect_true(inherits(object, "logistic6"))
  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 <- logistic6_new(x, y, w, start, max_iter, lower_bound, upper_bound)

  expect_true(inherits(object, "logistic6"))
  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, s)
  expect_equal(object$lower_bound, lb)
  expect_equal(object$upper_bound, ub)
})

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

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

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

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

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

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

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

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

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

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

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

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

  expect_error(
    logistic6_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(
    logistic6_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(
    logistic6_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(
    logistic6_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(
    logistic6_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(
    logistic6_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 <- ltd$stats_1[, 1]
  theta <- ltd$theta_6

  m <- length(x)

  true_value <- c(
    0.8, 0.79999984704883975, 0.77510647298564850, 0.61621786196251596,
    0.057664189665438460, 0.040170539009997540, 0.040164314348407453,
    0.040164314348407453
  )

  value <- logistic6_fn(x, theta)

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

  object <- structure(list(stats = ltd$stats_1), class = "logistic6")

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

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

  object <- structure(list(stats = ltd$stats_1), class = "logistic6_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 <- ltd$stats_1[, 1]
  theta <- ltd$theta_6

  m <- length(x)

  true_gradient <- matrix(
    c(
      # alpha
      rep(1, m),
      # delta
      9.6437492398195889e-23, 1.5295116025091289e-07, 0.024893527014351497,
      0.18378213803748404, 0.74233581033456154, 0.75982946099000246,
      0.75983568565159255, 0.75983568565159255,
      # eta
      2.4276486539439174e-21, 1.1736381182546699e-06, 0.041653962050552102,
      0.12331460190108316, -0.021582949777827438, -0.000033032686382156039,
      -1.8583072261039617e-42, -1.3923270509654377e-85,
      # phi
      4.8218746199097945e-23, 7.6475580125456447e-08, 0.012446749168025124,
      0.091576578355026419, 0.033030421569699108, 0.000012449068216533443,
      7.5377239464421580e-44, 2.8040905767036438e-87,
      # nu
      -1.2161526556915815e-21, -5.9051181979806446e-07, -0.021428022923202986,
      -0.066384858676451105, -0.051166597326927950, -0.052172373971476468,
      -0.052172801351587403, -0.052172801351587403,
      # xi
      0, 2.0926851264125609e-35, 2.3898584374049245e-09,
      0.000052415110619266708, 0.056356247266263610, 0.063317046904464116,
      0.063319640470966046, 0.063319640470966046
    ),
    nrow = m,
    ncol = 6
  )

  G <- logistic6_gradient(x, theta)

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

test_that("Hessian (1)", {
  x <- ltd$stats_1[, 1]
  theta <- ltd$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)
      -2.4276486539439174e-21, -1.1736381182546699e-06, -0.041653962050552102,
      -0.12331460190108316, 0.021582949777827438, 0.000033032686382156039,
      1.8583072261039617e-42, 1.3923270509654377e-85,
      # (delta, phi)
      -4.8218746199097945e-23, -7.6475580125456447e-08, -0.012446749168025124,
      -0.091576578355026419, -0.033030421569699108, -0.000012449068216533443,
      -7.5377239464421580e-44, -2.8040905767036438e-87,
      # (delta, nu)
      1.2161526556915815e-21, 5.9051181979806446e-07, 0.021428022923202986,
      0.066384858676451105, 0.051166597326927950, 0.052172373971476468,
      0.052172801351587403, 0.052172801351587403,
      # (delta, xi)
      -2.0853048558916222e-111, -2.0926851264125609e-35,
      -2.3898584374049245e-09, -0.000052415110619266708, -0.056356247266263610,
      -0.063317046904464116, -0.063319640470966046, -0.063319640470966046,
      # (eta, alpha)
      rep(0, m),
      # (eta, delta)
      -2.4276486539439174e-21, -1.1736381182546699e-06, -0.041653962050552102,
      -0.12331460190108316, 0.021582949777827438, 0.000033032686382156039,
      1.8583072261039617e-42, 1.3923270509654377e-85,
      # (eta, eta)
      -6.1111895803565778e-20, -9.0056618750765000e-06, -0.069698623184655625,
      -0.081605338544629843, 0.025068184884160287, 0.00017529242457613078,
      9.1627280890809950e-41, 1.3826761752674973e-83,
      # (eta, phi)
      -1.1897149538724097e-21, -5.4858126906460671e-07, -0.014603486473844729,
      -0.014813923084545845, -0.021848994443557445, -0.000059838132613523584,
      -3.6789258324757125e-42, -2.7706336490473572e-85,
      # (eta, nu)
      3.0007647424909336e-20, 4.2377570196669037e-06, 0.025441707980853873,
      0.013819938082184523, -0.0010074674090324482, -2.2678615960535156e-06,
      -1.2759744716990998e-43, -9.5601725502228708e-87,
      # (eta, xi)
      -2.6246988596328932e-109, -8.0288865734436783e-34,
      -1.9994569793277198e-08, -0.00017584811474708742, 0.0081926133528867864,
      0.000013763168318811387, 7.7429467754331736e-43, 5.8013627123559905e-86,
      # (phi, alpha)
      rep(0, m),
      # (phi, delta)
      -4.8218746199097945e-23, -7.6475580125456447e-08, -0.012446749168025124,
      -0.091576578355026419, -0.033030421569699108, -0.000012449068216533443,
      -7.5377239464421580e-44, -2.8040905767036438e-87,
      # (phi, eta)
      -1.1897149538724097e-21, -5.4858126906460671e-07, -0.014603486473844729,
      -0.014813923084545845, -0.021848994443557445, -0.000059838132613523584,
      -3.6789258324757125e-42, -2.7706336490473572e-85,
      # (phi, phi)
      -2.4109373099548972e-23, -3.8237790062728224e-08, -0.0062233387361772991,
      -0.045004753323180023, 0.058712357899407301, 0.000024897116603569502,
      1.5075447892884316e-43, 5.6081811534072877e-87,
      # (phi, nu)
      5.9602164129601628e-22, 2.7613701486766812e-07, 0.0076023154114251624,
      0.010263039600614143, 0.0015418222986489271, 8.5469172528994537e-07,
      5.1756475975932719e-45, 1.9253802288156686e-88,
      # (phi, xi)
      -5.2132621397290555e-111, -5.2317128160314022e-35,
      -5.9746392104900531e-09, -0.00013058930905553107, -0.012537928114042809,
      -5.1869417853060373e-06, -3.1407183110175658e-44, -1.1683710736265183e-87,
      # (nu, alpha)
      rep(0, m),
      # (nu, delta)
      1.2161526556915815e-21, 5.9051181979806446e-07, 0.021428022923202986,
      0.066384858676451105, 0.051166597326927950, 0.052172373971476468,
      0.052172801351587403, 0.052172801351587403,
      # (nu, eta)
      3.0007647424909336e-20, 4.2377570196669037e-06, 0.025441707980853873,
      0.013819938082184523, -0.0010074674090324482, -2.2678615960535156e-06,
      -1.2759744716990998e-43, -9.5601725502228708e-87,
      # (nu, phi)
      5.9602164129601628e-22, 2.7613701486766812e-07, 0.0076023154114251624,
      0.010263039600614143, 0.0015418222986489271, 8.5469172528994537e-07,
      5.1756475975932719e-45, 1.9253802288156686e-88,
      # (nu, nu)
      -1.4730072359242754e-20, -1.9869741882303387e-06, -0.0081199112657906938,
      0.0063612510169811997, 0.021964708947817155, 0.022503861272861439,
      0.022504045632100321, 0.022504045632100321,
      # (nu, xi)
      2.5254682512573550e-110, 7.0330685392331877e-35, 8.6223039892301619e-10,
      -7.2295861790158238e-06, -0.011458416654500857, -0.011482225054956227,
      -0.011482176671775895, -0.011482176671775895,
      # (xi, alpha)
      rep(0, m),
      # (xi, delta)
      -2.0853048558916222e-111, -2.0926851264125609e-35,
      -2.3898584374049245e-09, -0.000052415110619266708, -0.056356247266263610,
      -0.063317046904464116, -0.063319640470966046, -0.063319640470966046,
      # (xi, eta)
      -2.6246988596328932e-109, -8.0288865734436783e-34,
      -1.9994569793277198e-08, -0.00017584811474708742, 0.0081926133528867864,
      0.000013763168318811387, 7.7429467754331736e-43, 5.8013627123559905e-86,
      # (xi, phi)
      -5.2132621397290555e-111, -5.2317128160314022e-35,
      -5.9746392104900531e-09, -0.00013058930905553107, -0.012537928114042809,
      -5.1869417853060373e-06, -3.1407183110175658e-44, -1.1683710736265183e-87,
      # (xi, nu)
      2.5254682512573550e-110, 7.0330685392331877e-35, 8.6223039892301619e-10,
      -7.2295861790158238e-06, -0.011458416654500857, -0.011482225054956227,
      -0.011482176671775895, -0.011482176671775895,
      # (xi, xi)
      -2.2545673025434913e-199, -1.4316109244036342e-62,
      -1.1471703763678777e-15, -7.4744582105950359e-08, -0.021392115008602703,
      -0.026381238386562497, -0.026383183529569186, -0.026383183529569186
    ),
    dim = c(m, 6, 6)
  )

  H <- logistic6_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 <- ltd$stats_1[, 1]
  theta <- ltd$theta_6

  m <- length(x)

  true_gradient <- matrix(
    c(
      # alpha
      rep(1, m),
      # delta
      9.6437492398195889e-23, 1.5295116025091289e-07, 0.024893527014351497,
      0.18378213803748404, 0.74233581033456154, 0.75982946099000246,
      0.75983568565159255, 0.75983568565159255,
      # eta
      2.4276486539439174e-21, 1.1736381182546699e-06, 0.041653962050552102,
      0.12331460190108316, -0.021582949777827438, -0.000033032686382156039,
      -1.8583072261039617e-42, -1.3923270509654377e-85,
      # phi
      4.8218746199097945e-23, 7.6475580125456447e-08, 0.012446749168025124,
      0.091576578355026419, 0.033030421569699108, 0.000012449068216533443,
      7.5377239464421580e-44, 2.8040905767036438e-87,
      # nu
      -1.2161526556915815e-21, -5.9051181979806446e-07, -0.021428022923202986,
      -0.066384858676451105, -0.051166597326927950, -0.052172373971476468,
      -0.052172801351587403, -0.052172801351587403,
      # xi
      0, 2.0926851264125609e-35, 2.3898584374049245e-09,
      0.000052415110619266708, 0.056356247266263610, 0.063317046904464116,
      0.063319640470966046, 0.063319640470966046
    ),
    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)
      -2.4276486539439174e-21, -1.1736381182546699e-06, -0.041653962050552102,
      -0.12331460190108316, 0.021582949777827438, 0.000033032686382156039,
      1.8583072261039617e-42, 1.3923270509654377e-85,
      # (delta, phi)
      -4.8218746199097945e-23, -7.6475580125456447e-08, -0.012446749168025124,
      -0.091576578355026419, -0.033030421569699108, -0.000012449068216533443,
      -7.5377239464421580e-44, -2.8040905767036438e-87,
      # (delta, nu)
      1.2161526556915815e-21, 5.9051181979806446e-07, 0.021428022923202986,
      0.066384858676451105, 0.051166597326927950, 0.052172373971476468,
      0.052172801351587403, 0.052172801351587403,
      # (delta, xi)
      -2.0853048558916222e-111, -2.0926851264125609e-35,
      -2.3898584374049245e-09, -0.000052415110619266708, -0.056356247266263610,
      -0.063317046904464116, -0.063319640470966046, -0.063319640470966046,
      # (eta, alpha)
      rep(0, m),
      # (eta, delta)
      -2.4276486539439174e-21, -1.1736381182546699e-06, -0.041653962050552102,
      -0.12331460190108316, 0.021582949777827438, 0.000033032686382156039,
      1.8583072261039617e-42, 1.3923270509654377e-85,
      # (eta, eta)
      -6.1111895803565778e-20, -9.0056618750765000e-06, -0.069698623184655625,
      -0.081605338544629843, 0.025068184884160287, 0.00017529242457613078,
      9.1627280890809950e-41, 1.3826761752674973e-83,
      # (eta, phi)
      -1.1897149538724097e-21, -5.4858126906460671e-07, -0.014603486473844729,
      -0.014813923084545845, -0.021848994443557445, -0.000059838132613523584,
      -3.6789258324757125e-42, -2.7706336490473572e-85,
      # (eta, nu)
      3.0007647424909336e-20, 4.2377570196669037e-06, 0.025441707980853873,
      0.013819938082184523, -0.0010074674090324482, -2.2678615960535156e-06,
      -1.2759744716990998e-43, -9.5601725502228708e-87,
      # (eta, xi)
      -2.6246988596328932e-109, -8.0288865734436783e-34,
      -1.9994569793277198e-08, -0.00017584811474708742, 0.0081926133528867864,
      0.000013763168318811387, 7.7429467754331736e-43, 5.8013627123559905e-86,
      # (phi, alpha)
      rep(0, m),
      # (phi, delta)
      -4.8218746199097945e-23, -7.6475580125456447e-08, -0.012446749168025124,
      -0.091576578355026419, -0.033030421569699108, -0.000012449068216533443,
      -7.5377239464421580e-44, -2.8040905767036438e-87,
      # (phi, eta)
      -1.1897149538724097e-21, -5.4858126906460671e-07, -0.014603486473844729,
      -0.014813923084545845, -0.021848994443557445, -0.000059838132613523584,
      -3.6789258324757125e-42, -2.7706336490473572e-85,
      # (phi, phi)
      -2.4109373099548972e-23, -3.8237790062728224e-08, -0.0062233387361772991,
      -0.045004753323180023, 0.058712357899407301, 0.000024897116603569502,
      1.5075447892884316e-43, 5.6081811534072877e-87,
      # (phi, nu)
      5.9602164129601628e-22, 2.7613701486766812e-07, 0.0076023154114251624,
      0.010263039600614143, 0.0015418222986489271, 8.5469172528994537e-07,
      5.1756475975932719e-45, 1.9253802288156686e-88,
      # (phi, xi)
      -5.2132621397290555e-111, -5.2317128160314022e-35,
      -5.9746392104900531e-09, -0.00013058930905553107, -0.012537928114042809,
      -5.1869417853060373e-06, -3.1407183110175658e-44, -1.1683710736265183e-87,
      # (nu, alpha)
      rep(0, m),
      # (nu, delta)
      1.2161526556915815e-21, 5.9051181979806446e-07, 0.021428022923202986,
      0.066384858676451105, 0.051166597326927950, 0.052172373971476468,
      0.052172801351587403, 0.052172801351587403,
      # (nu, eta)
      3.0007647424909336e-20, 4.2377570196669037e-06, 0.025441707980853873,
      0.013819938082184523, -0.0010074674090324482, -2.2678615960535156e-06,
      -1.2759744716990998e-43, -9.5601725502228708e-87,
      # (nu, phi)
      5.9602164129601628e-22, 2.7613701486766812e-07, 0.0076023154114251624,
      0.010263039600614143, 0.0015418222986489271, 8.5469172528994537e-07,
      5.1756475975932719e-45, 1.9253802288156686e-88,
      # (nu, nu)
      -1.4730072359242754e-20, -1.9869741882303387e-06, -0.0081199112657906938,
      0.0063612510169811997, 0.021964708947817155, 0.022503861272861439,
      0.022504045632100321, 0.022504045632100321,
      # (nu, xi)
      2.5254682512573550e-110, 7.0330685392331877e-35, 8.6223039892301619e-10,
      -7.2295861790158238e-06, -0.011458416654500857, -0.011482225054956227,
      -0.011482176671775895, -0.011482176671775895,
      # (xi, alpha)
      rep(0, m),
      # (xi, delta)
      -2.0853048558916222e-111, -2.0926851264125609e-35,
      -2.3898584374049245e-09, -0.000052415110619266708, -0.056356247266263610,
      -0.063317046904464116, -0.063319640470966046, -0.063319640470966046,
      # (xi, eta)
      -2.6246988596328932e-109, -8.0288865734436783e-34,
      -1.9994569793277198e-08, -0.00017584811474708742, 0.0081926133528867864,
      0.000013763168318811387, 7.7429467754331736e-43, 5.8013627123559905e-86,
      # (xi, phi)
      -5.2132621397290555e-111, -5.2317128160314022e-35,
      -5.9746392104900531e-09, -0.00013058930905553107, -0.012537928114042809,
      -5.1869417853060373e-06, -3.1407183110175658e-44, -1.1683710736265183e-87,
      # (xi, nu)
      2.5254682512573550e-110, 7.0330685392331877e-35, 8.6223039892301619e-10,
      -7.2295861790158238e-06, -0.011458416654500857, -0.011482225054956227,
      -0.011482176671775895, -0.011482176671775895,
      # (xi, xi)
      -2.2545673025434913e-199, -1.4316109244036342e-62,
      -1.1471703763678777e-15, -7.4744582105950359e-08, -0.021392115008602703,
      -0.026381238386562497, -0.026383183529569186, -0.026383183529569186
    ),
    dim = c(m, 6, 6)
  )

  gh <- logistic6_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 <- ltd$stats_1[, 1]
  theta <- ltd$theta_6

  m <- length(x)

  true_gradient <- matrix(
    c(
      # alpha
      rep(1, m),
      # delta
      9.6437492398195889e-23, 1.5295116025091289e-07, 0.024893527014351497,
      0.18378213803748404, 0.74233581033456154, 0.75982946099000246,
      0.75983568565159255, 0.75983568565159255,
      # log_eta
      4.8552973078878348e-21, 2.3472762365093397e-06, 0.083307924101104204,
      0.24662920380216631, -0.043165899555654876, -0.000066065372764312077,
      -3.7166144522079233e-42, -2.7846541019308754e-85,
      # phi
      4.8218746199097945e-23, 7.6475580125456447e-08, 0.012446749168025124,
      0.091576578355026419, 0.033030421569699108, 0.000012449068216533443,
      7.5377239464421580e-44, 2.8040905767036438e-87,
      # log_nu
      -4.8646106227663261e-21, -2.3620472791922578e-06, -0.085712091692811944,
      -0.26553943470580442, -0.20466638930771180, -0.20868949588590587,
      -0.20869120540634961, -0.20869120540634961,
      # log_xi
      0, 6.2780553792376827e-35, 7.1695753122147735e-09, 0.00015724533185780012,
      0.16906874179879083, 0.18995114071339235, 0.18995892141289814,
      0.18995892141289814
    ),
    nrow = m,
    ncol = 6
  )

  G <- logistic6_gradient_2(x, theta)

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

test_that("Hessian (2)", {
  x <- ltd$stats_1[, 1]
  theta <- ltd$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, 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)
      -4.8552973078878348e-21, -2.3472762365093397e-06, -0.083307924101104204,
      -0.24662920380216631, 0.043165899555654876, 0.000066065372764312077,
      3.7166144522079233e-42, 2.7846541019308754e-85,
      # (delta, phi)
      -4.8218746199097945e-23, -7.6475580125456447e-08, -0.012446749168025124,
      -0.091576578355026419, -0.033030421569699108, -0.000012449068216533443,
      -7.5377239464421580e-44, -2.8040905767036438e-87,
      # (delta, log_nu)
      4.8646106227663261e-21, 2.3620472791922578e-06, 0.085712091692811944,
      0.26553943470580442, 0.20466638930771180, 0.20868949588590587,
      0.20869120540634961, 0.20869120540634961,
      # (delta, log_xi)
      0, -6.2780553792376827e-35, -7.1695753122147735e-09,
      -0.00015724533185780012, -0.16906874179879083, -0.18995114071339235,
      -0.18995892141289814, -0.18995892141289814,
      # (log_eta, alpha)
      rep(0, m),
      # (log_eta, delta)
      -4.8552973078878348e-21, -2.3472762365093397e-06, -0.083307924101104204,
      -0.24662920380216631, 0.043165899555654876, 0.000066065372764312077,
      3.7166144522079233e-42, 2.7846541019308754e-85,
      # (log_eta, log_eta)
      -2.3959228590637528e-19, -0.000033675371263796660, -0.19548656863751830,
      -0.079792150376353060, 0.057106839980986273, 0.00063510432554021104,
      3.6279250911103188e-40, 5.5028581600506806e-83,
      # (log_eta, phi)
      -2.3794299077448194e-21, -1.0971625381292134e-06, -0.029206972947689459,
      -0.029627846169091689, -0.043697988887114891, -0.00011967626522704717,
      -7.3578516649514251e-42, -5.5412672980947144e-85,
      # (log_eta, log_nu)
      2.4006117939927469e-19, 0.000033902056157335230, 0.20353366384683098,
      0.11055950465747618, -0.0080597392722595856, -0.000018142892768428125,
      -1.0207795773592799e-42, -7.6481380401782966e-86,
      # (log_eta, log_xi)
      0, -4.8173319440662070e-33, -1.1996741875966319e-07,
      -0.0010550886884825245, 0.049155680117320718, 0.000082579009912868325,
      4.6457680652599042e-42, 3.4808176274135943e-85,
      # (phi, alpha)
      rep(0, m),
      # (phi, delta)
      -4.8218746199097945e-23, -7.6475580125456447e-08, -0.012446749168025124,
      -0.091576578355026419, -0.033030421569699108, -0.000012449068216533443,
      -7.5377239464421580e-44, -2.8040905767036438e-87,
      # (phi, log_eta)
      -2.3794299077448194e-21, -1.0971625381292134e-06, -0.029206972947689459,
      -0.029627846169091689, -0.043697988887114891, -0.00011967626522704717,
      -7.3578516649514251e-42, -5.5412672980947144e-85,
      # (phi, phi)
      -2.4109373099548972e-23, -3.8237790062728224e-08, -0.0062233387361772991,
      -0.045004753323180023, 0.058712357899407301, 0.000024897116603569502,
      1.5075447892884316e-43, 5.6081811534072877e-87,
      # (phi, log_nu)
      2.3840865651840651e-21, 1.1045480594706725e-06, 0.030409261645700650,
      0.041052158402456571, 0.0061672891945957084, 3.4187669011597815e-06,
      2.0702590390373088e-44, 7.7015209152626744e-88,
      # (phi, log_xi)
      0, -1.5695138448094207e-34, -1.7923917631470159e-08,
      -0.00039176792716659321, -0.037613784342128428, -0.000015560825355918112,
      -9.4221549330526975e-44, -3.5051132208795548e-87,
      # (log_nu, alpha)
      rep(0, m),
      # (log_nu, delta)
      4.8646106227663261e-21, 2.3620472791922578e-06, 0.085712091692811944,
      0.26553943470580442, 0.20466638930771180, 0.20868949588590587,
      0.20869120540634961, 0.20869120540634961,
      # (log_nu, log_eta)
      2.4006117939927469e-19, 0.000033902056157335230, 0.20353366384683098,
      0.11055950465747618, -0.0080597392722595856, -0.000018142892768428125,
      -1.0207795773592799e-42, -7.6481380401782966e-86,
      # (log_nu, phi)
      2.3840865651840651e-21, 1.1045480594706725e-06, 0.030409261645700650,
      0.041052158402456571, 0.0061672891945957084, 3.4187669011597815e-06,
      2.0702590390373088e-44, 7.7015209152626744e-88,
      # (log_nu, log_nu)
      -2.4054576837065039e-19, -0.000034153634290877677, -0.21563067194546304,
      -0.16375941843410523, 0.14676895385736267, 0.15137228447987715,
      0.15137352470725553, 0.15137352470725553,
      # (log_nu, log_xi)
      0, 8.4396822470798252e-34, 1.0346764787076194e-08,
      -0.000086755034148189886, -0.13750099985401029, -0.13778670065947473,
      -0.13778612006131073, -0.13778612006131073,
      # (log_xi, alpha)
      rep(0, m),
      # (log_xi, delta)
      0, -6.2780553792376827e-35, -7.1695753122147735e-09,
      -0.00015724533185780012, -0.16906874179879083, -0.18995114071339235,
      -0.18995892141289814, -0.18995892141289814,
      # (log_xi, log_eta)
      0, -4.8173319440662070e-33, -1.1996741875966319e-07,
      -0.0010550886884825245, 0.049155680117320718, 0.000082579009912868325,
      4.6457680652599042e-42, 3.4808176274135943e-85,
      # (log_xi, phi)
      0, -1.5695138448094207e-34, -1.7923917631470159e-08,
      -0.00039176792716659321, -0.037613784342128428, -0.000015560825355918112,
      -9.4221549330526975e-44, -3.5051132208795548e-87,
      # (log_xi, log_nu)
      0, 8.4396822470798252e-34, 1.0346764787076194e-08,
      -0.000086755034148189886, -0.13750099985401029, -0.13778670065947473,
      -0.13778612006131073, -0.13778612006131073,
      # (log_xi, log_xi)
      0, 6.2780553792376827e-35, 7.1695649876813862e-09,
      0.00015657263061884657, -0.023460293278633494, -0.047480004765670128,
      -0.047489730353224534, -0.047489730353224534
    ),
    dim = c(m, 6, 6)
  )

  H <- logistic6_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 <- ltd$stats_1[, 1]
  theta <- ltd$theta_6

  m <- length(x)

  true_gradient <- matrix(
    c(
      # alpha
      rep(1, m),
      # delta
      9.6437492398195889e-23, 1.5295116025091289e-07, 0.024893527014351497,
      0.18378213803748404, 0.74233581033456154, 0.75982946099000246,
      0.75983568565159255, 0.75983568565159255,
      # log_eta
      4.8552973078878348e-21, 2.3472762365093397e-06, 0.083307924101104204,
      0.24662920380216631, -0.043165899555654876, -0.000066065372764312077,
      -3.7166144522079233e-42, -2.7846541019308754e-85,
      # phi
      4.8218746199097945e-23, 7.6475580125456447e-08, 0.012446749168025124,
      0.091576578355026419, 0.033030421569699108, 0.000012449068216533443,
      7.5377239464421580e-44, 2.8040905767036438e-87,
      # log_nu
      -4.8646106227663261e-21, -2.3620472791922578e-06, -0.085712091692811944,
      -0.26553943470580442, -0.20466638930771180, -0.20868949588590587,
      -0.20869120540634961, -0.20869120540634961,
      # log_xi
      0, 6.2780553792376827e-35, 7.1695753122147735e-09, 0.00015724533185780012,
      0.16906874179879083, 0.18995114071339235, 0.18995892141289814,
      0.18995892141289814
    ),
    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, 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)
      -4.8552973078878348e-21, -2.3472762365093397e-06, -0.083307924101104204,
      -0.24662920380216631, 0.043165899555654876, 0.000066065372764312077,
      3.7166144522079233e-42, 2.7846541019308754e-85,
      # (delta, phi)
      -4.8218746199097945e-23, -7.6475580125456447e-08, -0.012446749168025124,
      -0.091576578355026419, -0.033030421569699108, -0.000012449068216533443,
      -7.5377239464421580e-44, -2.8040905767036438e-87,
      # (delta, log_nu)
      4.8646106227663261e-21, 2.3620472791922578e-06, 0.085712091692811944,
      0.26553943470580442, 0.20466638930771180, 0.20868949588590587,
      0.20869120540634961, 0.20869120540634961,
      # (delta, log_xi)
      0, -6.2780553792376827e-35, -7.1695753122147735e-09,
      -0.00015724533185780012, -0.16906874179879083, -0.18995114071339235,
      -0.18995892141289814, -0.18995892141289814,
      # (log_eta, alpha)
      rep(0, m),
      # (log_eta, delta)
      -4.8552973078878348e-21, -2.3472762365093397e-06, -0.083307924101104204,
      -0.24662920380216631, 0.043165899555654876, 0.000066065372764312077,
      3.7166144522079233e-42, 2.7846541019308754e-85,
      # (log_eta, log_eta)
      -2.3959228590637528e-19, -0.000033675371263796660, -0.19548656863751830,
      -0.079792150376353060, 0.057106839980986273, 0.00063510432554021104,
      3.6279250911103188e-40, 5.5028581600506806e-83,
      # (log_eta, phi)
      -2.3794299077448194e-21, -1.0971625381292134e-06, -0.029206972947689459,
      -0.029627846169091689, -0.043697988887114891, -0.00011967626522704717,
      -7.3578516649514251e-42, -5.5412672980947144e-85,
      # (log_eta, log_nu)
      2.4006117939927469e-19, 0.000033902056157335230, 0.20353366384683098,
      0.11055950465747618, -0.0080597392722595856, -0.000018142892768428125,
      -1.0207795773592799e-42, -7.6481380401782966e-86,
      # (log_eta, log_xi)
      0, -4.8173319440662070e-33, -1.1996741875966319e-07,
      -0.0010550886884825245, 0.049155680117320718, 0.000082579009912868325,
      4.6457680652599042e-42, 3.4808176274135943e-85,
      # (phi, alpha)
      rep(0, m),
      # (phi, delta)
      -4.8218746199097945e-23, -7.6475580125456447e-08, -0.012446749168025124,
      -0.091576578355026419, -0.033030421569699108, -0.000012449068216533443,
      -7.5377239464421580e-44, -2.8040905767036438e-87,
      # (phi, log_eta)
      -2.3794299077448194e-21, -1.0971625381292134e-06, -0.029206972947689459,
      -0.029627846169091689, -0.043697988887114891, -0.00011967626522704717,
      -7.3578516649514251e-42, -5.5412672980947144e-85,
      # (phi, phi)
      -2.4109373099548972e-23, -3.8237790062728224e-08, -0.0062233387361772991,
      -0.045004753323180023, 0.058712357899407301, 0.000024897116603569502,
      1.5075447892884316e-43, 5.6081811534072877e-87,
      # (phi, log_nu)
      2.3840865651840651e-21, 1.1045480594706725e-06, 0.030409261645700650,
      0.041052158402456571, 0.0061672891945957084, 3.4187669011597815e-06,
      2.0702590390373088e-44, 7.7015209152626744e-88,
      # (phi, log_xi)
      0, -1.5695138448094207e-34, -1.7923917631470159e-08,
      -0.00039176792716659321, -0.037613784342128428, -0.000015560825355918112,
      -9.4221549330526975e-44, -3.5051132208795548e-87,
      # (log_nu, alpha)
      rep(0, m),
      # (log_nu, delta)
      4.8646106227663261e-21, 2.3620472791922578e-06, 0.085712091692811944,
      0.26553943470580442, 0.20466638930771180, 0.20868949588590587,
      0.20869120540634961, 0.20869120540634961,
      # (log_nu, log_eta)
      2.4006117939927469e-19, 0.000033902056157335230, 0.20353366384683098,
      0.11055950465747618, -0.0080597392722595856, -0.000018142892768428125,
      -1.0207795773592799e-42, -7.6481380401782966e-86,
      # (log_nu, phi)
      2.3840865651840651e-21, 1.1045480594706725e-06, 0.030409261645700650,
      0.041052158402456571, 0.0061672891945957084, 3.4187669011597815e-06,
      2.0702590390373088e-44, 7.7015209152626744e-88,
      # (log_nu, log_nu)
      -2.4054576837065039e-19, -0.000034153634290877677, -0.21563067194546304,
      -0.16375941843410523, 0.14676895385736267, 0.15137228447987715,
      0.15137352470725553, 0.15137352470725553,
      # (log_nu, log_xi)
      0, 8.4396822470798252e-34, 1.0346764787076194e-08,
      -0.000086755034148189886, -0.13750099985401029, -0.13778670065947473,
      -0.13778612006131073, -0.13778612006131073,
      # (log_xi, alpha)
      rep(0, m),
      # (log_xi, delta)
      0, -6.2780553792376827e-35, -7.1695753122147735e-09,
      -0.00015724533185780012, -0.16906874179879083, -0.18995114071339235,
      -0.18995892141289814, -0.18995892141289814,
      # (log_xi, log_eta)
      0, -4.8173319440662070e-33, -1.1996741875966319e-07,
      -0.0010550886884825245, 0.049155680117320718, 0.000082579009912868325,
      4.6457680652599042e-42, 3.4808176274135943e-85,
      # (log_xi, phi)
      0, -1.5695138448094207e-34, -1.7923917631470159e-08,
      -0.00039176792716659321, -0.037613784342128428, -0.000015560825355918112,
      -9.4221549330526975e-44, -3.5051132208795548e-87,
      # (log_xi, log_nu)
      0, 8.4396822470798252e-34, 1.0346764787076194e-08,
      -0.000086755034148189886, -0.13750099985401029, -0.13778670065947473,
      -0.13778612006131073, -0.13778612006131073,
      # (log_xi, log_xi)
      0, 6.2780553792376827e-35, 7.1695649876813862e-09,
      0.00015657263061884657, -0.023460293278633494, -0.047480004765670128,
      -0.047489730353224534, -0.047489730353224534
    ),
    dim = c(m, 6, 6)
  )

  gh <- logistic6_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 = ltd$stats_1), class = "logistic6")

  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 <- ltd$theta_6
  theta[c(3, 5, 6)] <- log(theta[c(3, 5, 6)])

  true_value <- 0.36255789619259319

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

  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 <- ltd$theta_6
  theta[c(3, 5, 6)] <- log(theta[c(3, 5, 6)])

  true_gradient <- c(
    -0.59376001498058462, -0.73330143767252627, 0.17369846729420963,
    0.017669970115745026, 0.080623092690263486, -0.19216034447600408
  )

  true_hessian <- matrix(
    c(
      # alpha
      19, 6.0512600139805846, 1.0235759411998543, 0.49030375496021016,
      -2.6823512883787693, 1.2676631260150411,
      # delta
      6.0512600139805846, 4.0989218556934594, -0.084427619158893432,
      0.12383872776827310, -1.3701711378917644, 1.1461364594573599,
      # log_eta
      1.0235759411998543, -0.084427619158893432, 0.12128671379906373,
      0.10305675773587517, -0.13332557577248706, -0.063800298599788432,
      # phi
      0.49030375496021016, 0.12383872776827310, 0.10305675773587517,
      -0.034948837041012630, -0.097276873662293284, 0.048434100462574590,
      # log_nu
      -2.6823512883787693, -1.3701711378917644, -0.13332557577248706,
      -0.097276873662293284, 0.29576678334915398, -0.11057898822366290,
      # log_xi
      1.2676631260150411, 1.1461364594573599, -0.063800298599788432,
      0.048434100462574590, -0.11057898822366290, 0.26231706463305267
    ),
    nrow = 6,
    ncol = 6
  )

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

  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 <- ltd$D$x
  y <- ltd$D$y
  w <- rep(1, length(y))

  max_iter <- 10000

  theta <- c(
    0, 1, -0.13948206931816815, 3.9201472292959156, 2.1396800314205292,
    0.54368155924756683
  )

  true_value <- c(
    0.84948649909664558, -0.80697277910347021, -0.13948206931816815,
    3.9201472292959156, 2.1396800314205292, 0.54368155924756683
  )

  object <- logistic6_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 <- ltd$D$x
  y <- ltd$D$y

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

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

  max_iter <- 10000

  # logistic6 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.050168602337768019

  fitted_values <- rep(
    c(
      0.84947145597244308, 0.83001285437491464, 0.62226246927410375,
      0.50732524066346644, 0.33637480874558723, 0.14349972605987800,
      0.092532338486404620, 0.092532338486404619
    ),
    k
  )

  residuals <- c(
    0.0032285440275569165, -0.092371455972443084, 0.087528544027556916,
    -0.0025128543749146426, -0.052212854374914643, 0.054587145625085357,
    -0.066162469274103751, 0.082537530725896249, 0.037974759336533557,
    -0.021325240663466443, 0.027774759336533557, -0.067425240663466443,
    0.018925191254412770, -0.023174808745587230, 0.013525191254412770,
    -0.0013997260598779998, -0.075732338486404620, 0.046067661513595380,
    0.030167661513595381
  )

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

  result <- fit(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_new(x, y, w, c(0, 1, 1, 1, 1, 1), max_iter, NULL, NULL)

  result <- fit(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- ltd$D$x
  y <- ltd$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.050168602337768019

  fitted_values <- rep(
    c(
      0.84947145597244308, 0.83001285437491464, 0.62226246927410375,
      0.50732524066346644, 0.33637480874558723, 0.14349972605987800,
      0.092532338486404620, 0.092532338486404619
    ),
    k
  )

  residuals <- c(
    0.0032285440275569165, -0.092371455972443084, 0.087528544027556916,
    -0.0025128543749146426, -0.052212854374914643, 0.054587145625085357,
    -0.066162469274103751, 0.082537530725896249, 0.037974759336533557,
    -0.021325240663466443, 0.027774759336533557, -0.067425240663466443,
    0.018925191254412770, -0.023174808745587230, 0.013525191254412770,
    -0.0013997260598779998, -0.075732338486404620, 0.046067661513595380,
    0.030167661513595381
  )

  object <- logistic6_new(
    x, y, w, NULL, max_iter,
    c(0.5, -1, 0.05, -5, 2, 0.5), c(1, -0.5, 5, 5, 10, 3)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_new(
    x, y, w, c(0.7, -0.6, 3, -1, 4, 1), max_iter,
    c(0.5, -1, 0.05, -5, 2, 0.5), c(1, -0.5, 5, 5, 10, 3)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_new(
    x, y, w, c(-2, 2, 7, -8, 1, 5), max_iter,
    c(0.5, -1, 0.05, -5, 2, 0.5), c(1, -0.5, 5, 5, 10, 3)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- ltd$D$x
  y <- ltd$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.062057350409983716

  fitted_values <- rep(
    c(
      0.79999935671542949, 0.79295099790131815, 0.62904933718751053,
      0.51070978631265628, 0.32891952067935980, 0.15665231157393833,
      0.090929983963630982, 0.090929983960543027
    ),
    k
  )

  residuals <- c(
    0.052700643284570514, -0.042899356715429486, 0.13700064328457051,
    0.034549002098681851, -0.015150997901318149, 0.091649002098681851,
    -0.072949337187510533, 0.075750662812489467, 0.034590213687343716,
    -0.024709786312656284, 0.024390213687343716, -0.070809786312656284,
    0.026380479320640199, -0.015719520679359801, 0.020980479320640199,
    -0.014552311573938331, -0.074129983963630982, 0.047670016036369018,
    0.031770016039456973
  )

  object <- logistic6_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, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_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, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_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, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- ltd$D$x
  y <- ltd$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.062057350409983716

  fitted_values <- rep(
    c(
      0.79999935671542949, 0.79295099790131815, 0.62904933718751053,
      0.51070978631265628, 0.32891952067935980, 0.15665231157393833,
      0.090929983963630982, 0.090929983960543027
    ),
    k
  )

  residuals <- c(
    0.052700643284570514, -0.042899356715429486, 0.13700064328457051,
    0.034549002098681851, -0.015150997901318149, 0.091649002098681851,
    -0.072949337187510533, 0.075750662812489467, 0.034590213687343716,
    -0.024709786312656284, 0.024390213687343716, -0.070809786312656284,
    0.026380479320640199, -0.015719520679359801, 0.020980479320640199,
    -0.014552311573938331, -0.074129983963630982, 0.047670016036369018,
    0.031770016039456973
  )

  object <- logistic6_new(
    x, y, w, NULL, max_iter,
    c(0.8, -0.9, 0.05, -5, 2, 1), c(0.8, -0.9, 3, 5, 6, 3)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_new(
    x, y, w, c(0.8, -0.9, 2, 0, 2.5, 1.2), max_iter,
    c(0.8, -0.9, 0.05, -5, 2, 1), c(0.8, -0.9, 3, 5, 6, 3)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_new(
    x, y, w, c(0, 1, 5, -8, 10, 5), max_iter,
    c(0.8, -0.9, 0.05, -5, 2, 1), c(0.8, -0.9, 3, 5, 6, 3)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- ltd$D$x
  y <- ltd$D$y
  w <- ltd$D$w

  n <- length(y)

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

  max_iter <- 10000

  # logistic6 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.026367789609414469

  fitted_values <- rep(
    c(
      0.85062856577440989, 0.85038705397619911, 0.68988971453996974,
      0.50866177107419315, 0.30821950267475678, 0.17833789292363378,
      0.087927106108376473, 0.087926555053940805
    ),
    k
  )

  residuals <- c(
    0.0020714342255901108, -0.093528565774409889, 0.086371434225590111,
    -0.022887053976199108, -0.072587053976199108, 0.034212946023800892,
    -0.13378971453996974, 0.014910285460030258, 0.036638228925806847,
    -0.022661771074193153, 0.026438228925806847, -0.068761771074193153,
    0.047080497325243221, 0.0049804973252432212, 0.041680497325243221,
    -0.036237892923633776, -0.071127106108376473, 0.050672893891623527,
    0.034773444946059195
  )

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

  result <- fit(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_new(
    x, y, w, c(1, -1, 1, 1, 1, 1), max_iter, NULL, NULL
  )

  result <- fit(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- ltd$D$x
  y <- ltd$D$y
  w <- ltd$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.026367789609414469

  fitted_values <- rep(
    c(
      0.85062856577440989, 0.85038705397619911, 0.68988971453996974,
      0.50866177107419315, 0.30821950267475678, 0.17833789292363378,
      0.087927106108376473, 0.087926555053940805
    ),
    k
  )

  residuals <- c(
    0.0020714342255901108, -0.093528565774409889, 0.086371434225590111,
    -0.022887053976199108, -0.072587053976199108, 0.034212946023800892,
    -0.13378971453996974, 0.014910285460030258, 0.036638228925806847,
    -0.022661771074193153, 0.026438228925806847, -0.068761771074193153,
    0.047080497325243221, 0.0049804973252432212, 0.041680497325243221,
    -0.036237892923633776, -0.071127106108376473, 0.050672893891623527,
    0.034773444946059195
  )

  object <- logistic6_new(
    x, y, w, NULL, max_iter,
    c(0.5, -1, 0.05, -5, 0.5, 0), c(1, -0.5, 3, 5, 3, 2)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_new(
    x, y, w, c(0.7, -0.6, 1, 2, 2.5, 1.8), max_iter,
    c(0.5, -1, 0.05, -5, 0.5, 0), c(1, -0.5, 3, 5, 3, 2)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_new(
    x, y, w, c(-2, -5, 4, -8, 4, 3), 10000,
    c(0.5, -1, 0.05, -5, 0.5, 0), c(1, -0.5, 3, 5, 3, 2)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- ltd$D$x
  y <- ltd$D$y
  w <- ltd$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.032325384005503981

  fitted_values <- rep(
    c(
      0.89996907976852931, 0.87386569158187010, 0.63652439966625265,
      0.51281254109481959, 0.33374773386909931, 0.14206152371471766,
      0.093890198677872793, 0.093890198677872789
    ),
    k
  )

  residuals <- c(
    -0.047269079768529308, -0.14286907976852931, 0.037030920231470692,
    -0.046365691581870097, -0.096065691581870097, 0.010734308418129903,
    -0.080424399666252655, 0.068275600333747345, 0.032487458905180414,
    -0.026812541094819586, 0.022287458905180414, -0.072912541094819586,
    0.021552266130900692, -0.020547733869099308, 0.016152266130900692,
    0.000038476285282337392, -0.077090198677872793, 0.044709801322127207,
    0.028809801322127211
  )

  object <- logistic6_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, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_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, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_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, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- ltd$D$x
  y <- ltd$D$y
  w <- ltd$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.032325384005503981

  fitted_values <- rep(
    c(
      0.89996907976852931, 0.87386569158187010, 0.63652439966625265,
      0.51281254109481959, 0.33374773386909931, 0.14206152371471766,
      0.093890198677872793, 0.093890198677872789
    ),
    k
  )

  residuals <- c(
    -0.047269079768529308, -0.14286907976852931, 0.037030920231470692,
    -0.046365691581870097, -0.096065691581870097, 0.010734308418129903,
    -0.080424399666252655, 0.068275600333747345, 0.032487458905180414,
    -0.026812541094819586, 0.022287458905180414, -0.072912541094819586,
    0.021552266130900692, -0.020547733869099308, 0.016152266130900692,
    0.000038476285282337392, -0.077090198677872793, 0.044709801322127207,
    0.028809801322127211
  )

  object <- logistic6_new(
    x, y, w, NULL, max_iter,
    c(0.9, -0.9, 0.05, -5, 5, 0.5), c(0.9, -0.9, 3, 5, 12, 4)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_new(
    x, y, w, c(0.9, -0.9, 2, 0, 11, 1), max_iter,
    c(0.9, -0.9, 0.05, -5, 5, 0.5), c(0.9, -0.9, 3, 5, 12, 4)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- logistic6_new(
    x, y, w, c(0, 1, 5, -8, 3, 6), max_iter,
    c(0.9, -0.9, 0.05, -5, 5, 0.5), c(0.9, -0.9, 3, 5, 12, 4)
  )

  result <- fit_constrained(object)

  expect_true(inherits(result, "logistic6_fit"))
  expect_true(inherits(result, "logistic"))
  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 <- ltd$D$x
  y <- ltd$D$y
  w <- ltd$D$w

  max_iter <- 10000

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

  sigma <- ltd$sigma

  true_value <- matrix(c(
      # alpha
      6206.96, 2093.9776752822887, 153.30946172983291, 153.11183995442809,
      -221.36872037331707, 148.97949961445271, 11595.531011291549,
      # delta
      2093.9776752822887, 1437.5265454670404, -14.586523790688970,
      40.797335330545767, -120.18243242866359, 134.67400631342478,
      10406.248038205941,
      # eta
      153.30946172983291, -14.586523790688970, -2.2648840073304669,
      16.727779439929406, -5.9742410009604452, -3.5836831170680006,
      -1060.9581981747777,
      # phi
      153.11183995442809, 40.797335330545767, 16.727779439929406,
      -12.239261546191805, -7.8866010818907720, 5.4396548685774482,
      -201.51541062927098,
      # nu
      -221.36872037331707, -120.18243242866359, -5.9742410009604452,
      -7.8866010818907720, 4.3346830493633020, -3.2808483355079093,
      -358.18422398482191,
      # xi
      148.97949961445271, 134.67400631342478, -3.5836831170680006,
      5.4396548685774482, -3.2808483355079093, 17.822668367095048,
      900.77323828870695,
      # sigma
      11595.531011291549, 10406.248038205941, -1060.9581981747777,
      -201.51541062927098, -358.18422398482191, 900.77323828870695,
      139735.83919814804
    ),
    nrow = 7,
    ncol = 7
  )

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

  object <- logistic6_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 <- ltd$D$x
  y <- ltd$D$y

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

  expect_error(
    drda(
      y ~ x, mean_function = "logistic6",
      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 = "logistic6",
      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 = "logistic6",
      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 = "logistic6",
      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 = "logistic6",
      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 <- ltd$D$x
  y <- ltd$D$y

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

  expect_error(
    drda(
      y ~ x, mean_function = "logistic6",
      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 = "logistic6",
      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 = "logistic6",
      lower_bound = rep(-Inf, 7),
      upper_bound = rep(Inf, 7)
    ),
    "'lower_bound' must be of length 6"
  )
})

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

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

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

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

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

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

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

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

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

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

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

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

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

  expect_equal(nauc(result), 0.42753516878440054)
  expect_equal(nauc(result, xlim = c(-2, 2)), 0.40530818958837455)
  expect_equal(nauc(result, ylim = c(0.3, 0.7)), 0.40513264312426645)
  expect_equal(nauc(result, xlim = c(-15, -10), ylim = c(0.3, 0.7)), 1.0)
  expect_equal(
    nauc(result, xlim = c(1, 5), ylim = c(0.3, 0.7)), 0.019678613960958503
  )
  expect_equal(nauc(result, xlim = c(10, 15), ylim = c(0.3, 0.7)), 0.0)
})

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

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

  expect_equal(naac(result), 1 - 0.42753516878440054)
  expect_equal(naac(result, xlim = c(-2, 2)), 1 - 0.40530818958837455)
  expect_equal(naac(result, ylim = c(0.3, 0.7)), 1 - 0.40513264312426645)
  expect_equal(naac(result, xlim = c(-15, -10), ylim = c(0.3, 0.7)), 0.0)
  expect_equal(
    naac(result, xlim = c(1, 5), ylim = c(0.3, 0.7)), 1 - 0.019678613960958503
  )
  expect_equal(naac(result, xlim = c(10, 15), ylim = c(0.3, 0.7)), 1.0)
})

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

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

  expect_equal(nauc(result), 0.65302073661957048)
  expect_equal(nauc(result, xlim = c(-2, 2)), 0.65733232476226463)
  expect_equal(nauc(result, ylim = c(0.3, 0.7)), 0.74700167012690062)
  expect_equal(nauc(result, xlim = c(-25, -18), ylim = c(0.3, 0.7)), 0.0)
  expect_equal(
    nauc(result, xlim = c(-5, -1), ylim = c(0.3, 0.7)), 0.62649140813097662
  )
  expect_equal(nauc(result, xlim = c(9, 12), ylim = c(0.3, 0.7)), 1.0)
})

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

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

  expect_equal(naac(result), 1 - 0.65302073661957048)
  expect_equal(naac(result, xlim = c(-2, 2)), 1 - 0.65733232476226463)
  expect_equal(naac(result, ylim = c(0.3, 0.7)), 1 - 0.74700167012690062)
  expect_equal(naac(result, xlim = c(-25, -18), ylim = c(0.3, 0.7)), 1.0)
  expect_equal(
    naac(result, xlim = c(-5, -1), ylim = c(0.3, 0.7)), 1 - 0.62649140813097662
  )
  expect_equal(naac(result, xlim = c(9, 12), ylim = c(0.3, 0.7)), 0.0)
})