R/conf_int_conditional.R

In adviksh/winfer:

#' Conditinoal Confidence Interval Using Truncated Normality
#' @export
#' @inheritParams ci_docs
#'
ci_conditional_unbiased <- function(estimates, st_errs, conf_level = 0.95) {

  # Defense -----------------------------------------------------------------
  check_estimates(estimates)
  check_st_errs(st_errs   = st_errs,
                estimates = estimates)
  check_conf_level(conf_level)

  # Function Body -----------------------------------------------------------
  half_alpha <- (1 - conf_level) / 2
  tail_tiles <- c(half_alpha, 1 - half_alpha)

  ordering <- order(estimates, decreasing = TRUE)

  estimates_sorted <- estimates[ordering]
  st_errs_sorted   <- st_errs[ordering]

  ci_conditional_unbiased_single(estimates_sorted = estimates_sorted,
                                 st_errs_sorted   = st_errs_sorted,
                                 tail_tiles = tail_tiles)
}

ci_conditional_unbiased_single <- function(estimates_sorted,
                                           st_errs_sorted,
                                           tail_tiles) {

  if (is.unsorted(rev(estimates_sorted))) {
    rlang::abort("estimates_unsorted is not sorted in decreasing order")
  }

  tp <- conditional_truncation_points(estimates_sorted)

  purrr::map_dbl(tail_tiles, mu_alpha_conditional,
                 y  = estimates_sorted[1],
                 se = st_errs_sorted[1],
                 trunc_lo = min(tp),
                 trunc_hi = max(tp))
}

mu_alpha_conditional <- function(alpha, y, se, trunc_lo = -Inf, trunc_hi = Inf) {

  ur_result <- stats::uniroot(f         = tcdf_minus_left_tail,
                              interval  = c(y +  3 * se, y - 3 * se),
                              y         = y,
                              se        = se,
                              alpha     = alpha,
                              trunc_lo  = trunc_lo,
                              trunc_hi  = trunc_hi,
                              extendInt = "yes")
  ur_result$root
}

tcdf_minus_left_tail <- function(mu, se, y, alpha,
                                 trunc_lo = -Inf,
                                 trunc_hi = Inf) {

  p <- TruncatedNormal::ptmvnorm(q     = y,
                                 mu    = mu,
                                 sigma = se,
                                 lb    = trunc_lo,
                                 ub    = trunc_hi,
                                 log   = FALSE,
                                 type  = "mc",
                                 B     = 1e5)

  if (is.na(p)) browser()

  p - (1 - alpha)
}