EVOLUTION-MANAGER

Edit File: _continuous_distns.py

# -*- coding: utf-8 -*-
#
# Author:  Travis Oliphant  2002-2011 with contributions from
#          SciPy Developers 2004-2011
#
from __future__ import division, print_function, absolute_import

import warnings
import functools

import numpy as np

from scipy._lib.doccer import (extend_notes_in_docstring,
                               replace_notes_in_docstring)
from scipy import optimize
from scipy import integrate
from scipy import interpolate
import scipy.special as sc
import scipy.special._ufuncs as scu
from scipy._lib._numpy_compat import broadcast_to
from scipy._lib._util import _lazyselect, _lazywhere
from . import _stats
from ._tukeylambda_stats import (tukeylambda_variance as _tlvar,
                                 tukeylambda_kurtosis as _tlkurt)
from ._distn_infrastructure import (get_distribution_names, _kurtosis,
                                    _ncx2_cdf, _ncx2_log_pdf, _ncx2_pdf,
                                    rv_continuous, _skew, valarray,
                                    _get_fixed_fit_value, _check_shape)
from ._constants import _XMIN, _EULER, _ZETA3, _XMAX, _LOGXMAX

# In numpy 1.12 and above, np.power refuses to raise integers to negative
# powers, and `np.float_power` is a new replacement.
try:
    float_power = np.float_power
except AttributeError:
    float_power = np.power

def _remove_optimizer_parameters(kwds):
    """
    Remove the optimizer-related keyword arguments 'loc', 'scale' and
    'optimizer' from `kwds`.  Then check that `kwds` is empty, and
    raise `TypeError("Unknown arguments: %s." % kwds)` if it is not.

This function is used in the fit method of distributions that override
    the default method and do not use the default optimization code.

`kwds` is modified in-place.
    """
    kwds.pop('loc', None)
    kwds.pop('scale', None)
    kwds.pop('optimizer', None)
    if kwds:
        raise TypeError("Unknown arguments: %s." % kwds)

## Kolmogorov-Smirnov one-sided and two-sided test statistics
class ksone_gen(rv_continuous):
    r"""General Kolmogorov-Smirnov one-sided test.

This is the distribution of the one-sided Kolmogorov-Smirnov (KS)
    statistics :math:`D_n^+` and :math:`D_n^-`
    for a finite sample size ``n`` (the shape parameter).

%(before_notes)s

Notes
    -----
    :math:`D_n^+` and :math:`D_n^-` are given by

.. math::

D_n^+ &= \text{sup}_x (F_n(x) - F(x)),\\
        D_n^- &= \text{sup}_x (F(x) - F_n(x)),\\

where :math:`F` is a CDF and :math:`F_n` is an empirical CDF. `ksone`
    describes the distribution under the null hypothesis of the KS test
    that the empirical CDF corresponds to :math:`n` i.i.d. random variates
    with CDF :math:`F`.

%(after_notes)s

See Also
    --------
    kstwobign, kstest

References
    ----------
    .. [1] Birnbaum, Z. W. and Tingey, F.H. "One-sided confidence contours
       for probability distribution functions", The Annals of Mathematical
       Statistics, 22(4), pp 592-596 (1951).

%(example)s

"""
    def _pdf(self, x, n):
        return -scu._smirnovp(n, x)

def _cdf(self, x, n):
        return scu._smirnovc(n, x)

def _sf(self, x, n):
        return sc.smirnov(n, x)

def _ppf(self, q, n):
        return scu._smirnovci(n, q)

def _isf(self, q, n):
        return sc.smirnovi(n, q)

ksone = ksone_gen(a=0.0, b=1.0, name='ksone')

class kstwobign_gen(rv_continuous):
    r"""Kolmogorov-Smirnov two-sided test for large N.

This is the asymptotic distribution of the two-sided Kolmogorov-Smirnov
    statistic :math:`\sqrt{n} D_n` that measures the maximum absolute
    distance of the theoretical CDF from the empirical CDF (see `kstest`).

%(before_notes)s

Notes
    -----
    :math:`\sqrt{n} D_n` is given by

.. math::

D_n = \text{sup}_x |F_n(x) - F(x)|

where :math:`F` is a CDF and :math:`F_n` is an empirical CDF. `kstwobign`
    describes the asymptotic distribution (i.e. the limit of
    :math:`\sqrt{n} D_n`) under the null hypothesis of the KS test that the
    empirical CDF corresponds to i.i.d. random variates with CDF :math:`F`.

%(after_notes)s

See Also
    --------
    levy, levy_l

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:

.. math::

\Phi = \begin{cases}
                \tan \left({\frac {\pi \alpha }{2}}\right)&\alpha \neq 1\\
                -{\frac {2}{\pi }}\log |t|&\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.

For evaluation of pdf we use either Zolotarev :math:`S_0` parameterization with integration,
    direct integration of standard parameterization of characteristic function or FFT of
    characteristic function. If set to other than None and if number of points is greater than
    ``levy_stable.pdf_fft_min_points_threshold`` (defaults to None) we use FFT otherwise we use one
    of the other methods.

The default method is 'best' which uses Zolotarev's method if alpha = 1 and integration of
    characteristic function otherwise. The default method can be changed by setting
    ``levy_stable.pdf_default_method`` to either 'zolotarev', 'quadrature' or 'best'.

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 a value that covers the
    input range * 4). Setting ``pdf_fft_n_points_two_power`` to 16 should be sufficiently accurate
    in most cases at the expense of CPU time.

For evaluation of cdf we use Zolatarev :math:`S_0` parameterization with integration or integral of
    the pdf FFT interpolated spline. The settings affecting FFT calculation are the same as
    for pdf calculation. Setting the threshold to ``None`` (default) will disable FFT. 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; so it's best that ``pdf_fft_min_points_threshold`` is left unset.

.. warning::

For pdf calculations implementation of Zolatarev is unstable for values where alpha = 1 and
        beta != 0. In this case the quadrature method is recommended. FFT calculation is also
        considered experimental.

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

%(after_notes)s

References
    ----------
    .. [MC] McCulloch, J., 1986. Simple consistent estimators of stable distribution parameters.
       Communications in Statistics - Simulation and Computation 15, 11091136.
    .. [MS] Mittnik, S.T. Rachev, T. Doganoglu, D. Chenyao, 1999. Maximum likelihood estimation
       of stable Paretian models, Mathematical and Computer Modelling, Volume 29, Issue 10,
       1999, Pages 275-293.
    .. [BS] Borak, S., Hardle, W., Rafal, W. 2005. Stable distributions, Economic Risk.

%(example)s

"""

def _rvs(self, alpha, beta):

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

