Package 'CompQuadForm'

Davies method

Description

Distribution function (survival function in fact) of quadratic forms in normal variables using Davies's method.

Usage

davies(q, lambda, h = rep(1, length(lambda)), delta = rep(0,
         length(lambda)), sigma = 0, lim = 10000, acc = 0.0001)

Arguments

q	value point at which distribution function is to be evaluated
lambda	the weights $\lambda_1, \lambda_2, ..., \lambda_n$, i.e. distinct non-zero characteristic roots of $A\Sigma$
h	respective orders of multiplicity $n_j$ of the $\lambda$s
delta	non-centrality parameters $\delta_j^2$ (should be positive)
sigma	coefficient $\sigma$ of the standard Gaussian
lim	maximum number of integration terms. Realistic values for 'lim' range from 1,000 if the procedure is to be called repeatedly up to 50,000 if it is to be called only occasionally
acc	error bound. Suitable values for 'acc' range from 0.001 to 0.00005 which should be adequate for most statistical purposes.

Details

Computes $P[Q>q]$ where $Q = \sum_{j=1}^r\lambda_jX_j+\sigma X_0$ where $X_j$ are independent random variables having a non-central $chi^2$ distribution with $n_j$ degrees of freedom and non-centrality parameter $delta_j^2$ for $j=1,...,r$ and $X_0$ having a standard Gaussian distribution.

Value

trace	vector, indicating performance of procedure, with the following components: 1: absolute value sum, 2: total number of integration terms, 3: number of integrations, 4: integration interval in main integration, 5: truncation point in initial integration, 6: standard deviation of convergence factor term, 7: number of cycles to locate integration parameters
ifault	fault indicator: 0: no error, 1: requested accuracy could not be obtained, 2: round-off error possibly significant, 3: invalid parameters, 4: unable to locate integration parameters
Qq	$P[Q>q]$

Author(s)

Pierre Lafaye de Micheaux ([email protected]) and Pierre Duchesne ([email protected])

References

P. Duchesne, P. Lafaye de Micheaux, Computing the distribution of quadratic forms: Further comparisons between the Liu-Tang-Zhang approximation and exact methods, Computational Statistics and Data Analysis, Volume 54, (2010), 858-862

Davies R.B., Algorithm AS 155: The Distribution of a Linear Combination of chi-2 Random Variables, Journal of the Royal Statistical Society. Series C (Applied Statistics), 29(3), p. 323-333, (1980)

Examples

# Some results from Table 3, p.327, Davies (1980)

 round(1 - davies(1, c(6, 3, 1), c(1, 1, 1))$Qq, 4)
 round(1 - davies(7, c(6, 3, 1), c(1, 1, 1))$Qq, 4)
 round(1 - davies(20, c(6, 3, 1), c(1, 1, 1))$Qq, 4)
 
 round(1 - davies(2, c(6, 3, 1), c(2, 2, 2))$Qq, 4)
 round(1 - davies(20, c(6, 3, 1), c(2, 2, 2))$Qq, 4)
 round(1 - davies(60, c(6, 3, 1), c(2, 2, 2))$Qq, 4)
 
 round(1 - davies(10, c(6, 3, 1), c(6, 4, 2))$Qq, 4)
 round(1 - davies(50, c(6, 3, 1), c(6, 4, 2))$Qq, 4)
 round(1 - davies(120, c(6, 3, 1), c(6, 4, 2))$Qq, 4)

 round(1 - davies(20, c(7, 3), c(6, 2), c(6, 2))$Qq, 4)
 round(1 - davies(100, c(7, 3), c(6, 2), c(6, 2))$Qq, 4)
 round(1 - davies(200, c(7, 3), c(6, 2), c(6, 2))$Qq, 4)

 round(1 - davies(10, c(7, 3), c(1, 1), c(6, 2))$Qq, 4)
 round(1 - davies(60, c(7, 3), c(1, 1), c(6, 2))$Qq, 4)
 round(1 - davies(150, c(7, 3), c(1, 1), c(6, 2))$Qq, 4)

 round(1 - davies(70, c(7, 3, 7, 3), c(6, 2, 1, 1), c(6, 2, 6, 2))$Qq, 4)
 round(1 - davies(160, c(7, 3, 7, 3), c(6, 2, 1, 1), c(6, 2, 6, 2))$Qq, 4)
 round(1 - davies(260, c(7, 3, 7, 3), c(6, 2, 1, 1), c(6, 2, 6, 2))$Qq, 4)

 round(1 - davies(-40, c(7, 3, -7, -3), c(6, 2, 1, 1), c(6, 2, 6,
 2))$Qq, 4)
 round(1 - davies(40, c(7, 3, -7, -3), c(6, 2, 1, 1), c(6, 2, 6, 2))$Qq,
 4)
 round(1 - davies(140, c(7, 3, -7, -3), c(6, 2, 1, 1), c(6, 2, 6,
 2))$Qq, 4)

