BUG: `stats.levy_stable.rvs`: "S0"-parameterization ignored about scipy

[Linfye]

@lucascolley I will do it later in my spare time.

[bewit]

@Linfye thank you for your response. Unfortunately, this code still leads to an AssertionError, even using a fresh setup.

I think the frozen property of scipy distributions leads to this behavior. The scipy.stats framework and the current implementation of using differently parametrized stable distributions are incompatible. As far as I found out, this affects any computation of levy_stable, not just the sampling. So everything is computed in "S1" and trying to set the distribution instance or class to "S0" is ignored all along.

I'm going to restate the issue and try to fix it in the future, but right now, I have no resources to do so.

[bewit]

@Linfye I'm sorry, I just didn't see your PR earlier. That should definitely do the trick for now. Thank you!

[fancidev]

From #21332 (comment), it seems the documented way to change the parameterization (by setting levy_stable.parameterization = …) only works if you also call the methods via the levy_stable object.

If you call the methods from a “frozen” distribution object, i.e. one fixed with shape, location, and scale parameters, then you’ll have to modify the dist attribute of the frozen object instead. For example, the following code works:

import numpy as np
from scipy.stats import levy_stable

stable_distribution = levy_stable(alpha=1.9, beta=0.1, loc=0.0, scale=1.0)
stable_distribution.dist.parameterization = "S0"
samples_s0 = stable_distribution.rvs(size=(1000,1), random_state=49)

stable_distribution.dist.parameterization = "S1"
samples_s1 = stable_distribution.rvs(size=(1000,1), random_state=49)

assert not np.all(np.isclose(samples_s0, samples_s1)) # passes without error

[mdhaber]

🤷
Could modify __setattr__ (tricky to get right; adds overhead). Could add a new setter set(attribute_name, value) method that the user must use if they want the attributes to work on frozen distributions (requires users to change behavior). Could document the implementation details (not great practice). Could close as wontfix (doesn't fix the problem). There's a downside to everything. Tough to decide what is "best". But since we're working with existing code that is very complicated, rigid, and fragile, the repetitive solution here might be the least worst. I won't stop it if others decide to go a different direction, though. Just trying to close the issue.

[Linfye]

I think I solve the pronlem. In the Scipy doc, it suggests to use levy_stable.parameterization to set the parameterization because it is class attribute instead of instance attribute. So the code should be like this:

import numpy as np
from scipy.stats import levy_stable

stable_distribution = levy_stable(alpha=1.9, beta=0.1, loc=0.0, scale=1.0)

levy_stable.parameterization = "S0"
samples_s0 = stable_distribution.rvs(size=(1000,1), random_state=49)

levy_stable.parameterization = "S1"
samples_s1 = stable_distribution.rvs(size=(1000,1), random_state=49)

assert not np.all(np.isclose(samples_s0, samples_s1))

In the Scipy source code, it uses wrong self.parameterization in the function so I correct it. And now it works fine.

[Linfye]

@bewit In my output, there's no AssertionError:

And the data in samples_s0 and samples_s1 are indeed different.

Could you please tell me more about your steps of running the modified code?

[bewit]

@Linfye Now I'm really puzzled 😄 Your screenshot looks just like I what I would expect the samples to look like, slightly shifted.
I opened up a fresh Google-Colab Notebook and copy-pasted your code. But here it still throws an AssertionError and the samples are identical. But I can try it again on a local machine (maybe using different versions).

Wondering if this might be due to setup differences? I'm not that knowledgeable about the scipy code base, just working a lot with the alpha-stable distribution lately. If anybody with more experience might have a look into this or could give a hint, I'm happy to do what I can.

[lucascolley]

One way to check is to add the check as a regression test in your PR @Linfye , and see what CI thinks. Mind doing that?

[Linfye]

@bewit You should modify the Scipy source code first like this:

And then paste my code:

import numpy as np
from scipy.stats import levy_stable

stable_distribution = levy_stable(alpha=1.9, beta=0.1, loc=0.0, scale=1.0)

levy_stable.parameterization = "S0"
samples_s0 = stable_distribution.rvs(size=(1000,1), random_state=49)

levy_stable.parameterization = "S1"
samples_s1 = stable_distribution.rvs(size=(1000,1), random_state=49)

assert not np.all(np.isclose(samples_s0, samples_s1))

In my code, I change your code stable_distribution.parameterization to levy_stable.parameterization.

Hope this would help.

[fancidev]

I could make some additions to the doc to clarify the usage of parameterization (and other attributes that modify the numerical aspects). There are three ways to set the attributes:

1. By modifying the class attribute levy_stable_gen.attr. This affects all objects except frozen objects created before that have dist.attr already set.
2. By modifying the instance attribute levy_stable.attr. This only affects methods called via this instance; in particular, it has no effect on frozen objects.
3. By modifying dist.attr of a frozen object. This affects only that frozen object.

On the other hand, I’m not super keen to advertise these “tricks”, in light of the upcoming reworked infrastructure which has native support for different parameterizations and algorithm parameters.

@mdhaber Do you think it’s worthwhile to modify the current doc as above? Or do you prefer to supersede it with the new infrastructure instead?

[mdhaber]

I don't think we need to modify the docs. I think we can just subclass rv_continuous_frozen for levy_stable and give it an @property named parameterization with a setter that modifies self.dist.parameterization. I'll open a PR.

[fancidev]

I don't think we need to modify the docs. I think we can just subclass rv_continuous_frozen for levy_stable and give it an @property named parameterization with a setter that modifies self.dist.parameterization. I'll open a PR.

That would work too. But levy_stable has a bunch of other parameters that work in the same fashion as parameterization, so it looks a bit tedious to also support them via rv_frozen…