#### /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))**