1
$\begingroup$

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

$\endgroup$
1
  • 1
    $\begingroup$ Sounds like you have a bug report for the authors of the package. $\endgroup$
    – Andrew M
    Jun 8 '15 at 16:43
0
$\begingroup$

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 
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.