Repository: scipy/scipy
Pull request: MAINT: stats: ensure that `rv_continuous._fitstart` shapes do not violate `_argcheck`
#### Reference issue
Closes gh-15689

#### What does this implement/fix?
gh-15689 reported that when using the `fit` method of `rv_continuous`, the `_fitstart` still needs to be overridden for several distributions because the default value of shape parameters (all `1`) violates `_argcheck`. This PR overrides `_fitstart` as needed and adds a test to ensure that new distributions do not have the same problem.

#### Additional information
Most of shape values produced by these overrides are arbitrary, but better than all `1`, and I think the need to improve distribution fitting on a case-by-case basis is adequately documented outside of gh-15689.

Problem statement:
BUG: default value of shape parameter in fit method of rv_continuous might not be admissible
### Describe your issue.

If no guess of the shape parameter is provided to `rv_continuous.fit`, the value of 1.0 is taken as a guess. However, it is not guaranteed that this value is allowed as a shape parameter, e.g., if `_argcheck(1.0)` is True

### Reproducing Code Example

```python
import numpy as np
from scipy.stats import rv_continuous
from scipy.optimize import fmin


class myDist_gen(rv_continuous):
    def _argcheck(self, a):
        return a > 2

    def _pdf(self, x, a):
        # we need a >= 2 to avoid division by 0 at x=Pi
        return (1 + np.cos(x))**(a-2)

myDist = myDist_gen(a=0.0, b=4.0, name='bla') 
data = myDist.rvs(a=3, size=20)
myDist.fit(data, floc=0, fscale=1)
```


### Error message

```shell
Exception: Optimization converged to parameters that are outside the range allowed by the distribution
```


### SciPy/NumPy/Python version information

1.7.1 1.21.2 sys.version_info(major=3, minor=9, micro=7, releaselevel='final', serial=0)

Hints from issue discussion:
To explain the root cause, I tried to replicate the behaviour of `fit` in `stats\_distn_infrastructure.py`:

```
args = ()
kwds = {'floc': 0.0, 'fscale': 1.0}

start = [None]*2
# fitstart will set the guess for a to 1.0
start = myDist._fitstart(data)
args += start[:-2]
loc = kwds.pop('loc', start[-2])
scale = kwds.pop('scale', start[-1])
args += (loc, scale)

x0, func, restore, args = myDist._reduce_func(args, kwds, data=data)
vals = fmin(func, x0, args=(np.ravel(data),), disp=1) # array([1.])

obj = func(vals, data) # np.inf

if not (np.all(self._argcheck(*shapes)) and scale > 0):
            raise Exception("Optimization converged to parameters that are "
                            "outside the range allowed by the distribution.")

```

From the above, it is clear that the optimizer did not converge...

