Coverage for /usr/lib/python3/dist-packages/scipy/stats/_discrete_distns.py: 33%
640 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1#
2# Author: Travis Oliphant 2002-2011 with contributions from
3# SciPy Developers 2004-2011
4#
5from functools import partial
7from scipy import special
8from scipy.special import entr, logsumexp, betaln, gammaln as gamln, zeta
9from scipy._lib._util import _lazywhere, rng_integers
10from scipy.interpolate import interp1d
12from numpy import floor, ceil, log, exp, sqrt, log1p, expm1, tanh, cosh, sinh
14import numpy as np
16from ._distn_infrastructure import (rv_discrete, get_distribution_names,
17 _check_shape, _ShapeInfo)
18import scipy.stats._boost as _boost
19from ._biasedurn import (_PyFishersNCHypergeometric,
20 _PyWalleniusNCHypergeometric,
21 _PyStochasticLib3)
24def _isintegral(x):
25 return x == np.round(x)
28class binom_gen(rv_discrete):
29 r"""A binomial discrete random variable.
31 %(before_notes)s
33 Notes
34 -----
35 The probability mass function for `binom` is:
37 .. math::
39 f(k) = \binom{n}{k} p^k (1-p)^{n-k}
41 for :math:`k \in \{0, 1, \dots, n\}`, :math:`0 \leq p \leq 1`
43 `binom` takes :math:`n` and :math:`p` as shape parameters,
44 where :math:`p` is the probability of a single success
45 and :math:`1-p` is the probability of a single failure.
47 %(after_notes)s
49 %(example)s
51 See Also
52 --------
53 hypergeom, nbinom, nhypergeom
55 """
56 def _shape_info(self):
57 return [_ShapeInfo("n", True, (0, np.inf), (True, False)),
58 _ShapeInfo("p", False, (0, 1), (True, True))]
60 def _rvs(self, n, p, size=None, random_state=None):
61 return random_state.binomial(n, p, size)
63 def _argcheck(self, n, p):
64 return (n >= 0) & _isintegral(n) & (p >= 0) & (p <= 1)
66 def _get_support(self, n, p):
67 return self.a, n
69 def _logpmf(self, x, n, p):
70 k = floor(x)
71 combiln = (gamln(n+1) - (gamln(k+1) + gamln(n-k+1)))
72 return combiln + special.xlogy(k, p) + special.xlog1py(n-k, -p)
74 def _pmf(self, x, n, p):
75 # binom.pmf(k) = choose(n, k) * p**k * (1-p)**(n-k)
76 return _boost._binom_pdf(x, n, p)
78 def _cdf(self, x, n, p):
79 k = floor(x)
80 return _boost._binom_cdf(k, n, p)
82 def _sf(self, x, n, p):
83 k = floor(x)
84 return _boost._binom_sf(k, n, p)
86 def _isf(self, x, n, p):
87 return _boost._binom_isf(x, n, p)
89 def _ppf(self, q, n, p):
90 return _boost._binom_ppf(q, n, p)
92 def _stats(self, n, p, moments='mv'):
93 mu = _boost._binom_mean(n, p)
94 var = _boost._binom_variance(n, p)
95 g1, g2 = None, None
96 if 's' in moments:
97 g1 = _boost._binom_skewness(n, p)
98 if 'k' in moments:
99 g2 = _boost._binom_kurtosis_excess(n, p)
100 return mu, var, g1, g2
102 def _entropy(self, n, p):
103 k = np.r_[0:n + 1]
104 vals = self._pmf(k, n, p)
105 return np.sum(entr(vals), axis=0)
108binom = binom_gen(name='binom')
111class bernoulli_gen(binom_gen):
112 r"""A Bernoulli discrete random variable.
114 %(before_notes)s
116 Notes
117 -----
118 The probability mass function for `bernoulli` is:
120 .. math::
122 f(k) = \begin{cases}1-p &\text{if } k = 0\\
123 p &\text{if } k = 1\end{cases}
125 for :math:`k` in :math:`\{0, 1\}`, :math:`0 \leq p \leq 1`
127 `bernoulli` takes :math:`p` as shape parameter,
128 where :math:`p` is the probability of a single success
129 and :math:`1-p` is the probability of a single failure.
131 %(after_notes)s
133 %(example)s
135 """
136 def _shape_info(self):
137 return [_ShapeInfo("p", False, (0, 1), (True, True))]
139 def _rvs(self, p, size=None, random_state=None):
140 return binom_gen._rvs(self, 1, p, size=size, random_state=random_state)
142 def _argcheck(self, p):
143 return (p >= 0) & (p <= 1)
145 def _get_support(self, p):
146 # Overrides binom_gen._get_support!x
147 return self.a, self.b
149 def _logpmf(self, x, p):
150 return binom._logpmf(x, 1, p)
152 def _pmf(self, x, p):
153 # bernoulli.pmf(k) = 1-p if k = 0
154 # = p if k = 1
155 return binom._pmf(x, 1, p)
157 def _cdf(self, x, p):
158 return binom._cdf(x, 1, p)
160 def _sf(self, x, p):
161 return binom._sf(x, 1, p)
163 def _isf(self, x, p):
164 return binom._isf(x, 1, p)
166 def _ppf(self, q, p):
167 return binom._ppf(q, 1, p)
169 def _stats(self, p):
170 return binom._stats(1, p)
172 def _entropy(self, p):
173 return entr(p) + entr(1-p)
176bernoulli = bernoulli_gen(b=1, name='bernoulli')
179class betabinom_gen(rv_discrete):
180 r"""A beta-binomial discrete random variable.
182 %(before_notes)s
184 Notes
185 -----
186 The beta-binomial distribution is a binomial distribution with a
187 probability of success `p` that follows a beta distribution.
189 The probability mass function for `betabinom` is:
191 .. math::
193 f(k) = \binom{n}{k} \frac{B(k + a, n - k + b)}{B(a, b)}
195 for :math:`k \in \{0, 1, \dots, n\}`, :math:`n \geq 0`, :math:`a > 0`,
196 :math:`b > 0`, where :math:`B(a, b)` is the beta function.
198 `betabinom` takes :math:`n`, :math:`a`, and :math:`b` as shape parameters.
200 References
201 ----------
202 .. [1] https://en.wikipedia.org/wiki/Beta-binomial_distribution
204 %(after_notes)s
206 .. versionadded:: 1.4.0
208 See Also
209 --------
210 beta, binom
212 %(example)s
214 """
215 def _shape_info(self):
216 return [_ShapeInfo("n", True, (0, np.inf), (True, False)),
217 _ShapeInfo("a", False, (0, np.inf), (False, False)),
218 _ShapeInfo("b", False, (0, np.inf), (False, False))]
220 def _rvs(self, n, a, b, size=None, random_state=None):
221 p = random_state.beta(a, b, size)
222 return random_state.binomial(n, p, size)
224 def _get_support(self, n, a, b):
225 return 0, n
227 def _argcheck(self, n, a, b):
228 return (n >= 0) & _isintegral(n) & (a > 0) & (b > 0)
230 def _logpmf(self, x, n, a, b):
231 k = floor(x)
232 combiln = -log(n + 1) - betaln(n - k + 1, k + 1)
233 return combiln + betaln(k + a, n - k + b) - betaln(a, b)
235 def _pmf(self, x, n, a, b):
236 return exp(self._logpmf(x, n, a, b))
238 def _stats(self, n, a, b, moments='mv'):
239 e_p = a / (a + b)
240 e_q = 1 - e_p
241 mu = n * e_p
242 var = n * (a + b + n) * e_p * e_q / (a + b + 1)
243 g1, g2 = None, None
244 if 's' in moments:
245 g1 = 1.0 / sqrt(var)
246 g1 *= (a + b + 2 * n) * (b - a)
247 g1 /= (a + b + 2) * (a + b)
248 if 'k' in moments:
249 g2 = (a + b).astype(e_p.dtype)
250 g2 *= (a + b - 1 + 6 * n)
251 g2 += 3 * a * b * (n - 2)
252 g2 += 6 * n ** 2
253 g2 -= 3 * e_p * b * n * (6 - n)
254 g2 -= 18 * e_p * e_q * n ** 2
255 g2 *= (a + b) ** 2 * (1 + a + b)
256 g2 /= (n * a * b * (a + b + 2) * (a + b + 3) * (a + b + n))
257 g2 -= 3
258 return mu, var, g1, g2
261betabinom = betabinom_gen(name='betabinom')
264class nbinom_gen(rv_discrete):
265 r"""A negative binomial discrete random variable.
267 %(before_notes)s
269 Notes
270 -----
271 Negative binomial distribution describes a sequence of i.i.d. Bernoulli
272 trials, repeated until a predefined, non-random number of successes occurs.
274 The probability mass function of the number of failures for `nbinom` is:
276 .. math::
278 f(k) = \binom{k+n-1}{n-1} p^n (1-p)^k
280 for :math:`k \ge 0`, :math:`0 < p \leq 1`
282 `nbinom` takes :math:`n` and :math:`p` as shape parameters where :math:`n`
283 is the number of successes, :math:`p` is the probability of a single
284 success, and :math:`1-p` is the probability of a single failure.
286 Another common parameterization of the negative binomial distribution is
287 in terms of the mean number of failures :math:`\mu` to achieve :math:`n`
288 successes. The mean :math:`\mu` is related to the probability of success
289 as
291 .. math::
293 p = \frac{n}{n + \mu}
295 The number of successes :math:`n` may also be specified in terms of a
296 "dispersion", "heterogeneity", or "aggregation" parameter :math:`\alpha`,
297 which relates the mean :math:`\mu` to the variance :math:`\sigma^2`,
298 e.g. :math:`\sigma^2 = \mu + \alpha \mu^2`. Regardless of the convention
299 used for :math:`\alpha`,
301 .. math::
303 p &= \frac{\mu}{\sigma^2} \\
304 n &= \frac{\mu^2}{\sigma^2 - \mu}
306 %(after_notes)s
308 %(example)s
310 See Also
311 --------
312 hypergeom, binom, nhypergeom
314 """
315 def _shape_info(self):
316 return [_ShapeInfo("n", True, (0, np.inf), (True, False)),
317 _ShapeInfo("p", False, (0, 1), (True, True))]
319 def _rvs(self, n, p, size=None, random_state=None):
320 return random_state.negative_binomial(n, p, size)
322 def _argcheck(self, n, p):
323 return (n > 0) & (p > 0) & (p <= 1)
325 def _pmf(self, x, n, p):
326 # nbinom.pmf(k) = choose(k+n-1, n-1) * p**n * (1-p)**k
327 return _boost._nbinom_pdf(x, n, p)
329 def _logpmf(self, x, n, p):
330 coeff = gamln(n+x) - gamln(x+1) - gamln(n)
331 return coeff + n*log(p) + special.xlog1py(x, -p)
333 def _cdf(self, x, n, p):
334 k = floor(x)
335 return _boost._nbinom_cdf(k, n, p)
337 def _logcdf(self, x, n, p):
338 k = floor(x)
339 cdf = self._cdf(k, n, p)
340 cond = cdf > 0.5
342 def f1(k, n, p):
343 return np.log1p(-special.betainc(k + 1, n, 1 - p))
345 # do calc in place
346 logcdf = cdf
347 with np.errstate(divide='ignore'):
348 logcdf[cond] = f1(k[cond], n[cond], p[cond])
349 logcdf[~cond] = np.log(cdf[~cond])
350 return logcdf
352 def _sf(self, x, n, p):
353 k = floor(x)
354 return _boost._nbinom_sf(k, n, p)
356 def _isf(self, x, n, p):
357 with np.errstate(over='ignore'): # see gh-17432
358 return _boost._nbinom_isf(x, n, p)
360 def _ppf(self, q, n, p):
361 with np.errstate(over='ignore'): # see gh-17432
362 return _boost._nbinom_ppf(q, n, p)
364 def _stats(self, n, p):
365 return (
366 _boost._nbinom_mean(n, p),
367 _boost._nbinom_variance(n, p),
368 _boost._nbinom_skewness(n, p),
369 _boost._nbinom_kurtosis_excess(n, p),
370 )
373nbinom = nbinom_gen(name='nbinom')
376class geom_gen(rv_discrete):
377 r"""A geometric discrete random variable.
379 %(before_notes)s
381 Notes
382 -----
383 The probability mass function for `geom` is:
385 .. math::
387 f(k) = (1-p)^{k-1} p
389 for :math:`k \ge 1`, :math:`0 < p \leq 1`
391 `geom` takes :math:`p` as shape parameter,
392 where :math:`p` is the probability of a single success
393 and :math:`1-p` is the probability of a single failure.
395 %(after_notes)s
397 See Also
398 --------
399 planck
401 %(example)s
403 """
405 def _shape_info(self):
406 return [_ShapeInfo("p", False, (0, 1), (True, True))]
408 def _rvs(self, p, size=None, random_state=None):
409 return random_state.geometric(p, size=size)
411 def _argcheck(self, p):
412 return (p <= 1) & (p > 0)
414 def _pmf(self, k, p):
415 return np.power(1-p, k-1) * p
417 def _logpmf(self, k, p):
418 return special.xlog1py(k - 1, -p) + log(p)
420 def _cdf(self, x, p):
421 k = floor(x)
422 return -expm1(log1p(-p)*k)
424 def _sf(self, x, p):
425 return np.exp(self._logsf(x, p))
427 def _logsf(self, x, p):
428 k = floor(x)
429 return k*log1p(-p)
431 def _ppf(self, q, p):
432 vals = ceil(log1p(-q) / log1p(-p))
433 temp = self._cdf(vals-1, p)
434 return np.where((temp >= q) & (vals > 0), vals-1, vals)
436 def _stats(self, p):
437 mu = 1.0/p
438 qr = 1.0-p
439 var = qr / p / p
440 g1 = (2.0-p) / sqrt(qr)
441 g2 = np.polyval([1, -6, 6], p)/(1.0-p)
442 return mu, var, g1, g2
444 def _entropy(self, p):
445 return -np.log(p) - np.log1p(-p) * (1.0-p) / p
448geom = geom_gen(a=1, name='geom', longname="A geometric")
451class hypergeom_gen(rv_discrete):
452 r"""A hypergeometric discrete random variable.
454 The hypergeometric distribution models drawing objects from a bin.
455 `M` is the total number of objects, `n` is total number of Type I objects.
456 The random variate represents the number of Type I objects in `N` drawn
457 without replacement from the total population.
459 %(before_notes)s
461 Notes
462 -----
463 The symbols used to denote the shape parameters (`M`, `n`, and `N`) are not
464 universally accepted. See the Examples for a clarification of the
465 definitions used here.
467 The probability mass function is defined as,
469 .. math:: p(k, M, n, N) = \frac{\binom{n}{k} \binom{M - n}{N - k}}
470 {\binom{M}{N}}
472 for :math:`k \in [\max(0, N - M + n), \min(n, N)]`, where the binomial
473 coefficients are defined as,
475 .. math:: \binom{n}{k} \equiv \frac{n!}{k! (n - k)!}.
477 %(after_notes)s
479 Examples
480 --------
481 >>> import numpy as np
482 >>> from scipy.stats import hypergeom
483 >>> import matplotlib.pyplot as plt
485 Suppose we have a collection of 20 animals, of which 7 are dogs. Then if
486 we want to know the probability of finding a given number of dogs if we
487 choose at random 12 of the 20 animals, we can initialize a frozen
488 distribution and plot the probability mass function:
490 >>> [M, n, N] = [20, 7, 12]
491 >>> rv = hypergeom(M, n, N)
492 >>> x = np.arange(0, n+1)
493 >>> pmf_dogs = rv.pmf(x)
495 >>> fig = plt.figure()
496 >>> ax = fig.add_subplot(111)
497 >>> ax.plot(x, pmf_dogs, 'bo')
498 >>> ax.vlines(x, 0, pmf_dogs, lw=2)
499 >>> ax.set_xlabel('# of dogs in our group of chosen animals')
500 >>> ax.set_ylabel('hypergeom PMF')
501 >>> plt.show()
503 Instead of using a frozen distribution we can also use `hypergeom`
504 methods directly. To for example obtain the cumulative distribution
505 function, use:
507 >>> prb = hypergeom.cdf(x, M, n, N)
509 And to generate random numbers:
511 >>> R = hypergeom.rvs(M, n, N, size=10)
513 See Also
514 --------
515 nhypergeom, binom, nbinom
517 """
518 def _shape_info(self):
519 return [_ShapeInfo("M", True, (0, np.inf), (True, False)),
520 _ShapeInfo("n", True, (0, np.inf), (True, False)),
521 _ShapeInfo("N", True, (0, np.inf), (True, False))]
523 def _rvs(self, M, n, N, size=None, random_state=None):
524 return random_state.hypergeometric(n, M-n, N, size=size)
526 def _get_support(self, M, n, N):
527 return np.maximum(N-(M-n), 0), np.minimum(n, N)
529 def _argcheck(self, M, n, N):
530 cond = (M > 0) & (n >= 0) & (N >= 0)
531 cond &= (n <= M) & (N <= M)
532 cond &= _isintegral(M) & _isintegral(n) & _isintegral(N)
533 return cond
535 def _logpmf(self, k, M, n, N):
536 tot, good = M, n
537 bad = tot - good
538 result = (betaln(good+1, 1) + betaln(bad+1, 1) + betaln(tot-N+1, N+1) -
539 betaln(k+1, good-k+1) - betaln(N-k+1, bad-N+k+1) -
540 betaln(tot+1, 1))
541 return result
543 def _pmf(self, k, M, n, N):
544 return _boost._hypergeom_pdf(k, n, N, M)
546 def _cdf(self, k, M, n, N):
547 return _boost._hypergeom_cdf(k, n, N, M)
549 def _stats(self, M, n, N):
550 M, n, N = 1. * M, 1. * n, 1. * N
551 m = M - n
553 # Boost kurtosis_excess doesn't return the same as the value
554 # computed here.
555 g2 = M * (M + 1) - 6. * N * (M - N) - 6. * n * m
556 g2 *= (M - 1) * M * M
557 g2 += 6. * n * N * (M - N) * m * (5. * M - 6)
558 g2 /= n * N * (M - N) * m * (M - 2.) * (M - 3.)
559 return (
560 _boost._hypergeom_mean(n, N, M),
561 _boost._hypergeom_variance(n, N, M),
562 _boost._hypergeom_skewness(n, N, M),
563 g2,
564 )
566 def _entropy(self, M, n, N):
567 k = np.r_[N - (M - n):min(n, N) + 1]
568 vals = self.pmf(k, M, n, N)
569 return np.sum(entr(vals), axis=0)
571 def _sf(self, k, M, n, N):
572 return _boost._hypergeom_sf(k, n, N, M)
574 def _logsf(self, k, M, n, N):
575 res = []
576 for quant, tot, good, draw in zip(*np.broadcast_arrays(k, M, n, N)):
577 if (quant + 0.5) * (tot + 0.5) < (good - 0.5) * (draw - 0.5):
578 # Less terms to sum if we calculate log(1-cdf)
579 res.append(log1p(-exp(self.logcdf(quant, tot, good, draw))))
580 else:
581 # Integration over probability mass function using logsumexp
582 k2 = np.arange(quant + 1, draw + 1)
583 res.append(logsumexp(self._logpmf(k2, tot, good, draw)))
584 return np.asarray(res)
586 def _logcdf(self, k, M, n, N):
587 res = []
588 for quant, tot, good, draw in zip(*np.broadcast_arrays(k, M, n, N)):
589 if (quant + 0.5) * (tot + 0.5) > (good - 0.5) * (draw - 0.5):
590 # Less terms to sum if we calculate log(1-sf)
591 res.append(log1p(-exp(self.logsf(quant, tot, good, draw))))
592 else:
593 # Integration over probability mass function using logsumexp
594 k2 = np.arange(0, quant + 1)
595 res.append(logsumexp(self._logpmf(k2, tot, good, draw)))
596 return np.asarray(res)
599hypergeom = hypergeom_gen(name='hypergeom')
602class nhypergeom_gen(rv_discrete):
603 r"""A negative hypergeometric discrete random variable.
605 Consider a box containing :math:`M` balls:, :math:`n` red and
606 :math:`M-n` blue. We randomly sample balls from the box, one
607 at a time and *without* replacement, until we have picked :math:`r`
608 blue balls. `nhypergeom` is the distribution of the number of
609 red balls :math:`k` we have picked.
611 %(before_notes)s
613 Notes
614 -----
615 The symbols used to denote the shape parameters (`M`, `n`, and `r`) are not
616 universally accepted. See the Examples for a clarification of the
617 definitions used here.
619 The probability mass function is defined as,
621 .. math:: f(k; M, n, r) = \frac{{{k+r-1}\choose{k}}{{M-r-k}\choose{n-k}}}
622 {{M \choose n}}
624 for :math:`k \in [0, n]`, :math:`n \in [0, M]`, :math:`r \in [0, M-n]`,
625 and the binomial coefficient is:
627 .. math:: \binom{n}{k} \equiv \frac{n!}{k! (n - k)!}.
629 It is equivalent to observing :math:`k` successes in :math:`k+r-1`
630 samples with :math:`k+r`'th sample being a failure. The former
631 can be modelled as a hypergeometric distribution. The probability
632 of the latter is simply the number of failures remaining
633 :math:`M-n-(r-1)` divided by the size of the remaining population
634 :math:`M-(k+r-1)`. This relationship can be shown as:
636 .. math:: NHG(k;M,n,r) = HG(k;M,n,k+r-1)\frac{(M-n-(r-1))}{(M-(k+r-1))}
638 where :math:`NHG` is probability mass function (PMF) of the
639 negative hypergeometric distribution and :math:`HG` is the
640 PMF of the hypergeometric distribution.
642 %(after_notes)s
644 Examples
645 --------
646 >>> import numpy as np
647 >>> from scipy.stats import nhypergeom
648 >>> import matplotlib.pyplot as plt
650 Suppose we have a collection of 20 animals, of which 7 are dogs.
651 Then if we want to know the probability of finding a given number
652 of dogs (successes) in a sample with exactly 12 animals that
653 aren't dogs (failures), we can initialize a frozen distribution
654 and plot the probability mass function:
656 >>> M, n, r = [20, 7, 12]
657 >>> rv = nhypergeom(M, n, r)
658 >>> x = np.arange(0, n+2)
659 >>> pmf_dogs = rv.pmf(x)
661 >>> fig = plt.figure()
662 >>> ax = fig.add_subplot(111)
663 >>> ax.plot(x, pmf_dogs, 'bo')
664 >>> ax.vlines(x, 0, pmf_dogs, lw=2)
665 >>> ax.set_xlabel('# of dogs in our group with given 12 failures')
666 >>> ax.set_ylabel('nhypergeom PMF')
667 >>> plt.show()
669 Instead of using a frozen distribution we can also use `nhypergeom`
670 methods directly. To for example obtain the probability mass
671 function, use:
673 >>> prb = nhypergeom.pmf(x, M, n, r)
675 And to generate random numbers:
677 >>> R = nhypergeom.rvs(M, n, r, size=10)
679 To verify the relationship between `hypergeom` and `nhypergeom`, use:
681 >>> from scipy.stats import hypergeom, nhypergeom
682 >>> M, n, r = 45, 13, 8
683 >>> k = 6
684 >>> nhypergeom.pmf(k, M, n, r)
685 0.06180776620271643
686 >>> hypergeom.pmf(k, M, n, k+r-1) * (M - n - (r-1)) / (M - (k+r-1))
687 0.06180776620271644
689 See Also
690 --------
691 hypergeom, binom, nbinom
693 References
694 ----------
695 .. [1] Negative Hypergeometric Distribution on Wikipedia
696 https://en.wikipedia.org/wiki/Negative_hypergeometric_distribution
698 .. [2] Negative Hypergeometric Distribution from
699 http://www.math.wm.edu/~leemis/chart/UDR/PDFs/Negativehypergeometric.pdf
701 """
703 def _shape_info(self):
704 return [_ShapeInfo("M", True, (0, np.inf), (True, False)),
705 _ShapeInfo("n", True, (0, np.inf), (True, False)),
706 _ShapeInfo("r", True, (0, np.inf), (True, False))]
708 def _get_support(self, M, n, r):
709 return 0, n
711 def _argcheck(self, M, n, r):
712 cond = (n >= 0) & (n <= M) & (r >= 0) & (r <= M-n)
713 cond &= _isintegral(M) & _isintegral(n) & _isintegral(r)
714 return cond
716 def _rvs(self, M, n, r, size=None, random_state=None):
718 @_vectorize_rvs_over_shapes
719 def _rvs1(M, n, r, size, random_state):
720 # invert cdf by calculating all values in support, scalar M, n, r
721 a, b = self.support(M, n, r)
722 ks = np.arange(a, b+1)
723 cdf = self.cdf(ks, M, n, r)
724 ppf = interp1d(cdf, ks, kind='next', fill_value='extrapolate')
725 rvs = ppf(random_state.uniform(size=size)).astype(int)
726 if size is None:
727 return rvs.item()
728 return rvs
730 return _rvs1(M, n, r, size=size, random_state=random_state)
732 def _logpmf(self, k, M, n, r):
733 cond = ((r == 0) & (k == 0))
734 result = _lazywhere(~cond, (k, M, n, r),
735 lambda k, M, n, r:
736 (-betaln(k+1, r) + betaln(k+r, 1) -
737 betaln(n-k+1, M-r-n+1) + betaln(M-r-k+1, 1) +
738 betaln(n+1, M-n+1) - betaln(M+1, 1)),
739 fillvalue=0.0)
740 return result
742 def _pmf(self, k, M, n, r):
743 # same as the following but numerically more precise
744 # return comb(k+r-1, k) * comb(M-r-k, n-k) / comb(M, n)
745 return exp(self._logpmf(k, M, n, r))
747 def _stats(self, M, n, r):
748 # Promote the datatype to at least float
749 # mu = rn / (M-n+1)
750 M, n, r = 1.*M, 1.*n, 1.*r
751 mu = r*n / (M-n+1)
753 var = r*(M+1)*n / ((M-n+1)*(M-n+2)) * (1 - r / (M-n+1))
755 # The skew and kurtosis are mathematically
756 # intractable so return `None`. See [2]_.
757 g1, g2 = None, None
758 return mu, var, g1, g2
761nhypergeom = nhypergeom_gen(name='nhypergeom')
764# FIXME: Fails _cdfvec
765class logser_gen(rv_discrete):
766 r"""A Logarithmic (Log-Series, Series) discrete random variable.
768 %(before_notes)s
770 Notes
771 -----
772 The probability mass function for `logser` is:
774 .. math::
776 f(k) = - \frac{p^k}{k \log(1-p)}
778 for :math:`k \ge 1`, :math:`0 < p < 1`
780 `logser` takes :math:`p` as shape parameter,
781 where :math:`p` is the probability of a single success
782 and :math:`1-p` is the probability of a single failure.
784 %(after_notes)s
786 %(example)s
788 """
790 def _shape_info(self):
791 return [_ShapeInfo("p", False, (0, 1), (True, True))]
793 def _rvs(self, p, size=None, random_state=None):
794 # looks wrong for p>0.5, too few k=1
795 # trying to use generic is worse, no k=1 at all
796 return random_state.logseries(p, size=size)
798 def _argcheck(self, p):
799 return (p > 0) & (p < 1)
801 def _pmf(self, k, p):
802 # logser.pmf(k) = - p**k / (k*log(1-p))
803 return -np.power(p, k) * 1.0 / k / special.log1p(-p)
805 def _stats(self, p):
806 r = special.log1p(-p)
807 mu = p / (p - 1.0) / r
808 mu2p = -p / r / (p - 1.0)**2
809 var = mu2p - mu*mu
810 mu3p = -p / r * (1.0+p) / (1.0 - p)**3
811 mu3 = mu3p - 3*mu*mu2p + 2*mu**3
812 g1 = mu3 / np.power(var, 1.5)
814 mu4p = -p / r * (
815 1.0 / (p-1)**2 - 6*p / (p - 1)**3 + 6*p*p / (p-1)**4)
816 mu4 = mu4p - 4*mu3p*mu + 6*mu2p*mu*mu - 3*mu**4
817 g2 = mu4 / var**2 - 3.0
818 return mu, var, g1, g2
821logser = logser_gen(a=1, name='logser', longname='A logarithmic')
824class poisson_gen(rv_discrete):
825 r"""A Poisson discrete random variable.
827 %(before_notes)s
829 Notes
830 -----
831 The probability mass function for `poisson` is:
833 .. math::
835 f(k) = \exp(-\mu) \frac{\mu^k}{k!}
837 for :math:`k \ge 0`.
839 `poisson` takes :math:`\mu \geq 0` as shape parameter.
840 When :math:`\mu = 0`, the ``pmf`` method
841 returns ``1.0`` at quantile :math:`k = 0`.
843 %(after_notes)s
845 %(example)s
847 """
849 def _shape_info(self):
850 return [_ShapeInfo("mu", False, (0, np.inf), (True, False))]
852 # Override rv_discrete._argcheck to allow mu=0.
853 def _argcheck(self, mu):
854 return mu >= 0
856 def _rvs(self, mu, size=None, random_state=None):
857 return random_state.poisson(mu, size)
859 def _logpmf(self, k, mu):
860 Pk = special.xlogy(k, mu) - gamln(k + 1) - mu
861 return Pk
863 def _pmf(self, k, mu):
864 # poisson.pmf(k) = exp(-mu) * mu**k / k!
865 return exp(self._logpmf(k, mu))
867 def _cdf(self, x, mu):
868 k = floor(x)
869 return special.pdtr(k, mu)
871 def _sf(self, x, mu):
872 k = floor(x)
873 return special.pdtrc(k, mu)
875 def _ppf(self, q, mu):
876 vals = ceil(special.pdtrik(q, mu))
877 vals1 = np.maximum(vals - 1, 0)
878 temp = special.pdtr(vals1, mu)
879 return np.where(temp >= q, vals1, vals)
881 def _stats(self, mu):
882 var = mu
883 tmp = np.asarray(mu)
884 mu_nonzero = tmp > 0
885 g1 = _lazywhere(mu_nonzero, (tmp,), lambda x: sqrt(1.0/x), np.inf)
886 g2 = _lazywhere(mu_nonzero, (tmp,), lambda x: 1.0/x, np.inf)
887 return mu, var, g1, g2
890poisson = poisson_gen(name="poisson", longname='A Poisson')
893class planck_gen(rv_discrete):
894 r"""A Planck discrete exponential random variable.
896 %(before_notes)s
898 Notes
899 -----
900 The probability mass function for `planck` is:
902 .. math::
904 f(k) = (1-\exp(-\lambda)) \exp(-\lambda k)
906 for :math:`k \ge 0` and :math:`\lambda > 0`.
908 `planck` takes :math:`\lambda` as shape parameter. The Planck distribution
909 can be written as a geometric distribution (`geom`) with
910 :math:`p = 1 - \exp(-\lambda)` shifted by ``loc = -1``.
912 %(after_notes)s
914 See Also
915 --------
916 geom
918 %(example)s
920 """
921 def _shape_info(self):
922 return [_ShapeInfo("lambda", False, (0, np.inf), (False, False))]
924 def _argcheck(self, lambda_):
925 return lambda_ > 0
927 def _pmf(self, k, lambda_):
928 return -expm1(-lambda_)*exp(-lambda_*k)
930 def _cdf(self, x, lambda_):
931 k = floor(x)
932 return -expm1(-lambda_*(k+1))
934 def _sf(self, x, lambda_):
935 return exp(self._logsf(x, lambda_))
937 def _logsf(self, x, lambda_):
938 k = floor(x)
939 return -lambda_*(k+1)
941 def _ppf(self, q, lambda_):
942 vals = ceil(-1.0/lambda_ * log1p(-q)-1)
943 vals1 = (vals-1).clip(*(self._get_support(lambda_)))
944 temp = self._cdf(vals1, lambda_)
945 return np.where(temp >= q, vals1, vals)
947 def _rvs(self, lambda_, size=None, random_state=None):
948 # use relation to geometric distribution for sampling
949 p = -expm1(-lambda_)
950 return random_state.geometric(p, size=size) - 1.0
952 def _stats(self, lambda_):
953 mu = 1/expm1(lambda_)
954 var = exp(-lambda_)/(expm1(-lambda_))**2
955 g1 = 2*cosh(lambda_/2.0)
956 g2 = 4+2*cosh(lambda_)
957 return mu, var, g1, g2
959 def _entropy(self, lambda_):
960 C = -expm1(-lambda_)
961 return lambda_*exp(-lambda_)/C - log(C)
964planck = planck_gen(a=0, name='planck', longname='A discrete exponential ')
967class boltzmann_gen(rv_discrete):
968 r"""A Boltzmann (Truncated Discrete Exponential) random variable.
970 %(before_notes)s
972 Notes
973 -----
974 The probability mass function for `boltzmann` is:
976 .. math::
978 f(k) = (1-\exp(-\lambda)) \exp(-\lambda k) / (1-\exp(-\lambda N))
980 for :math:`k = 0,..., N-1`.
982 `boltzmann` takes :math:`\lambda > 0` and :math:`N > 0` as shape parameters.
984 %(after_notes)s
986 %(example)s
988 """
989 def _shape_info(self):
990 return [_ShapeInfo("lambda_", False, (0, np.inf), (False, False)),
991 _ShapeInfo("N", True, (0, np.inf), (False, False))]
993 def _argcheck(self, lambda_, N):
994 return (lambda_ > 0) & (N > 0) & _isintegral(N)
996 def _get_support(self, lambda_, N):
997 return self.a, N - 1
999 def _pmf(self, k, lambda_, N):
1000 # boltzmann.pmf(k) =
1001 # (1-exp(-lambda_)*exp(-lambda_*k)/(1-exp(-lambda_*N))
1002 fact = (1-exp(-lambda_))/(1-exp(-lambda_*N))
1003 return fact*exp(-lambda_*k)
1005 def _cdf(self, x, lambda_, N):
1006 k = floor(x)
1007 return (1-exp(-lambda_*(k+1)))/(1-exp(-lambda_*N))
1009 def _ppf(self, q, lambda_, N):
1010 qnew = q*(1-exp(-lambda_*N))
1011 vals = ceil(-1.0/lambda_ * log(1-qnew)-1)
1012 vals1 = (vals-1).clip(0.0, np.inf)
1013 temp = self._cdf(vals1, lambda_, N)
1014 return np.where(temp >= q, vals1, vals)
1016 def _stats(self, lambda_, N):
1017 z = exp(-lambda_)
1018 zN = exp(-lambda_*N)
1019 mu = z/(1.0-z)-N*zN/(1-zN)
1020 var = z/(1.0-z)**2 - N*N*zN/(1-zN)**2
1021 trm = (1-zN)/(1-z)
1022 trm2 = (z*trm**2 - N*N*zN)
1023 g1 = z*(1+z)*trm**3 - N**3*zN*(1+zN)
1024 g1 = g1 / trm2**(1.5)
1025 g2 = z*(1+4*z+z*z)*trm**4 - N**4 * zN*(1+4*zN+zN*zN)
1026 g2 = g2 / trm2 / trm2
1027 return mu, var, g1, g2
1030boltzmann = boltzmann_gen(name='boltzmann', a=0,
1031 longname='A truncated discrete exponential ')
1034class randint_gen(rv_discrete):
1035 r"""A uniform discrete random variable.
1037 %(before_notes)s
1039 Notes
1040 -----
1041 The probability mass function for `randint` is:
1043 .. math::
1045 f(k) = \frac{1}{\texttt{high} - \texttt{low}}
1047 for :math:`k \in \{\texttt{low}, \dots, \texttt{high} - 1\}`.
1049 `randint` takes :math:`\texttt{low}` and :math:`\texttt{high}` as shape
1050 parameters.
1052 %(after_notes)s
1054 %(example)s
1056 """
1058 def _shape_info(self):
1059 return [_ShapeInfo("low", True, (-np.inf, np.inf), (False, False)),
1060 _ShapeInfo("high", True, (-np.inf, np.inf), (False, False))]
1062 def _argcheck(self, low, high):
1063 return (high > low) & _isintegral(low) & _isintegral(high)
1065 def _get_support(self, low, high):
1066 return low, high-1
1068 def _pmf(self, k, low, high):
1069 # randint.pmf(k) = 1./(high - low)
1070 p = np.ones_like(k) / (high - low)
1071 return np.where((k >= low) & (k < high), p, 0.)
1073 def _cdf(self, x, low, high):
1074 k = floor(x)
1075 return (k - low + 1.) / (high - low)
1077 def _ppf(self, q, low, high):
1078 vals = ceil(q * (high - low) + low) - 1
1079 vals1 = (vals - 1).clip(low, high)
1080 temp = self._cdf(vals1, low, high)
1081 return np.where(temp >= q, vals1, vals)
1083 def _stats(self, low, high):
1084 m2, m1 = np.asarray(high), np.asarray(low)
1085 mu = (m2 + m1 - 1.0) / 2
1086 d = m2 - m1
1087 var = (d*d - 1) / 12.0
1088 g1 = 0.0
1089 g2 = -6.0/5.0 * (d*d + 1.0) / (d*d - 1.0)
1090 return mu, var, g1, g2
1092 def _rvs(self, low, high, size=None, random_state=None):
1093 """An array of *size* random integers >= ``low`` and < ``high``."""
1094 if np.asarray(low).size == 1 and np.asarray(high).size == 1:
1095 # no need to vectorize in that case
1096 return rng_integers(random_state, low, high, size=size)
1098 if size is not None:
1099 # NumPy's RandomState.randint() doesn't broadcast its arguments.
1100 # Use `broadcast_to()` to extend the shapes of low and high
1101 # up to size. Then we can use the numpy.vectorize'd
1102 # randint without needing to pass it a `size` argument.
1103 low = np.broadcast_to(low, size)
1104 high = np.broadcast_to(high, size)
1105 randint = np.vectorize(partial(rng_integers, random_state),
1106 otypes=[np.int_])
1107 return randint(low, high)
1109 def _entropy(self, low, high):
1110 return log(high - low)
1113randint = randint_gen(name='randint', longname='A discrete uniform '
1114 '(random integer)')
1117# FIXME: problems sampling.
1118class zipf_gen(rv_discrete):
1119 r"""A Zipf (Zeta) discrete random variable.
1121 %(before_notes)s
1123 See Also
1124 --------
1125 zipfian
1127 Notes
1128 -----
1129 The probability mass function for `zipf` is:
1131 .. math::
1133 f(k, a) = \frac{1}{\zeta(a) k^a}
1135 for :math:`k \ge 1`, :math:`a > 1`.
1137 `zipf` takes :math:`a > 1` as shape parameter. :math:`\zeta` is the
1138 Riemann zeta function (`scipy.special.zeta`)
1140 The Zipf distribution is also known as the zeta distribution, which is
1141 a special case of the Zipfian distribution (`zipfian`).
1143 %(after_notes)s
1145 References
1146 ----------
1147 .. [1] "Zeta Distribution", Wikipedia,
1148 https://en.wikipedia.org/wiki/Zeta_distribution
1150 %(example)s
1152 Confirm that `zipf` is the large `n` limit of `zipfian`.
1154 >>> import numpy as np
1155 >>> from scipy.stats import zipfian
1156 >>> k = np.arange(11)
1157 >>> np.allclose(zipf.pmf(k, a), zipfian.pmf(k, a, n=10000000))
1158 True
1160 """
1162 def _shape_info(self):
1163 return [_ShapeInfo("a", False, (1, np.inf), (False, False))]
1165 def _rvs(self, a, size=None, random_state=None):
1166 return random_state.zipf(a, size=size)
1168 def _argcheck(self, a):
1169 return a > 1
1171 def _pmf(self, k, a):
1172 # zipf.pmf(k, a) = 1/(zeta(a) * k**a)
1173 Pk = 1.0 / special.zeta(a, 1) / k**a
1174 return Pk
1176 def _munp(self, n, a):
1177 return _lazywhere(
1178 a > n + 1, (a, n),
1179 lambda a, n: special.zeta(a - n, 1) / special.zeta(a, 1),
1180 np.inf)
1183zipf = zipf_gen(a=1, name='zipf', longname='A Zipf')
1186def _gen_harmonic_gt1(n, a):
1187 """Generalized harmonic number, a > 1"""
1188 # See https://en.wikipedia.org/wiki/Harmonic_number; search for "hurwitz"
1189 return zeta(a, 1) - zeta(a, n+1)
1192def _gen_harmonic_leq1(n, a):
1193 """Generalized harmonic number, a <= 1"""
1194 if not np.size(n):
1195 return n
1196 n_max = np.max(n) # loop starts at maximum of all n
1197 out = np.zeros_like(a, dtype=float)
1198 # add terms of harmonic series; starting from smallest to avoid roundoff
1199 for i in np.arange(n_max, 0, -1, dtype=float):
1200 mask = i <= n # don't add terms after nth
1201 out[mask] += 1/i**a[mask]
1202 return out
1205def _gen_harmonic(n, a):
1206 """Generalized harmonic number"""
1207 n, a = np.broadcast_arrays(n, a)
1208 return _lazywhere(a > 1, (n, a),
1209 f=_gen_harmonic_gt1, f2=_gen_harmonic_leq1)
1212class zipfian_gen(rv_discrete):
1213 r"""A Zipfian discrete random variable.
1215 %(before_notes)s
1217 See Also
1218 --------
1219 zipf
1221 Notes
1222 -----
1223 The probability mass function for `zipfian` is:
1225 .. math::
1227 f(k, a, n) = \frac{1}{H_{n,a} k^a}
1229 for :math:`k \in \{1, 2, \dots, n-1, n\}`, :math:`a \ge 0`,
1230 :math:`n \in \{1, 2, 3, \dots\}`.
1232 `zipfian` takes :math:`a` and :math:`n` as shape parameters.
1233 :math:`H_{n,a}` is the :math:`n`:sup:`th` generalized harmonic
1234 number of order :math:`a`.
1236 The Zipfian distribution reduces to the Zipf (zeta) distribution as
1237 :math:`n \rightarrow \infty`.
1239 %(after_notes)s
1241 References
1242 ----------
1243 .. [1] "Zipf's Law", Wikipedia, https://en.wikipedia.org/wiki/Zipf's_law
1244 .. [2] Larry Leemis, "Zipf Distribution", Univariate Distribution
1245 Relationships. http://www.math.wm.edu/~leemis/chart/UDR/PDFs/Zipf.pdf
1247 %(example)s
1249 Confirm that `zipfian` reduces to `zipf` for large `n`, `a > 1`.
1251 >>> import numpy as np
1252 >>> from scipy.stats import zipf
1253 >>> k = np.arange(11)
1254 >>> np.allclose(zipfian.pmf(k, a=3.5, n=10000000), zipf.pmf(k, a=3.5))
1255 True
1257 """
1259 def _shape_info(self):
1260 return [_ShapeInfo("a", False, (0, np.inf), (True, False)),
1261 _ShapeInfo("n", True, (0, np.inf), (False, False))]
1263 def _argcheck(self, a, n):
1264 # we need np.asarray here because moment (maybe others) don't convert
1265 return (a >= 0) & (n > 0) & (n == np.asarray(n, dtype=int))
1267 def _get_support(self, a, n):
1268 return 1, n
1270 def _pmf(self, k, a, n):
1271 return 1.0 / _gen_harmonic(n, a) / k**a
1273 def _cdf(self, k, a, n):
1274 return _gen_harmonic(k, a) / _gen_harmonic(n, a)
1276 def _sf(self, k, a, n):
1277 k = k + 1 # # to match SciPy convention
1278 # see http://www.math.wm.edu/~leemis/chart/UDR/PDFs/Zipf.pdf
1279 return ((k**a*(_gen_harmonic(n, a) - _gen_harmonic(k, a)) + 1)
1280 / (k**a*_gen_harmonic(n, a)))
1282 def _stats(self, a, n):
1283 # see # see http://www.math.wm.edu/~leemis/chart/UDR/PDFs/Zipf.pdf
1284 Hna = _gen_harmonic(n, a)
1285 Hna1 = _gen_harmonic(n, a-1)
1286 Hna2 = _gen_harmonic(n, a-2)
1287 Hna3 = _gen_harmonic(n, a-3)
1288 Hna4 = _gen_harmonic(n, a-4)
1289 mu1 = Hna1/Hna
1290 mu2n = (Hna2*Hna - Hna1**2)
1291 mu2d = Hna**2
1292 mu2 = mu2n / mu2d
1293 g1 = (Hna3/Hna - 3*Hna1*Hna2/Hna**2 + 2*Hna1**3/Hna**3)/mu2**(3/2)
1294 g2 = (Hna**3*Hna4 - 4*Hna**2*Hna1*Hna3 + 6*Hna*Hna1**2*Hna2
1295 - 3*Hna1**4) / mu2n**2
1296 g2 -= 3
1297 return mu1, mu2, g1, g2
1300zipfian = zipfian_gen(a=1, name='zipfian', longname='A Zipfian')
1303class dlaplace_gen(rv_discrete):
1304 r"""A Laplacian discrete random variable.
1306 %(before_notes)s
1308 Notes
1309 -----
1310 The probability mass function for `dlaplace` is:
1312 .. math::
1314 f(k) = \tanh(a/2) \exp(-a |k|)
1316 for integers :math:`k` and :math:`a > 0`.
1318 `dlaplace` takes :math:`a` as shape parameter.
1320 %(after_notes)s
1322 %(example)s
1324 """
1326 def _shape_info(self):
1327 return [_ShapeInfo("a", False, (0, np.inf), (False, False))]
1329 def _pmf(self, k, a):
1330 # dlaplace.pmf(k) = tanh(a/2) * exp(-a*abs(k))
1331 return tanh(a/2.0) * exp(-a * abs(k))
1333 def _cdf(self, x, a):
1334 k = floor(x)
1336 def f(k, a):
1337 return 1.0 - exp(-a * k) / (exp(a) + 1)
1339 def f2(k, a):
1340 return exp(a * (k + 1)) / (exp(a) + 1)
1342 return _lazywhere(k >= 0, (k, a), f=f, f2=f2)
1344 def _ppf(self, q, a):
1345 const = 1 + exp(a)
1346 vals = ceil(np.where(q < 1.0 / (1 + exp(-a)),
1347 log(q*const) / a - 1,
1348 -log((1-q) * const) / a))
1349 vals1 = vals - 1
1350 return np.where(self._cdf(vals1, a) >= q, vals1, vals)
1352 def _stats(self, a):
1353 ea = exp(a)
1354 mu2 = 2.*ea/(ea-1.)**2
1355 mu4 = 2.*ea*(ea**2+10.*ea+1.) / (ea-1.)**4
1356 return 0., mu2, 0., mu4/mu2**2 - 3.
1358 def _entropy(self, a):
1359 return a / sinh(a) - log(tanh(a/2.0))
1361 def _rvs(self, a, size=None, random_state=None):
1362 # The discrete Laplace is equivalent to the two-sided geometric
1363 # distribution with PMF:
1364 # f(k) = (1 - alpha)/(1 + alpha) * alpha^abs(k)
1365 # Reference:
1366 # https://www.sciencedirect.com/science/
1367 # article/abs/pii/S0378375804003519
1368 # Furthermore, the two-sided geometric distribution is
1369 # equivalent to the difference between two iid geometric
1370 # distributions.
1371 # Reference (page 179):
1372 # https://pdfs.semanticscholar.org/61b3/
1373 # b99f466815808fd0d03f5d2791eea8b541a1.pdf
1374 # Thus, we can leverage the following:
1375 # 1) alpha = e^-a
1376 # 2) probability_of_success = 1 - alpha (Bernoulli trial)
1377 probOfSuccess = -np.expm1(-np.asarray(a))
1378 x = random_state.geometric(probOfSuccess, size=size)
1379 y = random_state.geometric(probOfSuccess, size=size)
1380 return x - y
1383dlaplace = dlaplace_gen(a=-np.inf,
1384 name='dlaplace', longname='A discrete Laplacian')
1387class skellam_gen(rv_discrete):
1388 r"""A Skellam discrete random variable.
1390 %(before_notes)s
1392 Notes
1393 -----
1394 Probability distribution of the difference of two correlated or
1395 uncorrelated Poisson random variables.
1397 Let :math:`k_1` and :math:`k_2` be two Poisson-distributed r.v. with
1398 expected values :math:`\lambda_1` and :math:`\lambda_2`. Then,
1399 :math:`k_1 - k_2` follows a Skellam distribution with parameters
1400 :math:`\mu_1 = \lambda_1 - \rho \sqrt{\lambda_1 \lambda_2}` and
1401 :math:`\mu_2 = \lambda_2 - \rho \sqrt{\lambda_1 \lambda_2}`, where
1402 :math:`\rho` is the correlation coefficient between :math:`k_1` and
1403 :math:`k_2`. If the two Poisson-distributed r.v. are independent then
1404 :math:`\rho = 0`.
1406 Parameters :math:`\mu_1` and :math:`\mu_2` must be strictly positive.
1408 For details see: https://en.wikipedia.org/wiki/Skellam_distribution
1410 `skellam` takes :math:`\mu_1` and :math:`\mu_2` as shape parameters.
1412 %(after_notes)s
1414 %(example)s
1416 """
1417 def _shape_info(self):
1418 return [_ShapeInfo("mu1", False, (0, np.inf), (False, False)),
1419 _ShapeInfo("mu2", False, (0, np.inf), (False, False))]
1421 def _rvs(self, mu1, mu2, size=None, random_state=None):
1422 n = size
1423 return (random_state.poisson(mu1, n) -
1424 random_state.poisson(mu2, n))
1426 def _pmf(self, x, mu1, mu2):
1427 with np.errstate(over='ignore'): # see gh-17432
1428 px = np.where(x < 0,
1429 _boost._ncx2_pdf(2*mu2, 2*(1-x), 2*mu1)*2,
1430 _boost._ncx2_pdf(2*mu1, 2*(1+x), 2*mu2)*2)
1431 # ncx2.pdf() returns nan's for extremely low probabilities
1432 return px
1434 def _cdf(self, x, mu1, mu2):
1435 x = floor(x)
1436 with np.errstate(over='ignore'): # see gh-17432
1437 px = np.where(x < 0,
1438 _boost._ncx2_cdf(2*mu2, -2*x, 2*mu1),
1439 1 - _boost._ncx2_cdf(2*mu1, 2*(x+1), 2*mu2))
1440 return px
1442 def _stats(self, mu1, mu2):
1443 mean = mu1 - mu2
1444 var = mu1 + mu2
1445 g1 = mean / sqrt((var)**3)
1446 g2 = 1 / var
1447 return mean, var, g1, g2
1450skellam = skellam_gen(a=-np.inf, name="skellam", longname='A Skellam')
1453class yulesimon_gen(rv_discrete):
1454 r"""A Yule-Simon discrete random variable.
1456 %(before_notes)s
1458 Notes
1459 -----
1461 The probability mass function for the `yulesimon` is:
1463 .. math::
1465 f(k) = \alpha B(k, \alpha+1)
1467 for :math:`k=1,2,3,...`, where :math:`\alpha>0`.
1468 Here :math:`B` refers to the `scipy.special.beta` function.
1470 The sampling of random variates is based on pg 553, Section 6.3 of [1]_.
1471 Our notation maps to the referenced logic via :math:`\alpha=a-1`.
1473 For details see the wikipedia entry [2]_.
1475 References
1476 ----------
1477 .. [1] Devroye, Luc. "Non-uniform Random Variate Generation",
1478 (1986) Springer, New York.
1480 .. [2] https://en.wikipedia.org/wiki/Yule-Simon_distribution
1482 %(after_notes)s
1484 %(example)s
1486 """
1487 def _shape_info(self):
1488 return [_ShapeInfo("alpha", False, (0, np.inf), (False, False))]
1490 def _rvs(self, alpha, size=None, random_state=None):
1491 E1 = random_state.standard_exponential(size)
1492 E2 = random_state.standard_exponential(size)
1493 ans = ceil(-E1 / log1p(-exp(-E2 / alpha)))
1494 return ans
1496 def _pmf(self, x, alpha):
1497 return alpha * special.beta(x, alpha + 1)
1499 def _argcheck(self, alpha):
1500 return (alpha > 0)
1502 def _logpmf(self, x, alpha):
1503 return log(alpha) + special.betaln(x, alpha + 1)
1505 def _cdf(self, x, alpha):
1506 return 1 - x * special.beta(x, alpha + 1)
1508 def _sf(self, x, alpha):
1509 return x * special.beta(x, alpha + 1)
1511 def _logsf(self, x, alpha):
1512 return log(x) + special.betaln(x, alpha + 1)
1514 def _stats(self, alpha):
1515 mu = np.where(alpha <= 1, np.inf, alpha / (alpha - 1))
1516 mu2 = np.where(alpha > 2,
1517 alpha**2 / ((alpha - 2.0) * (alpha - 1)**2),
1518 np.inf)
1519 mu2 = np.where(alpha <= 1, np.nan, mu2)
1520 g1 = np.where(alpha > 3,
1521 sqrt(alpha - 2) * (alpha + 1)**2 / (alpha * (alpha - 3)),
1522 np.inf)
1523 g1 = np.where(alpha <= 2, np.nan, g1)
1524 g2 = np.where(alpha > 4,
1525 alpha + 3 + ((alpha**3 - 49 * alpha - 22) /
1526 (alpha * (alpha - 4) * (alpha - 3))),
1527 np.inf)
1528 g2 = np.where(alpha <= 2, np.nan, g2)
1529 return mu, mu2, g1, g2
1532yulesimon = yulesimon_gen(name='yulesimon', a=1)
1535def _vectorize_rvs_over_shapes(_rvs1):
1536 """Decorator that vectorizes _rvs method to work on ndarray shapes"""
1537 # _rvs1 must be a _function_ that accepts _scalar_ args as positional
1538 # arguments, `size` and `random_state` as keyword arguments.
1539 # _rvs1 must return a random variate array with shape `size`. If `size` is
1540 # None, _rvs1 must return a scalar.
1541 # When applied to _rvs1, this decorator broadcasts ndarray args
1542 # and loops over them, calling _rvs1 for each set of scalar args.
1543 # For usage example, see _nchypergeom_gen
1544 def _rvs(*args, size, random_state):
1545 _rvs1_size, _rvs1_indices = _check_shape(args[0].shape, size)
1547 size = np.array(size)
1548 _rvs1_size = np.array(_rvs1_size)
1549 _rvs1_indices = np.array(_rvs1_indices)
1551 if np.all(_rvs1_indices): # all args are scalars
1552 return _rvs1(*args, size, random_state)
1554 out = np.empty(size)
1556 # out.shape can mix dimensions associated with arg_shape and _rvs1_size
1557 # Sort them to arg_shape + _rvs1_size for easy indexing of dimensions
1558 # corresponding with the different sets of scalar args
1559 j0 = np.arange(out.ndim)
1560 j1 = np.hstack((j0[~_rvs1_indices], j0[_rvs1_indices]))
1561 out = np.moveaxis(out, j1, j0)
1563 for i in np.ndindex(*size[~_rvs1_indices]):
1564 # arg can be squeezed because singleton dimensions will be
1565 # associated with _rvs1_size, not arg_shape per _check_shape
1566 out[i] = _rvs1(*[np.squeeze(arg)[i] for arg in args],
1567 _rvs1_size, random_state)
1569 return np.moveaxis(out, j0, j1) # move axes back before returning
1570 return _rvs
1573class _nchypergeom_gen(rv_discrete):
1574 r"""A noncentral hypergeometric discrete random variable.
1576 For subclassing by nchypergeom_fisher_gen and nchypergeom_wallenius_gen.
1578 """
1580 rvs_name = None
1581 dist = None
1583 def _shape_info(self):
1584 return [_ShapeInfo("M", True, (0, np.inf), (True, False)),
1585 _ShapeInfo("n", True, (0, np.inf), (True, False)),
1586 _ShapeInfo("N", True, (0, np.inf), (True, False)),
1587 _ShapeInfo("odds", False, (0, np.inf), (False, False))]
1589 def _get_support(self, M, n, N, odds):
1590 N, m1, n = M, n, N # follow Wikipedia notation
1591 m2 = N - m1
1592 x_min = np.maximum(0, n - m2)
1593 x_max = np.minimum(n, m1)
1594 return x_min, x_max
1596 def _argcheck(self, M, n, N, odds):
1597 M, n = np.asarray(M), np.asarray(n),
1598 N, odds = np.asarray(N), np.asarray(odds)
1599 cond1 = (M.astype(int) == M) & (M >= 0)
1600 cond2 = (n.astype(int) == n) & (n >= 0)
1601 cond3 = (N.astype(int) == N) & (N >= 0)
1602 cond4 = odds > 0
1603 cond5 = N <= M
1604 cond6 = n <= M
1605 return cond1 & cond2 & cond3 & cond4 & cond5 & cond6
1607 def _rvs(self, M, n, N, odds, size=None, random_state=None):
1609 @_vectorize_rvs_over_shapes
1610 def _rvs1(M, n, N, odds, size, random_state):
1611 length = np.prod(size)
1612 urn = _PyStochasticLib3()
1613 rv_gen = getattr(urn, self.rvs_name)
1614 rvs = rv_gen(N, n, M, odds, length, random_state)
1615 rvs = rvs.reshape(size)
1616 return rvs
1618 return _rvs1(M, n, N, odds, size=size, random_state=random_state)
1620 def _pmf(self, x, M, n, N, odds):
1622 x, M, n, N, odds = np.broadcast_arrays(x, M, n, N, odds)
1623 if x.size == 0: # np.vectorize doesn't work with zero size input
1624 return np.empty_like(x)
1626 @np.vectorize
1627 def _pmf1(x, M, n, N, odds):
1628 urn = self.dist(N, n, M, odds, 1e-12)
1629 return urn.probability(x)
1631 return _pmf1(x, M, n, N, odds)
1633 def _stats(self, M, n, N, odds, moments):
1635 @np.vectorize
1636 def _moments1(M, n, N, odds):
1637 urn = self.dist(N, n, M, odds, 1e-12)
1638 return urn.moments()
1640 m, v = (_moments1(M, n, N, odds) if ("m" in moments or "v" in moments)
1641 else (None, None))
1642 s, k = None, None
1643 return m, v, s, k
1646class nchypergeom_fisher_gen(_nchypergeom_gen):
1647 r"""A Fisher's noncentral hypergeometric discrete random variable.
1649 Fisher's noncentral hypergeometric distribution models drawing objects of
1650 two types from a bin. `M` is the total number of objects, `n` is the
1651 number of Type I objects, and `odds` is the odds ratio: the odds of
1652 selecting a Type I object rather than a Type II object when there is only
1653 one object of each type.
1654 The random variate represents the number of Type I objects drawn if we
1655 take a handful of objects from the bin at once and find out afterwards
1656 that we took `N` objects.
1658 %(before_notes)s
1660 See Also
1661 --------
1662 nchypergeom_wallenius, hypergeom, nhypergeom
1664 Notes
1665 -----
1666 Let mathematical symbols :math:`N`, :math:`n`, and :math:`M` correspond
1667 with parameters `N`, `n`, and `M` (respectively) as defined above.
1669 The probability mass function is defined as
1671 .. math::
1673 p(x; M, n, N, \omega) =
1674 \frac{\binom{n}{x}\binom{M - n}{N-x}\omega^x}{P_0},
1676 for
1677 :math:`x \in [x_l, x_u]`,
1678 :math:`M \in {\mathbb N}`,
1679 :math:`n \in [0, M]`,
1680 :math:`N \in [0, M]`,
1681 :math:`\omega > 0`,
1682 where
1683 :math:`x_l = \max(0, N - (M - n))`,
1684 :math:`x_u = \min(N, n)`,
1686 .. math::
1688 P_0 = \sum_{y=x_l}^{x_u} \binom{n}{y}\binom{M - n}{N-y}\omega^y,
1690 and the binomial coefficients are defined as
1692 .. math:: \binom{n}{k} \equiv \frac{n!}{k! (n - k)!}.
1694 `nchypergeom_fisher` uses the BiasedUrn package by Agner Fog with
1695 permission for it to be distributed under SciPy's license.
1697 The symbols used to denote the shape parameters (`N`, `n`, and `M`) are not
1698 universally accepted; they are chosen for consistency with `hypergeom`.
1700 Note that Fisher's noncentral hypergeometric distribution is distinct
1701 from Wallenius' noncentral hypergeometric distribution, which models
1702 drawing a pre-determined `N` objects from a bin one by one.
1703 When the odds ratio is unity, however, both distributions reduce to the
1704 ordinary hypergeometric distribution.
1706 %(after_notes)s
1708 References
1709 ----------
1710 .. [1] Agner Fog, "Biased Urn Theory".
1711 https://cran.r-project.org/web/packages/BiasedUrn/vignettes/UrnTheory.pdf
1713 .. [2] "Fisher's noncentral hypergeometric distribution", Wikipedia,
1714 https://en.wikipedia.org/wiki/Fisher's_noncentral_hypergeometric_distribution
1716 %(example)s
1718 """
1720 rvs_name = "rvs_fisher"
1721 dist = _PyFishersNCHypergeometric
1724nchypergeom_fisher = nchypergeom_fisher_gen(
1725 name='nchypergeom_fisher',
1726 longname="A Fisher's noncentral hypergeometric")
1729class nchypergeom_wallenius_gen(_nchypergeom_gen):
1730 r"""A Wallenius' noncentral hypergeometric discrete random variable.
1732 Wallenius' noncentral hypergeometric distribution models drawing objects of
1733 two types from a bin. `M` is the total number of objects, `n` is the
1734 number of Type I objects, and `odds` is the odds ratio: the odds of
1735 selecting a Type I object rather than a Type II object when there is only
1736 one object of each type.
1737 The random variate represents the number of Type I objects drawn if we
1738 draw a pre-determined `N` objects from a bin one by one.
1740 %(before_notes)s
1742 See Also
1743 --------
1744 nchypergeom_fisher, hypergeom, nhypergeom
1746 Notes
1747 -----
1748 Let mathematical symbols :math:`N`, :math:`n`, and :math:`M` correspond
1749 with parameters `N`, `n`, and `M` (respectively) as defined above.
1751 The probability mass function is defined as
1753 .. math::
1755 p(x; N, n, M) = \binom{n}{x} \binom{M - n}{N-x}
1756 \int_0^1 \left(1-t^{\omega/D}\right)^x\left(1-t^{1/D}\right)^{N-x} dt
1758 for
1759 :math:`x \in [x_l, x_u]`,
1760 :math:`M \in {\mathbb N}`,
1761 :math:`n \in [0, M]`,
1762 :math:`N \in [0, M]`,
1763 :math:`\omega > 0`,
1764 where
1765 :math:`x_l = \max(0, N - (M - n))`,
1766 :math:`x_u = \min(N, n)`,
1768 .. math::
1770 D = \omega(n - x) + ((M - n)-(N-x)),
1772 and the binomial coefficients are defined as
1774 .. math:: \binom{n}{k} \equiv \frac{n!}{k! (n - k)!}.
1776 `nchypergeom_wallenius` uses the BiasedUrn package by Agner Fog with
1777 permission for it to be distributed under SciPy's license.
1779 The symbols used to denote the shape parameters (`N`, `n`, and `M`) are not
1780 universally accepted; they are chosen for consistency with `hypergeom`.
1782 Note that Wallenius' noncentral hypergeometric distribution is distinct
1783 from Fisher's noncentral hypergeometric distribution, which models
1784 take a handful of objects from the bin at once, finding out afterwards
1785 that `N` objects were taken.
1786 When the odds ratio is unity, however, both distributions reduce to the
1787 ordinary hypergeometric distribution.
1789 %(after_notes)s
1791 References
1792 ----------
1793 .. [1] Agner Fog, "Biased Urn Theory".
1794 https://cran.r-project.org/web/packages/BiasedUrn/vignettes/UrnTheory.pdf
1796 .. [2] "Wallenius' noncentral hypergeometric distribution", Wikipedia,
1797 https://en.wikipedia.org/wiki/Wallenius'_noncentral_hypergeometric_distribution
1799 %(example)s
1801 """
1803 rvs_name = "rvs_wallenius"
1804 dist = _PyWalleniusNCHypergeometric
1807nchypergeom_wallenius = nchypergeom_wallenius_gen(
1808 name='nchypergeom_wallenius',
1809 longname="A Wallenius' noncentral hypergeometric")
1812# Collect names of classes and objects in this module.
1813pairs = list(globals().copy().items())
1814_distn_names, _distn_gen_names = get_distribution_names(pairs, rv_discrete)
1816__all__ = _distn_names + _distn_gen_names