__init__.py

#

import warnings
from functools import partial

import numpy as np

from scipy import optimize
from scipy import integrate
from scipy.integrate._quadrature import _builtincoeffs
from scipy import interpolate
from scipy.interpolate import RectBivariateSpline
import scipy.special as sc
from scipy._lib._util import _lazywhere
from .._distn_infrastructure import rv_continuous, _ShapeInfo
from .._continuous_distns import uniform, expon, _norm_pdf, _norm_cdf
from .levyst import Nolan
from scipy._lib.doccer import inherit_docstring_from


__all__ = ["levy_stable", "levy_stable_gen", "pdf_from_cf_with_fft"]

# Stable distributions are known for various parameterisations
# some being advantageous for numerical considerations and others
# useful due to their location/scale awareness.
#
# Here we follow [NO] convention (see the references in the docstring
# for levy_stable_gen below).
#
# S0 / Z0 / x0 (aka Zoleterav's M)
# S1 / Z1 / x1
#
# Where S* denotes parameterisation, Z* denotes standardized
# version where gamma = 1, delta = 0 and x* denotes variable.
#
# Scipy's original Stable was a random variate generator. It
# uses S1 and unfortunately is not a location/scale aware.


# default numerical integration tolerance
# used for epsrel in piecewise and both epsrel and epsabs in dni
# (epsabs needed in dni since weighted quad requires epsabs > 0)
_QUAD_EPS = 1.2e-14


def _Phi_Z0(alpha, t):
    return (
        -np.tan(np.pi * alpha / 2) * (np.abs(t) ** (1 - alpha) - 1)
        if alpha != 1
        else -2.0 * np.log(np.abs(t)) / np.pi
    )


def _Phi_Z1(alpha, t):
    return (
        np.tan(np.pi * alpha / 2)
        if alpha != 1
        else -2.0 * np.log(np.abs(t)) / np.pi
    )


def _cf(Phi, t, alpha, beta):
    """Characteristic function."""
    return np.exp(
        -(np.abs(t) ** alpha) * (1 - 1j * beta * np.sign(t) * Phi(alpha, t))
    )


_cf_Z0 = partial(_cf, _Phi_Z0)
_cf_Z1 = partial(_cf, _Phi_Z1)


def _pdf_single_value_cf_integrate(Phi, x, alpha, beta, **kwds):
    """To improve DNI accuracy convert characteristic function in to real
    valued integral using Euler's formula, then exploit cosine symmetry to
    change limits to [0, inf). Finally use cosine addition formula to split
    into two parts that can be handled by weighted quad pack.
    """
    quad_eps = kwds.get("quad_eps", _QUAD_EPS)

    def integrand1(t):
        if t == 0:
            return 0
        return np.exp(-(t ** alpha)) * (
            np.cos(beta * (t ** alpha) * Phi(alpha, t))
        )

    def integrand2(t):
        if t == 0:
            return 0
        return np.exp(-(t ** alpha)) * (
            np.sin(beta * (t ** alpha) * Phi(alpha, t))
        )

    with np.errstate(invalid="ignore"):
        int1, *ret1 = integrate.quad(
            integrand1,
            0,
            np.inf,
            weight="cos",
            wvar=x,
            limit=1000,
            epsabs=quad_eps,
            epsrel=quad_eps,
            full_output=1,
        )

        int2, *ret2 = integrate.quad(
            integrand2,
            0,
            np.inf,
            weight="sin",
            wvar=x,
            limit=1000,
            epsabs=quad_eps,
            epsrel=quad_eps,
            full_output=1,
        )

    return (int1 + int2) / np.pi


_pdf_single_value_cf_integrate_Z0 = partial(
    _pdf_single_value_cf_integrate, _Phi_Z0
)
_pdf_single_value_cf_integrate_Z1 = partial(
    _pdf_single_value_cf_integrate, _Phi_Z1
)


def _nolan_round_x_near_zeta(x0, alpha, zeta, x_tol_near_zeta):
    """Round x close to zeta for Nolan's method in [NO]."""
    #   "8. When |x0-beta*tan(pi*alpha/2)| is small, the
    #   computations of the density and cumulative have numerical problems.
    #   The program works around this by setting
    #   z = beta*tan(pi*alpha/2) when
    #   |z-beta*tan(pi*alpha/2)| < tol(5)*alpha**(1/alpha).
    #   (The bound on the right is ad hoc, to get reasonable behavior
    #   when alpha is small)."
    # where tol(5) = 0.5e-2 by default.
    #
    # We seem to have partially addressed this through re-expression of
    # g(theta) here, but it still needs to be used in some extreme cases.
    # Perhaps tol(5) = 0.5e-2 could be reduced for our implementation.
    if np.abs(x0 - zeta) < x_tol_near_zeta * alpha ** (1 / alpha):
        x0 = zeta
    return x0


def _nolan_round_difficult_input(
    x0, alpha, beta, zeta, x_tol_near_zeta, alpha_tol_near_one
):
    """Round difficult input values for Nolan's method in [NO]."""

    # following Nolan's STABLE,
    #   "1. When 0 < |alpha-1| < 0.005, the program has numerical problems
    #   evaluating the pdf and cdf.  The current version of the program sets
    #   alpha=1 in these cases. This approximation is not bad in the S0
    #   parameterization."
    if np.abs(alpha - 1) < alpha_tol_near_one:
        alpha = 1.0

    #   "2. When alpha=1 and |beta| < 0.005, the program has numerical
    #   problems.  The current version sets beta=0."
    # We seem to have addressed this through re-expression of g(theta) here

    x0 = _nolan_round_x_near_zeta(x0, alpha, zeta, x_tol_near_zeta)
    return x0, alpha, beta