If one specifies a valid guess for a, it works:
```
myDist.fit(data, 2.2, floc=0, fscale=1) # I got (1028.5702001953148, 0, 1)
```
If we have `ShapeInfo` for all distributions, one could make sure that the guess of the shape parameter is valid
Yup. I [had asked](https://github.com/scipy/scipy/pull/15436#issuecomment-1043892084) whether you'd like me to open PRs for those (adding `ShapeInfo`) or whether you'd like to do it and I can review. The advantage of the latter is that it would give others a reason to test `stats.fit` out before 1.9.0 branches. 
Here is the start of a test that we can add to `test_fit.py` to check for when the `_fitstart` shape guess fails `_argcheck`.

```python3
def cases_test_fitstart():
    for distname, shapes in dict(distcont).items():
        if not isinstance(distname, str):
            continue
        yield distname, shapes


@pytest.mark.parametrize('distname, shapes', cases_test_fitstart())
def test_fitstart(distname, shapes):
    dist = getattr(stats, distname)
    data = np.random.rand(5)

    with np.errstate(invalid='ignore', divide='ignore'):  # irrelevant to test
        guess = dist._fitstart(data)

    assert dist._argcheck(*guess[:-2])
```
 
We're going to have to define `_fitstart` for at least some of these manually. For instance, `genhyperbolic` fails because it requires that `abs(b) < a` with strict inequality; `_shape_info` won't tell us that. For those that `_shape_info` is enough information to produce a valid guess, how were you thinking of using it? The maximum of `1` and some distance from the left bound of the parameter's domain? I'm thinking it might be easier to just manually override `_fitstart` of distributions that fail the test above.

<details>

```
=================================================================== FAILURES ====================================================================
______________________________________________________ test_fitstart[crystalball-shapes13] ______________________________________________________
..\..\scipy\stats\tests\test_fit.py:157: in test_fitstart
    assert dist._argcheck(*guess[:-2])
E   assert False
E    +  where False = <bound method crystalball_gen._argcheck of <scipy.stats._continuous_distns.crystalball_gen object at 0x0000011AE3DE5520>>(*(1.0, 1.0))
E    +    where <bound method crystalball_gen._argcheck of <scipy.stats._continuous_distns.crystalball_gen object at 0x0000011AE3DE5520>> = <scipy.stats._continuous_distns.crystalball_gen object at 0x0000011AE3DE5520>._argcheck
        data       = array([0.5384155 , 0.50860436, 0.18535913, 0.31765469, 0.06954131])
        dist       = <scipy.stats._continuous_distns.crystalball_gen object at 0x0000011AE3DE5520>
        distname   = 'crystalball'
        guess      = (1.0, 1.0, 0, 1)
        shapes     = (2.0, 3.0)
_____________________________________________________ test_fitstart[genhyperbolic-shapes32] _____________________________________________________
..\..\scipy\stats\tests\test_fit.py:157: in test_fitstart
    assert dist._argcheck(*guess[:-2])
E   assert False
E    +  where False = <bound method genhyperbolic_gen._argcheck of <scipy.stats._continuous_distns.genhyperbolic_gen object at 0x0000011AE3D859A0>>(*(1.0, 1.0, 1.0))
E    +    where <bound method genhyperbolic_gen._argcheck of <scipy.stats._continuous_distns.genhyperbolic_gen object at 0x0000011AE3D859A0>> = <scipy.stats._continuous_distns.genhyperbolic_gen object at 0x0000011AE3D859A0>._argcheck
        data       = array([0.56641824, 0.66658558, 0.79312473, 0.28752532, 0.28519154])
        dist       = <scipy.stats._continuous_distns.genhyperbolic_gen object at 0x0000011AE3D859A0>
        distname   = 'genhyperbolic'
        guess      = (1.0, 1.0, 1.0, 0, 1)
        shapes     = (0.5, 1.5, -0.5)
______________________________________________________ test_fitstart[loguniform-shapes65] _______________________________________________________
..\..\scipy\stats\tests\test_fit.py:157: in test_fitstart
    assert dist._argcheck(*guess[:-2])
E   assert False
E    +  where False = <bound method reciprocal_gen._argcheck of <scipy.stats._continuous_distns.reciprocal_gen object at 0x0000011AE3DE51C0>>(*(1.0, 1.0))
E    +    where <bound method reciprocal_gen._argcheck of <scipy.stats._continuous_distns.reciprocal_gen object at 0x0000011AE3DE51C0>> = <scipy.stats._continuous_distns.reciprocal_gen object at 0x0000011AE3DE51C0>._argcheck
        data       = array([0.7307671 , 0.10967873, 0.43040378, 0.39785707, 0.06877716])
        dist       = <scipy.stats._continuous_distns.reciprocal_gen object at 0x0000011AE3DE51C0>
        distname   = 'loguniform'
        guess      = (1.0, 1.0, 0, 1)
        shapes     = (0.01, 1.25)
_____________________________________________________ test_fitstart[norminvgauss-shapes75] ______________________________________________________
..\..\scipy\stats\tests\test_fit.py:157: in test_fitstart
    assert dist._argcheck(*guess[:-2])
E   assert False
E    +  where False = <bound method norminvgauss_gen._argcheck of <scipy.stats._continuous_distns.norminvgauss_gen object at 0x0000011AE3D95190>>(*(1.0, 1.0))
E    +    where <bound method norminvgauss_gen._argcheck of <scipy.stats._continuous_distns.norminvgauss_gen object at 0x0000011AE3D95190>> = <scipy.stats._continuous_distns.norminvgauss_gen object at 0x0000011AE3D95190>._argcheck
        data       = array([0.83545029, 0.51939515, 0.18484643, 0.16113219, 0.01119876])
        dist       = <scipy.stats._continuous_distns.norminvgauss_gen object at 0x0000011AE3D95190>
        distname   = 'norminvgauss'
        guess      = (1.0, 1.0, 0, 1)
        shapes     = (1.25, 0.5)
______________________________________________________ test_fitstart[reciprocal-shapes84] _______________________________________________________
..\..\scipy\stats\tests\test_fit.py:157: in test_fitstart
    assert dist._argcheck(*guess[:-2])
E   assert False
E    +  where False = <bound method reciprocal_gen._argcheck of <scipy.stats._continuous_distns.reciprocal_gen object at 0x0000011AE3DD9280>>(*(1.0, 1.0))
E    +    where <bound method reciprocal_gen._argcheck of <scipy.stats._continuous_distns.reciprocal_gen object at 0x0000011AE3DD9280>> = <scipy.stats._continuous_distns.reciprocal_gen object at 0x0000011AE3DD9280>._argcheck
        data       = array([0.1419588 , 0.97819743, 0.02484348, 0.49905794, 0.65659307])
        dist       = <scipy.stats._continuous_distns.reciprocal_gen object at 0x0000011AE3DD9280>
        distname   = 'reciprocal'
        guess      = (1.0, 1.0, 0, 1)
        shapes     = (0.01, 1.25)
___________________________________________________ test_fitstart[studentized_range-shapes89] ___________________________________________________
..\..\scipy\stats\tests\test_fit.py:155: in test_fitstart
    guess = dist._fitstart(data)
        data       = array([0.17973032, 0.78299111, 0.28032091, 0.6156475 , 0.74499654])
        dist       = <scipy.stats._continuous_distns.studentized_range_gen object at 0x0000011AE3DF69D0>
        distname   = 'studentized_range'
        shapes     = (3.0, 10.0)
..\..\scipy\stats\_continuous_distns.py:9313: in _fitstart
    return super(studentized_range_gen, self)._fitstart(data, args=(2, 1))
        __class__  = <class 'scipy.stats._continuous_distns.studentized_range_gen'>
        data       = array([0.17973032, 0.78299111, 0.28032091, 0.6156475 , 0.74499654])
        self       = <scipy.stats._continuous_distns.studentized_range_gen object at 0x0000011AE3DF69D0>
..\..\scipy\stats\_distn_infrastructure.py:2427: in _fitstart
    loc, scale = self._fit_loc_scale_support(data, *args)
        args       = (2, 1)
        data       = array([0.17973032, 0.78299111, 0.28032091, 0.6156475 , 0.74499654])
        self       = <scipy.stats._continuous_distns.studentized_range_gen object at 0x0000011AE3DF69D0>
..\..\scipy\stats\_distn_infrastructure.py:2720: in _fit_loc_scale_support
    loc_hat, scale_hat = self.fit_loc_scale(data, *args)
        args       = (2, 1)
        data       = array([0.17973032, 0.78299111, 0.28032091, 0.6156475 , 0.74499654])
        self       = <scipy.stats._continuous_distns.studentized_range_gen object at 0x0000011AE3DF69D0>
..\..\scipy\stats\_distn_infrastructure.py:2783: in fit_loc_scale
    mu, mu2 = self.stats(*args, **{'moments': 'mv'})
        args       = (2, 1)
        data       = array([0.17973032, 0.78299111, 0.28032091, 0.6156475 , 0.74499654])
        self       = <scipy.stats._continuous_distns.studentized_range_gen object at 0x0000011AE3DF69D0>
..\..\scipy\stats\_distn_infrastructure.py:1160: in stats
    mu = self._munp(1, *goodargs)
        args       = (array(2), array(1))
        cond       = True
        default    = array(nan)
        g1         = None
        g2         = None
        goodargs   = [array([2]), array([1])]
        kwds       = {'moments': 'mv'}
        loc        = array([0])
        moments    = 'mv'
        mu         = None
        mu2        = None
        output     = []
        scale      = array([1])
        self       = <scipy.stats._continuous_distns.studentized_range_gen object at 0x0000011AE3DF69D0>
..\..\scipy\stats\_continuous_distns.py:9334: in _munp
    return np.float64(ufunc(K, k, df))
        K          = 1
        _a         = 0
        _b         = inf
        _single_moment = <function studentized_range_gen._munp.<locals>._single_moment at 0x0000011AE5231040>
        cython_symbol = '_studentized_range_moment'
        df         = array([1])
        k          = array([2])
        self       = <scipy.stats._continuous_distns.studentized_range_gen object at 0x0000011AE3DF69D0>
        ufunc      = <ufunc '_single_moment (vectorized)'>
..\..\scipy\stats\_continuous_distns.py:9331: in _single_moment
    return integrate.nquad(llc, ranges=ranges, opts=opts)[0]
        K          = 1
        _a         = 0
        _b         = inf
        arg        = [1, 2, 1, 0.46735582791521785]
        cython_symbol = '_studentized_range_moment'
        df         = 1
        k          = 2
        llc        = LowLevelCallable(<capsule object "double (int, double *, void *)" at 0x0000011AE3C94840>, c_void_p(1214748133008))
        log_const  = 0.46735582791521785
        opts       = {'epsabs': 1e-11, 'epsrel': 1e-12}
        ranges     = [(-inf, inf), (0, inf), (0, inf)]
        usr_data   = c_void_p(1214748133008)
..\..\scipy\integrate\_quadpack_py.py:824: in nquad
    return _NQuad(func, ranges, opts, full_output).integrate(*args)
        args       = ()
        depth      = 3
        full_output = False
        func       = LowLevelCallable(<capsule object "double (int, double *, void *)" at 0x0000011AE3C94840>, c_void_p(1214748133008))
        opts       = [<scipy.integrate._quadpack_py._OptFunc object at 0x0000011AE5398910>, <scipy.integrate._quadpack_py._OptFunc object at 0x0000011AE5398910>, <scipy.integrate._quadpack_py._OptFunc object at 0x0000011AE5398910>]
        ranges     = [<scipy.integrate._quadpack_py._RangeFunc object at 0x0000011AE5398940>, <scipy.integrate._quadpack_py._RangeFunc object at 0x0000011AE5398D30>, <scipy.integrate._quadpack_py._RangeFunc object at 0x0000011AE5398BE0>]
..\..\scipy\integrate\_quadpack_py.py:878: in integrate
    quad_r = quad(f, low, high, args=args, full_output=self.full_output,
        args       = ()
        depth      = 0
        f          = functools.partial(<bound method _NQuad.integrate of <scipy.integrate._quadpack_py._NQuad object at 0x0000011AE518F730>>, depth=1)
        fn_opt     = <scipy.integrate._quadpack_py._OptFunc object at 0x0000011AE5398910>
        fn_range   = <scipy.integrate._quadpack_py._RangeFunc object at 0x0000011AE5398BE0>
        high       = inf
        ind        = -1
        kwargs     = {}
        low        = 0
        opt        = {'epsabs': 1e-11, 'epsrel': 1e-12}
        self       = <scipy.integrate._quadpack_py._NQuad object at 0x0000011AE518F730>
..\..\scipy\integrate\_quadpack_py.py:400: in quad
    warnings.warn(msg, IntegrationWarning, stacklevel=2)
E   scipy.integrate._quadpack_py.IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
E     If increasing the limit yields no improvement it is advised to analyze
E     the integrand in order to determine the difficulties.  If the position of a
E     local difficulty can be determined (singularity, discontinuity) one will
E     probably gain from splitting up the interval and calling the integrator
E     on the subranges.  Perhaps a special-purpose integrator should be used.
        a          = 0
        args       = ()
        b          = inf
        epsabs     = 1e-11
        epsrel     = 1e-12
        flip       = False
        full_output = False
        func       = functools.partial(<bound method _NQuad.integrate of <scipy.integrate._quadpack_py._NQuad object at 0x0000011AE518F730>>, depth=1)
        ier        = 1
        limit      = 50
        limlst     = 50
        maxp1      = 50
        msg        = 'The maximum number of subdivisions (50) has been achieved.\n  If increasing the limit yields no improvement it is adv...up the interval and calling the integrator \n  on the subranges.  Perhaps a special-purpose integrator should be used.'
        msgs       = {80: 'A Python error occurred possibly while calling the function.', 1: 'The maximum number of subdivisions (50) has b...\n  underestimated.', 3: 'Extremely bad integrand behavior occurs at some points of the\n  integration interval.', ...}
        points     = None
        retval     = (7.2434775656499895, 2.9875944049538376e-11, 1)
        weight     = None
        wopts      = None
        wvar       = None
_______________________________________________________ test_fitstart[truncnorm-shapes94] _______________________________________________________
..\..\scipy\stats\tests\test_fit.py:157: in test_fitstart
    assert dist._argcheck(*guess[:-2])
E   assert False
E    +  where False = <bound method truncnorm_gen._argcheck of <scipy.stats._continuous_distns.truncnorm_gen object at 0x0000011AE3DE5460>>(*(1.0, 1.0))
E    +    where <bound method truncnorm_gen._argcheck of <scipy.stats._continuous_distns.truncnorm_gen object at 0x0000011AE3DE5460>> = <scipy.stats._continuous_distns.truncnorm_gen object at 0x0000011AE3DE5460>._argcheck
        data       = array([0.79925833, 0.4462647 , 0.90538196, 0.51815275, 0.17035357])
        dist       = <scipy.stats._continuous_distns.truncnorm_gen object at 0x0000011AE3DE5460>
        distname   = 'truncnorm'
        guess      = (1.0, 1.0, 0, 1)
        shapes     = (0.1, 2.0)
___________________________________________________ test_fitstart[truncweibull_min-shapes95] ____________________________________________________
..\..\scipy\stats\tests\test_fit.py:157: in test_fitstart
    assert dist._argcheck(*guess[:-2])
E   assert False
E    +  where False = <bound method truncweibull_min_gen._argcheck of <scipy.stats._continuous_distns.truncweibull_min_gen object at 0x0000011AE3D78280>>(*(1.0, 1.0, 1.0))
E    +    where <bound method truncweibull_min_gen._argcheck of <scipy.stats._continuous_distns.truncweibull_min_gen object at 0x0000011AE3D78280>> = <scipy.stats._continuous_distns.truncweibull_min_gen object at 0x0000011AE3D78280>._argcheck
        data       = array([0.64502691, 0.03150199, 0.10353142, 0.37198364, 0.18305176])
        dist       = <scipy.stats._continuous_distns.truncweibull_min_gen object at 0x0000011AE3D78280>
        distname   = 'truncweibull_min'
        guess      = (1.0, 1.0, 1.0, 0, 1)
        shapes     = (2.5, 0.25, 1.75)
============================================================ short test summary info ============================================================
FAILED ..\..\scipy\stats\tests\test_fit.py::test_fitstart[crystalball-shapes13] - assert False
FAILED ..\..\scipy\stats\tests\test_fit.py::test_fitstart[genhyperbolic-shapes32] - assert False
FAILED ..\..\scipy\stats\tests\test_fit.py::test_fitstart[loguniform-shapes65] - assert False
FAILED ..\..\scipy\stats\tests\test_fit.py::test_fitstart[norminvgauss-shapes75] - assert False
FAILED ..\..\scipy\stats\tests\test_fit.py::test_fitstart[reciprocal-shapes84] - assert False
FAILED ..\..\scipy\stats\tests\test_fit.py::test_fitstart[studentized_range-shapes89] - scipy.integrate._quadpack_py.IntegrationWarning: The ma...
FAILED ..\..\scipy\stats\tests\test_fit.py::test_fitstart[truncnorm-shapes94] - assert False
FAILED ..\..\scipy\stats\tests\test_fit.py::test_fitstart[truncweibull_min-shapes95] - assert False
========================================================= 8 failed, 96 passed in 41.52s =========================================================
```

</details>

Commit message: MAINT: stats: ensure that _fitstart is not trivially invalid

Patch to review:
diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py
index 28793457feca..06f56cb5f557 100644
--- a/scipy/stats/_continuous_distns.py
+++ b/scipy/stats/_continuous_distns.py
@@ -2316,6 +2316,10 @@ def _shape_info(self):
         ib = _ShapeInfo("b", False, (0, np.inf), (False, False))
         return [ic, ia, ib]
 
+    def _fitstart(self, data):
+        # Arbitrary, but default a=b=c=1 is not valid
+        return super()._fitstart(data, args=(1, 0, 1))
+
     def _get_support(self, c, a, b):
         return a, b
 
@@ -3338,6 +3342,11 @@ def _shape_info(self):
         ib = _ShapeInfo("b", False, (-np.inf, np.inf), (False, False))
         return [ip, ia, ib]
 
+    def _fitstart(self, data):
+        # Arbitrary, but the default a=b=1 is not valid
+        return super()._fitstart(data, args=(1, 1, 0.5))
+
+
     def _logpdf(self, x, p, a, b):
         # kve instead of kv works better for large values of p
         # and smaller values of sqrt(a^2  - b^2)
@@ -4485,6 +4494,10 @@ def _shape_info(self):
         ib = _ShapeInfo("b", False, (-np.inf, np.inf), (False, False))
         return [ia, ib]
 
+    def _fitstart(self, data):
+        # Arbitrary, but the default a=b=1 is not valid
+        return super()._fitstart(data, args=(1, 0.5))
+
     def _pdf(self, x, a, b):
         gamma = np.sqrt(a**2 - b**2)
         fac1 = a / np.pi * np.exp(gamma)
@@ -7289,6 +7302,10 @@ def _shape_info(self):
         ib = _ShapeInfo("b", False, (0, np.inf), (False, False))
         return [ia, ib]
 
+    def _fitstart(self, data):
+        # Reasonable, since support is [a, b]
+        return super()._fitstart(data, args=(np.min(data), np.max(data)))
+
     def _get_support(self, a, b):
         return a, b
 
@@ -8275,6 +8292,10 @@ def _shape_info(self):
         ib = _ShapeInfo("b", False, (-np.inf, np.inf), (False, True))
         return [ia, ib]
 
+    def _fitstart(self, data):
+        # Reasonable, since support is [a, b]
+        return super()._fitstart(data, args=(np.min(data), np.max(data)))
+
     def _get_support(self, a, b):
         return a, b
 
@@ -9134,6 +9155,10 @@ def _shape_info(self):
         im = _ShapeInfo("m", False, (1, np.inf), (False, False))
         return [ibeta, im]
 
+    def _fitstart(self, data):
+        # Arbitrary, but the default m=1 is not valid
+        return super()._fitstart(data, args=(1, 1.5))
+
     def _pdf(self, x, beta, m):
         """
         Return PDF of the crystalball function.
diff --git a/scipy/stats/tests/test_fit.py b/scipy/stats/tests/test_fit.py
index 65f5cbd12499..143cc4d23daa 100644
--- a/scipy/stats/tests/test_fit.py
+++ b/scipy/stats/tests/test_fit.py
@@ -212,6 +212,25 @@ def cases_test_fit():
             yield dist
 
 
+def cases_test_fitstart():
+    for distname, shapes in dict(distcont).items():
+        if not isinstance(distname, str) or distname in {'studentized_range'}:
+            continue
+        yield distname, shapes
+
+
+@pytest.mark.parametrize('distname, shapes', cases_test_fitstart())
+def test_fitstart(distname, shapes):
+    dist = getattr(stats, distname)
+    rng = np.random.default_rng(216342614)
+    data = rng.random(10)
+
+    with np.errstate(invalid='ignore', divide='ignore'):  # irrelevant to test
+        guess = dist._fitstart(data)
+
+    assert dist._argcheck(*guess[:-2])
+
+
 def assert_nllf_less_or_close(dist, data, params1, params0, rtol=1e-7, atol=0):
     nllf1 = dist.nnlf(params1, data)
     nllf0 = dist.nnlf(params0, data)


Code review comments:
diff_hunk: @@ -3338,6 +3342,11 @@ def _shape_info(self):
         ib = _ShapeInfo("b", False, (-np.inf, np.inf), (False, False))
         return [ip, ia, ib]
 
+    def _fitstart(self, data):
+        # Arbitrary, but the default a=b=1 is not valid
+        return super()._fitstart(data, args=(1, 1, 0.5))
+
+
path: scipy/stats/_continuous_distns.py
text: @author:
```suggestion
```