R/ACI_est.R

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

#' Simple ACP Estimation
#'
#' Estimates average coverage probability using auxiliary randomization when \eqn{\lambda \ge 1/2}.
#'
#' @param chi a scalar half-length value \eqn{\chi}.
#' @param lam a scalar shrinkage factor \eqn{\lambda}.
#' @param xvec normalized observed outcome vector.
#' @param Z a vector of the length \code{length(xvec)},
#' to be used as auxiliary random variables; \code{Z} can be left unspecified.
#'
#' @return estimated average coverage probability value.
#' @export
#'
#' @examples xvec <- stats::rnorm(50)
#' gnZ(1.5, 0.7, xvec)
#' Z <- stats::rnorm(50)
#' gnZ(1.5, 0.7, xvec, Z)
gnZ <- function(chi, lam, xvec, Z = NULL){

  n <- length(xvec)

  if(is.null(Z)) Z <- stats::rnorm(n)

  new_rv <- ((lam - 1) / lam) * xvec + sqrt(1 - ((lam - 1) / lam)^2) * Z
  ind <- (new_rv <= chi / lam) & (new_rv >= -chi / lam)

  return(sum(ind) / n)
}

#' Simple ACP Estimation with Lower Variance
#'
#' Estimates average coverage probability using conditional mean of the
#' estimator with auxiliary randomization when \eqn{\lambda \ge 1/2}.
#'
#' @inheritParams gnZ
#'
#' @return estimated average coverage probability value.
#' @export
#'
#' @examples xvec <- stats::rnorm(50)
#' gnZ_cond(1.5, 0.7, xvec)
#' gnZ_cond(1, 0.5, xvec)
gnZ_cond <- function(chi, lam, xvec){

  n <- length(xvec)

  if(lam == 1/2){

    psi_i <- as.numeric((xvec <= 2 * chi) & (xvec >= -2 * chi))

  }else{

    sfun <- function(lam){

      (1 - ((lam - 1)/lam)^2)^(-1/2)
    }

    psi_i <- stats::pnorm(sfun(lam) * ((1 - lam)/lam * xvec + chi/lam)) -
      stats::pnorm(sfun(lam) * ((1 - lam)/lam * xvec - chi/lam))
  }

  res <- mean(psi_i)
  return(res)

}

#' ACP Estimation by Polynomial Approximation
#'
#' Estimates average coverage probability using polynomial approximation
#' of the ACP function when \eqn{\lambda < 1/2}.
#'
#' @param chi a scalar half-length value \eqn{\chi}.
#' @param lam a vector of shrinkage factors \eqn{\lambda}.
#' @param xvec normalized observed outcome vector,
#' corresponding to \eqn{z_i} in the paper.
#' @param Jn the order of polynomial approximation.
#'
#' @return vector of estimated average coverage probability values,
#' corresponding to each value of \eqn{\lambda} in \code{lam}.
#' @export
#'
#' @examples xvec <- rnorm(50)
#' gnB(1, c(0.3, 0.4), xvec, 2)
gnB <- function(chi, lam, xvec, Jn){

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

  coef <- cfun(chi, lam, Jn)
  x_pol <- H_cal(xvec, Jn)
  res <- as.numeric(colMeans(x_pol) %*% coef) # lamlen-dim vector

  return(res)

}


