BUG: levy_stable _fitstart fixes #14476

lutefiskhotdish · 2021-07-25T17:39:06Z

Reference issue

What does this implement/fix?

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

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.

Test scipy.stats.levy_stable._fitstart() with data for negative beta, negative loc (delta), and alpha==1.

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

tylerjereddy

I'll defer to a stats regular on the actual fix--I just added a few "mechanical" comments re: the tests.

scipy/stats/tests/test_levy_stable_fitstart.py

ragibson · 2021-07-26T14:26:54Z

Note that there is an ongoing overhaul of the levy_stable code at #9523.

This fix is partially correct, but is handled in a somewhat incorrect way (in terms of the current implementation structure).

Note that McCulloch's paper shows how to handle negative values in their interpolations.

Table 1:
Table 2:
Table 3:
Table 4:
Table 5:
Table 6:
Table 7:

Scipy's code is simply missing the negative nu_beta case in psi_1. Look at

# table III - alpha = psi_1(nu_alpha, nu_beta)
# table IV - beta = psi_2(nu_alpha, nu_beta)
# Table V - nu_c = psi_3(alpha, beta)
# Table VII - nu_zeta = psi_5(alpha, beta)

scipy/scipy/stats/_continuous_distns.py

Lines 5092 to 5099 in 88d4d94

    
           psi_1 = interpolate.interp2d(nu_beta_range, nu_alpha_range, alpha_table, kind='linear') 
        
           psi_2 = interpolate.interp2d(nu_beta_range, nu_alpha_range, beta_table, kind='linear') 
        
           psi_2_1 = lambda nu_beta, nu_alpha: psi_2(nu_beta, nu_alpha) if nu_beta > 0 else -psi_2(-nu_beta, nu_alpha) 
        
           phi_3 = interpolate.interp2d(beta_range, alpha_range, nu_c_table, kind='linear') 
        
           phi_3_1 = lambda beta, alpha: phi_3(beta, alpha) if beta > 0 else phi_3(-beta, alpha) 
        
           phi_5 = interpolate.interp2d(beta_range, alpha_range, nu_zeta_table, kind='linear') 
        
           phi_5_1 = lambda beta, alpha: phi_5(beta, alpha) if beta > 0 else -phi_5(-beta, alpha)

Really, I think you should add psi_1_1 similar to the existing definition of psi_2_1 and then use it in place of psi_1 at

scipy/scipy/stats/_continuous_distns.py

Lines 5111 to 5116 in 88d4d94

    
           if nu_alpha >= 2.439: 
        
               alpha = np.clip(psi_1(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)

I think you are also correct that the clipping of delta to [eps, infty) is incorrect here:

scipy/scipy/stats/_continuous_distns.py

Lines 5117 to 5121 in 88d4d94

    
           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) 
        
           return (alpha, beta, delta, c)

Of course, the location of the distribution can be negative and I see no reason why the code does this unless negative delta is handled elsewhere and I missed it while skimming.

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.

This is an unfortunate reality of the distribution itself in most common parameterizations (e.g. Zolotarev's / Nolan's S0). The mode of the distribution inherently diverges to infinity when delta is finite and alpha is close to (but not exactly) 1.0. McCulloch's paper has a few bits where this issue is mentioned, e.g.

Note: "if alpha is insignificantly different from unity, the true value of delta could lie anywhere".

That said, they do mention how to handle the case where alpha happens to be precisely unity and this is handled in the alpha == 1 check at

scipy/scipy/stats/_continuous_distns.py

Line 5119 in 88d4d94

    
           delta = np.clip(zeta-beta*c*np.tan(np.pi*alpha/2.) if alpha == 1. else zeta, np.finfo(float).eps, np.inf)

scipy/stats/_continuous_distns.py

Uses numpy.testing.assert_allclose() instead of generic assert with numpy.isclose(). Removes named constants that were only used once, for the sake of simplicity. Removes one test of dubious utility. Resolves NameError caused by incorrect module reference.

ragibson · 2021-07-27T21:41:21Z

@lutefiskhotdish I looked into this further and here are my suggested changes to the fix itself. I'm trying to keep the changes as minimal as reasonably possible so as to minimize the merging difficulties into the levy_stable overhaul in #9523.

diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py
index 7d779e22a..070738737 100644
--- a/scipy/stats/_continuous_distns.py
+++ b/scipy/stats/_continuous_distns.py
@@ -5090,6 +5090,7 @@ class levy_stable_gen(rv_continuous):
             [0, -0.061, -0.279, -0.659, -1.198]]
 
         psi_1 = interpolate.interp2d(nu_beta_range, nu_alpha_range, alpha_table, kind='linear')
+        psi_1_1 = lambda nu_beta, nu_alpha: psi_1(nu_beta, nu_alpha) if nu_beta > 0 else psi_1(-nu_beta, nu_alpha)
         psi_2 = interpolate.interp2d(nu_beta_range, nu_alpha_range, beta_table, kind='linear')
         psi_2_1 = lambda nu_beta, nu_alpha: psi_2(nu_beta, nu_alpha) if nu_beta > 0 else -psi_2(-nu_beta, nu_alpha)
 
@@ -5109,14 +5110,14 @@ class levy_stable_gen(rv_continuous):
         nu_beta = (p95 + p05 - 2*p50)/(p95 - p05)
 
         if nu_alpha >= 2.439:
-            alpha = np.clip(psi_1(abs(nu_beta), nu_alpha)[0], np.finfo(float).eps, 2.)
+            alpha = np.clip(psi_1_1(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 = zeta-beta*c*np.tan(np.pi*alpha/2.) if alpha == 1. else zeta
+        delta = zeta-beta*c*np.tan(np.pi*alpha/2.) if alpha != 1. else zeta
 
         return (alpha, beta, delta, c)

As mentioned in #14476 (comment), note that the previous code mixes the cases for the delta estimate.

@lutefiskhotdish With the above patch, your tests still pass

$ python runtests.py -v -m full -t scipy/stats/tests/test_levy_stable_fitstart.py 
Build OK (0:06:09.884302 elapsed)
================================================ test session starts =================================================
platform linux -- Python 3.9.6, pytest-6.2.4, py-1.10.0, pluggy-0.13.1 -- /usr/bin/python3.9
cachedir: .pytest_cache
rootdir: /home/ryan/PycharmProjects/scipy_stable/scipy, configfile: pytest.ini
collected 3 items                                                                                                    

scipy/stats/tests/test_levy_stable_fitstart.py::TestLevy_stable::test_neg_beta PASSED                          [ 33%]
scipy/stats/tests/test_levy_stable_fitstart.py::TestLevy_stable::test_neg_delta PASSED                         [ 66%]
scipy/stats/tests/test_levy_stable_fitstart.py::TestLevy_stable::test_delta_shift PASSED                       [100%]

================================================= 3 passed in 0.19s ==================================================

and the estimates are much more accurate

>>> from scipy.stats import levy_stable
>>> 
>>> SIZE = 5 * 10 ** 7
>>> rvs = levy_stable.rvs(alpha=1.75, beta=0.5, loc=2.0, scale=5.0, size=SIZE)
>>> print(levy_stable._fitstart(rvs))
(1.7533970685117124, 0.554005532981971, 2.058881157647261, 4.990065245947207)
>>> 
>>> rvs = levy_stable.rvs(alpha=1.75, beta=-0.5, loc=2.0, scale=5.0, size=SIZE)
>>> print(levy_stable._fitstart(rvs))
(1.7535113451206517, -0.5563912462386369, 1.937892365646908, 4.989910501051077)
>>> 
>>> rvs = levy_stable.rvs(alpha=1.75, beta=0.5, loc=-2.0, scale=5.0, size=SIZE)
>>> print(levy_stable._fitstart(rvs))
(1.7530215468936947, 0.5526121180489312, -1.9417311451832033, 4.989026437213537)
>>> 
>>> rvs = levy_stable.rvs(alpha=1.75, beta=-0.5, loc=-2.0, scale=5.0, size=SIZE)
>>> print(levy_stable._fitstart(rvs))
(1.7529422179586065, -0.5530799346384316, -2.05956703533164, 4.989271645373516)

These tests estimate loc ~ 0.9 in your code and scipy's current implementation.

@lutefiskhotdish You should also move your tests near the existing tests for levy_stable._fitstart() rather than making a completely new file. Somewhere around here:

scipy/scipy/stats/tests/test_distributions.py

Lines 3028 to 3052 in affdc4d

    
           class TestLevyStable: 
        
               def test_fit(self): 
        
                   # construct data to have percentiles that match 
        
                   # example in McCulloch 1986. 
        
                   x = [-.05413, -.05413, 
        
                        0., 0., 0., 0., 
        
                        .00533, .00533, .00533, .00533, .00533, 
        
                        .03354, .03354, .03354, .03354, .03354, 
        
                        .05309, .05309, .05309, .05309, .05309] 
        
                   alpha1, beta1, loc1, scale1 = stats.levy_stable._fitstart(x) 
        
                   assert_allclose(alpha1, 1.48, rtol=0, atol=0.01) 
        
                   assert_almost_equal(beta1, -.22, 2) 
        
                   assert_almost_equal(scale1, 0.01717, 4) 
        
                   # to 2 dps due to rounding error in McCulloch86 
        
                   assert_almost_equal(loc1, 0.00233, 2) 
        
                   # cover alpha=2 scenario 
        
                   x2 = x + [.05309, .05309, .05309, .05309, .05309] 
        
                   alpha2, beta2, loc2, scale2 = stats.levy_stable._fitstart(x2) 
        
                   assert_equal(alpha2, 2) 
        
                   assert_equal(beta2, -1) 
        
                   assert_almost_equal(scale2, .02503, 4) 
        
                   assert_almost_equal(loc2, .03354, 4)

perhaps with more specific names like test_fit_neg_beta(), test_fit_neg_delta(), and test_fit_delta_shift().

The tests in that location set random seeds with np.random.seed() rather than setting rv_continuous.random_state directly, so you'll probably have to change that to match.

I also think you might want to use numpy.testing.assert_almost_equal() for comparing scalars even if assert_allclose() technically works here.

Calculate alpha estimate when beta < 0, in style similar to the other parameters. Adjust delta estimate from zeta, when alpha is *not* 1, fixing apparent typo in comparator.

TestLevyStable.test_fit_beta_flip() confirms that sign of beta affects loc, not alpha or scale. TestLevyStable.test_fit_delta_shift() confirms that loc slides up and down if data shifts. TestLevyStable.test_fit_loc_extrap() confirms that loc goes out of sample for alpha close to 1.

Tests improved and relocated to test_distributions.py

lutefiskhotdish · 2021-07-30T06:12:18Z

@ragibson Thanks. I think I've taken your advice on board in the latest commits but let me know if I've missed anything or if you have additional suggestsions.

ragibson · 2021-07-31T17:01:26Z

@lutefiskhotdish This looks mostly good to me. The CI linter does complain about a few formatting items:

scipy/stats/_continuous_distns.py:5093:9: E731 do not assign a lambda expression, use a def
scipy/stats/_continuous_distns.py:5093:80: E501 line too long (114 > 79 characters)
scipy/stats/_continuous_distns.py:5113:80: E501 line too long (83 > 79 characters)
scipy/stats/tests/test_distributions.py:3052:1: W293 blank line contains whitespace
scipy/stats/tests/test_distributions.py:3064:1: W293 blank line contains whitespace
scipy/stats/tests/test_distributions.py:3075:1: W293 blank line contains whitespace

After this, it looks like all the CI tests will pass, so you might have to go ahead and format the legacy code in _fitstart() according to PEP8 (c.f. https://github.com/scipy/scipy/pull/9523/files/34b9f6f14f3fe6c74200c8074b311b704f31b99f#r670516701).

I wonder if the scipy maintainers would be okay with disabling the lint check for those few lines in _continuous_distns.py just to make it easier to resolve merge conflicts with #9523 (where we have to resolve these PEP8 violations anyway). I suspect the answer is "no", but I admit I don't know exactly how they would handle this situation.

Replace lambda expression with def. Wrap long line.

Removed whitespace on blank lines.

lutefiskhotdish · 2021-07-31T17:47:19Z

@ragibson I think the two commits I've just done make the bug fixes PEP8 compliant. Am I to understand I should fix PEP8 violations in legacy code before BUG fixes can be taken onboard?

ragibson · 2021-07-31T18:23:55Z

@ragibson I think the two commits I've just done make the bug fixes PEP8 compliant. Am I to understand I should fix PEP8 violations in legacy code before BUG fixes can be taken onboard?

@lutefiskhotdish Typically you try to get the CI tests to pass before a PR is merged. Scipy has a tool that checks the diff for PEP8 violations, so editing legacy code (that is not PEP8 compliant) means you inherit those formatting errors from the perspective of the CI tests.

I see a few remaining PEP8 violations from your edits:

scipy/stats/_continuous_distns.py:5094:40: W291 trailing whitespace
scipy/stats/_continuous_distns.py:5095:80: E501 line too long (89 > 79 characters)
scipy/stats/_continuous_distns.py:5116:59: W291 trailing whitespace

For reference, this is from the CI failure at https://dev.azure.com/scipy-org/SciPy/_build/results?buildId=13326&view=logs&j=fdd89725-f657-5341-7d77-3c5c44c4f5b6&t=325fc03e-e0b9-5339-adaf-099a49a14600 and you can run similar checks locally with e.g.

$ python runtests.py -v --pep8

which gives me the same three warnings. There are some details on setting up your local environment correctly for this, but if you're interested, see https://scipy.github.io/devdocs/dev/contributor/runtests.html and https://scipy.github.io/devdocs/dev/contributor/pep8.html.

Following scipy#9523. Replaced lambda expressions with def. Wrapped long lines. Removed trailing whitespace.

lutefiskhotdish · 2021-08-01T09:51:38Z

@ragibson Showing successful checks now so that must be progress, and now I know how to follow the trail of obscure links to find failing lint output....

mdhaber · 2022-02-20T23:20:10Z

@ragibson Was this one close enough to try to resolve the merge conflicts?

ragibson · 2022-02-21T00:04:26Z

@ragibson Was this one close enough to try to resolve the merge conflicts?

Yes, the edits to fitstart() in scipy/stats/_continuous_distns.py here should map directly to

scipy/scipy/stats/_levy_stable/__init__.py

Lines 488 to 629 in a574519

    
           def _fitstart_S1(data): 
        
               # We follow McCullock 1986 method - Simple Consistent Estimators 
        
               # of Stable Distribution Parameters 
        
               # fmt: off 
        
               # Table III and IV 
        
               nu_alpha_range = [2.439, 2.5, 2.6, 2.7, 2.8, 3, 3.2, 3.5, 4, 
        
                                 5, 6, 8, 10, 15, 25] 
        
               nu_beta_range = [0, 0.1, 0.2, 0.3, 0.5, 0.7, 1] 
        
               # table III - alpha = psi_1(nu_alpha, nu_beta) 
        
               alpha_table = [ 
        
                   [2.000, 2.000, 2.000, 2.000, 2.000, 2.000, 2.000], 
        
                   [1.916, 1.924, 1.924, 1.924, 1.924, 1.924, 1.924], 
        
                   [1.808, 1.813, 1.829, 1.829, 1.829, 1.829, 1.829], 
        
                   [1.729, 1.730, 1.737, 1.745, 1.745, 1.745, 1.745], 
        
                   [1.664, 1.663, 1.663, 1.668, 1.676, 1.676, 1.676], 
        
                   [1.563, 1.560, 1.553, 1.548, 1.547, 1.547, 1.547], 
        
                   [1.484, 1.480, 1.471, 1.460, 1.448, 1.438, 1.438], 
        
                   [1.391, 1.386, 1.378, 1.364, 1.337, 1.318, 1.318], 
        
                   [1.279, 1.273, 1.266, 1.250, 1.210, 1.184, 1.150], 
        
                   [1.128, 1.121, 1.114, 1.101, 1.067, 1.027, 0.973], 
        
                   [1.029, 1.021, 1.014, 1.004, 0.974, 0.935, 0.874], 
        
                   [0.896, 0.892, 0.884, 0.883, 0.855, 0.823, 0.769], 
        
                   [0.818, 0.812, 0.806, 0.801, 0.780, 0.756, 0.691], 
        
                   [0.698, 0.695, 0.692, 0.689, 0.676, 0.656, 0.597], 
        
                   [0.593, 0.590, 0.588, 0.586, 0.579, 0.563, 0.513]] 
        
               # table IV - beta = psi_2(nu_alpha, nu_beta) 
        
               beta_table = [ 
        
                   [0, 2.160, 1.000, 1.000, 1.000, 1.000, 1.000], 
        
                   [0, 1.592, 3.390, 1.000, 1.000, 1.000, 1.000], 
        
                   [0, 0.759, 1.800, 1.000, 1.000, 1.000, 1.000], 
        
                   [0, 0.482, 1.048, 1.694, 1.000, 1.000, 1.000], 
        
                   [0, 0.360, 0.760, 1.232, 2.229, 1.000, 1.000], 
        
                   [0, 0.253, 0.518, 0.823, 1.575, 1.000, 1.000], 
        
                   [0, 0.203, 0.410, 0.632, 1.244, 1.906, 1.000], 
        
                   [0, 0.165, 0.332, 0.499, 0.943, 1.560, 1.000], 
        
                   [0, 0.136, 0.271, 0.404, 0.689, 1.230, 2.195], 
        
                   [0, 0.109, 0.216, 0.323, 0.539, 0.827, 1.917], 
        
                   [0, 0.096, 0.190, 0.284, 0.472, 0.693, 1.759], 
        
                   [0, 0.082, 0.163, 0.243, 0.412, 0.601, 1.596], 
        
                   [0, 0.074, 0.147, 0.220, 0.377, 0.546, 1.482], 
        
                   [0, 0.064, 0.128, 0.191, 0.330, 0.478, 1.362], 
        
                   [0, 0.056, 0.112, 0.167, 0.285, 0.428, 1.274]] 
        
               # Table V and VII 
        
               alpha_range = [2, 1.9, 1.8, 1.7, 1.6, 1.5, 1.4, 1.3, 1.2, 1.1, 
        
                              1, 0.9, 0.8, 0.7, 0.6, 0.5] 
        
               beta_range = [0, 0.25, 0.5, 0.75, 1] 
        
               # Table V - nu_c = psi_3(alpha, beta) 
        
               nu_c_table = [ 
        
                   [1.908, 1.908, 1.908, 1.908, 1.908], 
        
                   [1.914, 1.915, 1.916, 1.918, 1.921], 
        
                   [1.921, 1.922, 1.927, 1.936, 1.947], 
        
                   [1.927, 1.930, 1.943, 1.961, 1.987], 
        
                   [1.933, 1.940, 1.962, 1.997, 2.043], 
        
                   [1.939, 1.952, 1.988, 2.045, 2.116], 
        
                   [1.946, 1.967, 2.022, 2.106, 2.211], 
        
                   [1.955, 1.984, 2.067, 2.188, 2.333], 
        
                   [1.965, 2.007, 2.125, 2.294, 2.491], 
        
                   [1.980, 2.040, 2.205, 2.435, 2.696], 
        
                   [2.000, 2.085, 2.311, 2.624, 2.973], 
        
                   [2.040, 2.149, 2.461, 2.886, 3.356], 
        
                   [2.098, 2.244, 2.676, 3.265, 3.912], 
        
                   [2.189, 2.392, 3.004, 3.844, 4.775], 
        
                   [2.337, 2.634, 3.542, 4.808, 6.247], 
        
                   [2.588, 3.073, 4.534, 6.636, 9.144]] 
        
               # Table VII - nu_zeta = psi_5(alpha, beta) 
        
               nu_zeta_table = [ 
        
                   [0, 0.000, 0.000, 0.000, 0.000], 
        
                   [0, -0.017, -0.032, -0.049, -0.064], 
        
                   [0, -0.030, -0.061, -0.092, -0.123], 
        
                   [0, -0.043, -0.088, -0.132, -0.179], 
        
                   [0, -0.056, -0.111, -0.170, -0.232], 
        
                   [0, -0.066, -0.134, -0.206, -0.283], 
        
                   [0, -0.075, -0.154, -0.241, -0.335], 
        
                   [0, -0.084, -0.173, -0.276, -0.390], 
        
                   [0, -0.090, -0.192, -0.310, -0.447], 
        
                   [0, -0.095, -0.208, -0.346, -0.508], 
        
                   [0, -0.098, -0.223, -0.380, -0.576], 
        
                   [0, -0.099, -0.237, -0.424, -0.652], 
        
                   [0, -0.096, -0.250, -0.469, -0.742], 
        
                   [0, -0.089, -0.262, -0.520, -0.853], 
        
                   [0, -0.078, -0.272, -0.581, -0.997], 
        
                   [0, -0.061, -0.279, -0.659, -1.198]] 
        
               # fmt: on 
        
               psi_1 = interpolate.interp2d( 
        
                   nu_beta_range, nu_alpha_range, alpha_table, kind="linear" 
        
               ) 
        
               psi_2 = interpolate.interp2d( 
        
                   nu_beta_range, nu_alpha_range, beta_table, kind="linear" 
        
               ) 
        
               def psi_2_1(nu_beta, nu_alpha): 
        
                   return psi_2(nu_beta, nu_alpha) \ 
        
                       if nu_beta > 0 else -psi_2(-nu_beta, nu_alpha) 
        
               phi_3 = interpolate.interp2d( 
        
                   beta_range, alpha_range, nu_c_table, kind="linear" 
        
               ) 
        
               def phi_3_1(beta, alpha): 
        
                   return phi_3(beta, alpha) if beta > 0 else phi_3(-beta, alpha) 
        
               phi_5 = interpolate.interp2d( 
        
                   beta_range, alpha_range, nu_zeta_table, kind="linear" 
        
               ) 
        
               def phi_5_1(beta, alpha): 
        
                   return phi_5(beta, alpha) if beta > 0 else -phi_5(-beta, alpha) 
        
               # quantiles 
        
               p05 = np.percentile(data, 5) 
        
               p50 = np.percentile(data, 50) 
        
               p95 = np.percentile(data, 95) 
        
               p25 = np.percentile(data, 25) 
        
               p75 = np.percentile(data, 75) 
        
               nu_alpha = (p95 - p05) / (p75 - p25) 
        
               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.0) 
        
                   beta = np.clip(psi_2_1(nu_beta, nu_alpha)[0], -1.0, 1.0) 
        
               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.0) 
        
                   if alpha == 1.0 
        
                   else zeta, 
        
                   np.finfo(float).eps, 
        
                   np.inf, 
        
               ) 
        
               return (alpha, beta, delta, c)

