qfratio (version 1.1.1)

qfmrm: Moment of multiple ratio of quadratic forms in normal variables

Description

qfmrm() is a front-end function to obtain the (compound) moment of a multiple ratio of quadratic forms in normal variables in the following special form: \( \mathrm{E} \left( \frac{(\mathbf{x^\mathit{T} A x})^p } {(\mathbf{x^\mathit{T} B x})^q (\mathbf{x^\mathit{T} D x})^r} \right) \), where \(\mathbf{x} \sim N_n(\bm{\mu}, \mathbf{\Sigma})\). Like qfrm(), this function calls one of the following “internal” functions for actual calculation, as appropriate.

qfmrm_ApBIqr_int(): For \(\mathbf{D} = \mathbf{I}_n\) and positive-integral \(p\)

qfmrm_ApBIqr_npi(): For \(\mathbf{D} = \mathbf{I}_n\) and non-integral \(p\)

qfmrm_IpBDqr_gen(): For \(\mathbf{A} = \mathbf{I}_n\)

qfmrm_ApBDqr_int(): For general \(\mathbf{A}\), \(\mathbf{B}\), and \(\mathbf{D}\), and positive-integral \(p\)

qfmrm_ApBDqr_npi(): For general \(\mathbf{A}\), \(\mathbf{B}\), and \(\mathbf{D}\), and non-integral \(p\)

Usage

qfmrm(
  A,
  B,
  D,
  p = 1,
  q = p/2,
  r = q,
  m = 100L,
  mu = rep.int(0, n),
  Sigma = diag(n),
  tol_zero = .Machine$double.eps * 100,
  tol_sing = tol_zero,
  ...
)
qfmrm_ApBIqr_int(
  A,
  B,
  p = 1,
  q = 1,
  r = 1,
  m = 100L,
  mu = rep.int(0, n),
  error_bound = TRUE,
  check_convergence = c("relative", "strict_relative", "absolute", "none"),
  use_cpp = TRUE,
  cpp_method = c("double", "long_double", "coef_wise"),
  nthreads = 0,
  alphaB = 1,
  tol_conv = .Machine$double.eps^(1/4),
  tol_zero = .Machine$double.eps * 100,
  tol_sing = tol_zero,
  thr_margin = 100
)
qfmrm_ApBIqr_npi(
  A,
  B,
  p = 1,
  q = 1,
  r = 1,
  m = 100L,
  mu = rep.int(0, n),
  check_convergence = c("relative", "strict_relative", "absolute", "none"),
  use_cpp = TRUE,
  cpp_method = c("double", "long_double", "coef_wise"),
  nthreads = 0,
  alphaA = 1,
  alphaB = 1,
  tol_conv = .Machine$double.eps^(1/4),
  tol_zero = .Machine$double.eps * 100,
  tol_sing = tol_zero,
  thr_margin = 100
)
qfmrm_IpBDqr_gen(
  B,
  D,
  p = 1,
  q = 1,
  r = 1,
  mu = rep.int(0, n),
  m = 100L,
  check_convergence = c("relative", "strict_relative", "absolute", "none"),
  use_cpp = TRUE,
  cpp_method = c("double", "long_double", "coef_wise"),
  nthreads = 0,
  alphaB = 1,
  alphaD = 1,
  tol_conv = .Machine$double.eps^(1/4),
  tol_zero = .Machine$double.eps * 100,
  tol_sing = tol_zero,
  thr_margin = 100
)
qfmrm_ApBDqr_int(
  A,
  B,
  D,
  p = 1,
  q = 1,
  r = 1,
  m = 100L,
  mu = rep.int(0, n),
  check_convergence = c("relative", "strict_relative", "absolute", "none"),
  use_cpp = TRUE,
  cpp_method = c("double", "long_double", "coef_wise"),
  nthreads = 0,
  alphaB = 1,
  alphaD = 1,
  tol_conv = .Machine$double.eps^(1/4),
  tol_zero = .Machine$double.eps * 100,
  tol_sing = tol_zero,
  thr_margin = 100
)
qfmrm_ApBDqr_npi(
  A,
  B,
  D,
  p = 1,
  q = 1,
  r = 1,
  m = 100L,
  mu = rep.int(0, n),
  check_convergence = c("relative", "strict_relative", "absolute", "none"),
  use_cpp = TRUE,
  cpp_method = c("double", "long_double", "coef_wise"),
  nthreads = 0,
  alphaA = 1,
  alphaB = 1,
  alphaD = 1,
  tol_conv = .Machine$double.eps^(1/4),
  tol_zero = .Machine$double.eps * 100,
  tol_sing = tol_zero,
  thr_margin = 100
)

Value

A qfrm object, as in qfrm() functions.

Arguments

A, B, D

Argument matrices. Should be square. Will be automatically symmetrized.

p, q, r

Exponents for \(\mathbf{A}\), \(\mathbf{B}\), and \(\mathbf{D}\), respectively. By default, q equals p/2 and r equals q. If unsure, specify all explicitly.

m

Order of polynomials at which the series expression is truncated. \(M\) in Hillier et al. (2009, 2014).

mu

Mean vector \(\bm{\mu}\) for \(\mathbf{x}\)

Sigma

Covariance matrix \(\mathbf{\Sigma}\) for \(\mathbf{x}\). Accommodated only by the front-end qfmrm(). See “Details” in qfrm.

tol_zero

Tolerance against which numerical zero is determined. Used to determine, e.g., whether mu is a zero vector, A or B equals the identity matrix, etc.

tol_sing

Tolerance against which matrix singularity and rank are determined. The eigenvalues smaller than this are considered zero.

...

Additional arguments in the front-end qfmrm() will be passed to the appropriate “internal” function.

