The rootWishartHD function family: backends, tails, p-values and critical values

Scope

The companion vignette “Exact and log-scale tail probabilities for Roy’s largest root” concentrates on the two work-horse functions doubleWishart() and doubleWishart_log(). rootWishartHD exports ten numerical functions (plus rootWishartHD_mpfr_enabled() for backend diagnostics), and this vignette is a guided tour of the rest: the alternative CDF backends, the robust-Pfaffian stress tool, the applied p-value and critical-value wrappers, and the complete single-Wishart family. The goal is to answer a practical question – which function should I call, and what does it cost?

library(rootWishartHD)

The double-Wishart functions use Chiani’s (s, m, n) parametrization. For a MANOVA-type problem with dimension p, hypothesis d.f. q_df and error d.f. encoded by m_df, we use the same mapping as the first vignette:

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)
}

A map of the family

Two design ideas recur throughout and are worth stating once:

• Log-scale tails. A double-precision CDF saturates: in moderate dimensions 1 - F underflows to 0, so the upper tail is unrecoverable. The *_log() functions return log F or log(1 - F) and stay finite arbitrarily far into the tail.
• Fast-then-exact. The applied wrappers (*_pvalue(), *_qalpha()) first try the cheap double CDF and only fall back to the slow arbitrary-precision (multiprecision) path for the points the double path cannot resolve. For large dimensions the wrappers bypass fragile fast-double passes and go directly to log-scale evaluation. Every result is tagged with the method actually used.

Double-Wishart: CDF backends

doubleWishart() and doubleWishart_scaled() both return the CDF F(theta) of theta = lambda / (1 + lambda). The first is the default user interface. It uses the historical fast path for small well-conditioned cases and routes to the scaled backend for large dimensions or arbitrary-precision runs. The second exposes the scaled Pfaffian backend directly for diagnostics.

On well-conditioned cases the two agree to roughly 1e-6 (the spread lives in the deep lower tail, where the values themselves are tiny):

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))
#>   theta doubleWishart       scaled     abs_diff
#> 1 0.500  1.208198e-13 1.207937e-13 2.604226e-17
#> 2 0.594  2.293764e-09 2.294117e-09 3.527681e-13
#> 3 0.688  5.311518e-06 5.312692e-06 1.174693e-09
#> 4 0.782  2.052932e-03 2.052813e-03 1.198231e-07
#> 5 0.876  1.345084e-01 1.345201e-01 1.171028e-05
#> 6 0.970  9.435361e-01 9.435115e-01 2.454725e-05

For everyday use, call doubleWishart(). Use doubleWishart_scaled() only when you want to inspect or force the scaled backend.

Saturation and log functions

The reason the tail functions exist is visible the moment the CDF approaches 1. The numbers below are measured at s = 40 (p = 80); they are shown as a static table because the exact evaluation takes ~40 s for three points and would slow the vignette down.

The robust-Pfaffian stress tool

doubleWishart_robustPfaffians() is a hardened entry point for large s: it fixes pf_method = "auto" (try several Pfaffian backends and take a conservative value), turns on adaptive precision, and always returns on the log scale. For well-behaved cases it is numerically identical to doubleWishart_log(); its value is as a fall-back when a single Pfaffian backend produces CDF > 1 overshoot artifacts that would collapse the tail.

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))
#>   theta         log robustPfaffians abs_diff
#> 1  0.85 -0.05257891     -0.05257891        0
#> 2  0.92 -0.63044036     -0.63044036        0
#> 3  0.97 -2.87380177     -2.87380177        0

Use doubleWishart_log() by default; reach for doubleWishart_robustPfaffians() only when you suspect a backend failure at large s.

Applied workflow: p-values and critical values

Most users never call the CDF directly. The two applied wrappers cover the common tasks.

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"))
#>   theta    p_value      method
#> 1  0.85 0.94877912 fast_double
#> 2  0.92 0.53234470 fast_double
#> 3  0.97 0.05648845 fast_double

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))
#>   alpha theta_crit      p_back      rel_err
#> 1  0.05  0.9712648 0.049999742 5.161086e-06
#> 2  0.01  0.9831993 0.009996837 3.163059e-04

Single-Wishart family

The single-Wishart functions describe the largest eigenvalue of a real Wishart matrix. The argument x is on the eigenvalue scale and the CDF is P(lambda_max <= x); the parameters are n_min = min(p, n) and n_max = max(p, n). The four exported functions mirror the double-Wishart applied set.

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)
#>    x         CDF CDF_from_logSF        logSF
#> 1 10 0.001675232    0.001675232 -0.001676636
#> 2 20 0.402072650    0.402072650 -0.514286021
#> 3 30 0.923259044    0.923259044 -2.567319737
#> 4 40 0.996243545    0.996243545 -5.584279582

CDF (fast double) and CDF_from_logSF (1 - exp(logSF), exact) agree, which validates the double path in the body of the distribution. The p-value and critical-value wrappers behave exactly like their double-Wishart counterparts, including the method tagging and a built-in sanity check that rejects a fast-double CDF if it disagrees wildly with a cheap multiprecision probe:

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))
)
#> $pvalues
#>    x           p      method
#> 1 10 0.998324768 fast_double
#> 2 20 0.597927350 fast_double
#> 3 30 0.076740956 fast_double
#> 4 40 0.003756455 fast_double
#> 
#> $qalpha
#>   alpha   x_crit p_back
#> 1  0.05 31.57641   0.05
#> 2  0.01 36.97514   0.01

The benchmark harness

The repository ships benchmark_family.R, a small sourceable harness that exercises this whole family and tags each call with its timing and the method used. It complements test_doubleWishartHD_sweep.R (which focuses on scaling the tail function to large p): this one answers “which function, at what cost?”.

source(system.file("validation", "benchmark_family.R", package = "rootWishartHD"))

res <- family_benchmark(which = c("small", "mid"))  # s = 10 and s = 20
family_report(res)                                  # compact usability tables

# individual probes
bench_double_backends(dsb(40, 26, 20))   # backend speed + agreement
bench_saturation(dsb(80, 54, 40))        # double-CDF saturation vs logSF (cached)
bench_single(5, 10)                      # single-Wishart family

family_report() prints four blocks: (1) CDF-backend speed and agreement, (2) the saturation table, (3) the p-value method and the qalpha round-trip, and (4) the single-Wishart family. Exact evaluations are cached under bench_cache/, so re-runs are fast.

When to call what

A rule of thumb on cost: the double backends are effectively free but saturate; the exact log path is accurate everywhere but grows quickly with s (tens of seconds per point near s = 40). The applied wrappers hide the switch, but the cost is inherited – at large s they run at exact-path speed.

References

• Chiani, M. (2014). Distribution of the largest root of a matrix for Roy’s test in multivariate analysis of variance. Journal of Multivariate Analysis, 131, 69-81.
• Chiani, M. (2016). Distribution of the largest eigenvalue for real Wishart and Gaussian random matrices and a simple approximation for the Tracy-Widom distribution. Journal of Multivariate Analysis, 143, 480-493.