test_that("Test shrinkage estimator", {
  # Parameters
  nSamples <- 500
  pTrue <- 2

  # True moments
  trueMean <- c(0, 0)
  trueSigma <- matrix(c(3, 2, 2, 2), nrow = 2)
  chol_trueSigma <- chol(trueSigma)

  # we run 100 shrinkage estimators
  lambdas <- rep(0, 100)
  for (i in seq(100)) {
    rr <- replicate(nSamples, trueMean) +
      t(chol_trueSigma) %*% matrix(stats::rnorm(pTrue * nSamples), nrow = pTrue, ncol = nSamples)
    # Estimate mean and covariance from samples
    mean_est <- rowMeans(rr)
    Sigma_est <- cov(t(rr))

    mean_est
    trueMean

    Sigma_est

    rr_centered <- rr - replicate(nSamples, mean_est)
    rowMeans(rr_centered)

    x <- t(rr)
    # x <- t(rr_centered)

    lambdas[i] <- schaferStrimmer_cov(x)$lambda_star
  }
  # The average over 100 runs must be within a certain range
  mean_lambdas <- mean(lambdas)

  # We expect this to be lower than the 99.99% quantile over 200 simulations of this scenario
  expect_lte(mean_lambdas, 0.005191634)

  # We expect this to be greater than the 0.01% quantile over 200 simulations of this scenario
  expect_gte(mean_lambdas, 0.004843328)
})

test_that("Test behavior for p=1", {
  # Parameters
  nSamples <- 500
  pTrue <- 1

  samples <- matrix(stats::rnorm(n = nSamples, sd = 2), ncol = pTrue)

  exp_sample_squared <- matrix(mean(samples^2), 1, 1)

  # The function should return E[x^2] and a warning
  expect_warning(expect_equal(schaferStrimmer_cov(samples)$shrink_cov, exp_sample_squared))
})