numerical_example <- function() {
  p <- matrix(c(0.024, 0.046, 0.021, 0.009,
                0.046, 0.096, 0.057, 0.001,
                0.021, 0.057, 0.216, 0.006,
                0.009, 0.001, 0.006, 0.384), nrow=4, byrow=TRUE)
  N <- 1000
  n <- p * N
  expected_N <- N
  expect_true(abs(sum(n) - expected_N) <= 0.001)
  expected_totals <- c(100, 200, 300, 400)
  rows <- rowSums(n)
  cols <- colSums(n)
  for (i in 1:4) {
    expect_true(abs(rows[i] - expected_totals[i]) <= 0.001)
    expect_true(abs(cols[i] - expected_totals[i]) <= 0.001)
  }
  n
}


set_up <- function() {
  delta <- 1.0e-6
  epsilon <- 100.0 * delta
  marginal_pi <- c(0.1, 0.2, 0.3, 0.4)
  r <- length(marginal_pi)
  kappa <- 0.28
  v <- matrix(0.0, nrow=r, ncol=r)
  v[1, 2] <- 2.3
  v[2, 1] <- v[1, 2]
  v[1, 3] <- 0.7
  v[3, 1] <- v[1, 3]
  v[2, 3] <- 0.95
  v[3, 2] <- v[2, 3]
  v[1, 4] <- 0.2250
  v[4, 1] <- v[1, 4]
  v[2, 4] <- 0.0125
  v[4, 2] <- v[2, 4]
  v[3, 4] <- 0.05
  v[4, 3] <- v[3, 4]

  pe <- sum(marginal_pi^2)
  v_tilde <- 1.0 + kappa * (1.0 - pe) / pe

  for (i in 1:r) {
    v[i, i] <- v_tilde
  }
  list(delta=delta, epsilon=epsilon, marginal_pi=marginal_pi,
       kappa=kappa, v=v)
}


test_that("derivative of pi wrt kappa is correct", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  r <- length(marginal_pi)
  v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  pi <- Schuster_compute_pi(marginal_pi, kappa, v)

  kappa <- kappa - delta
  v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  pi0 <- Schuster_compute_pi(marginal_pi, kappa, v0)

  kappa <- kappa + 2.0 * delta
  v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  pi1 <- Schuster_compute_pi(marginal_pi, kappa, v1)

  kappa <- kappa - delta
  Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)

  for (i in 1:r) {
    for (j in 1:r) {
      der <- ((pi[i, j] - pi0[i, j]) / delta + (pi1[i, j] - pi[i, j]) / delta) / 2.0
      deriv <- Schuster_derivative_pi_wrt_kappa(i, j, marginal_pi, kappa, v)
      expect_true(abs(der - deriv) <= epsilon, info=paste(i, j, der, deriv))
    }
  }
}
)


test_that("derivative wrt marginal_pi[k] is correct", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  r <- length(marginal_pi)
  pi <- Schuster_compute_pi(marginal_pi, kappa, v)
  v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  for (k in 1:(r - 1)) {
    marginal_pi[k] <- marginal_pi[k] - delta
    pi0 <- Schuster_compute_pi(marginal_pi, kappa, v)
    marginal_pi[k] <- marginal_pi[k] + 2.0 * delta
    pi1 <- Schuster_compute_pi(marginal_pi, kappa, v)
    marginal_pi[k] <- marginal_pi[k] - delta
    for (i in 1:r) {
      for (j in 1:r) {
        der <- ((pi[i, j] - pi0[i, j]) / delta + (pi1[i, j] - pi[i, j]) / delta) / 2.0
        deriv <- Schuster_derivative_pi_wrt_marginal_pi(i, j, k, marginal_pi, kappa, v)
        expect_true(abs(der - deriv) <= epsilon, info=paste(i, j, k, der, deriv))
      }
    }
  }
}
)


