Repository: scipy/scipy
Pull request: BUG: `stats.levy_stable`: fix ignored parameterization
<!-- 
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:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message.
However, please only include an issue number in the description, not the title,
and please ensure that any code names containing underscores are enclosed in backticks.

Depending on your changes, you can skip CI operations and save time and energy: 
http://scipy.github.io/devdocs/dev/contributor/continuous_integration.html#skipping

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.-->
Solves #20821
#### What does this implement/fix?
<!--Please explain your changes.-->
I modified the parameters used in the _parameterization function of the levy_stable_gen class. Previously, this function did not read the passed parameters because it used class attributes. By changing it to use instance attributes, this bug can be resolved.
#### Additional information
<!--Any additional information you think is important.-->
This is an example:
```
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))
```

Problem statement:
BUG: `stats.levy_stable.rvs`: "S0"-parameterization ignored
### Describe your issue.

When trying to sample from stats.levy_stable.rvs() with levy_stable.parametrization="S0", the returned samples are not transformed, but in the "S1"-parameterization. Line 850 of stats._levy_stable.__init__ does not get executed.
This is likely due to automatically resolving to the frozen distribution in stats._distn_infrastructure:rv_frozen.rvs() [line 493]. 


### Reproducing Code Example

```python
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.parameterization = "S0" #  we could also set this directly as levy_stable.parameterization, makes no difference
samples_s0 = stable_distribution.rvs(size=(1000,1), random_state=49)

stable_distribution.parameterization = "S1"
samples_s1 = stable_distribution.rvs(size=(1000,1), random_state=49)

# due to the needed transformation between "S0" and "S1" parameterizations, the samples should definitely not be identical, but they are
assert not np.all(np.isclose(samples_s0, samples_s1)) # this should hold!
```


### Error message

```shell
No error message
```


### SciPy/NumPy/Python version and system information

<details>

```shell
{
  "Compilers": {
    "c": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "cc"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "3.0.10",
      "commands": "cython"
    },
    "c++": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "c++"
    },
    "fortran": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "gfortran"
    },
    "pythran": {
      "version": "0.15.0",
      "include directory": "C:\\Users\\runneradmin\\AppData\\Local\\Temp\\pip-build-env-o281jlrt\\overlay\\Lib\\site-packages/pythran"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "windows"
    },
    "build": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "windows"
    },
    "cross-compiled": false
  },
  "Build Dependencies": {
    "blas": {
      "name": "openblas",
      "found": true,
      "version": "0.3.27",
      "detection method": "pkgconfig",
      "include directory": "/c/opt/64/include",
      "lib directory": "/c/opt/64/lib",
      "openblas configuration": "USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=24",
      "pc file directory": "c:/opt/64/lib/pkgconfig"
    },
    "lapack": {
      "name": "openblas",
      "found": true,
      "version": "0.3.27",
      "detection method": "pkgconfig",
      "include directory": "/c/opt/64/include",
      "lib directory": "/c/opt/64/lib",
      "openblas configuration": "USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=24",
      "pc file directory": "c:/opt/64/lib/pkgconfig"
    },
    "pybind11": {
      "name": "pybind11",
      "version": "2.12.0",
      "detection method": "config-tool",
      "include directory": "unknown"
    }
  },
  "Python Information": {
    "path": "C:\\Users\\runneradmin\\AppData\\Local\\Temp\\cibw-run-622t44xc\\cp39-win_amd64\\build\\venv\\Scripts\\python.exe",
    "version": "3.9"
  }
}
```

</details>