def _pdf_single_value_piecewise_Z1(x, alpha, beta, **kwds):
    # convert from Nolan's S_1 (aka S) to S_0 (aka Zolaterev M)
    # parameterization

    zeta = -beta * np.tan(np.pi * alpha / 2.0)
    x0 = x + zeta if alpha != 1 else x

    return _pdf_single_value_piecewise_Z0(x0, alpha, beta, **kwds)


def _pdf_single_value_piecewise_Z0(x0, alpha, beta, **kwds):

    quad_eps = kwds.get("quad_eps", _QUAD_EPS)
    x_tol_near_zeta = kwds.get("piecewise_x_tol_near_zeta", 0.005)
    alpha_tol_near_one = kwds.get("piecewise_alpha_tol_near_one", 0.005)

    zeta = -beta * np.tan(np.pi * alpha / 2.0)
    x0, alpha, beta = _nolan_round_difficult_input(
        x0, alpha, beta, zeta, x_tol_near_zeta, alpha_tol_near_one
    )

    # some other known distribution pdfs / analytical cases
    # TODO: add more where possible with test coverage,
    # eg https://en.wikipedia.org/wiki/Stable_distribution#Other_analytic_cases
    if alpha == 2.0:
        # normal
        return _norm_pdf(x0 / np.sqrt(2)) / np.sqrt(2)
    elif alpha == 0.5 and beta == 1.0:
        # levy
        # since S(1/2, 1, gamma, delta; <x>) ==
        # S(1/2, 1, gamma, gamma + delta; <x0>).
        _x = x0 + 1
        if _x <= 0:
            return 0

        return 1 / np.sqrt(2 * np.pi * _x) / _x * np.exp(-1 / (2 * _x))
    elif alpha == 0.5 and beta == 0.0 and x0 != 0:
        # analytical solution [HO]
        S, C = sc.fresnel([1 / np.sqrt(2 * np.pi * np.abs(x0))])
        arg = 1 / (4 * np.abs(x0))
        return (
            np.sin(arg) * (0.5 - S[0]) + np.cos(arg) * (0.5 - C[0])
        ) / np.sqrt(2 * np.pi * np.abs(x0) ** 3)
    elif alpha == 1.0 and beta == 0.0:
        # cauchy
        return 1 / (1 + x0 ** 2) / np.pi

    return _pdf_single_value_piecewise_post_rounding_Z0(
        x0, alpha, beta, quad_eps, x_tol_near_zeta
    )


def _pdf_single_value_piecewise_post_rounding_Z0(x0, alpha, beta, quad_eps,
                                                 x_tol_near_zeta):
    """Calculate pdf using Nolan's methods as detailed in [NO]."""

    _nolan = Nolan(alpha, beta, x0)
    zeta = _nolan.zeta
    xi = _nolan.xi
    c2 = _nolan.c2
    g = _nolan.g

    # round x0 to zeta again if needed. zeta was recomputed and may have
    # changed due to floating point differences.
    # See https://github.com/scipy/scipy/pull/18133
    x0 = _nolan_round_x_near_zeta(x0, alpha, zeta, x_tol_near_zeta)
    # handle Nolan's initial case logic
    if x0 == zeta:
        return (
            sc.gamma(1 + 1 / alpha)
            * np.cos(xi)
            / np.pi
            / ((1 + zeta ** 2) ** (1 / alpha / 2))
        )
    elif x0 < zeta:
        return _pdf_single_value_piecewise_post_rounding_Z0(
            -x0, alpha, -beta, quad_eps, x_tol_near_zeta
        )

    # following Nolan, we may now assume
    #   x0 > zeta when alpha != 1
    #   beta != 0 when alpha == 1

    # spare calculating integral on null set
    # use isclose as macos has fp differences
    if np.isclose(-xi, np.pi / 2, rtol=1e-014, atol=1e-014):
        return 0.0

    def integrand(theta):
        # limit any numerical issues leading to g_1 < 0 near theta limits
        g_1 = g(theta)
        if not np.isfinite(g_1) or g_1 < 0:
            g_1 = 0
        return g_1 * np.exp(-g_1)

    with np.errstate(all="ignore"):
        peak = optimize.bisect(
            lambda t: g(t) - 1, -xi, np.pi / 2, xtol=quad_eps
        )

        # this integrand can be very peaked, so we need to force
        # QUADPACK to evaluate the function inside its support
        #

        # lastly, we add additional samples at
        #   ~exp(-100), ~exp(-10), ~exp(-5), ~exp(-1)
        # to improve QUADPACK's detection of rapidly descending tail behavior
        # (this choice is fairly ad hoc)
        tail_points = [
            optimize.bisect(lambda t: g(t) - exp_height, -xi, np.pi / 2)
            for exp_height in [100, 10, 5]
            # exp_height = 1 is handled by peak
        ]
        intg_points = [0, peak] + tail_points
        intg, *ret = integrate.quad(
            integrand,
            -xi,
            np.pi / 2,
            points=intg_points,
            limit=100,
            epsrel=quad_eps,
            epsabs=0,
            full_output=1,
        )

    return c2 * intg