and the edits in scipy/stats/tests/test_distributions.py will probably have no conflicts since they can still be placed right after levy_stable's test_fit()

scipy/scipy/stats/tests/test_distributions.py

Lines 3293 to 3295 in 7818647

    
           assert_almost_equal(scale2, .02503, 4) 
        
           assert_almost_equal(loc2, .03354, 4)

…t_fixes

mdhaber · 2022-02-21T04:02:14Z

Hmmm. Many of the changes in this PR were already present. You can compare by going back to this commit and comparing against main. Maybe PEP8 made you make some of these changes when you moved to a new file?

I can make fixes if need be, but I'll let @steppi merge this one if it looks good.

steppi · 2022-02-21T05:07:51Z

I can make fixes if need be, but I'll let @steppi merge this one if it looks good.

Thanks. I should have time to take a look tomorrow.

This is scipy#14476 and the bug was not fixed.

steppi

This looks good to me. I pushed one commit to change the xfail reason for test_fit_rvs. I had also written the previous reason; it referred to this PR as a possible fix and was written before I looked into what this PR is actually doing. The test failure is an unrelated timeout. @ragibson, if you have no objections, I'll merge this tomorrow.

ragibson

Yes, this looks good to me. The main fixes were addressing the missing negative sign in psi_1 and removing the clipping of delta to [eps, inf) since that parameter can definitely be negative.