test_that("derivative pi[i, j] wrt v[i1, j1] is correct", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  r <- length(marginal_pi)
  pi <- Schuster_compute_pi(marginal_pi, kappa, v)
  for (i1 in 1:r) {
    for (j1 in 1:r) {
      v[i1, j1] <- v[i1, j1] - delta
      v[j1, i1] <- v[i1, j1]
      v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
      pi0 <- Schuster_compute_pi(marginal_pi, kappa, v0)

      v[i1, j1] <- v[i1, j1] + 2.0 * delta
      v[j1, i1] <- v[i1, j1]
      v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
      pi1 <- Schuster_compute_pi(marginal_pi, kappa, v1)

      v[i1, j1] <- v[i1, j1] - delta
      v[j1, i1] <- v[i1, j1]
      v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)

      for (i in 1:r) {
        for (j in 1:r) {
          der <- ((pi[i, j] - pi0[i, j]) / delta + (pi1[i, j] - pi[i, j]) / delta) / 2.0
          deriv <- Schuster_derivative_pi_wrt_v(i, j, i1, j1, marginal_pi, kappa, v)
          expect_true(abs(der - deriv) <= epsilon, info=paste(i, j, i1, j1, der, deriv))
        }
      }
    }
  }
}
)


test_that("derivative of log L wrt kappa is correct", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- nrow(n)

  v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  pi <- Schuster_compute_pi(marginal_pi, kappa, v)
  f <- log_likelihood(n, pi)

  kappa <- kappa - delta
  v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  pi0 <- Schuster_compute_pi(marginal_pi, kappa, v0)
  f0 <- log_likelihood(n, pi0)

  kappa <- kappa + 2.0 * delta
  v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  pi1 <- Schuster_compute_pi(marginal_pi, kappa, v1)
  f1 <- log_likelihood(n, pi1)

  kappa <- kappa - delta
  v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)

  der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
  deriv <- Schuster_derivative_log_l_wrt_kappa(n, marginal_pi, kappa, v)
  expect_true(abs(deriv - der) <= epsilon,
              info=paste("kappa", der, deriv))
}
)


test_that("derivative of log L wrt marginal_pi is correct", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- nrow(n)

  pi <- Schuster_compute_pi(marginal_pi, kappa, v)
  f <- log_likelihood(n, pi)
  for (k in 1:(r - 1)) {
    marginal_pi[k] <- marginal_pi[k] - delta
    marginal_pi[r] <- 1.0 - sum(marginal_pi[1:(r - 1)])
    pi <- Schuster_compute_pi(marginal_pi, kappa, v)
    f0 <- log_likelihood(n, pi)

    marginal_pi[k] <- marginal_pi[k] + 2.0 * delta
    marginal_pi[r] <- 1.0 - sum(marginal_pi[1:(r - 1)])
    pi <- Schuster_compute_pi(marginal_pi, kappa, v)
    f1 <- log_likelihood(n, pi)

    marginal_pi[k] <- marginal_pi[k] - delta
    marginal_pi[r] <- 1.0 - sum(marginal_pi[1:(r - 1)])

    der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
    deriv <- Schuster_derivative_log_l_wrt_marginal_pi(n, k, marginal_pi, kappa, v)
    expect_true(abs(deriv - der) <= epsilon,
                info=paste("marginal pi", k, der, deriv))
  }
}
)


test_that("constraint on v enforced correctly", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v_tilde <- Schuster_v_tilde(marginal_pi, kappa)
  v <- values$v
  r <- length(marginal_pi)
  marginal_pi[r] = 1.0 - sum(marginal_pi[1:(r - 1)])
  expect_true(0.0 < marginal_pi[r])
  v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  for (i in 1:(r - 1)) {
    expect_true(abs(v_tilde - v[i, i]) <= epsilon)
    expect_true(abs(1.0 - sum(v[i,] * marginal_pi)) <= epsilon,
                info=paste(i, sum(v[i, ] * marginal_pi)))
  }
}
)