#' ACP Estimation by Polynomial Approximation under Normality Assumption
#'
#' This is function is a version of \code{gnB} function where
#' higher moments of the true mean vector are estimated assuming the true mean vector is
#' normally distributed.
#'
#' \code{Jn} terms are estimated in the same way as \code{gnB}, while \code{Jn2 - Jn} terms
#' are estimated under the normality assumption.
#'
#' @inheritParams gnB
#' @param Jn the order of polynomial approximation to be used for
#' series estimation.
#' @param Jn2 the total order of polynomial approximation.
#' @param orcind oracle specification to be used; when \code{orcind = "m2"}, the value of
#' true second moment of the true mean vector is provided in \code{m2}; when \code{orcind = "th"},
#' the entire true mean vector is provided; the values other than \code{"default"} are used for
#' simulations.
#' @param m2 the value of the true second moment of the true mean vector;
#' used when \code{orcind = "m2"}.
#' @param th_vec the true mean vector; used when \code{orcind = "th"}.
#'
#' @return vector of estimated average coverage probability values,
#' corresponding to each value of \eqn{\lambda} in \code{lam}.
#' @export
#'
#' @examples th_vec <- stats::rnorm(50)
#' xvec <- stats::rnorm(50, th_vec)
#' gnB_norm(1, c(0.3, 0.4), xvec, 2, 4)
#' gnB_norm(1, c(0.3, 0.4), xvec, 0, 4)
#' m2 <- mean(th_vec^2)
#' gnB_norm(1, c(0.3, 0.4), xvec, 2, 4, orcind = "m2", m2)
#' gnB_norm(1, c(0.3, 0.4), xvec, 2, 4, orcind = "th", NULL, th_vec)
gnB_norm <- function(chi, lam, xvec, Jn, Jn2,
                     orcind = c("default", "m2", "theta"), m2, th_vec){

  # lam can be a vector // Jn2 > Jn // Jn >= 2

  orcind <- match.arg(orcind)

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

  coef_all <- cfun(chi, lam, Jn2)
  coef <- coef_all[1:(Jn + 1), ,drop = F]
  coef_rem <- coef_all[(Jn + 2):(Jn2 + 1), ,drop = F]

  x_pol <- H_cal(xvec, Jn)
  sum <- as.numeric(colMeans(x_pol) %*% coef) # lamlen-dim vector

  if(orcind == "m2" & !missing(m2)){

    mu2est = m2

  }else if(orcind == "default"){

    mu2est <- colMeans(x_pol)[3]

  }

  Jn_rem <- (Jn + 1):Jn2
  Jn_rem_len <- length(Jn_rem)

  for (i in 1:Jn_rem_len){

    coef_i <- coef_rem[i]
    j <- Jn_rem[i]
    coef_mom <- abs(hpol(j)[1]) * (j %% 2 == 0)

    if(orcind == "theta"){

      sum <- sum + coef_i * mean(th_vec^j)

    }else{

      sum <- sum + coef_i * coef_mom * mu2est^((j / 2) * (j %% 2 == 0))
    }

  }

  res <- sum
  return(res)

}


#' ACP Estimation by Polynomial Approximation under Normality Assumption
#'
#' This is function is a version of \code{gnB} function where
#' higher moments of the true mean vector are estimated assuming the true mean vector is
#' normally distributed.
#'
#' \code{Jn} terms are estimated in the same way as \code{gnB}, while \code{Jn2 - Jn} terms
#' are estimated under the normality assumption.
#'
#' @inheritParams gnB
#' @param Jn the order of polynomial approximation to be used for
#' series estimation.
#' @param Jn2 the total order of polynomial approximation.
#' @param orcind oracle specification to be used; when \code{orcind = "m2"}, the value of
#' true second moment of the true mean vector is provided in \code{m2}; when \code{orcind = "th"},
#' the entire true mean vector is provided; the values other than \code{"default"} are used for
#' simulations.
#' @param m2 the value of the true second moment of the true mean vector;
#' used when \code{orcind = "m2"}.
#' @param th_vec the true mean vector; used when \code{orcind = "th"}.
#'
#' @return vector of estimated average coverage probability values,
#' corresponding to each value of \eqn{\lambda} in \code{lam}.
#' @export
#'
#' @examples th_vec <- stats::rnorm(50)
#' xvec <- stats::rnorm(50, th_vec)
#' gnB_norm2(1, c(0.3, 0.4), xvec, 2, 4)
#' gnB_norm2(1, c(0.3, 0.4), xvec, 0, 4)
#' m2 <- mean(th_vec^2)
#' gnB_norm2(1, c(0.3, 0.4), xvec, 2, 4, orcind = "m2", m2)
#' gnB_norm2(1, c(0.3, 0.4), xvec, 2, 4, orcind = "th", NULL, th_vec)
gnB_norm2 <- function(chi, lam, xvec, Jn, Jn2,
                     orcind = c("default", "m2", "theta"), m2, th_vec){

  # lam can be a vector // Jn2 > Jn // Jn >= 2

  orcind <- match.arg(orcind)

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

  coef_all <- cfun(chi, lam, Jn2)
  coef <- coef_all[1:(Jn + 1), ,drop = F]
  coef_rem <- coef_all[(Jn + 2):(Jn2 + 1), ,drop = F]

  x_pol <- H_cal(xvec, Jn)
  sum <- as.numeric(colMeans(x_pol) %*% coef) # lamlen-dim vector

  if(orcind == "m2" & !missing(m2)){

    mu2est = m2

  }else if(orcind == "default"){

    mu2est <- colMeans(x_pol)[3]
    mu2est <- max(0, mu2est) # mu2est sometimes takes a negative value

  }

  Jn_rem <- (Jn + 1):Jn2
  Jn_rem_len <- length(Jn_rem)

  for (i in 1:Jn_rem_len){

    coef_i <- coef_rem[i]
    j <- Jn_rem[i]

    if(orcind == "theta"){

      sum <- sum + coef_i * mean(th_vec^j)

    }else{

      sum <- sum + coef_i * norm_moment(j, 0, sqrt(mu2est))
    }

  }

  res <- sum
  return(res)

}