poisson_conf_interval 'kraft-burrows-nousek' fails unexpectedly

[hamogu]

Description

Expected behavior

poisson_conf_interval 'kraft-burrows-nousek' calculated values or fails with an obvious error message.

Actual behavior

The error message comes form deep within scipy.optimize and it's not clear for a user how to react. It's particularly confusing here because, at least on my platform, the function could work if not for a hardcoded upper value of an interval. The problem is in line 1160 S_max = brentq(func, N - B, 100) (see pasted error output below). In the function, the upper end of the interval that brentq searches is hardcoded to 100. That's because (looking at the code comments and also speaking from memory as the person who originally implemented this) much larger value lead to numerical overflow error - scipy simple can't be used. However, for the numbers in my example, the calculation works when changing the hardcoded 100 to 200.

Notes for tackling this

So, this looks very simple - just increase the hardcoded 100. However, we want to be a bit careful not to increase it too much because nobody is helped by running into overflow's instead. So, one might instead have to decrease that "Don't use scipy with N > 100" rule to a smaller number - finding the balance here might require some experimentation and new numbers in the test.

Steps to Reproduce

With scipy installed:

from astropy.stats import poisson_conf_interval
poisson_conf_interval(n=int(98), background=3.8, confidence_level=0.95,
                                             interval='kraft-burrows-nousek')
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [48], in <cell line: 1>()
----> 1 poisson_conf_interval(n=int(98), background=3.8, confidence_level=0.95,
      2                                              interval='kraft-burrows-nousek')

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/astropy/stats/funcs.py:764, in poisson_conf_interval(n, interval, sigma, background, confidence_level)
    762     if np.any(background < 0):
    763         raise ValueError('Background must be >= 0.')
--> 764     conf_interval = np.vectorize(_kraft_burrows_nousek,
    765                                  cache=True)(n, background, confidence_level)
    766     conf_interval = np.vstack(conf_interval)
    767 else:

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/numpy/lib/function_base.py:2304, in vectorize.__call__(self, *args, **kwargs)
   2301     vargs = [args[_i] for _i in inds]
   2302     vargs.extend([kwargs[_n] for _n in names])
-> 2304 return self._vectorize_call(func=func, args=vargs)

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/numpy/lib/function_base.py:2382, in vectorize._vectorize_call(self, func, args)
   2380     res = func()
   2381 else:
-> 2382     ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args)
   2384     # Convert args to object arrays first
   2385     inputs = [asanyarray(a, dtype=object) for a in args]

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/numpy/lib/function_base.py:2342, in vectorize._get_ufunc_and_otypes(self, func, args)
   2338     raise ValueError('cannot call `vectorize` on size 0 inputs '
   2339                      'unless `otypes` is set')
   2341 inputs = [arg.flat[0] for arg in args]
-> 2342 outputs = func(*inputs)
   2344 # Performance note: profiling indicates that -- for simple
   2345 # functions at least -- this wrapping can almost double the
   2346 # execution time.
   2347 # Hence we make it optional.
   2348 if self.cache:

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/astropy/stats/funcs.py:1290, in _kraft_burrows_nousek(N, B, CL)
   1288 if HAS_SCIPY and N <= 100:
   1289     try:
-> 1290         return _scipy_kraft_burrows_nousek(N, B, CL)
   1291     except OverflowError:
   1292         if not HAS_MPMATH:

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/astropy/stats/funcs.py:1160, in _scipy_kraft_burrows_nousek(N, B, CL)
   1157     out = eqn9_left(s_min, s, N, B)
   1158     return out[0] - CL
-> 1160 S_max = brentq(func, N - B, 100)
   1161 S_min = find_s_min(S_max, N, B)
   1162 return S_min, S_max

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/scipy/optimize/_zeros_py.py:783, in brentq(f, a, b, args, xtol, rtol, maxiter, full_output, disp)
    781 if rtol < _rtol:
    782     raise ValueError("rtol too small (%g < %g)" % (rtol, _rtol))
--> 783 r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
    784 return results_c(full_output, r)

ValueError: f(a) and f(b) must have different signs

System Details

Linux-3.10.0-1062.12.1.el7.x86_64-x86_64-with-glibc2.17
Python 3.9.7 (default, Sep 16 2021, 13:09:58)
[GCC 7.5.0]
Numpy 1.21.2
pyerfa 2.0.0
astropy 5.1
Scipy 1.7.1

[bharathr98]

Hi, I would like to tackle this issue

[ohyos-isolt]

Simply increase the hardcoded value (100) may not help. In my case, I actually need to reduce the hardcoded value to 70, other wise I got the same error.

