Type: Package
Title: Distribution Function of Quadratic Forms in Normal Variables
Version: 1.4.3
Date: 2017-04-10
Author: P. Lafaye de Micheaux
Maintainer: P. Lafaye de Micheaux <lafaye@unsw.edu.au>
Description: Computes the distribution function of quadratic forms in normal variables using Imhof's method, Davies's algorithm, Farebrother's algorithm or Liu et al.'s algorithm.
License: GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
LazyLoad: yes
Packaged: 2017-04-10 02:26:43 UTC; lafaye
NeedsCompilation: yes
Repository: CRAN
Date/Publication: 2017-04-12 14:28:23 UTC

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 \lambdas

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 (lafaye@dms.umontreal.ca) and Pierre Duchesne (duchesne@dms.umontreal.ca)

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 \lambdas

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 (lafaye@dms.umontreal.ca) and Pierre Duchesne (duchesne@dms.umontreal.ca)

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 \lambdas

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 (lafaye@dms.umontreal.ca) and Pierre Duchesne (duchesne@dms.umontreal.ca)

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 (lafaye@dms.umontreal.ca) and Pierre Duchesne (duchesne@dms.umontreal.ca)

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)