rqfr {qfratio}

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)

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.

Value

Numeric vector of length nit.

See Also

qfrm and qfpm for analytic moments

dqfr for analytic distribution-related functions for simple ratios

Examples

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)