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

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
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[@]dogado.de.

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.
Health stats visible at Monitor.