New special functions: logarithm of incomplete gamma and beta functions

[dschmitz89]

The incomplete gamma and beta functions arise in a LOT of cumulative density functions. Very accurate implementations of them are quite common and also in Boost. What I have not seen yet are implementations of their logarithms. These would be very helpful to compute accurate logarithms of CDF and SFs. One example is an issue in SciPy where a user wants to compute the logcdf of the Poisson distribution for not even very extreme values: scipy/scipy#8424. Here the logarithm of the incomplete regularized gamma function would come very much in handy.

SciPy includes the logarithm of the normal distribution CDF already but incomplete beta and gamma functions are of course a larger beast to kill.

Is this something you would find worthwile for Boost?

Edit: CC ing the special function wizards @steppi @WarrenWeckesser

[mborland]

The incomplete gamma and beta functions arise in a LOT of cumulative density functions. Very accurate implementations of them are quite common and also in Boost. What I have not seen yet are implementations of their logarithms. These would be very helpful to compute accurate logarithms of CDF and SFs. One example is an issue in SciPy where a user wants to compute the logcdf of the Poisson distribution for not even very extreme values: scipy/scipy#8424. Here the logarithm of the incomplete regularized gamma function would come very much in handy.

SciPy includes the logarithm of the normal distribution CDF already but incomplete beta and gamma functions are of course a larger beast to kill.

Is this something you would find worthwile for Boost?

We certainly wouldn't turn them down. A while back @mdhaber asked about adding the logpdf and logcdfto the distribution so many were done (e.g. https://github.com/boostorg/math/blob/develop/include/boost/math/distributions/laplace.hpp#L235), but I think your proposal would help us both out. Right now the poisson dist simply uses: log(cdf(...))).

We have logsumexp here: https://github.com/boostorg/math/blob/develop/include/boost/math/special_functions/logsumexp.hpp, and continued fractions here: https://github.com/boostorg/math/blob/2a4351eb80a18c9abb5a40b96917df073f0939ae/include/boost/math/tools/simple_continued_fraction.hpp and here: https://github.com/boostorg/math/blob/2a4351eb80a18c9abb5a40b96917df073f0939ae/include/boost/math/tools/centered_continued_fraction.hpp depending on the route of implementation you're looking at.

Let us know what kind of assistance you need.

[jzmaddock]

Internally, the incomplete gamma and beta's are mostly computed as the product of an asymptotic part and a power term (I say mostly, because there are many code paths depending on the arguments). So in principle, calculating the log is relatively straightforward, it's just untangling all those code paths.

I might be tempted to do what I did for 1F1, and pass a scaling term down through all the calculations - the scaling term is an integer (so no calculation error at least in principle), and we then scale the final result by e^N, or of course, take logs if we're going that way.

BTW we already have the logpdf for this one correctly handled.

I also wonder if your original issue was motivated by the desire to get the complement of the CDF? ie he's asking for an accurate logcdf when the cdf is near 1. Why? To pass into expm1 to get the complement maybe? But I digress!

[dschmitz89]

Thanks for your feedback, sounds like you agree that these would be useful. I am personally not qualified to help out as I am not a C++ or special function expert though. Let's see if we can find a champion for this one :).

About the Poisson logpdf: do you list log methods in the distribution docs? For Poisson, I cannot find it and the general distribution page does not list it either. In general if log methods are implemented, even generically without high precision, I would expect to find that info somewhere.

[jzmaddock]

I think @mborland added those, did the docs get updated at the same time?

[mborland]

I think @mborland added those, did the docs get updated at the same time?

Apparently not so I'll get them merged in the next few days. The short of it is all distributions have logpdf and logcdf minimally through generic functions (e.g. log(cdf(...))) and some have one or both specialized.

[mborland]

Thanks for your feedback, sounds like you agree that these would be useful. I am personally not qualified to help out as I am not a C++ or special function expert though. Let's see if we can find a champion for this one :).

About the Poisson logpdf: do you list log methods in the distribution docs? For Poisson, I cannot find it and the general distribution page does not list it either. In general if log methods are implemented, even generically without high precision, I would expect to find that info somewhere.

See: #1176

[JacobHass8]

If this is still open, I wouldn't mind working on it. I'm not exactly sure how to tackle it though. Would taking the logarithm of the continued fraction be sufficient?

Apparently Fortran has an implementation. Could I translate this to C++ and use this?

[mdhaber]

Could I translate this to C++ and use this?

I'll let the Boost Math folks answer this directly, but the Fortran implementation is MIT-licensed, which is quite permissive.

Would taking the logarithm of the continued fraction be sufficient?

Maybe if you use the log of the terms and use the logsumexp trick to do the arithmetic. I wrote a private function for scipy.stats for evaluating generalized continued fractions. If you set log=True, then if your a and b callables evaluate the logarithm of the numerator and denominator coefficients, the result is the logarithm of the continued fraction, with all intermediate arithmetic performed in log-space. The function is vectorized and array API compatible (e.g. use CuPy backend), so you could use this to prototype something in Python and run extensive numerical tests to see how well it works.

[jzmaddock]

@JacobHass8 : take a look at my reply above. The existing code has a few different methods internally depending on what domain we're in, even taking logs we would need to emulate that to maintain accuracy and ensure we do do get convergence (eventually). Most of the methods though have a series or CF part which we can always calculate (doesn't overflow - at least probably not!) and the asymptotic part in regularised_gamma_prefix which computes (z^a)(e^-z)/tgamma(a) for the regularized form and (z^a)(e^-z) in full_igamma_prefix for the non-regularized one. Obviously these can be computed via logs to avoid overflow, but one word of warning: there is a danger, especially in the regularized prefix term, that the logs cancel out and the result is basically garbage. In short there are no free lunches here, which is why I asked what the motivation was...

