Repository: scipy/scipy
Pull request: BUG: levy_stable _fitstart fixes
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
https://docs.scipy.org/doc/numpy/dev/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
See #14470

#### What does this implement/fix?
<!--Please explain your changes.-->
Adds new tests for levy_stable._fitstart(), and fixes bugs so the tests now pass.

- A data sample and its mirror image (negative values) now return the same estimate of alpha.
- A data sample and its mirror image now return mirror image estimates of the location parameter (delta).
- A data sample with a negative location and a shifted sample with a slightly different negative location now return shifted estimates of the location parameter.

#### Additional information
<!--Any additional information you think is important.-->
There is still something flaky about the calculation of the location parameter (delta) when alpha==1, as it tries to calculate tan(pi/2) -> infinity.  In practice, it will be very rare that alpha==1. when estimated from data samples.

Problem statement:
levy_stable._fitstart() incorrect parameter estimates for negative beta, negative loc, or alpha==1
<!-- 
Thank you for taking the time to file a bug report. 
Please fill in the fields below, deleting the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
Note: This is a comment, and won't appear in the output.
-->

My issue is about scipy.levy_stable._fitstart().

- A data sample and its mirror image (negative values) should return the same estimate of alpha, but do not.
- A data sample and its mirror image should return mirror image estimates of the location parameter but do not.
- A data sample with a negative location and a shifted sample with a slightly different negative location should return shifted estimates of the location parameter but do not.
- A data sample with alpha==1 should not try to subtract infinity from a finite number but the calculation of the location parameter (delta) is trying to do that when it calculates tan(pi/2).  (My test code below induces an estimate of alpha==1 to exercise the relevant code, but does not currently check the value of the location parameter (delta) because I'm not sure what value is intended in this circumstance but I'm pretty sure it shouldn't involve trying to subtract tan(pi/2) from a finite value.)

#### Reproducing code example:
<!-- 
If you place your code between the triple backticks below, 
it will be rendered as a code block. 
-->

```
import pytest

# Test scipy.stats.levy_stable._fitstart().

import numpy as np
from scipy.stats import levy_stable
from scipy.stats import rv_continuous as rvc

class TestLevy_stable:
    def test_neg_beta(self):
        # Confirm that _fitstart() gets same estimate of alpha for +beta and -beta.
        rvc.random_state = 2021
        NPOINTS = 1001
        ALPHA = 1.5
        BETA_POS = 1
        
        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, size=NPOINTS)
        # Fit stable distribution to positive and negative mirror image data, examine alpha.
        alpha_pos, alpha_neg = [levy_stable._fitstart(data)[0] for data in [draw, -draw]]
        
        assert np.isclose(alpha_pos, alpha_neg)
        
    def test_neg_delta(self):
        # Confirm that _fitstart() returns symmetric +delta and -delta if beta flips.
        rvc.random_state = 2021
        NPOINTS = 1001
        ALPHA = 1.5
        BETA_POS = 1
        
        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, size=NPOINTS)
        # Fit stable distribution to positive and negative mirror image data, examine location.
        delta_pos, delta_neg = [levy_stable._fitstart(data)[2] for data in [draw, -draw]]
        
        assert np.isclose(delta_pos, -delta_neg)
        
    def test_delta_shift(self):
        # Confirm that _fitstart() returns delta that slides up and down if data shifts.
        rvc.random_state = 2021
        NPOINTS = 1001
        ALPHA = 1.5
        BETA_POS = 1
        LOC = -2
        SHIFT = 1
        
        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, loc=LOC, size=NPOINTS)
        # Fit stable distribution to data and shifted data, examine location parameter.
        delta0, delta_up = [levy_stable._fitstart(draw + x)[2] for x in [0, SHIFT]]
        
        assert np.isclose(delta0, delta_up-SHIFT)
        
    def test_alpha1(self):
        # Confirm that stable distribution can fit alpha==1 without blowing up.
        
        # Make a dataset whose quintiles lead to alpha==1.
        x5 = x25 = x50 = 5  # Percentiles 5, 10, 25 -> 5.
        x95 = 100 - x5  # Percentile 95.
        # Calculate value for percentile 75 that will interpolate alpha lookup to 1.
        #   We assume beta==1 so interpolation is between rows 4 and 5 of last column
        #   of _fitstart() alpha_table.
        x75 = x25 + (x95 - x5) / (4 + (1.150 - 1) / (1.150 - .973))
        #assert x5 <= x25 <= x50 <= x75 <= x95
        
        xx = np.array([x5, x25, x50, x75, x95])
        # np.percentile(xx, [5, 25, 50, 75, 95]) -> [5, 5, 5, ..., 95]
        
        # Place target percentiles either side.
        data = np.array([x5, x5,  # Percentiles 0, 10.
                         x25, x25,  # 20, 30.
                         x50, x50, x50,  # 40, 50, 60.
                         x75, x75,  # 70, 80.
                         x95, x95  # 90, 100.
                         ])
        #assert np.array_equal(np.percentile(data, [5, 25, 50, 75, 95]), xx)
        
        alpha, beta, loc, scale = levy_stable._fitstart(data)

        # Not sure what value of delta should be returned when alpha==1 so
        #  for now just verify that alpha==1.
        assert alpha == 1
        # Might want to check loc too--what value should it have when alpha==1?
        
        
        
        

#%%

if __name__ == "__main__":
    t = TestLevy_stable()
    t.test_neg_beta()
    t.test_neg_delta()
    t.test_delta_shift()
    t.test_alpha1()
```

#### pytest results:
```
% pytest
================================================== test session starts ===================================================
platform darwin -- Python 3.7.6, pytest-6.2.4, py-1.10.0, pluggy-0.13.1
rootdir: /Users/...
collected 4 items                                                                                                        

test_stats.py FFF.                                                                                                 [100%]

======================================================== FAILURES ========================================================
_____________________________________________ TestLevy_stable.test_neg_beta ______________________________________________

self = <test_stats.TestLevy_stable object at 0x7fc2f4c94510>

    def test_neg_beta(self):
        # Confirm that _fitstart() gets same estimate of alpha for +beta and -beta.
        rvc.random_state = 2021
        NPOINTS = 1001
        ALPHA = 1.5
        BETA_POS = 1
    
        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, size=NPOINTS)
        # Fit stable distribution to positive and negative mirror image data, examine alpha.
        alpha_pos, alpha_neg = [levy_stable._fitstart(data)[0] for data in [draw, -draw]]
    
>       assert np.isclose(alpha_pos, alpha_neg)
E       assert False
E        +  where False = <function isclose at 0x7fc2f3ccb5f0>(1.5529548477386994, 1.5672455848079985)
E        +    where <function isclose at 0x7fc2f3ccb5f0> = np.isclose

test_stats.py:21: AssertionError
_____________________________________________ TestLevy_stable.test_neg_delta _____________________________________________

self = <test_stats.TestLevy_stable object at 0x7fc2f4c94a10>

    def test_neg_delta(self):
        # Confirm that _fitstart() returns symmetric +delta and -delta if beta flips.
        rvc.random_state = 2021
        NPOINTS = 1001
        ALPHA = 1.5
        BETA_POS = 1
    
        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, size=NPOINTS)
        # Fit stable distribution to positive and negative mirror image data, examine location.
        delta_pos, delta_neg = [levy_stable._fitstart(data)[2] for data in [draw, -draw]]
    
>       assert np.isclose(delta_pos, -delta_neg)
E       assert False
E        +  where False = <function isclose at 0x7fc2f3ccb5f0>(2.220446049250313e-16, -1.0714006456656149)
E        +    where <function isclose at 0x7fc2f3ccb5f0> = np.isclose

test_stats.py:34: AssertionError
____________________________________________ TestLevy_stable.test_delta_shift ____________________________________________

self = <test_stats.TestLevy_stable object at 0x7fc2f4ce2250>

    def test_delta_shift(self):
        # Confirm that _fitstart() returns delta that slides up and down if data shifts.
        rvc.random_state = 2021
        NPOINTS = 1001
        ALPHA = 1.5
        BETA_POS = 1
        LOC = -2
        SHIFT = 1
    
        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, loc=LOC, size=NPOINTS)
        # Fit stable distribution to data and shifted data, examine location parameter.
        delta0, delta_up = [levy_stable._fitstart(draw + x)[2] for x in [0, SHIFT]]
    
>       assert np.isclose(delta0, delta_up-SHIFT)
E       assert False
E        +  where False = <function isclose at 0x7fc2f3ccb5f0>(2.220446049250313e-16, (2.220446049250313e-16 - 1))
E        +    where <function isclose at 0x7fc2f3ccb5f0> = np.isclose

test_stats.py:49: AssertionError
================================================ short test summary info =================================================
FAILED test_stats.py::TestLevy_stable::test_neg_beta - assert False
FAILED test_stats.py::TestLevy_stable::test_neg_delta - assert False
FAILED test_stats.py::TestLevy_stable::test_delta_shift - assert False
============================================== 3 failed, 1 passed in 5.47s ===============================================
```

#### Scipy/Numpy/Python version information:
<!-- You can simply run the following and paste the result in a code block
```
import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
```

-->

```
1.5.2 1.18.5 sys.version_info(major=3, minor=7, micro=6, releaselevel='final', serial=0)
```

Hints from issue discussion:
The issues all occur in the last few lines of _fitstart().  The first three issues are fixed in the code snippet below, though the issue remains with trying to return an arbitrarily large negative location parameter when alpha==1.  I haven't tracked down the relevant reference to determine what value should be returned for that case.

```
            ...
            if nu_alpha >= 2.439:
                # FIX: Use abs(nu_beta).
                alpha = np.clip(psi_1(abs(nu_beta), nu_alpha)[0], np.finfo(float).eps, 2.)
                beta = np.clip(psi_2_1(nu_beta, nu_alpha)[0], -1., 1.)
            else:
                alpha = 2.0
                beta = np.sign(nu_beta)
            # FIX NOT MADE: Could apply clipping to make sure, but nu_c_table guarantees positive values.
            c = (p75 - p25) / phi_3_1(beta, alpha)[0]  # Calculate scale.
            zeta = p50 + c*phi_5_1(beta, alpha)[0]  # Calculate loc (unless alpha==1).
            # FIX (partial?): Removed clipping, since loc (delta) range is (-inf, inf).
            delta = zeta-beta*c*np.tan(np.pi*alpha/2.) if alpha == 1. else zeta
    
            return (alpha, beta, delta, c)
```

Commit message: Update _continuous_distns.py

Fix estimate of alpha when beta < 0, and allow negative values of loc (delta)

Patch to review:
diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py
index 1a970393c6ae..7d779e22a5a7 100644
--- a/scipy/stats/_continuous_distns.py
+++ b/scipy/stats/_continuous_distns.py
@@ -5109,14 +5109,14 @@ def _fitstart(self, data):
         nu_beta = (p95 + p05 - 2*p50)/(p95 - p05)
 
         if nu_alpha >= 2.439:
-            alpha = np.clip(psi_1(nu_beta, nu_alpha)[0], np.finfo(float).eps, 2.)
+            alpha = np.clip(psi_1(abs(nu_beta), nu_alpha)[0], np.finfo(float).eps, 2.)
             beta = np.clip(psi_2_1(nu_beta, nu_alpha)[0], -1., 1.)
         else:
             alpha = 2.0
             beta = np.sign(nu_beta)
         c = (p75 - p25) / phi_3_1(beta, alpha)[0]
         zeta = p50 + c*phi_5_1(beta, alpha)[0]
-        delta = np.clip(zeta-beta*c*np.tan(np.pi*alpha/2.) if alpha == 1. else zeta, np.finfo(float).eps, np.inf)
+        delta = zeta-beta*c*np.tan(np.pi*alpha/2.) if alpha == 1. else zeta
 
         return (alpha, beta, delta, c)
 
diff --git a/scipy/stats/tests/test_levy_stable_fitstart.py b/scipy/stats/tests/test_levy_stable_fitstart.py
new file mode 100644
index 000000000000..59c0e68ea6e0
--- /dev/null
+++ b/scipy/stats/tests/test_levy_stable_fitstart.py
@@ -0,0 +1,91 @@
+# Test scipy.stats.levy_stable._fitstart().
+
+import numpy as np
+import pytest
+
+from scipy.stats import levy_stable
+from scipy.stats import rv_continuous as rvc
+
+class TestLevy_stable:
+    def test_neg_beta(self):
+        # Confirm that _fitstart() gets same estimate of alpha for +beta and -beta.
+        rvc.random_state = 2021
+        NPOINTS = 1001
+        ALPHA = 1.5
+        BETA_POS = 1
+        
+        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, size=NPOINTS)
+        # Fit stable distribution to positive and negative mirror image data, examine alpha.
+        alpha_pos, alpha_neg = [levy_stable._fitstart(data)[0] for data in [draw, -draw]]
+        
+        assert np.isclose(alpha_pos, alpha_neg)
+        
+    def test_neg_delta(self):
+        # Confirm that _fitstart() returns symmetric +delta and -delta if beta flips.
+        rvc.random_state = 2021
+        NPOINTS = 1001
+        ALPHA = 1.5
+        BETA_POS = 1
+        
+        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, size=NPOINTS)
+        # Fit stable distribution to positive and negative mirror image data, examine location.
+        delta_pos, delta_neg = [levy_stable._fitstart(data)[2] for data in [draw, -draw]]
+        
+        assert np.isclose(delta_pos, -delta_neg)
+        
+    def test_delta_shift(self):
+        # Confirm that _fitstart() returns delta that slides up and down if data shifts.
+        rvc.random_state = 2021
+        NPOINTS = 1001
+        ALPHA = 1.5
+        BETA_POS = 1
+        LOC = -2
+        SHIFT = 1
+        
+        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, loc=LOC, size=NPOINTS)
+        # Fit stable distribution to data and shifted data, examine location parameter.
+        delta0, delta_up = [levy_stable._fitstart(draw + x)[2] for x in [0, SHIFT]]
+        
+        assert np.isclose(delta0, delta_up-SHIFT)
+        
+    def test_alpha1(self):
+        # Confirm that stable distribution can fit alpha==1 without blowing up.
+        SHIFT = 0
+        SCALE = 1
+                
+        # Make a dataset whose quintiles lead to alpha==1.
+        x5 = x25 = x50 = 5  # Percentiles 5, 10, 25 -> 5.
+        x95 = 100 - x5  # Percentile 95.
+        # Calculate value for percentile 75 that will interpolate alpha lookup to 1.
+        #   We assume beta==1 so interpolation is between rows 4 and 5 of last column
+        #   of _fitstart() alpha_table.
+        x75 = x25 + (x95 - x5) / (4 + (1.150 - 1) / (1.150 - .973))
+        #assert x5 <= x25 <= x50 <= x75 <= x95
+        
+        xx = np.array([x5, x25, x50, x75, x95])
+        # np.percentile(xx, [5, 25, 50, 75, 95]) -> [5, 5, 5, ..., 95]
+        
+        # Place target percentiles either side.
+        data = np.array([x5, x5,  # Percentiles 0, 10.
+                         x25, x25,  # 20, 30.
+                         x50, x50, x50,  # 40, 50, 60.
+                         x75, x75,  # 70, 80.
+                         x95, x95  # 90, 100.
+                         ])
+        #assert np.array_equal(np.percentile(data, [5, 25, 50, 75, 95]), xx)
+        
+        alpha, beta, loc, scale = fitstart_fix._fitstart_fix(SCALE*data + SHIFT)
+        
+        # Not sure what value of delta should be returned when alpha==1 so
+        #  for now just verify that alpha==1.
+        assert alpha == 1
+        # Might want to check loc too--what value should it have when alpha==1?
+        
+#%%
+
+if __name__ == "__main__":
+    t = TestLevy_stable()
+    t.test_neg_beta()
+    t.test_neg_delta()
+    t.test_delta_shift()
+    t.test_alpha1()