The extra tests seem reasonable as well if they're passing.

mdhaber · 2022-03-03T14:43:47Z

Thanks @lutefiskhotdish and reviewers!

lutefiskhotdish added 2 commits July 25, 2021 17:02

Create test_levy_stable_fitstart.py

7677eb3

Test scipy.stats.levy_stable._fitstart() with data for negative beta, negative loc (delta), and alpha==1.

Update _continuous_distns.py

affdc4d

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

tylerjereddy added scipy.stats defect A clear bug or issue that prevents SciPy from being installed or used as expected labels Jul 25, 2021

tylerjereddy reviewed Jul 25, 2021

View reviewed changes

scipy/stats/tests/test_levy_stable_fitstart.py Outdated Show resolved Hide resolved

scipy/stats/tests/test_levy_stable_fitstart.py Outdated Show resolved Hide resolved

scipy/stats/tests/test_levy_stable_fitstart.py Outdated Show resolved Hide resolved

tylerjereddy reviewed Jul 25, 2021

View reviewed changes

scipy/stats/tests/test_levy_stable_fitstart.py Outdated Show resolved Hide resolved

tupui linked an issue Jul 26, 2021 that may be closed by this pull request

levy_stable._fitstart() incorrect parameter estimates for negative beta, negative loc, or alpha==1 #14470

