---
title: "Exact and log-scale tail probabilities for Roy's largest root with rootWishartHD"
author: "Stepan Grinek"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Exact and log-scale tail probabilities for Roy's largest root with rootWishartHD}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.6,
  fig.align = "center"
)
has_mc <- requireNamespace("rWishart", quietly = TRUE) &&
          requireNamespace("corpcor", quietly = TRUE)
```

## 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`](https://github.com/turgeonmaxime/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.

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

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

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

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

```{r multiprecision-option, eval=FALSE}
options(rootWishartHD.force_multiprecision = TRUE)
```

Local source builds can opt in to MPFR/GMP with

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

```{r mpfr-enabled}
rootWishartHD_mpfr_enabled()
```

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.

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

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:

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

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

```{r mc-helper, eval = has_mc}
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)
}
```

```{r mc-figure, eval = has_mc, fig.cap = "Empirical (Monte Carlo) vs analytic CDF for a small double-Wishart case."}
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)"))

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

```{r mc-skip, eval = !has_mc, echo = FALSE, results = "asis"}
cat("> _The Monte Carlo validation figure is skipped because `rWishart` and/or",
    "`corpcor` are not installed._\n")
```

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

```{r perf-table, echo = FALSE}
perf <- data.frame(
  setting   = c("p40", "p100", "p150", "p300", "p500"),
  p         = c(40, 100, 150, 300, 500),
  s         = c(20, 50, 35, 150, 498),
  `MC sim`  = c("4.1 s", "not run", "2.5 min", "2.5 min", "3.1 min"),
  `exact/pt`= c("2.5 s", "5.8 s", "4.1 s", "22 s", "3.7 min"),
  check.names = FALSE, stringsAsFactors = FALSE
)
knitr::kable(perf, caption = "Representative timings for one exact upper-tail logSF evaluation.")
```

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.

```{r validation-sweep, eval=FALSE}
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. \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}