PageRenderTime 113ms CodeModel.GetById 47ms app.highlight 55ms RepoModel.GetById 0ms app.codeStats 0ms

/scipy/stats/_continuous_distns.py

https://github.com/matthew-brett/scipy
Python | 1846 lines | 1832 code | 7 blank | 7 comment | 6 complexity | 790fa5e2faa5ee1f6aa083fd4f08f4ae MD5 | raw file
   1# -*- coding: utf-8 -*-
   2#
   3# Author:  Travis Oliphant  2002-2011 with contributions from
   4#          SciPy Developers 2004-2011
   5#
   6import warnings
   7from collections.abc import Iterable
   8from functools import wraps
   9import ctypes
  10
  11import numpy as np
  12
  13from scipy._lib.doccer import (extend_notes_in_docstring,
  14                               replace_notes_in_docstring)
  15from scipy._lib._ccallback import LowLevelCallable
  16from scipy import optimize
  17from scipy import integrate
  18from scipy import interpolate
  19import scipy.special as sc
  20
  21import scipy.special._ufuncs as scu
  22from scipy._lib._util import _lazyselect, _lazywhere
  23from . import _stats
  24from ._rvs_sampling import rvs_ratio_uniforms
  25from ._tukeylambda_stats import (tukeylambda_variance as _tlvar,
  26                                 tukeylambda_kurtosis as _tlkurt)
  27from ._distn_infrastructure import (
  28    get_distribution_names, _kurtosis, _ncx2_cdf, _ncx2_log_pdf, _ncx2_pdf,
  29    rv_continuous, _skew, _get_fixed_fit_value, _check_shape)
  30from ._ksstats import kolmogn, kolmognp, kolmogni
  31from ._constants import (_XMIN, _EULER, _ZETA3,
  32                         _SQRT_2_OVER_PI, _LOG_SQRT_2_OVER_PI)
  33import scipy.stats._boost as _boost
  34
  35
  36def _remove_optimizer_parameters(kwds):
  37    """
  38    Remove the optimizer-related keyword arguments 'loc', 'scale' and
  39    'optimizer' from `kwds`.  Then check that `kwds` is empty, and
  40    raise `TypeError("Unknown arguments: %s." % kwds)` if it is not.
  41
  42    This function is used in the fit method of distributions that override
  43    the default method and do not use the default optimization code.
  44
  45    `kwds` is modified in-place.
  46    """
  47    kwds.pop('loc', None)
  48    kwds.pop('scale', None)
  49    kwds.pop('optimizer', None)
  50    kwds.pop('method', None)
  51    if kwds:
  52        raise TypeError("Unknown arguments: %s." % kwds)
  53
  54
  55def _call_super_mom(fun):
  56    # if fit method is overridden only for MLE and doesn't specify what to do
  57    # if method == 'mm', this decorator calls generic implementation
  58    @wraps(fun)
  59    def wrapper(self, *args, **kwds):
  60        method = kwds.get('method', 'mle').lower()
  61        if method == 'mm':
  62            return super(type(self), self).fit(*args, **kwds)
  63        else:
  64            return fun(self, *args, **kwds)
  65    return wrapper
  66
  67
  68class ksone_gen(rv_continuous):
  69    r"""Kolmogorov-Smirnov one-sided test statistic distribution.
  70
  71    This is the distribution of the one-sided Kolmogorov-Smirnov (KS)
  72    statistics :math:`D_n^+` and :math:`D_n^-`
  73    for a finite sample size ``n`` (the shape parameter).
  74
  75    %(before_notes)s
  76
  77    See Also
  78    --------
  79    kstwobign, kstwo, kstest
  80
  81    Notes
  82    -----
  83    :math:`D_n^+` and :math:`D_n^-` are given by
  84
  85    .. math::
  86
  87        D_n^+ &= \text{sup}_x (F_n(x) - F(x)),\\
  88        D_n^- &= \text{sup}_x (F(x) - F_n(x)),\\
  89
  90    where :math:`F` is a continuous CDF and :math:`F_n` is an empirical CDF.
  91    `ksone` describes the distribution under the null hypothesis of the KS test
  92    that the empirical CDF corresponds to :math:`n` i.i.d. random variates
  93    with CDF :math:`F`.
  94
  95    %(after_notes)s
  96
  97    References
  98    ----------
  99    .. [1] Birnbaum, Z. W. and Tingey, F.H. "One-sided confidence contours
 100       for probability distribution functions", The Annals of Mathematical
 101       Statistics, 22(4), pp 592-596 (1951).
 102
 103    %(example)s
 104
 105    """
 106    def _pdf(self, x, n):
 107        return -scu._smirnovp(n, x)
 108
 109    def _cdf(self, x, n):
 110        return scu._smirnovc(n, x)
 111
 112    def _sf(self, x, n):
 113        return sc.smirnov(n, x)
 114
 115    def _ppf(self, q, n):
 116        return scu._smirnovci(n, q)
 117
 118    def _isf(self, q, n):
 119        return sc.smirnovi(n, q)
 120
 121
 122ksone = ksone_gen(a=0.0, b=1.0, name='ksone')
 123
 124
 125class kstwo_gen(rv_continuous):
 126    r"""Kolmogorov-Smirnov two-sided test statistic distribution.
 127
 128    This is the distribution of the two-sided Kolmogorov-Smirnov (KS)
 129    statistic :math:`D_n` for a finite sample size ``n``
 130    (the shape parameter).
 131
 132    %(before_notes)s
 133
 134    See Also
 135    --------
 136    kstwobign, ksone, kstest
 137
 138    Notes
 139    -----
 140    :math:`D_n` is given by
 141
 142    .. math::
 143
 144        D_n = \text{sup}_x |F_n(x) - F(x)|
 145
 146    where :math:`F` is a (continuous) CDF and :math:`F_n` is an empirical CDF.
 147    `kstwo` describes the distribution under the null hypothesis of the KS test
 148    that the empirical CDF corresponds to :math:`n` i.i.d. random variates
 149    with CDF :math:`F`.
 150
 151    %(after_notes)s
 152
 153    References
 154    ----------
 155    .. [1] Simard, R., L'Ecuyer, P. "Computing the Two-Sided
 156       Kolmogorov-Smirnov Distribution",  Journal of Statistical Software,
 157       Vol 39, 11, 1-18 (2011).
 158
 159    %(example)s
 160
 161    """
 162    def _get_support(self, n):
 163        return (0.5/(n if not isinstance(n, Iterable) else np.asanyarray(n)),
 164                1.0)
 165
 166    def _pdf(self, x, n):
 167        return kolmognp(n, x)
 168
 169    def _cdf(self, x, n):
 170        return kolmogn(n, x)
 171
 172    def _sf(self, x, n):
 173        return kolmogn(n, x, cdf=False)
 174
 175    def _ppf(self, q, n):
 176        return kolmogni(n, q, cdf=True)
 177
 178    def _isf(self, q, n):
 179        return kolmogni(n, q, cdf=False)
 180
 181
 182# Use the pdf, (not the ppf) to compute moments
 183kstwo = kstwo_gen(momtype=0, a=0.0, b=1.0, name='kstwo')
 184
 185
 186class kstwobign_gen(rv_continuous):
 187    r"""Limiting distribution of scaled Kolmogorov-Smirnov two-sided test statistic.
 188
 189    This is the asymptotic distribution of the two-sided Kolmogorov-Smirnov
 190    statistic :math:`\sqrt{n} D_n` that measures the maximum absolute
 191    distance of the theoretical (continuous) CDF from the empirical CDF.
 192    (see `kstest`).
 193
 194    %(before_notes)s
 195
 196    See Also
 197    --------
 198    ksone, kstwo, kstest
 199
 200    Notes
 201    -----
 202    :math:`\sqrt{n} D_n` is given by
 203
 204    .. math::
 205
 206        D_n = \text{sup}_x |F_n(x) - F(x)|
 207
 208    where :math:`F` is a continuous CDF and :math:`F_n` is an empirical CDF.
 209    `kstwobign`  describes the asymptotic distribution (i.e. the limit of
 210    :math:`\sqrt{n} D_n`) under the null hypothesis of the KS test that the
 211    empirical CDF corresponds to i.i.d. random variates with CDF :math:`F`.
 212
 213    %(after_notes)s
 214
 215    References
 216    ----------
 217    .. [1] Feller, W. "On the Kolmogorov-Smirnov Limit Theorems for Empirical
 218       Distributions",  Ann. Math. Statist. Vol 19, 177-189 (1948).
 219
 220    %(example)s
 221
 222    """
 223    def _pdf(self, x):
 224        return -scu._kolmogp(x)
 225
 226    def _cdf(self, x):
 227        return scu._kolmogc(x)
 228
 229    def _sf(self, x):
 230        return sc.kolmogorov(x)
 231
 232    def _ppf(self, q):
 233        return scu._kolmogci(q)
 234
 235    def _isf(self, q):
 236        return sc.kolmogi(q)
 237
 238
 239kstwobign = kstwobign_gen(a=0.0, name='kstwobign')
 240
 241
 242## Normal distribution
 243
 244# loc = mu, scale = std
 245# Keep these implementations out of the class definition so they can be reused
 246# by other distributions.
 247_norm_pdf_C = np.sqrt(2*np.pi)
 248_norm_pdf_logC = np.log(_norm_pdf_C)
 249
 250
 251def _norm_pdf(x):
 252    return np.exp(-x**2/2.0) / _norm_pdf_C
 253
 254
 255def _norm_logpdf(x):
 256    return -x**2 / 2.0 - _norm_pdf_logC
 257
 258
 259def _norm_cdf(x):
 260    return sc.ndtr(x)
 261
 262
 263def _norm_logcdf(x):
 264    return sc.log_ndtr(x)
 265
 266
 267def _norm_ppf(q):
 268    return sc.ndtri(q)
 269
 270
 271def _norm_sf(x):
 272    return _norm_cdf(-x)
 273
 274
 275def _norm_logsf(x):
 276    return _norm_logcdf(-x)
 277
 278
 279def _norm_isf(q):
 280    return -_norm_ppf(q)
 281
 282
 283class norm_gen(rv_continuous):
 284    r"""A normal continuous random variable.
 285
 286    The location (``loc``) keyword specifies the mean.
 287    The scale (``scale``) keyword specifies the standard deviation.
 288
 289    %(before_notes)s
 290
 291    Notes
 292    -----
 293    The probability density function for `norm` is:
 294
 295    .. math::
 296
 297        f(x) = \frac{\exp(-x^2/2)}{\sqrt{2\pi}}
 298
 299    for a real number :math:`x`.
 300
 301    %(after_notes)s
 302
 303    %(example)s
 304
 305    """
 306    def _rvs(self, size=None, random_state=None):
 307        return random_state.standard_normal(size)
 308
 309    def _pdf(self, x):
 310        # norm.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
 311        return _norm_pdf(x)
 312
 313    def _logpdf(self, x):
 314        return _norm_logpdf(x)
 315
 316    def _cdf(self, x):
 317        return _norm_cdf(x)
 318
 319    def _logcdf(self, x):
 320        return _norm_logcdf(x)
 321
 322    def _sf(self, x):
 323        return _norm_sf(x)
 324
 325    def _logsf(self, x):
 326        return _norm_logsf(x)
 327
 328    def _ppf(self, q):
 329        return _norm_ppf(q)
 330
 331    def _isf(self, q):
 332        return _norm_isf(q)
 333
 334    def _stats(self):
 335        return 0.0, 1.0, 0.0, 0.0
 336
 337    def _entropy(self):
 338        return 0.5*(np.log(2*np.pi)+1)
 339
 340    @_call_super_mom
 341    @replace_notes_in_docstring(rv_continuous, notes="""\
 342        For the normal distribution, method of moments and maximum likelihood
 343        estimation give identical fits, and explicit formulas for the estimates
 344        are available.
 345        This function uses these explicit formulas for the maximum likelihood
 346        estimation of the normal distribution parameters, so the
 347        `optimizer` and `method` arguments are ignored.\n\n""")
 348    def fit(self, data, **kwds):
 349
 350        floc = kwds.pop('floc', None)
 351        fscale = kwds.pop('fscale', None)
 352
 353        _remove_optimizer_parameters(kwds)
 354
 355        if floc is not None and fscale is not None:
 356            # This check is for consistency with `rv_continuous.fit`.
 357            # Without this check, this function would just return the
 358            # parameters that were given.
 359            raise ValueError("All parameters fixed. There is nothing to "
 360                             "optimize.")
 361
 362        data = np.asarray(data)
 363
 364        if not np.isfinite(data).all():
 365            raise RuntimeError("The data contains non-finite values.")
 366
 367        if floc is None:
 368            loc = data.mean()
 369        else:
 370            loc = floc
 371
 372        if fscale is None:
 373            scale = np.sqrt(((data - loc)**2).mean())
 374        else:
 375            scale = fscale
 376
 377        return loc, scale
 378
 379    def _munp(self, n):
 380        """
 381        @returns Moments of standard normal distribution for integer n >= 0
 382
 383        See eq. 16 of https://arxiv.org/abs/1209.4340v2
 384        """
 385        if n % 2 == 0:
 386            return sc.factorial2(n - 1)
 387        else:
 388            return 0.
 389
 390
 391norm = norm_gen(name='norm')
 392
 393
 394class alpha_gen(rv_continuous):
 395    r"""An alpha continuous random variable.
 396
 397    %(before_notes)s
 398
 399    Notes
 400    -----
 401    The probability density function for `alpha` ([1]_, [2]_) is:
 402
 403    .. math::
 404
 405        f(x, a) = \frac{1}{x^2 \Phi(a) \sqrt{2\pi}} *
 406                  \exp(-\frac{1}{2} (a-1/x)^2)
 407
 408    where :math:`\Phi` is the normal CDF, :math:`x > 0`, and :math:`a > 0`.
 409
 410    `alpha` takes ``a`` as a shape parameter.
 411
 412    %(after_notes)s
 413
 414    References
 415    ----------
 416    .. [1] Johnson, Kotz, and Balakrishnan, "Continuous Univariate
 417           Distributions, Volume 1", Second Edition, John Wiley and Sons,
 418           p. 173 (1994).
 419    .. [2] Anthony A. Salvia, "Reliability applications of the Alpha
 420           Distribution", IEEE Transactions on Reliability, Vol. R-34,
 421           No. 3, pp. 251-252 (1985).
 422
 423    %(example)s
 424
 425    """
 426    _support_mask = rv_continuous._open_support_mask
 427
 428    def _pdf(self, x, a):
 429        # alpha.pdf(x, a) = 1/(x**2*Phi(a)*sqrt(2*pi)) * exp(-1/2 * (a-1/x)**2)
 430        return 1.0/(x**2)/_norm_cdf(a)*_norm_pdf(a-1.0/x)
 431
 432    def _logpdf(self, x, a):
 433        return -2*np.log(x) + _norm_logpdf(a-1.0/x) - np.log(_norm_cdf(a))
 434
 435    def _cdf(self, x, a):
 436        return _norm_cdf(a-1.0/x) / _norm_cdf(a)
 437
 438    def _ppf(self, q, a):
 439        return 1.0/np.asarray(a-sc.ndtri(q*_norm_cdf(a)))
 440
 441    def _stats(self, a):
 442        return [np.inf]*2 + [np.nan]*2
 443
 444
 445alpha = alpha_gen(a=0.0, name='alpha')
 446
 447
 448class anglit_gen(rv_continuous):
 449    r"""An anglit continuous random variable.
 450
 451    %(before_notes)s
 452
 453    Notes
 454    -----
 455    The probability density function for `anglit` is:
 456
 457    .. math::
 458
 459        f(x) = \sin(2x + \pi/2) = \cos(2x)
 460
 461    for :math:`-\pi/4 \le x \le \pi/4`.
 462
 463    %(after_notes)s
 464
 465    %(example)s
 466
 467    """
 468    def _pdf(self, x):
 469        # anglit.pdf(x) = sin(2*x + \pi/2) = cos(2*x)
 470        return np.cos(2*x)
 471
 472    def _cdf(self, x):
 473        return np.sin(x+np.pi/4)**2.0
 474
 475    def _ppf(self, q):
 476        return np.arcsin(np.sqrt(q))-np.pi/4
 477
 478    def _stats(self):
 479        return 0.0, np.pi*np.pi/16-0.5, 0.0, -2*(np.pi**4 - 96)/(np.pi*np.pi-8)**2
 480
 481    def _entropy(self):
 482        return 1-np.log(2)
 483
 484
 485anglit = anglit_gen(a=-np.pi/4, b=np.pi/4, name='anglit')
 486
 487
 488class arcsine_gen(rv_continuous):
 489    r"""An arcsine continuous random variable.
 490
 491    %(before_notes)s
 492
 493    Notes
 494    -----
 495    The probability density function for `arcsine` is:
 496
 497    .. math::
 498
 499        f(x) = \frac{1}{\pi \sqrt{x (1-x)}}
 500
 501    for :math:`0 < x < 1`.
 502
 503    %(after_notes)s
 504
 505    %(example)s
 506
 507    """
 508    def _pdf(self, x):
 509        # arcsine.pdf(x) = 1/(pi*sqrt(x*(1-x)))
 510        with np.errstate(divide='ignore'):
 511            return 1.0/np.pi/np.sqrt(x*(1-x))
 512
 513    def _cdf(self, x):
 514        return 2.0/np.pi*np.arcsin(np.sqrt(x))
 515
 516    def _ppf(self, q):
 517        return np.sin(np.pi/2.0*q)**2.0
 518
 519    def _stats(self):
 520        mu = 0.5
 521        mu2 = 1.0/8
 522        g1 = 0
 523        g2 = -3.0/2.0
 524        return mu, mu2, g1, g2
 525
 526    def _entropy(self):
 527        return -0.24156447527049044468
 528
 529
 530arcsine = arcsine_gen(a=0.0, b=1.0, name='arcsine')
 531
 532
 533class FitDataError(ValueError):
 534    # This exception is raised by, for example, beta_gen.fit when both floc
 535    # and fscale are fixed and there are values in the data not in the open
 536    # interval (floc, floc+fscale).
 537    def __init__(self, distr, lower, upper):
 538        self.args = (
 539            "Invalid values in `data`.  Maximum likelihood "
 540            "estimation with {distr!r} requires that {lower!r} < "
 541            "(x - loc)/scale  < {upper!r} for each x in `data`.".format(
 542                distr=distr, lower=lower, upper=upper),
 543        )
 544
 545
 546class FitSolverError(RuntimeError):
 547    # This exception is raised by, for example, beta_gen.fit when
 548    # optimize.fsolve returns with ier != 1.
 549    def __init__(self, mesg):
 550        emsg = "Solver for the MLE equations failed to converge: "
 551        emsg += mesg.replace('\n', '')
 552        self.args = (emsg,)
 553
 554
 555def _beta_mle_a(a, b, n, s1):
 556    # The zeros of this function give the MLE for `a`, with
 557    # `b`, `n` and `s1` given.  `s1` is the sum of the logs of
 558    # the data. `n` is the number of data points.
 559    psiab = sc.psi(a + b)
 560    func = s1 - n * (-psiab + sc.psi(a))
 561    return func
 562
 563
 564def _beta_mle_ab(theta, n, s1, s2):
 565    # Zeros of this function are critical points of
 566    # the maximum likelihood function.  Solving this system
 567    # for theta (which contains a and b) gives the MLE for a and b
 568    # given `n`, `s1` and `s2`.  `s1` is the sum of the logs of the data,
 569    # and `s2` is the sum of the logs of 1 - data.  `n` is the number
 570    # of data points.
 571    a, b = theta
 572    psiab = sc.psi(a + b)
 573    func = [s1 - n * (-psiab + sc.psi(a)),
 574            s2 - n * (-psiab + sc.psi(b))]
 575    return func
 576
 577
 578class beta_gen(rv_continuous):
 579    r"""A beta continuous random variable.
 580
 581    %(before_notes)s
 582
 583    Notes
 584    -----
 585    The probability density function for `beta` is:
 586
 587    .. math::
 588
 589        f(x, a, b) = \frac{\Gamma(a+b) x^{a-1} (1-x)^{b-1}}
 590                          {\Gamma(a) \Gamma(b)}
 591
 592    for :math:`0 <= x <= 1`, :math:`a > 0`, :math:`b > 0`, where
 593    :math:`\Gamma` is the gamma function (`scipy.special.gamma`).
 594
 595    `beta` takes :math:`a` and :math:`b` as shape parameters.
 596
 597    %(after_notes)s
 598
 599    %(example)s
 600
 601    """
 602    def _rvs(self, a, b, size=None, random_state=None):
 603        return random_state.beta(a, b, size)
 604
 605    def _pdf(self, x, a, b):
 606        #                     gamma(a+b) * x**(a-1) * (1-x)**(b-1)
 607        # beta.pdf(x, a, b) = ------------------------------------
 608        #                              gamma(a)*gamma(b)
 609        return _boost._beta_pdf(x, a, b)
 610
 611    def _logpdf(self, x, a, b):
 612        lPx = sc.xlog1py(b - 1.0, -x) + sc.xlogy(a - 1.0, x)
 613        lPx -= sc.betaln(a, b)
 614        return lPx
 615
 616    def _cdf(self, x, a, b):
 617        return _boost._beta_cdf(x, a, b)
 618
 619    def _sf(self, x, a, b):
 620        return _boost._beta_sf(x, a, b)
 621
 622    def _isf(self, x, a, b):
 623        return _boost._beta_isf(x, a, b)
 624
 625    def _ppf(self, q, a, b):
 626        return _boost._beta_ppf(q, a, b)
 627
 628    def _stats(self, a, b):
 629        return(
 630            _boost._beta_mean(a, b),
 631            _boost._beta_variance(a, b),
 632            _boost._beta_skewness(a, b),
 633            _boost._beta_kurtosis_excess(a, b))
 634
 635    def _fitstart(self, data):
 636        g1 = _skew(data)
 637        g2 = _kurtosis(data)
 638
 639        def func(x):
 640            a, b = x
 641            sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)
 642            ku = a**3 - a**2*(2*b-1) + b**2*(b+1) - 2*a*b*(b+2)
 643            ku /= a*b*(a+b+2)*(a+b+3)
 644            ku *= 6
 645            return [sk-g1, ku-g2]
 646        a, b = optimize.fsolve(func, (1.0, 1.0))
 647        return super()._fitstart(data, args=(a, b))
 648
 649    @_call_super_mom
 650    @extend_notes_in_docstring(rv_continuous, notes="""\
 651        In the special case where `method="MLE"` and
 652        both `floc` and `fscale` are given, a
 653        `ValueError` is raised if any value `x` in `data` does not satisfy
 654        `floc < x < floc + fscale`.\n\n""")
 655    def fit(self, data, *args, **kwds):
 656        # Override rv_continuous.fit, so we can more efficiently handle the
 657        # case where floc and fscale are given.
 658
 659        floc = kwds.get('floc', None)
 660        fscale = kwds.get('fscale', None)
 661
 662        if floc is None or fscale is None:
 663            # do general fit
 664            return super().fit(data, *args, **kwds)
 665
 666        # We already got these from kwds, so just pop them.
 667        kwds.pop('floc', None)
 668        kwds.pop('fscale', None)
 669
 670        f0 = _get_fixed_fit_value(kwds, ['f0', 'fa', 'fix_a'])
 671        f1 = _get_fixed_fit_value(kwds, ['f1', 'fb', 'fix_b'])
 672
 673        _remove_optimizer_parameters(kwds)
 674
 675        if f0 is not None and f1 is not None:
 676            # This check is for consistency with `rv_continuous.fit`.
 677            raise ValueError("All parameters fixed. There is nothing to "
 678                             "optimize.")
 679
 680        # Special case: loc and scale are constrained, so we are fitting
 681        # just the shape parameters.  This can be done much more efficiently
 682        # than the method used in `rv_continuous.fit`.  (See the subsection
 683        # "Two unknown parameters" in the section "Maximum likelihood" of
 684        # the Wikipedia article on the Beta distribution for the formulas.)
 685
 686        if not np.isfinite(data).all():
 687            raise RuntimeError("The data contains non-finite values.")
 688
 689        # Normalize the data to the interval [0, 1].
 690        data = (np.ravel(data) - floc) / fscale
 691        if np.any(data <= 0) or np.any(data >= 1):
 692            raise FitDataError("beta", lower=floc, upper=floc + fscale)
 693
 694        xbar = data.mean()
 695
 696        if f0 is not None or f1 is not None:
 697            # One of the shape parameters is fixed.
 698
 699            if f0 is not None:
 700                # The shape parameter a is fixed, so swap the parameters
 701                # and flip the data.  We always solve for `a`.  The result
 702                # will be swapped back before returning.
 703                b = f0
 704                data = 1 - data
 705                xbar = 1 - xbar
 706            else:
 707                b = f1
 708
 709            # Initial guess for a.  Use the formula for the mean of the beta
 710            # distribution, E[x] = a / (a + b), to generate a reasonable
 711            # starting point based on the mean of the data and the given
 712            # value of b.
 713            a = b * xbar / (1 - xbar)
 714
 715            # Compute the MLE for `a` by solving _beta_mle_a.
 716            theta, info, ier, mesg = optimize.fsolve(
 717                _beta_mle_a, a,
 718                args=(b, len(data), np.log(data).sum()),
 719                full_output=True
 720            )
 721            if ier != 1:
 722                raise FitSolverError(mesg=mesg)
 723            a = theta[0]
 724
 725            if f0 is not None:
 726                # The shape parameter a was fixed, so swap back the
 727                # parameters.
 728                a, b = b, a
 729
 730        else:
 731            # Neither of the shape parameters is fixed.
 732
 733            # s1 and s2 are used in the extra arguments passed to _beta_mle_ab
 734            # by optimize.fsolve.
 735            s1 = np.log(data).sum()
 736            s2 = sc.log1p(-data).sum()
 737
 738            # Use the "method of moments" to estimate the initial
 739            # guess for a and b.
 740            fac = xbar * (1 - xbar) / data.var(ddof=0) - 1
 741            a = xbar * fac
 742            b = (1 - xbar) * fac
 743
 744            # Compute the MLE for a and b by solving _beta_mle_ab.
 745            theta, info, ier, mesg = optimize.fsolve(
 746                _beta_mle_ab, [a, b],
 747                args=(len(data), s1, s2),
 748                full_output=True
 749            )
 750            if ier != 1:
 751                raise FitSolverError(mesg=mesg)
 752            a, b = theta
 753
 754        return a, b, floc, fscale
 755
 756
 757beta = beta_gen(a=0.0, b=1.0, name='beta')
 758
 759
 760class betaprime_gen(rv_continuous):
 761    r"""A beta prime continuous random variable.
 762
 763    %(before_notes)s
 764
 765    Notes
 766    -----
 767    The probability density function for `betaprime` is:
 768
 769    .. math::
 770
 771        f(x, a, b) = \frac{x^{a-1} (1+x)^{-a-b}}{\beta(a, b)}
 772
 773    for :math:`x >= 0`, :math:`a > 0`, :math:`b > 0`, where
 774    :math:`\beta(a, b)` is the beta function (see `scipy.special.beta`).
 775
 776    `betaprime` takes ``a`` and ``b`` as shape parameters.
 777
 778    %(after_notes)s
 779
 780    %(example)s
 781
 782    """
 783    _support_mask = rv_continuous._open_support_mask
 784
 785    def _rvs(self, a, b, size=None, random_state=None):
 786        u1 = gamma.rvs(a, size=size, random_state=random_state)
 787        u2 = gamma.rvs(b, size=size, random_state=random_state)
 788        return u1 / u2
 789
 790    def _pdf(self, x, a, b):
 791        # betaprime.pdf(x, a, b) = x**(a-1) * (1+x)**(-a-b) / beta(a, b)
 792        return np.exp(self._logpdf(x, a, b))
 793
 794    def _logpdf(self, x, a, b):
 795        return sc.xlogy(a - 1.0, x) - sc.xlog1py(a + b, x) - sc.betaln(a, b)
 796
 797    def _cdf(self, x, a, b):
 798        return sc.betainc(a, b, x/(1.+x))
 799
 800    def _munp(self, n, a, b):
 801        if n == 1.0:
 802            return np.where(b > 1,
 803                            a/(b-1.0),
 804                            np.inf)
 805        elif n == 2.0:
 806            return np.where(b > 2,
 807                            a*(a+1.0)/((b-2.0)*(b-1.0)),
 808                            np.inf)
 809        elif n == 3.0:
 810            return np.where(b > 3,
 811                            a*(a+1.0)*(a+2.0)/((b-3.0)*(b-2.0)*(b-1.0)),
 812                            np.inf)
 813        elif n == 4.0:
 814            return np.where(b > 4,
 815                            (a*(a + 1.0)*(a + 2.0)*(a + 3.0) /
 816                             ((b - 4.0)*(b - 3.0)*(b - 2.0)*(b - 1.0))),
 817                            np.inf)
 818        else:
 819            raise NotImplementedError
 820
 821
 822betaprime = betaprime_gen(a=0.0, name='betaprime')
 823
 824
 825class bradford_gen(rv_continuous):
 826    r"""A Bradford continuous random variable.
 827
 828    %(before_notes)s
 829
 830    Notes
 831    -----
 832    The probability density function for `bradford` is:
 833
 834    .. math::
 835
 836        f(x, c) = \frac{c}{\log(1+c) (1+cx)}
 837
 838    for :math:`0 <= x <= 1` and :math:`c > 0`.
 839
 840    `bradford` takes ``c`` as a shape parameter for :math:`c`.
 841
 842    %(after_notes)s
 843
 844    %(example)s
 845
 846    """
 847    def _pdf(self, x, c):
 848        # bradford.pdf(x, c) = c / (k * (1+c*x))
 849        return c / (c*x + 1.0) / sc.log1p(c)
 850
 851    def _cdf(self, x, c):
 852        return sc.log1p(c*x) / sc.log1p(c)
 853
 854    def _ppf(self, q, c):
 855        return sc.expm1(q * sc.log1p(c)) / c
 856
 857    def _stats(self, c, moments='mv'):
 858        k = np.log(1.0+c)
 859        mu = (c-k)/(c*k)
 860        mu2 = ((c+2.0)*k-2.0*c)/(2*c*k*k)
 861        g1 = None
 862        g2 = None
 863        if 's' in moments:
 864            g1 = np.sqrt(2)*(12*c*c-9*c*k*(c+2)+2*k*k*(c*(c+3)+3))
 865            g1 /= np.sqrt(c*(c*(k-2)+2*k))*(3*c*(k-2)+6*k)
 866        if 'k' in moments:
 867            g2 = (c**3*(k-3)*(k*(3*k-16)+24)+12*k*c*c*(k-4)*(k-3) +
 868                  6*c*k*k*(3*k-14) + 12*k**3)
 869            g2 /= 3*c*(c*(k-2)+2*k)**2
 870        return mu, mu2, g1, g2
 871
 872    def _entropy(self, c):
 873        k = np.log(1+c)
 874        return k/2.0 - np.log(c/k)
 875
 876
 877bradford = bradford_gen(a=0.0, b=1.0, name='bradford')
 878
 879
 880class burr_gen(rv_continuous):
 881    r"""A Burr (Type III) continuous random variable.
 882
 883    %(before_notes)s
 884
 885    See Also
 886    --------
 887    fisk : a special case of either `burr` or `burr12` with ``d=1``
 888    burr12 : Burr Type XII distribution
 889    mielke : Mielke Beta-Kappa / Dagum distribution
 890
 891    Notes
 892    -----
 893    The probability density function for `burr` is:
 894
 895    .. math::
 896
 897        f(x, c, d) = c d x^{-c - 1} / (1 + x^{-c})^{d + 1}
 898
 899    for :math:`x >= 0` and :math:`c, d > 0`.
 900
 901    `burr` takes :math:`c` and :math:`d` as shape parameters.
 902
 903    This is the PDF corresponding to the third CDF given in Burr's list;
 904    specifically, it is equation (11) in Burr's paper [1]_. The distribution
 905    is also commonly referred to as the Dagum distribution [2]_. If the
 906    parameter :math:`c < 1` then the mean of the distribution does not
 907    exist and if :math:`c < 2` the variance does not exist [2]_.
 908    The PDF is finite at the left endpoint :math:`x = 0` if :math:`c * d >= 1`.
 909
 910    %(after_notes)s
 911
 912    References
 913    ----------
 914    .. [1] Burr, I. W. "Cumulative frequency functions", Annals of
 915       Mathematical Statistics, 13(2), pp 215-232 (1942).
 916    .. [2] https://en.wikipedia.org/wiki/Dagum_distribution
 917    .. [3] Kleiber, Christian. "A guide to the Dagum distributions."
 918       Modeling Income Distributions and Lorenz Curves  pp 97-117 (2008).
 919
 920    %(example)s
 921
 922    """
 923    # Do not set _support_mask to rv_continuous._open_support_mask
 924    # Whether the left-hand endpoint is suitable for pdf evaluation is dependent
 925    # on the values of c and d: if c*d >= 1, the pdf is finite, otherwise infinite.
 926
 927    def _pdf(self, x, c, d):
 928        # burr.pdf(x, c, d) = c * d * x**(-c-1) * (1+x**(-c))**(-d-1)
 929        output = _lazywhere(x == 0, [x, c, d],
 930                   lambda x_, c_, d_: c_ * d_ * (x_**(c_*d_-1)) / (1 + x_**c_),
 931                   f2 = lambda x_, c_, d_: (c_ * d_ * (x_ ** (-c_ - 1.0)) /
 932                                            ((1 + x_ ** (-c_)) ** (d_ + 1.0))))
 933        if output.ndim == 0:
 934            return output[()]
 935        return output
 936
 937    def _logpdf(self, x, c, d):
 938        output = _lazywhere(
 939            x == 0, [x, c, d],
 940            lambda x_, c_, d_: (np.log(c_) + np.log(d_) + sc.xlogy(c_*d_ - 1, x_)
 941                                - (d_+1) * sc.log1p(x_**(c_))),
 942            f2 = lambda x_, c_, d_: (np.log(c_) + np.log(d_)
 943                                     + sc.xlogy(-c_ - 1, x_)
 944                                     - sc.xlog1py(d_+1, x_**(-c_))))
 945        if output.ndim == 0:
 946            return output[()]
 947        return output
 948
 949    def _cdf(self, x, c, d):
 950        return (1 + x**(-c))**(-d)
 951
 952    def _logcdf(self, x, c, d):
 953        return sc.log1p(x**(-c)) * (-d)
 954
 955    def _sf(self, x, c, d):
 956        return np.exp(self._logsf(x, c, d))
 957
 958    def _logsf(self, x, c, d):
 959        return np.log1p(- (1 + x**(-c))**(-d))
 960
 961    def _ppf(self, q, c, d):
 962        return (q**(-1.0/d) - 1)**(-1.0/c)
 963
 964    def _stats(self, c, d):
 965        nc = np.arange(1, 5).reshape(4,1) / c
 966        #ek is the kth raw moment, e1 is the mean e2-e1**2 variance etc.
 967        e1, e2, e3, e4 = sc.beta(d + nc, 1. - nc) * d
 968        mu = np.where(c > 1.0, e1, np.nan)
 969        mu2_if_c = e2 - mu**2
 970        mu2 = np.where(c > 2.0, mu2_if_c, np.nan)
 971        g1 = _lazywhere(
 972            c > 3.0,
 973            (c, e1, e2, e3, mu2_if_c),
 974            lambda c, e1, e2, e3, mu2_if_c: (e3 - 3*e2*e1 + 2*e1**3) / np.sqrt((mu2_if_c)**3),
 975            fillvalue=np.nan)
 976        g2 = _lazywhere(
 977            c > 4.0,
 978            (c, e1, e2, e3, e4, mu2_if_c),
 979            lambda c, e1, e2, e3, e4, mu2_if_c: (
 980                ((e4 - 4*e3*e1 + 6*e2*e1**2 - 3*e1**4) / mu2_if_c**2) - 3),
 981            fillvalue=np.nan)
 982        if np.ndim(c) == 0:
 983            return mu.item(), mu2.item(), g1.item(), g2.item()
 984        return mu, mu2, g1, g2
 985
 986    def _munp(self, n, c, d):
 987        def __munp(n, c, d):
 988            nc = 1. * n / c
 989            return d * sc.beta(1.0 - nc, d + nc)
 990        n, c, d = np.asarray(n), np.asarray(c), np.asarray(d)
 991        return _lazywhere((c > n) & (n == n) & (d == d), (c, d, n),
 992                          lambda c, d, n: __munp(n, c, d),
 993                          np.nan)
 994
 995
 996burr = burr_gen(a=0.0, name='burr')
 997
 998
 999class burr12_gen(rv_continuous):
