qfpm: Moment of (product of) quadratic forms in normal variables

Description

Functions to obtain (compound) moments of a product of quadratic forms in normal variables, i.e., $ \mathrm{E} \left( (\mathbf{x^\mathit{T} A x})^p (\mathbf{x^\mathit{T} B x})^q (\mathbf{x^\mathit{T} D x})^r \right) $, where $\mathbf{x} \sim N_n(\bm{\mu}, \mathbf{\Sigma})$.

qfm_Ap_int() is for $q = r = 0$ (simple moment)

qfpm_ABpq_int() is for $r = 0$

qfpm_ABDpqr_int() is for the product of all three powers

Usage

qfm_Ap_int(
  A,
  p = 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
)
qfpm_ABpq_int(
  A,
  B,
  p = 1,
  q = 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
)
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
)

Value

A qfpm object which has the same elements as those returned by the qfrm functions. Use $statistic to access the value of the moment.

Arguments

A, B, D: Argument matrices. Should be square. Will be automatically symmetrized.
p, q, r: Exponents for $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{D}$, respectively. By default, these are set to the same value. If unsure, specify all explicitly.
mu: Mean vector $\bm{\mu}$ for $\mathbf{x}$
Sigma: Covariance matrix $\mathbf{\Sigma}$ for $\mathbf{x}$
use_cpp: Logical to specify whether the calculation is done with C++ functions via Rcpp. TRUE by default.
cpp_method: Variable type used in C++ calculations. In these functions this is ignored.
tol_zero: Tolerance against which numerical zero is determined. Used to determine, e.g., whether mu is a zero vector, A or B equals the identity matrix, etc.
tol_sing: Tolerance against which matrix singularity and rank are determined. The eigenvalues smaller than this are considered zero.

Details

These functions implement the super-short recursion algorithms described in Hillier et al. (2014: sec. 3.1--3.2 and 4). At present, only positive integers are accepted as exponents (negative exponents yield ratios, of course). All these yield exact results.

Examples

Run this code

## Some symmetric matrices and parameters
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 where x ~ N(0, I)
qfm_Ap_int(A, 2)

## This is the same but obviously less efficient
qfpm_ABpq_int(A, p = 2, q = 0)

## Expectation of (x^T A x) (x^T B x) (x^T D x) where x ~ N(0, I)
qfpm_ABDpqr_int(A, B, D, 1, 1, 1)

## Expectation of (x^T A x) (x^T B x) (x^T D x) where x ~ N(mu, Sigma)
qfpm_ABDpqr_int(A, B, D, 1, 1, 1, mu = mu, Sigma = Sigma)

## Expectations of (x^T x)^2 where x ~ N(0, I) and x ~ N(mu, I)
## i.e., roundabout way to obtain moments of
## central and noncentral chi-square variables
qfm_Ap_int(diag(nv), 2)
qfm_Ap_int(diag(nv), 2, mu = mu)