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.

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

Moment existence conditions

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

• B (and D) must be nonnegative definite (i.e., has no negative eigenvalues)
• when p is non-positive-integer, 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 B is positive definite, the moment exists if and only if $\frac{n}{2} + p > q$.
2. When B is positive semidefinite with rank l, consider a matrix of its eigenvectors P = (P₁, P₂), with P₁ and P₂ being n × l and n × (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 B is a subspace of that of D (or vice versa), the existence problem reduces to that of the simple ratio, and the above conditions apply by inserting the rank of B and q + r into l and q above (Watanabe, 2023). However, when both B and 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 Σ, it is always possible to transform the variables and argument matrices by using x_new = K⁻¹x, A_new = K^TAK, etc., where K is an n × n matrix that satisfies KK^T = Σ (Mathai and Provost, 1992, chap. 3). When Σ is singular, certain conditions need to be met by the argument matrices, Σ, and μ 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)) and for a positive-integral p, this can be simplified into (Bao and Kan, 2013, eq. (13), correcting a typo):

For the multiple ratio, the expression for an arbitrary p is (after Smith, 1989, eq. (5.2)) and for a positive-integral p (Watanabe, 2023, eq. (B28)),

In these expressions, β_⋅ are arbitrary scaling constants that satisfy $0 < \beta < 2 / \qfrlmax$, with $\qfrlmax$ being the largest eigenvalue of the argument matrix. In addition, h_k1, …ks and h̃_{k1; k2…ks} are the coefficients of t₁^k₁…t_s^k_s in the following power series expansions: (These notations follow Bao and Kan (2013); Hillier et al. (2014) used the symbol h for h̃ here.)

A common special case is when μ = 0_n, in which case the coefficients d in the following expansion is used:

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): 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 h̃ where A₁ 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 B = I_n. In this case, (2) simplifies into (Hillier et al., 2014, theorem 4) where ₁F₁(⋅; ⋅; ⋅) is the confluent hypergeometric function, and a_p, k are the coefficients of t^p in 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 μ = 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 I_n and μ = 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 ≤ 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, 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 B is positive definite. By denoting the expression (2) as M(A, B, p, q) and the partial sum of the same up to j = m as M̂_m(A, B, p, q), the error bound is (Hillier et al., 2014, theorem 7) where A⁺ is a symmetric matrix constructed from the eigenvectors and “positivized” eigenvalues of 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 d̃ and ĥ are coefficients in the following generating functions:

For the multiple ratio, a similar error bound is available when D = I_n, in addition to p being a positive integer. By denoting the expression (4) as M(A, B, D, p, q, r), and the partial sum of (4) up to j + k = m as M̂_m(A, B, D, p, q, r), the error bound is 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, h̃, and ĥ 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, where with the initial conditions h_0, 0, 0 = h̃_0, 0, 0 = ĥ_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 μ = 0_n (and hence omitting g altogether) in that for h, and those for the two-argument versions are by omitting any terms with A₃ 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): where with λ_i being the ith eigenvalue of A, δ_i = (h_i^Tμ)², h_i being the eigenvector of A corresponding to λ_i, and the initial conditions are d̃₀ = 1 and u_0, i = v_0, i = 0 for all i.

The recursion for a_p, k(A, μ) is (Hillier et al., 2014, eqs. (31), (32), noting that w should be indexed by k as well): where 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 A₂ or j − 1 as an index from that for the two-argument version of d above.

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}}$.

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

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:

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

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.

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)

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 (x^TAx)^p(x^TBx)^q(x^TDx)^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.

rqfr(nit = 1000L, A, B, p = 1, q = p, mu, Sigma, use_cpp = TRUE)
rqfmr(nit = 1000L, A, B, D, p = 1, q = p/2, r = q, mu, Sigma, use_cpp = TRUE)
rqfp(nit = 1000L, A, B, D, p = 1, q = 1, r = 1, mu, Sigma, use_cpp = TRUE)

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:

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

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.421011e-20 4.065758e-20 2.964615e-20 1.948176e-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 qfrm_ApBq_int(A = A, B = B, p = p, q = q, m = m, mu = mu, tol_zero = tol_zero, : 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 qfrm_ApBq_npi(A = A, B = B, p = p, q = q, m = m, mu = mu, tol_zero = tol_zero, : 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 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 qfmrm_ApBDqr_int(A = A, B = B, D = D, p = p, q = q, r = r, m = m, : 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