1000    r"""A Burr (Type XII) continuous random variable.
1001
1002    %(before_notes)s
1003
1004    See Also
1005    --------
1006    fisk : a special case of either `burr` or `burr12` with ``d=1``
1007    burr : Burr Type III distribution
1008
1009    Notes
1010    -----
1011    The probability density function for `burr` is:
1012
1013    .. math::
1014
1015        f(x, c, d) = c d x^{c-1} / (1 + x^c)^{d + 1}
1016
1017    for :math:`x >= 0` and :math:`c, d > 0`.
1018
1019    `burr12` takes ``c`` and ``d`` as shape parameters for :math:`c`
1020    and :math:`d`.
1021
1022    This is the PDF corresponding to the twelfth CDF given in Burr's list;
1023    specifically, it is equation (20) in Burr's paper [1]_.
1024
1025    %(after_notes)s
1026
1027    The Burr type 12 distribution is also sometimes referred to as
1028    the Singh-Maddala distribution from NIST [2]_.
1029
1030    References
1031    ----------
1032    .. [1] Burr, I. W. "Cumulative frequency functions", Annals of
1033       Mathematical Statistics, 13(2), pp 215-232 (1942).
1034
1035    .. [2] https://www.itl.nist.gov/div898/software/dataplot/refman2/auxillar/b12pdf.htm
1036
1037    .. [3] "Burr distribution",
1038       https://en.wikipedia.org/wiki/Burr_distribution
1039
1040    %(example)s
1041
1042    """
1043    def _pdf(self, x, c, d):
1044        # burr12.pdf(x, c, d) = c * d * x**(c-1) * (1+x**(c))**(-d-1)
1045        return np.exp(self._logpdf(x, c, d))
1046
1047    def _logpdf(self, x, c, d):
1048        return np.log(c) + np.log(d) + sc.xlogy(c - 1, x) + sc.xlog1py(-d-1, x**c)
1049
1050    def _cdf(self, x, c, d):
1051        return -sc.expm1(self._logsf(x, c, d))
1052
1053    def _logcdf(self, x, c, d):
1054        return sc.log1p(-(1 + x**c)**(-d))
1055
1056    def _sf(self, x, c, d):
1057        return np.exp(self._logsf(x, c, d))
1058
1059    def _logsf(self, x, c, d):
1060        return sc.xlog1py(-d, x**c)
1061
1062    def _ppf(self, q, c, d):
1063        # The following is an implementation of
1064        #   ((1 - q)**(-1.0/d) - 1)**(1.0/c)
1065        # that does a better job handling small values of q.
1066        return sc.expm1(-1/d * sc.log1p(-q))**(1/c)
1067
1068    def _munp(self, n, c, d):
1069        nc = 1. * n / c
1070        return d * sc.beta(1.0 + nc, d - nc)
1071
1072
1073burr12 = burr12_gen(a=0.0, name='burr12')
1074
1075
1076class fisk_gen(burr_gen):
1077    r"""A Fisk continuous random variable.
1078
1079    The Fisk distribution is also known as the log-logistic distribution.
1080
1081    %(before_notes)s
1082
1083    See Also
1084    --------
1085    burr
1086
1087    Notes
1088    -----
1089    The probability density function for `fisk` is:
1090
1091    .. math::
1092
1093        f(x, c) = c x^{-c-1} (1 + x^{-c})^{-2}
1094
1095    for :math:`x >= 0` and :math:`c > 0`.
1096
1097    `fisk` takes ``c`` as a shape parameter for :math:`c`.
1098
1099    `fisk` is a special case of `burr` or `burr12` with ``d=1``.
1100
1101    %(after_notes)s
1102
1103    %(example)s
1104
1105    """
1106    def _pdf(self, x, c):
1107        # fisk.pdf(x, c) = c * x**(-c-1) * (1 + x**(-c))**(-2)
1108        return burr._pdf(x, c, 1.0)
1109
1110    def _cdf(self, x, c):
1111        return burr._cdf(x, c, 1.0)
1112
1113    def _sf(self, x, c):
1114        return burr._sf(x, c, 1.0)
1115
1116    def _logpdf(self, x, c):
1117        # fisk.pdf(x, c) = c * x**(-c-1) * (1 + x**(-c))**(-2)
1118        return burr._logpdf(x, c, 1.0)
1119
1120    def _logcdf(self, x, c):
1121        return burr._logcdf(x, c, 1.0)
1122
1123    def _logsf(self, x, c):
1124        return burr._logsf(x, c, 1.0)
1125
1126    def _ppf(self, x, c):
1127        return burr._ppf(x, c, 1.0)
1128
1129    def _munp(self, n, c):
1130        return burr._munp(n, c, 1.0)
1131
1132    def _stats(self, c):
1133        return burr._stats(c, 1.0)
1134
1135    def _entropy(self, c):
1136        return 2 - np.log(c)
1137
1138
1139fisk = fisk_gen(a=0.0, name='fisk')
1140
1141
1142class cauchy_gen(rv_continuous):
1143    r"""A Cauchy continuous random variable.
1144
1145    %(before_notes)s
1146
1147    Notes
1148    -----
1149    The probability density function for `cauchy` is
1150
1151    .. math::
1152
1153        f(x) = \frac{1}{\pi (1 + x^2)}
1154
1155    for a real number :math:`x`.
1156
1157    %(after_notes)s
1158
1159    %(example)s
1160
1161    """
1162    def _pdf(self, x):
1163        # cauchy.pdf(x) = 1 / (pi * (1 + x**2))
1164        return 1.0/np.pi/(1.0+x*x)
1165
1166    def _cdf(self, x):
1167        return 0.5 + 1.0/np.pi*np.arctan(x)
1168
1169    def _ppf(self, q):
1170        return np.tan(np.pi*q-np.pi/2.0)
1171
1172    def _sf(self, x):
1173        return 0.5 - 1.0/np.pi*np.arctan(x)
1174
1175    def _isf(self, q):
1176        return np.tan(np.pi/2.0-np.pi*q)
1177
1178    def _stats(self):
1179        return np.nan, np.nan, np.nan, np.nan
1180
1181    def _entropy(self):
1182        return np.log(4*np.pi)
1183
1184    def _fitstart(self, data, args=None):
1185        # Initialize ML guesses using quartiles instead of moments.
1186        p25, p50, p75 = np.percentile(data, [25, 50, 75])
1187        return p50, (p75 - p25)/2
1188
1189
1190cauchy = cauchy_gen(name='cauchy')
1191
1192
1193class chi_gen(rv_continuous):
1194    r"""A chi continuous random variable.
1195
1196    %(before_notes)s
1197
1198    Notes
1199    -----
1200    The probability density function for `chi` is:
1201
1202    .. math::
1203
1204        f(x, k) = \frac{1}{2^{k/2-1} \Gamma \left( k/2 \right)}
1205                   x^{k-1} \exp \left( -x^2/2 \right)
1206
1207    for :math:`x >= 0` and :math:`k > 0` (degrees of freedom, denoted ``df``
1208    in the implementation). :math:`\Gamma` is the gamma function
1209    (`scipy.special.gamma`).
1210
1211    Special cases of `chi` are:
1212
1213        - ``chi(1, loc, scale)`` is equivalent to `halfnorm`
1214        - ``chi(2, 0, scale)`` is equivalent to `rayleigh`
1215        - ``chi(3, 0, scale)`` is equivalent to `maxwell`
1216
1217    `chi` takes ``df`` as a shape parameter.
1218
1219    %(after_notes)s
1220
1221    %(example)s
1222
1223    """
1224
1225    def _rvs(self, df, size=None, random_state=None):
1226        return np.sqrt(chi2.rvs(df, size=size, random_state=random_state))
1227
1228    def _pdf(self, x, df):
1229        #                   x**(df-1) * exp(-x**2/2)
1230        # chi.pdf(x, df) =  -------------------------
1231        #                   2**(df/2-1) * gamma(df/2)
1232        return np.exp(self._logpdf(x, df))
1233
1234    def _logpdf(self, x, df):
1235        l = np.log(2) - .5*np.log(2)*df - sc.gammaln(.5*df)
1236        return l + sc.xlogy(df - 1., x) - .5*x**2
1237
1238    def _cdf(self, x, df):
1239        return sc.gammainc(.5*df, .5*x**2)
1240
1241    def _sf(self, x, df):
1242        return sc.gammaincc(.5*df, .5*x**2)
1243
1244    def _ppf(self, q, df):
1245        return np.sqrt(2*sc.gammaincinv(.5*df, q))
1246
1247    def _isf(self, q, df):
1248        return np.sqrt(2*sc.gammainccinv(.5*df, q))
1249
1250    def _stats(self, df):
1251        mu = np.sqrt(2)*sc.gamma(df/2.0+0.5)/sc.gamma(df/2.0)
1252        mu2 = df - mu*mu
1253        g1 = (2*mu**3.0 + mu*(1-2*df))/np.asarray(np.power(mu2, 1.5))
1254        g2 = 2*df*(1.0-df)-6*mu**4 + 4*mu**2 * (2*df-1)
1255        g2 /= np.asarray(mu2**2.0)
1256        return mu, mu2, g1, g2
1257
1258
1259chi = chi_gen(a=0.0, name='chi')
1260
1261
1262class chi2_gen(rv_continuous):
1263    r"""A chi-squared continuous random variable.
1264
1265    For the noncentral chi-square distribution, see `ncx2`.
1266
1267    %(before_notes)s
1268
1269    See Also
1270    --------
1271    ncx2
1272
1273    Notes
1274    -----
1275    The probability density function for `chi2` is:
1276
1277    .. math::
1278
1279        f(x, k) = \frac{1}{2^{k/2} \Gamma \left( k/2 \right)}
1280                   x^{k/2-1} \exp \left( -x/2 \right)
1281
1282    for :math:`x > 0`  and :math:`k > 0` (degrees of freedom, denoted ``df``
1283    in the implementation).
1284
1285    `chi2` takes ``df`` as a shape parameter.
1286
1287    The chi-squared distribution is a special case of the gamma
1288    distribution, with gamma parameters ``a = df/2``, ``loc = 0`` and
1289    ``scale = 2``.
1290
1291    %(after_notes)s
1292
1293    %(example)s
1294
1295    """
1296    def _rvs(self, df, size=None, random_state=None):
1297        return random_state.chisquare(df, size)
1298
1299    def _pdf(self, x, df):
1300        # chi2.pdf(x, df) = 1 / (2*gamma(df/2)) * (x/2)**(df/2-1) * exp(-x/2)
1301        return np.exp(self._logpdf(x, df))
1302
1303    def _logpdf(self, x, df):
1304        return sc.xlogy(df/2.-1, x) - x/2. - sc.gammaln(df/2.) - (np.log(2)*df)/2.
1305
1306    def _cdf(self, x, df):
1307        return sc.chdtr(df, x)
1308
1309    def _sf(self, x, df):
1310        return sc.chdtrc(df, x)
1311
1312    def _isf(self, p, df):
1313        return sc.chdtri(df, p)
1314
1315    def _ppf(self, p, df):
1316        return 2*sc.gammaincinv(df/2, p)
1317
1318    def _stats(self, df):
1319        mu = df
1320        mu2 = 2*df
1321        g1 = 2*np.sqrt(2.0/df)
1322        g2 = 12.0/df
1323        return mu, mu2, g1, g2
1324
1325
1326chi2 = chi2_gen(a=0.0, name='chi2')
1327
1328
1329class cosine_gen(rv_continuous):
1330    r"""A cosine continuous random variable.
1331
1332    %(before_notes)s
1333
1334    Notes
1335    -----
1336    The cosine distribution is an approximation to the normal distribution.
1337    The probability density function for `cosine` is:
1338
1339    .. math::
1340
1341        f(x) = \frac{1}{2\pi} (1+\cos(x))
1342
1343    for :math:`-\pi \le x \le \pi`.
1344
1345    %(after_notes)s
1346
1347    %(example)s
1348
1349    """
1350    def _pdf(self, x):
1351        # cosine.pdf(x) = 1/(2*pi) * (1+cos(x))
1352        return 1.0/2/np.pi*(1+np.cos(x))
1353
1354    def _cdf(self, x):
1355        return scu._cosine_cdf(x)
1356
1357    def _sf(self, x):
1358        return scu._cosine_cdf(-x)
1359
1360    def _ppf(self, p):
1361        return scu._cosine_invcdf(p)
1362
1363    def _isf(self, p):
1364        return -scu._cosine_invcdf(p)
1365
1366    def _stats(self):
1367        return 0.0, np.pi*np.pi/3.0-2.0, 0.0, -6.0*(np.pi**4-90)/(5.0*(np.pi*np.pi-6)**2)
1368
1369    def _entropy(self):
1370        return np.log(4*np.pi)-1.0
1371
1372
1373cosine = cosine_gen(a=-np.pi, b=np.pi, name='cosine')
1374
1375
1376class dgamma_gen(rv_continuous):
1377    r"""A double gamma continuous random variable.
1378
1379    %(before_notes)s
1380
1381    Notes
1382    -----
1383    The probability density function for `dgamma` is:
1384
1385    .. math::
1386
1387        f(x, a) = \frac{1}{2\Gamma(a)} |x|^{a-1} \exp(-|x|)
1388
1389    for a real number :math:`x` and :math:`a > 0`. :math:`\Gamma` is the
1390    gamma function (`scipy.special.gamma`).
1391
1392    `dgamma` takes ``a`` as a shape parameter for :math:`a`.
1393
1394    %(after_notes)s
1395
1396    %(example)s
1397
1398    """
1399    def _rvs(self, a, size=None, random_state=None):
1400        u = random_state.uniform(size=size)
1401        gm = gamma.rvs(a, size=size, random_state=random_state)
1402        return gm * np.where(u >= 0.5, 1, -1)
1403
1404    def _pdf(self, x, a):
1405        # dgamma.pdf(x, a) = 1 / (2*gamma(a)) * abs(x)**(a-1) * exp(-abs(x))
1406        ax = abs(x)
1407        return 1.0/(2*sc.gamma(a))*ax**(a-1.0) * np.exp(-ax)
1408
1409    def _logpdf(self, x, a):
1410        ax = abs(x)
1411        return sc.xlogy(a - 1.0, ax) - ax - np.log(2) - sc.gammaln(a)
1412
1413    def _cdf(self, x, a):
1414        fac = 0.5*sc.gammainc(a, abs(x))
1415        return np.where(x > 0, 0.5 + fac, 0.5 - fac)
1416
1417    def _sf(self, x, a):
1418        fac = 0.5*sc.gammainc(a, abs(x))
1419        return np.where(x > 0, 0.5-fac, 0.5+fac)
1420
1421    def _ppf(self, q, a):
1422        fac = sc.gammainccinv(a, 1-abs(2*q-1))
1423        return np.where(q > 0.5, fac, -fac)
1424
1425    def _stats(self, a):
1426        mu2 = a*(a+1.0)
1427        return 0.0, mu2, 0.0, (a+2.0)*(a+3.0)/mu2-3.0
1428
1429
1430dgamma = dgamma_gen(name='dgamma')
1431
1432
1433class dweibull_gen(rv_continuous):
1434    r"""A double Weibull continuous random variable.
1435
1436    %(before_notes)s
1437
1438    Notes
1439    -----
1440    The probability density function for `dweibull` is given by
1441
1442    .. math::
1443
1444        f(x, c) = c / 2 |x|^{c-1} \exp(-|x|^c)
1445
1446    for a real number :math:`x` and :math:`c > 0`.
1447
1448    `dweibull` takes ``c`` as a shape parameter for :math:`c`.
1449
1450    %(after_notes)s
1451
1452    %(example)s
1453
1454    """
1455    def _rvs(self, c, size=None, random_state=None):
1456        u = random_state.uniform(size=size)
1457        w = weibull_min.rvs(c, size=size, random_state=random_state)
1458        return w * (np.where(u >= 0.5, 1, -1))
1459
1460    def _pdf(self, x, c):
1461        # dweibull.pdf(x, c) = c / 2 * abs(x)**(c-1) * exp(-abs(x)**c)
1462        ax = abs(x)
1463        Px = c / 2.0 * ax**(c-1.0) * np.exp(-ax**c)
1464        return Px
1465
1466    def _logpdf(self, x, c):
1467        ax = abs(x)
1468        return np.log(c) - np.log(2.0) + sc.xlogy(c - 1.0, ax) - ax**c
1469
1470    def _cdf(self, x, c):
1471        Cx1 = 0.5 * np.exp(-abs(x)**c)
1472        return np.where(x > 0, 1 - Cx1, Cx1)
1473
1474    def _ppf(self, q, c):
1475        fac = 2. * np.where(q <= 0.5, q, 1. - q)
1476        fac = np.power(-np.log(fac), 1.0 / c)
1477        return np.where(q > 0.5, fac, -fac)
1478
1479    def _munp(self, n, c):
1480        return (1 - (n % 2)) * sc.gamma(1.0 + 1.0 * n / c)
1481
1482    # since we know that all odd moments are zeros, return them at once.
1483    # returning Nones from _stats makes the public stats call _munp
1484    # so overall we're saving one or two gamma function evaluations here.
1485    def _stats(self, c):
1486        return 0, None, 0, None
1487
1488
1489dweibull = dweibull_gen(name='dweibull')
1490
1491
1492class expon_gen(rv_continuous):
1493    r"""An exponential continuous random variable.
1494
1495    %(before_notes)s
1496
1497    Notes
1498    -----
1499    The probability density function for `expon` is:
1500
1501    .. math::
1502
1503        f(x) = \exp(-x)
1504
1505    for :math:`x \ge 0`.
1506
1507    %(after_notes)s
1508
1509    A common parameterization for `expon` is in terms of the rate parameter
1510    ``lambda``, such that ``pdf = lambda * exp(-lambda * x)``. This
1511    parameterization corresponds to using ``scale = 1 / lambda``.
1512
1513    The exponential distribution is a special case of the gamma
1514    distributions, with gamma shape parameter ``a = 1``.
1515
1516    %(example)s
1517
1518    """
1519    def _rvs(self, size=None, random_state=None):
1520        return random_state.standard_exponential(size)
1521
1522    def _pdf(self, x):
1523        # expon.pdf(x) = exp(-x)
1524        return np.exp(-x)
1525
1526    def _logpdf(self, x):
1527        return -x
1528
1529    def _cdf(self, x):
1530        return -sc.expm1(-x)
1531
1532    def _ppf(self, q):
1533        return -sc.log1p(-q)
1534
1535    def _sf(self, x):
1536        return np.exp(-x)
1537
1538    def _logsf(self, x):
1539        return -x
1540
1541    def _isf(self, q):
1542        return -np.log(q)
1543
1544    def _stats(self):
1545        return 1.0, 1.0, 2.0, 6.0
1546
1547    def _entropy(self):
1548        return 1.0
1549
1550    @_call_super_mom
1551    @replace_notes_in_docstring(rv_continuous, notes="""\
1552        When `method='MLE'`,
1553        this function uses explicit formulas for the maximum likelihood
1554        estimation of the exponential distribution parameters, so the
1555        `optimizer`, `loc` and `scale` keyword arguments are
1556        ignored.\n\n""")
1557    def fit(self, data, *args, **kwds):
1558        if len(args) > 0:
1559            raise TypeError("Too many arguments.")
1560
1561        floc = kwds.pop('floc', None)
1562        fscale = kwds.pop('fscale', None)
1563
1564        _remove_optimizer_parameters(kwds)
1565
1566        if floc is not None and fscale is not None:
1567            # This check is for consistency with `rv_continuous.fit`.
1568            raise ValueError("All parameters fixed. There is nothing to "
1569                             "optimize.")
1570
1571        data = np.asarray(data)
1572
1573        if not np.isfinite(data).all():
1574            raise RuntimeError("The data contains non-finite values.")
1575
1576        data_min = data.min()
1577
1578        if floc is None:
1579            # ML estimate of the location is the minimum of the data.
1580            loc = data_min
1581        else:
1582            loc = floc
1583            if data_min < loc:
1584                # There are values that are less than the specified loc.
1585                raise FitDataError("expon", lower=floc, upper=np.inf)
1586
1587        if fscale is None:
1588            # ML estimate of the scale is the shifted mean.
1589            scale = data.mean() - loc
1590        else:
1591            scale = fscale
1592
1593        # We expect the return values to be floating point, so ensure it
1594        # by explicitly converting to float.
1595        return float(loc), float(scale)
1596
1597
1598expon = expon_gen(a=0.0, name='expon')
1599
1600
1601class exponnorm_gen(rv_continuous):
1602    r"""An exponentially modified Normal continuous random variable.
1603
1604    Also known as the exponentially modified Gaussian distribution [1]_.
1605
1606    %(before_notes)s
1607
1608    Notes
1609    -----
1610    The probability density function for `exponnorm` is:
1611
1612    .. math::
1613
1614        f(x, K) = \frac{1}{2K} \exp\left(\frac{1}{2 K^2} - x / K \right)
1615                  \text{erfc}\left(-\frac{x - 1/K}{\sqrt{2}}\right)
1616
1617    where :math:`x` is a real number and :math:`K > 0`.
1618
1619    It can be thought of as the sum of a standard normal random variable
1620    and an independent exponentially distributed random variable with rate
1621    ``1/K``.
1622
1623    %(after_notes)s
1624
1625    An alternative parameterization of this distribution (for example, in
1626    the Wikpedia article [1]_) involves three parameters, :math:`\mu`,
1627    :math:`\lambda` and :math:`\sigma`.
1628
1629    In the present parameterization this corresponds to having ``loc`` and
1630    ``scale`` equal to :math:`\mu` and :math:`\sigma`, respectively, and
1631    shape parameter :math:`K = 1/(\sigma\lambda)`.
1632
1633    .. versionadded:: 0.16.0
1634
1635    References
1636    ----------
1637    .. [1] Exponentially modified Gaussian distribution, Wikipedia,
1638           https://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution
1639
1640    %(example)s
1641
1642    """
1643    def _rvs(self, K, size=None, random_state=None):
1644        expval = random_state.standard_exponential(size) * K
1645        gval = random_state.standard_normal(size)
1646        return expval + gval
1647
1648    def _pdf(self, x, K):
1649        return np.exp(self._logpdf(x, K))
1650
1651    def _logpdf(self, x, K):
1652        invK = 1.0 / K
1653        exparg = invK * (0.5 * invK - x)
1654        return exparg + _norm_logcdf(x - invK) - np.log(K)
1655
1656    def _cdf(self, x, K):
1657        invK = 1.0 / K
1658        expval = invK * (0.5 * invK - x)
1659        logprod = expval + _norm_logcdf(x - invK)
1660        return _norm_cdf(x) - np.exp(logprod)
1661
1662    def _sf(self, x, K):
1663        invK = 1.0 / K
1664        expval = invK * (0.5 * invK - x)
1665        logprod = expval + _norm_logcdf(x - invK)
1666        return _norm_cdf(-x) + np.exp(logprod)
1667
1668    def _stats(self, K):
1669        K2 = K * K
1670        opK2 = 1.0 + K2
1671        skw = 2 * K**3 * opK2**(-1.5)
1672        krt = 6.0 * K2 * K2 * opK2**(-2)
1673        return K, opK2, skw, krt
1674
1675
1676exponnorm = exponnorm_gen(name='exponnorm')
1677
1678
1679class exponweib_gen(rv_continuous):
1680    r"""An exponentiated Weibull continuous random variable.
1681
1682    %(before_notes)s
1683
1684    See Also
1685    --------
1686    weibull_min, numpy.random.Generator.weibull
1687
1688    Notes
1689    -----
1690    The probability density function for `exponweib` is:
1691
1692    .. math::
1693
1694        f(x, a, c) = a c [1-\exp(-x^c)]^{a-1} \exp(-x^c) x^{c-1}
1695
1696    and its cumulative distribution function is:
1697
1698    .. math::
1699
1700        F(x, a, c) = [1-\exp(-x^c)]^a
1701
1702    for :math:`x > 0`, :math:`a > 0`, :math:`c > 0`.
1703
1704    `exponweib` takes :math:`a` and :math:`c` as shape parameters:
1705
1706    * :math:`a` is the exponentiation parameter,
1707      with the special case :math:`a=1` corresponding to the
1708      (non-exponentiated) Weibull distribution `weibull_min`.
1709    * :math:`c` is the shape parameter of the non-exponentiated Weibull law.
1710
1711    %(after_notes)s
1712
1713    References
1714    ----------
1715    https://en.wikipedia.org/wiki/Exponentiated_Weibull_distribution
1716
1717    %(example)s
1718
1719    """
1720    def _pdf(self, x, a, c):
1721        # exponweib.pdf(x, a, c) =
1722        #     a * c * (1-exp(-x**c))**(a-1) * exp(-x**c)*x**(c-1)
1723        return np.exp(self._logpdf(x, a, c))
1724
1725    def _logpdf(self, x, a, c):
1726        negxc = -x**c
1727        exm1c = -sc.expm1(negxc)
1728        logp = (np.log(a) + np.log(c) + sc.xlogy(a - 1.0, exm1c) +
1729                negxc + sc.xlogy(c - 1.0, x))
1730        return logp
1731
1732    def _cdf(self, x, a, c):
1733        exm1c = -sc.expm1(-x**c)
1734        return exm1c**a
1735
1736    def _ppf(self, q, a, c):
1737        return (-sc.log1p(-q**(1.0/a)))**np.asarray(1.0/c)
1738
1739
1740exponweib = exponweib_gen(a=0.0, name='exponweib')
1741
1742
1743class exponpow_gen(rv_continuous):
1744    r"""An exponential power continuous random variable.
1745
1746    %(before_notes)s
1747
1748    Notes
1749    -----
1750    The probability density function for `exponpow` is:
1751
1752    .. math::
1753
1754        f(x, b) = b x^{b-1} \exp(1 + x^b - \exp(x^b))
1755
1756    for :math:`x \ge 0`, :math:`b > 0`.  Note that this is a different
1757    distribution from the exponential power distribution that is also known
1758    under the names "generalized normal" or "generalized Gaussian".
1759
1760    `exponpow` takes ``b`` as a shape parameter for :math:`b`.
1761
1762    %(after_notes)s
1763
1764    References
1765    ----------
1766    http://www.math.wm.edu/~leemis/chart/UDR/PDFs/Exponentialpower.pdf
1767
1768    %(example)s
1769
1770    """
1771    def _pdf(self, x, b):
1772        # exponpow.pdf(x, b) = b * x**(b-1) * exp(1 + x**b - exp(x**b))
1773        return np.exp(self._logpdf(x, b))
1774
1775    def _logpdf(self, x, b):
1776        xb = x**b
1777        f = 1 + np.log(b) + sc.xlogy(b - 1.0, x) + xb - np.exp(xb)
1778        return f
1779
1780    def _cdf(self, x, b):
1781        return -sc.expm1(-sc.expm1(x**b))
1782
1783    def _sf(self, x, b):
1784        return np.exp(-sc.expm1(x**b))
1785
1786    def _isf(self, x, b):
1787        return (sc.log1p(-np.log(x)))**(1./b)
1788
1789    def _ppf(self, q, b):
1790        return pow(sc.log1p(-sc.log1p(-q)), 1.0/b)
1791
1792
1793exponpow = exponpow_gen(a=0.0, name='exponpow')
1794
1795
1796class fatiguelife_gen(rv_continuous):
1797    r"""A fatigue-life (Birnbaum-Saunders) continuous random variable.
1798
1799    %(before_notes)s
1800
1801    Notes
1802    -----
1803    The probability density function for `fatiguelife` is:
1804
1805    .. math::
1806
1807        f(x, c) = \frac{x+1}{2c\sqrt{2\pi x^3}} \exp(-\frac{(x-1)^2}{2x c^2})
1808
1809    for :math:`x >= 0` and :math:`c > 0`.
1810
1811    `fatiguelife` takes ``c`` as a shape parameter for :math:`c`.
1812
1813    %(after_notes)s
1814
1815    References
1816    ----------
1817    .. [1] "Birnbaum-Saunders distribution",
1818           https://en.wikipedia.org/wiki/Birnbaum-Saunders_distribution
1819
1820    %(example)s
1821
1822    """
1823    _support_mask = rv_continuous._open_support_mask
1824
1825    def _rvs(self, c, size=None, random_state=None):
1826        z = random_state.standard_normal(size)
1827        x = 0.5*c*z
1828        x2 = x*x
1829        t = 1.0 + 2*x2 + 2*x*np.sqrt(1 + x2)
1830        return t
1831
1832    def _pdf(self, x, c):
1833        # fatiguelife.pdf(x, c) =
1834        #     (x+1) / (2*c*sqrt(2*pi*x**3)) * exp(-(x-1)**2/(2*x*c**2))
1835        return np.exp(self._logpdf(x, c))
1836
1837    def _logpdf(self, x, c):
1838        return (np.log(x+1) - (x-1)**2 / (2.0*x*c**2) - np.log(2*c) -
1839                0.5*(np.log(2*np.pi) + 3*np.log(x)))
1840
1841    def _cdf(self, x, c):
1842        return _norm_cdf(1.0 / c * (np.sqrt(x) - 1.0/np.sqrt(x)))
1843
1844    def _ppf(self, q, c):
1845        tmp = c*sc.ndtri(q)
1846        return 0.25 * (tmp + np.sqrt(tmp**2 + 4))**