Distribution function

Forchini (2002, 2005) derived explicit expressions of F_Q using the top-order zonal and invariant polynomials.

Let Λ from above be arranged and partitioned such that where Λ₁(q) and −Λ₂(q) are n₊- and n₋-dimensional diagonal matrices of positive and negative eigenvalues of A − qB, respectively. By denoting with the partition of ν corresponding to that of the rows of Λ above, the expression of Forchini (2005) is, after correcting some errors, where $\qfrCid{\qfrbrc{k_1}}{\qfrbrc{k_2}}{\qfrbrc{k_1 + k_2}}{ \cdot }{ \cdot }$ are the top-order invariant polynomials of (k₁, k₂)-th degree (see above for other notations).

In the central case (μ = 0_n), the above expression simplifies into (Forchini, 2002) where $\qfrC[\qfrbrc{k}]{ \cdot }$ are the top-order zonal polynomials of kth degree.

These expressions can be numerically evaluated by truncating the infinite series at certain higher-order terms (i₁ + j₁ + i₂ + j₂ = m, say), and by using a recursive algorithm to calculate by Hillier et al. (2009, 2014) (see also the main vignette qfratio). The present package implements this algorithm in pqfr(..., method = "forchini") (see below).

The distribution function has points of nonanalyticity around the eigenvalues of B⁻¹A (assuming B⁻¹ is invertible) (Forchini, 2002). Practically speaking, around these points, the series expression can be slow to converge and evaluation of the hypergeometric function can fail because the argument $\qfrtr \left( {\boldsymbol{\Lambda}_1^*}^{-1} \right)$ becomes very close to 1. Otherwise, the series expression can be evaluated with high accuracy, although the computational cost of the recursive calculations can be substantial in large problems.

Density function

Apparently, the literature has an explicit expression for the density of Q only for the simple condition when B = I_n and μ = 0_n. In this case, the distribution of Q does not depend on the norm of x, so any spherically symmetric distribution of x yields the same distribution of Q (Hillier, 2001).

Let η₁ > … > η_s be s distinct eigenvalues of A, and n₁, …, n_s be the corresponding degrees of multiplicity ($\sum_{i=1}^{s} n_i = n$). Consider the transformed variable $V = \frac{Q - \eta_s}{\eta_1 - Q}$ and parameters $\psi_i = \frac{\eta_i - \eta_s}{\eta_1 - \eta_i}$ for i = 2, …s − 1, assuming s > 2 (see below for the case when s = 2). The density of V has different functional forms across its domain (from Hillier, 2001, lemmas 3 and 4): where $p_r = \sum_{i=1}^{r+1} n_i$, $\mathbf{D}_{r} = \qfrdiag \left( \psi_{2}^{-1} \mathbf{I}_{n_{2}} , \dots , \psi_{r+1}^{-1} \mathbf{I}_{n_{r+1}} \right)$, $\mathbf{D}_{r+1} = \qfrdiag \left( \psi_{r+2}^{-1} \mathbf{I}_{n_{r+2}} , \dots , \psi_{s-1}^{-1} \mathbf{I}_{n_{s-1}} \right)$, $\mathbf{D} = \qfrdiag \left( \psi_{2}^{-1} \mathbf{I}_{n_{2}} , \dots , \psi_{s-1}^{-1} \mathbf{I}_{n_{s-1}} \right)$, and c_r(j, k) are the coefficients defined as c_r(j, k) can be 0 when p_r or n − p_r is even, so that some terms in the above series disappear. Otherwise (whenever c_r(j, k) ≠ 0), it is possible to write to simplify calculation. The density is undefined at v = ψ_i (q = η_i) for any i.

From one of the above expressions, f_Q(q) can be obtained as It is seen that, when s = 2 (i.e., there are only two distinct eigenvalues), the above becomes which is the density of the (scaled) beta distribution in the interval [η_s, η₁] with the parameters n₁/2 and n_s/2. This result is expected from the basic relationship between the chi-square and beta distributions, i.e., $\qfrcchisq{n_1} / \left( \qfrcchisq{n_1} + \qfrcchisq{n_s} \right) = b \sim \qfrBtd{n_1 / 2}{n_s / 2}$, a standard beta distribution in [0, 1], with $\qfrcchisq{n_i}$ being independent chi-square variables; $Q = \left( \eta_1 \qfrcchisq{n_1} + \eta_s \qfrcchisq{n_s} \right) / \left( \qfrcchisq{n_1} + \qfrcchisq{n_s} \right) = \eta_1 b + \eta_s \left( 1 - b \right) = \left( \eta_1 - \eta_s \right) b + \eta_s$.

