test_that("elicitPower 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")

  # compound gamma example
  # central credible intervals with tertiles
  Z <- eglm(scenarios = X, link = "log.inverse",
    known.dispersion = FALSE, obs.name = "compoundPoisson",
    power = 1.1)

  muhat <- 1/2
  r1 <- 10
  s1 <- 2
  p1 <- 1.5
  w.size <- 1000
  rp1 <- r1*(2 - p1)/(2*muhat^(2 - p1))
  v1 <- r1*muhat^p1/(w.size*s1)

  mean.lower.bounds <- eglm:::qgt(c(1/3, 0.05), muhat, sqrt(v1), s1, s1)
  median.q0 <- eglm:::qUG(1/2, s1, rp1)

  Z <- elicitDispersion(Z, mu.hat = muhat, w = w.size,
                        lower.bound1 = mean.lower.bounds[1],
                        lower.bound2 = mean.lower.bounds[2])
  Z <- elicitPower(Z, median = median.q0)

  if (requireNamespace("tweedie")) {

    # Monte Carlo simulation
    set.seed(101)
    n.mc <- 2000
    lambda <- rgamma(n.mc, shape = s1/2, rate = r1/2)
    p.zero <- sample.means <- numeric(n.mc)
    for (i in 1:n.mc) {
      y <- tweedie::rtweedie(w.size, mu = muhat, phi = 1/lambda[i], power = p1)
      p.zero[i] <- sum(y == 0)/w.size
      sample.means[i] <- mean(y)
    }

    hist(p.zero, freq = FALSE, breaks = 20)
    curve(eglm:::dUG(x, Z$obs.model$s, Z$obs.model$rp), col = 'red', add = TRUE, n = 501)

    capture <- capture.output(result <- elicitPower(Z, median = median.q0,
                          evalCDF0 = c(0.8, 0.9, 0.95)))
    out <- strsplit(capture[4], " ")[[1]]
    nums <- as.numeric(out)
    fracs0 <- nums[!is.na(nums)]
    Fn0 <- ecdf(p.zero)
    expect_equal(fracs0,
                 Fn0(c(0.8, 0.9, 0.95)),
                 tolerance = 0.05)

  }
})