def _cdf_single_value_piecewise_Z1(x, alpha, beta, **kwds):
    # convert from Nolan's S_1 (aka S) to S_0 (aka Zolaterev M)
    # parameterization

    zeta = -beta * np.tan(np.pi * alpha / 2.0)
    x0 = x + zeta if alpha != 1 else x

    return _cdf_single_value_piecewise_Z0(x0, alpha, beta, **kwds)


def _cdf_single_value_piecewise_Z0(x0, alpha, beta, **kwds):

    quad_eps = kwds.get("quad_eps", _QUAD_EPS)
    x_tol_near_zeta = kwds.get("piecewise_x_tol_near_zeta", 0.005)
    alpha_tol_near_one = kwds.get("piecewise_alpha_tol_near_one", 0.005)

    zeta = -beta * np.tan(np.pi * alpha / 2.0)
    x0, alpha, beta = _nolan_round_difficult_input(
        x0, alpha, beta, zeta, x_tol_near_zeta, alpha_tol_near_one
    )

    # some other known distribution cdfs / analytical cases
    # TODO: add more where possible with test coverage,
    # eg https://en.wikipedia.org/wiki/Stable_distribution#Other_analytic_cases
    if alpha == 2.0:
        # normal
        return _norm_cdf(x0 / np.sqrt(2))
    elif alpha == 0.5 and beta == 1.0:
        # levy
        # since S(1/2, 1, gamma, delta; <x>) ==
        # S(1/2, 1, gamma, gamma + delta; <x0>).
        _x = x0 + 1
        if _x <= 0:
            return 0

        return sc.erfc(np.sqrt(0.5 / _x))
    elif alpha == 1.0 and beta == 0.0:
        # cauchy
        return 0.5 + np.arctan(x0) / np.pi

    return _cdf_single_value_piecewise_post_rounding_Z0(
        x0, alpha, beta, quad_eps, x_tol_near_zeta
    )


def _cdf_single_value_piecewise_post_rounding_Z0(x0, alpha, beta, quad_eps,
                                                 x_tol_near_zeta):
    """Calculate cdf using Nolan's methods as detailed in [NO]."""
    _nolan = Nolan(alpha, beta, x0)
    zeta = _nolan.zeta
    xi = _nolan.xi
    c1 = _nolan.c1
    # c2 = _nolan.c2
    c3 = _nolan.c3
    g = _nolan.g
    # round x0 to zeta again if needed. zeta was recomputed and may have
    # changed due to floating point differences.
    # See https://github.com/scipy/scipy/pull/18133
    x0 = _nolan_round_x_near_zeta(x0, alpha, zeta, x_tol_near_zeta)
    # handle Nolan's initial case logic
    if (alpha == 1 and beta < 0) or x0 < zeta:
        # NOTE: Nolan's paper has a typo here!
        # He states F(x) = 1 - F(x, alpha, -beta), but this is clearly
        # incorrect since F(-infty) would be 1.0 in this case
        # Indeed, the alpha != 1, x0 < zeta case is correct here.
        return 1 - _cdf_single_value_piecewise_post_rounding_Z0(
            -x0, alpha, -beta, quad_eps, x_tol_near_zeta
        )
    elif x0 == zeta:
        return 0.5 - xi / np.pi

    # following Nolan, we may now assume
    #   x0 > zeta when alpha != 1
    #   beta > 0 when alpha == 1

    # spare calculating integral on null set
    # use isclose as macos has fp differences
    if np.isclose(-xi, np.pi / 2, rtol=1e-014, atol=1e-014):
        return c1

    def integrand(theta):
        g_1 = g(theta)
        return np.exp(-g_1)

    with np.errstate(all="ignore"):
        # shrink supports where required
        left_support = -xi
        right_support = np.pi / 2
        if alpha > 1:
            # integrand(t) monotonic 0 to 1
            if integrand(-xi) != 0.0:
                res = optimize.minimize(
                    integrand,
                    (-xi,),
                    method="L-BFGS-B",
                    bounds=[(-xi, np.pi / 2)],
                )
                left_support = res.x[0]
        else:
            # integrand(t) monotonic 1 to 0
            if integrand(np.pi / 2) != 0.0:
                res = optimize.minimize(
                    integrand,
                    (np.pi / 2,),
                    method="L-BFGS-B",
                    bounds=[(-xi, np.pi / 2)],
                )
                right_support = res.x[0]

        intg, *ret = integrate.quad(
            integrand,
            left_support,
            right_support,
            points=[left_support, right_support],
            limit=100,
            epsrel=quad_eps,
            epsabs=0,
            full_output=1,
        )

    return c1 + c3 * intg