sz = self._size
        alpha = broadcast_to(alpha, sz)
        beta = broadcast_to(beta, sz)
        TH = uniform.rvs(loc=-np.pi/2.0, scale=np.pi, size=sz,
                         random_state=self._random_state)
        W = expon.rvs(size=sz, random_state=self._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 _argcheck(self, alpha, beta):
        return (alpha > 0) & (alpha <= 2) & (beta <= 1) & (beta >= -1)

@staticmethod
    def _cf(t, alpha, beta):
        Phi = lambda alpha, t: np.tan(np.pi*alpha/2) if alpha != 1 else -2.0*np.log(np.abs(t))/np.pi
        return np.exp(-(np.abs(t)**alpha)*(1-1j*beta*np.sign(t)*Phi(alpha, t)))

@staticmethod
    def _pdf_from_cf_with_fft(cf, h=0.01, q=9):
        """Calculates pdf from cf using fft. Using region around 0 with N=2**q points
        separated by distance h. As suggested by [MS].
        """
        N = 2**q
        n = np.arange(1,N+1)
        density = ((-1)**(n-1-N/2))*np.fft.fft(((-1)**(n-1))*cf(2*np.pi*(n-1-N/2)/h/N))/h/N
        x = (n-1-N/2)*h
        return (x, density)

@staticmethod
    def _pdf_single_value_best(x, alpha, beta):
        if alpha != 1. or (alpha == 1. and beta == 0.):
            return levy_stable_gen._pdf_single_value_zolotarev(x, alpha, beta)
        else:
            return levy_stable_gen._pdf_single_value_cf_integrate(x, alpha, beta)

@staticmethod
    def _pdf_single_value_cf_integrate(x, alpha, beta):
        cf = lambda t: levy_stable_gen._cf(t, alpha, beta)
        return integrate.quad(lambda t: np.real(np.exp(-1j*t*x)*cf(t)), -np.inf, np.inf, limit=1000)[0]/np.pi/2

@staticmethod
    def _pdf_single_value_zolotarev(x, alpha, beta):
        """Calculate pdf using Zolotarev's methods as detailed in [BS].
        """
        zeta = -beta*np.tan(np.pi*alpha/2.)
        if alpha != 1:
            x0 = x + zeta  # convert to S_0 parameterization
            xi = np.arctan(-zeta)/alpha

def V(theta):
                return np.cos(alpha*xi)**(1/(alpha-1)) * \
                                (np.cos(theta)/np.sin(alpha*(xi+theta)))**(alpha/(alpha-1)) * \
                                (np.cos(alpha*xi+(alpha-1)*theta)/np.cos(theta))
            if x0 > zeta:
                def g(theta):
                    return V(theta)*np.real(np.complex(x0-zeta)**(alpha/(alpha-1)))

def f(theta):
                    return g(theta) * np.exp(-g(theta))

# 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.

with np.errstate(all="ignore"):
                    intg_max = optimize.minimize_scalar(lambda theta: -f(theta), bounds=[-xi, np.pi/2])
                    intg_kwargs = {}
                    # windows quadpack less forgiving with points out of bounds
                    if intg_max.success and not np.isnan(intg_max.fun)\
                            and intg_max.x > -xi and intg_max.x < np.pi/2:
                        intg_kwargs["points"] = [intg_max.x]
                    intg = integrate.quad(f, -xi, np.pi/2, **intg_kwargs)[0]
                    return alpha * intg / np.pi / np.abs(alpha-1) / (x0-zeta)
            elif x0 == zeta:
                return sc.gamma(1+1/alpha)*np.cos(xi)/np.pi/((1+zeta**2)**(1/alpha/2))
            else:
                return levy_stable_gen._pdf_single_value_zolotarev(-x, alpha, -beta)
        else:
            # since location zero, no need to reposition x for S_0 parameterization
            xi = np.pi/2
            if beta != 0:
                warnings.warn('Density calculation unstable for alpha=1 and beta!=0.' +
                              ' Use quadrature method instead.', RuntimeWarning)

def V(theta):
                    expr_1 = np.pi/2+beta*theta
                    return 2. * expr_1 * np.exp(expr_1*np.tan(theta)/beta) / np.cos(theta) / np.pi

def g(theta):
                    return np.exp(-np.pi * x / 2. / beta) * V(theta)

def f(theta):
                    return g(theta) * np.exp(-g(theta))

with np.errstate(all="ignore"):
                    intg_max = optimize.minimize_scalar(lambda theta: -f(theta), bounds=[-np.pi/2, np.pi/2])
                    intg = integrate.fixed_quad(f, -np.pi/2, intg_max.x)[0] + integrate.fixed_quad(f, intg_max.x, np.pi/2)[0]
                    return intg / np.abs(beta) / 2.
            else:
                return 1/(1+x**2)/np.pi

@staticmethod
    def _cdf_single_value_zolotarev(x, alpha, beta):
        """Calculate cdf using Zolotarev's methods as detailed in [BS].
        """
        zeta = -beta*np.tan(np.pi*alpha/2.)
        if alpha != 1:
            x0 = x + zeta  # convert to S_0 parameterization
            xi = np.arctan(-zeta)/alpha

def f(theta):
                    return np.exp(-V(theta)*np.real(np.complex(x0-zeta)**(alpha/(alpha-1))))

with np.errstate(all="ignore"):
                    # 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):
                        intg = 0
                    else:
                        intg = integrate.quad(f, -xi, np.pi/2)[0]
                    return c_1 + np.sign(1-alpha) * intg / np.pi
            elif x0 == zeta:
                return .5 - xi/np.pi
            else:
                return 1 - levy_stable_gen._cdf_single_value_zolotarev(-x, alpha, -beta)

else:
            # since location zero, no need to reposition x for S_0 parameterization
            xi = np.pi/2
            if beta > 0:

def V(theta):
                    expr_1 = np.pi/2+beta*theta
                    return 2. * expr_1 * np.exp(expr_1*np.tan(theta)/beta) / np.cos(theta) / np.pi

with np.errstate(all="ignore"):
                    expr_1 = np.exp(-np.pi*x/beta/2.)
                    int_1 = integrate.quad(lambda theta: np.exp(-expr_1 * V(theta)), -np.pi/2, np.pi/2)[0]
                    return int_1 / np.pi
            elif beta == 0:
                return .5 + np.arctan(x)/np.pi
            else:
                return 1 - levy_stable_gen._cdf_single_value_zolotarev(-x, 1, -beta)

def _pdf(self, x, alpha, beta):

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 = getattr(self, 'pdf_default_method', 'best')
        if pdf_default_method_name == 'best':
            pdf_single_value_method = levy_stable_gen._pdf_single_value_best
        elif pdf_default_method_name == 'zolotarev':
            pdf_single_value_method = levy_stable_gen._pdf_single_value_zolotarev
        else:
            pdf_single_value_method = levy_stable_gen._pdf_single_value_cf_integrate

fft_min_points_threshold = getattr(self, 'pdf_fft_min_points_threshold', None)
        fft_grid_spacing = getattr(self, 'pdf_fft_grid_spacing', 0.001)
        fft_n_points_two_power = getattr(self, 'pdf_fft_n_points_two_power', None)

# group data in unique arrays of alpha, beta pairs
        uniq_param_pairs = np.vstack(list({tuple(row) for row in
                                           data_in[:, 1:]}))
        for pair in uniq_param_pairs:
            data_mask = np.all(data_in[:,1:] == pair, axis=-1)
            data_subset = data_in[data_mask]
            if fft_min_points_threshold is None or len(data_subset) < fft_min_points_threshold:
                data_out[data_mask] = np.array([pdf_single_value_method(_x, _alpha, _beta)
                            for _x, _alpha, _beta in data_subset]).reshape(len(data_subset), 1)
            else:
                warnings.warn('Density calculations experimental for FFT method.' +
                              ' Use combination of zolatarev and quadrature methods instead.', RuntimeWarning)
                _alpha, _beta = pair
                _x = data_subset[:,(0,)]

# need enough points to "cover" _x for interpolation
                h = fft_grid_spacing
                q = np.ceil(np.log(2*np.max(np.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 = levy_stable_gen._pdf_from_cf_with_fft(lambda t: levy_stable_gen._cf(t, _alpha, _beta), h=h, q=q)
                f = interpolate.interp1d(density_x, np.real(density))
                data_out[data_mask] = f(_x)

return data_out.T[0]

def _cdf(self, x, alpha, beta):

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))

# group data in unique arrays of alpha, beta pairs
        uniq_param_pairs = np.vstack(
            list({tuple(row) for row in data_in[:,1:]}))
        for pair in uniq_param_pairs:
            data_mask = np.all(data_in[:,1:] == pair, axis=-1)
            data_subset = data_in[data_mask]
            if fft_min_points_threshold is None or len(data_subset) < fft_min_points_threshold:
                data_out[data_mask] = np.array([levy_stable._cdf_single_value_zolotarev(_x, _alpha, _beta)
                            for _x, _alpha, _beta in data_subset]).reshape(len(data_subset), 1)
            else:
                warnings.warn(u'FFT method is considered experimental for ' +
                              u'cumulative distribution function ' +
                              u'evaluations. Use Zolotarev’s method instead).',
                              RuntimeWarning)
                _alpha, _beta = pair
                _x = data_subset[:,(0,)]

# need enough points to "cover" _x for interpolation
                h = fft_grid_spacing
                q = 16 if fft_n_points_two_power is None else int(fft_n_points_two_power)

density_x, density = levy_stable_gen._pdf_from_cf_with_fft(lambda t: levy_stable_gen._cf(t, _alpha, _beta), h=h, q=q)
                f = interpolate.InterpolatedUnivariateSpline(density_x, np.real(density))
                data_out[data_mask] = np.array([f.integral(self.a, x_1) for x_1 in _x]).reshape(data_out[data_mask].shape)

return data_out.T[0]

def _fitstart(self, data):
        # We follow McCullock 1986 method - Simple Consistent Estimators
        # of Stable Distribution Parameters

# 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]]

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)

# 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.)
            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 = 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)

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

levy_stable = levy_stable_gen(name='levy_stable')

class logistic_gen(rv_continuous):
    r"""A logistic (or Sech-squared) continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `logistic` is:

.. math::

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

`logistic` is a special case of `genlogistic` with ``c=1``.

%(after_notes)s

%(example)s

"""
    def _rvs(self):
        return self._random_state.logistic(size=self._size)

def _pdf(self, x):
        # logistic.pdf(x) = exp(-x) / (1+exp(-x))**2
        return np.exp(self._logpdf(x))

def _logpdf(self, x):
        return -x - 2. * sc.log1p(np.exp(-x))

def _cdf(self, x):
        return sc.expit(x)

def _ppf(self, q):
        return sc.logit(q)

def _sf(self, x):
        return sc.expit(-x)

def _isf(self, q):
        return -sc.logit(q)

def _stats(self):
        return 0, np.pi*np.pi/3.0, 0, 6.0/5.0

def _entropy(self):
        # https://en.wikipedia.org/wiki/Logistic_distribution
        return 2.0

logistic = logistic_gen(name='logistic')

class loggamma_gen(rv_continuous):
    r"""A log gamma continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `loggamma` is:

.. math::

f(x, c) = \frac{\exp(c x - \exp(x))}
                       {\Gamma(c)}

for all :math:`x, c > 0`. Here, :math:`\Gamma` is the
    gamma function (`scipy.special.gamma`).

`loggamma` takes ``c`` as a shape parameter for :math:`c`.

%(after_notes)s

%(example)s

"""
    def _rvs(self, c):
        return np.log(self._random_state.gamma(c, size=self._size))

def _pdf(self, x, c):
        # loggamma.pdf(x, c) = exp(c*x-exp(x)) / gamma(c)
        return np.exp(c*x-np.exp(x)-sc.gammaln(c))

def _cdf(self, x, c):
        return sc.gammainc(c, np.exp(x))

def _ppf(self, q, c):
        return np.log(sc.gammaincinv(c, q))

def _stats(self, c):
        # See, for example, "A Statistical Study of Log-Gamma Distribution", by
        # Ping Shing Chan (thesis, McMaster University, 1993).
        mean = sc.digamma(c)
        var = sc.polygamma(1, c)
        skewness = sc.polygamma(2, c) / np.power(var, 1.5)
        excess_kurtosis = sc.polygamma(3, c) / (var*var)
        return mean, var, skewness, excess_kurtosis

loggamma = loggamma_gen(name='loggamma')

class loglaplace_gen(rv_continuous):
    r"""A log-Laplace continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `loglaplace` is:

.. math::

f(x, c) = \begin{cases}\frac{c}{2} x^{ c-1}  &\text{for } 0 < x < 1\\
                               \frac{c}{2} x^{-c-1}  &\text{for } x \ge 1
                  \end{cases}

for :math:`c > 0`.

`loglaplace` takes ``c`` as a shape parameter for :math:`c`.

%(after_notes)s

References
    ----------
    T.J. Kozubowski and K. Podgorski, "A log-Laplace growth rate model",
    The Mathematical Scientist, vol. 28, pp. 49-60, 2003.

%(example)s

"""
    def _pdf(self, x, c):
        # loglaplace.pdf(x, c) = c / 2 * x**(c-1),   for 0 < x < 1
        #                      = c / 2 * x**(-c-1),  for x >= 1
        cd2 = c/2.0
        c = np.where(x < 1, c, -c)
        return cd2*x**(c-1)

