test_that("elicitDispersion matches density", {

  skip_on_cran()

  # scenarios matrix: two scenarios
  X <- matrix(c(1, 1, 0, 1), nrow = 2)
  rownames(X) <- c("scenario1", "scenario2")
  colnames(X) <- c("covariate1", "covariate2")

  # gamma example
  # central credible intervals with tertiles
  Z <- eglm(scenarios = X, link = "log.inverse", known.dispersion = FALSE, obs.name = "gamma")
  w.size <- 1000
  Z <- elicitDispersion(Z, mu.hat = 9, lower.bound1 = 8.5, lower.bound2 = 7, w = w.size)

  # Monte Carlo simulation
  set.seed(123)
  n.means <- 100000
  sample.means <- numeric(n.means)
  for (i in 1:n.means) {
    lambda <- rgamma(w.size, shape = Z$obs.model$s/2, rate = Z$obs.model$r/2)
    # lambda is shape parameter for gamma observation model
    # rate for observation model depends also on mu.hat
    r.obs <- lambda/9
    sample.means[i] <- mean(rgamma(w.size, shape = lambda, rate = r.obs))
  }
  q <- c(0.05, 1/3, 0.5, 2/3, 0.95)
  abline(v = quantile(sample.means, probs = q), col = 'red') # Monte Carlo estimates

  # numerical tests
  capture <- capture.output(result <- elicitDispersion(Z, mu.hat = 9, evalCDF = c(7, 8.5, 9), w = w.size))
  out <- strsplit(capture, " ")[[5]]
  nums <- as.numeric(out)
  fracs <- nums[!is.na(nums)]
  Fn <- ecdf(sample.means)
  abline(v = quantile(sample.means, q), col = 'red', lty = 2)
  expect_equal(fracs,
               Fn(c(7, 8.5, 9)),
               tolerance = 0.06)
})