test_that("sampling from univariate normal", {
  
  # Generate 1e4 samples from univariate Gaussian
  params <- list(mean=42, sd=1)
  distr <- "gaussian"
  n <- 1e4
  samples <- .distr_sample(params, distr, n)
  
  # Compute empirical mean and sd 
  sam_mean <- mean(samples)
  sam_sd <- sd(samples)
  
  # Check how close empirical values are to the truth
  m <- abs(sam_mean-42)/42
  s <- abs(sam_sd-1)
  
  expect_equal(m < 2e-3, TRUE)
  expect_equal(s < 4e-2, TRUE)
})

test_that("sampling from univariate nbinom", {
  
  # Generate 1e4 samples from negative binomial (size, prob)
  params <- list(size=12,prob=0.8)
  distr <- "nbinom"
  n <- 1e4
  samples <- .distr_sample(params, distr, n)
  
  # Compute empirical mean
  sam_mean <- mean(samples)
  true_mean <- params$size*(1-params$prob)/params$prob
  
  # Check how close empirical values are to the truth
  m <- abs(sam_mean-true_mean)/true_mean

  expect_equal(m < 3e-2, TRUE)
  
  # Generate 1e4 samples from negative binomial (size, mu)
  params <- list(size=12,mu=true_mean)
  distr <- "nbinom"
  n <- 1e4
  samples <- .distr_sample(params, distr, n)
  
  # Compute empirical mean
  sam_mean <- mean(samples)
  
  # Check how close empirical values are to the truth
  m <- abs(sam_mean-params$mu)/params$mu
  
  expect_equal(m < 3e-2, TRUE)
  
  # Check if it returns an error when all 3 parameters are specified
  params <- list(size=12,mu=true_mean,prob=0.8)
  distr <- "nbinom"
  n <- 1e4
  expect_error(.distr_sample(params, distr, n))
  
  # Check if it returns an error when size is not specified
  params <- list(mu=true_mean,prob=0.8)
  distr <- "nbinom"
  n <- 1e4
  expect_error(.distr_sample(params, distr, n))
  
  
  
})

test_that("sampling from univariate poisson", {
  
  # Generate 1e4 samples from poisson
  params <- list(lambda=10)
  distr <- "poisson"
  n <- 1e4
  samples <- .distr_sample(params, distr, n)
  
  # Compute empirical mean
  sam_mean <- mean(samples)
  
  # Check how close empirical values are to the truth
  m <- abs(sam_mean-10)/10
  
  expect_equal(m < 3e-2, TRUE)
})

test_that("sampling from multivariate normal", {

  # Generate 1e4 samples from bivariate normal
  mu=c(10,10)
  Sigma= matrix(c(1,0.7,0.7,1),nrow = 2)
  n <- 1e4
  samples <- .MVN_sample(n, mu, Sigma)
  
  # Compute empirical mean
  sam_mean <- colMeans(samples)
  
  # Check how close empirical values are to the truth
  m <- abs(sam_mean-10)/10
  
  expect_equal(all(m < 8e-3), TRUE)
})

test_that("MVN density works", {
  
  # Create 3x3 covariance matrix
  L <- matrix(0,nrow=3,ncol=3)
  L[lower.tri(L,diag=TRUE)] <- c(0.9,0.8,0.5,0.9,0.2,0.6)
  Sigma <- L%*%t(L)
  
  # create mean vector
  mu <- c(0,1,-1)
  
  # matrix where to evaluate the MVN
  xx <- matrix(c(0,2,1,
                 2,3,4,
                 0.5,0.5,0.5,
                 0,1,-1), ncol=3,byrow=TRUE)
  
  res <- .MVN_density(x=xx,mu=mu,Sigma=Sigma)
  
  true_val <- c(8.742644e-04, 1.375497e-11, 3.739985e-03, 1.306453e-01)
  expect_equal(res,true_val,tolerance = "3e")
  
  # Check if block-evaluation works
  xx <- matrix(runif(3*1e4),ncol=3,byrow=TRUE)
  
  res_chuncks <- .MVN_density(x=xx,mu=mu,Sigma=Sigma)
  res_all <- .MVN_density(x=xx,mu=mu,Sigma=Sigma,max_size_x = 1e4)
  
  expect_equal(res_chuncks,res_all)
  
})