def _cdf(self, x, c):
        return np.where(x < 1, 0.5*x**c, 1-0.5*x**(-c))

def _ppf(self, q, c):
        return np.where(q < 0.5, (2.0*q)**(1.0/c), (2*(1.0-q))**(-1.0/c))

def _munp(self, n, c):
        return c**2 / (c**2 - n**2)

def _entropy(self, c):
        return np.log(2.0/c) + 1.0

loglaplace = loglaplace_gen(a=0.0, name='loglaplace')

def _lognorm_logpdf(x, s):
    return _lazywhere(x != 0, (x, s),
                      lambda x, s: -np.log(x)**2 / (2*s**2) - np.log(s*x*np.sqrt(2*np.pi)),
                      -np.inf)

class lognorm_gen(rv_continuous):
    r"""A lognormal continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `lognorm` is:

.. math::

f(x, s) = \frac{1}{s x \sqrt{2\pi}}
                  \exp\left(-\frac{\log^2(x)}{2s^2}\right)

for :math:`x > 0`, :math:`s > 0`.

`lognorm` takes ``s`` as a shape parameter for :math:`s`.

%(after_notes)s

A common parametrization for a lognormal random variable ``Y`` is in
    terms of the mean, ``mu``, and standard deviation, ``sigma``, of the
    unique normally distributed random variable ``X`` such that exp(X) = Y.
    This parametrization corresponds to setting ``s = sigma`` and ``scale =
    exp(mu)``.

%(example)s

"""
    _support_mask = rv_continuous._open_support_mask

def _rvs(self, s):
        return np.exp(s * self._random_state.standard_normal(self._size))

def _pdf(self, x, s):
        # lognorm.pdf(x, s) = 1 / (s*x*sqrt(2*pi)) * exp(-1/2*(log(x)/s)**2)
        return np.exp(self._logpdf(x, s))

def _logpdf(self, x, s):
        return _lognorm_logpdf(x, s)

def _cdf(self, x, s):
        return _norm_cdf(np.log(x) / s)

def _logcdf(self, x, s):
        return _norm_logcdf(np.log(x) / s)

def _ppf(self, q, s):
        return np.exp(s * _norm_ppf(q))

def _sf(self, x, s):
        return _norm_sf(np.log(x) / s)

def _logsf(self, x, s):
        return _norm_logsf(np.log(x) / s)

def _stats(self, s):
        p = np.exp(s*s)
        mu = np.sqrt(p)
        mu2 = p*(p-1)
        g1 = np.sqrt((p-1))*(2+p)
        g2 = np.polyval([1, 2, 3, 0, -6.0], p)
        return mu, mu2, g1, g2

def _entropy(self, s):
        return 0.5 * (1 + np.log(2*np.pi) + 2 * np.log(s))

@extend_notes_in_docstring(rv_continuous, notes="""\
        When the location parameter is fixed by using the `floc` argument,
        this function uses explicit formulas for the maximum likelihood
        estimation of the log-normal shape and scale parameters, so the
        `optimizer`, `loc` and `scale` keyword arguments are ignored.\n\n""")
    def fit(self, data, *args, **kwds):
        floc = kwds.get('floc', None)
        if floc is None:
            # loc is not fixed.  Use the default fit method.
            return super(lognorm_gen, self).fit(data, *args, **kwds)

f0 = (kwds.get('f0', None) or kwds.get('fs', None) or
              kwds.get('fix_s', None))
        fscale = kwds.get('fscale', None)

if len(args) > 1:
            raise TypeError("Too many input arguments.")
        for name in ['f0', 'fs', 'fix_s', 'floc', 'fscale', 'loc', 'scale',
                     'optimizer']:
            kwds.pop(name, None)
        if kwds:
            raise TypeError("Unknown arguments: %s." % kwds)

# Special case: loc is fixed.  Use the maximum likelihood formulas
        # instead of the numerical solver.

if f0 is not None and fscale is not None:
            # This check is for consistency with `rv_continuous.fit`.
            raise ValueError("All parameters fixed. There is nothing to "
                             "optimize.")

data = np.asarray(data)

if not np.isfinite(data).all():
            raise RuntimeError("The data contains non-finite values.")

floc = float(floc)
        if floc != 0:
            # Shifting the data by floc. Don't do the subtraction in-place,
            # because `data` might be a view of the input array.
            data = data - floc
        if np.any(data <= 0):
            raise FitDataError("lognorm", lower=floc, upper=np.inf)
        lndata = np.log(data)

# Three cases to handle:
        # * shape and scale both free
        # * shape fixed, scale free
        # * shape free, scale fixed

if fscale is None:
            # scale is free.
            scale = np.exp(lndata.mean())
            if f0 is None:
                # shape is free.
                shape = lndata.std()
            else:
                # shape is fixed.
                shape = float(f0)
        else:
            # scale is fixed, shape is free
            scale = float(fscale)
            shape = np.sqrt(((lndata - np.log(scale))**2).mean())

return shape, floc, scale

lognorm = lognorm_gen(a=0.0, name='lognorm')

class gilbrat_gen(rv_continuous):
    r"""A Gilbrat continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `gilbrat` is:

.. math::

f(x) = \frac{1}{x \sqrt{2\pi}} \exp(-\frac{1}{2} (\log(x))^2)

`gilbrat` is a special case of `lognorm` with ``s=1``.

%(after_notes)s

%(example)s

"""
    _support_mask = rv_continuous._open_support_mask

def _rvs(self):
        return np.exp(self._random_state.standard_normal(self._size))

def _pdf(self, x):
        # gilbrat.pdf(x) = 1/(x*sqrt(2*pi)) * exp(-1/2*(log(x))**2)
        return np.exp(self._logpdf(x))

def _logpdf(self, x):
        return _lognorm_logpdf(x, 1.0)

def _cdf(self, x):
        return _norm_cdf(np.log(x))

def _ppf(self, q):
        return np.exp(_norm_ppf(q))

def _stats(self):
        p = np.e
        mu = np.sqrt(p)
        mu2 = p * (p - 1)
        g1 = np.sqrt((p - 1)) * (2 + p)
        g2 = np.polyval([1, 2, 3, 0, -6.0], p)
        return mu, mu2, g1, g2

def _entropy(self):
        return 0.5 * np.log(2 * np.pi) + 0.5

gilbrat = gilbrat_gen(a=0.0, name='gilbrat')

class maxwell_gen(rv_continuous):
    r"""A Maxwell continuous random variable.

%(before_notes)s

Notes
    -----
    A special case of a `chi` distribution,  with ``df=3``, ``loc=0.0``,
    and given ``scale = a``, where ``a`` is the parameter used in the
    Mathworld description [1]_.

The probability density function for `maxwell` is:

.. math::

f(x) = \sqrt{2/\pi}x^2 \exp(-x^2/2)

for :math:`x >= 0`.

%(after_notes)s

References
    ----------
    .. [1] http://mathworld.wolfram.com/MaxwellDistribution.html

%(example)s
    """
    def _rvs(self):
        return chi.rvs(3.0, size=self._size, random_state=self._random_state)

def _pdf(self, x):
        # maxwell.pdf(x) = sqrt(2/pi)x**2 * exp(-x**2/2)
        return np.sqrt(2.0/np.pi)*x*x*np.exp(-x*x/2.0)

def _cdf(self, x):
        return sc.gammainc(1.5, x*x/2.0)

def _ppf(self, q):
        return np.sqrt(2*sc.gammaincinv(1.5, q))

def _stats(self):
        val = 3*np.pi-8
        return (2*np.sqrt(2.0/np.pi),
                3-8/np.pi,
                np.sqrt(2)*(32-10*np.pi)/val**1.5,
                (-12*np.pi*np.pi + 160*np.pi - 384) / val**2.0)

def _entropy(self):
        return _EULER + 0.5*np.log(2*np.pi)-0.5

maxwell = maxwell_gen(a=0.0, name='maxwell')

class mielke_gen(rv_continuous):
    r"""A Mielke Beta-Kappa / Dagum continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `mielke` is:

.. math::

f(x, k, s) = \frac{k x^{k-1}}{(1+x^s)^{1+k/s}}

for :math:`x > 0` and :math:`k, s > 0`. The distribution is sometimes
    called Dagum distribution ([2]_). It was already defined in [3]_, called
    a Burr Type III distribution (`burr` with parameters ``c=s`` and
    ``d=k/s``).

`mielke` takes ``k`` and ``s`` as shape parameters.

%(after_notes)s

References
    ----------
    .. [1] Mielke, P.W., 1973 "Another Family of Distributions for Describing
           and Analyzing Precipitation Data." J. Appl. Meteor., 12, 275-280
    .. [2] Dagum, C., 1977 "A new model for personal income distribution."
           Economie Appliquee, 33, 327-367.
    .. [3] Burr, I. W. "Cumulative frequency functions", Annals of
           Mathematical Statistics, 13(2), pp 215-232 (1942).

%(example)s

"""
    def _argcheck(self, k, s):
        return (k > 0) & (s > 0)

def _pdf(self, x, k, s):
        return k*x**(k-1.0) / (1.0+x**s)**(1.0+k*1.0/s)

def _logpdf(self, x, k, s):
        return np.log(k) + np.log(x)*(k-1.0) - np.log1p(x**s)*(1.0+k*1.0/s)

def _cdf(self, x, k, s):
        return x**k / (1.0+x**s)**(k*1.0/s)

def _ppf(self, q, k, s):
        qsk = pow(q, s*1.0/k)
        return pow(qsk/(1.0-qsk), 1.0/s)

def _munp(self, n, k, s):
        def nth_moment(n, k, s):
            # n-th moment is defined for -k < n < s
            return sc.gamma((k+n)/s)*sc.gamma(1-n/s)/sc.gamma(k/s)

