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

#### Stepan Grinek

## 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

The main numerical functions exported by rootWishartHD.
function | ensemble | returns | role
--- | --- | --- | ---
doubleWishart() | double | CDF | default CDF interface; routes to scaled backend when needed
doubleWishart_scaled() | double | CDF | scaled Pfaffian CDF backend for ill-conditioned or large cases
doubleWishart_log() | double | log CDF / log SF | tail-safe log survival and log CDF
doubleWishart_robustPfaffians() | double | log CDF / log SF | auto Pfaffian backend, stress testing
doubleWishart_pvalue() | double | p-value | applied fast-then-exact upper-tail p-value
doubleWishart_qalpha() | double | critical value | invert the tail for a cutoff
singleWishart_cdf() | single | CDF | CDF, single-Wishart largest eigenvalue
singleWishart_log() | single | log CDF / log SF | tail-safe log survival, single-Wishart
singleWishart_pvalue() | single | p-value | applied p-value, single-Wishart
singleWishart_qalpha() | single | critical value | applied critical value, single-Wishart

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.

### `doubleWishart_pvalue()`

Given an observed statistic (as `theta`, or as `lambda` with `input="lambda"`), it returns the upper-tail p-value `P(Theta >= theta)`. It first tries the fast double CDF; if that has saturated it falls back to the exact log-survival path. The `method` attribute records which path produced each value.

    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

At small `s` the fast double path resolves everything (`method = "fast_double"`). The important caveat is dimensional: once `s` is large enough that the double CDF saturates even at moderate `theta`, every point falls back to `method = "exact"` and the call pays the full multiprecision cost. At `s = 40`, for instance, the fast path is already saturated at `theta = 0.90`, so a p-value there costs on the order of ten seconds, not milliseconds.

### `doubleWishart_qalpha()`

The inverse problem: find the critical value `theta_alpha` with `P(Theta >= theta_alpha) = alpha`. It brackets the root on a fast SF grid and then bisects, using the same fast-then-exact evaluator. Round-tripping the result back through `doubleWishart_pvalue()` recovers `alpha`:

    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 quick decision guide for the function family.
goal | call
--- | ---
CDF in the body of the distribution | doubleWishart()
Upper-tail probability that may be tiny | doubleWishart_log(tail = “upper”)
Observed statistic -> p-value | doubleWishart_pvalue()
Significance level -> critical value | doubleWishart_qalpha()
Suspected Pfaffian-backend failure at large s | doubleWishart_robustPfaffians()
Single-Wishart largest eigenvalue tail | singleWishart_pvalue() / singleWishart_qalpha()

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.