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)
  
})