Ah OK just looked at the original issue, so you want log(Q[k+1, lambda]), and the only time that taking the log of Q[k+1, lambda] fails is when the result is tiny, and in that special case everything goes through this function for integer k+1:

template <class T, class Policy>
BOOST_MATH_GPU_ENABLED inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0)
{
   //
   // Calculates normalised Q when a is an integer:
   //
   BOOST_MATH_STD_USING
   T e = exp(-x);
   T sum = e;
   if(sum != 0)
   {
      T term = sum;
      for(unsigned n = 1; n < a; ++n)
      {
         term /= n;
         term *= x;
         sum += term;
      }
   }
   if(pderivative)
   {
      *pderivative = e * pow(x, a) / boost::math::unchecked_factorial<T>(itrunc(T(a - 1), pol));
   }
   return sum;
}

Which should be trivial to evaluate as logs (ignore the derivative part in your case, and scale the sum by e^-x until the end and then take logs).

I still wonder why anyone would case about probabilities that small in the real world, but no matter!

Let me know if I'm missing anything, but a special case adaptation seems the best option rather than a fully general one here.

[JacobHass8]

Ah OK just looked at the original issue, so you want log(Q[k+1, lambda]), and the only time that taking the log of Q[k+1, lambda] fails is when the result is tiny, and in that special case everything goes through this function for integer k+1:

This is just for integer values of k+1, correct? I think that the idea was to generalize this to any values of $k \in \mathbb{R}$.

Which should be trivial to evaluate as logs (ignore the derivative part in your case, and scale the sum by e^-x until the end and then take logs).

I think for integer k+1 this would be an easy solution.

I still wonder why anyone would case about probabilities that small in the real world, but no matter!

I've actually used these small probabilities to study large deviations.

[JacobHass8]

Maybe if you use the log of the terms and use the logsumexp trick to do the arithmetic. I wrote a private function for scipy.stats for evaluating generalized continued fractions. If you set log=True, then if your a and b callables evaluate the logarithm of the numerator and denominator coefficients, the result is the logarithm of the continued fraction, with all intermediate arithmetic performed in log-space. The function is vectorized and array API compatible (e.g. use CuPy backend), so you could use this to prototype something in Python and run extensive numerical tests to see how well it works.

Using your code, I've been able to get some promising results. Here's an example script I made and the output:

Example

import numpy as np
from scipy.stats._continued_fraction import _continued_fraction as cf
from scipy.special import gammaincc, gamma
import warnings

def gammaincc_series(s, z, size=100):
    # Return series an, bn for incomplete gamma function at (s, z)
    an = np.zeros(size)
    bn = np.zeros(size)
    an[1] = 1
    bn[0] = 0
    bn[1:] = 2 * np.arange(0, size-1) + 1 + z - s
    an[2:] = np.arange(1, size-1) * (s - np.arange(1, size-1))
    return an, bn

s = 10
z = 1000
a, b = gammaincc_series(s, z, 100)

# Stop evaluation of series when an = 0
zero_idx = np.where(a == 0)[0][1]
a = a[:zero_idx + 1]
b = b[:zero_idx + 1]

# Calculate without using log arithmetic
res = cf(a=lambda n: a[n], b=lambda n: b[n], maxiter=len(b) - 1, log=False)
gammaincc_val = z ** s * np.exp(-z) * res.f
with np.errstate(divide='ignore'):
    print('Without using log arithmetic:', np.log(gammaincc_val))
    # Scipy implementation differs from wikipedia by factor of gamma(s)
    print('Current scipy.special implementation:', np.log(gammaincc(s, z) * gamma(s)))

# When log=True, get complex warning and divide by zero in log
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    a = np.log(a)
    b = np.log(b)
    res = cf(a=lambda n: a[n], b=lambda n: b[n], maxiter=len(b) - 1, log=True)
gammaincc_val = res.f + s * np.log(z) - z 
print('Using log arithmetic:', gammaincc_val)

where only the log arithmetic version returns a number:

Without using log arithmetic: -inf
Current scipy.special implementation: -inf
Using log arithmetic: -937.8211708900747

For reference, this matches Mathematica.

However, there are issues when $a_n$ or $b_n$ are negative. When this happens I can't take the log to put into the _continued_fraction function. I'm sure there's an easy fix, but I haven't dug into the code just yet. I'm thinking this might be the way to go.

[mdhaber]

I think _continued_fraction will "just work" with complex numbers. Use a complex dtype and the imaginary component of the log a negative number will be $\pi$ (or some multiple); e.g.

import numpy  as np
np.log(-1 + 0j)  # np.complex128(3.141592653589793j)

[JacobHass8]

I think _continued_fraction will "just work" with complex numbers. Use a complex dtype and the imaginary component of the log a negative number will be π (or some multiple); e.g.

Yeah, that's such a simple solution! Do we still need a C++ implementation of this or would this Python implementation be adequate?

[mdhaber]

We'd want C++. The goal would be for it to make its way into XSF, then into other backends besides SciPy, and then hopefully it could be part of the array API special function extension.

[JacobHass8]

Great, I'll work on translating your code to C++. Thanks for the help!