| qfmrm {qfratio} | R Documentation |
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
)
Arguments
A, B, D |
Argument matrices. Should be square. Will be automatically symmetrized. |
p, q, r |
Exponents for |
m |
Order of polynomials at which the series expression is truncated. |
mu |
Mean vector |
Sigma |
Covariance matrix |
tol_zero |
Tolerance against which numerical zero is determined. Used to determine,
e.g., whether |
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 |
error_bound |
Logical to specify whether an error bound is returned (if available). |
check_convergence |
Specifies how numerical convergence is checked (see “Details”). Options:
|
use_cpp |
Logical to specify whether the calculation is done with C++
functions via |
cpp_method |
Method used in C++ calculations to avoid numerical overflow/underflow (see “Details”). Options:
|
nthreads |
Number of threads used in OpenMP-enabled C++ functions. See
“Multithreading” in |
tol_conv |
Tolerance against which numerical convergence of series is checked. Used
with |
thr_margin |
Optional argument to adjust the threshold for scaling (see “Scaling”
in |
alphaA, alphaB, alphaD |
Factors for the scaling constants for |
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).
Value
A qfrm object, as in qfrm() functions.
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. 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. 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. 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)