---
title: "The rootWishartHD function family: backends, tails, p-values and critical values"
author: "Stepan Grinek"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{The rootWishartHD function family: backends, tails, p-values and critical values}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

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

```{r load}
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:

```{r 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)
}
```

## A map of the family

```{r 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.")
```

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

```{r 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))
```

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.

```{r 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))
```

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.

```{r 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"))
```

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`:

```{r 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-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.

```{r 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)
```

`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:

```{r 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))
)
```

## 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?".

```r
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

```{r 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.")
```

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. \doi{10.1016/j.jmva.2014.04.002}
* 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.
  \doi{10.1016/j.jmva.2015.10.007}