Code review comments:
diff_hunk: @@ -0,0 +1,91 @@
+# Test scipy.stats.levy_stable._fitstart().
+
+import numpy as np
+import pytest
+
+from scipy.stats import levy_stable
+from scipy.stats import rv_continuous as rvc
+
+class TestLevy_stable:
+    def test_neg_beta(self):
+        # Confirm that _fitstart() gets same estimate of alpha for +beta and -beta.
+        rvc.random_state = 2021
+        NPOINTS = 1001
+        ALPHA = 1.5
+        BETA_POS = 1
+        
+        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, size=NPOINTS)
+        # Fit stable distribution to positive and negative mirror image data, examine alpha.
+        alpha_pos, alpha_neg = [levy_stable._fitstart(data)[0] for data in [draw, -draw]]
+        
+        assert np.isclose(alpha_pos, alpha_neg)
path: scipy/stats/tests/test_levy_stable_fitstart.py
text: @user1:
Here and elsewhere, it is likely to be far more helpful to use `assert_allclose()` from the NumPy testing infrastructure than it is to use a plain `assert` on the boolean return value of the plain `isclose()` function, because `assert_allclose()` is specifically designed for testing purposes and clearly reporting the number and magnitude of differences.

@author:
"far more helpful to use `assert_allclose()`"