test_that("derivative of pi[i, r] with enforce constraints is right", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v_tilde <- Schuster_v_tilde(marginal_pi, kappa)
  v <- values$v
  r <- length(marginal_pi)
  pe <- sum(marginal_pi^2)
  marginal_pi[r] = 1.0 - sum(marginal_pi[1:(r - 1)])
  expect_true(0.0 < marginal_pi[r])
  v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  pi <- Schuster_compute_pi(marginal_pi, kappa, v)
  for (i in 1:(r - 1)) {
    for (k in 1:(r - 1)) {
      der_v_tilde_der_pk <- -2.0 * kappa * (marginal_pi[k] - marginal_pi[r]) / pe^2

      marginal_pi[k] <- marginal_pi[k] - delta
      v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
      pi0 <- Schuster_compute_pi(marginal_pi, kappa, v)

      marginal_pi[k] <- marginal_pi[k] + 2.0 * delta
      v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
      pi1 <- Schuster_compute_pi(marginal_pi, kappa, v)

      marginal_pi[k] <- marginal_pi[k] - delta

      der <- ((pi[i, r] - pi0[i, r]) / delta + (pi1[i, r] - pi[i, r]) / delta) / 2.0
      deriv <- marginal_pi[i] * marginal_pi[r] * (-v[i, k] - der_v_tilde_der_pk * marginal_pi[i]) / marginal_pi[r]
      # deriv <- deriv + marginal_pi[i] * marginal_pi[r] * v[i, r] / marginal_pi[r]
      # deriv <- deriv + marginal_pi[i] * v[i, r] * (-1)

      if (i == k) {
        deriv <- v[k, r] * marginal_pi[r]
        deriv <- deriv + marginal_pi[k] * marginal_pi[r] * (-v[k, k] - marginal_pi[k] * der_v_tilde_der_pk) / marginal_pi[r]
        # deriv <- deriv + marginal_pi[k] * marginal_pi[r] * v[k, r] / marginal_pi[r]
        # deriv <- deriv - marginal_pi[k] * v[k, r]
      }
      expect_true(abs(deriv - der) <= epsilon, info=paste(i, k, der, deriv))
    }
  }
}
)


test_that("derivative v[i1, j1] wrt v[i1, j2] is correct", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)
  v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  pi <- Schuster_compute_pi(marginal_pi, kappa, v)
  for (i2 in 1:r) {
    for (j2 in 1:r) {
      v[i2, j2] <- v[i2, j2] - delta
      v[j2, i2] <- v[i2, j2]
      v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)

      v[i2, j2] <- v[i2, j2] + 2.0 * delta
      v[j2, i2] <- v[i2, j2]
      v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)

      v[i2, j2] <- v[i2, j2] - delta
      v[j2, i2] <- v[i2, j2]
      v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)

      for (i1 in 1:r) {
        for (j1 in 1:r) {
          der <- ((v[i1, j1] - v0[i1, j1]) / delta + (v1[i1, j1] - v[i1, j1]) / delta) / 2.0
          deriv <- Schuster_derivative_v_wrt_v(i1, j1, i2, j2, marginal_pi, kappa, v)
          expect_true(abs(der - deriv) <= epsilon,
                      info=paste("der_v_wrt_v", i1, j1, i2, j2, der, deriv))
        }
      }
    }
  }
}
)


test_that("derivative of log L wrt v is correct", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)
  v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  pi <- Schuster_compute_pi(marginal_pi, kappa, v)
  f <- log_likelihood(n, pi)
  for (i1 in 1:(r - 1)) {
    for (j1 in 1:(r - 1)) {

      v[i1, j1] <- v[i1, j1] - delta
      v[j1, i1] <- v[i1, j1]
      pi0 <- Schuster_compute_pi(marginal_pi, kappa, v)
      f0 <- log_likelihood(n, pi0)

      v[i1, j1] <- v[i1, j1] + 2.0 * delta
      v[j1, i1] <- v[i1, j1]
      pi1 <- Schuster_compute_pi(marginal_pi, kappa, v)
      f1 <- log_likelihood(n, pi1)

      v[i1, j1] <- v[i1, j1] - delta
      v[j1, i1] <- v[i1, j1]

      deriv <- Schuster_derivative_log_l_wrt_v(n, i1, j1, marginal_pi, kappa, v)
      der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
      expect_true(abs(deriv - der) <= epsilon, info=paste("v log l", i1, j1, der, deriv))
    }
  }
}
)