As for the above series expression of F_Q, these expressions can be evaluated by taking a partial sum of the series and using the recursive algorithm for d. dqfr(..., method = "hillier") implements this algorithm (see below). This is reasonably quick and accurate in small problems, but the computational cost can be substantial in large problems.

Distribution function

From the famous formula of Imhof (1961) on the distribution of X_q, where

The above integral can be evaluated by using a regular numerical evaluation algorithm for infinite intervals. Alternatively, it can be evaluated by the trapezoidal integration algorithm of Davies (1973, 1980) which explicitly controls the numerical errors involved. The package CompQuadForm implements these methods in the function imhof() and davies(), respectively (strictly speaking, these are for 1 − F_Q as per Imhof’s original result). The present package utilizes those functions via pqfr(..., method = "imhof", use_cpp = FALSE) and pqfr(..., method = "davies"), respectively (see below). For the former method, the present package also has its own C++ implementation used via pqfr(..., method = "imhof", use_cpp = TRUE) (default). The numerical inversion can be evaluated fairly quickly on modern computers, and the accuracy will be sufficient for most practical purposes with slight error in numerical integration.

Density function

The density can be evaluated by numerical inversion of the characteristic function using Geary’s formula (Geary, 1944; Stuart and Ord, 1994, sec. 11.10). Broda and Paolella (2009) demonstrated that where, along with β and γ defined above, with H = P^TBP = (h_ij) and F = I_n + u²Λ².

The above expression can be evaluated with a regular numerical integration algorithm, and is implemented in dqfr(..., method = "broda") (see below).

Distribution function

A first-order saddlepoint approximation formula for the distribution function F_Q is (Butler and Paolella, 2007, 2008): where Φ(⋅) and ϕ(⋅) are the distribution and density functions, respectively, of the standard normal distribution, and The condition $\qfrE \left( X_q \right) = 0$ is equivalent to ŝ = 0, because of the elementary property of the cumulant generating function $\qfrE \left( X_q \right) = K_{X_q}' \left( 0 \right)$. Higher-order derivatives of K_Xq are (see also Paolella, 2007, chap. 10)

A more accurate second-order approximation is (Butler and Paolella, 2007) where κ̂_j = K_Xq^(j)(ŝ)/K_Xq″(ŝ)^j/2.

Evaluation of saddlepoint approximation is fairly quick, with the only potential complexity arising from the numerical root-finding. Empirically, the accuracy of this approximation seems to improve for large problems. This is expected since the distribution of X_q as a weighted sum approaches normality as n increases.

pqfr(..., method = "butler") implements this saddlepoint approximation. The second-order approximation is used by default (order_spa = 2) (see below).

Choosing a method

The evaluation method is specified by the argument method in these functions (by default, both functions choose a numerical inversion method). According to the choice, the actual calculations are done in one of the internal functions described below. Direct use of the internal functions are not recommended. These internal functions only accept a length-one quantile. To reduce computational time, they do not check argument structures or accommodate Sigma.

## Choice from alternative methods
A <- diag(1:3)
pqfr(1.5, A, method = "imhof")    # default
#> [1] 0.1978686
pqfr(1.5, A, method = "davies")   # similar
#> [1] 0.1979048
pqfr(1.5, A, method = "forchini") # series
#> [1] 0.1978686
pqfr(1.5, A, method = "butler")   # spa
#> [1] 0.1897189

dqfr(1.5, A, method = "broda")    # default
#> [1] 0.4506431
dqfr(1.5, A, method = "hillier")  # series
#> [1] 0.4506431
dqfr(1.5, A, method = "butler")   # spa
#> [1] 0.4523631

## Not recommended; for diagnostic use only
qfratio:::pqfr_imhof(1.5, A)
#> $p
#> [1] 0.1978686
#> 
#> $abserr
#> [1] 2.031334e-09
qfratio:::pqfr_A1B1(1.5, A, m = 9, check_convergence = FALSE)
#> $p
#> [1] 0.1978641
#> 
#> $terms
#>  [1] 1.495393e-01 3.427348e-02 9.516357e-03 2.974620e-03 1.002253e-03
#>  [6] 3.544878e-04 1.295434e-04 4.844575e-05 1.842967e-05 7.103869e-06


## This is okay
x <- c(1.5, 2.5, 3.5)
pqfr(x, A)
#> [1] 0.1978686 0.8021314 1.0000000