Thanks.  That would be sound advice for sure in a context that could involve multiple values.  In `_fitstart()` the parameter lookups are assumed to be scalars so there is arguably less gain from the extra output of `assert_allclose()`, and perhaps potential for incorrectly implying that multiple values are potentially involved.  Thoughts?  Would you still advice using `assert_allclose()`?

@author:
I've moved the tests to a better neighbourhood in test_distributions.py around line 3052 and rewritten them more in the style of what was there already.  The tests now mostly use numpy testing functions, except for a few plain asserts involving inequalities--is there a better way to do those?

diff_hunk: @@ -0,0 +1,91 @@
+# Test scipy.stats.levy_stable._fitstart().
+
+import numpy as np
+import pytest
+
+from scipy.stats import levy_stable
+from scipy.stats import rv_continuous as rvc
+
+class TestLevy_stable:
+    def test_neg_beta(self):
+        # Confirm that _fitstart() gets same estimate of alpha for +beta and -beta.
+        rvc.random_state = 2021
+        NPOINTS = 1001
+        ALPHA = 1.5
+        BETA_POS = 1
+        
+        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, size=NPOINTS)
+        # Fit stable distribution to positive and negative mirror image data, examine alpha.
+        alpha_pos, alpha_neg = [levy_stable._fitstart(data)[0] for data in [draw, -draw]]
+        
+        assert np.isclose(alpha_pos, alpha_neg)
+        
+    def test_neg_delta(self):
+        # Confirm that _fitstart() returns symmetric +delta and -delta if beta flips.
+        rvc.random_state = 2021
+        NPOINTS = 1001
+        ALPHA = 1.5
+        BETA_POS = 1
path: scipy/stats/tests/test_levy_stable_fitstart.py
text: @user1:
here and elsewhere, I'm not sure much is gained by defining the variables like this and then using them only once vs. just using the values directly in the keyword arguments