test_that("second order derivative pi wrt kappa^2 is 0", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)

  for (i in 1:r) {
    for (j in 1:r) {
      v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
      f <- Schuster_derivative_pi_wrt_kappa(i, j, marginal_pi, kappa, v)

      kappa <- kappa - delta
      v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
      f0 <- Schuster_derivative_pi_wrt_kappa(i, j, marginal_pi, kappa, v)

      kappa <- kappa + 2.0 * delta
      v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
      f1 <- Schuster_derivative_pi_wrt_kappa(i, j, marginal_pi, kappa, v)

      kappa <- kappa - delta

      der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
      deriv <- Schuster_second_deriv_pi_wrt_kappa_2(i, j, marginal_pi, kappa, v)
      expect_true(abs(der - deriv) < epsilon, info=paste(i, j, der, deriv))
    }
  }
}
)


test_that("second order derivative pi wrt kappa, v", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)
  for (i in 1:r) {
    for (j in 1:r) {
      for (i1 in 1:(r - 1)) {
        for (j1 in 1:(r - 1)) {
          v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
          pi <- Schuster_compute_pi(marginal_pi, kappa, v)
          f <- Schuster_derivative_pi_wrt_v(i, j, i1, j1, marginal_pi, kappa, v)

          kappa <- kappa - delta
          v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
          pi0 <- Schuster_compute_pi(marginal_pi, kappa, v0)
          f0 <- Schuster_derivative_pi_wrt_v(i, j, i1, j1, marginal_pi, kappa, v)

          kappa <- kappa + 2.0 * delta
          v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
          pi1 <- Schuster_compute_pi(marginal_pi, kappa, v1)
          f1 <- Schuster_derivative_pi_wrt_v(i, j, i1, j1, marginal_pi, kappa, v)

          kappa <- kappa - delta

          der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
          deriv <- Schuster_second_deriv_pi_wrt_kappa_v(i, j, i1, j1, marginal_pi, kappa, v)
          expect_true(abs(der - deriv) <= epsilon,
                      info=paste(i, j, i1, j1, der, deriv))
        }
      }
    }
  }
}
)


test_that("second order derivative pi wrt v^2", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)
  for (i in 1:r) {
    for (j in 1:r) {
      for (i1 in 1:(r - 1)) {
        for (j1 in 1:(r - 1)) {
          v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
          f <- Schuster_derivative_pi_wrt_v(i, j, i1, j1, marginal_pi, kappa, v)
          if (j1 == i1) next
          for (i2 in 1:(r - 1)) {
            for (j2 in 1:(r - 1)) {
              if (j2 == i2) next
              v[i2, j2] <- v[i2, j2] - delta
              v[j2, i2] <- v[i2, j2]
              v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
              f0 <- Schuster_derivative_pi_wrt_v(i, j, i1, j1, marginal_pi, kappa, v0)

              v[i2, j2] <- v[i2, j2] + 2.0 * delta
              v[j2, i1] <- v[i2, j2]
              v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
              f1 <- Schuster_derivative_pi_wrt_v(i, j, i1, j1, marginal_pi, kappa, v1)

              v[i2, j2] <- v[i2, j2] - delta
              v[j2, i2] <- v[i2, j2]
              v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)

              der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
              deriv <- Schuster_second_deriv_pi_wrt_v_2(i, j, i1, j1, i2, j2, marginal_pi, kappa, v)
              expect_true(abs(der - deriv) <= epsilon,
                          info=paste("der_v^2", i, j, i1, j2, der, deriv))
            }
          }
        }
      }
    }
  }
}
)


test_that("second order derivative pi wrt marginal_pi and kappa", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)
  for (i in 1:r) {
    for (j in 1:r) {
      for (k in 1:(r - 1)) {
        v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
        f <- Schuster_derivative_pi_wrt_marginal_pi(i, j, k, marginal_pi, kappa, v)

        kappa <- kappa - delta
        v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
        f0 <- Schuster_derivative_pi_wrt_marginal_pi(i, j, k, marginal_pi, kappa, v0)

        kappa <- kappa + 2.0 * delta
        v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
        f1 <- Schuster_derivative_pi_wrt_marginal_pi(i, j, k, marginal_pi, kappa, v1)

        kappa <- kappa - delta
        v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)

        der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
        deriv <- Schuster_second_deriv_pi_wrt_marginal_pi_kappa(i, j, k, marginal_pi, kappa, v)
        expect_true(abs(der - deriv) <= epsilon,
                    info=paste("der marginal_pi der kappa", i, j, k, der, deriv, abs(der - deriv)))
      }
    }
  }
}
)


