Repository: scipy/scipy
Pull request: MAINT: update boost_math
Closes #18302

Problem statement:
BUG: beta.pdf is broken on main (1.11.0.dev0)
### Describe your issue.

beta.pdf has changes from 1.10 to the 1.11.0.dev0 build and now returns incorrect results.




### Reproducing Code Example

```python
import numpy as np
from scipy import stats

xx = np.array([0.        , 0.11111111, 0.22222222, 0.33333333, 0.44444444, 0.55555556, 0.66666667, 0.77777778, 0.88888889, 1.        ])
k = np.array([0., 1., 2., 3., 4.])
n=4.0

pdf = stats.beta.pdf(xx[..., None], k + 1, n - k + 1)


np.testing.assert_allclose(pdf[-1,-1],5.0)

```
Passes on 1.10.0, fails on 1.11.0.dev0.



### Error message

```shell
AssertionError: 
Not equal to tolerance rtol=1e-07, atol=0
Mismatched elements: 1 / 1 (100%)
Max absolute difference: 5.
Max relative difference: 1.
 x: array(0.)
 y: array(5.)
```


### SciPy/NumPy/Python version and system information

```shell
1.11.0.dev0+1817.98c51a6 1.23.5 sys.version_info(major=3, minor=11, micro=2, releaselevel='final', serial=0)
{
  "Compilers": {
    "c": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "cc"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "0.29.34",
      "commands": "cython"
    },
    "c++": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "c++"
    },
    "fortran": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "gfortran"
    },
    "pythran": {
      "version": "0.12.1",
      "include directory": "C:\\Users\\runneradmin\\AppData\\Local\\Temp\\pip-build-env-024pp_ot\\overlay\\Lib\\site-packages/pythran"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "windows"
    },
    "build": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "windows"
    },
    "cross-compiled": false
  },
  "Build Dependencies": {
    "blas": {
      "name": "openblas",
      "found": true,
      "version": "0.3.21.dev",
      "detection method": "pkgconfig",
      "include directory": "/c/opt/64/include",
      "lib directory": "/c/opt/64/lib",
      "openblas configuration": "USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= SKYLAKEX MAX_THREADS=2",
      "pc file directory": "c:/opt/64/lib/pkgconfig"
    },
    "lapack": {
      "name": "openblas",
      "found": true,
      "version": "0.3.21.dev",
      "detection method": "pkgconfig",
      "include directory": "/c/opt/64/include",
      "lib directory": "/c/opt/64/lib",
      "openblas configuration": "USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= SKYLAKEX MAX_THREADS=2",
      "pc file directory": "c:/opt/64/lib/pkgconfig"
    }
  },
  "Python Information": {
    "path": "C:\\Users\\runneradmin\\AppData\\Local\\Temp\\cibw-run-cgdr6yw7\\cp311-win_amd64\\build\\venv\\Scripts\\python.exe",
    "version": "3.11"
  }
}
```
```

Hints from issue discussion:
Creating a failure on statsmodels pre-release testing.

Thanks @bashtage I confirm the issue. A smaller reproducer is

```python
stats.beta.pdf(1e-30, a=1, b=5)
# 5.000000000000035
stats.beta.pdf(0, a=1, b=5)
# 0.0
```
I can get to the actually error playing with the float

```python
stats.beta.pdf(1e-310, a=1, b=5)
Traceback (most recent call last):
  File ".../envs/scipy-dev-311/lib/python3.11/site-packages/IPython/core/interactiveshell.py", line 3505, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-50-530e9aaceec0>", line 1, in <module>
    stats.beta.pdf(1e-310, a=1, b=5)
  File ".../scipy/build-install/lib/python3.11/site-packages/scipy/stats/_distn_infrastructure.py", line 1992, in pdf
    place(output, cond, self._pdf(*goodargs) / scale)
                        ^^^^^^^^^^^^^^^^^^^^
  File ".../scipy/build-install/lib/python3.11/site-packages/scipy/stats/_continuous_distns.py", line 681, in _pdf
    return _boost._beta_pdf(x, a, b)
           ^^^^^^^^^^^^^^^^^^^^^^^^^