@author:
Fair point, thanks.

diff_hunk: @@ -0,0 +1,91 @@
+# Test scipy.stats.levy_stable._fitstart().
+
+import numpy as np
+import pytest
+
+from scipy.stats import levy_stable
+from scipy.stats import rv_continuous as rvc
+
+class TestLevy_stable:
+    def test_neg_beta(self):
+        # Confirm that _fitstart() gets same estimate of alpha for +beta and -beta.
+        rvc.random_state = 2021
+        NPOINTS = 1001
+        ALPHA = 1.5
+        BETA_POS = 1
+        
+        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, size=NPOINTS)
+        # Fit stable distribution to positive and negative mirror image data, examine alpha.
+        alpha_pos, alpha_neg = [levy_stable._fitstart(data)[0] for data in [draw, -draw]]
+        
+        assert np.isclose(alpha_pos, alpha_neg)
+        
+    def test_neg_delta(self):
+        # Confirm that _fitstart() returns symmetric +delta and -delta if beta flips.
+        rvc.random_state = 2021
+        NPOINTS = 1001
+        ALPHA = 1.5
+        BETA_POS = 1
+        
+        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, size=NPOINTS)
+        # Fit stable distribution to positive and negative mirror image data, examine location.
+        delta_pos, delta_neg = [levy_stable._fitstart(data)[2] for data in [draw, -draw]]
+        
+        assert np.isclose(delta_pos, -delta_neg)
+        
+    def test_delta_shift(self):
+        # Confirm that _fitstart() returns delta that slides up and down if data shifts.
+        rvc.random_state = 2021
+        NPOINTS = 1001
+        ALPHA = 1.5
+        BETA_POS = 1
+        LOC = -2
+        SHIFT = 1
+        
+        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, loc=LOC, size=NPOINTS)
+        # Fit stable distribution to data and shifted data, examine location parameter.
+        delta0, delta_up = [levy_stable._fitstart(draw + x)[2] for x in [0, SHIFT]]
+        
+        assert np.isclose(delta0, delta_up-SHIFT)
+        
+    def test_alpha1(self):
+        # Confirm that stable distribution can fit alpha==1 without blowing up.
+        SHIFT = 0
+        SCALE = 1
+                
+        # Make a dataset whose quintiles lead to alpha==1.
+        x5 = x25 = x50 = 5  # Percentiles 5, 10, 25 -> 5.
+        x95 = 100 - x5  # Percentile 95.
+        # Calculate value for percentile 75 that will interpolate alpha lookup to 1.
+        #   We assume beta==1 so interpolation is between rows 4 and 5 of last column
+        #   of _fitstart() alpha_table.
+        x75 = x25 + (x95 - x5) / (4 + (1.150 - 1) / (1.150 - .973))
+        #assert x5 <= x25 <= x50 <= x75 <= x95
+        
+        xx = np.array([x5, x25, x50, x75, x95])
+        # np.percentile(xx, [5, 25, 50, 75, 95]) -> [5, 5, 5, ..., 95]
+        
+        # Place target percentiles either side.
+        data = np.array([x5, x5,  # Percentiles 0, 10.
+                         x25, x25,  # 20, 30.
+                         x50, x50, x50,  # 40, 50, 60.
+                         x75, x75,  # 70, 80.
+                         x95, x95  # 90, 100.
+                         ])
+        #assert np.array_equal(np.percentile(data, [5, 25, 50, 75, 95]), xx)
+        
+        alpha, beta, loc, scale = fitstart_fix._fitstart_fix(SCALE*data + SHIFT)
+        
+        # Not sure what value of delta should be returned when alpha==1 so
+        #  for now just verify that alpha==1.
+        assert alpha == 1
+        # Might want to check loc too--what value should it have when alpha==1?
+        
+#%%
+
+if __name__ == "__main__":
+    t = TestLevy_stable()
+    t.test_neg_beta()
+    t.test_neg_delta()
+    t.test_delta_shift()
+    t.test_alpha1()
path: scipy/stats/tests/test_levy_stable_fitstart.py
text: @user1:
I don't think we normally add blocks like this--instead, leaving the work to `pytest` to collect and run the tests automatically

