rqfr: Monte Carlo sampling of ratio/product of quadratic forms

Description

rqfr(), rqfmr(), and rqfp() calculate a random sample of a simple ratio, multiple ratio (of special form), and product, respectively, of quadratic forms in normal variables of specified mean and covariance (standard multivariate normal by default). These functions are primarily for empirical verification of the analytic results provided in this package.

Usage

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)

Value

Numeric vector of length nit.

Arguments

nit: Number of iteration or sample size. Should be an integer-alike of length 1.
A, B, D: Argument matrices (see “Details”). Assumed to be square matrices of the same order. When missing, set to the identity matrix. At least one of these must be specified.
p, q, r: Exponents for A, B, D, respectively (see “Details”). Assumed to be numeric of length 1 each. See “Details” for default values.
mu: Mean vector \(\bm{\mu}\) for \(\mathbf{x}\). Default zero vector.
Sigma: Covariance matrix \(\mathbf{\Sigma}\) for \(\mathbf{x}\). Default identity matrix. mu and Sigma are assumed to be of the same order as the argument matrices.
use_cpp: Logical to specify whether an C++ version is called or not. TRUE by default.

Details

These functions generate a random sample of \( \frac{(\mathbf{x^\mathit{T} A x})^p}{(\mathbf{x^\mathit{T} B x})^q} \) (rqfr()), \( \frac{(\mathbf{x^\mathit{T} A x})^p} {(\mathbf{x^\mathit{T} B x})^q (\mathbf{x^\mathit{T} Dx})^r} \) (rqfmr()), and \( (\mathbf{x^\mathit{T} A x})^p (\mathbf{x^\mathit{T} B x})^q (\mathbf{x^\mathit{T} D x})^r \) (rqfp()), where \(\mathbf{x} \sim N_n(\bm{\mu}, \mathbf{\Sigma}) \). (Internally, rqfr() and rqfmr() just call rqfp() with negative exponents.)

When only one of p and q is provided in rqfr(), the other (missing) one is set to the same value.

In rqfmr(), q and r are set to p/2 when both missing, and set to the same value when only one is missing. When p is missing, this is set to be q + r. If unsure, specify all these explicitly.

In rqfp(), p, q and r are 1 by default, provided that the corresponding argument matrices are given. If both an argument matrix and its exponent (e.g., D and r) are missing, the exponent is set to 0 so that the factor be unity.

Examples

Run this code

p <- 4
A <- diag(1:p)
B <- diag(p:1)
D <- diag(sqrt(1:p))

## By default B = I, p = q = 1;
## i.e., (x^T A x) / (x^T x), x ~ N(0, I)
rqfr(5, A)

## (x^T A x) / ((x^T B x)(x^T D x))^(1/2), x ~ N(0, I)
rqfmr(5, A, B, D, 1, 1/2, 1/2)

## (x^T A x), x ~ N(0, I)
rqfp(5, A)

## (x^T A x) (x^T B x), x ~ N(0, I)
rqfp(5, A, B)

## (x^T A x) (x^T B x) (x^T D x), x ~ N(0, I)
rqfp(5, A, B, D)

## Example with non-standard normal
mu <- 1:p / p
Sigma <- matrix(0.5, p, p)
diag(Sigma) <- 1
rqfr(5, A, mu = 1:p / p, Sigma = Sigma)

## Compare Monte Carlo sample and analytic expression
set.seed(3)
mcres <- rqfr(1000, A, p = 2)
mean(mcres)
(anres <- qfrm(A, p = 2))
stats::t.test(mcres, mu = anres$statistic)