return _lazywhere(n < s, (n, k, s), nth_moment, np.inf)

mielke = mielke_gen(a=0.0, name='mielke')

class kappa4_gen(rv_continuous):
    r"""Kappa 4 parameter distribution.

%(before_notes)s

Notes
    -----
    The probability density function for kappa4 is:

.. math::

f(x, h, k) = (1 - k x)^{1/k - 1} (1 - h (1 - k x)^{1/k})^{1/h-1}

if :math:`h` and :math:`k` are not equal to 0.

If :math:`h` or :math:`k` are zero then the pdf can be simplified:

h = 0 and k != 0::

kappa4.pdf(x, h, k) = (1.0 - k*x)**(1.0/k - 1.0)*
                              exp(-(1.0 - k*x)**(1.0/k))

h != 0 and k = 0::

kappa4.pdf(x, h, k) = exp(-x)*(1.0 - h*exp(-x))**(1.0/h - 1.0)

h = 0 and k = 0::

kappa4.pdf(x, h, k) = exp(-x)*exp(-exp(-x))

kappa4 takes :math:`h` and :math:`k` as shape parameters.

The kappa4 distribution returns other distributions when certain
    :math:`h` and :math:`k` values are used.

(1) There are at least five generalized logistic distributions.
        Four are described here:
        https://en.wikipedia.org/wiki/Generalized_logistic_distribution
        The "fifth" one is the one kappa4 should match which currently
        isn't implemented in scipy:
        https://en.wikipedia.org/wiki/Talk:Generalized_logistic_distribution
        https://www.mathwave.com/help/easyfit/html/analyses/distributions/gen_logistic.html
    (2) This distribution is currently not in scipy.

References
    ----------
    J.C. Finney, "Optimization of a Skewed Logistic Distribution With Respect
    to the Kolmogorov-Smirnov Test", A Dissertation Submitted to the Graduate
    Faculty of the Louisiana State University and Agricultural and Mechanical
    College, (August, 2004),
    https://digitalcommons.lsu.edu/gradschool_dissertations/3672

J.R.M. Hosking, "The four-parameter kappa distribution". IBM J. Res.
    Develop. 38 (3), 25 1-258 (1994).

B. Kumphon, A. Kaew-Man, P. Seenoi, "A Rainfall Distribution for the Lampao
    Site in the Chi River Basin, Thailand", Journal of Water Resource and
    Protection, vol. 4, 866-869, (2012).
    https://doi.org/10.4236/jwarp.2012.410101

C. Winchester, "On Estimation of the Four-Parameter Kappa Distribution", A
    Thesis Submitted to Dalhousie University, Halifax, Nova Scotia, (March
    2000).
    http://www.nlc-bnc.ca/obj/s4/f2/dsk2/ftp01/MQ57336.pdf

%(after_notes)s

%(example)s

"""
    def _argcheck(self, h, k):
        return h == h

def _get_support(self, h, k):
        condlist = [np.logical_and(h > 0, k > 0),
                    np.logical_and(h > 0, k == 0),
                    np.logical_and(h > 0, k < 0),
                    np.logical_and(h <= 0, k > 0),
                    np.logical_and(h <= 0, k == 0),
                    np.logical_and(h <= 0, k < 0)]

def f0(h, k):
            return (1.0 - float_power(h, -k))/k

def f1(h, k):
            return np.log(h)

def f3(h, k):
            a = np.empty(np.shape(h))
            a[:] = -np.inf
            return a

def f5(h, k):
            return 1.0/k

_a = _lazyselect(condlist,
                             [f0, f1, f0, f3, f3, f5],
                             [h, k],
                             default=np.nan)

def f0(h, k):
            return 1.0/k

def f1(h, k):
            a = np.empty(np.shape(h))
            a[:] = np.inf
            return a

_b = _lazyselect(condlist,
                             [f0, f1, f1, f0, f1, f1],
                             [h, k],
                             default=np.nan)
        return _a, _b

def _pdf(self, x, h, k):
        # kappa4.pdf(x, h, k) = (1.0 - k*x)**(1.0/k - 1.0)*
        #                       (1.0 - h*(1.0 - k*x)**(1.0/k))**(1.0/h-1)
        return np.exp(self._logpdf(x, h, k))

def _logpdf(self, x, h, k):
        condlist = [np.logical_and(h != 0, k != 0),
                    np.logical_and(h == 0, k != 0),
                    np.logical_and(h != 0, k == 0),
                    np.logical_and(h == 0, k == 0)]

def f0(x, h, k):
            '''pdf = (1.0 - k*x)**(1.0/k - 1.0)*(
                      1.0 - h*(1.0 - k*x)**(1.0/k))**(1.0/h-1.0)
               logpdf = ...
            '''
            return (sc.xlog1py(1.0/k - 1.0, -k*x) +
                    sc.xlog1py(1.0/h - 1.0, -h*(1.0 - k*x)**(1.0/k)))

def f1(x, h, k):
            '''pdf = (1.0 - k*x)**(1.0/k - 1.0)*np.exp(-(
                      1.0 - k*x)**(1.0/k))
               logpdf = ...
            '''
            return sc.xlog1py(1.0/k - 1.0, -k*x) - (1.0 - k*x)**(1.0/k)

def f2(x, h, k):
            '''pdf = np.exp(-x)*(1.0 - h*np.exp(-x))**(1.0/h - 1.0)
               logpdf = ...
            '''
            return -x + sc.xlog1py(1.0/h - 1.0, -h*np.exp(-x))

def f3(x, h, k):
            '''pdf = np.exp(-x-np.exp(-x))
               logpdf = ...
            '''
            return -x - np.exp(-x)

return _lazyselect(condlist,
                           [f0, f1, f2, f3],
                           [x, h, k],
                           default=np.nan)

def _cdf(self, x, h, k):
        return np.exp(self._logcdf(x, h, k))

def _logcdf(self, x, h, k):
        condlist = [np.logical_and(h != 0, k != 0),
                    np.logical_and(h == 0, k != 0),
                    np.logical_and(h != 0, k == 0),
                    np.logical_and(h == 0, k == 0)]

def f0(x, h, k):
            '''cdf = (1.0 - h*(1.0 - k*x)**(1.0/k))**(1.0/h)
               logcdf = ...
            '''
            return (1.0/h)*sc.log1p(-h*(1.0 - k*x)**(1.0/k))

def f1(x, h, k):
            '''cdf = np.exp(-(1.0 - k*x)**(1.0/k))
               logcdf = ...
            '''
            return -(1.0 - k*x)**(1.0/k)

def f2(x, h, k):
            '''cdf = (1.0 - h*np.exp(-x))**(1.0/h)
               logcdf = ...
            '''
            return (1.0/h)*sc.log1p(-h*np.exp(-x))

def f3(x, h, k):
            '''cdf = np.exp(-np.exp(-x))
               logcdf = ...
            '''
            return -np.exp(-x)

return _lazyselect(condlist,
                           [f0, f1, f2, f3],
                           [x, h, k],
                           default=np.nan)

def _ppf(self, q, h, k):
        condlist = [np.logical_and(h != 0, k != 0),
                    np.logical_and(h == 0, k != 0),
                    np.logical_and(h != 0, k == 0),
                    np.logical_and(h == 0, k == 0)]

def f0(q, h, k):
            return 1.0/k*(1.0 - ((1.0 - (q**h))/h)**k)

def f1(q, h, k):
            return 1.0/k*(1.0 - (-np.log(q))**k)

def f2(q, h, k):
            '''ppf = -np.log((1.0 - (q**h))/h)
            '''
            return -sc.log1p(-(q**h)) + np.log(h)

def f3(q, h, k):
            return -np.log(-np.log(q))

return _lazyselect(condlist,
                           [f0, f1, f2, f3],
                           [q, h, k],
                           default=np.nan)

def _stats(self, h, k):
        if h >= 0 and k >= 0:
            maxr = 5
        elif h < 0 and k >= 0:
            maxr = int(-1.0/h*k)
        elif k < 0:
            maxr = int(-1.0/k)
        else:
            maxr = 5

outputs = [None if r < maxr else np.nan for r in range(1, 5)]
        return outputs[:]

kappa4 = kappa4_gen(name='kappa4')

class kappa3_gen(rv_continuous):
    r"""Kappa 3 parameter distribution.

%(before_notes)s

Notes
    -----
    The probability density function for `kappa3` is:

.. math::

f(x, a) = a (a + x^a)^{-(a + 1)/a}

for :math:`x > 0` and :math:`a > 0`.

`kappa3` takes ``a`` as a shape parameter for :math:`a`.

References
    ----------
    P.W. Mielke and E.S. Johnson, "Three-Parameter Kappa Distribution Maximum
    Likelihood and Likelihood Ratio Tests", Methods in Weather Research,
    701-707, (September, 1973),
    https://doi.org/10.1175/1520-0493(1973)101<0701:TKDMLE>2.3.CO;2

B. Kumphon, "Maximum Entropy and Maximum Likelihood Estimation for the
    Three-Parameter Kappa Distribution", Open Journal of Statistics, vol 2,
    415-419 (2012), https://doi.org/10.4236/ojs.2012.24050

%(after_notes)s

%(example)s

"""
    def _argcheck(self, a):
        return a > 0

def _pdf(self, x, a):
        # kappa3.pdf(x, a) = a*(a + x**a)**(-(a + 1)/a),     for x > 0
        return a*(a + x**a)**(-1.0/a-1)

def _cdf(self, x, a):
        return x*(a + x**a)**(-1.0/a)

def _ppf(self, q, a):
        return (a/(q**-a - 1.0))**(1.0/a)

def _stats(self, a):
        outputs = [None if i < a else np.nan for i in range(1, 5)]
        return outputs[:]

kappa3 = kappa3_gen(a=0.0, name='kappa3')