@author:
This code block is not included in the new tests now.

diff_hunk: @@ -0,0 +1,91 @@
+# Test scipy.stats.levy_stable._fitstart().
+
+import numpy as np
+import pytest
+
+from scipy.stats import levy_stable
+from scipy.stats import rv_continuous as rvc
+
+class TestLevy_stable:
+    def test_neg_beta(self):
+        # Confirm that _fitstart() gets same estimate of alpha for +beta and -beta.
+        rvc.random_state = 2021
+        NPOINTS = 1001
+        ALPHA = 1.5
+        BETA_POS = 1
+        
+        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, size=NPOINTS)
+        # Fit stable distribution to positive and negative mirror image data, examine alpha.
+        alpha_pos, alpha_neg = [levy_stable._fitstart(data)[0] for data in [draw, -draw]]
+        
+        assert np.isclose(alpha_pos, alpha_neg)
+        
+    def test_neg_delta(self):
+        # Confirm that _fitstart() returns symmetric +delta and -delta if beta flips.
+        rvc.random_state = 2021
+        NPOINTS = 1001
+        ALPHA = 1.5
+        BETA_POS = 1
+        
+        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, size=NPOINTS)
+        # Fit stable distribution to positive and negative mirror image data, examine location.
+        delta_pos, delta_neg = [levy_stable._fitstart(data)[2] for data in [draw, -draw]]
+        
+        assert np.isclose(delta_pos, -delta_neg)
+        
+    def test_delta_shift(self):
+        # Confirm that _fitstart() returns delta that slides up and down if data shifts.
+        rvc.random_state = 2021
+        NPOINTS = 1001
+        ALPHA = 1.5
+        BETA_POS = 1
+        LOC = -2
+        SHIFT = 1
+        
+        draw = levy_stable.rvs(alpha=ALPHA, beta=BETA_POS, loc=LOC, size=NPOINTS)
+        # Fit stable distribution to data and shifted data, examine location parameter.
+        delta0, delta_up = [levy_stable._fitstart(draw + x)[2] for x in [0, SHIFT]]
+        
+        assert np.isclose(delta0, delta_up-SHIFT)
+        
+    def test_alpha1(self):
+        # Confirm that stable distribution can fit alpha==1 without blowing up.
+        SHIFT = 0
+        SCALE = 1
+                
+        # Make a dataset whose quintiles lead to alpha==1.
+        x5 = x25 = x50 = 5  # Percentiles 5, 10, 25 -> 5.
+        x95 = 100 - x5  # Percentile 95.
+        # Calculate value for percentile 75 that will interpolate alpha lookup to 1.
+        #   We assume beta==1 so interpolation is between rows 4 and 5 of last column
+        #   of _fitstart() alpha_table.
+        x75 = x25 + (x95 - x5) / (4 + (1.150 - 1) / (1.150 - .973))
+        #assert x5 <= x25 <= x50 <= x75 <= x95
+        
+        xx = np.array([x5, x25, x50, x75, x95])
+        # np.percentile(xx, [5, 25, 50, 75, 95]) -> [5, 5, 5, ..., 95]
+        
+        # Place target percentiles either side.
+        data = np.array([x5, x5,  # Percentiles 0, 10.
+                         x25, x25,  # 20, 30.
+                         x50, x50, x50,  # 40, 50, 60.
+                         x75, x75,  # 70, 80.
+                         x95, x95  # 90, 100.
+                         ])
+        #assert np.array_equal(np.percentile(data, [5, 25, 50, 75, 95]), xx)
+        
+        alpha, beta, loc, scale = fitstart_fix._fitstart_fix(SCALE*data + SHIFT)
path: scipy/stats/tests/test_levy_stable_fitstart.py
text: @user1:
looks like the CI is complaining about this one: `E   NameError: name 'fitstart_fix' is not defined`