from astropy.stats import poisson_conf_interval
poisson_conf_interval(n=int(79), background=0.54, confidence_level=0.998,
                                             interval='kraft-burrows-nousek')
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [48], in <cell line: 1>()
----> 1 poisson_conf_interval(n=int(98), background=3.8, confidence_level=0.95,
      2                                              interval='kraft-burrows-nousek')

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/astropy/stats/funcs.py:764, in poisson_conf_interval(n, interval, sigma, background, confidence_level)
    762     if np.any(background < 0):
    763         raise ValueError('Background must be >= 0.')
--> 764     conf_interval = np.vectorize(_kraft_burrows_nousek,
    765                                  cache=True)(n, background, confidence_level)
    766     conf_interval = np.vstack(conf_interval)
    767 else:

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/numpy/lib/function_base.py:2304, in vectorize.__call__(self, *args, **kwargs)
   2301     vargs = [args[_i] for _i in inds]
   2302     vargs.extend([kwargs[_n] for _n in names])
-> 2304 return self._vectorize_call(func=func, args=vargs)

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/numpy/lib/function_base.py:2382, in vectorize._vectorize_call(self, func, args)
   2380     res = func()
   2381 else:
-> 2382     ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args)
   2384     # Convert args to object arrays first
   2385     inputs = [asanyarray(a, dtype=object) for a in args]

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/numpy/lib/function_base.py:2342, in vectorize._get_ufunc_and_otypes(self, func, args)
   2338     raise ValueError('cannot call `vectorize` on size 0 inputs '
   2339                      'unless `otypes` is set')
   2341 inputs = [arg.flat[0] for arg in args]
-> 2342 outputs = func(*inputs)
   2344 # Performance note: profiling indicates that -- for simple
   2345 # functions at least -- this wrapping can almost double the
   2346 # execution time.
   2347 # Hence we make it optional.
   2348 if self.cache:

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/astropy/stats/funcs.py:1290, in _kraft_burrows_nousek(N, B, CL)
   1288 if HAS_SCIPY and N <= 100:
   1289     try:
-> 1290         return _scipy_kraft_burrows_nousek(N, B, CL)
   1291     except OverflowError:
   1292         if not HAS_MPMATH:

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/astropy/stats/funcs.py:1160, in _scipy_kraft_burrows_nousek(N, B, CL)
   1157     out = eqn9_left(s_min, s, N, B)
   1158     return out[0] - CL
-> 1160 S_max = brentq(func, N - B, 100)
   1161 S_min = find_s_min(S_max, N, B)
   1162 return S_min, S_max

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/scipy/optimize/_zeros_py.py:783, in brentq(f, a, b, args, xtol, rtol, maxiter, full_output, disp)
    781 if rtol < _rtol:
    782     raise ValueError("rtol too small (%g < %g)" % (rtol, _rtol))
--> 783 r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
    784 return results_c(full_output, r)

ValueError: f(a) and f(b) must have different signs

System detail
macOS 13.1
Python 3.10.9 (main, Dec 15 2022, 18:25:35) [Clang 14.0.0 (clang-1400.0.29.202)]
scipy 1.10.0
astropy 5.2.1
numpy 1.23.5

[efrain99]

Working on a solution for this issue but running into pre-commit errors. Hopefully I can resolve it soon so I can resolve this issue.

[pllim]

@efrain99 , would https://docs.astropy.org/en/latest/development/workflow/development_workflow.html#pre-commit help?

[efrain99]

@efrain99 , would https://docs.astropy.org/en/latest/development/workflow/development_workflow.html#pre-commit help?

Its apparently not an issue with pre-commit but an issue with an ssl certificate. So when I make a git commit pre-commit returns an ssl error. Apparently this seems to be an issue with windows... sigh

[pllim]

@efrain99 , did you install https://pypi.org/project/certifi/ ?

[efrain99]

@pllim , yes I just tried that but same error message unfortunately. I am going to keep looking thank you for the help!

[pllim]

@efrain99 , have you tried on WSL2 (assuming you are able to install it on that Windows machine)?

[salmanhiro]

Hi, I tried to change something in this PR #15208, but I still not fully understand about the reason behind S_max upper limit should be a certain value. Would be helpful if anyone can review the PR

Note: Well, still need to fix the pre commit error

Edit: I use the wrong branch to submit PR #15207. I move it to #15208

[datajungler]

@pllim , yes I just tried that but same error message unfortunately. I am going to keep looking thank you for the help!

During the pre-commit process, I encountered a similar issue with the SSL certificate. Despite attempting to install different versions of pip and conda, the error persisted. However, I found a solution by deactivating the virtual environment, which resolved the problem.