class moyal_gen(rv_continuous):
    r"""A Moyal continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `moyal` is:

.. math::

f(x) = \exp(-(x + \exp(-x))/2) / \sqrt{2\pi}

for a real number :math:`x`.

%(after_notes)s

This distribution has utility in high-energy physics and radiation
    detection. It describes the energy loss of a charged relativistic
    particle due to ionization of the medium [1]_. It also provides an
    approximation for the Landau distribution. For an in depth description
    see [2]_. For additional description, see [3]_.

References
    ----------
    .. [1] J.E. Moyal, "XXX. Theory of ionization fluctuations",
           The London, Edinburgh, and Dublin Philosophical Magazine
           and Journal of Science, vol 46, 263-280, (1955).
           :doi:`10.1080/14786440308521076` (gated)
    .. [2] G. Cordeiro et al., "The beta Moyal: a useful skew distribution",
           International Journal of Research and Reviews in Applied Sciences,
           vol 10, 171-192, (2012).
           http://www.arpapress.com/Volumes/Vol10Issue2/IJRRAS_10_2_02.pdf
    .. [3] C. Walck, "Handbook on Statistical Distributions for
           Experimentalists; International Report SUF-PFY/96-01", Chapter 26,
           University of Stockholm: Stockholm, Sweden, (2007).
           http://www.stat.rice.edu/~dobelman/textfiles/DistributionsHandbook.pdf

.. versionadded:: 1.1.0

%(example)s

"""
    def _rvs(self):
        sz, rndm = self._size, self._random_state
        u1 = gamma.rvs(a = 0.5, scale = 2, size=sz, random_state=rndm)
        return -np.log(u1)

def _pdf(self, x):
        return np.exp(-0.5 * (x + np.exp(-x))) / np.sqrt(2*np.pi)

def _cdf(self, x):
        return sc.erfc(np.exp(-0.5 * x) / np.sqrt(2))

def _sf(self, x):
        return sc.erf(np.exp(-0.5 * x) / np.sqrt(2))

def _ppf(self, x):
        return -np.log(2 * sc.erfcinv(x)**2)

def _stats(self):
        mu = np.log(2) + np.euler_gamma
        mu2 = np.pi**2 / 2
        g1 = 28 * np.sqrt(2) * sc.zeta(3) / np.pi**3
        g2 = 4.
        return mu, mu2, g1, g2

def _munp(self, n):
        if n == 1.0:
            return np.log(2) + np.euler_gamma
        elif n == 2.0:
            return np.pi**2 / 2 + (np.log(2) + np.euler_gamma)**2
        elif n == 3.0:
            tmp1 = 1.5 * np.pi**2 * (np.log(2)+np.euler_gamma)
            tmp2 = (np.log(2)+np.euler_gamma)**3
            tmp3 = 14 * sc.zeta(3)
            return tmp1 + tmp2 + tmp3
        elif n == 4.0:
            tmp1 = 4 * 14 * sc.zeta(3) * (np.log(2) + np.euler_gamma)
            tmp2 = 3 * np.pi**2 * (np.log(2) + np.euler_gamma)**2
            tmp3 = (np.log(2) + np.euler_gamma)**4
            tmp4 = 7 * np.pi**4 / 4
            return tmp1 + tmp2 + tmp3 + tmp4
        else:
            # return generic for higher moments
            # return rv_continuous._mom1_sc(self, n, b)
            return self._mom1_sc(n)

moyal = moyal_gen(name="moyal")

class nakagami_gen(rv_continuous):
    r"""A Nakagami continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `nakagami` is:

.. math::

f(x, \nu) = \frac{2 \nu^\nu}{\Gamma(\nu)} x^{2\nu-1} \exp(-\nu x^2)

for :math:`x >= 0`, :math:`\nu > 0`.

`nakagami` takes ``nu`` as a shape parameter for :math:`\nu`.

%(after_notes)s

%(example)s

"""
    def _pdf(self, x, nu):
        # nakagami.pdf(x, nu) = 2 * nu**nu / gamma(nu) *
        #                       x**(2*nu-1) * exp(-nu*x**2)
        return 2*nu**nu/sc.gamma(nu)*(x**(2*nu-1.0))*np.exp(-nu*x*x)

def _cdf(self, x, nu):
        return sc.gammainc(nu, nu*x*x)

def _ppf(self, q, nu):
        return np.sqrt(1.0/nu*sc.gammaincinv(nu, q))

def _stats(self, nu):
        mu = sc.gamma(nu+0.5)/sc.gamma(nu)/np.sqrt(nu)
        mu2 = 1.0-mu*mu
        g1 = mu * (1 - 4*nu*mu2) / 2.0 / nu / np.power(mu2, 1.5)
        g2 = -6*mu**4*nu + (8*nu-2)*mu**2-2*nu + 1
        g2 /= nu*mu2**2.0
        return mu, mu2, g1, g2

nakagami = nakagami_gen(a=0.0, name="nakagami")

class ncx2_gen(rv_continuous):
    r"""A non-central chi-squared continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `ncx2` is:

.. math::

f(x, k, \lambda) = \frac{1}{2} \exp(-(\lambda+x)/2)
            (x/\lambda)^{(k-2)/4}  I_{(k-2)/2}(\sqrt{\lambda x})

for :math:`x >= 0` and :math:`k, \lambda > 0`. :math:`k` specifies the
    degrees of freedom (denoted ``df`` in the implementation) and
    :math:`\lambda` is the non-centrality parameter (denoted ``nc`` in the
    implementation). :math:`I_\nu` denotes the modified Bessel function of
    first order of degree :math:`\nu` (`scipy.special.iv`).

`ncx2` takes ``df`` and ``nc`` as shape parameters.

%(after_notes)s

%(example)s

"""
    def _argcheck(self, df, nc):
        return (df > 0) & (nc >= 0)

def _rvs(self, df, nc):
        return self._random_state.noncentral_chisquare(df, nc, self._size)

def _logpdf(self, x, df, nc):
        cond = np.ones_like(x, dtype=bool) & (nc != 0)
        return _lazywhere(cond, (x, df, nc), f=_ncx2_log_pdf, f2=chi2.logpdf)

def _pdf(self, x, df, nc):
        # ncx2.pdf(x, df, nc) = exp(-(nc+x)/2) * 1/2 * (x/nc)**((df-2)/4)
        #                       * I[(df-2)/2](sqrt(nc*x))
        cond = np.ones_like(x, dtype=bool) & (nc != 0)
        return _lazywhere(cond, (x, df, nc), f=_ncx2_pdf, f2=chi2.pdf)

def _cdf(self, x, df, nc):
        cond = np.ones_like(x, dtype=bool) & (nc != 0)
        return _lazywhere(cond, (x, df, nc), f=_ncx2_cdf, f2=chi2.cdf)

def _ppf(self, q, df, nc):
        cond = np.ones_like(q, dtype=bool) & (nc != 0)
        return _lazywhere(cond, (q, df, nc), f=sc.chndtrix, f2=chi2.ppf)

def _stats(self, df, nc):
        val = df + 2.0*nc
        return (df + nc,
                2*val,
                np.sqrt(8)*(val+nc)/val**1.5,
                12.0*(val+2*nc)/val**2.0)

ncx2 = ncx2_gen(a=0.0, name='ncx2')

class ncf_gen(rv_continuous):
    r"""A non-central F distribution continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `ncf` is:

.. math::

f(x, n_1, n_2, \lambda) =
                          \exp(\frac{\lambda}{2} + \lambda n_1 \frac{x}{2(n_1 x+n_2)})
                          n_1^{n_1/2} n_2^{n_2/2} x^{n_1/2 - 1} \\
                          (n_2+n_1 x)^{-(n_1+n_2)/2}
                          \gamma(n_1/2) \gamma(1+n_2/2) \\
                         \frac{L^{\frac{v_1}{2}-1}_{v_2/2}
                               (-\lambda v_1 \frac{x}{2(v_1 x+v_2)})}
                              {B(v_1/2, v_2/2)  \gamma(\frac{v_1+v_2}{2})}

for :math:`n_1 > 1`, :math:`n_2, \lambda > 0`.  Here :math:`n_1` is the
    degrees of freedom in the numerator, :math:`n_2` the degrees of freedom in
    the denominator, :math:`\lambda` the non-centrality parameter,
    :math:`\gamma` is the logarithm of the Gamma function, :math:`L_n^k` is a
    generalized Laguerre polynomial and :math:`B` is the beta function.

`ncf` takes ``df1``, ``df2`` and ``nc`` as shape parameters.

%(after_notes)s

%(example)s

"""
    def _rvs(self, dfn, dfd, nc):
        return self._random_state.noncentral_f(dfn, dfd, nc, self._size)

def _pdf_skip(self, x, dfn, dfd, nc):
        # ncf.pdf(x, df1, df2, nc) = exp(nc/2 + nc*df1*x/(2*(df1*x+df2))) *
        #             df1**(df1/2) * df2**(df2/2) * x**(df1/2-1) *
        #             (df2+df1*x)**(-(df1+df2)/2) *
        #             gamma(df1/2)*gamma(1+df2/2) *
        #             L^{v1/2-1}^{v2/2}(-nc*v1*x/(2*(v1*x+v2))) /
        #             (B(v1/2, v2/2) * gamma((v1+v2)/2))
        n1, n2 = dfn, dfd
        term = -nc/2+nc*n1*x/(2*(n2+n1*x)) + sc.gammaln(n1/2.)+sc.gammaln(1+n2/2.)
        term -= sc.gammaln((n1+n2)/2.0)
        Px = np.exp(term)
        Px *= n1**(n1/2) * n2**(n2/2) * x**(n1/2-1)
        Px *= (n2+n1*x)**(-(n1+n2)/2)
        Px *= sc.assoc_laguerre(-nc*n1*x/(2.0*(n2+n1*x)), n2/2, n1/2-1)
        Px /= sc.beta(n1/2, n2/2)
        # This function does not have a return.  Drop it for now, the generic
        # function seems to work OK.