def _rvs_Z1(alpha, beta, size=None, random_state=None):
    """Simulate random variables using Nolan's methods as detailed in [NO].
    """

    def alpha1func(alpha, beta, TH, aTH, bTH, cosTH, tanTH, W):
        return (
            2
            / np.pi
            * (
                (np.pi / 2 + bTH) * tanTH
                - beta * np.log((np.pi / 2 * W * cosTH) / (np.pi / 2 + bTH))
            )
        )

    def beta0func(alpha, beta, TH, aTH, bTH, cosTH, tanTH, W):
        return (
            W
            / (cosTH / np.tan(aTH) + np.sin(TH))
            * ((np.cos(aTH) + np.sin(aTH) * tanTH) / W) ** (1.0 / alpha)
        )

    def otherwise(alpha, beta, TH, aTH, bTH, cosTH, tanTH, W):
        # alpha is not 1 and beta is not 0
        val0 = beta * np.tan(np.pi * alpha / 2)
        th0 = np.arctan(val0) / alpha
        val3 = W / (cosTH / np.tan(alpha * (th0 + TH)) + np.sin(TH))
        res3 = val3 * (
            (
                np.cos(aTH)
                + np.sin(aTH) * tanTH
                - val0 * (np.sin(aTH) - np.cos(aTH) * tanTH)
            )
            / W
        ) ** (1.0 / alpha)
        return res3

    def alphanot1func(alpha, beta, TH, aTH, bTH, cosTH, tanTH, W):
        res = _lazywhere(
            beta == 0,
            (alpha, beta, TH, aTH, bTH, cosTH, tanTH, W),
            beta0func,
            f2=otherwise,
        )
        return res

    alpha = np.broadcast_to(alpha, size)
    beta = np.broadcast_to(beta, size)
    TH = uniform.rvs(
        loc=-np.pi / 2.0, scale=np.pi, size=size, random_state=random_state
    )
    W = expon.rvs(size=size, random_state=random_state)
    aTH = alpha * TH
    bTH = beta * TH
    cosTH = np.cos(TH)
    tanTH = np.tan(TH)
    res = _lazywhere(
        alpha == 1,
        (alpha, beta, TH, aTH, bTH, cosTH, tanTH, W),
        alpha1func,
        f2=alphanot1func,
    )
    return res


def _fitstart_S0(data):
    alpha, beta, delta1, gamma = _fitstart_S1(data)

    # Formulas for mapping parameters in S1 parameterization to
    # those in S0 parameterization can be found in [NO]. Note that
    # only delta changes.
    if alpha != 1:
        delta0 = delta1 + beta * gamma * np.tan(np.pi * alpha / 2.0)
    else:
        delta0 = delta1 + 2 * beta * gamma * np.log(gamma) / np.pi

    return alpha, beta, delta0, gamma


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 = np.array([
        [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]]).T
    # transpose because interpolation with `RectBivariateSpline` is with
    # `nu_beta` as `x` and `nu_alpha` as `y`

    # table IV - beta = psi_2(nu_alpha, nu_beta)
    beta_table = np.array([
        [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]]).T

    # Table V and VII
    # These are ordered with decreasing `alpha_range`; so we will need to
    # reverse them as required by RectBivariateSpline.
    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][::-1]
    beta_range = [0, 0.25, 0.5, 0.75, 1]

    # Table V - nu_c = psi_3(alpha, beta)
    nu_c_table = np.array([
        [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]])[::-1].T
    # transpose because interpolation with `RectBivariateSpline` is with
    # `beta` as `x` and `alpha` as `y`

    # Table VII - nu_zeta = psi_5(alpha, beta)
    nu_zeta_table = np.array([
        [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]])[::-1].T
    # fmt: on

    psi_1 = RectBivariateSpline(nu_beta_range, nu_alpha_range,
                                alpha_table, kx=1, ky=1, s=0)

    def psi_1_1(nu_beta, nu_alpha):
        return psi_1(nu_beta, nu_alpha) \
            if nu_beta > 0 else psi_1(-nu_beta, nu_alpha)

    psi_2 = RectBivariateSpline(nu_beta_range, nu_alpha_range,
                                beta_table, kx=1, ky=1, s=0)

    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 = RectBivariateSpline(beta_range, alpha_range, nu_c_table,
                                kx=1, ky=1, s=0)

    def phi_3_1(beta, alpha):
        return phi_3(beta, alpha) if beta > 0 else phi_3(-beta, alpha)

    phi_5 = RectBivariateSpline(beta_range, alpha_range, nu_zeta_table,
                                kx=1, ky=1, s=0)

    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:
        eps = np.finfo(float).eps
        alpha = np.clip(psi_1_1(nu_beta, nu_alpha)[0, 0], eps, 2.)
        beta = np.clip(psi_2_1(nu_beta, nu_alpha)[0, 0], -1.0, 1.0)
    else:
        alpha = 2.0
        beta = np.sign(nu_beta)
    c = (p75 - p25) / phi_3_1(beta, alpha)[0, 0]
    zeta = p50 + c * phi_5_1(beta, alpha)[0, 0]
    delta = zeta-beta*c*np.tan(np.pi*alpha/2.) if alpha != 1. else zeta

    return (alpha, beta, delta, c)