# You should sometimes play with the 'lim' parameter:
davies(0.00001,lambda=0.2)
imhof(0.00001,lambda=0.2)$Qq
davies(0.00001,lambda=0.2, lim=20000)

---

Ruben/Farebrother method

Description

Distribution function (survival function in fact) of quadratic forms in normal variables using Farebrother's algorithm.

Usage

farebrother(q, lambda, h = rep(1, length(lambda)),
            delta = rep(0, length(lambda)), maxit = 100000,
            eps = 10^(-10), mode = 1)

Arguments

q	value point at which distribution function is to be evaluated
lambda	the weights $\lambda_1, \lambda_2, ..., \lambda_n$, i.e. the distinct non-zero characteristic roots of $A\Sigma$
h	vector of the respective orders of multiplicity $m_i$ of the $\lambda$s
delta	the non-centrality parameters $\delta_i$ (should be positive)
maxit	the maximum number of term K in equation below
eps	the desired level of accuracy
mode	if 'mode' > 0 then $\beta=mode*\lambda_{min}$ otherwise $\beta=\beta_B=2/(1/\lambda_{min}+1/\lambda_{max})$

Details

Computes P[Q>q] where $Q=\sum_{j=1}^n\lambda_j\chi^2(m_j,\delta_j^2)$. P[Q<q] is approximated by $\sum_k=0^{K-1} a_k P[\chi^2(m+2k)<q/\beta]$ where $m=\sum_{j=1}^n m_j$ and $\beta$ is an arbitrary constant (as given by argument mode).

Value

dnsty	the density of the linear form
ifault	the fault indicator. -i: one or more of the constraints $\lambda_i>0$

, $m_i>0$ and $\delta_i^2\geq0$ is not satisfied. 1: non-fatal underflow of $a_0$. 2: one or more of the constraints $n>0$, $q>0$, $maxit>0$ and $eps>0$ is not satisfied. 3: the current estimate of the probability is greater than 2. 4: the required accuracy could not be obtained in 'maxit' iterations. 5: the value returned by the procedure does not satisfy $0\leq RUBEN\leq 1$. 6: 'dnsty' is negative. 9: faults 4 and 5. 10: faults 4 and 6. 0: otherwise.

Qq

$P[Q>q]$

Author(s)

Pierre Lafaye de Micheaux ([email protected]) and Pierre Duchesne ([email protected])

References

P. Duchesne, P. Lafaye de Micheaux, Computing the distribution of quadratic forms: Further comparisons between the Liu-Tang-Zhang approximation and exact methods, Computational Statistics and Data Analysis, Volume 54, (2010), 858-862

Farebrother R.W., Algorithm AS 204: The distribution of a Positive Linear Combination of chi-squared random variables, Journal of the Royal Statistical Society, Series C (applied Statistics), Vol. 33, No. 3 (1984), p. 332-339

Examples

# Some results from Table 3, p.327, Davies (1980)

 1 - farebrother(1, c(6, 3, 1), c(1, 1, 1), c(0, 0, 0))$Qq

---

Imhof method.

Description

Distribution function (survival function in fact) of quadratic forms in normal variables using Imhof's method.

Usage

imhof(q, lambda, h = rep(1, length(lambda)),
      delta = rep(0, length(lambda)),
      epsabs = 10^(-6), epsrel = 10^(-6), limit = 10000)

Arguments

q	value point at which the survival function is to be evaluated
lambda	distinct non-zero characteristic roots of $A\Sigma$
h	respective orders of multiplicity of the $\lambda$s
delta	non-centrality parameters (should be positive)
epsabs	absolute accuracy requested
epsrel	relative accuracy requested
limit	determines the maximum number of subintervals in the partition of the given integration interval

Details

Let $\boldsymbol{X}=(X_1,\ldots,X_n)'$ be a column random vector which follows a multidimensional normal law with mean vector $\boldsymbol{0}$ and non-singular covariance matrix $\boldsymbol{\Sigma}$. Let $\boldsymbol{\mu}=(\mu_1,\ldots,\mu_n)'$ be a constant vector, and consider the quadratic form

$Q=(\boldsymbol{x}+\boldsymbol{\mu})'\boldsymbol{A}(\boldsymbol{x}+\boldsymbol{\mu})=\sum_{r=1}^m\lambda_r\chi^2_{h_r;\delta_r}.$

The function imhof computes $P[Q>q]$.