OverflowError: Error in function ibeta_derivative<d>(%1%,%1%,%1%): Overflow Error
```

Seems like it's a boost thing. cc @mckib2 
Ping @jzmaddock @mborland 
Interesting; I will look at it. Does Scipy track the current Boost.Math develop, or is pinned to a particular commit?
Thanks @mborland!

It's tracking `7203fa2def6347b0d5f8fe1e8522d5b0a618db9d` but I could reproduce on `109a814e89f77ff8a3fc8f0391f6b35a12640669` (develop). If I did this correctly, I updated and rebuilt from scratch 😅
There are two issues: the return value at the boundary when the shape parameter associated with that boundary is 1, and the overflow when `x` is subnormal (e.g. 1e-310).  @mborland has a fix in progress at https://github.com/boostorg/math/pull/977.

In SciPy 1.10.1, `beta.pdf(0, 1, 5)` returns 5.  That is the correct limit as $x \rightarrow 0^+$.  Currently, `beta.pdf(0, 1, 5)` returns 0, and that appears to be what the beta distribution in boost/math has always done (I haven't checked the history thoroughly).  The justification there is that Wolfram Alpha returns 0.  But the "correct" value of the PDF at the boundaries is more often a matter of convention than mathematical rigor (@rkern wrote a nice summary of this, but I don't recall where).  We can restore the old convention in SciPy by adding an additional check for this special case here:

https://github.com/scipy/scipy/blob/b07ed536bdd1bee735e0ded1dbe5bf32ac209604/scipy/stats/_boost/include/func_defs.hpp#L71-L88

E.g. this should do it:

```C++
template<template <typename, typename> class Dst, class RealType, class...Args>
RealType
boost_pdf_beta(const RealType x, const RealType a, const RealType b)
{
    if (std::isfinite(x)) {
        if ((x >= 1) && (b < 1)) {
            // x>1 should really be 0, but rv_continuous will do that for us
            return INFINITY;
        }
        else if ((x <= 0) && (a < 1)) {
            return INFINITY;
        }
        else if ((x == 0) && (a == 1)) {
            return b;
        }
        else if ((x == 1) && (b == 1)) {
            return a;
        }
        return boost::math::pdf(
	    boost::math::beta_distribution<RealType, Policy>(a, b), x);
    }
    return NAN;
}
```
For Wolfram alpha, not sure how you checked. I personally just looked at the graph here https://www.wolframalpha.com/input?i=Beta%281%2C+5%29&assumption=%7B%22C%22%2C+%22Beta%22%7D+-%3E+%7B%22ProbabilityDistribution%22%7D
You can use [`PDF[BetaDistribution[1, 5], 0]`](https://www.wolframalpha.com/input?i=PDF%5BBetaDistribution%5B1%2C+5%5D%2C+0%5D).
What do you get with scipy main for the upper bound?

```
beta.pdf(0, 1, 5), beta.pdf(1, 5, 1)
(5.0, 5.0)
scipy.__version__
'1.7.3'
```

The application in statsmodels is for bernstein polynomials.
With this old behavior of beta.pdf at the boundaries, a theoretically implied round trip test for bernstein polynomials also works at the 0, 1 boundary values. This test fails now at the upper boundary with scipy main.

And it looks like those are the correct values by continuity.

Commit message: ENH/TST: filter overflow warning and better tests

Patch to review:
diff --git a/scipy/_lib/boost_math b/scipy/_lib/boost_math
index 109a814e89f7..298a243ccd36 160000
--- a/scipy/_lib/boost_math
+++ b/scipy/_lib/boost_math
@@ -1 +1 @@
-Subproject commit 109a814e89f77ff8a3fc8f0391f6b35a12640669
+Subproject commit 298a243ccd3639b6eaa59bcdab7ab9d5f008fb36
diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py
index 20e75371e6e8..853dfaeb586f 100644
--- a/scipy/stats/_continuous_distns.py
+++ b/scipy/stats/_continuous_distns.py
@@ -678,6 +678,8 @@ def _pdf(self, x, a, b):
         #                     gamma(a+b) * x**(a-1) * (1-x)**(b-1)
         # beta.pdf(x, a, b) = ------------------------------------
         #                              gamma(a)*gamma(b)
+        message = "overflow encountered in _beta_pdf"
+        warnings.filterwarnings('ignore', message=message)
         return _boost._beta_pdf(x, a, b)
 
     def _logpdf(self, x, a, b):
diff --git a/scipy/stats/tests/test_boost_ufuncs.py b/scipy/stats/tests/test_boost_ufuncs.py
index 55797953c8ca..e3b3950896a6 100644
--- a/scipy/stats/tests/test_boost_ufuncs.py
+++ b/scipy/stats/tests/test_boost_ufuncs.py
@@ -20,6 +20,8 @@
     (_boost._beta_cdf, (0.5, 2, 3), 0.6875),
     (_boost._beta_ppf, (0.6875, 2, 3), 0.5),
     (_boost._beta_pdf, (0.5, 2, 3), 1.5),
+    (_boost._beta_pdf, (0, 1, 5), 5.0),
+    (_boost._beta_pdf, (1, 5, 1), 5.0),
     (_boost._beta_sf, (0.5, 2, 1), 0.75),
     (_boost._beta_isf, (0.75, 2, 1), 0.5),
     (_boost._binom_cdf, (1, 3, 0.5), 0.5),
diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py
index c6764577f410..b4bab6682498 100644
--- a/scipy/stats/tests/test_distributions.py
+++ b/scipy/stats/tests/test_distributions.py
@@ -4037,6 +4037,18 @@ def test_endpoints(self):
         a, b = 0.2, 3
         assert_equal(stats.beta.pdf(0, a, b), np.inf)
 
+        # Confirm that boost's beta distribution returns 5 at x=0
+        # when a=1, b=5
+        a, b = 1, 5
+        assert_equal(stats.beta.pdf(0, a, b), 5)
+        assert_equal(stats.beta.pdf(1e-310, a, b), 5)
+
+        # Confirm that boost's beta distribution returns 5 at x=1
+        # when a=5, b=1
+        a, b = 5, 1
+        assert_equal(stats.beta.pdf(1, a, b), 5)
+        assert_equal(stats.beta.pdf(1-1e-310, a, b), 5)
+
     @pytest.mark.xfail(IS_PYPY, reason="Does not convert boost warning")
     def test_boost_eval_issue_14606(self):
         q, a, b = 0.995, 1.0e11, 1.0e13


Code review comments:
diff_hunk: @@ -678,6 +678,8 @@ def _pdf(self, x, a, b):
         #                     gamma(a+b) * x**(a-1) * (1-x)**(b-1)
         # beta.pdf(x, a, b) = ------------------------------------
         #                              gamma(a)*gamma(b)
+        message = "overflow encountered in _beta_pdf"
+        warnings.filterwarnings('ignore', message=message)
path: scipy/stats/_continuous_distns.py
text: @user1:
If this is necessary, it should be handled with [`numpy.errstate`](https://numpy.org/doc/stable/reference/generated/numpy.errstate.html) instead of using the `warnings` module.

@author:
Why?

@author:
Tried it locally and it passes, I will commit.

Is `np.errstate` doing the same and just a convenience for consistency??