# Evaluate the c.d.f. of a weighted sum of chi-squared deviates 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 [REMOVE_ME]$$ P(q< \sum_{i=1}^r \lambda_i X_i + \sigma_z Z)P(q < sum_j lb[j] X_j + sigz Z) [REMOVE_ME_2]$$ 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. ## 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)P(q < sum_j lb[j] X_j + sigz 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. ```r 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 $lb[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`. ## 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 ## Author(s) Simon N. Wood simon.wood@r-project.org ## Examples ```r 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) ```

psum.chisq function