Repository: scipy/scipy
File: scipy/stats/_continuous_distns.py (Python)

Code before:
            -0.00000000006514,
            -0.00000000000536,
            +0.00000000000603,
            -0.00000000000233,
            +0.00000000000069,
            -0.00000000000018,
            +0.00000000000003,
        )
        w0 = np.real(sc.lambertw(np.exp(x)))
        # For "very large values of x", sc.lambertw fails, instead the
        # equation provided in the paper (see docstring) works better.
        LIMIT_WHERE_TO_STOP_USING_LAMBERT_EQUATION_SOLVER = 400
        large_x = x[x>LIMIT_WHERE_TO_STOP_USING_LAMBERT_EQUATION_SOLVER]
        w0[x>LIMIT_WHERE_TO_STOP_USING_LAMBERT_EQUATION_SOLVER] = large_x-large_x*np.log(large_x)/(large_x+1)

        w4 = 2.9262710624435
        phi = np.polynomial.chebyshev.Chebyshev(dn)(2*w4/w0-1)/w0**2
        phi[(x < 4)] = float('NaN')
        return phi

    x_is_iterable = hasattr(x, '__iter__')
    if not x_is_iterable:
        x = [x]
    x = np.asarray(x)

    with np.errstate(over='ignore', invalid='ignore', divide='ignore'):
        phi = np.nansum(np.stack([f(x) for f in (phi_0, phi_1, phi_2, phi_3, phi_4)]), axis=0)

    if x_is_iterable:
        return phi
    else:
        return phi[0]


def _landau_cdf(x):
    """Fast calculation of Landau cumulative distribution function.

    Calculation formula from K.S. Kölbig, B. Schorr, A program package for
    the Landau distribution, Computer Physics Communications, Volume 31,
    Issue 1, 1984, Pages 97-111, ISSN 0010-4655, https://doi.org/10.1016/0010-4655(84)90085-7.
    Public preprint freely accessible in CERN Document Server, https://cds.cern.ch/record/145105

    Parameters
    ----------
    x: float, numpy array
        Point in which to calculate the Landau distribution.

    Returns
    -------
    landau_pdf: float, numpy array
        The value of the Landau distribution.

Review comment (performance):
Turns out there's an efficient and relatively simple way to [generate random variates from stable distributions](https://www.jstor.org/stable/pdf/2285309.pdf?casa_token=XMi92GzUp9sAAAAA:VwyR-RlzSUZbnXYlFHJ_H6LlIvxJe_Vk-CSP7ba9uFVkY9zq7ylkRWBuATIhMC7sAhgv9JasHo-FV2QH6sSbR-ELOZdtWTXpdqFzAg7a8OygbUXbijkx). 

For Landau ($\alpha=\beta=1$):
<img width="487" alt="image" src="https://github.com/user-attachments/assets/9f81648c-bb87-47e0-9059-a3053ffd0018">
($W$ is standard exponential.)

```python3
from scipy.special import _ufuncs
import matplotlib.pyplot as plt

import numpy as np
rng = np.random.default_rng(2348923492)
size = 100000
pi_2 = np.pi/2

U = rng.uniform(-np.pi/2, np.pi/2, size=size)
W = rng.standard_exponential(size=size)
S = 2/np.pi * ((pi_2 + U) * np.tan(U) - np.log((pi_2*W*np.cos(U))/(pi_2 + U)))

x = np.linspace(-3, 20, 300)
pdf = _ufuncs._landau_pdf(x, 0, 1)
plt.plot(x, pdf)
plt.hist(S, density=True, bins=x[::5])
plt.show()
```

<img width="785" alt="image" src="https://github.com/user-attachments/assets/f064bedf-52bf-49af-8830-43cb2dc0f44d">

Thanks for finding this @tk-yoshimura!




Code after:
    def _cdf(self, x):
        return scu._landau_cdf(x, 0, 1)

    def _sf(self, x):
        return scu._landau_sf(x, 0, 1)

    def _ppf(self, p):
        return scu._landau_ppf(p, 0, 1)

    def _isf(self, p):
        return scu._landau_isf(p, 0, 1)

    def _stats(self):
        return np.nan, np.nan, np.nan, np.nan

    def _munp(self, n):
        return np.nan

    def _fitstart(self, data, args=None):
        # Initialize ML guesses using quartiles instead of moments.
        if isinstance(data, CensoredData):
            data = data._uncensor()
        p25, p50, p75 = np.percentile(data, [25, 50, 75])
        return p50, (p75 - p25)/2

    def _rvs(self, size=None, random_state=None):
        # Method from https://www.jstor.org/stable/2285309 Eq. 2.4
        pi_2 = np.pi / 2
        U = random_state.uniform(-np.pi / 2, np.pi / 2, size=size)
        W = random_state.standard_exponential(size=size)
        S = 2 / np.pi * ((pi_2 + U) * np.tan(U)
                         - np.log((pi_2 * W * np.cos(U)) / (pi_2 + U)))
        return S


landau = landau_gen(name='landau')


class laplace_gen(rv_continuous):
    r"""A Laplace continuous random variable.

    %(before_notes)s

    Notes
    -----
    The probability density function for `laplace` is

    .. math::

        f(x) = \frac{1}{2} \exp(-|x|)