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] 1

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

[1] 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)

[1] NaN
Warning message:
In sqrt(raw.cumulants[2]) : 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
 A: 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)
> library(sadists)
> 
> # 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))
[1] 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))
[1] 0.998275 1.000000
> 
> # approximation with psumchisqpow -----------------------------------------
> 
> psumchisqpow(
+   c(-2000, 1000),
+   weights,
+   degrees_of_freedoms,
+   noncentrality_parameters,
+   order.max = 9
+ )
Loading required package: polynom
[1] 0.9987426 1.0000000
Warning message:
package ‘polynom’ was built under R version 3.5.2 