def _cdf(self, x, dfn, dfd, nc):
        return sc.ncfdtr(dfn, dfd, nc, x)

def _ppf(self, q, dfn, dfd, nc):
        return sc.ncfdtri(dfn, dfd, nc, q)

def _munp(self, n, dfn, dfd, nc):
        val = (dfn * 1.0/dfd)**n
        term = sc.gammaln(n+0.5*dfn) + sc.gammaln(0.5*dfd-n) - sc.gammaln(dfd*0.5)
        val *= np.exp(-nc / 2.0+term)
        val *= sc.hyp1f1(n+0.5*dfn, 0.5*dfn, 0.5*nc)
        return val

def _stats(self, dfn, dfd, nc):
        # Note: the rv_continuous class ensures that dfn > 0 when this function
        # is called, so we don't have  to check for division by zero with dfn
        # in the following.
        mu_num = dfd * (dfn + nc)
        mu_den = dfn * (dfd - 2)
        mu = np.full_like(mu_num, dtype=np.float64, fill_value=np.inf)
        np.true_divide(mu_num, mu_den, where=dfd > 2, out=mu)

mu2_num = 2*((dfn + nc)**2 + (dfn + 2*nc)*(dfd - 2))*(dfd/dfn)**2
        mu2_den = (dfd - 2)**2 * (dfd - 4)
        mu2 = np.full_like(mu2_num, dtype=np.float64, fill_value=np.inf)
        np.true_divide(mu2_num, mu2_den, where=dfd > 4, out=mu2)

return mu, mu2, None, None

ncf = ncf_gen(a=0.0, name='ncf')

class t_gen(rv_continuous):
    r"""A Student's t continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `t` is:

.. math::

f(x, \nu) = \frac{\Gamma((\nu+1)/2)}
                        {\sqrt{\pi \nu} \Gamma(\nu/2)}
                    (1+x^2/\nu)^{-(\nu+1)/2}

where :math:`x` is a real number and the degrees of freedom parameter
    :math:`\nu` (denoted ``df`` in the implementation) satisfies
    :math:`\nu > 0`. :math:`\Gamma` is the gamma function
    (`scipy.special.gamma`).

%(after_notes)s

%(example)s

"""
    def _argcheck(self, df):
        return df > 0

def _rvs(self, df):
        return self._random_state.standard_t(df, size=self._size)

def _pdf(self, x, df):
        #                                gamma((df+1)/2)
        # t.pdf(x, df) = ---------------------------------------------------
        #                sqrt(pi*df) * gamma(df/2) * (1+x**2/df)**((df+1)/2)
        r = np.asarray(df*1.0)
        Px = np.exp(sc.gammaln((r+1)/2)-sc.gammaln(r/2))
        Px /= np.sqrt(r*np.pi)*(1+(x**2)/r)**((r+1)/2)
        return Px

def _logpdf(self, x, df):
        r = df*1.0
        lPx = sc.gammaln((r+1)/2)-sc.gammaln(r/2)
        lPx -= 0.5*np.log(r*np.pi) + (r+1)/2*np.log(1+(x**2)/r)
        return lPx

def _cdf(self, x, df):
        return sc.stdtr(df, x)

def _sf(self, x, df):
        return sc.stdtr(df, -x)

def _ppf(self, q, df):
        return sc.stdtrit(df, q)

def _isf(self, q, df):
        return -sc.stdtrit(df, q)

def _stats(self, df):
        mu = np.where(df > 1, 0.0, np.inf)
        mu2 = _lazywhere(df > 2, (df,),
                         lambda df: df / (df-2.0),
                         np.inf)
        mu2 = np.where(df <= 1, np.nan, mu2)
        g1 = np.where(df > 3, 0.0, np.nan)
        g2 = _lazywhere(df > 4, (df,),
                        lambda df: 6.0 / (df-4.0),
                        np.inf)
        g2 = np.where(df <= 2, np.nan, g2)
        return mu, mu2, g1, g2

t = t_gen(name='t')

class nct_gen(rv_continuous):
    r"""A non-central Student's t continuous random variable.

%(before_notes)s

Notes
    -----
    If :math:`Y` is a standard normal random variable and :math:`V` is
    an independent chi-square random variable (`chi2`) with :math:`k` degrees
    of freedom, then

.. math::

X = \frac{Y + c}{\sqrt{V/k}}

has a non-central Student's t distribution on the real line.
    The degrees of freedom parameter :math:`k` (denoted ``df`` in the
    implementation) satisfies :math:`k > 0` and the noncentrality parameter
    :math:`c` (denoted ``nc`` in the implementation) is a real number.

%(after_notes)s

%(example)s

"""
    def _argcheck(self, df, nc):
        return (df > 0) & (nc == nc)

def _rvs(self, df, nc):
        sz, rndm = self._size, self._random_state
        n = norm.rvs(loc=nc, size=sz, random_state=rndm)
        c2 = chi2.rvs(df, size=sz, random_state=rndm)
        return n * np.sqrt(df) / np.sqrt(c2)

def _pdf(self, x, df, nc):
        n = df*1.0
        nc = nc*1.0
        x2 = x*x
        ncx2 = nc*nc*x2
        fac1 = n + x2
        trm1 = n/2.*np.log(n) + sc.gammaln(n+1)
        trm1 -= n*np.log(2)+nc*nc/2.+(n/2.)*np.log(fac1)+sc.gammaln(n/2.)
        Px = np.exp(trm1)
        valF = ncx2 / (2*fac1)
        trm1 = np.sqrt(2)*nc*x*sc.hyp1f1(n/2+1, 1.5, valF)
        trm1 /= np.asarray(fac1*sc.gamma((n+1)/2))
        trm2 = sc.hyp1f1((n+1)/2, 0.5, valF)
        trm2 /= np.asarray(np.sqrt(fac1)*sc.gamma(n/2+1))
        Px *= trm1+trm2
        return Px

def _cdf(self, x, df, nc):
        return sc.nctdtr(df, nc, x)

def _ppf(self, q, df, nc):
        return sc.nctdtrit(df, nc, q)

def _stats(self, df, nc, moments='mv'):
        #
        # See D. Hogben, R.S. Pinkham, and M.B. Wilk,
        # 'The moments of the non-central t-distribution'
        # Biometrika 48, p. 465 (2961).
        # e.g. https://www.jstor.org/stable/2332772 (gated)
        #
        mu, mu2, g1, g2 = None, None, None, None

gfac = sc.gamma(df/2.-0.5) / sc.gamma(df/2.)
        c11 = np.sqrt(df/2.) * gfac
        c20 = df / (df-2.)
        c22 = c20 - c11*c11
        mu = np.where(df > 1, nc*c11, np.inf)
        mu2 = np.where(df > 2, c22*nc*nc + c20, np.inf)
        if 's' in moments:
            c33t = df * (7.-2.*df) / (df-2.) / (df-3.) + 2.*c11*c11
            c31t = 3.*df / (df-2.) / (df-3.)
            mu3 = (c33t*nc*nc + c31t) * c11*nc
            g1 = np.where(df > 3, mu3 / np.power(mu2, 1.5), np.nan)
        # kurtosis
        if 'k' in moments:
            c44 = df*df / (df-2.) / (df-4.)
            c44 -= c11*c11 * 2.*df*(5.-df) / (df-2.) / (df-3.)
            c44 -= 3.*c11**4
            c42 = df / (df-4.) - c11*c11 * (df-1.) / (df-3.)
            c42 *= 6.*df / (df-2.)
            c40 = 3.*df*df / (df-2.) / (df-4.)

mu4 = c44 * nc**4 + c42*nc**2 + c40
            g2 = np.where(df > 4, mu4/mu2**2 - 3., np.nan)
        return mu, mu2, g1, g2

nct = nct_gen(name="nct")

class pareto_gen(rv_continuous):
    r"""A Pareto continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `pareto` is:

.. math::

f(x, b) = \frac{b}{x^{b+1}}

for :math:`x \ge 1`, :math:`b > 0`.

`pareto` takes ``b`` as a shape parameter for :math:`b`.

%(after_notes)s

%(example)s

"""
    def _pdf(self, x, b):
        # pareto.pdf(x, b) = b / x**(b+1)
        return b * x**(-b-1)

def _cdf(self, x, b):
        return 1 - x**(-b)

def _ppf(self, q, b):
        return pow(1-q, -1.0/b)

def _sf(self, x, b):
        return x**(-b)

def _stats(self, b, moments='mv'):
        mu, mu2, g1, g2 = None, None, None, None
        if 'm' in moments:
            mask = b > 1
            bt = np.extract(mask, b)
            mu = valarray(np.shape(b), value=np.inf)
            np.place(mu, mask, bt / (bt-1.0))
        if 'v' in moments:
            mask = b > 2
            bt = np.extract(mask, b)
            mu2 = valarray(np.shape(b), value=np.inf)
            np.place(mu2, mask, bt / (bt-2.0) / (bt-1.0)**2)
        if 's' in moments:
            mask = b > 3
            bt = np.extract(mask, b)
            g1 = valarray(np.shape(b), value=np.nan)
            vals = 2 * (bt + 1.0) * np.sqrt(bt - 2.0) / ((bt - 3.0) * np.sqrt(bt))
            np.place(g1, mask, vals)
        if 'k' in moments:
            mask = b > 4
            bt = np.extract(mask, b)
            g2 = valarray(np.shape(b), value=np.nan)
            vals = (6.0*np.polyval([1.0, 1.0, -6, -2], bt) /
                    np.polyval([1.0, -7.0, 12.0, 0.0], bt))
            np.place(g2, mask, vals)
        return mu, mu2, g1, g2

def _entropy(self, c):
        return 1 + 1.0/c - np.log(c)

pareto = pareto_gen(a=1.0, name="pareto")