Hints from issue discussion:
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 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 In my output, there's no AssertionError:
![8c04b19421a163bac79c7963bb3c505](https://github.com/user-attachments/assets/1cef1d2a-6fe1-4e78-bd06-1993f05d6fd6)
And the data in samples_s0 and samples_s1 are indeed different.
![947c59bc58899b8da8c0184bec6e1f9](https://github.com/user-attachments/assets/bdd57c64-7015-42e4-b3f6-3a74b5048d40)
Could you please tell me more about your steps of running the modified code?
@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.
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?
@bewit You should modify the Scipy source code **first** like this:
![PixPin_2024-08-07_22-50-45](https://github.com/user-attachments/assets/7a49b187-19e9-47a4-b3c5-8bfd30395bc0)
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.
@lucascolley I will do it later in my spare time.
@Linfye I'm sorry, I just didn't see your PR earlier. That should definitely do the trick for now. Thank you!
From https://github.com/scipy/scipy/pull/21332#issuecomment-2278882307, 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:

```Python
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
```
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](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.levy_stable.html) as above? Or do you prefer to supersede it with the new infrastructure instead?

Commit message: MAINT: stats.levy_stable: adjust approach

Patch to review:
diff --git a/scipy/stats/_levy_stable/__init__.py b/scipy/stats/_levy_stable/__init__.py
index 13d143e3dd56..9803dbf609e5 100644
--- a/scipy/stats/_levy_stable/__init__.py
+++ b/scipy/stats/_levy_stable/__init__.py
@@ -12,7 +12,7 @@
 from scipy.interpolate import RectBivariateSpline
 import scipy.special as sc
 from scipy._lib._util import _lazywhere
-from .._distn_infrastructure import rv_continuous, _ShapeInfo
+from .._distn_infrastructure import rv_continuous, _ShapeInfo, rv_continuous_frozen
 from .._continuous_distns import uniform, expon, _norm_pdf, _norm_cdf
 from .levyst import Nolan
 from scipy._lib.doccer import inherit_docstring_from
@@ -816,6 +816,11 @@ class levy_stable_gen(rv_continuous):
     pdf_fft_interpolation_level = 3
     pdf_fft_interpolation_degree = 3
 
+    def freeze(self, *args, **params):
+        dist = levy_stable_frozen(self, *args, **params)
+        dist.parameterization = self.parameterization
+        return dist
+
     def _argcheck(self, alpha, beta):
         return (alpha > 0) & (alpha <= 2) & (beta <= 1) & (beta >= -1)
 
@@ -1222,3 +1227,13 @@ def pdf_from_cf_with_fft(cf, h=0.01, q=9, level=3):
 
 
 levy_stable = levy_stable_gen(name="levy_stable")
+
+
+class levy_stable_frozen(rv_continuous_frozen):
+    @property
+    def parameterization(self):
+        return self.dist.parameterization
+
+    @parameterization.setter
+    def parameterization(self, value):
+        self.dist.parameterization = value
diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py
index 61235b81b27c..cfe0ecf1341e 100644
--- a/scipy/stats/tests/test_distributions.py
+++ b/scipy/stats/tests/test_distributions.py
@@ -5864,6 +5864,47 @@ def test_x_near_zeta(
             expected,
         )
 
+    def test_frozen_parameterization_gh20821(self):
+        # gh-20821 reported that frozen distributions ignore the parameterization.
+        # Check that this is resolved and that the frozen distribution's
+        # parameterization can be changed independently of stats.levy_stable
+        rng = np.random.default_rng
+        shapes = dict(alpha=1.9, beta=0.1, loc=0.0, scale=1.0)
+        unfrozen = stats.levy_stable
+        frozen = stats.levy_stable(**shapes)
+
+        unfrozen.parameterization = "S0"
+        frozen.parameterization = "S1"
+        unfrozen_a = unfrozen.rvs(**shapes, size=10, random_state=rng(329823498))
+        frozen_a = frozen.rvs(size=10, random_state=rng(329823498))
+        assert not np.any(frozen_a == unfrozen_a)
+
+        unfrozen.parameterization = "S1"
+        frozen.parameterization = "S0"
+        unfrozen_b = unfrozen.rvs(**shapes, size=10, random_state=rng(329823498))
+        frozen_b = frozen.rvs(size=10, random_state=rng(329823498))
+        assert_equal(frozen_b, unfrozen_a)
+        assert_equal(unfrozen_b, frozen_a)
+
+    def test_frozen_parameterization_gh20821b(self):
+        # Check that the parameterization of the frozen distribution is that of
+        # the unfrozen distribution at the time of freezing
+        rng = np.random.default_rng
+        shapes = dict(alpha=1.9, beta=0.1, loc=0.0, scale=1.0)
+        unfrozen = stats.levy_stable
+
+        unfrozen.parameterization = "S0"
+        frozen = stats.levy_stable(**shapes)
+        unfrozen_a = unfrozen.rvs(**shapes, size=10, random_state=rng(329823498))
+        frozen_a = frozen.rvs(size=10, random_state=rng(329823498))
+        assert_equal(frozen_a, unfrozen_a)
+
+        unfrozen.parameterization = "S1"
+        frozen = stats.levy_stable(**shapes)
+        unfrozen_b = unfrozen.rvs(**shapes, size=10, random_state=rng(329823498))
+        frozen_b = frozen.rvs(size=10, random_state=rng(329823498))
+        assert_equal(frozen_b, unfrozen_b)
+
 
 class TestArrayArgument:  # test for ticket:992
     def setup_method(self):


Code review comments:
diff_hunk: @@ -816,6 +816,11 @@ class levy_stable_gen(rv_continuous):
     pdf_fft_interpolation_level = 3
     pdf_fft_interpolation_degree = 3
 
+    def freeze(self, *args, **params):
path: scipy/stats/_levy_stable/__init__.py
text: @user1:
```suggestion
    def __call__(self, *args, **params):
```
This should fix the refguide check failure.