error_bound

Logical to specify whether an error bound is returned (if available).

check_convergence

Specifies how numerical convergence is checked (see “Details”). Options:

"relative"

default; magnitude of the last term of the series relative to the sum is compared with tol_conv

"strict_relative" or TRUE

same, but stricter than default by setting tol_conv = .Machine$double.eps (unless a smaller value is specified by the user)

"absolute"

absolute magnitude of the last term is compared with tol_conv

"none" or FALSE

skips convergence check

use_cpp

Logical to specify whether the calculation is done with C++ functions via Rcpp. TRUE by default.

cpp_method

Method used in C++ calculations to avoid numerical overflow/underflow (see “Details”). Options:

"double"

default; fastest but prone to underflow in some conditions

"long_double"

same algorithm but using the long double variable type; robust but slow and memory-inefficient

"coef_wise"

coefficient-wise scaling algorithm; most robust but variably slow

nthreads

Number of threads used in OpenMP-enabled C++ functions. See “Multithreading” in qfrm.

tol_conv

Tolerance against which numerical convergence of series is checked. Used with check_convergence.

thr_margin

Optional argument to adjust the threshold for scaling (see “Scaling” in d1_i). Passed to internal functions (d1_i, d2_ij, d3_ijk) or their C++ equivalents.

alphaA, alphaB, alphaD

Factors for the scaling constants for \(\mathbf{A}\), \(\mathbf{B}\), and \(\mathbf{D}\), respectively. See “Details” in qfrm.

Details

The usage of these functions is similar to qfrm, to which the user is referred for documentation. It is assumed that \(\mathbf{B} \neq \mathbf{D}\) (otherwise, the problem reduces to a simple ratio).

When B is identity or missing, this and its exponent q will be swapped with D and r, respectively, before qfmrm_ApBIqr_***() is called.

The existence conditions for the moments of this multiple ratio can be reduced to those for a simple ratio, provided that one of the null spaces of \(\mathbf{B}\) and \(\mathbf{D}\) is a subspace of the other (including the case they are null). The conditions of Bao and Kan (2013: proposition 1) can then be applied by replacing \(q\) and \(m\) there by \(q + r\) and \(\min{( \mathrm{rank}(\mathbf{B}), \mathrm{rank}(\mathbf{D}) )} \), respectively (see also Smith 1989: p. 258 for nonsingular \(\mathbf{B}\), \(\mathbf{D}\)). An error is thrown if these conditions are not met in this case. Otherwise (i.e., \(\mathbf{B}\) and \(\mathbf{D}\) are both singular and neither of their null spaces is a subspace of the other), it seems difficult to define general moment existence conditions. A sufficient condition can be obtained by applying the same proposition with a new denominator matrix whose null space is union of those of \(\mathbf{B}\) and \(\mathbf{D}\) (Watanabe, 2023). A warning is thrown if that condition is not met in this case.

Most of these functions, excepting qfmrm_ApBIqr_int() with zero mu, involve evaluation of multiple series, which can suffer from numerical overflow and underflow (see “Scaling” in d1_i and “Details” in qfrm). To avoid this, cpp_method = "long_double" or "coef_wise" options can be used (see “Details” in qfrm).

References

Bao, Y. and Kan, R. (2013) On the moments of ratios of quadratic forms in normal random variables. Journal of Multivariate Analysis, 117, 229--245. tools:::Rd_expr_doi("10.1016/j.jmva.2013.03.002").

Smith, M. D. (1989). On the expectation of a ratio of quadratic forms in normal variables. Journal of Multivariate Analysis, 31, 244--257. tools:::Rd_expr_doi("10.1016/0047-259X(89)90065-1").

Watanabe, J. (2023) Exact expressions and numerical evaluation of average evolvability measures for characterizing and comparing G matrices. Journal of Mathematical Biology, 86, 95. tools:::Rd_expr_doi("10.1007/s00285-023-01930-8").

See Also

qfrm for simple ratio

Examples

## Some symmetric matrices and parameters
nv <- 4
A <- diag(nv:1)
B <- diag(sqrt(1:nv))
D <- diag((1:nv)^2 / nv)
mu <- nv:1 / nv
Sigma <- matrix(0.5, nv, nv)
diag(Sigma) <- 1

## Expectation of (x^T A x)^2 / (x^T B x) (x^T x) where x ~ N(0, I)
(res1 <- qfmrm(A, B, p = 2, q = 1, r = 1))
plot(res1)

# The above internally calls the following:
qfmrm_ApBIqr_int(A, B, p = 2, q = 1, r = 1) ## The same

# Similar result with different expression
# This is a suboptimal option and throws a warning
qfmrm_ApBIqr_npi(A, B, p = 2, q = 1, r = 1)

## Expectation of (x^T A x) / (x^T B x)^(1/2) (x^T D x)^(1/2) where x ~ N(0, I)
(res2 <- qfmrm(A, B, D, p = 1, q = 1/2, r = 1/2))
plot(res2)

# The above internally calls the following:
qfmrm_ApBDqr_int(A, B, D, p = 1, q = 1/2, r = 1/2) ## The same

## Average response correlation between A and B
(res3 <- qfmrm(crossprod(A, B), crossprod(A), crossprod(B),
               p = 1, q = 1/2, r = 1/2))
plot(res3)

## Same, but with x ~ N(mu, Sigma)
(res4 <- qfmrm(crossprod(A, B), crossprod(A), crossprod(B),
               p = 1, q = 1/2, r = 1/2, mu = mu, Sigma = Sigma))
plot(res4)

## Average autonomy of D
(res5 <- qfmrm(B = D, D = solve(D), p = 2, q = 1, r = 1))
plot(res5)