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
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.
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)q |
is the vector of quantile values at which to evaluate. |
lb |
contains lb[i], the weight for deviate |
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 |
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.
Simon N. Wood simon.wood@r-project.org
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
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)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.