Closed

ragibson reviewed Jul 26, 2021

View reviewed changes

scipy/stats/_continuous_distns.py Outdated Show resolved Hide resolved

lutefiskhotdish added 2 commits July 26, 2021 19:47

Merge branch 'scipy:master' into levy_stable_fitstart_fixes

918bffb

lutefiskhotdish added 4 commits July 28, 2021 06:59

BUG: Fix levy_stable._fitstart()

c98b1a6

Calculate alpha estimate when beta < 0, in style similar to the other parameters. Adjust delta estimate from zeta, when alpha is *not* 1, fixing apparent typo in comparator.

Merge branch 'scipy:master' into levy_stable_fitstart_fixes

8e70105

TST: delete test_levy_stable_fitstart.py

772f5f3

Tests improved and relocated to test_distributions.py

ragibson mentioned this pull request Jul 31, 2021

ENH: improvements to the Stable distribution #9523

Merged

lutefiskhotdish added 2 commits July 31, 2021 18:27

STY: Resolve linter complaints

71ba6f9

Replace lambda expression with def. Wrap long line.

STY: Resolve linter complaints

51b6006

Removed whitespace on blank lines.

lutefiskhotdish added 3 commits August 1, 2021 08:45

Merge branch 'scipy:master' into levy_stable_fitstart_fixes

