beta(-n,n) should return ±1/n, not NaN

[stevengj]

In #14165, @CarloLucibello noticed a bug in our beta function: beta(-n,n) is returning NaN for integers n, whereas the correct answer is (-1)^n / n, I think.

[stevengj]

scipy.special.beta gives the correct answer here.

[stevengj]

It would be easy to simply check for integer arguments. However, it is important to make sure that we also correctly handle any arguments close to (-n,n).

[stevengj]

The scipy beta implementation is from Cephes, and includes some special-case code that we are missing:

• It has special handling for negative-integer arguments to avoid NaN.
• It has special handling of beta(a,b) for |a| >> |b| to avoid cancellation errors.
• For smallish arguments, it calls gamma directly rather than exponentiating lgamma; I suspect this is mostly a performance optimization.

[andreasnoack]

What is the defintion for the beta function with negative (integer) first argument? It appears to me that the usual integral definition doesn't converge for e.g. (-2,2). Negative integer arguments for x are allowed in Mathematica but I couldn't find any explanation or definition in their documentation. Abramowitz and Stegun as well as Wikipedia require positive real part of both arguments.

[CarloLucibello]

It is defined through the gamma functions I think, for negative x and y the integral representation doesn't hold anymore.

If you agree on this I could create a PR addressing this bug mimicking the scipy implementation

[stevengj]

If you take the limit of the gamma function definition as (z,n) approaches (-n,n), you get a finite and well-defined result. Besides being unambiguous, this seems to be totally standard; e.g. it is what Mathematica returns.

We should definitely implement the tricks from the Cephes implementation, which is well-respected.

[andreasnoack]

I have no problem with extending the beta function, but

...this seems to be totally standard

is not a fair description. Cephes doesn't support negative arguments. Support for that was a relatively recent modification. Neither does the Naval Surface Warfare Center library and, as mentioned, Abramowitz and Stegun doesn't define the function for negative arguments.

[stevengj]

@andreasnoack, thanks for the clarification, I didn't realize it was a recent modification to Cephes.

Maple also defines beta for negative integers by the limit.

[andreasnoack]

Thanks for the link. In contrast to Mathematica, Maple's documentation states the definition for negative integers. It would be great if we included a similar comment about negative integer x in our documentation.

[simonbyrne]

I mentioned the |a| >> |b| problem here.

I remember looking at the Cephes code before, and was a bit underwhelmed, but there probably are some good ideas there. Also, I never managed to find its licence status (though presumably whomever has copyright presumably doesn't care that much).

[stevengj]

@simonbyrne, apparently the Cephes author gave permission to distribute it under 3-clause BSD; see e.g. here. Also, apparently Debian inquired about this 11 years ago and eventually got 3-clause BSD permission: the /usr/share/doc/python-scipy/copyright file on my Debian machine lists:

Files: scipy/special/cephes/*
Copyright: 1984, 1987, 1989, 2000 Stephen L. Moshier
License: BSD-3-clause

So, as long as we use the Cephes code from SciPy, we should be okay.

[simonbyrne]

Good to know, thanks.