psum.chisq: Evaluate the c.d.f. of a weighted sum of chi-squared deviates

In mgcv: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation

Evaluate the c.d.f. of a weighted sum of chi-squared deviates

Description

Evaluates the c.d.f. of a weighted sum of chi-squared random variables by the method of Davies (1973, 1980). That is it computes

P(q< \sum_{i=1}^r \lambda_i X_i + \sigma_z Z)

where X_j is a chi-squared random variable with df[j] (integer) degrees of freedom and non-centrality parameter nc[j], while Z is a standard normal deviate.

Usage

psum.chisq(q,lb,df=rep(1,length(lb)),nc=rep(0,length(lb)),sigz=0,
           lower.tail=FALSE,tol=2e-5,nlim=100000,trace=FALSE)

Arguments

q	is the vector of quantile values at which to evaluate.
lb	contains \lambda_i, the weight for deviate i. Weights can be positive and/or negative.
df	is the integer vector of chi-squared degrees of freedom.
nc	is the vector of non-centrality parameters for the chi-squared deviates.
sigz	is the multiplier for the standard normal deviate. Non- positive to exclude this term.
lower.tail	indicates whether lower of upper tail probabilities are required.
tol	is the numerical tolerance to work to.
nlim	is the maximum number of integration steps to allow
trace	can be set to TRUE to return some trace information and a fault code as attributes.

Details

This calls a C translation of the original Algol60 code from Davies (1980), which numerically inverts the characteristic function of the distribution (see Davies, 1973). Some modifications have been made to remove goto statements and global variables, to use a slightly more efficient sorting of lb and to use R functions for log(1+x). In addition the integral and associated error are accumulated in single terms, rather than each being split into 2, since only their sums are ever used. If q is a vector then psum.chisq calls the algorithm separately for each q[i].

If the Davies algorithm returns an error then an attempt will be made to use the approximation of Liu et al (2009) and a warning will be issued. If that is not possible then an NA is returned. A warning will also be issued if the algorithm detects that round off errors may be significant.

If trace is set to TRUE then the result will have two attributes. "ifault" is 0 for no problem, 1 if the desired accuracy can not be obtained, 2 if round-off error may be significant, 3 is invalid parameters have been supplied or 4 if integration parameters can not be located. "trace" is a 7 element vector: 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. sd of convergence factor term; 7. number of cycles to locate integration parameters. See Davies (1980) for more details. Note that for vector q these attributes relate to the final element of q.

Author(s)

Simon N. Wood [e-mail]

References

Davies, R. B. (1973). Numerical inversion of a characteristic function. Biometrika, 60(2), 415-417.

Davies, R. B. (1980) Algorithm AS 155: The Distribution of a Linear Combination of Chi-squared Random Variables. J. R. Statist. Soc. C 29, 323-333

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

Examples

  require(mgcv)
  lb <- c(4.1,1.2,1e-3,-1) ## weights
  df <- c(2,1,1,1) ## degrees of freedom
  nc <- c(1,1.5,4,1) ## non-centrality parameter
  q <- c(1,6,20) ## quantiles to evaluate

  psum.chisq(q,lb,df,nc)

  ## same by simulation...
  
  psc.sim <- function(q,lb,df=lb*0+1,nc=df*0,ns=10000) {
    r <- length(lb);p <- q
    X <- rowSums(rep(lb,each=ns) *
         matrix(rchisq(r*ns,rep(df,each=ns),rep(nc,each=ns)),ns,r))
    apply(matrix(q),1,function(q) mean(X>q))	 
  } ## psc.sim
  
  psum.chisq(q,lb,df,nc)
  psc.sim(q,lb,df,nc,100000)

mgcv documentation built on May 29, 2024, 4:34 a.m.