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

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

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.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.