# Exact and log-scale tail probabilities for Roy’s largest root with rootWishartHD

#### Stepan Grinek

## Overview

`rootWishartHD` computes distribution functions and log-scale tail probabilities for Roy’s largest root in single- and double-Wishart (Jacobi ensemble) settings. It is derived from the `rootWishart` and extends it with:

- numerically robust log-CDF and log-survival evaluation that does not saturate to 0 or 1 in extreme tails;
- tail-aware arbitrary precision, so probabilities remain accurate far into the upper tail;
- scaling that pushes reliable evaluation from `p ~ 30-40` (the practical limit of the original package) up to `p ~ 500`.

The methods follow Chiani (2014, 2016): the CDF of the largest root is a Pfaffian of a skew-symmetric matrix built from incomplete-beta entries.

    library(rootWishartHD)

The double-Wishart distribution is parametrized by `(s, m, n)` (Chiani’s notation). For a MANOVA-type problem with dimension `p`, hypothesis degrees of freedom `q_df` and error degrees of freedom related to `m_df`, the mapping used throughout this vignette is:

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

## Basic usage: the CDF

`doubleWishart()` returns the CDF `F(theta)` of the largest root statistic `theta = lambda / (1 + lambda)` on `(0, 1)`. `type = "double"` uses fast double precision; `type = "arbitrary"` uses the arbitrary-precision backend.

    par0 <- dsb_params(p = 20, m_df = 14, q_df = 10)   # s = 10
    theta <- c(0.3, 0.5, 0.7, 0.9)
    doubleWishart(theta, s = par0$s, m = par0$m, n = par0$n,
                  type = "double", verbose = FALSE)
    #> [1] 2.220446e-16 1.208198e-13 1.251410e-05 2.832402e-01

## Multiprecision backend

The portable default build uses Boost’s header-only `cpp_dec_float` backend from `BH` (`DW_USE_MPFR=0`). This is the CRAN-safe setting and does not require system MPFR/GMP libraries. To force the arbitrary-precision path at runtime, use

    options(rootWishartHD.force_multiprecision = TRUE)

Local source builds can opt in to MPFR/GMP with

    DW_USE_MPFR=1 R CMD INSTALL rootWishartHD_0.95.1.tar.gz

Runtime adaptive precision (`adaptive = TRUE`) requires an MPFR/GMP build. With the default `DW_USE_MPFR=0` build, adaptive requests are downgraded to fixed `cpp_dec_float` precision with a warning; increase fixed precision at install time, for example with `DW_MP_DIGITS=300 R CMD INSTALL .`.

Check the compiled backend with:

    rootWishartHD_mpfr_enabled()
    #> [1] FALSE

If the deprecated `force_mpfr` interface is used but the package was built with `DW_USE_MPFR=0`, `rootWishartHD` warns once and falls back to Boost `cpp_dec_float`.

## Log scale to avoid tail saturation

In moderate to high dimensions the CDF is extremely close to 0 or 1 over most of its support, so a plain double-precision CDF saturates: `1 - F` underflows to exactly 0 and the upper tail is lost. The log-survival function keeps the tail resolvable.

    # Upper-tail log-survival log(1 - F). 'lower' gives log F.
    logSF <- doubleWishart_log(
      c(0.85, 0.92, 0.97), s = par0$s, m = par0$m, n = par0$n,
      type = "arbitrary", tail = "upper", verbose = FALSE)
    data.frame(theta = c(0.85, 0.92, 0.97),
               logSF = logSF,
               SF    = exp(logSF))
    #>   theta       logSF         SF
    #> 1  0.85 -0.05257891 0.94877945
    #> 2  0.92 -0.63044036 0.53235732
    #> 3  0.97 -2.87380177 0.05648378

The plain CDF would report `1` (survival `0`) for these points, whereas the log-scale survival values remain finite and accurate.

## Regulating precision

The exact path exposes several precision knobs. The most important:

argument | meaning
--- | ---
`type` | `"double"` (fast) or `"arbitrary"` (multiprecision)
`adaptive` | grow precision at runtime until the result converges
`start_digits10` | starting decimal precision for the adaptive search
`max_digits10` | cap on decimal precision
`tol` | convergence tolerance in log-space
`pf_method` | Pfaffian backend: `"gauss"`, `"lu"`, `"svd"`, `"schur"`, `"auto"`
`scale_iter` | symmetric equilibration iterations before the Pfaffian

A practical two-stage strategy is to evaluate every point at moderate precision first, and only re-evaluate the points that underflow (return `-Inf`) at a much higher precision:

    logsf_two_stage <- function(theta, s, m, n,
                                stage1_max = 600L, stage2_max = 20000L) {
      v <- doubleWishart_log(
        theta, s = s, m = m, n = n, type = "arbitrary", tail = "upper",
        adaptive = TRUE, start_digits10 = 200L, max_digits10 = stage1_max,
        tol = 1e-8, pf_method = "gauss", verbose = FALSE)
      unresolved <- !is.finite(v)
      if (any(unresolved)) {
        v[unresolved] <- doubleWishart_log(
          theta[unresolved], s = s, m = m, n = n, type = "arbitrary", tail = "upper",
          adaptive = TRUE, start_digits10 = 200L, max_digits10 = stage2_max,
          tol = 1e-12, pf_method = "gauss", verbose = FALSE)
      }
      v
    }

    logsf_two_stage(c(0.9, 0.97), s = par0$s, m = par0$m, n = par0$n)
    #> Warning: adaptive=TRUE requires an MPFR/GMP build (DW_USE_MPFR=1). Using fixed
    #> Boost cpp_dec_float precision instead. For more fixed digits, reinstall with
    #> environment variable DW_MP_DIGITS=<digits>; for runtime adaptive precision,
    #> also set DW_USE_MPFR=1.
    #> [1] -0.3330177 -2.8738018

