# qfratio: Moments of Ratios of Quadratic Forms

\[ \DeclareMathOperator{\qfrE}{E} \DeclareMathOperator{\qfrtr}{tr} \newcommand{\qfrGmf}[1]{\Gamma \! \left( #1 \right)} \newcommand{\qfrbrc}[1]{\left[ {#1} \right]} \newcommand{\qfrCit}[7]{ C^{#1, #2, #3}_{#4} \! \left( #5, #6, #7 \right) } \newcommand{\qfrrf}[2][k]{\left( #2 \right)_{#1}} \newcommand{\qfrdk}[2][k]{ d_{#1} \! \left( #2 \right) } \newcommand{\qfrdijk}[4][k]{ d_{#1} \! \left( #2, #3, #4 \right) } \newcommand{\qfrdijkl}[5][k]{ d_{#1} \! \left( #2, #3, #4, #5 \right) } \newcommand{\qfrdtk}[3][k]{ \tilde{d}_{#1} \! \left( #2, #3 \right) } \newcommand{\qfrhij}[3][k]{ h_{#1} \! \left( #2, #3 \right) } \newcommand{\qfrhijk}[4][k]{ h_{#1} \! \left( #2, #3, #4 \right) } \newcommand{\qfrhtij}[3][k]{ \tilde{h}_{#1} \! \left( #2; #3 \right) } \newcommand{\qfrhtijk}[4][k]{ \tilde{h}_{#1} \! \left( #2; #3, #4 \right) } \newcommand{\qfrhhij}[3][k]{ \hat{h}_{#1} \! \left( #2; #3 \right) } \newcommand{\qfrhhijk}[4][k]{ \hat{h}_{#1} \! \left( #2; #3, #4 \right) } \renewcommand{\det}[1]{\left\lvert #1 \right\rvert} \newcommand{\qfrlmax}{\lambda_{\max}} \newcommand{\qfrmvnorm}[3][n]{N_{#1} \! \left( #2 , #3 \right)} \]

This document is to provide an overview of the moment-related functionalities of the `R` package `qfratio` (CRAN; GitHub). Distribution-related functionalities are documented in a separate vignette: `vignette("qfratio_distr")`.

## Symbols used

- \(n\): number of variables
- \(m\): order of evaluation at which infinite series is truncated
- \(p, q, r\): exponents of quadratic forms (nonnegative real numbers)
- \(t_1, \dots, t_s\): parameters in generating functions; \(s\) here denotes the number of arguments (typically \(s = 2\) or \(3\))
- \(a\), \(d\), \(\tilde{d}\), \(h\), \(\tilde{h}\), \(\hat{h}\): coefficients in generating functions defined in text
- \(C\): (top-order) zonal and invariant polynomials of matrix arguments (Chikuse, 1980, 1987; Davis, 1980; Muirhead, 1982)
- \(\mathbf{x}\): \(n\)-variate random vector (in most cases normal)
- \(\qfrmvnorm{\boldsymbol{\mu}}{\boldsymbol{\Sigma}}\): \(n\)-variate normal distribution with mean vector \(\boldsymbol{\mu}\) and covariance matrix \(\boldsymbol{\Sigma}\)
- \(\mathbf{A}, \mathbf{B}, \mathbf{D}\): \(n \times n\) argument matrices
- \(\mathbf{I}_n\): \(n\)-dimensional identity matrix
- \(\mathbf{0}_n\): \(n\)-variate vector of \(0\)’s
- \(\mathbf{0}_{n \times n}\): \(n \times n\) matrix of \(0\)’s
- \(\qfrE \left( \cdot \right)\): expectation/mean
- \(\qfrGmf{\cdot}\): gamma function
- \((x)_a\): Pochhammer symbol, that is, \((x)_a = x (x + 1) \dots (x + a - 1)\) (with the convention \((x)_0 = 1\))
- \(\mathbf{A}^T\): matrix transposition
- \(\mathbf{A}^{-1}\): matrix inverse
- \(\mathbf{A}^{-\frac{1}{2}}\): (symmetric) matrix satisfying \(\mathbf{A}^{-\frac{1}{2}} \mathbf{A} \mathbf{A}^{-\frac{1}{2}} = \mathbf{I}_n\), assuming positive definite \(\mathbf{A}\)
- \(\det{\mathbf{A}}\): matrix determinant
- \(\qfrtr{\mathbf{A}}\): matrix trace
- \(\lambda_i\): \(i\)th eigenvalue of the matrix concerned
- \(\qfrlmax\): largest eigenvalue of the matrix concerned
- \(\mathbf{h}_i\): eigenvector corresponding to \(\lambda_i\)
- \(\beta_{\mathbf{A}}\): scaling constant that satisfies \(0 < \beta < 2 / \qfrlmax\) for \(\mathbf{A}\)
- \(\mathbf{A}^+\): matrix with “positivized” eigenvalues, i.e., \(\mathbf{A}^+ = \sum_{i=1}^n \lvert \lambda_i \rvert \mathbf{h}_i \mathbf{h}_i^T\)

The use of \(n\) for the number of variables and \(p, q, \dots\) for exponents is to conform with Bao and Kan (2013).

# Target

The package primarily concerns moments of ratios of quadratic forms of the following forms: \[\begin{equation} \frac{\left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p} {\left( \mathbf{x}^T \mathbf{B} \mathbf{x} \right)^q}, \quad \frac{\left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p} {\left( \mathbf{x}^T \mathbf{B} \mathbf{x} \right)^q \left( \mathbf{x}^T \mathbf{D} \mathbf{x} \right)^r}. \end{equation}\] In this package, these are called a simple and multiple ratio, respectively, following Smith (1989) (to be precise, the latter is a very specific case of the multiple ratios treated therein).

Quantities of the above forms frequently occur in the statistical and econometric literature (Mathai and Provost, 1992; Hillier et al., 2009, 2014). In evolutionary biology, the quantities known as evolvability measures have the same general form, and evaluation of their expectation (average or mean evolvability measures) can be of practical interest (Hansen and Houle, 2008; Melo et al., 2016; Watanabe, 2023).

The package is mainly for evaluating the moments of the above quantity assuming the multivariate normality of \(\mathbf{x}\), that is, \[\begin{equation} \qfrE \left( \frac{\left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p} {\left( \mathbf{x}^T \mathbf{B} \mathbf{x} \right)^q} \right), \quad \qfrE \left( \frac{\left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p} {\left( \mathbf{x}^T \mathbf{B} \mathbf{x} \right)^q \left( \mathbf{x}^T \mathbf{D} \mathbf{x} \right)^r} \right), \end{equation}\] where \(\mathbf{x} \sim \qfrmvnorm{\boldsymbol{\mu}}{\boldsymbol{\Sigma}}\).

## First examples

Here are simple example use cases:

    ## Simple matrices
    nv <- 4
    A <- diag(1:nv)
    B <- diag(sqrt(nv:1))
    D <- diag((nv:1) ^ 2)

    ## Expectation of (x^T A x)^2 / (x^T x)^1 where x ~ N(0, I)
    qfrm(A, p = 2, q = 1) # By default, B = I
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 26.66667
    #> This value is exact

    ## Expectation of (x^T A x)^2 / (x^T B x)^2 where x ~ N(0, I)
    qfrm(A, B, p = 2) # By default, q = p
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 3.467871, Error = 0 (one-sided)
    #> Possible range:
    #>  3.46787143 3.46787143

    ## Expectation of (x^T A x)^1/2 / (x^T B x)^1 where x ~ N(0, I)
    qfrm(A, B, p = 1/2, q = 1)
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 0.6652398
    #> Error bound unavailable; recommended to inspect plot() of this object

    ## Expectation of (x^T A x)^2 / ((x^T B x)^1 (x^T D x)^1) where x ~ N(0, I)
    qfmrm(A, B, D, p = 2, q = 1, r = 1)
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 1.135126
    #> Error bound unavailable; recommended to inspect plot() of this object

    ## Expectation of (x^T A x)^1 / (x^T B x)^1 where x ~ N(mu, Sigma)
    ## with arbitrary mu and Sigma
    mu <- 1:nv / nv
    Sigma <- diag(runif(nv) * 3)
    qfrm(A, B, p = 1, q = 1, mu = mu, Sigma = Sigma, m = 300)
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 2.210412, Error = 0 (two-sided)
    #> Possible range:
    #>  2.21041201 2.21041201

In each of the above examples, the function `qfrm()` (quadratic form ratio moment) or `qfmrm()` (quadratic form multiple ratio moment) is called to calculate the moment. The function calls are followed by outputs from the `print` method for the results, which show the values of the moments.

For the simple ratio, there is also a wrapper function `mqfr()` which returns the result from `qfrm()` as a numeric vector. This function assumes a special case \(p = q\), whose value is specified by its first argument `power` (against which it is vectorized, unlike `qfrm()`). This syntax is for consistency with other distribution-related functions in the package, `dqfr()`, `pqfr()`, `qqfr()`, and `rqfr()`, for which see `vignette("qfratio_distr")`.

    ## Expectation of (x^T A x)^2 / (x^T x)^2 where x ~ N(0, I)
    mqfr(2, A) # By default, B = I
    #> [1] 6.666667

    ## Expectation of (x^T A x)^p / (x^T B x)^p where x ~ N(0, I), p = 0,0.5,1,2
    mqfr(c(0, 0.5, 1, 2), A, B)
    #> [1] 1.000000 1.286306 1.724404 3.467871

    ## Expectation of (x^T A x)^1 / (x^T B x)^1 where x ~ N(mu, Sigma)
    ## with arbitrary mu and Sigma
    mqfr(1, A, B, mu = mu, Sigma = Sigma, m = 300)
    #> [1] 2.210412

You might have wondered what it means by results being “exact” or having “error”, or about the unexplained argument `m = 300` in the above examples. Seeing this requires some excursion into the math involved.

# Mathematical details

## Moment existence conditions

For the above quantities to be well-defined and have finite moments:

- \(\mathbf{B}\) (and \(\mathbf{D}\)) must be nonnegative definite (i.e., has no negative eigenvalues)
- when \(p\) is non-positive-integer, \(\mathbf{A}\) must be nonnegative definite as well

For the simple ratio, the necessary and sufficient conditions for the moment to exist is stated in Bao and Kan (2013, proposition 1):

1. When \(\mathbf{B}\) is positive definite, the moment exists if and only if \(\frac{n}{2} + p > q\).
2. When \(\mathbf{B}\) is positive semidefinite with rank \(l\), consider a matrix of its eigenvectors \(\mathbf{P} = (\mathbf{P}_1, \mathbf{P}_2)\), with \(\mathbf{P}_1\) and \(\mathbf{P}_2\) being \(n \times l\) and \(n \times (n - l)\) submatrices corresponding to the nonzero and zero eigenvalues, respectively, and let \(\tilde{\mathbf{A}}_{12} = \mathbf{P}_1^T \mathbf{A} \mathbf{P}_2\) and \(\tilde{\mathbf{A}}_{22} = \mathbf{P}_2^T \mathbf{A} \mathbf{P}_2\). Then,
   - when \(\tilde{\mathbf{A}}_{12} = \mathbf{0}_{l \times (n - l)}\) and \(\tilde{\mathbf{A}}_{22} = \mathbf{0}_{(n - l) \times (n - l)}\), the moment exists if and only if \(\frac{l}{2} + p > q\);
   - when \(\tilde{\mathbf{A}}_{12} \neq \mathbf{0}_{l \times (n - l)}\) and \(\tilde{\mathbf{A}}_{22} = \mathbf{0}_{(n - l) \times (n - l)}\), the moment exists if and only if \(\frac{l + p}{2} > q\); and
   - when \(\tilde{\mathbf{A}}_{22} \neq \mathbf{0}_{(n - l) \times (n - l)}\), the moment exists if and only if \(\frac{l}{2} > q\).

For the multiple ratio, the existence conditions are more difficult to define. When the null space of \(\mathbf{B}\) is a subspace of that of \(\mathbf{D}\) (or vice versa), the existence problem reduces to that of the simple ratio, and the above conditions apply by inserting the rank of \(\mathbf{B}\) and \(q + r\) into \(l\) and \(q\) above (Watanabe, 2023). However, when both \(\mathbf{B}\) and \(\mathbf{D}\) are singular and neither of their ranges is a subspace of each other, no general existence conditions seem to be known in the literature. See the documentations for `qfrm` and `qfmrm` for the behaviors when the existence conditions are not satisfied.

## Expressions for moments

There are several equivalent expressions for moments of ratios of quadratic forms in normal variables, only some of which are amenable to straightforward numerical evaluation. The package uses one class of such expressions, based on a series expansion and term-by-term integration of a ratio (Smith, 1989, 1993; Hillier, 2001; Hillier et al., 2009, 2014; Bao and Kan, 2013).

The following expressions assume \(\mathbf{x} \sim \qfrmvnorm{\boldsymbol{\mu}}{\mathbf{I}_n}\). For arbitrary but nonsingular \(\boldsymbol{\Sigma}\), it is always possible to transform the variables and argument matrices by using \(\mathbf{x}_{\mathrm{new}} = \mathbf{K}^{-1} \mathbf{x}\), \(\mathbf{A}_{\mathrm{new}} = \mathbf{K}^T \mathbf{A} \mathbf{K}\), etc., where \(\mathbf{K}\) is an \(n \times n\) matrix that satisfies \(\mathbf{K} \mathbf{K}^T = \boldsymbol{\Sigma}\) (Mathai and Provost, 1992, chap. 3). When \(\boldsymbol{\Sigma}\) is singular, certain conditions need to be met by the argument matrices, \(\boldsymbol{\Sigma}\), and \(\boldsymbol{\mu}\) for this transformation, and hence the following expressions, to be valid (Watanabe, 2023, appendix C).

For the simple ratio with an arbitrary (positive real) \(p\), when the moment exists, the expression is (Bao and Kan, 2013, eq. (12)) \[\begin{equation} \begin{aligned} \qfrE \left( \frac{(\mathbf{x}^T \mathbf{A} \mathbf{x})^p}{(\mathbf{x}^T \mathbf{B} \mathbf{x})^q} \right) = & 2^{p - q} \beta_{\mathbf{A}}^{-p} \beta_{\mathbf{B}}^{q} \qfrGmf{ \frac{n}{2} + p - q } \\ & \quad \cdot \sum_{i=0}^{\infty} \sum_{j=0}^{\infty} \frac{ \qfrrf[i]{-p} \qfrrf[j]{q} }{ \qfrGmf{\frac{n}{2} + i + j } } \\ & \qquad \cdot \qfrhij[i,j]{ \mathbf{I}_n - \beta_{\mathbf{A}} \mathbf{A} }{ \mathbf{I}_n - \beta_{\mathbf{B}} \mathbf{B} } , \end{aligned} \tag{1} \end{equation}\] and for a positive-integral \(p\), this can be simplified into (Bao and Kan, 2013, eq. (13), correcting a typo): \[\begin{equation} \begin{aligned} \qfrE \left( \frac{(\mathbf{x}^T \mathbf{A} \mathbf{x})^p}{(\mathbf{x}^T \mathbf{B} \mathbf{x})^q} \right) = & 2^{p - q} \beta_{\mathbf{B}}^{q} p! \qfrGmf{ \frac{n}{2} + p - q } \sum_{j=0}^{\infty} \frac{ \qfrrf[j]{q} }{ \qfrGmf{ \frac{n}{2} + p + j } } \qfrhtij[p;j]{ \mathbf{A} }{ \mathbf{I}_n - \beta_{\mathbf{B}} \mathbf{B} } . \end{aligned} \tag{2} \end{equation}\]

For the multiple ratio, the expression for an arbitrary \(p\) is (after Smith, 1989, eq. (5.2)) \[\begin{equation} \begin{aligned} \qfrE \! \left( \frac{(\mathbf{x}^T \mathbf{A} \mathbf{x})^p}{(\mathbf{x}^T \mathbf{B} \mathbf{x})^q (\mathbf{x}^T \mathbf{D} \mathbf{x})^r} \right) = & 2^{p - q - r} \beta_{\mathbf{A}}^{-p} \beta_{\mathbf{B}}^{q} \beta_{\mathbf{D}}^{r} \qfrGmf{ \frac{n}{2} + p - q - r } \\ & \quad \cdot \sum_{i=0}^{\infty} \sum_{j=0}^{\infty} \sum_{k=0}^{\infty} \frac{ \qfrrf[i]{-p} \qfrrf[j]{q} \qfrrf[k]{r} }{ \qfrGmf{ \frac{n}{2} + i + j + k } } \\ & \quad \quad \cdot \qfrhijk[i, j, k]{ \mathbf{I}_n - \beta_{\mathbf{A}} \mathbf{A} }{ \mathbf{I}_n - \beta_{\mathbf{B}} \mathbf{B} }{ \mathbf{I}_n - \beta_{\mathbf{D}} \mathbf{D} }, \end{aligned} \tag{3} \end{equation}\] and for a positive-integral \(p\) (Watanabe, 2023, eq. (B28)), \[\begin{equation} \begin{aligned} \qfrE \! \left( \frac{(\mathbf{x}^T \mathbf{A} \mathbf{x})^p}{(\mathbf{x}^T \mathbf{B} \mathbf{x})^q (\mathbf{x}^T \mathbf{D} \mathbf{x})^r} \right) = & 2^{p - q - r} \beta_{\mathbf{B}}^{q} \beta_{\mathbf{D}}^{r} p! \qfrGmf{ \frac{n}{2} + p - q - r } \\ & \quad \cdot \sum_{j=0}^{\infty} \sum_{k=0}^{\infty} \frac{ \qfrrf[j]{q} \qfrrf[k]{r} }{ \qfrGmf{ \frac{n}{2} + p + j + k } } \\ & \quad \quad \cdot \qfrhtijk[p; j,k]{ \mathbf{A} }{ \mathbf{I}_n - \beta_{\mathbf{B}} \mathbf{B} }{ \mathbf{I}_n - \beta_{\mathbf{D}} \mathbf{D} }. \end{aligned} \tag{4} \end{equation}\]

In these expressions, \(\beta_{\cdot}\) are arbitrary scaling constants that satisfy \(0 < \beta < 2 / \qfrlmax\), with \(\qfrlmax\) being the largest eigenvalue of the argument matrix. In addition, \(h_{k_1, \dots k_s}\) and \(\tilde{h}_{k_1; k_2 \dots k_s}\) are the coefficients of \(t_1^{k_1} \dots t_s^{k_s}\) in the following power series expansions: \[\begin{multline} \det{ \mathbf{I}_n - t_1 \mathbf{A}_1 - \dots - t_s \mathbf{A}_s }^{-\frac{1}{2}} \\ \cdot \exp \left( \frac{ (1 - t_1 - \dots - t_s) \boldsymbol{\mu}^T \left( \mathbf{I}_n - t_1 \mathbf{A}_1 - \dots - t_s \mathbf{A}_s \right)^{-1} \boldsymbol{\mu} - \boldsymbol{\mu}^T \boldsymbol{\mu} }{2} \right) \\ = \sum_{k_1 = 0}^{\infty} \dots \sum_{k_s = 0}^{\infty} \qfrhijk[k_1, \dots, k_s]{ \mathbf{A}_1 }{\dots}{ \mathbf{A}_s } t_1^{k_1} \dots t_s^{k_s}, \label{gfun_hij} \\ \end{multline}\] \[\begin{multline} \det{ \mathbf{I}_n - t_1 \mathbf{A}_1 - \dots - t_s \mathbf{A}_s }^{-\frac{1}{2}} \\ \cdot \exp \left( \frac{ (1 - t_2 - \dots - t_s) \boldsymbol{\mu}^T \left( \mathbf{I}_n - t_1 \mathbf{A}_1 - \dots - t_s \mathbf{A}_s \right)^{-1} \boldsymbol{\mu} - \boldsymbol{\mu}^T \boldsymbol{\mu} }{2} \right) \\ = \sum_{k_1 = 0}^{\infty} \dots \sum_{k_s = 0}^{\infty} \qfrhtijk[k_1; k_2, \dots, k_s]{ \mathbf{A}_1 }{ \mathbf{A}_2, \dots }{ \mathbf{A}_s } t_1^{k_1} t_2^{k_2} \dots t_s^{k_s}. \end{multline}\] (These notations follow Bao and Kan (2013); Hillier et al. (2014) used the symbol \(h\) for \(\tilde{h}\) here.)

A common special case is when \(\boldsymbol{\mu} = \mathbf{0}_n\), in which case the coefficients \(d\) in the following expansion is used: \[\begin{equation} \det{ \mathbf{I}_n - t_1 \mathbf{A}_1 - \dots - t_s \mathbf{A}_s }^{-\frac{1}{2}} = \sum_{k_1 = 0}^{\infty} \dots \sum_{k_s = 0}^{\infty} \qfrdijk[k_1, \dots, k_s]{ \mathbf{A}_1 }{\dots}{ \mathbf{A}_s } t_1^{k_1} \dots t_s^{k_s} . \end{equation}\]

Actually, the above expressions (1)–(4) can be written in terms of \(d\) with an additional argument (Hillier et al., 2009; Bao and Kan, 2013), but those alternative expressions are not utilized in the package.

A theoretically important result is that this \(d\) is related to the top-order zonal and invariant polynomials of matrix arguments, \(C\) (Chikuse, 1987; Hillier et al., 2009): \[\begin{align} \qfrdijk[k_1, \dots, k_s]{ \mathbf{A}_1 }{\dots}{ \mathbf{A}_s } = \frac{ \qfrrf[k_1 + \dots + k_s]{\frac{1}{2}} }{ k_1 ! \dots k_s !} \qfrCit{\qfrbrc{k_1}}{\dots}{\qfrbrc{k_s}}{\qfrbrc{k_1 + \dots + k_s}}{ \mathbf{A}_1 }{\dots}{ \mathbf{A}_s } , \end{align}\] where \(\qfrbrc{k_i}\) and the like denote the top-order partitions of \(k_i\), in the notation of Smith (1989).

Note that the order of the argument matrices is arbitrary in \(d\) or \(h\), but not in \(\tilde{h}\) where \(\mathbf{A}_1\) is not interchangeable with the other arguments.

## Special cases

A notable special case for moments of simple ratios is when \(p\) is a positive integer and \(\mathbf{B} = \mathbf{I}_n\). In this case, (2) simplifies into (Hillier et al., 2014, theorem 4) \[\begin{equation} \qfrE \left( \frac{(\mathbf{x}^T \mathbf{A} \mathbf{x})^p}{(\mathbf{x}^T \mathbf{x})^q} \right) = 2^{p - q} p! \sum_{k=0}^{p} \frac{ \qfrGmf{ \frac{n}{2} + p - q + k } }{ 2^k k! \qfrGmf{ \frac{n}{2} + p + k } } {}_1 F_1 \left( q ; \frac{n}{2} + p + k ; - \frac{\boldsymbol{\mu}^T \boldsymbol{\mu}}{2} \right) a_{p, k} \left( \mathbf{A}, \boldsymbol{\mu} \right) , \tag{5} \end{equation}\] where \({}_1 F_1 \left( \cdot ; \cdot ; \cdot \right)\) is the confluent hypergeometric function, and \(a_{p, k}\) are the coefficients of \(t^p\) in \[\begin{align} & \det{ \mathbf{I}_n - t \mathbf{A} }^{-\frac{1}{2}} \left( \boldsymbol{\mu}^T \left( \mathbf{I}_n - t \mathbf{A} \right)^{-1} \boldsymbol{\mu} - \boldsymbol{\mu}^T \boldsymbol{\mu} \right) ^ k = \sum_{i = k}^{\infty} a_{p, i} \left( \mathbf{A}, \boldsymbol{\mu} \right) t^{i}. \end{align}\] As the lowest-order term in this last expansion is that of \(t^k\), \(a_{p, k} = 0\) for \(k > p\), so that the series in (5) truncates at \(k = p\). This truncation enables quasi-exact evaluation of the moment, given accurate evaluation of the hypergeometric function.

When \(\boldsymbol{\mu} = \mathbf{0}_n\) in addition, (5) further simplifies as \(a_{p, k} = 0\) for \(k > 0\), so that the series essentially vanishes leaving only the term for \(k = 0\). Similar simplification happens in (1), (3), and (4) as well when any of the argument matrices is \(\mathbf{I}_n\) and \(\boldsymbol{\mu} = \mathbf{0}_n\), although they will still involve a single or double infinite series.

## Numerical evaluation

The above expressions (1)–(4) are exact, but their exact values cannot generally be evaluated as they involve infinite series, except when the conditions for (5) apply. Fortunately, the series are absolutely convergent, so that their partial sums up to certain orders typically yield a sufficiently accurate result. This is what is done by the `qfrm()` and `qfmrm()` functions. The order of evaluation is determined by the argument `m`; the functions calculate the terms from \(i = j = k = 0\) to \(i + j + k \leq m\) (as \(i,j,k\) appear in the above summations) and then return an \(m + 1\) vector that contains sums of the same-order terms, along with other things.

The default is `m = 100`, which sufficed for most of the above examples, but for large problems, the user would want to use a larger value (which can be hundreds of thousands for very large problems; see below).

The coefficients \(d\), \(h\), \(\tilde{h}\) are evaluated by recursive algorithms described below. At the core of this package are internal functions to implement these recursive algorithms, which can be computationally intensive as the problem becomes larger.

## Truncation error

A great advantage in the above expressions is that an error bound is available for a partial (truncated) sum for the simple ratio, provided that \(p\) is a positive integer and \(\mathbf{B}\) is positive definite. By denoting the expression (2) as \(M \left( \mathbf{A}, \mathbf{B}, p, q \right)\) and the partial sum of the same up to \(j = m\) as \(\hat{M}_m \left( \mathbf{A}, \mathbf{B}, p, q \right)\), the error bound is (Hillier et al., 2014, theorem 7) \[\begin{multline} \lvert M \left( \mathbf{A}, \mathbf{B}, p, q \right) - \hat{M}_m \left( \mathbf{A}, \mathbf{B}, p, q \right) \rvert \\ \leq \frac{ 2^{p - q} \beta_{\mathbf{B}}^{q} p! \qfrGmf{ \frac{n}{2} + p - q } \qfrrf[m+1]{q} } { \qfrGmf{ \frac{n}{2} + p + m + 1 } } \left[ \frac{ \exp \left( \frac{ \bar{\boldsymbol{\mu}}^T \bar{\boldsymbol{\mu}} - \boldsymbol{\mu}^T \boldsymbol{\mu} }{2} \right) \qfrdtk[p]{ \bar{A} }{ \bar{\boldsymbol{\mu}} } }{ \det{\beta_{\mathbf{B}} \mathbf{B}}^{\frac{1}{2}} } - \sum_{j=0}^{m} \qfrhhij[p;j]{ \mathbf{A}^+ }{ \mathbf{I}_n - \beta_{\mathbf{B}} \mathbf{B} } \right] , \end{multline}\] where \(\mathbf{A}^+\) is a symmetric matrix constructed from the eigenvectors and “positivized” eigenvalues of \(\mathbf{A}\) (above), \(\bar{\boldsymbol{\mu}} = \sqrt{2}\left( \beta_\mathbf{B} \mathbf{B} \right)^{-\frac{1}{2}} \boldsymbol{\mu}\), \(\bar{\mathbf{A}} = \beta_{\mathbf{B}}^{-1} \mathbf{B}^{-\frac{1}{2}} \mathbf{A}^+ \mathbf{B}^{-\frac{1}{2}}\), and \(\tilde{d}\) and \(\hat{h}\) are coefficients in the following generating functions: \[\begin{equation} \det{ \mathbf{I}_n - t \mathbf{A} }^{-\frac{1}{2}} \exp \left( \frac{ \boldsymbol{\mu}^T \left( \mathbf{I}_n - t \mathbf{A} \right)^{-1} \boldsymbol{\mu} - \boldsymbol{\mu}^T \boldsymbol{\mu} }{2} \right) = \sum_{k = 0}^{\infty} \qfrdtk[k]{ \mathbf{A} }{ \boldsymbol{\mu} } t^{k} , \end{equation}\] \[\begin{multline} \det{ \mathbf{I}_n - t_1 \mathbf{A}_1 - \dots - t_s \mathbf{A}_s }^{-\frac{1}{2}} \\ \quad \cdot \exp \left( \frac{ (1 + t_2 + \dots + t_s) \boldsymbol{\mu}^T \left( \mathbf{I}_n - t_1 \mathbf{A}_1 - \dots - t_s \mathbf{A}_s \right)^{-1} \boldsymbol{\mu} - \boldsymbol{\mu}^T \boldsymbol{\mu} }{2} \right) \\ = \sum_{k_1 = 0}^{\infty} \dots \sum_{k_s = 0}^{\infty} \qfrhhijk[k_1; k_2, \dots, k_s]{ \mathbf{A}_1 }{ \mathbf{A}_2, \dots }{ \mathbf{A}_s } t_1^{k_1} t_2^{k_2} \dots t_s^{k_s} . \end{multline}\]

For the multiple ratio, a similar error bound is available when \(\mathbf{D} = \mathbf{I}_n\), in addition to \(p\) being a positive integer. By denoting the expression (4) as \(M \left( \mathbf{A}, \mathbf{B}, \mathbf{D}, p, q, r \right)\), and the partial sum of (4) up to \(j + k = m\) as \(\hat{M}_m \left( \mathbf{A}, \mathbf{B}, \mathbf{D}, p, q, r \right)\), the error bound is \[\begin{equation} \begin{aligned} & \lvert M \left( \mathbf{A}, \mathbf{B}, \mathbf{I}_n, p, q, r \right) - \hat{M}_m \left( \mathbf{A}, \mathbf{B}, \mathbf{I}_n, p, q, r \right) \rvert \\ & \quad \leq \frac{ 2^{p - q - r} \beta_{\mathbf{B}}^{q} p! \qfrGmf{ \frac{n}{2} + p - q - r } \qfrrf[m+1]{\max \left( q, r \right)} } { \qfrGmf{ \frac{n}{2} + p + m + 1 } } \\ & \quad \quad \cdot \left[ \frac{ \exp \left( \frac{ \tilde{\boldsymbol{\mu}}^T \tilde{\boldsymbol{\mu}} - \boldsymbol{\mu}^T \boldsymbol{\mu} }{2} \right) \qfrdtk[p]{ \bar{A} }{ \tilde{\boldsymbol{\mu}} } }{ \det{\beta_{\mathbf{B}} \mathbf{B}}^{\frac{1}{2}} } - \sum_{k=0}^{m} \sum_{j=0}^{k} \qfrhhijk[p;k-j,k]{ \mathbf{A}^+ }{ \mathbf{I}_n - \beta_{\mathbf{B}} \mathbf{B} }{ \mathbf{0}_{n \times n} } \right] , \end{aligned} \end{equation}\] where \(\tilde{\boldsymbol{\mu}} = \sqrt{3}\left( \beta_\mathbf{B} \mathbf{B} \right)^{-\frac{1}{2}} \boldsymbol{\mu}\) and other notations are as above.

## Recursions

The coefficients \(h\), \(\tilde{h}\), and \(\hat{h}\) in the above series expressions are evaluated by the “super-short” recursive algorithm developed by Hillier and colleagues (Bao and Kan, 2013; Hillier et al., 2014). That is, \[\begin{align} h_{i,j,k} = & \frac{\qfrtr{\mathbf{G}_{i,j,k}} + \boldsymbol{\mu}^T \mathbf{g}_{i,j,k}}{2(i + j + k)}, \\ \tilde{h}_{i,j,k} = & \frac{\qfrtr{\tilde{\mathbf{G}}_{i,j,k}} + \boldsymbol{\mu}^T \tilde{\mathbf{g}}_{i,j,k}}{2(i + j + k)}, \\ \hat{h}_{i,j,k} = & \frac{\qfrtr{\hat{\mathbf{G}}_{i,j,k}} + \boldsymbol{\mu}^T \hat{\mathbf{g}}_{i,j,k}}{2(i + j + k)}, \end{align}\] where \[\begin{align} \mathbf{G}_{i,j,k} = & \mathbf{A}_1 \left( h_{i-1,j,k} \mathbf{I}_n + \mathbf{G}_{i-1,j,k} \right) \\ & + \mathbf{A}_2 \left( h_{i,j-1,k} \mathbf{I}_n + \mathbf{G}_{i,j-1,k} \right) \\ & + \mathbf{A}_3 \left( h_{i,j,k-1} \mathbf{I}_n + \mathbf{G}_{i,j,k-1} \right) , \\ \mathbf{g}_{i,j,k} = & \left( \mathbf{G}_{i,j,k} - \mathbf{G}_{i-1,j,k} - \mathbf{G}_{i,j-1,k} - \mathbf{G}_{i,j,k-1} \right) \boldsymbol{\mu} \\ & - \left( h_{i-1,j,k} + h_{i,j-1,k} + h_{i,j,k-1} \right) \boldsymbol{\mu} \\ & + \mathbf{A}_1 \mathbf{g}_{i-1,j,k} + \mathbf{A}_2 \mathbf{g}_{i,j-1,k} + \mathbf{A}_3 \mathbf{g}_{i,j,k-1} , \\ \tilde{\mathbf{G}}_{i,j,k} = & \mathbf{A}_1 \left( \tilde{h}_{i-1,j,k} \mathbf{I}_n + \tilde{\mathbf{G}}_{i-1,j,k} \right) \\ & + \mathbf{A}_2 \left( \tilde{h}_{i,j-1,k} \mathbf{I}_n + \tilde{\mathbf{G}}_{i,j-1,k} \right) \\ & + \mathbf{A}_3 \left( \tilde{h}_{i,j,k-1} \mathbf{I}_n + \tilde{\mathbf{G}}_{i,j,k-1} \right) , \\ \tilde{\mathbf{g}}_{i,j,k} = & \left( \tilde{\mathbf{G}}_{i,j,k} - \tilde{\mathbf{G}}_{i,j-1,k} - \tilde{\mathbf{G}}_{i,j,k-1} \right) \boldsymbol{\mu} \\ & - \left( \tilde{h}_{i,j-1,k} + \tilde{h}_{i,j,k-1} \right) \boldsymbol{\mu} \\ & + \mathbf{A}_1 \tilde{\mathbf{g}}_{i-1,j,k} + \mathbf{A}_2 \tilde{\mathbf{g}}_{i,j-1,k} + \mathbf{A}_3 \tilde{\mathbf{g}}_{i,j,k-1} , \\ \hat{\mathbf{G}}_{i,j,k} = & \mathbf{A}_1 \left( \hat{h}_{i-1,j,k} \mathbf{I}_n + \hat{\mathbf{G}}_{i-1,j,k} \right) \\ & + \mathbf{A}_2 \left( \hat{h}_{i,j-1,k} \mathbf{I}_n + \hat{\mathbf{G}}_{i,j-1,k} \right) \\ & + \mathbf{A}_3 \left( \hat{h}_{i,j,k-1} \mathbf{I}_n + \hat{\mathbf{G}}_{i,j,k-1} \right) , \\ \hat{\mathbf{g}}_{i,j,k} = & \left( \hat{\mathbf{G}}_{i,j,k} + \hat{\mathbf{G}}_{i,j-1,k} + \hat{\mathbf{G}}_{i,j,k-1} \right) \boldsymbol{\mu} \\ & + \left( \hat{h}_{i,j-1,k} + \hat{h}_{i,j,k-1} \right) \boldsymbol{\mu} \\ & + \mathbf{A}_1 \hat{\mathbf{g}}_{i-1,j,k} + \mathbf{A}_2 \hat{\mathbf{g}}_{i,j-1,k} + \mathbf{A}_3 \hat{\mathbf{g}}_{i,j,k-1} , \\ \end{align}\] with the initial conditions \(h_{0,0,0} = \tilde{h}_{0,0,0} = \hat{h}_{0,0,0} = 1\) (as per Hillier et al. (2014), but not \(0\), as incorrectly printed in Bao and Kan (2013)), \(\mathbf{G}_{0,0,0} = \tilde{\mathbf{G}}_{0,0,0} = \hat{\mathbf{G}}_{0,0,0} = \mathbf{0}_{n \times n}\), and \(\mathbf{g}_{0,0,0} = \tilde{\mathbf{g}}_{0,0,0} = \hat{\mathbf{g}}_{0,0,0} = \mathbf{0}_n\). If any of the indices is negative, the term is treated as \(0\).

The recursion for \(d\) can be obtained by just inserting \(\boldsymbol{\mu} = \mathbf{0}_n\) (and hence omitting \(\mathbf{g}\) altogether) in that for \(h\), and those for the two-argument versions are by omitting any terms with \(\mathbf{A}_3\) or \(k - 1\) as an index in the corresponding three-argument versions.

The recursion for \(\qfrdtk{\mathbf{A}}{\boldsymbol{\mu}}\) is (Bao and Kan, 2013; Hillier et al., 2014): \[\begin{equation} \tilde{d}_{k} = \frac{1}{2k} \sum_{i=1}^n \left(u_{k,i} + v_{k,i} \right), \end{equation}\] where \[\begin{align} u_{k,i} = & \lambda_i \left( \tilde{d}_{k-1} + u_{k-1,i} \right), \\ v_{k,i} = & \delta_i u_{k,i} + \lambda_i v_{k-1,i}, \end{align}\] with \(\lambda_i\) being the \(i\)th eigenvalue of \(\mathbf{A}\), \(\delta_i = \left( \mathbf{h}_i^T \boldsymbol{\mu} \right)^2\), \(\mathbf{h}_i\) being the eigenvector of \(\mathbf{A}\) corresponding to \(\lambda_i\), and the initial conditions are \(\tilde{d}_0 = 1\) and \(u_{0,i} = v_{0,i} = 0\) for all \(i\).

The recursion for \(a_{p,k} \left( \mathbf{A}, \boldsymbol{\mu} \right)\) is (Hillier et al., 2014, eqs. (31), (32), noting that \(w\) should be indexed by \(k\) as well): \[\begin{equation} a_{p,k} = \sum_{i=1}^n \delta_i w_{p,k,i}, \end{equation}\] where \[\begin{equation} w_{p,k,i} = \lambda_i \left( a_{p-1,k-1} + w_{p-1,k,i} \right), \end{equation}\] with \(a_{p,0} = \qfrdk[p]{\mathbf{A}}\) and \(w_{0,k,i} = 0\) for any \(k, i\).

The recursion for \(\qfrdk{\mathbf{A}}\) can be obtained either as that for \(\qfrdtk{\mathbf{A}}{\mathbf{0}_n}\) (where all \(v\)’s are \(0\)) or by dropping any terms with \(\mathbf{A}_2\) or \(j-1\) as an index from that for the two-argument version of \(d\) above.

# Functions

The package is designed for evaluation of moments of ratios of quadratic forms using the recursive algorithms described above. The primary front-end functions are `qfrm()` and `qfmrm()`, which return a `qfrm` object containing terms in the truncated series and error bound when available. Simple `print` and `plot` methods are defined for the `qfrm` class, for quick inspection of evaluation results.

In addition, there is a convenience wrapper `mqfr()`, which extracts and returns the main evaluation result from `qfrm()` as a numeric vector (see `vignette("qfratio_distr")`).

## `qfrm()`: quadratic form ratio moment

This function is for evaluating moments of the simple ratios of quadratic forms, \(\qfrE \left( \frac{\left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p} {\left( \mathbf{x}^T \mathbf{B} \mathbf{x} \right)^q} \right)\), where \(\mathbf{x} \sim \qfrmvnorm{\boldsymbol{\mu}}{\boldsymbol{\Sigma}}\).

#### Usage

    qfrm(A, B, p = 1, q = p, m = 100L,
         mu = rep.int(0, n), Sigma = diag(n),
         tol_zero = .Machine$double.eps * 100,
         tol_sing = tol_zero, simplify_ratio = TRUE, ...)

`A` and `B` are the argument matrices; when one is missing, \(\mathbf{I}_n\) is inserted.

`m` determines the order of evaluation, so will be one of the most important arguments specified by the user.

`mu` and `Sigma` specify the mean and covariance of \(\mathbf{x}\), so by default \(\mathbf{x} \sim \qfrmvnorm{\mathbf{0}_n}{\mathbf{I}_n}\).

`tol_zero` and `tol_sing` just determine the numerical tolerance for assessing identity and singularity, respectively, so usually need not be changed except when handling edgy cases.

The optional argument `simplify_ratio` is to specify the behavior in special conditions when `A == B` or `n == 1`. When `TRUE` and `A == B` or `n == 1`, the function tries to simplify the problem by partially cancelling the exponents `p` and `q`, so that the new `p` is the smallest nonnegative integer that ensures `q` to be nonnegative. Such a combination of `p` and `q` usually allows for efficient evaluation of the series expression and calculation of error bound (below).

#### Inside the function

When `Sigma` is provided and not \(\mathbf{I}_n\), this function calls itself recursively with `A`, `B`, and `mu` transformed as described above (and dropping `Sigma`).

As the evaluation is based on different series expressions depending on whether \(\mathbf{B} = \mathbf{I}_n\) or not and whether \(p\) is an integer or not, the actual calculation is done by one of the following four “internal” functions, for which `qfrm()` acts as a switcher:

Function | Argument matrices | \(p\) | Expression | Error bound | Number of recursions
--- | --- | --- | --- | --- | ---
`qfrm_ApIq_int()` | \(\mathbf{B} = \mathbf{I}_n\) | integer | (5) | Exact* | \(\sim p\)
`qfrm_ApIq_npi()` | \(\mathbf{B} = \mathbf{I}_n\) | non-positive-integer | (1) (simplified when \(\boldsymbol{\mu} = \mathbf{0}_n\)) | Available only when \(\boldsymbol{\mu} = \mathbf{0}_n\) | \(\sim m\) when \(\boldsymbol{\mu} = \mathbf{0}_n\); \(\sim m^2 / 2\) otherwise
`qfrm_ApBq_int()` | \(\mathbf{B}\) arbitrary | integer | (2) | Available | \(\sim (p + 1)m\)
`qfrm_ApBq_npi()` | \(\mathbf{B}\) arbitrary | non-positive-integer | (1) | Unavailable | \(\sim m^2 / 2\)

*: assuming accurate evaluation of the confluent hypergeometric function when \(\boldsymbol{\mu} \neq \mathbf{0}_n\)

These functions are exported into the package namespace so can be directly used if necessary.

They do the following operations:

1. Check structures of arguments (e.g., `A`, `B` being square, `p`, `q` being scalar)
2. Symmetrize the argument matrices (`A <- (A + t(A)) / 2`; same for `B`)
3. Rotate `A`, `B`, and `mu` with the eigenvectors of `B`, to
   - check for the moment existence conditions (above)
   - see whether `A` is co-diagonalized, in which case computation can be simplified
4. Calculate \(d\), \(h\), or \(\tilde{h}\) up to the order `m` with an internal function
5. Multiply the result with the appropriate coefficients in (1), (2) or (5) with an internal function
6. In `qfrm_ApBq_npi()`, the same-order terms in the double series are summed up
7. Conduct a rough check of numerical convergence
8. Calculate the sequence of error bound as in the steps 4–6 when applicable
9. Return a structured object including the series and error bound

By default, the steps 4–8 are handled by internal `C++` functions. The steps 7 and 8 can be suppressed if desired (by using `check_convergence = "none"` or `FALSE` and `error_bound = FALSE`, respectively).

## `qfmrm()`: quadratic form multiple ratio moment

This function is for evaluating moments of the multiple ratios of quadratic forms, \(\qfrE \left( \frac{\left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p} {\left( \mathbf{x}^T \mathbf{B} \mathbf{x} \right)^q\left( \mathbf{x}^T \mathbf{D} \mathbf{x} \right)^r} \right)\), where \(\mathbf{x} \sim \qfrmvnorm{\boldsymbol{\mu}}{\boldsymbol{\Sigma}}\).

It works in a way very similar to `qfrm()`, just with the additional arguments `D` and `r`. The corresponding “internal” functions are:

Function | Argument matrices | \(p\) | Expression | Error bound | Number of recursions
--- | --- | --- | --- | --- | ---
`qfmrm_ApBIqr_int()` | \(\mathbf{D} = \mathbf{I}_n\) | integer | (4) (simplified when \(\boldsymbol{\mu} = \mathbf{0}_n\)) | Available | \(\sim (p + 1) m\) when \(\boldsymbol{\mu} = \mathbf{0}_n\); \(\sim (p + 1) m^2 / 2\) otherwise
`qfmrm_ApBIqr_npi()` | \(\mathbf{D} = \mathbf{I}_n\) | non-positive-integer | (3) (simplified when \(\boldsymbol{\mu} = \mathbf{0}_n\)) | Unavailable | \(\sim m^2 / 2\) when \(\boldsymbol{\mu} = \mathbf{0}_n\); \(\sim m^3 / 6\) otherwise
`qfmrm_IpBDqr_gen()` | \(\mathbf{A} = \mathbf{I}_n\) | general | (4) simplified (when \(\boldsymbol{\mu} = \mathbf{0}_n\)) or (3) (otherwise) | Unavailable | \(\sim m^2 / 2\) when \(\boldsymbol{\mu} = \mathbf{0}_n\); \(\sim m^3 / 6\) otherwise
`qfmrm_ApBDqr_int()` | Arbitrary | integer | (4) | Unavailable | \(\sim (p + 1) m^2 / 2\)
`qfmrm_ApBDqr_npi()` | Arbitrary | non-positive-integer | (3) | Unavailable | \(\sim m^3 / 6\)

#### 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, simplify_ratio = TRUE, ...)

Essentially the same as `qfrm()`.

## `qfm_Ap_int()`, `qfpm_ABpq_int()`, `qfpm_ABDpqr_int()`: quadratic form (product) moment

These functions are for evaluating moments of (products of) quadratic forms, \(\qfrE \left( \left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p \right)\), \(\qfrE \left( \left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p \left( \mathbf{x}^T \mathbf{B} \mathbf{x} \right)^q \right)\), and \(\qfrE \left( \left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p \left( \mathbf{x}^T \mathbf{B} \mathbf{x} \right)^q \left( \mathbf{x}^T \mathbf{D} \mathbf{x} \right)^r \right)\), where \(\mathbf{x} \sim \qfrmvnorm{\boldsymbol{\mu}}{\boldsymbol{\Sigma}}\), using the algorithm described in Hillier et al. (2014, secs. 3–4). The exponents must be positive integers.

In principle, the same algorithm can be applied to a product of more than three quadratic forms, but this package has not implemented functions for them.

As these expressions only involve finite series, all the results are exact. The functions return a `qfpm` object, for which a simple `print` method is defined.

#### Usage

    qfpm_ABDpqr_int(A, B, D, p = 1, q = 1, r = 1,
                    mu = rep.int(0, n), Sigma = diag(n),
                    use_cpp = TRUE, cpp_method = "double",
                    tol_zero = .Machine$double.eps * 100,
                    tol_sing = tol_zero)

The other two functions have the same arguments, just dropping the inapplicable argument matrices and exponents.

`cpp_method` is a dummy argument and ignored in these functions.

## `rqfr()`, `rqfmr()`, `rqfp()`: random number generation

These functions generate a random sample of the simple ratio \(\frac{\left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p} {\left( \mathbf{x}^T \mathbf{B} \mathbf{x} \right)^q}\), multiple ratio \(\frac{\left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p} {\left( \mathbf{x}^T \mathbf{B} \mathbf{x} \right)^q\left( \mathbf{x}^T \mathbf{D} \mathbf{x} \right)^r}\), and product \(\left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p \left( \mathbf{x}^T \mathbf{B} \mathbf{x} \right)^q \left( \mathbf{x}^T \mathbf{D} \mathbf{x} \right)^r\), respectively, of quadratic forms in normal variables \(\mathbf{x} \sim \qfrmvnorm{\boldsymbol{\mu}}{\boldsymbol{\Sigma}}\).

They can be used for verifying results from the above functions, or as a Monte Carlo surrogate when numerical convergence cannot be obtained within reasonable computational memory and time.

#### Usage

    rqfr(nit, A, B, power = 1, q = power, mu = rep.int(0, n), Sigma = diag(n), use_cpp = TRUE)
    rqfmr(nit, A, B, D, p = 1, q = p/2, r = q, mu = rep.int(0, n), Sigma = diag(n), use_cpp = TRUE)
    rqfp(nit, A, B, D, p = 1, q = 1, r = 1, mu = rep.int(0, n), Sigma = diag(n), use_cpp = TRUE)

Similar to the regular random number generation functions, e.g., `stats::rnorm()`; the first argument `nit` is the desired size of the random sample (or number of iteration), followed by other parameters.

## `d1_i()`, `d2_*j_*()`, `d3_*jk_*()`, etc.: recursion algorithms

The package has a number of internal functions that implement the recursive algorithms described above. They are usually not accessed by the user, but are documented herein to aid understanding of the package structure. Note also that, by default, the actual calculations are handled by `C++` equivalents, which are not exposed to the user.

Arranged in the target variable and the number of argument matrices, these are:

Target | One matrix | Two matrices | Three matrices
--- | --- | --- | ---
\(d\) | `d1_i()` | `d2_ij_*()`, `d2_pj_*()` | `d3_ijk_*()`, `d3_pjk_*()`
\(\tilde{d}\) | `dtil1_i_*()` | `dtil2_pq_*()` | `dtil3_pqr_*()`
\(h\) | — | `h2_ij_*()` | `h3_ijk_*()`
\(\tilde{h}\) | — | `htil2_pj_*()` | `htil3_pjk_*()`
\(\hat{h}\) | — | `hhat2_pj_*()` | `hhat3_pjk_*()`
\(a\) | `a1_pk()` | — | —

The functions with `*_ij_*` or `*_ijk_*` are to calculate terms in double and triple infinite series, with the dimensions of the output determined by the argument `m` (typically used with non-integer \(p\), (1) or (3) above), whereas `*_pj_*` and `*_pjk_*` are for the terms in series truncated along the first argument matrix (used with integral \(p\), (2) or (4)). The one matrix functions assume either of these two cases, and `dtil2_pq_*()` and `dtil3_pqr_*()` assume finite series along all argument matrices.

All these functions, excepting `d1_i()` and `a1_pk()`, have two versions ending in `*_m` or `*_v`, which stand for the matrix and vector “methods”. The latter is for the situations where all the argument matrices are co-diagonalized by the same eigenvectors, in which case the above recursive algorithms reduce from full matrix arithmetic to element-wise operations of the eigenvalues (and the rotated mean vector), which is much less computationally intensive.

Most of these functions automatically scale the result to avoid numerical overflow, and (log of) the scaling factors are appended as an attribute of the return value, so that the front-end functions can appropriately scale the results back to the original one (see below).

# Using the functionality

## Recommended workflow

This package implements the recursive algorithms with `C++` functions with the aid of the `RcppEigen` package to accelerate calculations. The evaluation of moments can still be computationally intensive and take non-trivial computational time, especially when a double or triple series is involved (see (1), (3) and (4)).

Unless an exact solution is available (i.e., when using (5)), it is always recommended to check for numerical convergence by inspecting the error bound and the `plot` method. Thus, a recommended workflow is:

1. Use `qfrm()` (or `qfmrm()`) with the desired arguments with a small `m` (perhaps the default `m = 100`), assigning the result into an object
   - This may yield a warning when non-convergence is evident
2. Inspect the result with the `print` method (or list element access), paying attention to the error bound whenever available
3. Use `plot()` of the object to see if numerical convergence has been reached
4. If unsatisfied, repeat the cycle with a larger `m`

Here are some examples:

    ## Preparing simple objects
    nv <- 10
    A <- diag(1:nv)
    B <- diag(sqrt(nv:1))

    #####
    ## I: Simple case, (x^T A t)^2 / (x^T x)^2, where x ~ N(0, I)
    ##    Exact solution available
    res <- qfrm(A, p = 2, q = 2)
    res
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 31.625
    #> This value is exact
    ## In this case, plot is moot because the result is exact

    #####
    ## II: Simple case, (x^T A t)^2 / (x^T B x)^2, where x ~ N(0, I)
    ##     Error bound available
    res <- qfrm(A, B, 2, 2)
    res
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 7.319712, Error = 8.470329e-21 (one-sided)
    #> Possible range:
    #>  7.31971208 7.31971208

    ## Displaying more digits
    print(res, digits = 10)
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 7.319712084, Error = 8.470329473e-21 (one-sided)
    #> Possible range:
    #>  7.31971208375 7.31971208375

    ## Direct access to the results:
    names(res)
    #> [1] "statistic"   "error_bound" "terms"       "seq_error"
    res$statistic       # Evaluation result (partial sum)
    #> [1] 7.319712
    res$error_bound     # Error bound for $statistic
    #> [1] 8.470329e-21
    tail(res$terms)     # Tail of truncated series; sum($terms) == $statistic
    #> [1] 3.543891e-20 2.339201e-20 1.544571e-20 1.020231e-20 6.741185e-21
    #> [6] 4.455725e-21
    tail(res$seq_error) # Tail of sequence of error bounds; last one is $error_bound
    #> [1] 7.877406e-20 5.251604e-20 4.150461e-20 2.879912e-20 1.863472e-20
    #> [6] 8.470329e-21

    ## Inspect plot
    plot(res)

    ## Note that profile for error bound is shown as
    ## upper bound of the true value in the plot

    #####
    ## III: Simple case, (x^T A t)^(1/2) / (x^T x)^(1/2), where x ~ N(0, I)
    ##      Error bound can be negative
    res <- qfrm(A, p = 1/2, q = 1/2)
    res # Note negative (though one-sided) error bound
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 2.3313, Error = -8.091943e-14 (one-sided)
    #> Possible range:
    #>  2.33129983 2.33129983

    ## Inspect plot
    plot(res)

    #####
    ## IV: Non-central case, (x^T A t)^2 / (x^T B x)^2, where x ~ N(mu, I)
    ##     Two-sided error bound
    mu <- 1:nv / sum(1:nv)
    res <- qfrm(A, B, p = 2, q = 2, mu = mu)
    res # Note two-sided error bound
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 7.428816, Error = 3.489776e-19 (two-sided)
    #> Possible range:
    #>  7.42881636 7.42881636

    ## plot() automatically handles two-sided error bound
    plot(res)

    #####
    ## V: Difficult case with highly dispersed eigenvalues
    ##    Error bound available
    A <- diag((1:nv) ^ 2)
    B <- diag((nv:1) ^ 3)

    res <- qfrm(A, B, 2, 2) # Note the warning
    #> Warning in .run_check_conv(ansseq, check_convergence, tol_conv): Last term is >1.2e-04 times as large as the series,
    #>   suggesting non-convergence. Consider using larger m
    res
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 0.06799312, Error = 2.254503 (one-sided)
    #> Possible range:
    #>  0.0679931219 2.3224957273

    plot(res) # Non-convergence evident

    ## Try larger m
    res <- qfrm(A, B, 2, 2, m = 2000)
    res
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 0.06976012, Error = 4.865235e-07 (one-sided)
    #> Possible range:
    #>  0.0697601234 0.0697606100

    plot(res) # Much better

    #####
    ## VI: Difficult case with highly dispersed eigenvalues
    ##     Error bound unavailable

    ## It is recommended to use a tol_conv stricter than
    ## the default value, which is .Machine$double.eps ^ (1/4)
    res <- qfrm(A, B, 1/2, 1/2, tol_conv = sqrt(.Machine$double.eps))
    #> Warning in .run_check_conv(ansseq, check_convergence, tol_conv): Last term is >1.5e-08 times as large as the series,
    #>   suggesting non-convergence. Consider using larger m
    res
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 0.3954463
    #> Error bound unavailable; recommended to inspect plot() of this object

    plot(res) # Close to convergence, but still growing

    ## Try larger m
    res <- qfrm(A, B, 1/2, 1/2, m = 500, tol_conv = sqrt(.Machine$double.eps))
    res
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 0.3954721
    #> Error bound unavailable; recommended to inspect plot() of this object

    plot(res) # Better

## Checking for numerical convergence

In the last example above, the argument `tol_conv` was used to adjust the tolerance in assessing numerical convergence; by default (`check_convergence = "relative"`), a warning is thrown when the ratio of the last term in the series to the partial sum is larger than this tolerance, whose default is `.Machine$double.eps ^ (1/4)` (= `~1.2e-4`). Using a smaller value for this argument as done above can safeguard against non-convergence, although it is up to the user’s discretion at which point the result is deemed satisfactory.

Other options:

- `check_convergence = "strict_relative"` or `TRUE`: sets `tol_conv = .Machine$double.eps` (unless a smaller `tol_conv` is specified by the user), which is rather strict
- `check_convergence = "absolute"`: check is done in terms of absolute value rather than relative magnitude
- `check_convergence = "none"` or `FALSE`: skips convergence check

## Computational cost

The computational cost for the recursive algorithms increases drastically with the order of evaluation \(m\). An approximate number of recursions required in each function is listed in the tables for `qfrm` and `qfmrm` functions above; the computational memory and time increase at least proportionally to this value. Of course the computational cost is also influenced by the number of variables \(n\).

In extremely large problems, therefore, it may be impossible to obtain a satisfactory numerical convergence within a reasonable amount of time (and memory). A potential practical alternative is a Monte Carlo approach using `rqfr()`.

    ## A large problem
    nv <- 200
    large_A <- diag(c(1000, rep.int(1, nv - 1)))
    large_B <- diag(c(rep.int(1, nv - 1), 1000))
    large_D <- diag((nv:1) ^ 2)

    res <- qfmrm(large_A, large_B, large_D, 1, 1/2, 1/2, m = 500)
    #> Warning in .run_check_conv(ansseq, check_convergence, tol_conv): Last term is >1.2e-04 times as large as the series,
    #>   suggesting non-convergence. Consider using larger m
    res
    #>
    #>  Moment of ratio of quadratic forms
    #>
    #> Moment = 0.02519373
    #> Error bound unavailable; recommended to inspect plot() of this object

    plot(res) # Far to convergence

    ## This problem needs much larger m (>5000)
    ## taking computational time + memory
    ## (though still manageable with a regular machine in this particular case)

    ## Monte Carlo sample
    MCres <- rqfmr(10000, large_A, large_B, large_D, 1, 1/2, 1/2)

    ## Monte Carlo 95% CI of the moment
    mean(MCres) + sd(MCres) / sqrt(10000) * qt(c(0.025, 0.975), 10000 - 1)
    #> [1] 0.02958917 0.03117620

# 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.
Chikuse, Y. (1980) Invariant polynomials with matrix arguments and their applications. In: Gupta, R. P., (ed.), Multivariate Statistical Analysis. Amsterdam: North-Holland. pp. 53–68.
Chikuse, Y. (1987) Methods for constructing top order invariant polynomials. Econometric Theory, 3, 195–207. doi:10.1017/S026646660001029X.
Davis, A. W. (1980) Invariant polynomials with two matrix arguments, extending the zonal polynomials. In: Krishnaiah, P. R., (ed.), Multivariate Analysis—V. Amsterdam: North-Holland. pp. 287–299.
Hansen, T. F. and Houle, D. (2008) Measuring and comparing evolvability and constraint in multivariate characters. Journal of Evolutionary Biology, 21, 1201–1219. doi:10.1111/j.1420-9101.2008.01573.x.
Hillier, G. (2001) The density of a quadratic form in a vector uniformly distributed on the n-sphere. Econometric Theory, 17, 1–28. doi:10.1017/S026646660117101X.
Hillier, G., Kan, R. and Wang, X. (2009) Computationally efficient recursions for top-order invariant polynomials with applications. Econometric Theory, 25, 211–242. 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. doi:10.1017/S0266466613000364.
Mathai, A. M. and Provost, S. B. (1992) Quadratic Forms in Random Variables: Theory and Applications. New York, New York: Marcel Dekker.
Melo, D., Garcia, G., Hubbe, A., Assis, A. P. and Marroig, G. (2016) EvolQG—an R package for evolutionary quantitative genetics [version 3; referees: 2 approved, 1 approved with reservations]. F1000Research, 4, 925. doi:10.12688/f1000research.7082.3.
Muirhead, R. J. (1982) Aspects of Multivariate Statistical Theory. Hoboken, New Jersey: John Wiley & Sons. doi:10.1002/9780470316559.
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.
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. 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. doi:10.1007/s00285-023-01930-8.