class lomax_gen(rv_continuous):
    r"""A Lomax (Pareto of the second kind) continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `lomax` is:

.. math::

f(x, c) = \frac{c}{(1+x)^{c+1}}

for :math:`x \ge 0`, :math:`c > 0`.

`lomax` takes ``c`` as a shape parameter for :math:`c`.

`lomax` is a special case of `pareto` with ``loc=-1.0``.

%(after_notes)s

%(example)s

"""
    def _pdf(self, x, c):
        # lomax.pdf(x, c) = c / (1+x)**(c+1)
        return c*1.0/(1.0+x)**(c+1.0)

def _logpdf(self, x, c):
        return np.log(c) - (c+1)*sc.log1p(x)

def _cdf(self, x, c):
        return -sc.expm1(-c*sc.log1p(x))

def _sf(self, x, c):
        return np.exp(-c*sc.log1p(x))

def _logsf(self, x, c):
        return -c*sc.log1p(x)

def _ppf(self, q, c):
        return sc.expm1(-sc.log1p(-q)/c)

def _stats(self, c):
        mu, mu2, g1, g2 = pareto.stats(c, loc=-1.0, moments='mvsk')
        return mu, mu2, g1, g2

def _entropy(self, c):
        return 1+1.0/c-np.log(c)

lomax = lomax_gen(a=0.0, name="lomax")

class pearson3_gen(rv_continuous):
    r"""A pearson type III continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `pearson3` is:

.. math::

f(x, skew) = \frac{|\beta|}{\Gamma(\alpha)}
                     (\beta (x - \zeta))^{\alpha - 1}
                     \exp(-\beta (x - \zeta))

where:

.. math::

\beta = \frac{2}{skew  stddev}
            \alpha = (stddev \beta)^2
            \zeta = loc - \frac{\alpha}{\beta}

:math:`\Gamma` is the gamma function (`scipy.special.gamma`).
    `pearson3` takes ``skew`` as a shape parameter for :math:`skew`.

%(after_notes)s

%(example)s

References
    ----------
    R.W. Vogel and D.E. McMartin, "Probability Plot Goodness-of-Fit and
    Skewness Estimation Procedures for the Pearson Type 3 Distribution", Water
    Resources Research, Vol.27, 3149-3158 (1991).

L.R. Salvosa, "Tables of Pearson's Type III Function", Ann. Math. Statist.,
    Vol.1, 191-198 (1930).

"Using Modern Computing Tools to Fit the Pearson Type III Distribution to
    Aviation Loads Data", Office of Aviation Research (2003).

"""
    def _preprocess(self, x, skew):
        # The real 'loc' and 'scale' are handled in the calling pdf(...). The
        # local variables 'loc' and 'scale' within pearson3._pdf are set to
        # the defaults just to keep them as part of the equations for
        # documentation.
        loc = 0.0
        scale = 1.0

# If skew is small, return _norm_pdf. The divide between pearson3
        # and norm was found by brute force and is approximately a skew of
        # 0.000016.  No one, I hope, would actually use a skew value even
        # close to this small.
        norm2pearson_transition = 0.000016

ans, x, skew = np.broadcast_arrays([1.0], x, skew)
        ans = ans.copy()

# mask is True where skew is small enough to use the normal approx.
        mask = np.absolute(skew) < norm2pearson_transition
        invmask = ~mask

beta = 2.0 / (skew[invmask] * scale)
        alpha = (scale * beta)**2
        zeta = loc - alpha / beta

transx = beta * (x[invmask] - zeta)
        return ans, x, transx, mask, invmask, beta, alpha, zeta

def _argcheck(self, skew):
        # The _argcheck function in rv_continuous only allows positive
        # arguments.  The skew argument for pearson3 can be zero (which I want
        # to handle inside pearson3._pdf) or negative.  So just return True
        # for all skew args.
        return np.ones(np.shape(skew), dtype=bool)

def _stats(self, skew):
        _, _, _, _, _, beta, alpha, zeta = (
            self._preprocess([1], skew))
        m = zeta + alpha / beta
        v = alpha / (beta**2)
        s = 2.0 / (alpha**0.5) * np.sign(beta)
        k = 6.0 / alpha
        return m, v, s, k

def _pdf(self, x, skew):
        # pearson3.pdf(x, skew) = abs(beta) / gamma(alpha) *
        #     (beta * (x - zeta))**(alpha - 1) * exp(-beta*(x - zeta))
        # Do the calculation in _logpdf since helps to limit
        # overflow/underflow problems
        ans = np.exp(self._logpdf(x, skew))
        if ans.ndim == 0:
            if np.isnan(ans):
                return 0.0
            return ans
        ans[np.isnan(ans)] = 0.0
        return ans

def _logpdf(self, x, skew):
        #   PEARSON3 logpdf                           GAMMA logpdf
        #   np.log(abs(beta))
        # + (alpha - 1)*np.log(beta*(x - zeta))          + (a - 1)*np.log(x)
        # - beta*(x - zeta)                           - x
        # - sc.gammalnalpha)                              - sc.gammalna)
        ans, x, transx, mask, invmask, beta, alpha, _ = (
            self._preprocess(x, skew))

ans[mask] = np.log(_norm_pdf(x[mask]))
        ans[invmask] = np.log(abs(beta)) + gamma._logpdf(transx, alpha)
        return ans

def _cdf(self, x, skew):
        ans, x, transx, mask, invmask, _, alpha, _ = (
            self._preprocess(x, skew))

ans[mask] = _norm_cdf(x[mask])
        ans[invmask] = gamma._cdf(transx, alpha)
        return ans

def _rvs(self, skew):
        skew = broadcast_to(skew, self._size)
        ans, _, _, mask, invmask, beta, alpha, zeta = (
            self._preprocess([0], skew))

nsmall = mask.sum()
        nbig = mask.size - nsmall
        ans[mask] = self._random_state.standard_normal(nsmall)
        ans[invmask] = (self._random_state.standard_gamma(alpha, nbig)/beta +
                        zeta)

if self._size == ():
            ans = ans[0]
        return ans

def _ppf(self, q, skew):
        ans, q, _, mask, invmask, beta, alpha, zeta = (
            self._preprocess(q, skew))
        ans[mask] = _norm_ppf(q[mask])
        ans[invmask] = sc.gammaincinv(alpha, q[invmask])/beta + zeta
        return ans

pearson3 = pearson3_gen(name="pearson3")

class powerlaw_gen(rv_continuous):
    r"""A power-function continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `powerlaw` is:

.. math::

f(x, a) = a x^{a-1}

for :math:`0 \le x \le 1`, :math:`a > 0`.

`powerlaw` takes ``a`` as a shape parameter for :math:`a`.

%(after_notes)s

`powerlaw` is a special case of `beta` with ``b=1``.

%(example)s

"""
    def _pdf(self, x, a):
        # powerlaw.pdf(x, a) = a * x**(a-1)
        return a*x**(a-1.0)

def _logpdf(self, x, a):
        return np.log(a) + sc.xlogy(a - 1, x)

def _cdf(self, x, a):
        return x**(a*1.0)

def _logcdf(self, x, a):
        return a*np.log(x)

def _ppf(self, q, a):
        return pow(q, 1.0/a)

def _stats(self, a):
        return (a / (a + 1.0),
                a / (a + 2.0) / (a + 1.0) ** 2,
                -2.0 * ((a - 1.0) / (a + 3.0)) * np.sqrt((a + 2.0) / a),
                6 * np.polyval([1, -1, -6, 2], a) / (a * (a + 3.0) * (a + 4)))

def _entropy(self, a):
        return 1 - 1.0/a - np.log(a)

powerlaw = powerlaw_gen(a=0.0, b=1.0, name="powerlaw")

class powerlognorm_gen(rv_continuous):
    r"""A power log-normal continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `powerlognorm` is:

.. math::

f(x, c, s) = \frac{c}{x s} \phi(\log(x)/s)
                     (\Phi(-\log(x)/s))^{c-1}

where :math:`\phi` is the normal pdf, and :math:`\Phi` is the normal cdf,
    and :math:`x > 0`, :math:`s, c > 0`.

`powerlognorm` takes :math:`c` and :math:`s` as shape parameters.

%(after_notes)s

%(example)s

"""
    _support_mask = rv_continuous._open_support_mask

def _pdf(self, x, c, s):
        # powerlognorm.pdf(x, c, s) = c / (x*s) * phi(log(x)/s) *
        #                                         (Phi(-log(x)/s))**(c-1),
        return (c/(x*s) * _norm_pdf(np.log(x)/s) *
                pow(_norm_cdf(-np.log(x)/s), c*1.0-1.0))

def _cdf(self, x, c, s):
        return 1.0 - pow(_norm_cdf(-np.log(x)/s), c*1.0)

def _ppf(self, q, c, s):
        return np.exp(-s * _norm_ppf(pow(1.0 - q, 1.0 / c)))

powerlognorm = powerlognorm_gen(a=0.0, name="powerlognorm")

class powernorm_gen(rv_continuous):
    r"""A power normal continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `powernorm` is:

.. math::

f(x, c) = c \phi(x) (\Phi(-x))^{c-1}

where :math:`\phi` is the normal pdf, and :math:`\Phi` is the normal cdf,
    and :math:`x >= 0`, :math:`c > 0`.

`powernorm` takes ``c`` as a shape parameter for :math:`c`.

%(after_notes)s

%(example)s

"""
    def _pdf(self, x, c):
        # powernorm.pdf(x, c) = c * phi(x) * (Phi(-x))**(c-1)
        return c*_norm_pdf(x) * (_norm_cdf(-x)**(c-1.0))

def _logpdf(self, x, c):
        return np.log(c) + _norm_logpdf(x) + (c-1)*_norm_logcdf(-x)

def _cdf(self, x, c):
        return 1.0-_norm_cdf(-x)**(c*1.0)

def _ppf(self, q, c):
        return -_norm_ppf(pow(1.0 - q, 1.0 / c))

powernorm = powernorm_gen(name='powernorm')

