R/aux_fcn.R

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

#' Coefficient Calculation Helper Function for Series Expansion
#'
#' Calculates \eqn{h_j} function in the paper.
#'
#' @param x a value where \eqn{H_j} function is evaluated at; it can be a vector as well.
#' @param j the corresponding order.
#'
#' @return the coefficient value, with length of \code{len(x)}.
#' @export
#'
#' @examples hfun(1, 2)
#' hfun(c(-1, 1), 3)
hfun <- function(x, j){ # x can be a vector

  if(j == 0){

    res <- stats::pnorm(x)

  }else{

    res <- stats::dnorm(x) * (-1)^(j - 1) * EQL::hermite(x, j - 1)
  }

  return(res)
}

#' Coefficient Calculation for Series Expansion
#'
#' Calculates \eqn{c_j} function in the paper.
#'
#' @param chi a scalar half-length value \eqn{\chi}.
#' @param lam a vector of shrinkage factors \eqn{\lambda}.
#' @param Jn the order of polynomial approximation.
#'
#' @return a matrix with the row length \code{Jn + 1}
#'  and the column length \code{length(lam)}.
#' @export
#'
#' @examples cfun(1, c(0.2, 0.4), 2)
cfun <- function(chi, lam, Jn){ # lam can be a vector

  lamlen <- length(lam)
  res <- matrix(0, nrow = Jn + 1, ncol = lamlen)

  for(j in 0 : Jn){

    res[j + 1, ] <- (1 / factorial(j)) * ((1 - lam) / lam)^j *
      (hfun(chi / lam, j) - hfun(-chi / lam, j))
  }

  return(res)
}

#' Hermite Polynomial
#'
#' @param xvec vector of observations
#' @param J order of the polynomial
#'
#' @return matrix with the row length \code{length(xvec)}
#' and the column length \code{J + 1}.
#' @export
#'
#' @examples H_cal(c(1, 2, 3), 2)
H_cal <- function(xvec, J){

  n = length(xvec)
  xvec_pol = matrix(0, n, J + 1)
  coef_mat <- matrix(0, J + 1, J + 1)
  for(i in 1:(J + 1)){
    j = i - 1
    coefs = orthopolynom::hermite.h.polynomials(j)[[j + 1]]
    coefs = as.numeric(coefs)
    scale_const <- sqrt(2)^(-j) * sqrt(2)^(-(0:j))
    coefs = coefs * scale_const
    coefs = c(coefs, rep(0, J + 1 - i))
    coef_mat[i, ] <- coefs
    xvec_pol[, i] <- xvec^j
  }

  coef_mat <- t(coef_mat)

  res <- xvec_pol %*% coef_mat

  return(res)

}

#' Probabilist's Hermite polynomial Coefficients
#'
#' Returns probabilist's Hermite polynomial coefficients as a vector.
#'
#' @param J the order of polynomials.
#'
#' @return a vector with the length \code{J + 1}.
#' @export
#'
#' @examples hpol(5)
hpol <- function(J){

  res <- orthopolynom::hermite.he.polynomials(J)
  res <- as.numeric(res[[J + 1]])

  return(res)

}

#' Cutoff Value Finder for the AKP Procedure
#'
#' @param n a sample size.
#' @param cutoff_len a grid length for cutoff values to be tested.
#' @param m_len a grid length for spread values to be tested.
#' @param w_len a grid length for skewness values to be tested.
#' @param m_max the largest spread value; the default is 1.
#' @param m_min the smallest spread value, the default is 0.1.
#' @param cutoff_max the largest cutoff value; the default is 0.1.
#' @param cutoff_min the smallest cutoff value; the default is 0.01.
#' @param w_min the smallest skeweness value; the default is 0.1.
#' @param alpha the desired average non-coverage probability; the default is 0.05.
#' @param eps fixed error term; added for a simulation purpose.
#'
#' @return a \code{m_len * w_len * cutoff_len} dimensional vector, containing average
#' coverage probability for each specification.
#' @export
#'
#' @examples ebci_cutoff(500, 2, 2, 2)
#' ebci_cutoff(500, 3, 4, 1)
ebci_cutoff <- function(n, cutoff_len, m_len, w_len, m_max = 1, m_min = 0.1, cutoff_max = 0.1,
                        cutoff_min = 0.01, w_min = 0.1, alpha = 0.05, eps = NULL){

  cutoff_vec <- seq(from = cutoff_min, to = cutoff_max, length.out = cutoff_len)
  m_vec <- seq(from = m_min, to = m_max, length.out = m_len)

  if(w_len == 1){

    w_vec <- 0.5

  }else{

    w_vec <- seq(from = w_min, to = 0.5, length.out = w_len)
  }

  res <- array(0, c(m_len, w_len, cutoff_len))

  for(i in 1:m_len){
    for(j in 1:w_len){
      for(k in 1:cutoff_len){

        m <- m_vec[i]
        w <- w_vec[j]
        cutoff <- cutoff_vec[k]

        n1 <- as.integer(n * w)
        n2 <- n - n1
        th <- c(rep(sqrt(m), n1), rep(-sqrt(m), n2))
        th <- th - mean(th)

        if(is.null(eps)){

          y <- stats::rnorm(n, th)

        }else{

          y <- th + eps
        }

        CIres <- ACI_AKP(y, rep(1, n), alpha, cutoff)

        CI_l <- CIres$CI_l
        CI_u <- CIres$CI_u

        cover_l <- th > CI_l
        cover_u <- th < CI_u
        cover_ind <- cover_l * cover_u

        res[i, j, k] <- mean(cover_ind)
      }
    }
  }

  res_vec <- as.vector(apply(res, 3, c))
  return(res_vec)

}

#' Double Factorial
#'
#' @param n a even or odd number
#'
#' @return the value of n!!
#' @export
#'
#' @examples d.factorial(6)
#' d.factorial(7)
d.factorial <- function(n){

  if(n %% 2 == 0){

    res <- 2^{n/2} * factorial(n/2)

  }else{

    k <- (n + 1)/2
    res <- factorial(2 * k) / (2^k * factorial(k))
  }

  return(res)
}

#' Gaussian Moments
#'
#' @param n the power of moment to be calculated
#' @param mu mean of the normal distribution
#' @param sig standard deviation of the normal distribution
#'
#' @return nth moment of N(\code{mu}, \code{sig}^2)
#' @export
#'
#' @examples norm_moment(6, 1, 2)
norm_moment <- function(n, mu, sig){

  # from https://www.johndcook.com/blog/2012/11/06/general-formula-for-normal-moments/

  sumend <- floor(n/2)
  sum <- 0

  for(j in 0:sumend){

    term1 <- factorial(n) / (factorial(n - 2 * j) * factorial(2 * j))
    term2 <- d.factorial(2*j - 1)
    term3 <- sig^(2 * j) * mu^(n - 2 * j)

    sum = sum + term1 * term2 * term3
  }

  return(sum)
}