R/ACI_orc.R

In koohyun-kwon/OptACI: Methods to Construct Optimal Average Coverage Confidence Intervals

#' Average Coverage Probability
#'
#' Function for the average coverage probability calculation under a known value of the mean vector.
#'
#' This function was created for the purpose of simulations.
#'
#' @param chi a scalar half-length value \eqn{\chi}.
#' @param lam a scalar shrinkage factor \eqn{\lambda}.
#' @param th_vec a normalized mean vector.
#'
#' @return average coverage probability.
#' @export
#'
#' @examples cov_prob(1.5, 0.5, seq(from = -1, to = 1, by = 0.01))
cov_prob <- function(chi, lam, th_vec){

  probs <- stats::pnorm((1 - lam)/lam * th_vec + chi / lam) - stats::pnorm((1 - lam)/lam * th_vec - chi / lam)
  res <- mean(probs)

  return(res)
}


#' Average Coverage Probability: Truncated
#'
#' Function for the average coverage probability calculation under a known value of the mean vector,
#' when truncated series approximation is used.
#'
#' This function was created for the purpose of simulations.
#'
#' @inheritParams gnB
#' @param th_vec a normalized mean vector.
#'
#' @return vector of estimated average coverage probability values,
#' corresponding to each value of \eqn{\lambda} in \code{lam}.
#' @export
#'
#' @examples cov_prob_tr(1, 0.5, seq(from = -1, to = 1, by = 0.01), 4)
cov_prob_tr <- function(chi, lam, th_vec, Jn){ # lam can be a vector

  n <- length(th_vec)
  lamlen <- length(lam)

  sum <- numeric(lamlen)

  coef <- cfun(chi, lam, Jn)

  for(j in 0:Jn){

    newsum <- coef[j + 1, ] * mean(th_vec^j)
    sum <- sum + newsum
  }

  res <- sum
  return(res)

}

#' Unbiased Half-oracle Average Coverage Probability Estimator
#'
#' Estimates average coverage probability with the same series approximation
#' procedure as \code{gnB}, and then removes the truncation bias.
#'
#' @inheritParams cov_prob_tr
#' @param lam a scalar shrinkage factor \eqn{\lambda}.
#' @param xvec normalized observed outcome vector,
#' corresponding to \eqn{z_i} in the paper.
#'
#' @return average coverage probability.
#' @export
#'
#' @examples th_vec <- stats::rnorm(50)
#' xvec <- stats::rnorm(50, th_vec)
#' cov_prob_nobias(1, 0.5, th_vec, 2, xvec)
cov_prob_nobias <- function(chi, lam, th_vec, Jn, xvec){

    rem <- cov_prob(chi, lam, th_vec) - cov_prob_tr(chi, lam, th_vec, Jn)
    res <- gnB(chi, lam, xvec, Jn) + rem
    return(res)
}