qfpm {qfratio}

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
)

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.

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.

See Also

qfrm and qfmrm for moments of ratios

Examples

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