Calculating miniscule numbers for chi-squared distribution -- numerical precision

5
votes

I am using the pchisq function in R to calculate the cumulative distribution function for the chi-squared distribution. I would like to calculate very small values, such that 1-pchisq(...) can have a value smaller than 2.2e-16 (which is the numerical precision limit for R's numeric format). Right now, these very small values simply become 0.

Things I've tried:

  • Setting the digits display option to 22 (max)

  • Using the Rmpfr package for increased precision, but that number format doesn't work with the pchisq function

  • Breaking the CDF function into its component gamma functions, but that results in similar precision limits. Any ideas on how I can calculate what I want?

Background: I'm using Fisher's method to combine a bunch of p-values. Yes, I know these p-values are minuscule, but it is actually useful for what I'm analyzing.

1
rstatisticsprecisionchi-squared
this is probably a wrong place to ask about stats. – Mox
could you compute log(1-pchisq) ? – Severin Pappadeux

1 Answers

4
votes

A couple of things.

  • 2.2e-16 is not the lower precision limit of values in R; it's just the way that R prints very small p-values by default, using format.pval:
format.pval(1e-20)
## [1] "< 2.22e-16"
  • values smaller than approximately 1e-320 do round down to zero:
1e-330
## [1] 0

@SeverinPappadeux's suggestion is exactly right:

pchisq(121231,1,lower.tail=FALSE,log.p=TRUE)
## [1] -60621.58

This is equivalent to 10^(-26327):

-60621.58/log(10)
## -26327.62

Check for a less extreme value:

log10(pchisq(100,1,lower.tail=FALSE) )
## [1] -22.81702
pchisq(100,1,lower.tail=FALSE,log.p=TRUE)/log(10)
## [1] -22.81702

Furthermore, log(p) is exactly what you need to use for Fisher's method.