R/solver.R

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

#' Half-length Solver
#'
#' Solves for the half-length corresponding to a shrinkage factor.
#'
#' @param lam a scalar shrinkage factor \eqn{\lambda}.
#' @param xvec normalized observed outcome vector,
#' corresponding to \eqn{z_i} in the paper.
#' @param Jn the order of polynomial approximation.
#' @param alpha desired level of non-coverage probability.
#'
#' @return the estimated value of the half-length \eqn{\chi} which gives the correct coverage
#' given the value of \code{lam}.
#' @export
#'
#' @examples thvec <- rnorm(50)
#' xvec <- rnorm(50, thvec)
#' chilam(0.6, xvec, 0.05)
#' chilam(0.4, xvec, 0.05, 2)
chilam <- function(lam, xvec, alpha, Jn){

  chirange <- c(0, stats::qnorm(1 - alpha/2))
  n <- length(xvec)

  if(lam >= 1/2){

    # A function returning gnZ_cond - (1 - alpha)
    gnZ_0 <- function(chi, lam, xvec, alpha){

      res <- gnZ_cond(chi, lam, xvec) - (1 - alpha)
      return(res)
    }

    if(gnZ_0(chirange[2], lam, xvec, alpha) <= 0){

      res <- chirange[2]

    }else{

      sol = stats::uniroot(f = gnZ_0, interval = chirange, lam = lam, xvec = xvec, alpha = alpha,
                           extendInt = "upX")
      res <- sol$root
    }

  }else{

    # A function returning gnB - (1 - alpha)
    gnB_0 <- function(chi, lam, xvec, Jn, alpha){

      res <- gnB(chi, lam, xvec, Jn) - (1 - alpha)
      return(res)
    }

    if(gnB_0(chirange[2], lam, xvec, Jn, alpha) <= 0){

      res <- chirange[2]

    }else{

      sol <- stats::uniroot(f = gnB_0, interval = chirange, lam = lam, xvec = xvec, Jn = Jn, alpha = alpha,
                            extendInt = "upX")
      res <- sol$root
    }
  }

  return(res)
}


#' Oracle Half-length Solver
#'
#' Solves for the half-length corresponding to a shrinkage factor,
#' when the true mean vector is known.
#'
#' This function was created for the purpose of simulations.
#'
#' @param lam a scalar shrinkage factor \eqn{\lambda}.
#' @param th_vec a normalized mean vector.
#' @param alpha desired level of non-coverage probability.
#'
#' @return the value of the half-length \eqn{\chi} which gives the correct coverage
#' given the value of \code{lam}.
#' @export
#'
#' @examples th_vec <- rnorm(50)
#' chilam_orc(0.4, th_vec, 0.05)
chilam_orc <- function(lam, th_vec, alpha){

  # returns cov_prob - (1 - alpha)
  cov_prob0 <- function(chival, lamval, th_vec, alpha){

    res <- cov_prob(chival, lamval, th_vec) - (1 - alpha)
    return(res)
  }

  chirange <- c(0, stats::qnorm(1 - alpha/2))

  sol <- stats::uniroot(f = cov_prob0, interval = chirange, lamval = lam, th_vec = th_vec, alpha = alpha,
                        extendInt = "upX")

  res <- sol$root


  return(res)
}


#'  Half-length Solver under Normality Assumption
#'
#'  Solves for the half-length corresponding to a shrinkage factor, using the function
#'  \code{gnB_norm}.
#'
#' @inheritParams chilam
#' @param Jn the order of polynomial approximation to be used for
#' series estimation.
#' @param Jn2 the total order of polynomial approximation.
#'
#' @return the estimated value of the half-length \eqn{\chi} which gives the correct coverage
#' given the value of \code{lam}.
#' @export
#'
#' @examples thvec <- rnorm(50)
#' xvec <- rnorm(50, thvec)
#' chilam_norm(0.6, xvec, 0.05, 2, 4)
#' chilam_norm(0.4, xvec, 0.05, 2, 4)
chilam_norm <- function(lam, xvec, alpha, Jn, Jn2){

  chirange <- c(0, stats::qnorm(1 - alpha/2))
  n <- length(xvec)

  if(lam >= 1/2){

    # A function returning gnZ_cond - (1 - alpha)
    gnZ_0 <- function(chi, lam, xvec, alpha){

      res <- gnZ_cond(chi, lam, xvec) - (1 - alpha)
      return(res)
    }

    if(gnZ_0(chirange[2], lam, xvec, alpha) <= 0){

      res <- chirange[2]

    }else{

      sol = stats::uniroot(f = gnZ_0, interval = chirange, lam = lam, xvec = xvec, alpha = alpha,
                           extendInt = "upX")
      res <- sol$root
    }

  }else{

    # A function returning gnB_norm - (1 - alpha)
    gnB_norm_0 <- function(chi, lam, xvec, Jn, Jn2, alpha){

      res <- gnB_norm2(chi, lam, xvec, Jn, Jn2) - (1 - alpha)
      return(res)
    }

    if(gnB_norm_0(chirange[2], lam, xvec, Jn, Jn2, alpha) <= 0){

      res <- chirange[2]

    }else{

      sol <- stats::uniroot(f = gnB_norm_0, interval = chirange, lam = lam, xvec = xvec,
                            Jn = Jn, Jn2 = Jn2, alpha = alpha, extendInt = "upX")
      res <- sol$root
    }
  }

  return(res)
}