test_that("second order derivative pi wrt marginal_pi^2", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)
  for (i in 1:r) {
    for (j in 1:r) {
      for (k in 1:(r - 1)) {
        v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
        f <- Schuster_derivative_pi_wrt_marginal_pi(i, j, k, marginal_pi, kappa, v)
        for (k2 in 1:(r - 1)) {
          marginal_pi[k2] <- marginal_pi[k2] - delta
          v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
          f0 <- Schuster_derivative_pi_wrt_marginal_pi(i, j, k, marginal_pi, kappa, v0)

          marginal_pi[k2] <- marginal_pi[k2] + 2.0 * delta
          v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
          f1 <- Schuster_derivative_pi_wrt_marginal_pi(i, j, k, marginal_pi, kappa, v1)

          marginal_pi[k2] <- marginal_pi[k2] - delta
          v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)

          der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
          deriv <- Schuster_second_deriv_pi_wrt_marginal_pi_2(i, j, k, k2, marginal_pi, kappa, v)
          expect_true(abs(der - deriv) <= epsilon,
                      info=paste("der marginal_pi^2", i, j, k, k2, der, deriv))
        }
      }
    }
  }
}
)


test_that("second order derivative pi wrt marginal_pi and v", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)
  # i & j idex pi
  for (i in 1:r) {
    for (j in 1:r) {
      # i1 & j1 index v
      for (i1 in 1:(r - 1)) {
        for (j1 in 1:(r - 1)) {
          if (j1 == i1) {
            next
          }
          # k indexes marginal_pi
          for (k in 1:(r - 1)) {
            v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
            f <- Schuster_derivative_pi_wrt_v(i, j, i1, j1, marginal_pi, kappa, v)

            marginal_pi[k] <- marginal_pi[k] - delta
            v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
            f0 <- Schuster_derivative_pi_wrt_v(i, j, i1, j1, marginal_pi, kappa, v0)

            marginal_pi[k] <- marginal_pi[k] + 2.0 * delta
            v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
            f1 <- Schuster_derivative_pi_wrt_v(i, j, i1, j1, marginal_pi, kappa, v1)

            marginal_pi[k] <- marginal_pi[k] - delta
            v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)

            der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
            deriv <- Schuster_second_deriv_pi_wrt_marginal_pi_v(i, j, k, i1, j1, marginal_pi, kappa, v)
            expect_true(abs(der - deriv) <= epsilon,
                        info=paste("der marginal_pi der v", i, j, k, i1, j1, der, deriv))
          }
        }
      }
    }
  }
}
)


test_that("second order log l wrt marginal_pi^2 is correct", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)
  deriv_mat <- Schuster_second_deriv_log_l_wrt_marginal_pi_2(n, marginal_pi, kappa, v)

  for (k in 1:(r - 1)) {
    for (k1 in (1:(r - 1))) {
      v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
      pi <- Schuster_compute_pi(marginal_pi, kappa, v)
      f <- Schuster_derivative_log_l_wrt_marginal_pi(n, k1, marginal_pi, kappa, v)

      marginal_pi[k] <- marginal_pi[k] - delta
      v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
      pi <- Schuster_compute_pi(marginal_pi, kappa, v0)
      f0 <- Schuster_derivative_log_l_wrt_marginal_pi(n, k1, marginal_pi, kappa, v0)

      marginal_pi[k] <- marginal_pi[k] + 2.0 * delta
      v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
      pi <- Schuster_compute_pi(marginal_pi, kappa, v1)
      f1 <- Schuster_derivative_log_l_wrt_marginal_pi(n, k1, marginal_pi, kappa, v1)

      marginal_pi[k] <- marginal_pi[k] - delta

      der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
      expect_true(abs(deriv_mat[k, k1] - deriv_mat[k1, k]) <= epsilon)
      expect_true(abs(der - deriv_mat[k, k1]) <= epsilon,
                  info=paste("der marginal_pi^2", k, k1, der, deriv))
    }
  }
}
)