2f6dbeb

STY: Resolve PEP8 violations

87a3ec0

Following scipy#9523. Replaced lambda expressions with def. Wrapped long lines. Removed trailing whitespace.

STY: Wrap long lines

724f2d1

ragibson mentioned this pull request Dec 7, 2021

MAINT: Wrapping some things up bsdz/scipy#11

Merged

Merge remote-tracking branch 'upstream/main' into levy_stable_fitstar…

df8b034

…t_fixes

steppi pushed a commit to steppi/scipy that referenced this pull request Feb 22, 2022

TST: Update reason for test failure

fcf1ca8

This is scipy#14476 and the bug was not fixed.

TST: Update reason for test failure

bea64a0

This is scipy#14476 and the bug was not fixed.

steppi force-pushed the levy_stable_fitstart_fixes branch from fcf1ca8 to bea64a0 Compare February 22, 2022 02:25

steppi approved these changes Feb 22, 2022

View reviewed changes

ragibson approved these changes Feb 22, 2022

View reviewed changes

steppi merged commit 4a856f1 into scipy:main Feb 22, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

BUG: levy_stable _fitstart fixes #14476

BUG: levy_stable _fitstart fixes #14476

lutefiskhotdish commented Jul 25, 2021

tylerjereddy left a comment

ragibson commented Jul 26, 2021 •

