## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width = 7,
  fig.height = 4.4,
  fig.align = "center"
)

## ----load---------------------------------------------------------------------
library(rootWishartHD)

## ----dsb----------------------------------------------------------------------
dsb <- function(p, m_df, q_df) {
  s   <- q_df
  mC  <- 0.5 * (abs(m_df - s) - 1)
  df2 <- p - m_df + q_df
  nC  <- 0.5 * (abs(df2 - s) - 1)
  list(s = s, m = mC, n = nC)
}

## ----family-map, echo = FALSE-------------------------------------------------
fam <- data.frame(
  `function` = c(
    "doubleWishart()", "doubleWishart_scaled()", "doubleWishart_log()",
    "doubleWishart_robustPfaffians()", "doubleWishart_pvalue()",
    "doubleWishart_qalpha()", "singleWishart_cdf()",
    "singleWishart_log()", "singleWishart_pvalue()",
    "singleWishart_qalpha()"),
  ensemble = c(rep("double", 6), rep("single", 4)),
  returns  = c("CDF", "CDF", "log CDF / log SF", "log CDF / log SF",
               "p-value", "critical value", "CDF", "log CDF / log SF",
               "p-value", "critical value"),
  role = c(
    "default CDF interface; routes to scaled backend when needed",
    "scaled Pfaffian CDF backend for ill-conditioned or large cases",
    "tail-safe log survival and log CDF",
    "auto Pfaffian backend, stress testing",
    "applied fast-then-exact upper-tail p-value",
    "invert the tail for a cutoff",
    "CDF, single-Wishart largest eigenvalue",
    "tail-safe log survival, single-Wishart",
    "applied p-value, single-Wishart",
    "applied critical value, single-Wishart"),
  check.names = FALSE, stringsAsFactors = FALSE)
knitr::kable(fam, caption = "The main numerical functions exported by rootWishartHD.")

## ----backends-----------------------------------------------------------------
pr <- dsb(20, 14, 10)                    # s = 10
theta <- seq(0.5, 0.97, length.out = 6)
v0 <- doubleWishart(theta, s = pr$s, m = pr$m, n = pr$n,
                    type = "double", verbose = FALSE)
v2 <- doubleWishart_scaled(theta, s = pr$s, m = pr$m, n = pr$n,
                            type = "double", verbose = FALSE)
data.frame(theta, doubleWishart = v0, scaled = v2,
           abs_diff = abs(v0 - v2))

## ----robust-------------------------------------------------------------------
th <- c(0.85, 0.92, 0.97)
a <- doubleWishart_log(th, s = pr$s, m = pr$m, n = pr$n,
                                    type = "arbitrary", tail = "upper",
                                    verbose = FALSE)
b <- doubleWishart_robustPfaffians(th, s = pr$s, m = pr$m, n = pr$n,
                                   type = "arbitrary", tail = "upper",
                                   verbose = FALSE)
data.frame(theta = th, log = a, robustPfaffians = b,
           abs_diff = abs(a - b))

## ----pvalue-------------------------------------------------------------------
pv <- doubleWishart_pvalue(c(0.85, 0.92, 0.97),
                           s = pr$s, m = pr$m, n = pr$n,
                           input = "theta", verbose = FALSE)
data.frame(theta = c(0.85, 0.92, 0.97),
           p_value = as.numeric(pv),
           method  = attr(pv, "method"))

## ----qalpha-------------------------------------------------------------------
crit <- sapply(c(0.05, 0.01), function(a)
  doubleWishart_qalpha(a, s = pr$s, m = pr$m, n = pr$n,
                       return = "theta", verbose = FALSE))
back <- doubleWishart_pvalue(crit, s = pr$s, m = pr$m, n = pr$n,
                             input = "theta", verbose = FALSE)
data.frame(alpha = c(0.05, 0.01), theta_crit = crit,
           p_back = as.numeric(back),
           rel_err = abs(as.numeric(back) - c(0.05, 0.01)) / c(0.05, 0.01))

## ----single-cdf---------------------------------------------------------------
nmin <- 5L; nmax <- 10L
x <- c(10, 20, 30, 40)
F2  <- singleWishart_cdf(x, n_min = nmin, n_max = nmax,
                                   type = "double", verbose = FALSE)
lsf <- singleWishart_log(x, n_min = nmin, n_max = nmax,
                                      type = "arbitrary", tail = "upper",
                                      verbose = FALSE)
data.frame(x, CDF = F2, CDF_from_logSF = -expm1(lsf), logSF = lsf)

## ----single-applied-----------------------------------------------------------
pv <- singleWishart_pvalue(x, n_min = nmin, n_max = nmax, verbose = FALSE)
xa <- sapply(c(0.05, 0.01), function(a)
  singleWishart_qalpha(a, n_min = nmin, n_max = nmax, verbose = FALSE))
xb <- singleWishart_pvalue(xa, n_min = nmin, n_max = nmax, verbose = FALSE)
list(
  pvalues = data.frame(x, p = as.numeric(pv), method = attr(pv, "method")),
  qalpha  = data.frame(alpha = c(0.05, 0.01), x_crit = xa,
                       p_back = as.numeric(xb))
)

## ----decision, echo = FALSE---------------------------------------------------
dec <- data.frame(
  goal = c(
    "CDF in the body of the distribution",
    "Upper-tail probability that may be tiny",
    "Observed statistic -> p-value",
    "Significance level -> critical value",
    "Suspected Pfaffian-backend failure at large s",
    "Single-Wishart largest eigenvalue tail"),
  call = c(
    "doubleWishart()",
    "doubleWishart_log(tail = \"upper\")",
    "doubleWishart_pvalue()",
    "doubleWishart_qalpha()",
    "doubleWishart_robustPfaffians()",
    "singleWishart_pvalue() / singleWishart_qalpha()"),
  check.names = FALSE, stringsAsFactors = FALSE)
knitr::kable(dec, caption = "A quick decision guide for the function family.")