## This is not
qfratio:::pqfr_imhof(x, A)
#> Error in qfratio:::pqfr_imhof(x, A): In pqfr_imhof, quantile must be length-one

Use with ks.test()

In principle, pqfr() is compatible with stats:::ks.test(), but care must be exercised as evaluation result may involve non-trivial error. It is recommended to inspect error bounds beforehand, using pqfr(..., method = "imhof", return_abserr_attr = TRUE). In addition, the argument B is pre-occupied by the same-named argument in ks.test(), so it cannot be passed via ...; this means that a typical syntax with non-default B should be something like ks.test(x, function(q) pqfr(q, A, B, ...)) rather than ks.test(x, pqfr, A = A, B = B, ...)). For example,

## Small Monte Carlo sample
A <- diag(1:3)
B <- diag(sqrt(1:3))
x <- rqfr(10, A, B)

## Calculate p-values
pseq <- pqfr(x, A, B, return_abserr_attr = TRUE)

## Maximum error when evaluated at x;
## looks small enough
max(attr(pseq, "abserr"))
#> [1] 8.318564e-07

## Correct syntax, expected outcome
## \(q) syntax could also be used in recent versions of R
ks.test(x, function(q) pqfr(q, A, B))
#> 
#>  Exact one-sample Kolmogorov-Smirnov test
#> 
#> data:  x
#> D = 0.29324, p-value = 0.2946
#> alternative hypothesis: two-sided

rather than

## Incorrect; no error/warning because
## B is passed to ks.test rather than to pqfr
ks.test(x, pqfr, A = A, B = B)
#> 
#>  Exact one-sample Kolmogorov-Smirnov test
#> 
#> data:  x
#> D = 0.78335, p-value = 5.182e-07
#> alternative hypothesis: two-sided

Specifying integration error

The above integration functions try to find an absolute error bound e_I that is bounded by the user-specified tolerance for absolute ϵ_abs and relative ϵ_rel errors: e_I ≤ ϵ_abs + |I|ϵ_rel, where I is the result of integration.

Internally, ϵ_abs is calculated from the user-specified arguments epsabs and epsrel to appropriately constrain the density or distribution function (whereas ϵ_rel is always specified by epsrel). In dqfr_broda(), pi * epsabs is used as ϵ_abs, and the resultant error bound abserr is subsequently divided by pi, so is the integration result itself to yield the density: f_Q = I/π (see above).

Situation is more complicated for pqfr_imhof(), because the relative error in I cannot in general be directly transformed to that of the distribution function, which is F_Q = 1/2 − I/π (see above). In this function, pi * (epsabs * epsrel / 2) is passed as ϵ_abs, and the resultant error bound abserr is divided by pi. This procedure ensures an equivalent of the above inequality to hold for F_Q, provided I ≤ 0 (F_Q ≥ 1/2) or ϵ_rel = 0. Otherwise, an error bound calculated in the same way can only be conservative; pqfr_imhof() returns this value, but it can violate the user-specified relative tolerance epsrel.

A <- diag(1:4)

## This error bound satisfies "abserr < value * epsrel"
pqfr(3.9, A, method = "imhof", return_abserr_attr = TRUE,
     epsabs = 0, epsrel = 1e-6)
#> [1] 0.9944167
#> attr(,"abserr")
#> [1] 3.325725e-08

## This one violates "abserr < value * epsrel",
## although abserr is a valid error bound
pqfr(1.2, A, method = "imhof", return_abserr_attr = TRUE,
     epsabs = 0, epsrel = 1e-6)
#> [1] 0.01611023
#> attr(,"abserr")
#> [1] 4.165996e-07

autoscale_args

Numerical integration involved in these functions typically fail when the magnitude of eigenvalues is too small or too large, whence the integrand functions can decrease too slowly (i.e., divergent-looking) or too quickly (i.e., looks constant 0) with respect to the integration parameter (u above). To avoid such failures, these functions internally scale the eigenvalues by default, so that max λ_i − min λ_i is equal to the argument autoscale_args (default 1); remember that min λ_i is negative, so this quantity is sum of the absolute values.

A <- diag(1:3)
B <- diag(sqrt(1:3))

## Without autoscale_args
## We know these are equal
pqfr(1.5, A, B, autoscale_args = FALSE)
#> [1] 0.6376791
pqfr(1.5, A * 1e-10, B * 1e-10, autoscale_args = FALSE)
#> [1] 0.5
## The latter failed because of numerically small eigenvalues

## With autoscale_args = 1 (default)
pqfr(1.5, A * 1e-10, B * 1e-10)
#> [1] 0.6376791