class levy_stable_gen(rv_continuous):
    r"""A Levy-stable continuous random variable.

    %(before_notes)s

    See Also
    --------
    levy, levy_l, cauchy, norm

    Notes
    -----
    The distribution for `levy_stable` has characteristic function:

    .. math::

        \varphi(t, \alpha, \beta, c, \mu) =
        e^{it\mu -|ct|^{\alpha}(1-i\beta\operatorname{sign}(t)\Phi(\alpha, t))}

    where two different parameterizations are supported. The first :math:`S_1`:

    .. math::

        \Phi = \begin{cases}
                \tan \left({\frac {\pi \alpha }{2}}\right)&\alpha \neq 1\\
                -{\frac {2}{\pi }}\log |t|&\alpha =1
                \end{cases}

    The second :math:`S_0`:

    .. math::

        \Phi = \begin{cases}
                -\tan \left({\frac {\pi \alpha }{2}}\right)(|ct|^{1-\alpha}-1)
                &\alpha \neq 1\\
                -{\frac {2}{\pi }}\log |ct|&\alpha =1
                \end{cases}


    The probability density function for `levy_stable` is:

    .. math::

        f(x) = \frac{1}{2\pi}\int_{-\infty}^\infty \varphi(t)e^{-ixt}\,dt

    where :math:`-\infty < t < \infty`. This integral does not have a known
    closed form.

    `levy_stable` generalizes several distributions.  Where possible, they
    should be used instead.  Specifically, when the shape parameters
    assume the values in the table below, the corresponding equivalent
    distribution should be used.

    =========  ========  ===========
    ``alpha``  ``beta``   Equivalent
    =========  ========  ===========
     1/2       -1        `levy_l`
     1/2       1         `levy`
     1         0         `cauchy`
     2         any       `norm` (with ``scale=sqrt(2)``)
    =========  ========  ===========

    Evaluation of the pdf uses Nolan's piecewise integration approach with the
    Zolotarev :math:`M` parameterization by default. There is also the option
    to use direct numerical integration of the standard parameterization of the
    characteristic function or to evaluate by taking the FFT of the
    characteristic function.

    The default method can changed by setting the class variable
    ``levy_stable.pdf_default_method`` to one of 'piecewise' for Nolan's
    approach, 'dni' for direct numerical integration, or 'fft-simpson' for the
    FFT based approach. For the sake of backwards compatibility, the methods
    'best' and 'zolotarev' are equivalent to 'piecewise' and the method
    'quadrature' is equivalent to 'dni'.

    The parameterization can be changed  by setting the class variable
    ``levy_stable.parameterization`` to either 'S0' or 'S1'.
    The default is 'S1'.

    To improve performance of piecewise and direct numerical integration one
    can specify ``levy_stable.quad_eps`` (defaults to 1.2e-14). This is used
    as both the absolute and relative quadrature tolerance for direct numerical
    integration and as the relative quadrature tolerance for the piecewise
    method. One can also specify ``levy_stable.piecewise_x_tol_near_zeta``
    (defaults to 0.005) for how close x is to zeta before it is considered the
    same as x [NO]. The exact check is
    ``abs(x0 - zeta) < piecewise_x_tol_near_zeta*alpha**(1/alpha)``. One can
    also specify ``levy_stable.piecewise_alpha_tol_near_one`` (defaults to
    0.005) for how close alpha is to 1 before being considered equal to 1.

    To increase accuracy of FFT calculation one can specify
    ``levy_stable.pdf_fft_grid_spacing`` (defaults to 0.001) and
    ``pdf_fft_n_points_two_power`` (defaults to None which means a value is
    calculated that sufficiently covers the input range).

    Further control over FFT calculation is available by setting
    ``pdf_fft_interpolation_degree`` (defaults to 3) for spline order and
    ``pdf_fft_interpolation_level`` for determining the number of points to use
    in the Newton-Cotes formula when approximating the characteristic function
    (considered experimental).

    Evaluation of the cdf uses Nolan's piecewise integration approach with the
    Zolatarev :math:`S_0` parameterization by default. There is also the option
    to evaluate through integration of an interpolated spline of the pdf
    calculated by means of the FFT method. The settings affecting FFT
    calculation are the same as for pdf calculation. The default cdf method can
    be changed by setting ``levy_stable.cdf_default_method`` to either
    'piecewise' or 'fft-simpson'.  For cdf calculations the Zolatarev method is
    superior in accuracy, so FFT is disabled by default.

    Fitting estimate uses quantile estimation method in [MC]. MLE estimation of
    parameters in fit method uses this quantile estimate initially. Note that
    MLE doesn't always converge if using FFT for pdf calculations; this will be
    the case if alpha <= 1 where the FFT approach doesn't give good
    approximations.

    Any non-missing value for the attribute
    ``levy_stable.pdf_fft_min_points_threshold`` will set
    ``levy_stable.pdf_default_method`` to 'fft-simpson' if a valid
    default method is not otherwise set.



    .. warning::

        For pdf calculations FFT calculation is considered experimental.

        For cdf calculations FFT calculation is considered experimental. Use
        Zolatarev's method instead (default).

    The probability density above is defined in the "standardized" form. To
    shift and/or scale the distribution use the ``loc`` and ``scale``
    parameters.
    Generally ``%(name)s.pdf(x, %(shapes)s, loc, scale)`` is identically
    equivalent to ``%(name)s.pdf(y, %(shapes)s) / scale`` with
    ``y = (x - loc) / scale``, except in the ``S1`` parameterization if
    ``alpha == 1``.  In that case ``%(name)s.pdf(x, %(shapes)s, loc, scale)``
    is identically equivalent to ``%(name)s.pdf(y, %(shapes)s) / scale`` with
    ``y = (x - loc - 2 * beta * scale * np.log(scale) / np.pi) / scale``.
    See [NO2]_ Definition 1.8 for more information.
    Note that shifting the location of a distribution
    does not make it a "noncentral" distribution.

    References
    ----------
    .. [MC] McCulloch, J., 1986. Simple consistent estimators of stable
        distribution parameters. Communications in Statistics - Simulation and
        Computation 15, 11091136.
    .. [WZ] Wang, Li and Zhang, Ji-Hong, 2008. Simpson's rule based FFT method
        to compute densities of stable distribution.
    .. [NO] Nolan, J., 1997. Numerical Calculation of Stable Densities and
        distributions Functions.
    .. [NO2] Nolan, J., 2018. Stable Distributions: Models for Heavy Tailed
        Data.
    .. [HO] Hopcraft, K. I., Jakeman, E., Tanner, R. M. J., 1999. Lévy random
        walks with fluctuating step number and multiscale behavior.

    %(example)s

    """
    # Configurable options as class variables
    # (accessible from self by attribute lookup).
    parameterization = "S1"
    pdf_default_method = "piecewise"
    cdf_default_method = "piecewise"
    quad_eps = _QUAD_EPS
    piecewise_x_tol_near_zeta = 0.005
    piecewise_alpha_tol_near_one = 0.005
    pdf_fft_min_points_threshold = None
    pdf_fft_grid_spacing = 0.001
    pdf_fft_n_points_two_power = None
    pdf_fft_interpolation_level = 3
    pdf_fft_interpolation_degree = 3

    def _argcheck(self, alpha, beta):
        return (alpha > 0) & (alpha <= 2) & (beta <= 1) & (beta >= -1)

    def _shape_info(self):
        ialpha = _ShapeInfo("alpha", False, (0, 2), (False, True))
        ibeta = _ShapeInfo("beta", False, (-1, 1), (True, True))
        return [ialpha, ibeta]

    def _parameterization(self):
        allowed = ("S0", "S1")
        pz = self.parameterization
        if pz not in allowed:
            raise RuntimeError(
                f"Parameterization '{pz}' in supported list: {allowed}"
            )
        return pz

    @inherit_docstring_from(rv_continuous)
    def rvs(self, *args, **kwds):
        X1 = super().rvs(*args, **kwds)

        kwds.pop("discrete", None)
        kwds.pop("random_state", None)
        (alpha, beta), delta, gamma, size = self._parse_args_rvs(*args, **kwds)

        # shift location for this parameterisation (S1)
        X1 = np.where(
            alpha == 1.0, X1 + 2 * beta * gamma * np.log(gamma) / np.pi, X1
        )

        if self._parameterization() == "S0":
            return np.where(
                alpha == 1.0,
                X1 - (beta * 2 * gamma * np.log(gamma) / np.pi),
                X1 - gamma * beta * np.tan(np.pi * alpha / 2.0),
            )
        elif self._parameterization() == "S1":
            return X1

    def _rvs(self, alpha, beta, size=None, random_state=None):
        return _rvs_Z1(alpha, beta, size, random_state)

    @inherit_docstring_from(rv_continuous)
    def pdf(self, x, *args, **kwds):
        # override base class version to correct
        # location for S1 parameterization
        if self._parameterization() == "S0":
            return super().pdf(x, *args, **kwds)
        elif self._parameterization() == "S1":
            (alpha, beta), delta, gamma = self._parse_args(*args, **kwds)
            if np.all(np.reshape(alpha, (1, -1))[0, :] != 1):
                return super().pdf(x, *args, **kwds)
            else:
                # correct location for this parameterisation
                x = np.reshape(x, (1, -1))[0, :]
                x, alpha, beta = np.broadcast_arrays(x, alpha, beta)

                data_in = np.dstack((x, alpha, beta))[0]
                data_out = np.empty(shape=(len(data_in), 1))
                # group data in unique arrays of alpha, beta pairs
                uniq_param_pairs = np.unique(data_in[:, 1:], axis=0)
                for pair in uniq_param_pairs:
                    _alpha, _beta = pair
                    _delta = (
                        delta + 2 * _beta * gamma * np.log(gamma) / np.pi
                        if _alpha == 1.0
                        else delta
                    )
                    data_mask = np.all(data_in[:, 1:] == pair, axis=-1)
                    _x = data_in[data_mask, 0]
                    data_out[data_mask] = (
                        super()
                        .pdf(_x, _alpha, _beta, loc=_delta, scale=gamma)
                        .reshape(len(_x), 1)
                    )
                output = data_out.T[0]
                if output.shape == (1,):
                    return output[0]
                return output

    def _pdf(self, x, alpha, beta):
        if self._parameterization() == "S0":
            _pdf_single_value_piecewise = _pdf_single_value_piecewise_Z0
            _pdf_single_value_cf_integrate = _pdf_single_value_cf_integrate_Z0
            _cf = _cf_Z0
        elif self._parameterization() == "S1":
            _pdf_single_value_piecewise = _pdf_single_value_piecewise_Z1
            _pdf_single_value_cf_integrate = _pdf_single_value_cf_integrate_Z1
            _cf = _cf_Z1

        x = np.asarray(x).reshape(1, -1)[0, :]

        x, alpha, beta = np.broadcast_arrays(x, alpha, beta)

        data_in = np.dstack((x, alpha, beta))[0]
        data_out = np.empty(shape=(len(data_in), 1))

        pdf_default_method_name = self.pdf_default_method
        if pdf_default_method_name in ("piecewise", "best", "zolotarev"):
            pdf_single_value_method = _pdf_single_value_piecewise
        elif pdf_default_method_name in ("dni", "quadrature"):
            pdf_single_value_method = _pdf_single_value_cf_integrate
        elif (
            pdf_default_method_name == "fft-simpson"
            or self.pdf_fft_min_points_threshold is not None
        ):
            pdf_single_value_method = None

        pdf_single_value_kwds = {
            "quad_eps": self.quad_eps,
            "piecewise_x_tol_near_zeta": self.piecewise_x_tol_near_zeta,
            "piecewise_alpha_tol_near_one": self.piecewise_alpha_tol_near_one,
        }

        fft_grid_spacing = self.pdf_fft_grid_spacing
        fft_n_points_two_power = self.pdf_fft_n_points_two_power
        fft_interpolation_level = self.pdf_fft_interpolation_level
        fft_interpolation_degree = self.pdf_fft_interpolation_degree

        # group data in unique arrays of alpha, beta pairs
        uniq_param_pairs = np.unique(data_in[:, 1:], axis=0)
        for pair in uniq_param_pairs:
            data_mask = np.all(data_in[:, 1:] == pair, axis=-1)
            data_subset = data_in[data_mask]
            if pdf_single_value_method is not None:
                data_out[data_mask] = np.array(
                    [
                        pdf_single_value_method(
                            _x, _alpha, _beta, **pdf_single_value_kwds
                        )
                        for _x, _alpha, _beta in data_subset
                    ]
                ).reshape(len(data_subset), 1)
            else:
                warnings.warn(
                    "Density calculations experimental for FFT method."
                    + " Use combination of piecewise and dni methods instead.",
                    RuntimeWarning, stacklevel=3,
                )
                _alpha, _beta = pair
                _x = data_subset[:, (0,)]

                if _alpha < 1.0:
                    raise RuntimeError(
                        "FFT method does not work well for alpha less than 1."
                    )

                # need enough points to "cover" _x for interpolation
                if fft_grid_spacing is None and fft_n_points_two_power is None:
                    raise ValueError(
                        "One of fft_grid_spacing or fft_n_points_two_power "
                        + "needs to be set."
                    )
                max_abs_x = np.max(np.abs(_x))
                h = (
                    2 ** (3 - fft_n_points_two_power) * max_abs_x
                    if fft_grid_spacing is None
                    else fft_grid_spacing
                )
                q = (
                    np.ceil(np.log(2 * max_abs_x / h) / np.log(2)) + 2
                    if fft_n_points_two_power is None
                    else int(fft_n_points_two_power)
                )

                # for some parameters, the range of x can be quite
                # large, let's choose an arbitrary cut off (8GB) to save on
                # computer memory.
                MAX_Q = 30
                if q > MAX_Q:
                    raise RuntimeError(
                        "fft_n_points_two_power has a maximum "
                        + f"value of {MAX_Q}"
                    )

                density_x, density = pdf_from_cf_with_fft(
                    lambda t: _cf(t, _alpha, _beta),
                    h=h,
                    q=q,
                    level=fft_interpolation_level,
                )
                f = interpolate.InterpolatedUnivariateSpline(
                    density_x, np.real(density), k=fft_interpolation_degree
                )  # patch FFT to use cubic
                data_out[data_mask] = f(_x)

        return data_out.T[0]

    @inherit_docstring_from(rv_continuous)
    def cdf(self, x, *args, **kwds):
        # override base class version to correct
        # location for S1 parameterization
        # NOTE: this is near identical to pdf() above
        if self._parameterization() == "S0":
            return super().cdf(x, *args, **kwds)
        elif self._parameterization() == "S1":
            (alpha, beta), delta, gamma = self._parse_args(*args, **kwds)
            if np.all(np.reshape(alpha, (1, -1))[0, :] != 1):
                return super().cdf(x, *args, **kwds)
            else:
                # correct location for this parameterisation
                x = np.reshape(x, (1, -1))[0, :]
                x, alpha, beta = np.broadcast_arrays(x, alpha, beta)

                data_in = np.dstack((x, alpha, beta))[0]
                data_out = np.empty(shape=(len(data_in), 1))
                # group data in unique arrays of alpha, beta pairs
                uniq_param_pairs = np.unique(data_in[:, 1:], axis=0)
                for pair in uniq_param_pairs:
                    _alpha, _beta = pair
                    _delta = (
                        delta + 2 * _beta * gamma * np.log(gamma) / np.pi
                        if _alpha == 1.0
                        else delta
                    )
                    data_mask = np.all(data_in[:, 1:] == pair, axis=-1)
                    _x = data_in[data_mask, 0]
                    data_out[data_mask] = (
                        super()
                        .cdf(_x, _alpha, _beta, loc=_delta, scale=gamma)
                        .reshape(len(_x), 1)
                    )
                output = data_out.T[0]
                if output.shape == (1,):
                    return output[0]
                return output

    def _cdf(self, x, alpha, beta):
        if self._parameterization() == "S0":
            _cdf_single_value_piecewise = _cdf_single_value_piecewise_Z0
            _cf = _cf_Z0
        elif self._parameterization() == "S1":
            _cdf_single_value_piecewise = _cdf_single_value_piecewise_Z1
            _cf = _cf_Z1

        x = np.asarray(x).reshape(1, -1)[0, :]

        x, alpha, beta = np.broadcast_arrays(x, alpha, beta)

        data_in = np.dstack((x, alpha, beta))[0]
        data_out = np.empty(shape=(len(data_in), 1))

        cdf_default_method_name = self.cdf_default_method
        if cdf_default_method_name == "piecewise":
            cdf_single_value_method = _cdf_single_value_piecewise
        elif cdf_default_method_name == "fft-simpson":
            cdf_single_value_method = None

        cdf_single_value_kwds = {
            "quad_eps": self.quad_eps,
            "piecewise_x_tol_near_zeta": self.piecewise_x_tol_near_zeta,
            "piecewise_alpha_tol_near_one": self.piecewise_alpha_tol_near_one,
        }

        fft_grid_spacing = self.pdf_fft_grid_spacing
        fft_n_points_two_power = self.pdf_fft_n_points_two_power
        fft_interpolation_level = self.pdf_fft_interpolation_level
        fft_interpolation_degree = self.pdf_fft_interpolation_degree

        # group data in unique arrays of alpha, beta pairs
        uniq_param_pairs = np.unique(data_in[:, 1:], axis=0)
        for pair in uniq_param_pairs:
            data_mask = np.all(data_in[:, 1:] == pair, axis=-1)
            data_subset = data_in[data_mask]
            if cdf_single_value_method is not None:
                data_out[data_mask] = np.array(
                    [
                        cdf_single_value_method(
                            _x, _alpha, _beta, **cdf_single_value_kwds
                        )
                        for _x, _alpha, _beta in data_subset
                    ]
                ).reshape(len(data_subset), 1)
            else:
                warnings.warn(
                    "Cumulative density calculations experimental for FFT"
                    + " method. Use piecewise method instead.",
                    RuntimeWarning, stacklevel=3,
                )
                _alpha, _beta = pair
                _x = data_subset[:, (0,)]

                # need enough points to "cover" _x for interpolation
                if fft_grid_spacing is None and fft_n_points_two_power is None:
                    raise ValueError(
                        "One of fft_grid_spacing or fft_n_points_two_power "
                        + "needs to be set."
                    )
                max_abs_x = np.max(np.abs(_x))
                h = (
                    2 ** (3 - fft_n_points_two_power) * max_abs_x
                    if fft_grid_spacing is None
                    else fft_grid_spacing
                )
                q = (
                    np.ceil(np.log(2 * max_abs_x / h) / np.log(2)) + 2
                    if fft_n_points_two_power is None
                    else int(fft_n_points_two_power)
                )

                density_x, density = pdf_from_cf_with_fft(
                    lambda t: _cf(t, _alpha, _beta),
                    h=h,
                    q=q,
                    level=fft_interpolation_level,
                )
                f = interpolate.InterpolatedUnivariateSpline(
                    density_x, np.real(density), k=fft_interpolation_degree
                )
                data_out[data_mask] = np.array(
                    [f.integral(self.a, float(x_1.squeeze())) for x_1 in _x]
                ).reshape(data_out[data_mask].shape)

        return data_out.T[0]

    def _fitstart(self, data):
        if self._parameterization() == "S0":
            _fitstart = _fitstart_S0
        elif self._parameterization() == "S1":
            _fitstart = _fitstart_S1
        return _fitstart(data)

    def _stats(self, alpha, beta):
        mu = 0 if alpha > 1 else np.nan
        mu2 = 2 if alpha == 2 else np.inf
        g1 = 0.0 if alpha == 2.0 else np.nan
        g2 = 0.0 if alpha == 2.0 else np.nan
        return mu, mu2, g1, g2