test_that("second order log l wrt marginal_pi and v is correct", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)
  deriv_mat <- Schuster_second_deriv_log_l_wrt_marginal_pi_v(n, marginal_pi, kappa, v)

  for (k in 1:(r - 1)) {
    v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
    pi <- Schuster_compute_pi(marginal_pi, kappa, v)

    index <- 1
    for (i1 in (1:(r - 2))) {
      for (j1 in (i1 + 1): (r - 1)) {
        for (k1 in 1:(r-1)) {
          f <- Schuster_derivative_log_l_wrt_v(n, i1, j1, marginal_pi, kappa, v)

          marginal_pi[k] <- marginal_pi[k] - delta
          v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
          pi <- Schuster_compute_pi(marginal_pi, kappa, v0)
          f0 <- Schuster_derivative_log_l_wrt_v(n, i1, j1, marginal_pi, kappa, v0)

          marginal_pi[k] <- marginal_pi[k] + 2.0 * delta
          v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
          pi <- Schuster_compute_pi(marginal_pi, kappa, v1)
          f1 <- Schuster_derivative_log_l_wrt_v(n, i1, j1, marginal_pi, kappa, v1)

          marginal_pi[k] <- marginal_pi[k] - delta

          der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
          expect_true(abs(der - deriv_mat[k, index]) <= epsilon,
                    info=paste("der marginal_pi^2", k, k1, der, deriv_mat[k, index]))
        }
        index <- index + 1
      }
    }
  }
}
)


test_that("second order log l wrt kappa^2 is correct", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)
  deriv <- Schuster_second_deriv_log_l_wrt_kappa_2(n, marginal_pi, kappa, v)
  v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  pi <- Schuster_compute_pi(marginal_pi, kappa, v)
  f <- Schuster_derivative_log_l_wrt_kappa(n, marginal_pi, kappa, v)

  kappa <- kappa - delta
  v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  pi <- Schuster_compute_pi(marginal_pi, kappa, v0)
  f0 <- Schuster_derivative_log_l_wrt_kappa(n, marginal_pi, kappa, v0)

  kappa <- kappa + 2.0 * delta
  v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  pi <- Schuster_compute_pi(marginal_pi, kappa, v1)
  f1 <- Schuster_derivative_log_l_wrt_kappa(n, marginal_pi, kappa, v1)

  kappa <- kappa - delta

  der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
  expect_true(abs(der - deriv) <= epsilon,
              info=paste("der marginal_pi^2", k, der, deriv))
}
)


test_that("second order log l wrt v^2 is correct", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)

  deriv_mat <- Schuster_second_deriv_log_l_wrt_v_2(n, marginal_pi, kappa, v)

  index <- 1
  for (i1 in 1:(r - 2)) {
    for (j1 in (i1 + 1):(r - 1)) {
      index2 <- 1
      for (i2 in 1:(r - 2)) {
        for (j2 in (i2 + 1): (r - 1)) {
          v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
          f <- Schuster_derivative_log_l_wrt_v(n, i1, j1, marginal_pi, kappa, v)

          v[i2, j2] <- v[i2, j2] - delta
          v[j2, i2] <- v[i2, j2]
          v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
          f0 <- Schuster_derivative_log_l_wrt_v(n, i1, j1, marginal_pi, kappa, v0)

          v[i2, j2] <- v[i2, j2] + 2.0 * delta
          v[j2, i2] <- v[i2, j2]
          v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
          f1 <- Schuster_derivative_log_l_wrt_v(n, i1, j1, marginal_pi, kappa, v1)

          v[i2, j2] <- v[i2, j2] - delta
          v[j2, i2] <- v[i2, j2]

          der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
          expect_true(abs(deriv_mat[index, index2] - deriv_mat[index2, index]) <= epsilon)
          expect_true(abs(der - deriv_mat[index, index2]) <= epsilon,
                      info=paste("der marginal_pi^2", i1, j1, i2, j2, der, deriv_mat[index, index2]))
          index2 <- index2 + 1
        }
      }
      index <- index + 1
    }
  }
}
)