trim_values

Numerical integration can yield spurious results that are outside the mathematically permissible supports; [0, ∞) and [0, 1] for the density and distribution functions, respectively. By default (trim_values = TRUE), the external functions dqfr() and pqfr() trim those values into the permissible range by using tol_zero as a margin; e.g., negative p-values are replaced by ~2.2e-14 (default tol_zero). A warning is thrown if this happens, because it usually means that numerically accurate evaluation was impossible, at least with the given parameters. Turn trim_values = FALSE to skip these trimming and warning, although pqfr_imhof() and pqfr_davies() can still throw a warning from CompQuadForm functions. Note that, on the other hand, all these functions try to return exact 0 or 1 when q is outside the possible range of Q (as numerically determined).

## Result without trimming;
## (typically) negative density, which is absurd
## In this case, error interval typically spans across 0
dqfr(1.2, diag(1:30), return_abserr_attr = TRUE,
     trim_values = FALSE)
#> [1] -4.859181e-17
#> attr(,"abserr")
#> [1] 8.043729e-08

## Result with trimming (default)
dqfr(1.2, diag(1:30), return_abserr_attr = TRUE)
#> Warning in dqfr(1.2, diag(1:30), return_abserr_attr = TRUE): values < 0 trimmed
#> up to tol_zero
#> [1] 2.220446e-14
#> attr(,"abserr")
#> [1] 8.043726e-08
## Note that the actual value is only bounded by
## 0 and abserr

Options

The saddlepoint approximation density does not integrate to unity, but can be normalized by setting normalize_spa = TRUE in dqfr() (note that this is done in the external function). The normalized density can be more accurate (although it is usually a matter of empiricism). However, this is usually slower than the numerical inversion method for a small number of quantiles.

The second-order approximation is used by default (order_spa = 2) (internally, any value > 1 calls this option). The first-order approximation can be used by setting order_spa = 1, but this is usually less accurate and only slightly faster than the second-order approximation.

A <- diag(1:3)

## Default for spa distribution function
pqfr(1.2, A, method = "butler", order_spa = 2)
#> [1] 0.07183068

## First-order spa
pqfr(1.2, A, method = "butler", order_spa = 1)
#> [1] 0.0790331

## More accurate numerical inversion
pqfr(1.2, A)
#> [1] 0.07359703


## Default for density
dqfr(1.2, A, method = "butler",
     order_spa = 2, normalize_spa = FALSE)
#> [1] 0.3716931

## First-order
dqfr(1.2, A, method = "butler",
     order_spa = 1, normalize_spa = FALSE)
#> [1] 0.4577787

## Normalized density, second-order
dqfr(1.2, A, method = "butler",
     order_spa = 2, normalize_spa = TRUE)
#> [1] 0.3913688

## Normalized density, first-order
dqfr(1.2, A, method = "butler",
     order_spa = 1, normalize_spa = TRUE)
#> [1] 0.4412349

## More accurate numerical inversion
dqfr(1.2, A)
#> [1] 0.3837318

Paolella (2007, program listing 10.4) noted that the second-order approximation for the distribution function can be “problematic”, which presumably means that the evaluation result can be unstable. In development of this package, some instability in the second-order approximation was encountered, but experiments suggest that this was due to sensitivity of the result to the numerically found root ŝ. This instability is rarely encountered with the present default setting, but the user may want to adjust root-finding-related parameters when any doubt exists.

pqfr() and dqfr()

Return values from pqfr_imhof() and dqfr_broda() have an error bound abserr for numerical integration, along with the evaluation result itself. Technically, the error bound from the integration algorithm is divided by pi before returned, as the evaluation result itself is. This can be passed to the external functions and returned as an attribute by setting return_abserr_attr = TRUE (as already used in above examples):

A <- diag(1:4)

pqfr(1.5, A, return_abserr_attr = TRUE)
#> [1] 0.06819534
#> attr(,"abserr")
#> [1] 9.576418e-08

dqfr(1.5, A, return_abserr_attr = TRUE)
#> [1] 0.22202
#> attr(,"abserr")
#> [1] 5.738788e-08

This error bound tries to accommodate the effect of trim_values. If the integration result is outside the permissible support (e.g., negative density), the possible error bound is only on the direction toward the support (assuming things are calculated accurately). The returned abserr is truncated accordingly, unless trimming is beyond the original abserr (in which case it is replaced by tol_zero). See this in examples:

## Without trimming, result is (typically) negative
## But note that value + abserr is positive
dqfr(1.2, diag(1:35), return_abserr_attr = TRUE,
     epsabs = 1e-10, trim_values = FALSE)
#> [1] -6.626156e-18
#> attr(,"abserr")
#> [1] 2.378925e-12

## With trimming, value is replaced by tol_zero
## Note slightly shortened abserr
dqfr(1.2, diag(1:35), return_abserr_attr = TRUE,
     epsabs = 1e-10)
#> Warning in dqfr(1.2, diag(1:35), return_abserr_attr = TRUE, epsabs = 1e-10):
#> values < 0 trimmed up to tol_zero
#> [1] 2.220446e-14
#> attr(,"abserr")
#> [1] 2.356714e-12


## When untrimmed value + abserr < tol_zero
dqfr(1.1, diag(1:35), return_abserr_attr = TRUE,
     epsabs = 1e-15, trim_values = FALSE)
#> [1] -3.313078e-18
#> attr(,"abserr")
#> [1] 7.760891e-16
## True value is somewhere between 0 and value + abserr
## (assuming these are reliable)

## When trimmed, abserr reflects tol_zero
## because the true value is between 0 and tol_zero
dqfr(1.1, diag(1:35), return_abserr_attr = TRUE,
     epsabs = 1e-15)
#> Warning in dqfr(1.1, diag(1:35), return_abserr_attr = TRUE, epsabs = 1e-15):
#> values < 0 trimmed up to tol_zero
#> [1] 2.220446e-14
#> attr(,"abserr")
#> [1] 2.220446e-14

When log/log.p = TRUE, abserr is transformed so that it is a conservative absolute error bound on the log scale. That is, if the original value and its error bound is denoted by x̂ and δx̂, respectively, and the log-transformed value and its error bound is by log x̂ and δ(log x̂), the latter error bound is set so that log x̂ − δ(log x̂) = log (x̂ − δx̂), i.e., $\delta (\log \hat{x}) = - \log \left( 1 - \frac{\delta \hat{x}}{\hat{x}} \right)$. Note that the upper error bound $\log \left( 1 + \frac{\delta \hat{x}}{\hat{x}} \right)$ is narrower than this unless δx̂ > x̂ (i.e., x̂ − δx̂ < 0), in which case it should be taken as δ(log x̂) = ∞. In summary, the new error bound is calculated as ifelse(abserr > ans, Inf, -log1p(-abserr/ans)).

qqfr()

The option return_abserr_attr = TRUE is available in qqfr() as well:

A <- diag(1:4)

qqfr(0.95, A, return_abserr_attr = TRUE)
#> [1] 3.587557
#> attr(,"abserr")
#> [1] 6.648746e-06

In qqfr(), numerical errors arise from the root-finding with stats::uniroot() as well as in propagation from pqfr(). When return_abserr_attr = TRUE, it tries to evaluate a conservative error bound as follows:

1. Store the estimated error in root-finding uniroot()$estim.prec as δq_rf
2. Store that in pqfr() at the root as δF
3. If δF ≠ 0, calculate the density f and its error bound δf at the root using dqfr(..., method = "broda"), so that the conservative slope of the tangent there is b = max (f − δf, 0). The error δq_p in the root arising from pqfr() is at most b⁻¹δF. If δF = 0, δq_p = 0 regardless of b.
4. The total error in the quantile is δq_rf + δq_p

If log.p = TRUE, the root-finding is done on log F, so the slope used in 3 is replaced by b = max (f − δf, 0)/F.

For probability = 0 or 1, the quantile corresponds to the minimum or maximum of the ratio, and the above error bound does not apply. At present, an arbitrary value of .Machine$double.eps * 100 (~2.2e-14) is returned as an error bound for a finite minimum/maximum, although the actual error in calculation can be larger.

qqfr(0, A, return_abserr_attr = TRUE)
#> [1] 1
#> attr(,"abserr")
#> [1] 2.220446e-14

When A is nonnegative definite or p is an odd integer

In this case, the relationship between Q and R is one-to-one, so that Thus, the result can be obtained by a single recursive call of pqfr(..., p = 1) or dqfr(..., p = 1) with transformed quantile.

When A is indefinite and p is an even integer

In this case, R is an even function of Q, so that Thus, for r > 0, the result is obtained from two recursive calls of pqfr(..., p = 1) or dqfr(..., p = 1) with transformed quantile.

When A is indefinite and p is non-integer

In this case, R can be undefined, so pqfr() and dqfr() return an error, "A must be nonnegative definite when p is non-integer", regardless of the value of quantile.