## check that predictive distributions hold what they promise,
## i.e. that they describe the sum of n new data points.

set.seed(42343)

## precision at which reference and tested must match
eps <- 1e-2

## percentiles to test for
p_quants <- seq(0.1, 0.9, by = 0.1)

## number of samples used for sampling method
Nsamp <- 1e5

## define the different test cases
beta <- mixbeta(c(1, 11, 4))
betaMix <- mixbeta(c(0.8, 11, 4), c(0.2, 1, 1))

gamma <- mixgamma(c(1, 65, 50))
gammaMix <- mixgamma(rob = c(0.5, 8, 0.5), inf = c(0.5, 9, 2), param = "ms")

norm <- mixnorm(c(1, 0, 0.5), sigma = 1)
normMix <- mixnorm(c(0.2, 0, 2), c(0.8, 2, 1), sigma = 1)

n <- 25

preddist_cmp <- function(mix, n, n_rng, N = Nsamp, qntls = p_quants, stat = c("sum", "mean"), Teps = eps) {
  skip_on_cran()

  ## sample for each draw a single hyper-parameter which is then
  ## used n times in the rng function to return n samples from the
  ## sampling distribution
  stat <- match.arg(stat)
  test <- replicate(N, n_rng(rmix(mix, 1)))
  if (stat == "sum") test_stat <- colSums(test)
  if (stat == "mean") test_stat <- colMeans(test)

  quants_stest <- quantile(test_stat, qntls)
  ## note: in particular for the discrete/counting distributions,
  ## sampling gives better estimates
  quants_sref <- qmix(preddist(mix, n = n), qntls)
  res_sum <- abs(quants_sref - quants_stest)
  ## note: errors are accumulating with n, hence to check the mean,
  ## we scale eps with n
  if (stat == "sum") expect_true(all(res_sum < n * Teps))
  if (stat == "mean") expect_true(all(res_sum < Teps))

  quants_test <- quantile(test[1, ], qntls)
  quants_ref <- qmix(preddist(mix, n = 1), qntls)
  res <- abs(quants_ref - quants_test)
  expect_true(all(res < Teps))

  if (inherits(mix, "betaMix")) {
    ## specifically test BetaBinomial (which is from Stan)
    predmix <- preddist(mix, n = 30)
    dens <- dmix(predmix, 0:30)
    cdens <- cumsum(dens)
    pc <- pmix(predmix, 0:30)
    upc <- pmix(predmix, 0:30, FALSE)
    expect_true(all(abs(cdens - pc) < Teps))
    expect_true(all(abs(1 - cdens - upc) < Teps))
  }

  ## test that output length of the predictive is the same as the
  ## input data vector, even for negative values for values outside
  ## the valid support of the postive only distribtions
  cprob_ltest <- pmix(preddist(mix, n = 1), c(-1, quants_ref))
  expect_true(length(cprob_ltest) == length(c(-1, quants_ref)))
}


test_that("Predictive for a beta evaluates correctly (binary)", {
  preddist_cmp(beta, n, Curry(rbinom, n = n, size = 1))
})
test_that("Predictive for a beta mixture evaluates correctly (binary)", {
  preddist_cmp(betaMix, n, Curry(rbinom, n = n, size = 1))
})

test_that("Predictive for a gamma evaluates correctly (poisson)", {
  preddist_cmp(gamma, n, Curry(rpois, n = n))
})
test_that("Predictive for a gamma mixture evaluates correctly (poisson)", {
  preddist_cmp(gammaMix, n, Curry(rpois, n = n), Teps = 1E-1)
})

test_that("Predictive for a normal evaluates correctly (normal)", {
  preddist_cmp(norm, n, Curry(rnorm, n = n, sd = sigma(norm)), stat = "mean")
})
test_that("Predictive for a normal mixture evaluates correctly (normal)", {
  preddist_cmp(normMix, n, Curry(rnorm, n = n, sd = sigma(normMix)), stat = "mean", Teps = 1E-1)
})