test_that("second order log l wrt kappa and v is correct", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)

  deriv_vec <- Schuster_second_deriv_log_l_wrt_kappa_v(n, marginal_pi, kappa, v)
  index <- 1
  for (i1 in 1:(r - 2)) {
    for (j1 in (i1 + 1): (r - 1)) {
      v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
      f <- Schuster_derivative_log_l_wrt_v(n, i1, j1, marginal_pi, kappa, v)

      kappa <- kappa - delta
      v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
      f0 <- Schuster_derivative_log_l_wrt_v(n, i1, j1, marginal_pi, kappa, v0)

      kappa <- kappa + 2.0 * delta
      v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
      f1 <- Schuster_derivative_log_l_wrt_v(n, i1, j1, marginal_pi, kappa, v1)

      kappa <- kappa - delta

      der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
      expect_true(abs(der - deriv_vec[index]) <= epsilon,
                  info=paste(i1, j1, der, deriv_vec[index]))
      index <- index + 1
    }
  }
}
)


test_that("second order log l wrt marginal_pi kappa is correct", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)

  v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  deriv_vec <- Schuster_second_deriv_log_l_wrt_marginal_pi_kappa(n, marginal_pi, kappa, v)

  for (k in 1:(r - 1)) {
    v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
    f <- Schuster_derivative_log_l_wrt_marginal_pi(n, k, marginal_pi, kappa, v)

    kappa <- kappa - delta
    v0 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
    f0 <- Schuster_derivative_log_l_wrt_marginal_pi(n, k, marginal_pi, kappa, v0)

    kappa <- kappa + 2.0 * delta
    v1 <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
    f1 <- Schuster_derivative_log_l_wrt_marginal_pi(n, k, marginal_pi, kappa, v1)

    kappa <- kappa - delta

    der <- ((f - f0) / delta + (f1 - f) / delta) / 2.0
    expect_true(abs(der - deriv_vec[k]) <= epsilon,
                info=paste("der marginal_pi k", k, der, deriv_vec[k]))
  }
}
)


test_that("hessian is correct",  {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)

  hessian <- Schuster_hessian(n, marginal_pi, kappa, v)
  expect_equal(nrow(hessian), ncol(hessian))
  for (i in 1:nrow(hessian)) {
    for (j in 1:ncol(hessian)) {
      expect_true(abs(hessian[i, j] - hessian[j, i]) <= epsilon)
    }
  }

  mat_marginal_pi <- Schuster_second_deriv_log_l_wrt_marginal_pi_2(n, marginal_pi, kappa, v)
  for (k in 1:(r - 1)) {
    for (k1 in 1:(r - 1)) {
      expect_true(abs(mat_marginal_pi[k, k1] - hessian[k, k1]) <= epsilon)
    }
  }

  vec_marginal_pi_kappa <- Schuster_second_deriv_log_l_wrt_marginal_pi_kappa(n, marginal_pi, kappa, v)
  for (i in 1:length(vec_marginal_pi_kappa)) {
    expect_true(abs(hessian[i, r] - vec_marginal_pi_kappa[i]) <= epsilon)
    expect_true(abs(hessian[r, i] - vec_marginal_pi_kappa[i]) <= epsilon)
  }

  deriv_kappa <- Schuster_second_deriv_log_l_wrt_kappa_2(n, marginal_pi, kappa, v)
  expect_true(abs(hessian[r, r] - deriv_kappa) <= epsilon,
              info=paste("deriv kappa", r, deriv_kappa, hessian[r, r]))

  mat_marginal_pi_v <- Schuster_second_deriv_log_l_wrt_marginal_pi_v(n, marginal_pi, kappa, v)
  for (i in 1:nrow(mat_marginal_pi_v)) {
    for (j in 1:ncol(mat_marginal_pi_v)) {
      expect_true(abs(hessian[i, r + j] - mat_marginal_pi_v[i, j]) <= epsilon,
                  info=paste("marginal_pi, v", i, j, hessian[i, r + j], mat_marginal_pi_v[i, j]))
      expect_true(abs(hessian[r + j, i] - mat_marginal_pi_v[i, j]) <= epsilon,
                  info=paste("marginal_pi, v", i, j, hessian[r + j, i], mat_marginal_pi_v[j, i]))
    }
  }

  vec_kappa_v <- Schuster_second_deriv_log_l_wrt_kappa_v(n, marginal_pi, kappa, v)
  for (i in 1:length(vec_kappa_v)) {
    expect_true(abs(hessian[r, r + i] - vec_kappa_v[i]) <= epsilon)
    expect_true(abs(hessian[r + i, r] - vec_kappa_v[i]) <= epsilon)
  }

  mat_v_v <- Schuster_second_deriv_log_l_wrt_v_2(n, marginal_pi, kappa, v)
  for (i in 1:nrow(mat_v_v)) {
    for (j in 1:ncol(mat_v_v)) {
      expect_true(abs(hessian[r + i, r + j] - mat_v_v[i, j]) <= epsilon)
      expect_true(abs(hessian[r + j, r + i] - mat_v_v[i, j]) <= epsilon)
    }
  }
}
)


