/scipy/stats/_continuous_distns.py
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))**