# sadists package/ Sum of (non-central) chi-squares to a power

I try use "sadists" package in R to compute quantiles and probabilities on sum of non-central chi-squares distribution, but there are some issues in this package. I give an example:

wts = c(0.5000, -103.5088, -853.8842)

df = c(1, 10, 8)

ncp = c(0.0000000, 3.6615283, 0.9199399)

pow = c(1, 1, 1)

psumchisqpow(-2000, wts, df, ncp, pow)

 1

psumchisqpow(1000, wts, df, ncp, pow)

 0.9955206

So Prob(X <= -2000) = 1 and Prob(X <= 1000) = 0.995 !! How this behavior can occur? That means there are negative values in the PDF in the range [-2000, 1000].

The second issue is when I add lower.tail = FALSE, psumchisqpow gives NaN:

psumchisqpow(1000, wts, df, ncp, pow, lower.tail = FALSE)

 NaN

Warning message: In sqrt(raw.cumulants) : NaNs produced

How can I deal with these issues (I miss something?)? Is there an alternative way to estimate quantiles and probabilities on weighted sum of non-central chi squares distribution? Thanks

• Sounds like you have a bug report for the authors of the package. Jun 8 '15 at 16:43

You may want to look into the approximations underlying the function you're using. It seems within your control to resolve using the order.max argument, given the below examples.

> # libraries and options ---------------------------------------------------
>
> library(pbapply)
>
> # setup -------------------------------------------------------------------
>
> degrees_of_freedoms <- c(1, 10, 8)
> noncentrality_parameters <- c(0.0000000, 3.6615283, 0.9199399)
> weights <- c(0.5000, -103.5088, -853.8842)
>
> # approximation through simulation ----------------------------------------
>
> sim_sum_noncentral_chi_squares <- function(
+   degrees_of_freedoms,
+   noncentrality_parameters,
+   weights
+ ) {
+   variates <- sapply(1:length(degrees_of_freedoms), function(index) {
+     rchisq(1, degrees_of_freedoms[index], noncentrality_parameters[index])
+   })
+   return(sum(weights * variates))
+ }
>
> random_draws <- pbreplicate(
+   1000000,
+   sim_sum_noncentral_chi_squares(
+     degrees_of_freedoms,
+     noncentrality_parameters,
+     weights
+   )
+ )
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=21s
>
> ecdf(random_draws)(c(-2000, 1000))
 0.99834 1.00000
>
> # approximation with sadists ----------------------------------------------
>
> random_draws_2 <- rsumchisqpow(
+   1000000,
+   weights,
+   degrees_of_freedoms,
+   noncentrality_parameters
+ )
>
> ecdf(random_draws_2)(c(-2000, 1000))
 0.998275 1.000000
>
> # approximation with psumchisqpow -----------------------------------------
>
> psumchisqpow(
+   c(-2000, 1000),
+   weights,
+   degrees_of_freedoms,
+   noncentrality_parameters,
+   order.max = 9
+ )