Coverage for /usr/lib/python3/dist-packages/scipy/stats/_distn_infrastructure.py: 26%
1573 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 scipy._lib._util import getfullargspec_no_self as _getfullargspec
7import sys
8import keyword
9import re
10import types
11import warnings
12from itertools import zip_longest
14from scipy._lib import doccer
15from ._distr_params import distcont, distdiscrete
16from scipy._lib._util import check_random_state
18from scipy.special import comb, entr
21# for root finding for continuous distribution ppf, and maximum likelihood
22# estimation
23from scipy import optimize
25# for functions of continuous distributions (e.g. moments, entropy, cdf)
26from scipy import integrate
28# to approximate the pdf of a continuous distribution given its cdf
29from scipy._lib._finite_differences import _derivative
31# for scipy.stats.entropy. Attempts to import just that function or file
32# have cause import problems
33from scipy import stats
35from numpy import (arange, putmask, ones, shape, ndarray, zeros, floor,
36 logical_and, log, sqrt, place, argmax, vectorize, asarray,
37 nan, inf, isinf, NINF, empty)
39import numpy as np
40from ._constants import _XMAX, _LOGXMAX
41from ._censored_data import CensoredData
42from scipy.stats._warnings_errors import FitError
44# These are the docstring parts used for substitution in specific
45# distribution docstrings
47docheaders = {'methods': """\nMethods\n-------\n""",
48 'notes': """\nNotes\n-----\n""",
49 'examples': """\nExamples\n--------\n"""}
51_doc_rvs = """\
52rvs(%(shapes)s, loc=0, scale=1, size=1, random_state=None)
53 Random variates.
54"""
55_doc_pdf = """\
56pdf(x, %(shapes)s, loc=0, scale=1)
57 Probability density function.
58"""
59_doc_logpdf = """\
60logpdf(x, %(shapes)s, loc=0, scale=1)
61 Log of the probability density function.
62"""
63_doc_pmf = """\
64pmf(k, %(shapes)s, loc=0, scale=1)
65 Probability mass function.
66"""
67_doc_logpmf = """\
68logpmf(k, %(shapes)s, loc=0, scale=1)
69 Log of the probability mass function.
70"""
71_doc_cdf = """\
72cdf(x, %(shapes)s, loc=0, scale=1)
73 Cumulative distribution function.
74"""
75_doc_logcdf = """\
76logcdf(x, %(shapes)s, loc=0, scale=1)
77 Log of the cumulative distribution function.
78"""
79_doc_sf = """\
80sf(x, %(shapes)s, loc=0, scale=1)
81 Survival function (also defined as ``1 - cdf``, but `sf` is sometimes more accurate).
82"""
83_doc_logsf = """\
84logsf(x, %(shapes)s, loc=0, scale=1)
85 Log of the survival function.
86"""
87_doc_ppf = """\
88ppf(q, %(shapes)s, loc=0, scale=1)
89 Percent point function (inverse of ``cdf`` --- percentiles).
90"""
91_doc_isf = """\
92isf(q, %(shapes)s, loc=0, scale=1)
93 Inverse survival function (inverse of ``sf``).
94"""
95_doc_moment = """\
96moment(order, %(shapes)s, loc=0, scale=1)
97 Non-central moment of the specified order.
98"""
99_doc_stats = """\
100stats(%(shapes)s, loc=0, scale=1, moments='mv')
101 Mean('m'), variance('v'), skew('s'), and/or kurtosis('k').
102"""
103_doc_entropy = """\
104entropy(%(shapes)s, loc=0, scale=1)
105 (Differential) entropy of the RV.
106"""
107_doc_fit = """\
108fit(data)
109 Parameter estimates for generic data.
110 See `scipy.stats.rv_continuous.fit <https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.fit.html#scipy.stats.rv_continuous.fit>`__ for detailed documentation of the
111 keyword arguments.
112"""
113_doc_expect = """\
114expect(func, args=(%(shapes_)s), loc=0, scale=1, lb=None, ub=None, conditional=False, **kwds)
115 Expected value of a function (of one argument) with respect to the distribution.
116"""
117_doc_expect_discrete = """\
118expect(func, args=(%(shapes_)s), loc=0, lb=None, ub=None, conditional=False)
119 Expected value of a function (of one argument) with respect to the distribution.
120"""
121_doc_median = """\
122median(%(shapes)s, loc=0, scale=1)
123 Median of the distribution.
124"""
125_doc_mean = """\
126mean(%(shapes)s, loc=0, scale=1)
127 Mean of the distribution.
128"""
129_doc_var = """\
130var(%(shapes)s, loc=0, scale=1)
131 Variance of the distribution.
132"""
133_doc_std = """\
134std(%(shapes)s, loc=0, scale=1)
135 Standard deviation of the distribution.
136"""
137_doc_interval = """\
138interval(confidence, %(shapes)s, loc=0, scale=1)
139 Confidence interval with equal areas around the median.
140"""
141_doc_allmethods = ''.join([docheaders['methods'], _doc_rvs, _doc_pdf,
142 _doc_logpdf, _doc_cdf, _doc_logcdf, _doc_sf,
143 _doc_logsf, _doc_ppf, _doc_isf, _doc_moment,
144 _doc_stats, _doc_entropy, _doc_fit,
145 _doc_expect, _doc_median,
146 _doc_mean, _doc_var, _doc_std, _doc_interval])
148_doc_default_longsummary = """\
149As an instance of the `rv_continuous` class, `%(name)s` object inherits from it
150a collection of generic methods (see below for the full list),
151and completes them with details specific for this particular distribution.
152"""
154_doc_default_frozen_note = """
155Alternatively, the object may be called (as a function) to fix the shape,
156location, and scale parameters returning a "frozen" continuous RV object:
158rv = %(name)s(%(shapes)s, loc=0, scale=1)
159 - Frozen RV object with the same methods but holding the given shape,
160 location, and scale fixed.
161"""
162_doc_default_example = """\
163Examples
164--------
165>>> import numpy as np
166>>> from scipy.stats import %(name)s
167>>> import matplotlib.pyplot as plt
168>>> fig, ax = plt.subplots(1, 1)
170Calculate the first four moments:
172%(set_vals_stmt)s
173>>> mean, var, skew, kurt = %(name)s.stats(%(shapes)s, moments='mvsk')
175Display the probability density function (``pdf``):
177>>> x = np.linspace(%(name)s.ppf(0.01, %(shapes)s),
178... %(name)s.ppf(0.99, %(shapes)s), 100)
179>>> ax.plot(x, %(name)s.pdf(x, %(shapes)s),
180... 'r-', lw=5, alpha=0.6, label='%(name)s pdf')
182Alternatively, the distribution object can be called (as a function)
183to fix the shape, location and scale parameters. This returns a "frozen"
184RV object holding the given parameters fixed.
186Freeze the distribution and display the frozen ``pdf``:
188>>> rv = %(name)s(%(shapes)s)
189>>> ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')
191Check accuracy of ``cdf`` and ``ppf``:
193>>> vals = %(name)s.ppf([0.001, 0.5, 0.999], %(shapes)s)
194>>> np.allclose([0.001, 0.5, 0.999], %(name)s.cdf(vals, %(shapes)s))
195True
197Generate random numbers:
199>>> r = %(name)s.rvs(%(shapes)s, size=1000)
201And compare the histogram:
203>>> ax.hist(r, density=True, bins='auto', histtype='stepfilled', alpha=0.2)
204>>> ax.set_xlim([x[0], x[-1]])
205>>> ax.legend(loc='best', frameon=False)
206>>> plt.show()
208"""
210_doc_default_locscale = """\
211The probability density above is defined in the "standardized" form. To shift
212and/or scale the distribution use the ``loc`` and ``scale`` parameters.
213Specifically, ``%(name)s.pdf(x, %(shapes)s, loc, scale)`` is identically
214equivalent to ``%(name)s.pdf(y, %(shapes)s) / scale`` with
215``y = (x - loc) / scale``. Note that shifting the location of a distribution
216does not make it a "noncentral" distribution; noncentral generalizations of
217some distributions are available in separate classes.
218"""
220_doc_default = ''.join([_doc_default_longsummary,
221 _doc_allmethods,
222 '\n',
223 _doc_default_example])
225_doc_default_before_notes = ''.join([_doc_default_longsummary,
226 _doc_allmethods])
228docdict = {
229 'rvs': _doc_rvs,
230 'pdf': _doc_pdf,
231 'logpdf': _doc_logpdf,
232 'cdf': _doc_cdf,
233 'logcdf': _doc_logcdf,
234 'sf': _doc_sf,
235 'logsf': _doc_logsf,
236 'ppf': _doc_ppf,
237 'isf': _doc_isf,
238 'stats': _doc_stats,
239 'entropy': _doc_entropy,
240 'fit': _doc_fit,
241 'moment': _doc_moment,
242 'expect': _doc_expect,
243 'interval': _doc_interval,
244 'mean': _doc_mean,
245 'std': _doc_std,
246 'var': _doc_var,
247 'median': _doc_median,
248 'allmethods': _doc_allmethods,
249 'longsummary': _doc_default_longsummary,
250 'frozennote': _doc_default_frozen_note,
251 'example': _doc_default_example,
252 'default': _doc_default,
253 'before_notes': _doc_default_before_notes,
254 'after_notes': _doc_default_locscale
255}
257# Reuse common content between continuous and discrete docs, change some
258# minor bits.
259docdict_discrete = docdict.copy()
261docdict_discrete['pmf'] = _doc_pmf
262docdict_discrete['logpmf'] = _doc_logpmf
263docdict_discrete['expect'] = _doc_expect_discrete
264_doc_disc_methods = ['rvs', 'pmf', 'logpmf', 'cdf', 'logcdf', 'sf', 'logsf',
265 'ppf', 'isf', 'stats', 'entropy', 'expect', 'median',
266 'mean', 'var', 'std', 'interval']
267for obj in _doc_disc_methods:
268 docdict_discrete[obj] = docdict_discrete[obj].replace(', scale=1', '')
270_doc_disc_methods_err_varname = ['cdf', 'logcdf', 'sf', 'logsf']
271for obj in _doc_disc_methods_err_varname:
272 docdict_discrete[obj] = docdict_discrete[obj].replace('(x, ', '(k, ')
274docdict_discrete.pop('pdf')
275docdict_discrete.pop('logpdf')
277_doc_allmethods = ''.join([docdict_discrete[obj] for obj in _doc_disc_methods])
278docdict_discrete['allmethods'] = docheaders['methods'] + _doc_allmethods
280docdict_discrete['longsummary'] = _doc_default_longsummary.replace(
281 'rv_continuous', 'rv_discrete')
283_doc_default_frozen_note = """
284Alternatively, the object may be called (as a function) to fix the shape and
285location parameters returning a "frozen" discrete RV object:
287rv = %(name)s(%(shapes)s, loc=0)
288 - Frozen RV object with the same methods but holding the given shape and
289 location fixed.
290"""
291docdict_discrete['frozennote'] = _doc_default_frozen_note
293_doc_default_discrete_example = """\
294Examples
295--------
296>>> import numpy as np
297>>> from scipy.stats import %(name)s
298>>> import matplotlib.pyplot as plt
299>>> fig, ax = plt.subplots(1, 1)
301Calculate the first four moments:
303%(set_vals_stmt)s
304>>> mean, var, skew, kurt = %(name)s.stats(%(shapes)s, moments='mvsk')
306Display the probability mass function (``pmf``):
308>>> x = np.arange(%(name)s.ppf(0.01, %(shapes)s),
309... %(name)s.ppf(0.99, %(shapes)s))
310>>> ax.plot(x, %(name)s.pmf(x, %(shapes)s), 'bo', ms=8, label='%(name)s pmf')
311>>> ax.vlines(x, 0, %(name)s.pmf(x, %(shapes)s), colors='b', lw=5, alpha=0.5)
313Alternatively, the distribution object can be called (as a function)
314to fix the shape and location. This returns a "frozen" RV object holding
315the given parameters fixed.
317Freeze the distribution and display the frozen ``pmf``:
319>>> rv = %(name)s(%(shapes)s)
320>>> ax.vlines(x, 0, rv.pmf(x), colors='k', linestyles='-', lw=1,
321... label='frozen pmf')
322>>> ax.legend(loc='best', frameon=False)
323>>> plt.show()
325Check accuracy of ``cdf`` and ``ppf``:
327>>> prob = %(name)s.cdf(x, %(shapes)s)
328>>> np.allclose(x, %(name)s.ppf(prob, %(shapes)s))
329True
331Generate random numbers:
333>>> r = %(name)s.rvs(%(shapes)s, size=1000)
334"""
337_doc_default_discrete_locscale = """\
338The probability mass function above is defined in the "standardized" form.
339To shift distribution use the ``loc`` parameter.
340Specifically, ``%(name)s.pmf(k, %(shapes)s, loc)`` is identically
341equivalent to ``%(name)s.pmf(k - loc, %(shapes)s)``.
342"""
344docdict_discrete['example'] = _doc_default_discrete_example
345docdict_discrete['after_notes'] = _doc_default_discrete_locscale
347_doc_default_before_notes = ''.join([docdict_discrete['longsummary'],
348 docdict_discrete['allmethods']])
349docdict_discrete['before_notes'] = _doc_default_before_notes
351_doc_default_disc = ''.join([docdict_discrete['longsummary'],
352 docdict_discrete['allmethods'],
353 docdict_discrete['frozennote'],
354 docdict_discrete['example']])
355docdict_discrete['default'] = _doc_default_disc
357# clean up all the separate docstring elements, we do not need them anymore
358for obj in [s for s in dir() if s.startswith('_doc_')]:
359 exec('del ' + obj)
360del obj
363def _moment(data, n, mu=None):
364 if mu is None:
365 mu = data.mean()
366 return ((data - mu)**n).mean()
369def _moment_from_stats(n, mu, mu2, g1, g2, moment_func, args):
370 if (n == 0):
371 return 1.0
372 elif (n == 1):
373 if mu is None:
374 val = moment_func(1, *args)
375 else:
376 val = mu
377 elif (n == 2):
378 if mu2 is None or mu is None:
379 val = moment_func(2, *args)
380 else:
381 val = mu2 + mu*mu
382 elif (n == 3):
383 if g1 is None or mu2 is None or mu is None:
384 val = moment_func(3, *args)
385 else:
386 mu3 = g1 * np.power(mu2, 1.5) # 3rd central moment
387 val = mu3+3*mu*mu2+mu*mu*mu # 3rd non-central moment
388 elif (n == 4):
389 if g1 is None or g2 is None or mu2 is None or mu is None:
390 val = moment_func(4, *args)
391 else:
392 mu4 = (g2+3.0)*(mu2**2.0) # 4th central moment
393 mu3 = g1*np.power(mu2, 1.5) # 3rd central moment
394 val = mu4+4*mu*mu3+6*mu*mu*mu2+mu*mu*mu*mu
395 else:
396 val = moment_func(n, *args)
398 return val
401def _skew(data):
402 """
403 skew is third central moment / variance**(1.5)
404 """
405 data = np.ravel(data)
406 mu = data.mean()
407 m2 = ((data - mu)**2).mean()
408 m3 = ((data - mu)**3).mean()
409 return m3 / np.power(m2, 1.5)
412def _kurtosis(data):
413 """kurtosis is fourth central moment / variance**2 - 3."""
414 data = np.ravel(data)
415 mu = data.mean()
416 m2 = ((data - mu)**2).mean()
417 m4 = ((data - mu)**4).mean()
418 return m4 / m2**2 - 3
421def _fit_determine_optimizer(optimizer):
422 if not callable(optimizer) and isinstance(optimizer, str):
423 if not optimizer.startswith('fmin_'):
424 optimizer = "fmin_"+optimizer
425 if optimizer == 'fmin_':
426 optimizer = 'fmin'
427 try:
428 optimizer = getattr(optimize, optimizer)
429 except AttributeError as e:
430 raise ValueError("%s is not a valid optimizer" % optimizer) from e
431 return optimizer
434def _sum_finite(x):
435 """
436 For a 1D array x, return a tuple containing the sum of the
437 finite values of x and the number of nonfinite values.
439 This is a utility function used when evaluating the negative
440 loglikelihood for a distribution and an array of samples.
442 Examples
443 --------
444 >>> tot, nbad = _sum_finite(np.array([-2, -np.inf, 5, 1]))
445 >>> tot
446 4.0
447 >>> nbad
448 1
449 """
450 finite_x = np.isfinite(x)
451 bad_count = finite_x.size - np.count_nonzero(finite_x)
452 return np.sum(x[finite_x]), bad_count
455# Frozen RV class
456class rv_frozen:
458 def __init__(self, dist, *args, **kwds):
459 self.args = args
460 self.kwds = kwds
462 # create a new instance
463 self.dist = dist.__class__(**dist._updated_ctor_param())
465 shapes, _, _ = self.dist._parse_args(*args, **kwds)
466 self.a, self.b = self.dist._get_support(*shapes)
468 @property
469 def random_state(self):
470 return self.dist._random_state
472 @random_state.setter
473 def random_state(self, seed):
474 self.dist._random_state = check_random_state(seed)
476 def cdf(self, x):
477 return self.dist.cdf(x, *self.args, **self.kwds)
479 def logcdf(self, x):
480 return self.dist.logcdf(x, *self.args, **self.kwds)
482 def ppf(self, q):
483 return self.dist.ppf(q, *self.args, **self.kwds)
485 def isf(self, q):
486 return self.dist.isf(q, *self.args, **self.kwds)
488 def rvs(self, size=None, random_state=None):
489 kwds = self.kwds.copy()
490 kwds.update({'size': size, 'random_state': random_state})
491 return self.dist.rvs(*self.args, **kwds)
493 def sf(self, x):
494 return self.dist.sf(x, *self.args, **self.kwds)
496 def logsf(self, x):
497 return self.dist.logsf(x, *self.args, **self.kwds)
499 def stats(self, moments='mv'):
500 kwds = self.kwds.copy()
501 kwds.update({'moments': moments})
502 return self.dist.stats(*self.args, **kwds)
504 def median(self):
505 return self.dist.median(*self.args, **self.kwds)
507 def mean(self):
508 return self.dist.mean(*self.args, **self.kwds)
510 def var(self):
511 return self.dist.var(*self.args, **self.kwds)
513 def std(self):
514 return self.dist.std(*self.args, **self.kwds)
516 def moment(self, order=None):
517 return self.dist.moment(order, *self.args, **self.kwds)
519 def entropy(self):
520 return self.dist.entropy(*self.args, **self.kwds)
522 def interval(self, confidence=None):
523 return self.dist.interval(confidence, *self.args, **self.kwds)
525 def expect(self, func=None, lb=None, ub=None, conditional=False, **kwds):
526 # expect method only accepts shape parameters as positional args
527 # hence convert self.args, self.kwds, also loc/scale
528 # See the .expect method docstrings for the meaning of
529 # other parameters.
530 a, loc, scale = self.dist._parse_args(*self.args, **self.kwds)
531 if isinstance(self.dist, rv_discrete):
532 return self.dist.expect(func, a, loc, lb, ub, conditional, **kwds)
533 else:
534 return self.dist.expect(func, a, loc, scale, lb, ub,
535 conditional, **kwds)
537 def support(self):
538 return self.dist.support(*self.args, **self.kwds)
541class rv_discrete_frozen(rv_frozen):
543 def pmf(self, k):
544 return self.dist.pmf(k, *self.args, **self.kwds)
546 def logpmf(self, k): # No error
547 return self.dist.logpmf(k, *self.args, **self.kwds)
550class rv_continuous_frozen(rv_frozen):
552 def pdf(self, x):
553 return self.dist.pdf(x, *self.args, **self.kwds)
555 def logpdf(self, x):
556 return self.dist.logpdf(x, *self.args, **self.kwds)
559def argsreduce(cond, *args):
560 """Clean arguments to:
562 1. Ensure all arguments are iterable (arrays of dimension at least one
563 2. If cond != True and size > 1, ravel(args[i]) where ravel(condition) is
564 True, in 1D.
566 Return list of processed arguments.
568 Examples
569 --------
570 >>> import numpy as np
571 >>> rng = np.random.default_rng()
572 >>> A = rng.random((4, 5))
573 >>> B = 2
574 >>> C = rng.random((1, 5))
575 >>> cond = np.ones(A.shape)
576 >>> [A1, B1, C1] = argsreduce(cond, A, B, C)
577 >>> A1.shape
578 (4, 5)
579 >>> B1.shape
580 (1,)
581 >>> C1.shape
582 (1, 5)
583 >>> cond[2,:] = 0
584 >>> [A1, B1, C1] = argsreduce(cond, A, B, C)
585 >>> A1.shape
586 (15,)
587 >>> B1.shape
588 (1,)
589 >>> C1.shape
590 (15,)
592 """
593 # some distributions assume arguments are iterable.
594 newargs = np.atleast_1d(*args)
596 # np.atleast_1d returns an array if only one argument, or a list of arrays
597 # if more than one argument.
598 if not isinstance(newargs, list):
599 newargs = [newargs, ]
601 if np.all(cond):
602 # broadcast arrays with cond
603 *newargs, cond = np.broadcast_arrays(*newargs, cond)
604 return [arg.ravel() for arg in newargs]
606 s = cond.shape
607 # np.extract returns flattened arrays, which are not broadcastable together
608 # unless they are either the same size or size == 1.
609 return [(arg if np.size(arg) == 1
610 else np.extract(cond, np.broadcast_to(arg, s)))
611 for arg in newargs]
614parse_arg_template = """
615def _parse_args(self, %(shape_arg_str)s %(locscale_in)s):
616 return (%(shape_arg_str)s), %(locscale_out)s
618def _parse_args_rvs(self, %(shape_arg_str)s %(locscale_in)s, size=None):
619 return self._argcheck_rvs(%(shape_arg_str)s %(locscale_out)s, size=size)
621def _parse_args_stats(self, %(shape_arg_str)s %(locscale_in)s, moments='mv'):
622 return (%(shape_arg_str)s), %(locscale_out)s, moments
623"""
626class rv_generic:
627 """Class which encapsulates common functionality between rv_discrete
628 and rv_continuous.
630 """
632 def __init__(self, seed=None):
633 super().__init__()
635 # figure out if _stats signature has 'moments' keyword
636 sig = _getfullargspec(self._stats)
637 self._stats_has_moments = ((sig.varkw is not None) or
638 ('moments' in sig.args) or
639 ('moments' in sig.kwonlyargs))
640 self._random_state = check_random_state(seed)
642 @property
643 def random_state(self):
644 """Get or set the generator object for generating random variates.
646 If `random_state` is None (or `np.random`), the
647 `numpy.random.RandomState` singleton is used.
648 If `random_state` is an int, a new ``RandomState`` instance is used,
649 seeded with `random_state`.
650 If `random_state` is already a ``Generator`` or ``RandomState``
651 instance, that instance is used.
653 """
654 return self._random_state
656 @random_state.setter
657 def random_state(self, seed):
658 self._random_state = check_random_state(seed)
660 def __setstate__(self, state):
661 try:
662 self.__dict__.update(state)
663 # attaches the dynamically created methods on each instance.
664 # if a subclass overrides rv_generic.__setstate__, or implements
665 # it's own _attach_methods, then it must make sure that
666 # _attach_argparser_methods is called.
667 self._attach_methods()
668 except ValueError:
669 # reconstitute an old pickle scipy<1.6, that contains
670 # (_ctor_param, random_state) as state
671 self._ctor_param = state[0]
672 self._random_state = state[1]
673 self.__init__()
675 def _attach_methods(self):
676 """Attaches dynamically created methods to the rv_* instance.
678 This method must be overridden by subclasses, and must itself call
679 _attach_argparser_methods. This method is called in __init__ in
680 subclasses, and in __setstate__
681 """
682 raise NotImplementedError
684 def _attach_argparser_methods(self):
685 """
686 Generates the argument-parsing functions dynamically and attaches
687 them to the instance.
689 Should be called from `_attach_methods`, typically in __init__ and
690 during unpickling (__setstate__)
691 """
692 ns = {}
693 exec(self._parse_arg_template, ns)
694 # NB: attach to the instance, not class
695 for name in ['_parse_args', '_parse_args_stats', '_parse_args_rvs']:
696 setattr(self, name, types.MethodType(ns[name], self))
698 def _construct_argparser(
699 self, meths_to_inspect, locscale_in, locscale_out):
700 """Construct the parser string for the shape arguments.
702 This method should be called in __init__ of a class for each
703 distribution. It creates the `_parse_arg_template` attribute that is
704 then used by `_attach_argparser_methods` to dynamically create and
705 attach the `_parse_args`, `_parse_args_stats`, `_parse_args_rvs`
706 methods to the instance.
708 If self.shapes is a non-empty string, interprets it as a
709 comma-separated list of shape parameters.
711 Otherwise inspects the call signatures of `meths_to_inspect`
712 and constructs the argument-parsing functions from these.
713 In this case also sets `shapes` and `numargs`.
714 """
716 if self.shapes:
717 # sanitize the user-supplied shapes
718 if not isinstance(self.shapes, str):
719 raise TypeError('shapes must be a string.')
721 shapes = self.shapes.replace(',', ' ').split()
723 for field in shapes:
724 if keyword.iskeyword(field):
725 raise SyntaxError('keywords cannot be used as shapes.')
726 if not re.match('^[_a-zA-Z][_a-zA-Z0-9]*$', field):
727 raise SyntaxError(
728 'shapes must be valid python identifiers')
729 else:
730 # find out the call signatures (_pdf, _cdf etc), deduce shape
731 # arguments. Generic methods only have 'self, x', any further args
732 # are shapes.
733 shapes_list = []
734 for meth in meths_to_inspect:
735 shapes_args = _getfullargspec(meth) # NB does not contain self
736 args = shapes_args.args[1:] # peel off 'x', too
738 if args:
739 shapes_list.append(args)
741 # *args or **kwargs are not allowed w/automatic shapes
742 if shapes_args.varargs is not None:
743 raise TypeError(
744 '*args are not allowed w/out explicit shapes')
745 if shapes_args.varkw is not None:
746 raise TypeError(
747 '**kwds are not allowed w/out explicit shapes')
748 if shapes_args.kwonlyargs:
749 raise TypeError(
750 'kwonly args are not allowed w/out explicit shapes')
751 if shapes_args.defaults is not None:
752 raise TypeError('defaults are not allowed for shapes')
754 if shapes_list:
755 shapes = shapes_list[0]
757 # make sure the signatures are consistent
758 for item in shapes_list:
759 if item != shapes:
760 raise TypeError('Shape arguments are inconsistent.')
761 else:
762 shapes = []
764 # have the arguments, construct the method from template
765 shapes_str = ', '.join(shapes) + ', ' if shapes else '' # NB: not None
766 dct = dict(shape_arg_str=shapes_str,
767 locscale_in=locscale_in,
768 locscale_out=locscale_out,
769 )
771 # this string is used by _attach_argparser_methods
772 self._parse_arg_template = parse_arg_template % dct
774 self.shapes = ', '.join(shapes) if shapes else None
775 if not hasattr(self, 'numargs'):
776 # allows more general subclassing with *args
777 self.numargs = len(shapes)
779 def _construct_doc(self, docdict, shapes_vals=None):
780 """Construct the instance docstring with string substitutions."""
781 tempdict = docdict.copy()
782 tempdict['name'] = self.name or 'distname'
783 tempdict['shapes'] = self.shapes or ''
785 if shapes_vals is None:
786 shapes_vals = ()
787 vals = ', '.join('%.3g' % val for val in shapes_vals)
788 tempdict['vals'] = vals
790 tempdict['shapes_'] = self.shapes or ''
791 if self.shapes and self.numargs == 1:
792 tempdict['shapes_'] += ','
794 if self.shapes:
795 tempdict['set_vals_stmt'] = f'>>> {self.shapes} = {vals}'
796 else:
797 tempdict['set_vals_stmt'] = ''
799 if self.shapes is None:
800 # remove shapes from call parameters if there are none
801 for item in ['default', 'before_notes']:
802 tempdict[item] = tempdict[item].replace(
803 "\n%(shapes)s : array_like\n shape parameters", "")
804 for i in range(2):
805 if self.shapes is None:
806 # necessary because we use %(shapes)s in two forms (w w/o ", ")
807 self.__doc__ = self.__doc__.replace("%(shapes)s, ", "")
808 try:
809 self.__doc__ = doccer.docformat(self.__doc__, tempdict)
810 except TypeError as e:
811 raise Exception("Unable to construct docstring for "
812 "distribution \"%s\": %s" %
813 (self.name, repr(e))) from e
815 # correct for empty shapes
816 self.__doc__ = self.__doc__.replace('(, ', '(').replace(', )', ')')
818 def _construct_default_doc(self, longname=None,
819 docdict=None, discrete='continuous'):
820 """Construct instance docstring from the default template."""
821 if longname is None:
822 longname = 'A'
823 self.__doc__ = ''.join([f'{longname} {discrete} random variable.',
824 '\n\n%(before_notes)s\n', docheaders['notes'],
825 '\n%(example)s'])
826 self._construct_doc(docdict)
828 def freeze(self, *args, **kwds):
829 """Freeze the distribution for the given arguments.
831 Parameters
832 ----------
833 arg1, arg2, arg3,... : array_like
834 The shape parameter(s) for the distribution. Should include all
835 the non-optional arguments, may include ``loc`` and ``scale``.
837 Returns
838 -------
839 rv_frozen : rv_frozen instance
840 The frozen distribution.
842 """
843 if isinstance(self, rv_continuous):
844 return rv_continuous_frozen(self, *args, **kwds)
845 else:
846 return rv_discrete_frozen(self, *args, **kwds)
848 def __call__(self, *args, **kwds):
849 return self.freeze(*args, **kwds)
850 __call__.__doc__ = freeze.__doc__
852 # The actual calculation functions (no basic checking need be done)
853 # If these are defined, the others won't be looked at.
854 # Otherwise, the other set can be defined.
855 def _stats(self, *args, **kwds):
856 return None, None, None, None
858 # Noncentral moments (also known as the moment about the origin).
859 # Expressed in LaTeX, munp would be $\mu'_{n}$, i.e. "mu-sub-n-prime".
860 # The primed mu is a widely used notation for the noncentral moment.
861 def _munp(self, n, *args):
862 # Silence floating point warnings from integration.
863 with np.errstate(all='ignore'):
864 vals = self.generic_moment(n, *args)
865 return vals
867 def _argcheck_rvs(self, *args, **kwargs):
868 # Handle broadcasting and size validation of the rvs method.
869 # Subclasses should not have to override this method.
870 # The rule is that if `size` is not None, then `size` gives the
871 # shape of the result (integer values of `size` are treated as
872 # tuples with length 1; i.e. `size=3` is the same as `size=(3,)`.)
873 #
874 # `args` is expected to contain the shape parameters (if any), the
875 # location and the scale in a flat tuple (e.g. if there are two
876 # shape parameters `a` and `b`, `args` will be `(a, b, loc, scale)`).
877 # The only keyword argument expected is 'size'.
878 size = kwargs.get('size', None)
879 all_bcast = np.broadcast_arrays(*args)
881 def squeeze_left(a):
882 while a.ndim > 0 and a.shape[0] == 1:
883 a = a[0]
884 return a
886 # Eliminate trivial leading dimensions. In the convention
887 # used by numpy's random variate generators, trivial leading
888 # dimensions are effectively ignored. In other words, when `size`
889 # is given, trivial leading dimensions of the broadcast parameters
890 # in excess of the number of dimensions in size are ignored, e.g.
891 # >>> np.random.normal([[1, 3, 5]], [[[[0.01]]]], size=3)
892 # array([ 1.00104267, 3.00422496, 4.99799278])
893 # If `size` is not given, the exact broadcast shape is preserved:
894 # >>> np.random.normal([[1, 3, 5]], [[[[0.01]]]])
895 # array([[[[ 1.00862899, 3.00061431, 4.99867122]]]])
896 #
897 all_bcast = [squeeze_left(a) for a in all_bcast]
898 bcast_shape = all_bcast[0].shape
899 bcast_ndim = all_bcast[0].ndim
901 if size is None:
902 size_ = bcast_shape
903 else:
904 size_ = tuple(np.atleast_1d(size))
906 # Check compatibility of size_ with the broadcast shape of all
907 # the parameters. This check is intended to be consistent with
908 # how the numpy random variate generators (e.g. np.random.normal,
909 # np.random.beta) handle their arguments. The rule is that, if size
910 # is given, it determines the shape of the output. Broadcasting
911 # can't change the output size.
913 # This is the standard broadcasting convention of extending the
914 # shape with fewer dimensions with enough dimensions of length 1
915 # so that the two shapes have the same number of dimensions.
916 ndiff = bcast_ndim - len(size_)
917 if ndiff < 0:
918 bcast_shape = (1,)*(-ndiff) + bcast_shape
919 elif ndiff > 0:
920 size_ = (1,)*ndiff + size_
922 # This compatibility test is not standard. In "regular" broadcasting,
923 # two shapes are compatible if for each dimension, the lengths are the
924 # same or one of the lengths is 1. Here, the length of a dimension in
925 # size_ must not be less than the corresponding length in bcast_shape.
926 ok = all([bcdim == 1 or bcdim == szdim
927 for (bcdim, szdim) in zip(bcast_shape, size_)])
928 if not ok:
929 raise ValueError("size does not match the broadcast shape of "
930 f"the parameters. {size}, {size_}, {bcast_shape}")
932 param_bcast = all_bcast[:-2]
933 loc_bcast = all_bcast[-2]
934 scale_bcast = all_bcast[-1]
936 return param_bcast, loc_bcast, scale_bcast, size_
938 # These are the methods you must define (standard form functions)
939 # NB: generic _pdf, _logpdf, _cdf are different for
940 # rv_continuous and rv_discrete hence are defined in there
941 def _argcheck(self, *args):
942 """Default check for correct values on args and keywords.
944 Returns condition array of 1's where arguments are correct and
945 0's where they are not.
947 """
948 cond = 1
949 for arg in args:
950 cond = logical_and(cond, (asarray(arg) > 0))
951 return cond
953 def _get_support(self, *args, **kwargs):
954 """Return the support of the (unscaled, unshifted) distribution.
956 *Must* be overridden by distributions which have support dependent
957 upon the shape parameters of the distribution. Any such override
958 *must not* set or change any of the class members, as these members
959 are shared amongst all instances of the distribution.
961 Parameters
962 ----------
963 arg1, arg2, ... : array_like
964 The shape parameter(s) for the distribution (see docstring of the
965 instance object for more information).
967 Returns
968 -------
969 a, b : numeric (float, or int or +/-np.inf)
970 end-points of the distribution's support for the specified
971 shape parameters.
972 """
973 return self.a, self.b
975 def _support_mask(self, x, *args):
976 a, b = self._get_support(*args)
977 with np.errstate(invalid='ignore'):
978 return (a <= x) & (x <= b)
980 def _open_support_mask(self, x, *args):
981 a, b = self._get_support(*args)
982 with np.errstate(invalid='ignore'):
983 return (a < x) & (x < b)
985 def _rvs(self, *args, size=None, random_state=None):
986 # This method must handle size being a tuple, and it must
987 # properly broadcast *args and size. size might be
988 # an empty tuple, which means a scalar random variate is to be
989 # generated.
991 # Use basic inverse cdf algorithm for RV generation as default.
992 U = random_state.uniform(size=size)
993 Y = self._ppf(U, *args)
994 return Y
996 def _logcdf(self, x, *args):
997 with np.errstate(divide='ignore'):
998 return log(self._cdf(x, *args))
1000 def _sf(self, x, *args):
1001 return 1.0-self._cdf(x, *args)
1003 def _logsf(self, x, *args):
1004 with np.errstate(divide='ignore'):
1005 return log(self._sf(x, *args))
1007 def _ppf(self, q, *args):
1008 return self._ppfvec(q, *args)
1010 def _isf(self, q, *args):
1011 return self._ppf(1.0-q, *args) # use correct _ppf for subclasses
1013 # These are actually called, and should not be overwritten if you
1014 # want to keep error checking.
1015 def rvs(self, *args, **kwds):
1016 """Random variates of given type.
1018 Parameters
1019 ----------
1020 arg1, arg2, arg3,... : array_like
1021 The shape parameter(s) for the distribution (see docstring of the
1022 instance object for more information).
1023 loc : array_like, optional
1024 Location parameter (default=0).
1025 scale : array_like, optional
1026 Scale parameter (default=1).
1027 size : int or tuple of ints, optional
1028 Defining number of random variates (default is 1).
1029 random_state : {None, int, `numpy.random.Generator`,
1030 `numpy.random.RandomState`}, optional
1032 If `random_state` is None (or `np.random`), the
1033 `numpy.random.RandomState` singleton is used.
1034 If `random_state` is an int, a new ``RandomState`` instance is
1035 used, seeded with `random_state`.
1036 If `random_state` is already a ``Generator`` or ``RandomState``
1037 instance, that instance is used.
1039 Returns
1040 -------
1041 rvs : ndarray or scalar
1042 Random variates of given `size`.
1044 """
1045 discrete = kwds.pop('discrete', None)
1046 rndm = kwds.pop('random_state', None)
1047 args, loc, scale, size = self._parse_args_rvs(*args, **kwds)
1048 cond = logical_and(self._argcheck(*args), (scale >= 0))
1049 if not np.all(cond):
1050 message = ("Domain error in arguments. The `scale` parameter must "
1051 "be positive for all distributions, and many "
1052 "distributions have restrictions on shape parameters. "
1053 f"Please see the `scipy.stats.{self.name}` "
1054 "documentation for details.")
1055 raise ValueError(message)
1057 if np.all(scale == 0):
1058 return loc*ones(size, 'd')
1060 # extra gymnastics needed for a custom random_state
1061 if rndm is not None:
1062 random_state_saved = self._random_state
1063 random_state = check_random_state(rndm)
1064 else:
1065 random_state = self._random_state
1067 vals = self._rvs(*args, size=size, random_state=random_state)
1069 vals = vals * scale + loc
1071 # do not forget to restore the _random_state
1072 if rndm is not None:
1073 self._random_state = random_state_saved
1075 # Cast to int if discrete
1076 if discrete and not isinstance(self, rv_sample):
1077 if size == ():
1078 vals = int(vals)
1079 else:
1080 vals = vals.astype(np.int64)
1082 return vals
1084 def stats(self, *args, **kwds):
1085 """Some statistics of the given RV.
1087 Parameters
1088 ----------
1089 arg1, arg2, arg3,... : array_like
1090 The shape parameter(s) for the distribution (see docstring of the
1091 instance object for more information)
1092 loc : array_like, optional
1093 location parameter (default=0)
1094 scale : array_like, optional (continuous RVs only)
1095 scale parameter (default=1)
1096 moments : str, optional
1097 composed of letters ['mvsk'] defining which moments to compute:
1098 'm' = mean,
1099 'v' = variance,
1100 's' = (Fisher's) skew,
1101 'k' = (Fisher's) kurtosis.
1102 (default is 'mv')
1104 Returns
1105 -------
1106 stats : sequence
1107 of requested moments.
1109 """
1110 args, loc, scale, moments = self._parse_args_stats(*args, **kwds)
1111 # scale = 1 by construction for discrete RVs
1112 loc, scale = map(asarray, (loc, scale))
1113 args = tuple(map(asarray, args))
1114 cond = self._argcheck(*args) & (scale > 0) & (loc == loc)
1115 output = []
1116 default = np.full(shape(cond), fill_value=self.badvalue)
1118 # Use only entries that are valid in calculation
1119 if np.any(cond):
1120 goodargs = argsreduce(cond, *(args+(scale, loc)))
1121 scale, loc, goodargs = goodargs[-2], goodargs[-1], goodargs[:-2]
1123 if self._stats_has_moments:
1124 mu, mu2, g1, g2 = self._stats(*goodargs,
1125 **{'moments': moments})
1126 else:
1127 mu, mu2, g1, g2 = self._stats(*goodargs)
1129 if 'm' in moments:
1130 if mu is None:
1131 mu = self._munp(1, *goodargs)
1132 out0 = default.copy()
1133 place(out0, cond, mu * scale + loc)
1134 output.append(out0)
1136 if 'v' in moments:
1137 if mu2 is None:
1138 mu2p = self._munp(2, *goodargs)
1139 if mu is None:
1140 mu = self._munp(1, *goodargs)
1141 # if mean is inf then var is also inf
1142 with np.errstate(invalid='ignore'):
1143 mu2 = np.where(~np.isinf(mu), mu2p - mu**2, np.inf)
1144 out0 = default.copy()
1145 place(out0, cond, mu2 * scale * scale)
1146 output.append(out0)
1148 if 's' in moments:
1149 if g1 is None:
1150 mu3p = self._munp(3, *goodargs)
1151 if mu is None:
1152 mu = self._munp(1, *goodargs)
1153 if mu2 is None:
1154 mu2p = self._munp(2, *goodargs)
1155 mu2 = mu2p - mu * mu
1156 with np.errstate(invalid='ignore'):
1157 mu3 = (-mu*mu - 3*mu2)*mu + mu3p
1158 g1 = mu3 / np.power(mu2, 1.5)
1159 out0 = default.copy()
1160 place(out0, cond, g1)
1161 output.append(out0)
1163 if 'k' in moments:
1164 if g2 is None:
1165 mu4p = self._munp(4, *goodargs)
1166 if mu is None:
1167 mu = self._munp(1, *goodargs)
1168 if mu2 is None:
1169 mu2p = self._munp(2, *goodargs)
1170 mu2 = mu2p - mu * mu
1171 if g1 is None:
1172 mu3 = None
1173 else:
1174 # (mu2**1.5) breaks down for nan and inf
1175 mu3 = g1 * np.power(mu2, 1.5)
1176 if mu3 is None:
1177 mu3p = self._munp(3, *goodargs)
1178 with np.errstate(invalid='ignore'):
1179 mu3 = (-mu * mu - 3 * mu2) * mu + mu3p
1180 with np.errstate(invalid='ignore'):
1181 mu4 = ((-mu**2 - 6*mu2) * mu - 4*mu3)*mu + mu4p
1182 g2 = mu4 / mu2**2.0 - 3.0
1183 out0 = default.copy()
1184 place(out0, cond, g2)
1185 output.append(out0)
1186 else: # no valid args
1187 output = [default.copy() for _ in moments]
1189 output = [out[()] for out in output]
1190 if len(output) == 1:
1191 return output[0]
1192 else:
1193 return tuple(output)
1195 def entropy(self, *args, **kwds):
1196 """Differential entropy of the RV.
1198 Parameters
1199 ----------
1200 arg1, arg2, arg3,... : array_like
1201 The shape parameter(s) for the distribution (see docstring of the
1202 instance object for more information).
1203 loc : array_like, optional
1204 Location parameter (default=0).
1205 scale : array_like, optional (continuous distributions only).
1206 Scale parameter (default=1).
1208 Notes
1209 -----
1210 Entropy is defined base `e`:
1212 >>> import numpy as np
1213 >>> drv = rv_discrete(values=((0, 1), (0.5, 0.5)))
1214 >>> np.allclose(drv.entropy(), np.log(2.0))
1215 True
1217 """
1218 args, loc, scale = self._parse_args(*args, **kwds)
1219 # NB: for discrete distributions scale=1 by construction in _parse_args
1220 loc, scale = map(asarray, (loc, scale))
1221 args = tuple(map(asarray, args))
1222 cond0 = self._argcheck(*args) & (scale > 0) & (loc == loc)
1223 output = zeros(shape(cond0), 'd')
1224 place(output, (1-cond0), self.badvalue)
1225 goodargs = argsreduce(cond0, scale, *args)
1226 goodscale = goodargs[0]
1227 goodargs = goodargs[1:]
1228 place(output, cond0, self.vecentropy(*goodargs) + log(goodscale))
1229 return output[()]
1231 def moment(self, order, *args, **kwds):
1232 """non-central moment of distribution of specified order.
1234 Parameters
1235 ----------
1236 order : int, order >= 1
1237 Order of moment.
1238 arg1, arg2, arg3,... : float
1239 The shape parameter(s) for the distribution (see docstring of the
1240 instance object for more information).
1241 loc : array_like, optional
1242 location parameter (default=0)
1243 scale : array_like, optional
1244 scale parameter (default=1)
1246 """
1247 n = order
1248 shapes, loc, scale = self._parse_args(*args, **kwds)
1249 args = np.broadcast_arrays(*(*shapes, loc, scale))
1250 *shapes, loc, scale = args
1252 i0 = np.logical_and(self._argcheck(*shapes), scale > 0)
1253 i1 = np.logical_and(i0, loc == 0)
1254 i2 = np.logical_and(i0, loc != 0)
1256 args = argsreduce(i0, *shapes, loc, scale)
1257 *shapes, loc, scale = args
1259 if (floor(n) != n):
1260 raise ValueError("Moment must be an integer.")
1261 if (n < 0):
1262 raise ValueError("Moment must be positive.")
1263 mu, mu2, g1, g2 = None, None, None, None
1264 if (n > 0) and (n < 5):
1265 if self._stats_has_moments:
1266 mdict = {'moments': {1: 'm', 2: 'v', 3: 'vs', 4: 'mvsk'}[n]}
1267 else:
1268 mdict = {}
1269 mu, mu2, g1, g2 = self._stats(*shapes, **mdict)
1270 val = np.empty(loc.shape) # val needs to be indexed by loc
1271 val[...] = _moment_from_stats(n, mu, mu2, g1, g2, self._munp, shapes)
1273 # Convert to transformed X = L + S*Y
1274 # E[X^n] = E[(L+S*Y)^n] = L^n sum(comb(n, k)*(S/L)^k E[Y^k], k=0...n)
1275 result = zeros(i0.shape)
1276 place(result, ~i0, self.badvalue)
1278 if i1.any():
1279 res1 = scale[loc == 0]**n * val[loc == 0]
1280 place(result, i1, res1)
1282 if i2.any():
1283 mom = [mu, mu2, g1, g2]
1284 arrs = [i for i in mom if i is not None]
1285 idx = [i for i in range(4) if mom[i] is not None]
1286 if any(idx):
1287 arrs = argsreduce(loc != 0, *arrs)
1288 j = 0
1289 for i in idx:
1290 mom[i] = arrs[j]
1291 j += 1
1292 mu, mu2, g1, g2 = mom
1293 args = argsreduce(loc != 0, *shapes, loc, scale, val)
1294 *shapes, loc, scale, val = args
1296 res2 = zeros(loc.shape, dtype='d')
1297 fac = scale / loc
1298 for k in range(n):
1299 valk = _moment_from_stats(k, mu, mu2, g1, g2, self._munp,
1300 shapes)
1301 res2 += comb(n, k, exact=True)*fac**k * valk
1302 res2 += fac**n * val
1303 res2 *= loc**n
1304 place(result, i2, res2)
1306 return result[()]
1308 def median(self, *args, **kwds):
1309 """Median of the distribution.
1311 Parameters
1312 ----------
1313 arg1, arg2, arg3,... : array_like
1314 The shape parameter(s) for the distribution (see docstring of the
1315 instance object for more information)
1316 loc : array_like, optional
1317 Location parameter, Default is 0.
1318 scale : array_like, optional
1319 Scale parameter, Default is 1.
1321 Returns
1322 -------
1323 median : float
1324 The median of the distribution.
1326 See Also
1327 --------
1328 rv_discrete.ppf
1329 Inverse of the CDF
1331 """
1332 return self.ppf(0.5, *args, **kwds)
1334 def mean(self, *args, **kwds):
1335 """Mean of the distribution.
1337 Parameters
1338 ----------
1339 arg1, arg2, arg3,... : array_like
1340 The shape parameter(s) for the distribution (see docstring of the
1341 instance object for more information)
1342 loc : array_like, optional
1343 location parameter (default=0)
1344 scale : array_like, optional
1345 scale parameter (default=1)
1347 Returns
1348 -------
1349 mean : float
1350 the mean of the distribution
1352 """
1353 kwds['moments'] = 'm'
1354 res = self.stats(*args, **kwds)
1355 if isinstance(res, ndarray) and res.ndim == 0:
1356 return res[()]
1357 return res
1359 def var(self, *args, **kwds):
1360 """Variance of the distribution.
1362 Parameters
1363 ----------
1364 arg1, arg2, arg3,... : array_like
1365 The shape parameter(s) for the distribution (see docstring of the
1366 instance object for more information)
1367 loc : array_like, optional
1368 location parameter (default=0)
1369 scale : array_like, optional
1370 scale parameter (default=1)
1372 Returns
1373 -------
1374 var : float
1375 the variance of the distribution
1377 """
1378 kwds['moments'] = 'v'
1379 res = self.stats(*args, **kwds)
1380 if isinstance(res, ndarray) and res.ndim == 0:
1381 return res[()]
1382 return res
1384 def std(self, *args, **kwds):
1385 """Standard deviation of the distribution.
1387 Parameters
1388 ----------
1389 arg1, arg2, arg3,... : array_like
1390 The shape parameter(s) for the distribution (see docstring of the
1391 instance object for more information)
1392 loc : array_like, optional
1393 location parameter (default=0)
1394 scale : array_like, optional
1395 scale parameter (default=1)
1397 Returns
1398 -------
1399 std : float
1400 standard deviation of the distribution
1402 """
1403 kwds['moments'] = 'v'
1404 res = sqrt(self.stats(*args, **kwds))
1405 return res
1407 def interval(self, confidence, *args, **kwds):
1408 """Confidence interval with equal areas around the median.
1410 Parameters
1411 ----------
1412 confidence : array_like of float
1413 Probability that an rv will be drawn from the returned range.
1414 Each value should be in the range [0, 1].
1415 arg1, arg2, ... : array_like
1416 The shape parameter(s) for the distribution (see docstring of the
1417 instance object for more information).
1418 loc : array_like, optional
1419 location parameter, Default is 0.
1420 scale : array_like, optional
1421 scale parameter, Default is 1.
1423 Returns
1424 -------
1425 a, b : ndarray of float
1426 end-points of range that contain ``100 * alpha %`` of the rv's
1427 possible values.
1429 Notes
1430 -----
1431 This is implemented as ``ppf([p_tail, 1-p_tail])``, where
1432 ``ppf`` is the inverse cumulative distribution function and
1433 ``p_tail = (1-confidence)/2``. Suppose ``[c, d]`` is the support of a
1434 discrete distribution; then ``ppf([0, 1]) == (c-1, d)``. Therefore,
1435 when ``confidence=1`` and the distribution is discrete, the left end
1436 of the interval will be beyond the support of the distribution.
1437 For discrete distributions, the interval will limit the probability
1438 in each tail to be less than or equal to ``p_tail`` (usually
1439 strictly less).
1441 """
1442 alpha = confidence
1444 alpha = asarray(alpha)
1445 if np.any((alpha > 1) | (alpha < 0)):
1446 raise ValueError("alpha must be between 0 and 1 inclusive")
1447 q1 = (1.0-alpha)/2
1448 q2 = (1.0+alpha)/2
1449 a = self.ppf(q1, *args, **kwds)
1450 b = self.ppf(q2, *args, **kwds)
1451 return a, b
1453 def support(self, *args, **kwargs):
1454 """Support of the distribution.
1456 Parameters
1457 ----------
1458 arg1, arg2, ... : array_like
1459 The shape parameter(s) for the distribution (see docstring of the
1460 instance object for more information).
1461 loc : array_like, optional
1462 location parameter, Default is 0.
1463 scale : array_like, optional
1464 scale parameter, Default is 1.
1466 Returns
1467 -------
1468 a, b : array_like
1469 end-points of the distribution's support.
1471 """
1472 args, loc, scale = self._parse_args(*args, **kwargs)
1473 arrs = np.broadcast_arrays(*args, loc, scale)
1474 args, loc, scale = arrs[:-2], arrs[-2], arrs[-1]
1475 cond = self._argcheck(*args) & (scale > 0)
1476 _a, _b = self._get_support(*args)
1477 if cond.all():
1478 return _a * scale + loc, _b * scale + loc
1479 elif cond.ndim == 0:
1480 return self.badvalue, self.badvalue
1481 # promote bounds to at least float to fill in the badvalue
1482 _a, _b = np.asarray(_a).astype('d'), np.asarray(_b).astype('d')
1483 out_a, out_b = _a * scale + loc, _b * scale + loc
1484 place(out_a, 1-cond, self.badvalue)
1485 place(out_b, 1-cond, self.badvalue)
1486 return out_a, out_b
1488 def nnlf(self, theta, x):
1489 """Negative loglikelihood function.
1490 Notes
1491 -----
1492 This is ``-sum(log pdf(x, theta), axis=0)`` where `theta` are the
1493 parameters (including loc and scale).
1494 """
1495 loc, scale, args = self._unpack_loc_scale(theta)
1496 if not self._argcheck(*args) or scale <= 0:
1497 return inf
1498 x = (asarray(x)-loc) / scale
1499 n_log_scale = len(x) * log(scale)
1500 if np.any(~self._support_mask(x, *args)):
1501 return inf
1502 return self._nnlf(x, *args) + n_log_scale
1504 def _nnlf(self, x, *args):
1505 return -np.sum(self._logpxf(x, *args), axis=0)
1507 def _nlff_and_penalty(self, x, args, log_fitfun):
1508 # negative log fit function
1509 cond0 = ~self._support_mask(x, *args)
1510 n_bad = np.count_nonzero(cond0, axis=0)
1511 if n_bad > 0:
1512 x = argsreduce(~cond0, x)[0]
1513 logff = log_fitfun(x, *args)
1514 finite_logff = np.isfinite(logff)
1515 n_bad += np.sum(~finite_logff, axis=0)
1516 if n_bad > 0:
1517 penalty = n_bad * log(_XMAX) * 100
1518 return -np.sum(logff[finite_logff], axis=0) + penalty
1519 return -np.sum(logff, axis=0)
1521 def _penalized_nnlf(self, theta, x):
1522 """Penalized negative loglikelihood function.
1523 i.e., - sum (log pdf(x, theta), axis=0) + penalty
1524 where theta are the parameters (including loc and scale)
1525 """
1526 loc, scale, args = self._unpack_loc_scale(theta)
1527 if not self._argcheck(*args) or scale <= 0:
1528 return inf
1529 x = asarray((x-loc) / scale)
1530 n_log_scale = len(x) * log(scale)
1531 return self._nlff_and_penalty(x, args, self._logpxf) + n_log_scale
1533 def _penalized_nlpsf(self, theta, x):
1534 """Penalized negative log product spacing function.
1535 i.e., - sum (log (diff (cdf (x, theta))), axis=0) + penalty
1536 where theta are the parameters (including loc and scale)
1537 Follows reference [1] of scipy.stats.fit
1538 """
1539 loc, scale, args = self._unpack_loc_scale(theta)
1540 if not self._argcheck(*args) or scale <= 0:
1541 return inf
1542 x = (np.sort(x) - loc)/scale
1544 def log_psf(x, *args):
1545 x, lj = np.unique(x, return_counts=True) # fast for sorted x
1546 cdf_data = self._cdf(x, *args) if x.size else []
1547 if not (x.size and 1 - cdf_data[-1] <= 0):
1548 cdf = np.concatenate(([0], cdf_data, [1]))
1549 lj = np.concatenate((lj, [1]))
1550 else:
1551 cdf = np.concatenate(([0], cdf_data))
1552 # here we could use logcdf w/ logsumexp trick to take differences,
1553 # but in the context of the method, it seems unlikely to matter
1554 return lj * np.log(np.diff(cdf) / lj)
1556 return self._nlff_and_penalty(x, args, log_psf)
1559class _ShapeInfo:
1560 def __init__(self, name, integrality=False, domain=(-np.inf, np.inf),
1561 inclusive=(True, True)):
1562 self.name = name
1563 self.integrality = integrality
1565 domain = list(domain)
1566 if np.isfinite(domain[0]) and not inclusive[0]:
1567 domain[0] = np.nextafter(domain[0], np.inf)
1568 if np.isfinite(domain[1]) and not inclusive[1]:
1569 domain[1] = np.nextafter(domain[1], -np.inf)
1570 self.domain = domain
1573def _get_fixed_fit_value(kwds, names):
1574 """
1575 Given names such as `['f0', 'fa', 'fix_a']`, check that there is
1576 at most one non-None value in `kwds` associaed with those names.
1577 Return that value, or None if none of the names occur in `kwds`.
1578 As a side effect, all occurrences of those names in `kwds` are
1579 removed.
1580 """
1581 vals = [(name, kwds.pop(name)) for name in names if name in kwds]
1582 if len(vals) > 1:
1583 repeated = [name for name, val in vals]
1584 raise ValueError("fit method got multiple keyword arguments to "
1585 "specify the same fixed parameter: " +
1586 ', '.join(repeated))
1587 return vals[0][1] if vals else None
1590# continuous random variables: implement maybe later
1591#
1592# hf --- Hazard Function (PDF / SF)
1593# chf --- Cumulative hazard function (-log(SF))
1594# psf --- Probability sparsity function (reciprocal of the pdf) in
1595# units of percent-point-function (as a function of q).
1596# Also, the derivative of the percent-point function.
1599class rv_continuous(rv_generic):
1600 """A generic continuous random variable class meant for subclassing.
1602 `rv_continuous` is a base class to construct specific distribution classes
1603 and instances for continuous random variables. It cannot be used
1604 directly as a distribution.
1606 Parameters
1607 ----------
1608 momtype : int, optional
1609 The type of generic moment calculation to use: 0 for pdf, 1 (default)
1610 for ppf.
1611 a : float, optional
1612 Lower bound of the support of the distribution, default is minus
1613 infinity.
1614 b : float, optional
1615 Upper bound of the support of the distribution, default is plus
1616 infinity.
1617 xtol : float, optional
1618 The tolerance for fixed point calculation for generic ppf.
1619 badvalue : float, optional
1620 The value in a result arrays that indicates a value that for which
1621 some argument restriction is violated, default is np.nan.
1622 name : str, optional
1623 The name of the instance. This string is used to construct the default
1624 example for distributions.
1625 longname : str, optional
1626 This string is used as part of the first line of the docstring returned
1627 when a subclass has no docstring of its own. Note: `longname` exists
1628 for backwards compatibility, do not use for new subclasses.
1629 shapes : str, optional
1630 The shape of the distribution. For example ``"m, n"`` for a
1631 distribution that takes two integers as the two shape arguments for all
1632 its methods. If not provided, shape parameters will be inferred from
1633 the signature of the private methods, ``_pdf`` and ``_cdf`` of the
1634 instance.
1635 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
1636 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
1637 singleton is used.
1638 If `seed` is an int, a new ``RandomState`` instance is used,
1639 seeded with `seed`.
1640 If `seed` is already a ``Generator`` or ``RandomState`` instance then
1641 that instance is used.
1643 Methods
1644 -------
1645 rvs
1646 pdf
1647 logpdf
1648 cdf
1649 logcdf
1650 sf
1651 logsf
1652 ppf
1653 isf
1654 moment
1655 stats
1656 entropy
1657 expect
1658 median
1659 mean
1660 std
1661 var
1662 interval
1663 __call__
1664 fit
1665 fit_loc_scale
1666 nnlf
1667 support
1669 Notes
1670 -----
1671 Public methods of an instance of a distribution class (e.g., ``pdf``,
1672 ``cdf``) check their arguments and pass valid arguments to private,
1673 computational methods (``_pdf``, ``_cdf``). For ``pdf(x)``, ``x`` is valid
1674 if it is within the support of the distribution.
1675 Whether a shape parameter is valid is decided by an ``_argcheck`` method
1676 (which defaults to checking that its arguments are strictly positive.)
1678 **Subclassing**
1680 New random variables can be defined by subclassing the `rv_continuous` class
1681 and re-defining at least the ``_pdf`` or the ``_cdf`` method (normalized
1682 to location 0 and scale 1).
1684 If positive argument checking is not correct for your RV
1685 then you will also need to re-define the ``_argcheck`` method.
1687 For most of the scipy.stats distributions, the support interval doesn't
1688 depend on the shape parameters. ``x`` being in the support interval is
1689 equivalent to ``self.a <= x <= self.b``. If either of the endpoints of
1690 the support do depend on the shape parameters, then
1691 i) the distribution must implement the ``_get_support`` method; and
1692 ii) those dependent endpoints must be omitted from the distribution's
1693 call to the ``rv_continuous`` initializer.
1695 Correct, but potentially slow defaults exist for the remaining
1696 methods but for speed and/or accuracy you can over-ride::
1698 _logpdf, _cdf, _logcdf, _ppf, _rvs, _isf, _sf, _logsf
1700 The default method ``_rvs`` relies on the inverse of the cdf, ``_ppf``,
1701 applied to a uniform random variate. In order to generate random variates
1702 efficiently, either the default ``_ppf`` needs to be overwritten (e.g.
1703 if the inverse cdf can expressed in an explicit form) or a sampling
1704 method needs to be implemented in a custom ``_rvs`` method.
1706 If possible, you should override ``_isf``, ``_sf`` or ``_logsf``.
1707 The main reason would be to improve numerical accuracy: for example,
1708 the survival function ``_sf`` is computed as ``1 - _cdf`` which can
1709 result in loss of precision if ``_cdf(x)`` is close to one.
1711 **Methods that can be overwritten by subclasses**
1712 ::
1714 _rvs
1715 _pdf
1716 _cdf
1717 _sf
1718 _ppf
1719 _isf
1720 _stats
1721 _munp
1722 _entropy
1723 _argcheck
1724 _get_support
1726 There are additional (internal and private) generic methods that can
1727 be useful for cross-checking and for debugging, but might work in all
1728 cases when directly called.
1730 A note on ``shapes``: subclasses need not specify them explicitly. In this
1731 case, `shapes` will be automatically deduced from the signatures of the
1732 overridden methods (`pdf`, `cdf` etc).
1733 If, for some reason, you prefer to avoid relying on introspection, you can
1734 specify ``shapes`` explicitly as an argument to the instance constructor.
1737 **Frozen Distributions**
1739 Normally, you must provide shape parameters (and, optionally, location and
1740 scale parameters to each call of a method of a distribution.
1742 Alternatively, the object may be called (as a function) to fix the shape,
1743 location, and scale parameters returning a "frozen" continuous RV object:
1745 rv = generic(<shape(s)>, loc=0, scale=1)
1746 `rv_frozen` object with the same methods but holding the given shape,
1747 location, and scale fixed
1749 **Statistics**
1751 Statistics are computed using numerical integration by default.
1752 For speed you can redefine this using ``_stats``:
1754 - take shape parameters and return mu, mu2, g1, g2
1755 - If you can't compute one of these, return it as None
1756 - Can also be defined with a keyword argument ``moments``, which is a
1757 string composed of "m", "v", "s", and/or "k".
1758 Only the components appearing in string should be computed and
1759 returned in the order "m", "v", "s", or "k" with missing values
1760 returned as None.
1762 Alternatively, you can override ``_munp``, which takes ``n`` and shape
1763 parameters and returns the n-th non-central moment of the distribution.
1765 **Deepcopying / Pickling**
1767 If a distribution or frozen distribution is deepcopied (pickled/unpickled,
1768 etc.), any underlying random number generator is deepcopied with it. An
1769 implication is that if a distribution relies on the singleton RandomState
1770 before copying, it will rely on a copy of that random state after copying,
1771 and ``np.random.seed`` will no longer control the state.
1773 Examples
1774 --------
1775 To create a new Gaussian distribution, we would do the following:
1777 >>> from scipy.stats import rv_continuous
1778 >>> class gaussian_gen(rv_continuous):
1779 ... "Gaussian distribution"
1780 ... def _pdf(self, x):
1781 ... return np.exp(-x**2 / 2.) / np.sqrt(2.0 * np.pi)
1782 >>> gaussian = gaussian_gen(name='gaussian')
1784 ``scipy.stats`` distributions are *instances*, so here we subclass
1785 `rv_continuous` and create an instance. With this, we now have
1786 a fully functional distribution with all relevant methods automagically
1787 generated by the framework.
1789 Note that above we defined a standard normal distribution, with zero mean
1790 and unit variance. Shifting and scaling of the distribution can be done
1791 by using ``loc`` and ``scale`` parameters: ``gaussian.pdf(x, loc, scale)``
1792 essentially computes ``y = (x - loc) / scale`` and
1793 ``gaussian._pdf(y) / scale``.
1795 """
1797 def __init__(self, momtype=1, a=None, b=None, xtol=1e-14,
1798 badvalue=None, name=None, longname=None,
1799 shapes=None, seed=None):
1801 super().__init__(seed)
1803 # save the ctor parameters, cf generic freeze
1804 self._ctor_param = dict(
1805 momtype=momtype, a=a, b=b, xtol=xtol,
1806 badvalue=badvalue, name=name, longname=longname,
1807 shapes=shapes, seed=seed)
1809 if badvalue is None:
1810 badvalue = nan
1811 if name is None:
1812 name = 'Distribution'
1813 self.badvalue = badvalue
1814 self.name = name
1815 self.a = a
1816 self.b = b
1817 if a is None:
1818 self.a = -inf
1819 if b is None:
1820 self.b = inf
1821 self.xtol = xtol
1822 self.moment_type = momtype
1823 self.shapes = shapes
1825 self._construct_argparser(meths_to_inspect=[self._pdf, self._cdf],
1826 locscale_in='loc=0, scale=1',
1827 locscale_out='loc, scale')
1828 self._attach_methods()
1830 if longname is None:
1831 if name[0] in ['aeiouAEIOU']:
1832 hstr = "An "
1833 else:
1834 hstr = "A "
1835 longname = hstr + name
1837 if sys.flags.optimize < 2:
1838 # Skip adding docstrings if interpreter is run with -OO
1839 if self.__doc__ is None:
1840 self._construct_default_doc(longname=longname,
1841 docdict=docdict,
1842 discrete='continuous')
1843 else:
1844 dct = dict(distcont)
1845 self._construct_doc(docdict, dct.get(self.name))
1847 def __getstate__(self):
1848 dct = self.__dict__.copy()
1850 # these methods will be remade in __setstate__
1851 # _random_state attribute is taken care of by rv_generic
1852 attrs = ["_parse_args", "_parse_args_stats", "_parse_args_rvs",
1853 "_cdfvec", "_ppfvec", "vecentropy", "generic_moment"]
1854 [dct.pop(attr, None) for attr in attrs]
1855 return dct
1857 def _attach_methods(self):
1858 """
1859 Attaches dynamically created methods to the rv_continuous instance.
1860 """
1861 # _attach_methods is responsible for calling _attach_argparser_methods
1862 self._attach_argparser_methods()
1864 # nin correction
1865 self._ppfvec = vectorize(self._ppf_single, otypes='d')
1866 self._ppfvec.nin = self.numargs + 1
1867 self.vecentropy = vectorize(self._entropy, otypes='d')
1868 self._cdfvec = vectorize(self._cdf_single, otypes='d')
1869 self._cdfvec.nin = self.numargs + 1
1871 if self.moment_type == 0:
1872 self.generic_moment = vectorize(self._mom0_sc, otypes='d')
1873 else:
1874 self.generic_moment = vectorize(self._mom1_sc, otypes='d')
1875 # Because of the *args argument of _mom0_sc, vectorize cannot count the
1876 # number of arguments correctly.
1877 self.generic_moment.nin = self.numargs + 1
1879 def _updated_ctor_param(self):
1880 """Return the current version of _ctor_param, possibly updated by user.
1882 Used by freezing.
1883 Keep this in sync with the signature of __init__.
1884 """
1885 dct = self._ctor_param.copy()
1886 dct['a'] = self.a
1887 dct['b'] = self.b
1888 dct['xtol'] = self.xtol
1889 dct['badvalue'] = self.badvalue
1890 dct['name'] = self.name
1891 dct['shapes'] = self.shapes
1892 return dct
1894 def _ppf_to_solve(self, x, q, *args):
1895 return self.cdf(*(x, )+args)-q
1897 def _ppf_single(self, q, *args):
1898 factor = 10.
1899 left, right = self._get_support(*args)
1901 if np.isinf(left):
1902 left = min(-factor, right)
1903 while self._ppf_to_solve(left, q, *args) > 0.:
1904 left, right = left * factor, left
1905 # left is now such that cdf(left) <= q
1906 # if right has changed, then cdf(right) > q
1908 if np.isinf(right):
1909 right = max(factor, left)
1910 while self._ppf_to_solve(right, q, *args) < 0.:
1911 left, right = right, right * factor
1912 # right is now such that cdf(right) >= q
1914 return optimize.brentq(self._ppf_to_solve,
1915 left, right, args=(q,)+args, xtol=self.xtol)
1917 # moment from definition
1918 def _mom_integ0(self, x, m, *args):
1919 return x**m * self.pdf(x, *args)
1921 def _mom0_sc(self, m, *args):
1922 _a, _b = self._get_support(*args)
1923 return integrate.quad(self._mom_integ0, _a, _b,
1924 args=(m,)+args)[0]
1926 # moment calculated using ppf
1927 def _mom_integ1(self, q, m, *args):
1928 return (self.ppf(q, *args))**m
1930 def _mom1_sc(self, m, *args):
1931 return integrate.quad(self._mom_integ1, 0, 1, args=(m,)+args)[0]
1933 def _pdf(self, x, *args):
1934 return _derivative(self._cdf, x, dx=1e-5, args=args, order=5)
1936 # Could also define any of these
1937 def _logpdf(self, x, *args):
1938 p = self._pdf(x, *args)
1939 with np.errstate(divide='ignore'):
1940 return log(p)
1942 def _logpxf(self, x, *args):
1943 # continuous distributions have PDF, discrete have PMF, but sometimes
1944 # the distinction doesn't matter. This lets us use `_logpxf` for both
1945 # discrete and continuous distributions.
1946 return self._logpdf(x, *args)
1948 def _cdf_single(self, x, *args):
1949 _a, _b = self._get_support(*args)
1950 return integrate.quad(self._pdf, _a, x, args=args)[0]
1952 def _cdf(self, x, *args):
1953 return self._cdfvec(x, *args)
1955 # generic _argcheck, _logcdf, _sf, _logsf, _ppf, _isf, _rvs are defined
1956 # in rv_generic
1958 def pdf(self, x, *args, **kwds):
1959 """Probability density function at x of the given RV.
1961 Parameters
1962 ----------
1963 x : array_like
1964 quantiles
1965 arg1, arg2, arg3,... : array_like
1966 The shape parameter(s) for the distribution (see docstring of the
1967 instance object for more information)
1968 loc : array_like, optional
1969 location parameter (default=0)
1970 scale : array_like, optional
1971 scale parameter (default=1)
1973 Returns
1974 -------
1975 pdf : ndarray
1976 Probability density function evaluated at x
1978 """
1979 args, loc, scale = self._parse_args(*args, **kwds)
1980 x, loc, scale = map(asarray, (x, loc, scale))
1981 args = tuple(map(asarray, args))
1982 dtyp = np.promote_types(x.dtype, np.float64)
1983 x = np.asarray((x - loc)/scale, dtype=dtyp)
1984 cond0 = self._argcheck(*args) & (scale > 0)
1985 cond1 = self._support_mask(x, *args) & (scale > 0)
1986 cond = cond0 & cond1
1987 output = zeros(shape(cond), dtyp)
1988 putmask(output, (1-cond0)+np.isnan(x), self.badvalue)
1989 if np.any(cond):
1990 goodargs = argsreduce(cond, *((x,)+args+(scale,)))
1991 scale, goodargs = goodargs[-1], goodargs[:-1]
1992 place(output, cond, self._pdf(*goodargs) / scale)
1993 if output.ndim == 0:
1994 return output[()]
1995 return output
1997 def logpdf(self, x, *args, **kwds):
1998 """Log of the probability density function at x of the given RV.
2000 This uses a more numerically accurate calculation if available.
2002 Parameters
2003 ----------
2004 x : array_like
2005 quantiles
2006 arg1, arg2, arg3,... : array_like
2007 The shape parameter(s) for the distribution (see docstring of the
2008 instance object for more information)
2009 loc : array_like, optional
2010 location parameter (default=0)
2011 scale : array_like, optional
2012 scale parameter (default=1)
2014 Returns
2015 -------
2016 logpdf : array_like
2017 Log of the probability density function evaluated at x
2019 """
2020 args, loc, scale = self._parse_args(*args, **kwds)
2021 x, loc, scale = map(asarray, (x, loc, scale))
2022 args = tuple(map(asarray, args))
2023 dtyp = np.promote_types(x.dtype, np.float64)
2024 x = np.asarray((x - loc)/scale, dtype=dtyp)
2025 cond0 = self._argcheck(*args) & (scale > 0)
2026 cond1 = self._support_mask(x, *args) & (scale > 0)
2027 cond = cond0 & cond1
2028 output = empty(shape(cond), dtyp)
2029 output.fill(NINF)
2030 putmask(output, (1-cond0)+np.isnan(x), self.badvalue)
2031 if np.any(cond):
2032 goodargs = argsreduce(cond, *((x,)+args+(scale,)))
2033 scale, goodargs = goodargs[-1], goodargs[:-1]
2034 place(output, cond, self._logpdf(*goodargs) - log(scale))
2035 if output.ndim == 0:
2036 return output[()]
2037 return output
2039 def cdf(self, x, *args, **kwds):
2040 """
2041 Cumulative distribution function of the given RV.
2043 Parameters
2044 ----------
2045 x : array_like
2046 quantiles
2047 arg1, arg2, arg3,... : array_like
2048 The shape parameter(s) for the distribution (see docstring of the
2049 instance object for more information)
2050 loc : array_like, optional
2051 location parameter (default=0)
2052 scale : array_like, optional
2053 scale parameter (default=1)
2055 Returns
2056 -------
2057 cdf : ndarray
2058 Cumulative distribution function evaluated at `x`
2060 """
2061 args, loc, scale = self._parse_args(*args, **kwds)
2062 x, loc, scale = map(asarray, (x, loc, scale))
2063 args = tuple(map(asarray, args))
2064 _a, _b = self._get_support(*args)
2065 dtyp = np.promote_types(x.dtype, np.float64)
2066 x = np.asarray((x - loc)/scale, dtype=dtyp)
2067 cond0 = self._argcheck(*args) & (scale > 0)
2068 cond1 = self._open_support_mask(x, *args) & (scale > 0)
2069 cond2 = (x >= np.asarray(_b)) & cond0
2070 cond = cond0 & cond1
2071 output = zeros(shape(cond), dtyp)
2072 place(output, (1-cond0)+np.isnan(x), self.badvalue)
2073 place(output, cond2, 1.0)
2074 if np.any(cond): # call only if at least 1 entry
2075 goodargs = argsreduce(cond, *((x,)+args))
2076 place(output, cond, self._cdf(*goodargs))
2077 if output.ndim == 0:
2078 return output[()]
2079 return output
2081 def logcdf(self, x, *args, **kwds):
2082 """Log of the cumulative distribution function at x of the given RV.
2084 Parameters
2085 ----------
2086 x : array_like
2087 quantiles
2088 arg1, arg2, arg3,... : array_like
2089 The shape parameter(s) for the distribution (see docstring of the
2090 instance object for more information)
2091 loc : array_like, optional
2092 location parameter (default=0)
2093 scale : array_like, optional
2094 scale parameter (default=1)
2096 Returns
2097 -------
2098 logcdf : array_like
2099 Log of the cumulative distribution function evaluated at x
2101 """
2102 args, loc, scale = self._parse_args(*args, **kwds)
2103 x, loc, scale = map(asarray, (x, loc, scale))
2104 args = tuple(map(asarray, args))
2105 _a, _b = self._get_support(*args)
2106 dtyp = np.promote_types(x.dtype, np.float64)
2107 x = np.asarray((x - loc)/scale, dtype=dtyp)
2108 cond0 = self._argcheck(*args) & (scale > 0)
2109 cond1 = self._open_support_mask(x, *args) & (scale > 0)
2110 cond2 = (x >= _b) & cond0
2111 cond = cond0 & cond1
2112 output = empty(shape(cond), dtyp)
2113 output.fill(NINF)
2114 place(output, (1-cond0)*(cond1 == cond1)+np.isnan(x), self.badvalue)
2115 place(output, cond2, 0.0)
2116 if np.any(cond): # call only if at least 1 entry
2117 goodargs = argsreduce(cond, *((x,)+args))
2118 place(output, cond, self._logcdf(*goodargs))
2119 if output.ndim == 0:
2120 return output[()]
2121 return output
2123 def sf(self, x, *args, **kwds):
2124 """Survival function (1 - `cdf`) at x of the given RV.
2126 Parameters
2127 ----------
2128 x : array_like
2129 quantiles
2130 arg1, arg2, arg3,... : array_like
2131 The shape parameter(s) for the distribution (see docstring of the
2132 instance object for more information)
2133 loc : array_like, optional
2134 location parameter (default=0)
2135 scale : array_like, optional
2136 scale parameter (default=1)
2138 Returns
2139 -------
2140 sf : array_like
2141 Survival function evaluated at x
2143 """
2144 args, loc, scale = self._parse_args(*args, **kwds)
2145 x, loc, scale = map(asarray, (x, loc, scale))
2146 args = tuple(map(asarray, args))
2147 _a, _b = self._get_support(*args)
2148 dtyp = np.promote_types(x.dtype, np.float64)
2149 x = np.asarray((x - loc)/scale, dtype=dtyp)
2150 cond0 = self._argcheck(*args) & (scale > 0)
2151 cond1 = self._open_support_mask(x, *args) & (scale > 0)
2152 cond2 = cond0 & (x <= _a)
2153 cond = cond0 & cond1
2154 output = zeros(shape(cond), dtyp)
2155 place(output, (1-cond0)+np.isnan(x), self.badvalue)
2156 place(output, cond2, 1.0)
2157 if np.any(cond):
2158 goodargs = argsreduce(cond, *((x,)+args))
2159 place(output, cond, self._sf(*goodargs))
2160 if output.ndim == 0:
2161 return output[()]
2162 return output
2164 def logsf(self, x, *args, **kwds):
2165 """Log of the survival function of the given RV.
2167 Returns the log of the "survival function," defined as (1 - `cdf`),
2168 evaluated at `x`.
2170 Parameters
2171 ----------
2172 x : array_like
2173 quantiles
2174 arg1, arg2, arg3,... : array_like
2175 The shape parameter(s) for the distribution (see docstring of the
2176 instance object for more information)
2177 loc : array_like, optional
2178 location parameter (default=0)
2179 scale : array_like, optional
2180 scale parameter (default=1)
2182 Returns
2183 -------
2184 logsf : ndarray
2185 Log of the survival function evaluated at `x`.
2187 """
2188 args, loc, scale = self._parse_args(*args, **kwds)
2189 x, loc, scale = map(asarray, (x, loc, scale))
2190 args = tuple(map(asarray, args))
2191 _a, _b = self._get_support(*args)
2192 dtyp = np.promote_types(x.dtype, np.float64)
2193 x = np.asarray((x - loc)/scale, dtype=dtyp)
2194 cond0 = self._argcheck(*args) & (scale > 0)
2195 cond1 = self._open_support_mask(x, *args) & (scale > 0)
2196 cond2 = cond0 & (x <= _a)
2197 cond = cond0 & cond1
2198 output = empty(shape(cond), dtyp)
2199 output.fill(NINF)
2200 place(output, (1-cond0)+np.isnan(x), self.badvalue)
2201 place(output, cond2, 0.0)
2202 if np.any(cond):
2203 goodargs = argsreduce(cond, *((x,)+args))
2204 place(output, cond, self._logsf(*goodargs))
2205 if output.ndim == 0:
2206 return output[()]
2207 return output
2209 def ppf(self, q, *args, **kwds):
2210 """Percent point function (inverse of `cdf`) at q of the given RV.
2212 Parameters
2213 ----------
2214 q : array_like
2215 lower tail probability
2216 arg1, arg2, arg3,... : array_like
2217 The shape parameter(s) for the distribution (see docstring of the
2218 instance object for more information)
2219 loc : array_like, optional
2220 location parameter (default=0)
2221 scale : array_like, optional
2222 scale parameter (default=1)
2224 Returns
2225 -------
2226 x : array_like
2227 quantile corresponding to the lower tail probability q.
2229 """
2230 args, loc, scale = self._parse_args(*args, **kwds)
2231 q, loc, scale = map(asarray, (q, loc, scale))
2232 args = tuple(map(asarray, args))
2233 _a, _b = self._get_support(*args)
2234 cond0 = self._argcheck(*args) & (scale > 0) & (loc == loc)
2235 cond1 = (0 < q) & (q < 1)
2236 cond2 = cond0 & (q == 0)
2237 cond3 = cond0 & (q == 1)
2238 cond = cond0 & cond1
2239 output = np.full(shape(cond), fill_value=self.badvalue)
2241 lower_bound = _a * scale + loc
2242 upper_bound = _b * scale + loc
2243 place(output, cond2, argsreduce(cond2, lower_bound)[0])
2244 place(output, cond3, argsreduce(cond3, upper_bound)[0])
2246 if np.any(cond): # call only if at least 1 entry
2247 goodargs = argsreduce(cond, *((q,)+args+(scale, loc)))
2248 scale, loc, goodargs = goodargs[-2], goodargs[-1], goodargs[:-2]
2249 place(output, cond, self._ppf(*goodargs) * scale + loc)
2250 if output.ndim == 0:
2251 return output[()]
2252 return output
2254 def isf(self, q, *args, **kwds):
2255 """Inverse survival function (inverse of `sf`) at q of the given RV.
2257 Parameters
2258 ----------
2259 q : array_like
2260 upper tail probability
2261 arg1, arg2, arg3,... : array_like
2262 The shape parameter(s) for the distribution (see docstring of the
2263 instance object for more information)
2264 loc : array_like, optional
2265 location parameter (default=0)
2266 scale : array_like, optional
2267 scale parameter (default=1)
2269 Returns
2270 -------
2271 x : ndarray or scalar
2272 Quantile corresponding to the upper tail probability q.
2274 """
2275 args, loc, scale = self._parse_args(*args, **kwds)
2276 q, loc, scale = map(asarray, (q, loc, scale))
2277 args = tuple(map(asarray, args))
2278 _a, _b = self._get_support(*args)
2279 cond0 = self._argcheck(*args) & (scale > 0) & (loc == loc)
2280 cond1 = (0 < q) & (q < 1)
2281 cond2 = cond0 & (q == 1)
2282 cond3 = cond0 & (q == 0)
2283 cond = cond0 & cond1
2284 output = np.full(shape(cond), fill_value=self.badvalue)
2286 lower_bound = _a * scale + loc
2287 upper_bound = _b * scale + loc
2288 place(output, cond2, argsreduce(cond2, lower_bound)[0])
2289 place(output, cond3, argsreduce(cond3, upper_bound)[0])
2291 if np.any(cond):
2292 goodargs = argsreduce(cond, *((q,)+args+(scale, loc)))
2293 scale, loc, goodargs = goodargs[-2], goodargs[-1], goodargs[:-2]
2294 place(output, cond, self._isf(*goodargs) * scale + loc)
2295 if output.ndim == 0:
2296 return output[()]
2297 return output
2299 def _unpack_loc_scale(self, theta):
2300 try:
2301 loc = theta[-2]
2302 scale = theta[-1]
2303 args = tuple(theta[:-2])
2304 except IndexError as e:
2305 raise ValueError("Not enough input arguments.") from e
2306 return loc, scale, args
2308 def _nnlf_and_penalty(self, x, args):
2309 """
2310 Compute the penalized negative log-likelihood for the
2311 "standardized" data (i.e. already shifted by loc and
2312 scaled by scale) for the shape parameters in `args`.
2314 `x` can be a 1D numpy array or a CensoredData instance.
2315 """
2316 if isinstance(x, CensoredData):
2317 # Filter out the data that is not in the support.
2318 xs = x._supported(*self._get_support(*args))
2319 n_bad = len(x) - len(xs)
2320 i1, i2 = xs._interval.T
2321 terms = [
2322 # logpdf of the noncensored data.
2323 self._logpdf(xs._uncensored, *args),
2324 # logcdf of the left-censored data.
2325 self._logcdf(xs._left, *args),
2326 # logsf of the right-censored data.
2327 self._logsf(xs._right, *args),
2328 # log of probability of the interval-censored data.
2329 np.log(self._delta_cdf(i1, i2, *args)),
2330 ]
2331 else:
2332 cond0 = ~self._support_mask(x, *args)
2333 n_bad = np.count_nonzero(cond0)
2334 if n_bad > 0:
2335 x = argsreduce(~cond0, x)[0]
2336 terms = [self._logpdf(x, *args)]
2338 totals, bad_counts = zip(*[_sum_finite(term) for term in terms])
2339 total = sum(totals)
2340 n_bad += sum(bad_counts)
2342 return -total + n_bad * _LOGXMAX * 100
2344 def _penalized_nnlf(self, theta, x):
2345 """Penalized negative loglikelihood function.
2347 i.e., - sum (log pdf(x, theta), axis=0) + penalty
2348 where theta are the parameters (including loc and scale)
2349 """
2350 loc, scale, args = self._unpack_loc_scale(theta)
2351 if not self._argcheck(*args) or scale <= 0:
2352 return inf
2353 if isinstance(x, CensoredData):
2354 x = (x - loc) / scale
2355 n_log_scale = (len(x) - x.num_censored()) * log(scale)
2356 else:
2357 x = (x - loc) / scale
2358 n_log_scale = len(x) * log(scale)
2360 return self._nnlf_and_penalty(x, args) + n_log_scale
2362 def _fitstart(self, data, args=None):
2363 """Starting point for fit (shape arguments + loc + scale)."""
2364 if args is None:
2365 args = (1.0,)*self.numargs
2366 loc, scale = self._fit_loc_scale_support(data, *args)
2367 return args + (loc, scale)
2369 def _reduce_func(self, args, kwds, data=None):
2370 """
2371 Return the (possibly reduced) function to optimize in order to find MLE
2372 estimates for the .fit method.
2373 """
2374 # Convert fixed shape parameters to the standard numeric form: e.g. for
2375 # stats.beta, shapes='a, b'. To fix `a`, the caller can give a value
2376 # for `f0`, `fa` or 'fix_a'. The following converts the latter two
2377 # into the first (numeric) form.
2378 shapes = []
2379 if self.shapes:
2380 shapes = self.shapes.replace(',', ' ').split()
2381 for j, s in enumerate(shapes):
2382 key = 'f' + str(j)
2383 names = [key, 'f' + s, 'fix_' + s]
2384 val = _get_fixed_fit_value(kwds, names)
2385 if val is not None:
2386 kwds[key] = val
2388 args = list(args)
2389 Nargs = len(args)
2390 fixedn = []
2391 names = ['f%d' % n for n in range(Nargs - 2)] + ['floc', 'fscale']
2392 x0 = []
2393 for n, key in enumerate(names):
2394 if key in kwds:
2395 fixedn.append(n)
2396 args[n] = kwds.pop(key)
2397 else:
2398 x0.append(args[n])
2400 methods = {"mle", "mm"}
2401 method = kwds.pop('method', "mle").lower()
2402 if method == "mm":
2403 n_params = len(shapes) + 2 - len(fixedn)
2404 exponents = (np.arange(1, n_params+1))[:, np.newaxis]
2405 data_moments = np.sum(data[None, :]**exponents/len(data), axis=1)
2407 def objective(theta, x):
2408 return self._moment_error(theta, x, data_moments)
2410 elif method == "mle":
2411 objective = self._penalized_nnlf
2412 else:
2413 raise ValueError("Method '{}' not available; must be one of {}"
2414 .format(method, methods))
2416 if len(fixedn) == 0:
2417 func = objective
2418 restore = None
2419 else:
2420 if len(fixedn) == Nargs:
2421 raise ValueError(
2422 "All parameters fixed. There is nothing to optimize.")
2424 def restore(args, theta):
2425 # Replace with theta for all numbers not in fixedn
2426 # This allows the non-fixed values to vary, but
2427 # we still call self.nnlf with all parameters.
2428 i = 0
2429 for n in range(Nargs):
2430 if n not in fixedn:
2431 args[n] = theta[i]
2432 i += 1
2433 return args
2435 def func(theta, x):
2436 newtheta = restore(args[:], theta)
2437 return objective(newtheta, x)
2439 return x0, func, restore, args
2441 def _moment_error(self, theta, x, data_moments):
2442 loc, scale, args = self._unpack_loc_scale(theta)
2443 if not self._argcheck(*args) or scale <= 0:
2444 return inf
2446 dist_moments = np.array([self.moment(i+1, *args, loc=loc, scale=scale)
2447 for i in range(len(data_moments))])
2448 if np.any(np.isnan(dist_moments)):
2449 raise ValueError("Method of moments encountered a non-finite "
2450 "distribution moment and cannot continue. "
2451 "Consider trying method='MLE'.")
2453 return (((data_moments - dist_moments) /
2454 np.maximum(np.abs(data_moments), 1e-8))**2).sum()
2456 def fit(self, data, *args, **kwds):
2457 """
2458 Return estimates of shape (if applicable), location, and scale
2459 parameters from data. The default estimation method is Maximum
2460 Likelihood Estimation (MLE), but Method of Moments (MM)
2461 is also available.
2463 Starting estimates for the fit are given by input arguments;
2464 for any arguments not provided with starting estimates,
2465 ``self._fitstart(data)`` is called to generate such.
2467 One can hold some parameters fixed to specific values by passing in
2468 keyword arguments ``f0``, ``f1``, ..., ``fn`` (for shape parameters)
2469 and ``floc`` and ``fscale`` (for location and scale parameters,
2470 respectively).
2472 Parameters
2473 ----------
2474 data : array_like or `CensoredData` instance
2475 Data to use in estimating the distribution parameters.
2476 arg1, arg2, arg3,... : floats, optional
2477 Starting value(s) for any shape-characterizing arguments (those not
2478 provided will be determined by a call to ``_fitstart(data)``).
2479 No default value.
2480 **kwds : floats, optional
2481 - `loc`: initial guess of the distribution's location parameter.
2482 - `scale`: initial guess of the distribution's scale parameter.
2484 Special keyword arguments are recognized as holding certain
2485 parameters fixed:
2487 - f0...fn : hold respective shape parameters fixed.
2488 Alternatively, shape parameters to fix can be specified by name.
2489 For example, if ``self.shapes == "a, b"``, ``fa`` and ``fix_a``
2490 are equivalent to ``f0``, and ``fb`` and ``fix_b`` are
2491 equivalent to ``f1``.
2493 - floc : hold location parameter fixed to specified value.
2495 - fscale : hold scale parameter fixed to specified value.
2497 - optimizer : The optimizer to use. The optimizer must take
2498 ``func`` and starting position as the first two arguments,
2499 plus ``args`` (for extra arguments to pass to the
2500 function to be optimized) and ``disp=0`` to suppress
2501 output as keyword arguments.
2503 - method : The method to use. The default is "MLE" (Maximum
2504 Likelihood Estimate); "MM" (Method of Moments)
2505 is also available.
2507 Raises
2508 ------
2509 TypeError, ValueError
2510 If an input is invalid
2511 `~scipy.stats.FitError`
2512 If fitting fails or the fit produced would be invalid
2514 Returns
2515 -------
2516 parameter_tuple : tuple of floats
2517 Estimates for any shape parameters (if applicable), followed by
2518 those for location and scale. For most random variables, shape
2519 statistics will be returned, but there are exceptions (e.g.
2520 ``norm``).
2522 Notes
2523 -----
2524 With ``method="MLE"`` (default), the fit is computed by minimizing
2525 the negative log-likelihood function. A large, finite penalty
2526 (rather than infinite negative log-likelihood) is applied for
2527 observations beyond the support of the distribution.
2529 With ``method="MM"``, the fit is computed by minimizing the L2 norm
2530 of the relative errors between the first *k* raw (about zero) data
2531 moments and the corresponding distribution moments, where *k* is the
2532 number of non-fixed parameters.
2533 More precisely, the objective function is::
2535 (((data_moments - dist_moments)
2536 / np.maximum(np.abs(data_moments), 1e-8))**2).sum()
2538 where the constant ``1e-8`` avoids division by zero in case of
2539 vanishing data moments. Typically, this error norm can be reduced to
2540 zero.
2541 Note that the standard method of moments can produce parameters for
2542 which some data are outside the support of the fitted distribution;
2543 this implementation does nothing to prevent this.
2545 For either method,
2546 the returned answer is not guaranteed to be globally optimal; it
2547 may only be locally optimal, or the optimization may fail altogether.
2548 If the data contain any of ``np.nan``, ``np.inf``, or ``-np.inf``,
2549 the `fit` method will raise a ``RuntimeError``.
2551 Examples
2552 --------
2554 Generate some data to fit: draw random variates from the `beta`
2555 distribution
2557 >>> from scipy.stats import beta
2558 >>> a, b = 1., 2.
2559 >>> x = beta.rvs(a, b, size=1000)
2561 Now we can fit all four parameters (``a``, ``b``, ``loc`` and
2562 ``scale``):
2564 >>> a1, b1, loc1, scale1 = beta.fit(x)
2566 We can also use some prior knowledge about the dataset: let's keep
2567 ``loc`` and ``scale`` fixed:
2569 >>> a1, b1, loc1, scale1 = beta.fit(x, floc=0, fscale=1)
2570 >>> loc1, scale1
2571 (0, 1)
2573 We can also keep shape parameters fixed by using ``f``-keywords. To
2574 keep the zero-th shape parameter ``a`` equal 1, use ``f0=1`` or,
2575 equivalently, ``fa=1``:
2577 >>> a1, b1, loc1, scale1 = beta.fit(x, fa=1, floc=0, fscale=1)
2578 >>> a1
2579 1
2581 Not all distributions return estimates for the shape parameters.
2582 ``norm`` for example just returns estimates for location and scale:
2584 >>> from scipy.stats import norm
2585 >>> x = norm.rvs(a, b, size=1000, random_state=123)
2586 >>> loc1, scale1 = norm.fit(x)
2587 >>> loc1, scale1
2588 (0.92087172783841631, 2.0015750750324668)
2589 """
2590 method = kwds.get('method', "mle").lower()
2592 censored = isinstance(data, CensoredData)
2593 if censored:
2594 if method != 'mle':
2595 raise ValueError('For censored data, the method must'
2596 ' be "MLE".')
2597 if data.num_censored() == 0:
2598 # There are no censored values in data, so replace the
2599 # CensoredData instance with a regular array.
2600 data = data._uncensored
2601 censored = False
2603 Narg = len(args)
2604 if Narg > self.numargs:
2605 raise TypeError("Too many input arguments.")
2607 # Check the finiteness of data only if data is not an instance of
2608 # CensoredData. The arrays in a CensoredData instance have already
2609 # been validated.
2610 if not censored:
2611 # Note: `ravel()` is called for backwards compatibility.
2612 data = np.asarray(data).ravel()
2613 if not np.isfinite(data).all():
2614 raise ValueError("The data contains non-finite values.")
2616 start = [None]*2
2617 if (Narg < self.numargs) or not ('loc' in kwds and
2618 'scale' in kwds):
2619 # get distribution specific starting locations
2620 start = self._fitstart(data)
2621 args += start[Narg:-2]
2622 loc = kwds.pop('loc', start[-2])
2623 scale = kwds.pop('scale', start[-1])
2624 args += (loc, scale)
2625 x0, func, restore, args = self._reduce_func(args, kwds, data=data)
2626 optimizer = kwds.pop('optimizer', optimize.fmin)
2627 # convert string to function in scipy.optimize
2628 optimizer = _fit_determine_optimizer(optimizer)
2629 # by now kwds must be empty, since everybody took what they needed
2630 if kwds:
2631 raise TypeError("Unknown arguments: %s." % kwds)
2633 # In some cases, method of moments can be done with fsolve/root
2634 # instead of an optimizer, but sometimes no solution exists,
2635 # especially when the user fixes parameters. Minimizing the sum
2636 # of squares of the error generalizes to these cases.
2637 vals = optimizer(func, x0, args=(data,), disp=0)
2638 obj = func(vals, data)
2640 if restore is not None:
2641 vals = restore(args, vals)
2642 vals = tuple(vals)
2644 loc, scale, shapes = self._unpack_loc_scale(vals)
2645 if not (np.all(self._argcheck(*shapes)) and scale > 0):
2646 raise FitError("Optimization converged to parameters that are "
2647 "outside the range allowed by the distribution.")
2649 if method == 'mm':
2650 if not np.isfinite(obj):
2651 raise FitError("Optimization failed: either a data moment "
2652 "or fitted distribution moment is "
2653 "non-finite.")
2655 return vals
2657 def _fit_loc_scale_support(self, data, *args):
2658 """Estimate loc and scale parameters from data accounting for support.
2660 Parameters
2661 ----------
2662 data : array_like
2663 Data to fit.
2664 arg1, arg2, arg3,... : array_like
2665 The shape parameter(s) for the distribution (see docstring of the
2666 instance object for more information).
2668 Returns
2669 -------
2670 Lhat : float
2671 Estimated location parameter for the data.
2672 Shat : float
2673 Estimated scale parameter for the data.
2675 """
2676 if isinstance(data, CensoredData):
2677 # For this estimate, "uncensor" the data by taking the
2678 # given endpoints as the data for the left- or right-censored
2679 # data, and the mean for the interval-censored data.
2680 data = data._uncensor()
2681 else:
2682 data = np.asarray(data)
2684 # Estimate location and scale according to the method of moments.
2685 loc_hat, scale_hat = self.fit_loc_scale(data, *args)
2687 # Compute the support according to the shape parameters.
2688 self._argcheck(*args)
2689 _a, _b = self._get_support(*args)
2690 a, b = _a, _b
2691 support_width = b - a
2693 # If the support is empty then return the moment-based estimates.
2694 if support_width <= 0:
2695 return loc_hat, scale_hat
2697 # Compute the proposed support according to the loc and scale
2698 # estimates.
2699 a_hat = loc_hat + a * scale_hat
2700 b_hat = loc_hat + b * scale_hat
2702 # Use the moment-based estimates if they are compatible with the data.
2703 data_a = np.min(data)
2704 data_b = np.max(data)
2705 if a_hat < data_a and data_b < b_hat:
2706 return loc_hat, scale_hat
2708 # Otherwise find other estimates that are compatible with the data.
2709 data_width = data_b - data_a
2710 rel_margin = 0.1
2711 margin = data_width * rel_margin
2713 # For a finite interval, both the location and scale
2714 # should have interesting values.
2715 if support_width < np.inf:
2716 loc_hat = (data_a - a) - margin
2717 scale_hat = (data_width + 2 * margin) / support_width
2718 return loc_hat, scale_hat
2720 # For a one-sided interval, use only an interesting location parameter.
2721 if a > -np.inf:
2722 return (data_a - a) - margin, 1
2723 elif b < np.inf:
2724 return (data_b - b) + margin, 1
2725 else:
2726 raise RuntimeError
2728 def fit_loc_scale(self, data, *args):
2729 """
2730 Estimate loc and scale parameters from data using 1st and 2nd moments.
2732 Parameters
2733 ----------
2734 data : array_like
2735 Data to fit.
2736 arg1, arg2, arg3,... : array_like
2737 The shape parameter(s) for the distribution (see docstring of the
2738 instance object for more information).
2740 Returns
2741 -------
2742 Lhat : float
2743 Estimated location parameter for the data.
2744 Shat : float
2745 Estimated scale parameter for the data.
2747 """
2748 mu, mu2 = self.stats(*args, **{'moments': 'mv'})
2749 tmp = asarray(data)
2750 muhat = tmp.mean()
2751 mu2hat = tmp.var()
2752 Shat = sqrt(mu2hat / mu2)
2753 with np.errstate(invalid='ignore'):
2754 Lhat = muhat - Shat*mu
2755 if not np.isfinite(Lhat):
2756 Lhat = 0
2757 if not (np.isfinite(Shat) and (0 < Shat)):
2758 Shat = 1
2759 return Lhat, Shat
2761 def _entropy(self, *args):
2762 def integ(x):
2763 val = self._pdf(x, *args)
2764 return entr(val)
2766 # upper limit is often inf, so suppress warnings when integrating
2767 _a, _b = self._get_support(*args)
2768 with np.errstate(over='ignore'):
2769 h = integrate.quad(integ, _a, _b)[0]
2771 if not np.isnan(h):
2772 return h
2773 else:
2774 # try with different limits if integration problems
2775 low, upp = self.ppf([1e-10, 1. - 1e-10], *args)
2776 if np.isinf(_b):
2777 upper = upp
2778 else:
2779 upper = _b
2780 if np.isinf(_a):
2781 lower = low
2782 else:
2783 lower = _a
2784 return integrate.quad(integ, lower, upper)[0]
2786 def expect(self, func=None, args=(), loc=0, scale=1, lb=None, ub=None,
2787 conditional=False, **kwds):
2788 """Calculate expected value of a function with respect to the
2789 distribution by numerical integration.
2791 The expected value of a function ``f(x)`` with respect to a
2792 distribution ``dist`` is defined as::
2794 ub
2795 E[f(x)] = Integral(f(x) * dist.pdf(x)),
2796 lb
2798 where ``ub`` and ``lb`` are arguments and ``x`` has the ``dist.pdf(x)``
2799 distribution. If the bounds ``lb`` and ``ub`` correspond to the
2800 support of the distribution, e.g. ``[-inf, inf]`` in the default
2801 case, then the integral is the unrestricted expectation of ``f(x)``.
2802 Also, the function ``f(x)`` may be defined such that ``f(x)`` is ``0``
2803 outside a finite interval in which case the expectation is
2804 calculated within the finite range ``[lb, ub]``.
2806 Parameters
2807 ----------
2808 func : callable, optional
2809 Function for which integral is calculated. Takes only one argument.
2810 The default is the identity mapping f(x) = x.
2811 args : tuple, optional
2812 Shape parameters of the distribution.
2813 loc : float, optional
2814 Location parameter (default=0).
2815 scale : float, optional
2816 Scale parameter (default=1).
2817 lb, ub : scalar, optional
2818 Lower and upper bound for integration. Default is set to the
2819 support of the distribution.
2820 conditional : bool, optional
2821 If True, the integral is corrected by the conditional probability
2822 of the integration interval. The return value is the expectation
2823 of the function, conditional on being in the given interval.
2824 Default is False.
2826 Additional keyword arguments are passed to the integration routine.
2828 Returns
2829 -------
2830 expect : float
2831 The calculated expected value.
2833 Notes
2834 -----
2835 The integration behavior of this function is inherited from
2836 `scipy.integrate.quad`. Neither this function nor
2837 `scipy.integrate.quad` can verify whether the integral exists or is
2838 finite. For example ``cauchy(0).mean()`` returns ``np.nan`` and
2839 ``cauchy(0).expect()`` returns ``0.0``.
2841 Likewise, the accuracy of results is not verified by the function.
2842 `scipy.integrate.quad` is typically reliable for integrals that are
2843 numerically favorable, but it is not guaranteed to converge
2844 to a correct value for all possible intervals and integrands. This
2845 function is provided for convenience; for critical applications,
2846 check results against other integration methods.
2848 The function is not vectorized.
2850 Examples
2851 --------
2853 To understand the effect of the bounds of integration consider
2855 >>> from scipy.stats import expon
2856 >>> expon(1).expect(lambda x: 1, lb=0.0, ub=2.0)
2857 0.6321205588285578
2859 This is close to
2861 >>> expon(1).cdf(2.0) - expon(1).cdf(0.0)
2862 0.6321205588285577
2864 If ``conditional=True``
2866 >>> expon(1).expect(lambda x: 1, lb=0.0, ub=2.0, conditional=True)
2867 1.0000000000000002
2869 The slight deviation from 1 is due to numerical integration.
2871 The integrand can be treated as a complex-valued function
2872 by passing ``complex_func=True`` to `scipy.integrate.quad` .
2874 >>> import numpy as np
2875 >>> from scipy.stats import vonmises
2876 >>> res = vonmises(loc=2, kappa=1).expect(lambda x: np.exp(1j*x),
2877 ... complex_func=True)
2878 >>> res
2879 (-0.18576377217422957+0.40590124735052263j)
2881 >>> np.angle(res) # location of the (circular) distribution
2882 2.0
2884 """
2885 lockwds = {'loc': loc,
2886 'scale': scale}
2887 self._argcheck(*args)
2888 _a, _b = self._get_support(*args)
2889 if func is None:
2890 def fun(x, *args):
2891 return x * self.pdf(x, *args, **lockwds)
2892 else:
2893 def fun(x, *args):
2894 return func(x) * self.pdf(x, *args, **lockwds)
2895 if lb is None:
2896 lb = loc + _a * scale
2897 if ub is None:
2898 ub = loc + _b * scale
2900 cdf_bounds = self.cdf([lb, ub], *args, **lockwds)
2901 invfac = cdf_bounds[1] - cdf_bounds[0]
2903 kwds['args'] = args
2905 # split interval to help integrator w/ infinite support; see gh-8928
2906 alpha = 0.05 # split body from tails at probability mass `alpha`
2907 inner_bounds = np.array([alpha, 1-alpha])
2908 cdf_inner_bounds = cdf_bounds[0] + invfac * inner_bounds
2909 c, d = loc + self._ppf(cdf_inner_bounds, *args) * scale
2911 # Do not silence warnings from integration.
2912 lbc = integrate.quad(fun, lb, c, **kwds)[0]
2913 cd = integrate.quad(fun, c, d, **kwds)[0]
2914 dub = integrate.quad(fun, d, ub, **kwds)[0]
2915 vals = (lbc + cd + dub)
2917 if conditional:
2918 vals /= invfac
2919 return np.array(vals)[()] # make it a numpy scalar like other methods
2921 def _param_info(self):
2922 shape_info = self._shape_info()
2923 loc_info = _ShapeInfo("loc", False, (-np.inf, np.inf), (False, False))
2924 scale_info = _ShapeInfo("scale", False, (0, np.inf), (False, False))
2925 param_info = shape_info + [loc_info, scale_info]
2926 return param_info
2928 # For now, _delta_cdf is a private method.
2929 def _delta_cdf(self, x1, x2, *args, loc=0, scale=1):
2930 """
2931 Compute CDF(x2) - CDF(x1).
2933 Where x1 is greater than the median, compute SF(x1) - SF(x2),
2934 otherwise compute CDF(x2) - CDF(x1).
2936 This function is only useful if `dist.sf(x, ...)` has an implementation
2937 that is numerically more accurate than `1 - dist.cdf(x, ...)`.
2938 """
2939 cdf1 = self.cdf(x1, *args, loc=loc, scale=scale)
2940 # Possible optimizations (needs investigation-these might not be
2941 # better):
2942 # * Use _lazywhere instead of np.where
2943 # * Instead of cdf1 > 0.5, compare x1 to the median.
2944 result = np.where(cdf1 > 0.5,
2945 (self.sf(x1, *args, loc=loc, scale=scale)
2946 - self.sf(x2, *args, loc=loc, scale=scale)),
2947 self.cdf(x2, *args, loc=loc, scale=scale) - cdf1)
2948 if result.ndim == 0:
2949 result = result[()]
2950 return result
2953# Helpers for the discrete distributions
2954def _drv2_moment(self, n, *args):
2955 """Non-central moment of discrete distribution."""
2956 def fun(x):
2957 return np.power(x, n) * self._pmf(x, *args)
2959 _a, _b = self._get_support(*args)
2960 return _expect(fun, _a, _b, self.ppf(0.5, *args), self.inc)
2963def _drv2_ppfsingle(self, q, *args): # Use basic bisection algorithm
2964 _a, _b = self._get_support(*args)
2965 b = _b
2966 a = _a
2967 if isinf(b): # Be sure ending point is > q
2968 b = int(max(100*q, 10))
2969 while 1:
2970 if b >= _b:
2971 qb = 1.0
2972 break
2973 qb = self._cdf(b, *args)
2974 if (qb < q):
2975 b += 10
2976 else:
2977 break
2978 else:
2979 qb = 1.0
2980 if isinf(a): # be sure starting point < q
2981 a = int(min(-100*q, -10))
2982 while 1:
2983 if a <= _a:
2984 qb = 0.0
2985 break
2986 qa = self._cdf(a, *args)
2987 if (qa > q):
2988 a -= 10
2989 else:
2990 break
2991 else:
2992 qa = self._cdf(a, *args)
2994 while 1:
2995 if (qa == q):
2996 return a
2997 if (qb == q):
2998 return b
2999 if b <= a+1:
3000 if qa > q:
3001 return a
3002 else:
3003 return b
3004 c = int((a+b)/2.0)
3005 qc = self._cdf(c, *args)
3006 if (qc < q):
3007 if a != c:
3008 a = c
3009 else:
3010 raise RuntimeError('updating stopped, endless loop')
3011 qa = qc
3012 elif (qc > q):
3013 if b != c:
3014 b = c
3015 else:
3016 raise RuntimeError('updating stopped, endless loop')
3017 qb = qc
3018 else:
3019 return c
3022# Must over-ride one of _pmf or _cdf or pass in
3023# x_k, p(x_k) lists in initialization
3026class rv_discrete(rv_generic):
3027 """A generic discrete random variable class meant for subclassing.
3029 `rv_discrete` is a base class to construct specific distribution classes
3030 and instances for discrete random variables. It can also be used
3031 to construct an arbitrary distribution defined by a list of support
3032 points and corresponding probabilities.
3034 Parameters
3035 ----------
3036 a : float, optional
3037 Lower bound of the support of the distribution, default: 0
3038 b : float, optional
3039 Upper bound of the support of the distribution, default: plus infinity
3040 moment_tol : float, optional
3041 The tolerance for the generic calculation of moments.
3042 values : tuple of two array_like, optional
3043 ``(xk, pk)`` where ``xk`` are integers and ``pk`` are the non-zero
3044 probabilities between 0 and 1 with ``sum(pk) = 1``. ``xk``
3045 and ``pk`` must have the same shape.
3046 inc : integer, optional
3047 Increment for the support of the distribution.
3048 Default is 1. (other values have not been tested)
3049 badvalue : float, optional
3050 The value in a result arrays that indicates a value that for which
3051 some argument restriction is violated, default is np.nan.
3052 name : str, optional
3053 The name of the instance. This string is used to construct the default
3054 example for distributions.
3055 longname : str, optional
3056 This string is used as part of the first line of the docstring returned
3057 when a subclass has no docstring of its own. Note: `longname` exists
3058 for backwards compatibility, do not use for new subclasses.
3059 shapes : str, optional
3060 The shape of the distribution. For example "m, n" for a distribution
3061 that takes two integers as the two shape arguments for all its methods
3062 If not provided, shape parameters will be inferred from
3063 the signatures of the private methods, ``_pmf`` and ``_cdf`` of
3064 the instance.
3065 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
3066 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
3067 singleton is used.
3068 If `seed` is an int, a new ``RandomState`` instance is used,
3069 seeded with `seed`.
3070 If `seed` is already a ``Generator`` or ``RandomState`` instance then
3071 that instance is used.
3073 Methods
3074 -------
3075 rvs
3076 pmf
3077 logpmf
3078 cdf
3079 logcdf
3080 sf
3081 logsf
3082 ppf
3083 isf
3084 moment
3085 stats
3086 entropy
3087 expect
3088 median
3089 mean
3090 std
3091 var
3092 interval
3093 __call__
3094 support
3096 Notes
3097 -----
3098 This class is similar to `rv_continuous`. Whether a shape parameter is
3099 valid is decided by an ``_argcheck`` method (which defaults to checking
3100 that its arguments are strictly positive.)
3101 The main differences are:
3103 - the support of the distribution is a set of integers
3104 - instead of the probability density function, ``pdf`` (and the
3105 corresponding private ``_pdf``), this class defines the
3106 *probability mass function*, `pmf` (and the corresponding
3107 private ``_pmf``.)
3108 - scale parameter is not defined.
3110 To create a new discrete distribution, we would do the following:
3112 >>> from scipy.stats import rv_discrete
3113 >>> class poisson_gen(rv_discrete):
3114 ... "Poisson distribution"
3115 ... def _pmf(self, k, mu):
3116 ... return exp(-mu) * mu**k / factorial(k)
3118 and create an instance::
3120 >>> poisson = poisson_gen(name="poisson")
3122 Note that above we defined the Poisson distribution in the standard form.
3123 Shifting the distribution can be done by providing the ``loc`` parameter
3124 to the methods of the instance. For example, ``poisson.pmf(x, mu, loc)``
3125 delegates the work to ``poisson._pmf(x-loc, mu)``.
3127 **Discrete distributions from a list of probabilities**
3129 Alternatively, you can construct an arbitrary discrete rv defined
3130 on a finite set of values ``xk`` with ``Prob{X=xk} = pk`` by using the
3131 ``values`` keyword argument to the `rv_discrete` constructor.
3133 **Deepcopying / Pickling**
3135 If a distribution or frozen distribution is deepcopied (pickled/unpickled,
3136 etc.), any underlying random number generator is deepcopied with it. An
3137 implication is that if a distribution relies on the singleton RandomState
3138 before copying, it will rely on a copy of that random state after copying,
3139 and ``np.random.seed`` will no longer control the state.
3141 Examples
3142 --------
3143 Custom made discrete distribution:
3145 >>> import numpy as np
3146 >>> from scipy import stats
3147 >>> xk = np.arange(7)
3148 >>> pk = (0.1, 0.2, 0.3, 0.1, 0.1, 0.0, 0.2)
3149 >>> custm = stats.rv_discrete(name='custm', values=(xk, pk))
3150 >>>
3151 >>> import matplotlib.pyplot as plt
3152 >>> fig, ax = plt.subplots(1, 1)
3153 >>> ax.plot(xk, custm.pmf(xk), 'ro', ms=12, mec='r')
3154 >>> ax.vlines(xk, 0, custm.pmf(xk), colors='r', lw=4)
3155 >>> plt.show()
3157 Random number generation:
3159 >>> R = custm.rvs(size=100)
3161 """
3162 def __new__(cls, a=0, b=inf, name=None, badvalue=None,
3163 moment_tol=1e-8, values=None, inc=1, longname=None,
3164 shapes=None, seed=None):
3166 if values is not None:
3167 # dispatch to a subclass
3168 return super().__new__(rv_sample)
3169 else:
3170 # business as usual
3171 return super().__new__(cls)
3173 def __init__(self, a=0, b=inf, name=None, badvalue=None,
3174 moment_tol=1e-8, values=None, inc=1, longname=None,
3175 shapes=None, seed=None):
3177 super().__init__(seed)
3179 # cf generic freeze
3180 self._ctor_param = dict(
3181 a=a, b=b, name=name, badvalue=badvalue,
3182 moment_tol=moment_tol, values=values, inc=inc,
3183 longname=longname, shapes=shapes, seed=seed)
3185 if badvalue is None:
3186 badvalue = nan
3187 self.badvalue = badvalue
3188 self.a = a
3189 self.b = b
3190 self.moment_tol = moment_tol
3191 self.inc = inc
3192 self.shapes = shapes
3194 if values is not None:
3195 raise ValueError("rv_discrete.__init__(..., values != None, ...)")
3197 self._construct_argparser(meths_to_inspect=[self._pmf, self._cdf],
3198 locscale_in='loc=0',
3199 # scale=1 for discrete RVs
3200 locscale_out='loc, 1')
3201 self._attach_methods()
3202 self._construct_docstrings(name, longname)
3204 def __getstate__(self):
3205 dct = self.__dict__.copy()
3206 # these methods will be remade in __setstate__
3207 attrs = ["_parse_args", "_parse_args_stats", "_parse_args_rvs",
3208 "_cdfvec", "_ppfvec", "generic_moment"]
3209 [dct.pop(attr, None) for attr in attrs]
3210 return dct
3212 def _attach_methods(self):
3213 """Attaches dynamically created methods to the rv_discrete instance."""
3214 self._cdfvec = vectorize(self._cdf_single, otypes='d')
3215 self.vecentropy = vectorize(self._entropy)
3217 # _attach_methods is responsible for calling _attach_argparser_methods
3218 self._attach_argparser_methods()
3220 # nin correction needs to be after we know numargs
3221 # correct nin for generic moment vectorization
3222 _vec_generic_moment = vectorize(_drv2_moment, otypes='d')
3223 _vec_generic_moment.nin = self.numargs + 2
3224 self.generic_moment = types.MethodType(_vec_generic_moment, self)
3226 # correct nin for ppf vectorization
3227 _vppf = vectorize(_drv2_ppfsingle, otypes='d')
3228 _vppf.nin = self.numargs + 2
3229 self._ppfvec = types.MethodType(_vppf, self)
3231 # now that self.numargs is defined, we can adjust nin
3232 self._cdfvec.nin = self.numargs + 1
3234 def _construct_docstrings(self, name, longname):
3235 if name is None:
3236 name = 'Distribution'
3237 self.name = name
3239 # generate docstring for subclass instances
3240 if longname is None:
3241 if name[0] in ['aeiouAEIOU']:
3242 hstr = "An "
3243 else:
3244 hstr = "A "
3245 longname = hstr + name
3247 if sys.flags.optimize < 2:
3248 # Skip adding docstrings if interpreter is run with -OO
3249 if self.__doc__ is None:
3250 self._construct_default_doc(longname=longname,
3251 docdict=docdict_discrete,
3252 discrete='discrete')
3253 else:
3254 dct = dict(distdiscrete)
3255 self._construct_doc(docdict_discrete, dct.get(self.name))
3257 # discrete RV do not have the scale parameter, remove it
3258 self.__doc__ = self.__doc__.replace(
3259 '\n scale : array_like, '
3260 'optional\n scale parameter (default=1)', '')
3262 def _updated_ctor_param(self):
3263 """Return the current version of _ctor_param, possibly updated by user.
3265 Used by freezing.
3266 Keep this in sync with the signature of __init__.
3267 """
3268 dct = self._ctor_param.copy()
3269 dct['a'] = self.a
3270 dct['b'] = self.b
3271 dct['badvalue'] = self.badvalue
3272 dct['moment_tol'] = self.moment_tol
3273 dct['inc'] = self.inc
3274 dct['name'] = self.name
3275 dct['shapes'] = self.shapes
3276 return dct
3278 def _nonzero(self, k, *args):
3279 return floor(k) == k
3281 def _pmf(self, k, *args):
3282 return self._cdf(k, *args) - self._cdf(k-1, *args)
3284 def _logpmf(self, k, *args):
3285 return log(self._pmf(k, *args))
3287 def _logpxf(self, k, *args):
3288 # continuous distributions have PDF, discrete have PMF, but sometimes
3289 # the distinction doesn't matter. This lets us use `_logpxf` for both
3290 # discrete and continuous distributions.
3291 return self._logpmf(k, *args)
3293 def _unpack_loc_scale(self, theta):
3294 try:
3295 loc = theta[-1]
3296 scale = 1
3297 args = tuple(theta[:-1])
3298 except IndexError as e:
3299 raise ValueError("Not enough input arguments.") from e
3300 return loc, scale, args
3302 def _cdf_single(self, k, *args):
3303 _a, _b = self._get_support(*args)
3304 m = arange(int(_a), k+1)
3305 return np.sum(self._pmf(m, *args), axis=0)
3307 def _cdf(self, x, *args):
3308 k = floor(x)
3309 return self._cdfvec(k, *args)
3311 # generic _logcdf, _sf, _logsf, _ppf, _isf, _rvs defined in rv_generic
3313 def rvs(self, *args, **kwargs):
3314 """Random variates of given type.
3316 Parameters
3317 ----------
3318 arg1, arg2, arg3,... : array_like
3319 The shape parameter(s) for the distribution (see docstring of the
3320 instance object for more information).
3321 loc : array_like, optional
3322 Location parameter (default=0).
3323 size : int or tuple of ints, optional
3324 Defining number of random variates (Default is 1). Note that `size`
3325 has to be given as keyword, not as positional argument.
3326 random_state : {None, int, `numpy.random.Generator`,
3327 `numpy.random.RandomState`}, optional
3329 If `random_state` is None (or `np.random`), the
3330 `numpy.random.RandomState` singleton is used.
3331 If `random_state` is an int, a new ``RandomState`` instance is
3332 used, seeded with `random_state`.
3333 If `random_state` is already a ``Generator`` or ``RandomState``
3334 instance, that instance is used.
3336 Returns
3337 -------
3338 rvs : ndarray or scalar
3339 Random variates of given `size`.
3341 """
3342 kwargs['discrete'] = True
3343 return super().rvs(*args, **kwargs)
3345 def pmf(self, k, *args, **kwds):
3346 """Probability mass function at k of the given RV.
3348 Parameters
3349 ----------
3350 k : array_like
3351 Quantiles.
3352 arg1, arg2, arg3,... : array_like
3353 The shape parameter(s) for the distribution (see docstring of the
3354 instance object for more information)
3355 loc : array_like, optional
3356 Location parameter (default=0).
3358 Returns
3359 -------
3360 pmf : array_like
3361 Probability mass function evaluated at k
3363 """
3364 args, loc, _ = self._parse_args(*args, **kwds)
3365 k, loc = map(asarray, (k, loc))
3366 args = tuple(map(asarray, args))
3367 _a, _b = self._get_support(*args)
3368 k = asarray(k-loc)
3369 cond0 = self._argcheck(*args)
3370 cond1 = (k >= _a) & (k <= _b)
3371 if not isinstance(self, rv_sample):
3372 cond1 = cond1 & self._nonzero(k, *args)
3373 cond = cond0 & cond1
3374 output = zeros(shape(cond), 'd')
3375 place(output, (1-cond0) + np.isnan(k), self.badvalue)
3376 if np.any(cond):
3377 goodargs = argsreduce(cond, *((k,)+args))
3378 place(output, cond, np.clip(self._pmf(*goodargs), 0, 1))
3379 if output.ndim == 0:
3380 return output[()]
3381 return output
3383 def logpmf(self, k, *args, **kwds):
3384 """Log of the probability mass function at k of the given RV.
3386 Parameters
3387 ----------
3388 k : array_like
3389 Quantiles.
3390 arg1, arg2, arg3,... : array_like
3391 The shape parameter(s) for the distribution (see docstring of the
3392 instance object for more information).
3393 loc : array_like, optional
3394 Location parameter. Default is 0.
3396 Returns
3397 -------
3398 logpmf : array_like
3399 Log of the probability mass function evaluated at k.
3401 """
3402 args, loc, _ = self._parse_args(*args, **kwds)
3403 k, loc = map(asarray, (k, loc))
3404 args = tuple(map(asarray, args))
3405 _a, _b = self._get_support(*args)
3406 k = asarray(k-loc)
3407 cond0 = self._argcheck(*args)
3408 cond1 = (k >= _a) & (k <= _b)
3409 if not isinstance(self, rv_sample):
3410 cond1 = cond1 & self._nonzero(k, *args)
3411 cond = cond0 & cond1
3412 output = empty(shape(cond), 'd')
3413 output.fill(NINF)
3414 place(output, (1-cond0) + np.isnan(k), self.badvalue)
3415 if np.any(cond):
3416 goodargs = argsreduce(cond, *((k,)+args))
3417 place(output, cond, self._logpmf(*goodargs))
3418 if output.ndim == 0:
3419 return output[()]
3420 return output
3422 def cdf(self, k, *args, **kwds):
3423 """Cumulative distribution function of the given RV.
3425 Parameters
3426 ----------
3427 k : array_like, int
3428 Quantiles.
3429 arg1, arg2, arg3,... : array_like
3430 The shape parameter(s) for the distribution (see docstring of the
3431 instance object for more information).
3432 loc : array_like, optional
3433 Location parameter (default=0).
3435 Returns
3436 -------
3437 cdf : ndarray
3438 Cumulative distribution function evaluated at `k`.
3440 """
3441 args, loc, _ = self._parse_args(*args, **kwds)
3442 k, loc = map(asarray, (k, loc))
3443 args = tuple(map(asarray, args))
3444 _a, _b = self._get_support(*args)
3445 k = asarray(k-loc)
3446 cond0 = self._argcheck(*args)
3447 cond1 = (k >= _a) & (k < _b)
3448 cond2 = (k >= _b)
3449 cond3 = np.isneginf(k)
3450 cond = cond0 & cond1 & np.isfinite(k)
3452 output = zeros(shape(cond), 'd')
3453 place(output, cond2*(cond0 == cond0), 1.0)
3454 place(output, cond3*(cond0 == cond0), 0.0)
3455 place(output, (1-cond0) + np.isnan(k), self.badvalue)
3457 if np.any(cond):
3458 goodargs = argsreduce(cond, *((k,)+args))
3459 place(output, cond, np.clip(self._cdf(*goodargs), 0, 1))
3460 if output.ndim == 0:
3461 return output[()]
3462 return output
3464 def logcdf(self, k, *args, **kwds):
3465 """Log of the cumulative distribution function at k of the given RV.
3467 Parameters
3468 ----------
3469 k : array_like, int
3470 Quantiles.
3471 arg1, arg2, arg3,... : array_like
3472 The shape parameter(s) for the distribution (see docstring of the
3473 instance object for more information).
3474 loc : array_like, optional
3475 Location parameter (default=0).
3477 Returns
3478 -------
3479 logcdf : array_like
3480 Log of the cumulative distribution function evaluated at k.
3482 """
3483 args, loc, _ = self._parse_args(*args, **kwds)
3484 k, loc = map(asarray, (k, loc))
3485 args = tuple(map(asarray, args))
3486 _a, _b = self._get_support(*args)
3487 k = asarray(k-loc)
3488 cond0 = self._argcheck(*args)
3489 cond1 = (k >= _a) & (k < _b)
3490 cond2 = (k >= _b)
3491 cond = cond0 & cond1
3492 output = empty(shape(cond), 'd')
3493 output.fill(NINF)
3494 place(output, (1-cond0) + np.isnan(k), self.badvalue)
3495 place(output, cond2*(cond0 == cond0), 0.0)
3497 if np.any(cond):
3498 goodargs = argsreduce(cond, *((k,)+args))
3499 place(output, cond, self._logcdf(*goodargs))
3500 if output.ndim == 0:
3501 return output[()]
3502 return output
3504 def sf(self, k, *args, **kwds):
3505 """Survival function (1 - `cdf`) at k of the given RV.
3507 Parameters
3508 ----------
3509 k : array_like
3510 Quantiles.
3511 arg1, arg2, arg3,... : array_like
3512 The shape parameter(s) for the distribution (see docstring of the
3513 instance object for more information).
3514 loc : array_like, optional
3515 Location parameter (default=0).
3517 Returns
3518 -------
3519 sf : array_like
3520 Survival function evaluated at k.
3522 """
3523 args, loc, _ = self._parse_args(*args, **kwds)
3524 k, loc = map(asarray, (k, loc))
3525 args = tuple(map(asarray, args))
3526 _a, _b = self._get_support(*args)
3527 k = asarray(k-loc)
3528 cond0 = self._argcheck(*args)
3529 cond1 = (k >= _a) & (k < _b)
3530 cond2 = ((k < _a) | np.isneginf(k)) & cond0
3531 cond = cond0 & cond1 & np.isfinite(k)
3532 output = zeros(shape(cond), 'd')
3533 place(output, (1-cond0) + np.isnan(k), self.badvalue)
3534 place(output, cond2, 1.0)
3535 if np.any(cond):
3536 goodargs = argsreduce(cond, *((k,)+args))
3537 place(output, cond, np.clip(self._sf(*goodargs), 0, 1))
3538 if output.ndim == 0:
3539 return output[()]
3540 return output
3542 def logsf(self, k, *args, **kwds):
3543 """Log of the survival function of the given RV.
3545 Returns the log of the "survival function," defined as 1 - `cdf`,
3546 evaluated at `k`.
3548 Parameters
3549 ----------
3550 k : array_like
3551 Quantiles.
3552 arg1, arg2, arg3,... : array_like
3553 The shape parameter(s) for the distribution (see docstring of the
3554 instance object for more information).
3555 loc : array_like, optional
3556 Location parameter (default=0).
3558 Returns
3559 -------
3560 logsf : ndarray
3561 Log of the survival function evaluated at `k`.
3563 """
3564 args, loc, _ = self._parse_args(*args, **kwds)
3565 k, loc = map(asarray, (k, loc))
3566 args = tuple(map(asarray, args))
3567 _a, _b = self._get_support(*args)
3568 k = asarray(k-loc)
3569 cond0 = self._argcheck(*args)
3570 cond1 = (k >= _a) & (k < _b)
3571 cond2 = (k < _a) & cond0
3572 cond = cond0 & cond1
3573 output = empty(shape(cond), 'd')
3574 output.fill(NINF)
3575 place(output, (1-cond0) + np.isnan(k), self.badvalue)
3576 place(output, cond2, 0.0)
3577 if np.any(cond):
3578 goodargs = argsreduce(cond, *((k,)+args))
3579 place(output, cond, self._logsf(*goodargs))
3580 if output.ndim == 0:
3581 return output[()]
3582 return output
3584 def ppf(self, q, *args, **kwds):
3585 """Percent point function (inverse of `cdf`) at q of the given RV.
3587 Parameters
3588 ----------
3589 q : array_like
3590 Lower tail probability.
3591 arg1, arg2, arg3,... : array_like
3592 The shape parameter(s) for the distribution (see docstring of the
3593 instance object for more information).
3594 loc : array_like, optional
3595 Location parameter (default=0).
3597 Returns
3598 -------
3599 k : array_like
3600 Quantile corresponding to the lower tail probability, q.
3602 """
3603 args, loc, _ = self._parse_args(*args, **kwds)
3604 q, loc = map(asarray, (q, loc))
3605 args = tuple(map(asarray, args))
3606 _a, _b = self._get_support(*args)
3607 cond0 = self._argcheck(*args) & (loc == loc)
3608 cond1 = (q > 0) & (q < 1)
3609 cond2 = (q == 1) & cond0
3610 cond = cond0 & cond1
3611 output = np.full(shape(cond), fill_value=self.badvalue, dtype='d')
3612 # output type 'd' to handle nin and inf
3613 place(output, (q == 0)*(cond == cond), _a-1 + loc)
3614 place(output, cond2, _b + loc)
3615 if np.any(cond):
3616 goodargs = argsreduce(cond, *((q,)+args+(loc,)))
3617 loc, goodargs = goodargs[-1], goodargs[:-1]
3618 place(output, cond, self._ppf(*goodargs) + loc)
3620 if output.ndim == 0:
3621 return output[()]
3622 return output
3624 def isf(self, q, *args, **kwds):
3625 """Inverse survival function (inverse of `sf`) at q of the given RV.
3627 Parameters
3628 ----------
3629 q : array_like
3630 Upper tail probability.
3631 arg1, arg2, arg3,... : array_like
3632 The shape parameter(s) for the distribution (see docstring of the
3633 instance object for more information).
3634 loc : array_like, optional
3635 Location parameter (default=0).
3637 Returns
3638 -------
3639 k : ndarray or scalar
3640 Quantile corresponding to the upper tail probability, q.
3642 """
3643 args, loc, _ = self._parse_args(*args, **kwds)
3644 q, loc = map(asarray, (q, loc))
3645 args = tuple(map(asarray, args))
3646 _a, _b = self._get_support(*args)
3647 cond0 = self._argcheck(*args) & (loc == loc)
3648 cond1 = (q > 0) & (q < 1)
3649 cond2 = (q == 1) & cond0
3650 cond3 = (q == 0) & cond0
3651 cond = cond0 & cond1
3653 # same problem as with ppf; copied from ppf and changed
3654 output = np.full(shape(cond), fill_value=self.badvalue, dtype='d')
3655 # output type 'd' to handle nin and inf
3656 lower_bound = _a - 1 + loc
3657 upper_bound = _b + loc
3658 place(output, cond2*(cond == cond), lower_bound)
3659 place(output, cond3*(cond == cond), upper_bound)
3661 # call place only if at least 1 valid argument
3662 if np.any(cond):
3663 goodargs = argsreduce(cond, *((q,)+args+(loc,)))
3664 loc, goodargs = goodargs[-1], goodargs[:-1]
3665 # PB same as ticket 766
3666 place(output, cond, self._isf(*goodargs) + loc)
3668 if output.ndim == 0:
3669 return output[()]
3670 return output
3672 def _entropy(self, *args):
3673 if hasattr(self, 'pk'):
3674 return stats.entropy(self.pk)
3675 else:
3676 _a, _b = self._get_support(*args)
3677 return _expect(lambda x: entr(self.pmf(x, *args)),
3678 _a, _b, self.ppf(0.5, *args), self.inc)
3680 def expect(self, func=None, args=(), loc=0, lb=None, ub=None,
3681 conditional=False, maxcount=1000, tolerance=1e-10, chunksize=32):
3682 """
3683 Calculate expected value of a function with respect to the distribution
3684 for discrete distribution by numerical summation.
3686 Parameters
3687 ----------
3688 func : callable, optional
3689 Function for which the expectation value is calculated.
3690 Takes only one argument.
3691 The default is the identity mapping f(k) = k.
3692 args : tuple, optional
3693 Shape parameters of the distribution.
3694 loc : float, optional
3695 Location parameter.
3696 Default is 0.
3697 lb, ub : int, optional
3698 Lower and upper bound for the summation, default is set to the
3699 support of the distribution, inclusive (``lb <= k <= ub``).
3700 conditional : bool, optional
3701 If true then the expectation is corrected by the conditional
3702 probability of the summation interval. The return value is the
3703 expectation of the function, `func`, conditional on being in
3704 the given interval (k such that ``lb <= k <= ub``).
3705 Default is False.
3706 maxcount : int, optional
3707 Maximal number of terms to evaluate (to avoid an endless loop for
3708 an infinite sum). Default is 1000.
3709 tolerance : float, optional
3710 Absolute tolerance for the summation. Default is 1e-10.
3711 chunksize : int, optional
3712 Iterate over the support of a distributions in chunks of this size.
3713 Default is 32.
3715 Returns
3716 -------
3717 expect : float
3718 Expected value.
3720 Notes
3721 -----
3722 For heavy-tailed distributions, the expected value may or
3723 may not exist,
3724 depending on the function, `func`. If it does exist, but the
3725 sum converges
3726 slowly, the accuracy of the result may be rather low. For instance, for
3727 ``zipf(4)``, accuracy for mean, variance in example is only 1e-5.
3728 increasing `maxcount` and/or `chunksize` may improve the result,
3729 but may also make zipf very slow.
3731 The function is not vectorized.
3733 """
3734 if func is None:
3735 def fun(x):
3736 # loc and args from outer scope
3737 return (x+loc)*self._pmf(x, *args)
3738 else:
3739 def fun(x):
3740 # loc and args from outer scope
3741 return func(x+loc)*self._pmf(x, *args)
3742 # used pmf because _pmf does not check support in randint and there
3743 # might be problems(?) with correct self.a, self.b at this stage maybe
3744 # not anymore, seems to work now with _pmf
3746 _a, _b = self._get_support(*args)
3747 if lb is None:
3748 lb = _a
3749 else:
3750 lb = lb - loc # convert bound for standardized distribution
3751 if ub is None:
3752 ub = _b
3753 else:
3754 ub = ub - loc # convert bound for standardized distribution
3755 if conditional:
3756 invfac = self.sf(lb-1, *args) - self.sf(ub, *args)
3757 else:
3758 invfac = 1.0
3760 if isinstance(self, rv_sample):
3761 res = self._expect(fun, lb, ub)
3762 return res / invfac
3764 # iterate over the support, starting from the median
3765 x0 = self.ppf(0.5, *args)
3766 res = _expect(fun, lb, ub, x0, self.inc, maxcount, tolerance, chunksize)
3767 return res / invfac
3769 def _param_info(self):
3770 shape_info = self._shape_info()
3771 loc_info = _ShapeInfo("loc", True, (-np.inf, np.inf), (False, False))
3772 param_info = shape_info + [loc_info]
3773 return param_info
3776def _expect(fun, lb, ub, x0, inc, maxcount=1000, tolerance=1e-10,
3777 chunksize=32):
3778 """Helper for computing the expectation value of `fun`."""
3779 # short-circuit if the support size is small enough
3780 if (ub - lb) <= chunksize:
3781 supp = np.arange(lb, ub+1, inc)
3782 vals = fun(supp)
3783 return np.sum(vals)
3785 # otherwise, iterate starting from x0
3786 if x0 < lb:
3787 x0 = lb
3788 if x0 > ub:
3789 x0 = ub
3791 count, tot = 0, 0.
3792 # iterate over [x0, ub] inclusive
3793 for x in _iter_chunked(x0, ub+1, chunksize=chunksize, inc=inc):
3794 count += x.size
3795 delta = np.sum(fun(x))
3796 tot += delta
3797 if abs(delta) < tolerance * x.size:
3798 break
3799 if count > maxcount:
3800 warnings.warn('expect(): sum did not converge', RuntimeWarning)
3801 return tot
3803 # iterate over [lb, x0)
3804 for x in _iter_chunked(x0-1, lb-1, chunksize=chunksize, inc=-inc):
3805 count += x.size
3806 delta = np.sum(fun(x))
3807 tot += delta
3808 if abs(delta) < tolerance * x.size:
3809 break
3810 if count > maxcount:
3811 warnings.warn('expect(): sum did not converge', RuntimeWarning)
3812 break
3814 return tot
3817def _iter_chunked(x0, x1, chunksize=4, inc=1):
3818 """Iterate from x0 to x1 in chunks of chunksize and steps inc.
3820 x0 must be finite, x1 need not be. In the latter case, the iterator is
3821 infinite.
3822 Handles both x0 < x1 and x0 > x1. In the latter case, iterates downwards
3823 (make sure to set inc < 0.)
3825 >>> [x for x in _iter_chunked(2, 5, inc=2)]
3826 [array([2, 4])]
3827 >>> [x for x in _iter_chunked(2, 11, inc=2)]
3828 [array([2, 4, 6, 8]), array([10])]
3829 >>> [x for x in _iter_chunked(2, -5, inc=-2)]
3830 [array([ 2, 0, -2, -4])]
3831 >>> [x for x in _iter_chunked(2, -9, inc=-2)]
3832 [array([ 2, 0, -2, -4]), array([-6, -8])]
3834 """
3835 if inc == 0:
3836 raise ValueError('Cannot increment by zero.')
3837 if chunksize <= 0:
3838 raise ValueError('Chunk size must be positive; got %s.' % chunksize)
3840 s = 1 if inc > 0 else -1
3841 stepsize = abs(chunksize * inc)
3843 x = x0
3844 while (x - x1) * inc < 0:
3845 delta = min(stepsize, abs(x - x1))
3846 step = delta * s
3847 supp = np.arange(x, x + step, inc)
3848 x += step
3849 yield supp
3852class rv_sample(rv_discrete):
3853 """A 'sample' discrete distribution defined by the support and values.
3855 The ctor ignores most of the arguments, only needs the `values` argument.
3856 """
3858 def __init__(self, a=0, b=inf, name=None, badvalue=None,
3859 moment_tol=1e-8, values=None, inc=1, longname=None,
3860 shapes=None, seed=None):
3862 super(rv_discrete, self).__init__(seed)
3864 if values is None:
3865 raise ValueError("rv_sample.__init__(..., values=None,...)")
3867 # cf generic freeze
3868 self._ctor_param = dict(
3869 a=a, b=b, name=name, badvalue=badvalue,
3870 moment_tol=moment_tol, values=values, inc=inc,
3871 longname=longname, shapes=shapes, seed=seed)
3873 if badvalue is None:
3874 badvalue = nan
3875 self.badvalue = badvalue
3876 self.moment_tol = moment_tol
3877 self.inc = inc
3878 self.shapes = shapes
3879 self.vecentropy = self._entropy
3881 xk, pk = values
3883 if np.shape(xk) != np.shape(pk):
3884 raise ValueError("xk and pk must have the same shape.")
3885 if np.less(pk, 0.0).any():
3886 raise ValueError("All elements of pk must be non-negative.")
3887 if not np.allclose(np.sum(pk), 1):
3888 raise ValueError("The sum of provided pk is not 1.")
3890 indx = np.argsort(np.ravel(xk))
3891 self.xk = np.take(np.ravel(xk), indx, 0)
3892 self.pk = np.take(np.ravel(pk), indx, 0)
3893 self.a = self.xk[0]
3894 self.b = self.xk[-1]
3896 self.qvals = np.cumsum(self.pk, axis=0)
3898 self.shapes = ' ' # bypass inspection
3900 self._construct_argparser(meths_to_inspect=[self._pmf],
3901 locscale_in='loc=0',
3902 # scale=1 for discrete RVs
3903 locscale_out='loc, 1')
3905 self._attach_methods()
3907 self._construct_docstrings(name, longname)
3909 def __getstate__(self):
3910 dct = self.__dict__.copy()
3912 # these methods will be remade in rv_generic.__setstate__,
3913 # which calls rv_generic._attach_methods
3914 attrs = ["_parse_args", "_parse_args_stats", "_parse_args_rvs"]
3915 [dct.pop(attr, None) for attr in attrs]
3917 return dct
3919 def _attach_methods(self):
3920 """Attaches dynamically created argparser methods."""
3921 self._attach_argparser_methods()
3923 def _get_support(self, *args):
3924 """Return the support of the (unscaled, unshifted) distribution.
3926 Parameters
3927 ----------
3928 arg1, arg2, ... : array_like
3929 The shape parameter(s) for the distribution (see docstring of the
3930 instance object for more information).
3932 Returns
3933 -------
3934 a, b : numeric (float, or int or +/-np.inf)
3935 end-points of the distribution's support.
3936 """
3937 return self.a, self.b
3939 def _pmf(self, x):
3940 return np.select([x == k for k in self.xk],
3941 [np.broadcast_arrays(p, x)[0] for p in self.pk], 0)
3943 def _cdf(self, x):
3944 xx, xxk = np.broadcast_arrays(x[:, None], self.xk)
3945 indx = np.argmax(xxk > xx, axis=-1) - 1
3946 return self.qvals[indx]
3948 def _ppf(self, q):
3949 qq, sqq = np.broadcast_arrays(q[..., None], self.qvals)
3950 indx = argmax(sqq >= qq, axis=-1)
3951 return self.xk[indx]
3953 def _rvs(self, size=None, random_state=None):
3954 # Need to define it explicitly, otherwise .rvs() with size=None
3955 # fails due to explicit broadcasting in _ppf
3956 U = random_state.uniform(size=size)
3957 if size is None:
3958 U = np.array(U, ndmin=1)
3959 Y = self._ppf(U)[0]
3960 else:
3961 Y = self._ppf(U)
3962 return Y
3964 def _entropy(self):
3965 return stats.entropy(self.pk)
3967 def generic_moment(self, n):
3968 n = asarray(n)
3969 return np.sum(self.xk**n[np.newaxis, ...] * self.pk, axis=0)
3971 def _expect(self, fun, lb, ub, *args, **kwds):
3972 # ignore all args, just do a brute force summation
3973 supp = self.xk[(lb <= self.xk) & (self.xk <= ub)]
3974 vals = fun(supp)
3975 return np.sum(vals)
3978def _check_shape(argshape, size):
3979 """
3980 This is a utility function used by `_rvs()` in the class geninvgauss_gen.
3981 It compares the tuple argshape to the tuple size.
3983 Parameters
3984 ----------
3985 argshape : tuple of integers
3986 Shape of the arguments.
3987 size : tuple of integers or integer
3988 Size argument of rvs().
3990 Returns
3991 -------
3992 The function returns two tuples, scalar_shape and bc.
3994 scalar_shape : tuple
3995 Shape to which the 1-d array of random variates returned by
3996 _rvs_scalar() is converted when it is copied into the
3997 output array of _rvs().
3999 bc : tuple of booleans
4000 bc is an tuple the same length as size. bc[j] is True if the data
4001 associated with that index is generated in one call of _rvs_scalar().
4003 """
4004 scalar_shape = []
4005 bc = []
4006 for argdim, sizedim in zip_longest(argshape[::-1], size[::-1],
4007 fillvalue=1):
4008 if sizedim > argdim or (argdim == sizedim == 1):
4009 scalar_shape.append(sizedim)
4010 bc.append(True)
4011 else:
4012 bc.append(False)
4013 return tuple(scalar_shape[::-1]), tuple(bc[::-1])
4016def get_distribution_names(namespace_pairs, rv_base_class):
4017 """Collect names of statistical distributions and their generators.
4019 Parameters
4020 ----------
4021 namespace_pairs : sequence
4022 A snapshot of (name, value) pairs in the namespace of a module.
4023 rv_base_class : class
4024 The base class of random variable generator classes in a module.
4026 Returns
4027 -------
4028 distn_names : list of strings
4029 Names of the statistical distributions.
4030 distn_gen_names : list of strings
4031 Names of the generators of the statistical distributions.
4032 Note that these are not simply the names of the statistical
4033 distributions, with a _gen suffix added.
4035 """
4036 distn_names = []
4037 distn_gen_names = []
4038 for name, value in namespace_pairs:
4039 if name.startswith('_'):
4040 continue
4041 if name.endswith('_gen') and issubclass(value, rv_base_class):
4042 distn_gen_names.append(name)
4043 if isinstance(value, rv_base_class):
4044 distn_names.append(name)
4045 return distn_names, distn_gen_names