qfratio-package: qfratio: Moments and Distributions of Ratios of Quadratic Forms

Description

This package is for evaluating moments of ratios (and products) of quadratic forms in normal variables, specifically using recursive algorithms developed by Bao et al. (2013) and Hillier et al. (2014) (see also Smith, 1989, 1993; Hillier et al., 2009). It also provides some functions to evaluate distribution, quantile, and probability density functions of simple ratios of quadratic forms in normal variables using several algorithms. It was originally developed as a supplement to Watanabe (2023) for evaluating average evolvability measures in evolutionary quantitative genetics, but can be used for a broader class of statistics.

Arguments

Author/Maintainer

Junya Watanabe Junya.Watanabe@uab.cat

Details

The primary front-end functions of this package are qfrm() and qfmrm() for evaluating moments of ratios of quadratic forms. These pass arguments to one of the several “internal” (though exported) functions which do actual calculations, depending on the argument matrices and exponents. In addition, there are a few functions to calculate moments of products of quadratic forms (integer exponents only; qfpm).

There are many internal functions for calculating coefficients in power-series expansion of generating functions for these moments (d1_i, d2_ij, d3_ijk, dtil2_pq) using “super-short” recursions (Bao and Kan, 2013; Hillier et al. 2014). Some of these coefficients are related to the top-order zonal and invariant polynomials of matrix arguments.

The package also has some functions to evaluate distribution, quantile, and density functions of simple ratios of quadratic forms: pqfr(), qqfr(), and dqfr().

See package vignettes (vignette("qfratio") and vignette("qfratio_distr")) for more details.

The DESCRIPTION file: tools:::Rd_package_DESCRIPTION("qfratio") tools:::Rd_package_indices("qfratio")

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

Hillier, G., Kan, R. and Wang, X. (2009) Computationally efficient recursions for top-order invariant polynomials with applications. Econometric Theory, 25, 211--242. tools:::Rd_expr_doi("10.1017/S0266466608090075").

Hillier, G., Kan, R. and Wang, X. (2014) Generating functions and short recursions, with applications to the moments of quadratic forms in noncentral normal vectors. Econometric Theory, 30, 436--473. tools:::Rd_expr_doi("10.1017/S0266466613000364").

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

Smith, M. D. (1993) Expectations of ratios of quadratic forms in normal variables: evaluating some top-order invariant polynomials. Australian Journal of Statistics, 35, 271--282. tools:::Rd_expr_doi("10.1111/j.1467-842X.1993.tb01335.x").

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

Examples

Run this code

## Symmetric matrices
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 x)^2 where x ~ N(0, I)
qfrm(A, p = 2)
## And a Monte Carlo mean of the same
mean(rqfr(1000, A = A, p = 2))

## Expectation of (x^T A x)^1/2 / (x^T x)^1/2 where x ~ N(0, I)
(res1 <- qfrm(A, p = 1/2))
plot(res1)
## A Monte Carlo mean
mean(rqfr(1000, A = A, p = 1/2))

## (x^T A x)^2 / (x^T B x)^3 where x ~ N(mu, Sigma)
(res2 <- qfrm(A, B, p = 2, q = 3, mu = mu, Sigma = Sigma))
plot(res2)
## A Monte Carlo mean
mean(rqfr(1000, A = A, B = B, p = 2, q = 3, mu = mu, Sigma = Sigma))

## Expectation of (x^T A x)^2 / (x^T B x) (x^T x) where x ~ N(0, I)
(res3 <- qfmrm(A, B, p = 2, q = 1, r = 1))
plot(res3)
## A Monte Carlo mean
mean(rqfmr(1000, A = A, B = B, p = 2, q = 1, r = 1))

## Expectation of (x^T A x)^2 where x ~ N(0, I)
qfm_Ap_int(A, 2)
## A Monte Carlo mean
mean(rqfp(1000, A = A, p = 2, q = 0, r = 0))

## Expectation of (x^T A x) (x^T B x) (x^T D x) where x ~ N(mu, I)
qfpm_ABDpqr_int(A, B, D, 1, 1, 1, mu = mu)
## A Monte Carlo mean
mean(rqfp(1000, A = A, B = B, D = D, p = 1, q = 1, r = 1, mu = mu))

## Distribution and quantile functions,
## and density of (x^T A x) / (x^T B x)
quantiles <- 0:nv + 0.5
(probs <- pqfr(quantiles, A, B))
qqfr(probs, A, B)     # p = 1 yields maximum of ratio
dqfr(quantiles, A, B)