# cotes numbers - see sequence from http://oeis.org/A100642
Cotes_table = np.array(
    [[], [1]] + [v[2] for v in _builtincoeffs.values()], dtype=object
)
Cotes = np.array(
    [
        np.pad(r, (0, len(Cotes_table) - 1 - len(r)), mode='constant')
        for r in Cotes_table
    ]
)


def pdf_from_cf_with_fft(cf, h=0.01, q=9, level=3):
    """Calculates pdf from characteristic function.

    Uses fast Fourier transform with Newton-Cotes integration following [WZ].
    Defaults to using Simpson's method (3-point Newton-Cotes integration).

    Parameters
    ----------
    cf : callable
        Single argument function from float -> complex expressing a
        characteristic function for some distribution.
    h : Optional[float]
        Step size for Newton-Cotes integration. Default: 0.01
    q : Optional[int]
        Use 2**q steps when performing Newton-Cotes integration.
        The infinite integral in the inverse Fourier transform will then
        be restricted to the interval [-2**q * h / 2, 2**q * h / 2]. Setting
        the number of steps equal to a power of 2 allows the fft to be
        calculated in O(n*log(n)) time rather than O(n**2).
        Default: 9
    level : Optional[int]
        Calculate integral using n-point Newton-Cotes integration for
        n = level. The 3-point Newton-Cotes formula corresponds to Simpson's
        rule. Default: 3

    Returns
    -------
    x_l : ndarray
        Array of points x at which pdf is estimated. 2**q equally spaced
        points from -pi/h up to but not including pi/h.
    density : ndarray
        Estimated values of pdf corresponding to cf at points in x_l.

    References
    ----------
    .. [WZ] Wang, Li and Zhang, Ji-Hong, 2008. Simpson's rule based FFT method
        to compute densities of stable distribution.
    """
    n = level
    N = 2**q
    steps = np.arange(0, N)
    L = N * h / 2
    x_l = np.pi * (steps - N / 2) / L
    if level > 1:
        indices = np.arange(n).reshape(n, 1)
        s1 = np.sum(
            (-1) ** steps * Cotes[n, indices] * np.fft.fft(
                (-1)**steps * cf(-L + h * steps + h * indices / (n - 1))
            ) * np.exp(
                1j * np.pi * indices / (n - 1)
                - 2 * 1j * np.pi * indices * steps /
                (N * (n - 1))
            ),
            axis=0
        )
    else:
        s1 = (-1) ** steps * Cotes[n, 0] * np.fft.fft(
            (-1) ** steps * cf(-L + h * steps)
        )
    density = h * s1 / (2 * np.pi * np.sum(Cotes[n]))
    return (x_l, density)


levy_stable = levy_stable_gen(name="levy_stable")