pymare.stats.log_chi2_sf

log_chi2_sf(q, df)[source]

Natural logarithm of the chi-squared upper tail, without underflowing.

Parameters:
Returns:

log(P(chi^2_df > q)).

Return type:

numpy.ndarray

Notes

scipy.stats.chi2.logsf computes the survival function first and takes its logarithm afterwards, so it returns -inf for every q whose tail falls below the smallest positive double – which Cochran’s Q reaches on a few hundred heterogeneous observations, well inside the range of an ordinary meta-analysis.

The tail is the regularized upper incomplete gamma function Q(a, x) at a = df/2, x = q/2, which factors as

Q(a, x) = \frac{e^{-x} x^{a}}{\Gamma(a)} \cdot F(a, x),

where F is a continued fraction of order one, evaluated here by the modified Lentz algorithm. Everything that underflows lives in the prefactor, and the prefactor is an exponential, so evaluating the fraction on its own and adding -x + a log x - lnGamma(a) to its logarithm gives the answer with no intermediate that can flush to zero.

The fraction is used only where SciPy has already failed, and SciPy is used everywhere else. That is not a preference for one over the other but the split that makes both fast: the fraction’s convergence slows as x approaches a from above, needing hundreds of iterations at x = a + 1 with a large a, and that is precisely the region where the tail is around one half, nothing underflows and SciPy is exact. Where the tail is small enough to have underflowed, x is far enough beyond a that the fraction converges in a handful of iterations.