The $\lambda_r$'s are the distinct non-zero characteristic roots of $A\Sigma$, the $h_r$'s their respective orders of multiplicity, the $\delta_r$'s are certain linear combinations of $\mu_1,\ldots,\mu_n$ and the $\chi^2_{h_r;\delta_r}$ are independent $\chi^2$-variables with $h_r$ degrees of freedom and non-centrality parameter $\delta_r$. The variable $\chi^2_{h,\delta}$ is defined here by the relation $\chi^2_{h,\delta}=(X_1 + \delta)^2+\sum_{i=2}^hX_i^2$, where $X_1,\ldots,X_h$ are independent unit normal deviates.

Value

Qq	$P[Q>q]$
abserr	estimate of the modulus of the absolute error, which should equal or exceed abs(i - result)

Author(s)

Pierre Lafaye de Micheaux ([email protected]) and Pierre Duchesne ([email protected])

References

P. Duchesne, P. Lafaye de Micheaux, Computing the distribution of quadratic forms: Further comparisons between the Liu-Tang-Zhang approximation and exact methods, Computational Statistics and Data Analysis, Volume 54, (2010), 858-862

J. P. Imhof, Computing the Distribution of Quadratic Forms in Normal Variables, Biometrika, Volume 48, Issue 3/4 (Dec., 1961), 419-426

Examples

# Some results from Table 1, p.424, Imhof (1961)

# Q1 with x = 2
round(imhof(2, c(0.6, 0.3, 0.1))$Qq, 4)

# Q2 with x = 6
round(imhof(6, c(0.6, 0.3, 0.1), c(2, 2, 2))$Qq, 4)

# Q6 with x = 15
round(imhof(15, c(0.7, 0.3), c(1, 1), c(6, 2))$Qq, 4)

---

Liu's method

Description

Distribution function (survival function in fact) of quadratic forms in normal variables using Liu et al.'s method.

Usage

liu(q, lambda, h = rep(1, length(lambda)),
    delta = rep(0, length(lambda)))

Arguments

q	value point at which the survival function is to be evaluated
lambda	distinct non-zero characteristic roots of $A\Sigma$, i.e. the $\lambda_i$'s
h	respective orders of multiplicity $h_i$'s of the $\lambda$'s
delta	non-centrality parameters $\delta_i$'s (should be positive)

Details

New chi-square approximation to the distribution of non-negative definite quadratic forms in non-central normal variables.

Computes $P[Q>q]$ where $Q=\sum_{j=1}^n\lambda_j\chi^2(h_j,\delta_j)$.

This method does not work as good as the Imhof's method. Thus Imhof's method should be recommended.

Value

Qq

$P[Q>q]$

Author(s)

Pierre Lafaye de Micheaux ([email protected]) and Pierre Duchesne ([email protected])

References

P. Duchesne, P. Lafaye de Micheaux, Computing the distribution of quadratic forms: Further comparisons between the Liu-Tang-Zhang approximation and exact methods, Computational Statistics and Data Analysis, Volume 54, (2010), 858-862

H. Liu, Y. Tang, H.H. Zhang, A new chi-square approximation to the distribution of non-negative definite quadratic forms in non-central normal variables, Computational Statistics and Data Analysis, Volume 53, (2009), 853-856

Examples

# Some results from Liu et al. (2009)
# Q1 from Liu et al.
round(liu(2, c(0.5, 0.4, 0.1), c(1, 2, 1), c(1, 0.6, 0.8)), 6)
round(liu(6, c(0.5, 0.4, 0.1), c(1, 2, 1), c(1, 0.6, 0.8)), 6)
round(liu(8, c(0.5, 0.4, 0.1), c(1, 2, 1), c(1, 0.6, 0.8)), 6)

# Q2 from Liu et al.
round(liu(1, c(0.7, 0.3), c(1, 1), c(6, 2)), 6)
round(liu(6, c(0.7, 0.3), c(1, 1), c(6, 2)), 6)
round(liu(15, c(0.7, 0.3), c(1, 1), c(6, 2)), 6)

# Q3 from Liu et al.
round(liu(2, c(0.995, 0.005), c(1, 2), c(1, 1)), 6)
round(liu(8, c(0.995, 0.005), c(1, 2), c(1, 1)), 6)
round(liu(12, c(0.995, 0.005), c(1, 2), c(1, 1)), 6)

# Q4 from Liu et al.
round(liu(3.5, c(0.35, 0.15, 0.35, 0.15), c(1, 1, 6, 2), c(6, 2, 6, 2)),
6)
round(liu(8, c(0.35, 0.15, 0.35, 0.15), c(1, 1, 6, 2), c(6, 2, 6, 2)), 6)
round(liu(13, c(0.35, 0.15, 0.35, 0.15), c(1, 1, 6, 2), c(6, 2, 6, 2)), 6)

---