R/pairwise_conf_region.R

In carter-allen/mvnormR: Multivariate normal plotting

#' Pairwise confidence ellipse plotter for mean vector
#'
#' This function allows you to plot (1-alpha)*100% confidence ellipses for pairs of variables
#' @param X data matrix
#' @param var1 An integer corresponding to the column of variable 1
#' @param var2 An integer corresponding to the column of variable 2
#' @export
#'

pairwise_conf_region <- function(X, var1 = 1, var2 = 2, alpha = 0.05)
{
  n = nrow(X)
  p = ncol(X)
  X_sub = X[, c(var1, var2)]
  X_bar = matrix(colMeans(X_sub), nrow = 2, byrow = TRUE)
  S = cov(X_sub)
  S.inv = solve(S)
  ret_df = NULL
  for (i in seq(min(X_sub[, 1]), max(X_sub[, 1]), length.out = 100))
    {
      for (j in seq(min(X_sub[, 2]), max(X_sub[, 2]), length.out = 100))
        {
          mu_ij = matrix(c(i, j), nrow = 2, byrow = TRUE)
          q1 = n * t(X_bar - mu_ij) %*% S.inv %*% (X_bar - mu_ij)
          q2 = (p * (n - 1))/(n - p) * qf(1 - alpha, p, n - p)
          if (q1 <= q2)
          {
            ret_df = rbind(ret_df, t(mu_ij))
          }
        }
      }
      colnames(ret_df) = c(paste("mu", var1, sep = ""), paste("mu", var2, sep = ""))
      return(as.data.frame(ret_df))
}