@author:
Yikes!  Clearly juggling too many new tools at once and taking the wrong shortcuts.

diff_hunk: @@ -5109,14 +5109,14 @@ def _fitstart(self, data):
         nu_beta = (p95 + p05 - 2*p50)/(p95 - p05)
 
         if nu_alpha >= 2.439:
-            alpha = np.clip(psi_1(nu_beta, nu_alpha)[0], np.finfo(float).eps, 2.)
+            alpha = np.clip(psi_1(abs(nu_beta), nu_alpha)[0], np.finfo(float).eps, 2.)
path: scipy/stats/_continuous_distns.py
text: @user1:
For consistency, I'd suggest handling this in the same way that the current code does by adding a `psi_1_1` function:
https://github.com/scipy/scipy/blob/affdc4d54142cdcfedcfa06e1185c2ed0441c36c/scipy/stats/_continuous_distns.py#L5092-L5099

McCulloch's paper has the negative `nu_beta` case, but it was somehow missed in scipy's implementation.
![image](https://user-images.githubusercontent.com/14023456/127035848-a775baf3-2b14-48d2-8dc1-ded584d5a6c0.png)

The current code essentially treats negative `nu_beta` as `nu_beta == 0.0` in the interpolation, which is incorrect.


@author:
"I'd suggest handling this in the same way that the current code does by adding a psi_1_1 function"

In terms of consistency, you're probably right.

Any appetite for a slightly more ambitious revision with new helper functions `_rectified()` and `_reflected()` that would result in something like this:

```
def _rectified(func):
    return lambda xval, yval: func(val if xval > 0 else -xval, yval)

def _reflected(func):
    rectified_func = _rectified(func)
    def lookup(xval, yval):
        rectified_result = rectified_func(val,  yval)
        return rectified_result if xval > 0 else -rectified_result
    return lookup

psi_1_1 = _rectified(interpolate.interp2d(nu_beta_range, nu_alpha_range, alpha_table))
psi_2_1 = _reflected(interpolate.interp2d(nu_beta_range, nu_alpha_range, beta_table))

phi_3_1 = _rectified(interpolate.interp2d(beta_range, alpha_range, nu_c_table))
phi_5_1 = _reflected(interpolate.interp2d(beta_range, alpha_range, nu_zeta_table))

```

Or is it better to go with the minimal fix--less risk of breaking something that currently works?

@author:
Thanks for the excerpts from McCulloch--I haven't tracked down a copy yet.  McCulloch's equation 4.6 for the delta estimator suggests that the adjustment to zeta should be applied when alpha is **not** 1.  This is exactly the opposite of what the code is doing (I show it here with spurious clipping removed):

`delta = zeta-beta*c*np.tan(np.pi*alpha/2.) if alpha == 1. else zeta`

So I think the `==` should be `!=`.  Is that right?

If the true value could lie anywhere when alpha is exactly 1, there is arguably no point in checking for this event which is vanishingly unlikely anyway.  So maybe the code should just be:

`delta = zeta-beta*c*np.tan(np.pi*alpha/2.)`

Thoughts?

@user1:
I think you're right, but I'll have to check more closely tonight. There are many different parameterizations of stable laws and the literature has historically been very inconsistent in this regard.

I know scipy's current code uses Nolan's S1 parameterization and much of the theoretical work in papers used to use his S0 parameterization, which requires a shift of exactly the form used in the delta estimate here. It's unlikely, but that *could* explain the difference in the code.

For instance, Nolan notes the S0 <-> S1 transformation requires the following shift
![image](https://user-images.githubusercontent.com/14023456/127168112-758500a3-d06f-474f-ad3d-b788817c47e9.png)
(Also note that the field never really converged on consistent notation. Above, gamma plays a role similar to McCulloch's "c" in your code).

@user1:
It looks like swapping `alpha == 1` to `alpha != 1` is all that's needed for the delta estimation issue. More details are in https://github.com/scipy/scipy/pull/14476#issuecomment-887853522.