edited

Loading

ragibson commented Jul 27, 2021 •

edited

Loading

lutefiskhotdish commented Jul 30, 2021

ragibson commented Jul 31, 2021

lutefiskhotdish commented Jul 31, 2021 •

edited

Loading

ragibson commented Jul 31, 2021

lutefiskhotdish commented Aug 1, 2021

mdhaber commented Feb 20, 2022

ragibson commented Feb 21, 2022 •

edited

Loading

mdhaber commented Feb 21, 2022

steppi commented Feb 21, 2022

steppi left a comment

ragibson left a comment

mdhaber commented Mar 3, 2022

BUG: levy_stable _fitstart fixes #14476

BUG: levy_stable _fitstart fixes #14476

Conversation

lutefiskhotdish commented Jul 25, 2021

Reference issue

What does this implement/fix?

Additional information

tylerjereddy left a comment

Choose a reason for hiding this comment

ragibson commented Jul 26, 2021 • edited Loading

ragibson commented Jul 27, 2021 • edited Loading

lutefiskhotdish commented Jul 30, 2021

ragibson commented Jul 31, 2021

lutefiskhotdish commented Jul 31, 2021 • edited Loading

ragibson commented Jul 31, 2021

lutefiskhotdish commented Aug 1, 2021

mdhaber commented Feb 20, 2022

ragibson commented Feb 21, 2022 • edited Loading

mdhaber commented Feb 21, 2022

steppi commented Feb 21, 2022

steppi left a comment

Choose a reason for hiding this comment

ragibson left a comment

Choose a reason for hiding this comment

mdhaber commented Mar 3, 2022

ragibson commented Jul 26, 2021 •

edited

Loading

ragibson commented Jul 27, 2021 •

edited

Loading

lutefiskhotdish commented Jul 31, 2021 •

edited

Loading

ragibson commented Feb 21, 2022 •

edited

Loading