class rdist_gen(rv_continuous):
    r"""An R-distributed (symmetric beta) continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `rdist` is:

.. math::

f(x, c) = \frac{(1-x^2)^{c/2-1}}{B(1/2, c/2)}

for :math:`-1 \le x \le 1`, :math:`c > 0`. `rdist` is also called the
    symmetric beta distribution: if B has a `beta` distribution with
    parameters (c/2, c/2), then X = 2*B - 1 follows a R-distribution with
    parameter c.

`rdist` takes ``c`` as a shape parameter for :math:`c`.

This distribution includes the following distribution kernels as
    special cases::

c = 2:  uniform
        c = 3:  `semicircular`
        c = 4:  Epanechnikov (parabolic)
        c = 6:  quartic (biweight)
        c = 8:  triweight

%(after_notes)s

%(example)s

"""
    # use relation to the beta distribution for pdf, cdf, etc
    def _pdf(self, x, c):
        return 0.5*beta._pdf((x + 1)/2, c/2, c/2)

def _logpdf(self, x, c):
        return -np.log(2) + beta._logpdf((x + 1)/2, c/2, c/2)

def _cdf(self, x, c):
        return beta._cdf((x + 1)/2, c/2, c/2)

def _ppf(self, q, c):
        return 2*beta._ppf(q, c/2, c/2) - 1

def _rvs(self, c):
        return 2 * self._random_state.beta(c/2, c/2, self._size) - 1

def _munp(self, n, c):
        numerator = (1 - (n % 2)) * sc.beta((n + 1.0) / 2, c / 2.0)
        return numerator / sc.beta(1. / 2, c / 2.)

rdist = rdist_gen(a=-1.0, b=1.0, name="rdist")

class rayleigh_gen(rv_continuous):
    r"""A Rayleigh continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `rayleigh` is:

.. math::

f(x) = x \exp(-x^2/2)

for :math:`x \ge 0`.

`rayleigh` is a special case of `chi` with ``df=2``.

%(after_notes)s

%(example)s

"""
    _support_mask = rv_continuous._open_support_mask

def _rvs(self):
        return chi.rvs(2, size=self._size, random_state=self._random_state)

def _pdf(self, r):
        # rayleigh.pdf(r) = r * exp(-r**2/2)
        return np.exp(self._logpdf(r))

def _logpdf(self, r):
        return np.log(r) - 0.5 * r * r

def _cdf(self, r):
        return -sc.expm1(-0.5 * r**2)

def _ppf(self, q):
        return np.sqrt(-2 * sc.log1p(-q))

def _sf(self, r):
        return np.exp(self._logsf(r))

def _logsf(self, r):
        return -0.5 * r * r

def _isf(self, q):
        return np.sqrt(-2 * np.log(q))

def _stats(self):
        val = 4 - np.pi
        return (np.sqrt(np.pi/2),
                val/2,
                2*(np.pi-3)*np.sqrt(np.pi)/val**1.5,
                6*np.pi/val-16/val**2)

def _entropy(self):
        return _EULER/2.0 + 1 - 0.5*np.log(2)

rayleigh = rayleigh_gen(a=0.0, name="rayleigh")

class reciprocal_gen(rv_continuous):
    r"""A loguniform or reciprocal continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for this class is:

.. math::

f(x, a, b) = \frac{1}{x \log(b/a)}

for :math:`a \le x \le b`, :math:`b > a > 0`. This class takes
    :math:`a` and :math:`b` as shape parameters. %(after_notes)s

%(example)s

This doesn't show the equal probability of ``0.01``, ``0.1`` and
    ``1``. This is best when the x-axis is log-scaled:

>>> import numpy as np
    >>> fig, ax = plt.subplots(1, 1)
    >>> ax.hist(np.log10(r))
    >>> ax.set_ylabel("Frequency")
    >>> ax.set_xlabel("Value of random variable")
    >>> ax.xaxis.set_major_locator(plt.FixedLocator([-2, -1, 0]))
    >>> ticks = ["$10^{{ {} }}$".format(i) for i in [-2, -1, 0]]
    >>> ax.set_xticklabels(ticks)  # doctest: +SKIP
    >>> plt.show()

This random variable will be log-uniform regardless of the base chosen for
    ``a`` and ``b``. Let's specify with base ``2`` instead:

>>> rvs = %(name)s(2**-2, 2**0).rvs(size=1000)

Values of ``1/4``, ``1/2`` and ``1`` are equally likely with this random
    variable.  Here's the histogram:

>>> fig, ax = plt.subplots(1, 1)
    >>> ax.hist(np.log2(rvs))
    >>> ax.set_ylabel("Frequency")
    >>> ax.set_xlabel("Value of random variable")
    >>> ax.xaxis.set_major_locator(plt.FixedLocator([-2, -1, 0]))
    >>> ticks = ["$2^{{ {} }}$".format(i) for i in [-2, -1, 0]]
    >>> ax.set_xticklabels(ticks)  # doctest: +SKIP
    >>> plt.show()

"""
    def _argcheck(self, a, b):
        return (a > 0) & (b > a)

def _get_support(self, a, b):
        return a, b

def _pdf(self, x, a, b):
        # reciprocal.pdf(x, a, b) = 1 / (x*log(b/a))
        return 1.0 / (x * np.log(b * 1.0 / a))

def _logpdf(self, x, a, b):
        return -np.log(x) - np.log(np.log(b * 1.0 / a))

def _cdf(self, x, a, b):
        return (np.log(x)-np.log(a)) / np.log(b * 1.0 / a)

def _ppf(self, q, a, b):
        return a*pow(b*1.0/a, q)

def _munp(self, n, a, b):
        return 1.0/np.log(b*1.0/a) / n * (pow(b*1.0, n) - pow(a*1.0, n))

def _entropy(self, a, b):
        return 0.5*np.log(a*b)+np.log(np.log(b*1.0/a))

loguniform = reciprocal_gen(name="loguniform")
reciprocal = reciprocal_gen(name="reciprocal")

class rice_gen(rv_continuous):
    r"""A Rice continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `rice` is:

.. math::

f(x, b) = x \exp(- \frac{x^2 + b^2}{2}) I_0(x b)

for :math:`x >= 0`, :math:`b > 0`. :math:`I_0` is the modified Bessel
    function of order zero (`scipy.special.i0`).

`rice` takes ``b`` as a shape parameter for :math:`b`.

%(after_notes)s

The Rice distribution describes the length, :math:`r`, of a 2-D vector with
    components :math:`(U+u, V+v)`, where :math:`U, V` are constant, :math:`u,
    v` are independent Gaussian random variables with standard deviation
    :math:`s`.  Let :math:`R = \sqrt{U^2 + V^2}`. Then the pdf of :math:`r` is
    ``rice.pdf(x, R/s, scale=s)``.

%(example)s

"""
    def _argcheck(self, b):
        return b >= 0

def _rvs(self, b):
        # https://en.wikipedia.org/wiki/Rice_distribution
        t = b/np.sqrt(2) + self._random_state.standard_normal(size=(2,) +
                                                              self._size)
        return np.sqrt((t*t).sum(axis=0))

def _cdf(self, x, b):
        return sc.chndtr(np.square(x), 2, np.square(b))

def _ppf(self, q, b):
        return np.sqrt(sc.chndtrix(q, 2, np.square(b)))

def _pdf(self, x, b):
        # rice.pdf(x, b) = x * exp(-(x**2+b**2)/2) * I[0](x*b)
        #
        # We use (x**2 + b**2)/2 = ((x-b)**2)/2 + xb.
        # The factor of np.exp(-xb) is then included in the i0e function
        # in place of the modified Bessel function, i0, improving
        # numerical stability for large values of xb.
        return x * np.exp(-(x-b)*(x-b)/2.0) * sc.i0e(x*b)

def _munp(self, n, b):
        nd2 = n/2.0
        n1 = 1 + nd2
        b2 = b*b/2.0
        return (2.0**(nd2) * np.exp(-b2) * sc.gamma(n1) *
                sc.hyp1f1(n1, 1, b2))

rice = rice_gen(a=0.0, name="rice")

# FIXME: PPF does not work.
class recipinvgauss_gen(rv_continuous):
    r"""A reciprocal inverse Gaussian continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `recipinvgauss` is:

.. math::

f(x, \mu) = \frac{1}{\sqrt{2\pi x}}
                    \exp\left(\frac{-(1-\mu x)^2}{2\mu^2x}\right)

for :math:`x \ge 0`.

`recipinvgauss` takes ``mu`` as a shape parameter for :math:`\mu`.

%(after_notes)s

%(example)s

"""

def _pdf(self, x, mu):
        # recipinvgauss.pdf(x, mu) =
        #                     1/sqrt(2*pi*x) * exp(-(1-mu*x)**2/(2*x*mu**2))
        return 1.0/np.sqrt(2*np.pi*x)*np.exp(-(1-mu*x)**2.0 / (2*x*mu**2.0))

def _logpdf(self, x, mu):
        return -(1-mu*x)**2.0 / (2*x*mu**2.0) - 0.5*np.log(2*np.pi*x)

def _cdf(self, x, mu):
        trm1 = 1.0/mu - x
        trm2 = 1.0/mu + x
        isqx = 1.0/np.sqrt(x)
        return 1.0-_norm_cdf(isqx*trm1)-np.exp(2.0/mu)*_norm_cdf(-isqx*trm2)

def _rvs(self, mu):
        return 1.0/self._random_state.wald(mu, 1.0, size=self._size)

recipinvgauss = recipinvgauss_gen(a=0.0, name='recipinvgauss')

class semicircular_gen(rv_continuous):
    r"""A semicircular continuous random variable.

%(before_notes)s

Notes
    -----
    The probability density function for `semicircular` is:

.. math::

f(x) = \frac{2}{\pi} \sqrt{1-x^2}

for :math:`-1 \le x \le 1`.

The distribution is a special case of `rdist` with `c = 3`.

%(after_notes)s