test_that("Newton-Raphson steps improve log(likelihoo)", {
  values <- set_up()
  delta <- values$delta
  epsilon <- values$epsilon
  marginal_pi <- values$marginal_pi
  kappa <- values$kappa
  v <- values$v
  n <- numerical_example()
  r <- length(marginal_pi)
  v <- Schuster_enforce_constraints_on_v(marginal_pi, kappa, v)
  pi <- Schuster_compute_pi(marginal_pi, kappa, v)
  pi <- pi / sum(pi)
  original_logL <- log_likelihood(n, pi)
  result <- Schuster_symmetric_rater_agreement_model(n, max_iter=200, verbose=FALSE)
  marginal_pi <- result$marginal_pi
  kappa <- result$kappa
  v <- result$v
  pi <- Schuster_compute_pi(marginal_pi, kappa, v)
  new_logL <- log_likelihood(n, pi)
  expect_true(new_logL > original_logL)
}
)


test_that("actual budget data gives similar answer to article", {
  n <- budget_actual
  result <- Schuster_symmetric_rater_agreement_model(n, max_iter=250, verbose=FALSE)
  expect_true(abs(result$g_squared - 3.105) <= 0.2,
              info=paste("G^2", 3.105, result$g_squared))
  expect_true(abs(result$kappa - 0.308) <= 0.025,
              info=paste("kappa", 0.308, result$kappa))
  expect_true(abs(result$v[1,1] - 1.59) <= 0.05,
              info=paste("v-tilde", 1.59, result$v[1, 1]))
  expect_equal(5, result$df)
}
)


test_that("expected budget data gives similar answer to article", {
  n <- budget_expected
  r <- nrow(budget_expected)
  N <- sum(n)
  # can't match values in paper. This has better fit than reported there
  result <- Schuster_symmetric_rater_agreement_model(n, max_iter=250, verbose=FALSE)
  expect_true(abs(result$g_squared - 9.384) <= 0.0005,
              info=paste("G^2", 9.384, result$g_squared))
  expect_true(abs(result$kappa - 0.319) <= 0.0005,
              info=paste("kappa", 0.319, result$kappa))
  expect_true(abs(result$v[1,1] - 1.592) <= 0.0005,
              info=paste("v-tilde", 1.592, result$v[1, 1]))
  expect_equal(5, result$df)

  logL <- result$logL

  # checking fit of numbers in article.
  marginal_pi <- vector("double", r)
  marginal_pi[1] <- 36.9 / N
  marginal_pi[2] <- 36.9 / N
  marginal_pi[3] <- 57.2 / N
  kappa <- 0.254
  v_tilde <- 1.47
  v <- result$v
  for (i in 1:r) {
    v[i, i] <- v_tilde
  }
  pi <- Schuster_compute_pi(marginal_pi, kappa, v)
  pi <- pi / sum(pi)
  # assert that this solution is better than the one published
  g2 <- likelihood_ratio_chisq(n, pi)
  expect_true(result$g_squared < g2)
  x2 <- pearson_chisq(n, pi)
  expect_true(result$chisq < x2)
}
)