## Validation: empirical vs analytic CDF

We can check the analytic CDF against a Monte Carlo reference. The helper below draws the largest-root statistic by simulation (requires the `rWishart` and `corpcor` packages).

    one_draw_theta <- function(dfA, dfB, S) {
      A <- rWishart::rSingularWishart(1L, dfA, S)[, , 1]
      B <- rWishart::rSingularWishart(1L, dfB, S)[, , 1]
      eigA   <- eigen(A, symmetric = TRUE)
      rankVr <- min(dfA, nrow(A))
      V      <- eigA$vectors[, 1:rankVr, drop = FALSE]
      vals   <- pmax(eigA$values[1:rankVr], 1e-12 * max(eigA$values[1:rankVr]))
      Xp     <- sweep(V, 2, sqrt(1 / vals), `*`)
      C      <- crossprod(Xp, B %*% Xp)
      W      <- Xp %*% corpcor::fast.svd(C)$u
      lmax   <- max(crossprod(W, B %*% W))
      lmax / (1 + lmax)
    }

    set.seed(1)
    p <- 20; m_df <- 14; q_df <- 10
    par1 <- dsb_params(p, m_df, q_df)            # s = 10
    S <- diag(1, p)
    theta_mc <- replicate(100, one_draw_theta(m_df, q_df, S))
    theta_mc <- pmin(pmax(theta_mc, 1e-12), 1 - 1e-12)

    grid  <- as.numeric(quantile(theta_mc, probs = seq(0.02, 0.98, length.out = 30)))
    F_ana <- doubleWishart(grid, s = par1$s, m = par1$m, n = par1$n,
                           type = "double", verbose = FALSE)
    F_emp <- ecdf(theta_mc)(grid)

    plot(grid, F_emp, type = "s", col = "steelblue", lwd = 2, ylim = c(0, 1),
         xlab = expression(theta), ylab = expression(F(theta)),
         main = sprintf("s = %d: empirical vs analytic CDF", par1$s))
    lines(grid, F_ana, col = "firebrick", lwd = 2, lty = 2)
    points(grid, F_ana, col = "firebrick", pch = 19, cex = 0.5)
    legend("bottomright", bty = "n", lwd = 2, lty = c(1, 2),
           col = c("steelblue", "firebrick"),
           legend = c("empirical (Monte Carlo)", "analytic (rootWishartHD)"))

Empirical (Monte Carlo) vs analytic CDF for a small double-Wishart case.

    cat(sprintf("max |F_emp - F_ana| = %.4f\n", max(abs(F_emp - F_ana))))
    #> max |F_emp - F_ana| = 0.0686

## Performance across dimensions

The table below summarises measured timings on a multi-core Linux machine (R 4.6, 16 PSOCK workers). `exact/pt` is the wall time for one adaptive upper-tail `logSF` evaluation; `MC sim` is the one-off cost of generating the Monte Carlo reference (cached on disk for reuse). The repository ships a benchmarking harness, `test_doubleWishartHD_sweep.R`, that reproduces these numbers and the figure above across a grid of settings.

Representative timings for one exact upper-tail logSF evaluation.
setting | p | s | MC sim | exact/pt
--- | --- | --- | --- | ---
p40 | 40 | 20 | 4.1 s | 2.5 s
p100 | 100 | 50 | not run | 5.8 s
p150 | 150 | 35 | 2.5 min | 4.1 s
p300 | 300 | 150 | 2.5 min | 22 s
p500 | 500 | 498 | 3.1 min | 3.7 min

A key internal optimisation makes the matrix construction `O(s)` incomplete-beta evaluations instead of `O(s^2)` (the `a_{ij}` entries depend only on `i + j`). For `s = 50` this reduced one exact evaluation from about 36 s to under 6 s with bit-identical results.

## Reproducing the benchmarks

The full numerical validation sweep is included with the package but is not run automatically, because it can be computationally expensive and may use arbitrary-precision arithmetic.

    sweep_file <- system.file(
      "validation", "test_doubleWishartHD_sweep.R",
      package = "rootWishartHD"
    )

    if (!nzchar(sweep_file)) {
      stop("Validation script not found. Reinstall rootWishartHD with inst/validation included.")
    }

    source(sweep_file)

    # performance sweep over selected settings, 8 workers
    res <- run_sweep(c("p40", "p150"), n_exact = 6, n_cores = 8)
    perf_table(res, "kable")               # markdown performance table

    # empirical-vs-analytic CDF comparison + figure
    cmp <- cdf_compare("p40", n_grid = 25)
    plot_cdf_compare(cmp)

    # regulate precision: two-stage arbitrary-precision evaluation
    ex <- make_exact(start_digits10 = 150, stage1_max = 800, stage2_max = 20000)
    run_sweep("p200", exact = ex, plot = TRUE)   # writes CDF PNGs to diag_figs/

Monte Carlo draws are cached under `diag_cache/`.

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