R/tester.R

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

#' Test of Monotonicity
#'
#' Tests whether the average coverage probability function
#' estimated by the series approximation is monotone in \eqn{\chi}.
#'
#' @inheritParams gnB
#' @param alpha the desired level of non-coverage probability
#' @param chilen the number of grid points of \eqn{\chi} to be used in the test;
#' the default is 50.
#'
#' @return a vector of logical values, each value taking \code{TRUE} if monotone,
#' with each value corresponding to each value of \eqn{\lambda} in \code{lam}.
#' @export
#'
#' @examples xvec <- rnorm(50)
#' lam <- c(0.3, 0.35, 0.4, 0.45)
#' mono_chi(lam, xvec, 2, 0.05)
mono_chi <- function(lam, xvec, Jn, alpha, chilen = 50){

  chivec <- seq(from = 0, to = stats::qnorm(1 - alpha / 2), length.out = chilen)
  lamlen <- length(lam)
  gres <- matrix(0, nrow = chilen, ncol = lamlen)

  for(i in 1:chilen){

    chi <- chivec[i]
    gres[i, ] <- gnB(chi, lam, xvec, Jn)
  }

  res <- numeric(lamlen)

  for(l in 1:lamlen){
    res[l] <- !is.unsorted(gres[ , l])
  }

  return(res)
}

#' Bias Upper Bound Estimation
#'
#' Calculates the bias upper bound estimate for the average coverage probability
#' function estimator using the series approximation.
#'
#' Even if an even \code{Jn} is used in the actual estimation, we can plug in
#' \code{Jn + 1} in this function, as described in the paper.
#'
#' @inheritParams gnB
#' @param Jn an odd number.
#'
#' @return the value of the bias upper bound estimate.
#' @export
#'
#' @examples xvec <- rnorm(50)
#' lam <- c(0.3, 0.35, 0.4, 0.45)
#' gnB_bias(lam, xvec, 3)
gnB_bias <- function(lam, xvec, Jn){

  # Bias UB estimator for gnB: Jn should be an odd number
  # lam can be a vector

  n <- length(xvec)
  K <- (1 - lam) / lam

  const1 <- sqrt(2 / pi) / factorial(Jn + 1) *  K^(Jn + 1)

  hroots <- polyroot(hpol(Jn + 1))
  const2 <- max(abs(EQL::hermite(hroots, Jn, prob = TRUE) * exp(-hroots^2 / 2)))

  est <- mean(EQL::hermite(xvec, Jn + 1, prob = TRUE))

  res <- const1 * const2 * est

  return(res)

}

#' Test of Large Bias
#'
#' Tests whether the bias of the average coverage probability function
#' estimated by the series approximation is too large.
#'
#' The bias is deemed large if the upper bound bias estimate computed by
#' \code{gnB_bias} is large than \code{tol}.
#'
#' @inheritParams mono_chi
#' @param tol a threshold value to test for the large bias;
#' the default is \code{alpha/2}.
#'
#' @return a vector of logical values, each value taking \code{TRUE} if estimated
#' bias is smaller than \code{tol}, with each value corresponding to
#' each value of \eqn{\lambda} in \code{lam}.
#' @export
#'
#' @examples xvec <- rnorm(50)
#' lam <- c(0.3, 0.35, 0.4, 0.45)
#' bias_test(lam, xvec, 2, 0.05)
bias_test <- function(lam, xvec, Jn, alpha, tol = alpha/2){

  bias_bd <- gnB_bias(lam, xvec, Jn)
  res <- bias_bd < tol

  return(res)

}

#' Shrinkage Factor Adjustment
#'
#' Determines the minimum shrinkage factor to be used using
#' tests for monotonicity and large bias. If no other value of
#' shrinkage factor passes the tests, returns \eqn{\lambda = 1/2}.
#'
#' The tests are conducted using functions \code{mono_chi} and \code{bias_test}
#' in the package \code{OptACI}.
#'
#' @inheritParams mono_chi
#' @inheritParams bias_test
#' @param lam a grid of shrinkage factors, smaller than 1/2.
#'
#' @return the value of adjusted minimum value of shrinkage factor.
#' @export
#'
#' @examples xvec <- rnorm(50)
#' lam <- c(0.3, 0.35, 0.4, 0.45)
#' lam_adj(lam, xvec, 2, 0.05)
lam_adj <- function(lam, xvec, Jn, alpha, chilen = 50, tol = alpha/2){

  monoind <- mono_chi(lam, xvec, Jn, alpha, chilen)
  biasind <- bias_test(lam, xvec, Jn, tol)
  lamind <- as.logical(monoind * biasind)
  lammin_new <- min(c(lam[lamind], 0.5))

  return(lammin_new)
}