|
# # Author: Travis Oliphant 2002-2011 with contributions from # SciPy Developers 2004-2011 #
# for root finding for discrete distribution ppf, and max likelihood estimation
# for functions of continuous distributions (e.g. moments, entropy, cdf)
# to approximate the pdf of a continuous distribution given its cdf
product, reshape, zeros, floor, logical_and, log, sqrt, exp)
asarray, nan, inf, isinf, NINF, empty)
else: instancemethod = types.MethodType
# These are the docstring parts used for substitution in specific # distribution docstrings
'notes': """\nNotes\n-----\n""", 'examples': """\nExamples\n--------\n"""}
rvs(%(shapes)s, loc=0, scale=1, size=1, random_state=None) Random variates. """ pdf(x, %(shapes)s, loc=0, scale=1) Probability density function. """ logpdf(x, %(shapes)s, loc=0, scale=1) Log of the probability density function. """ pmf(k, %(shapes)s, loc=0, scale=1) Probability mass function. """ logpmf(k, %(shapes)s, loc=0, scale=1) Log of the probability mass function. """ cdf(x, %(shapes)s, loc=0, scale=1) Cumulative distribution function. """ logcdf(x, %(shapes)s, loc=0, scale=1) Log of the cumulative distribution function. """ sf(x, %(shapes)s, loc=0, scale=1) Survival function (also defined as ``1 - cdf``, but `sf` is sometimes more accurate). """ logsf(x, %(shapes)s, loc=0, scale=1) Log of the survival function. """ ppf(q, %(shapes)s, loc=0, scale=1) Percent point function (inverse of ``cdf`` --- percentiles). """ isf(q, %(shapes)s, loc=0, scale=1) Inverse survival function (inverse of ``sf``). """ moment(n, %(shapes)s, loc=0, scale=1) Non-central moment of order n """ stats(%(shapes)s, loc=0, scale=1, moments='mv') Mean('m'), variance('v'), skew('s'), and/or kurtosis('k'). """ entropy(%(shapes)s, loc=0, scale=1) (Differential) entropy of the RV. """ fit(data, %(shapes)s, loc=0, scale=1) Parameter estimates for generic data. """ expect(func, args=(%(shapes_)s), loc=0, scale=1, lb=None, ub=None, conditional=False, **kwds) Expected value of a function (of one argument) with respect to the distribution. """ expect(func, args=(%(shapes_)s), loc=0, lb=None, ub=None, conditional=False) Expected value of a function (of one argument) with respect to the distribution. """ median(%(shapes)s, loc=0, scale=1) Median of the distribution. """ mean(%(shapes)s, loc=0, scale=1) Mean of the distribution. """ var(%(shapes)s, loc=0, scale=1) Variance of the distribution. """ std(%(shapes)s, loc=0, scale=1) Standard deviation of the distribution. """ interval(alpha, %(shapes)s, loc=0, scale=1) Endpoints of the range that contains alpha percent of the distribution """ _doc_logpdf, _doc_cdf, _doc_logcdf, _doc_sf, _doc_logsf, _doc_ppf, _doc_isf, _doc_moment, _doc_stats, _doc_entropy, _doc_fit, _doc_expect, _doc_median, _doc_mean, _doc_var, _doc_std, _doc_interval])
As an instance of the `rv_continuous` class, `%(name)s` object inherits from it a collection of generic methods (see below for the full list), and completes them with details specific for this particular distribution. """
Alternatively, the object may be called (as a function) to fix the shape, location, and scale parameters returning a "frozen" continuous RV object:
rv = %(name)s(%(shapes)s, loc=0, scale=1) - Frozen RV object with the same methods but holding the given shape, location, and scale fixed. """ Examples -------- >>> from scipy.stats import %(name)s >>> import matplotlib.pyplot as plt >>> fig, ax = plt.subplots(1, 1)
Calculate a few first moments:
%(set_vals_stmt)s >>> mean, var, skew, kurt = %(name)s.stats(%(shapes)s, moments='mvsk')
Display the probability density function (``pdf``):
>>> x = np.linspace(%(name)s.ppf(0.01, %(shapes)s), ... %(name)s.ppf(0.99, %(shapes)s), 100) >>> ax.plot(x, %(name)s.pdf(x, %(shapes)s), ... 'r-', lw=5, alpha=0.6, label='%(name)s pdf')
Alternatively, the distribution object can be called (as a function) to fix the shape, location and scale parameters. This returns a "frozen" RV object holding the given parameters fixed.
Freeze the distribution and display the frozen ``pdf``:
>>> rv = %(name)s(%(shapes)s) >>> ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')
Check accuracy of ``cdf`` and ``ppf``:
>>> vals = %(name)s.ppf([0.001, 0.5, 0.999], %(shapes)s) >>> np.allclose([0.001, 0.5, 0.999], %(name)s.cdf(vals, %(shapes)s)) True
Generate random numbers:
>>> r = %(name)s.rvs(%(shapes)s, size=1000)
And compare the histogram:
>>> ax.hist(r, density=True, histtype='stepfilled', alpha=0.2) >>> ax.legend(loc='best', frameon=False) >>> plt.show()
"""
The probability density above is defined in the "standardized" form. To shift and/or scale the distribution use the ``loc`` and ``scale`` parameters. Specifically, ``%(name)s.pdf(x, %(shapes)s, loc, scale)`` is identically equivalent to ``%(name)s.pdf(y, %(shapes)s) / scale`` with ``y = (x - loc) / scale``. """
_doc_allmethods, '\n', _doc_default_example])
_doc_allmethods])
'rvs': _doc_rvs, 'pdf': _doc_pdf, 'logpdf': _doc_logpdf, 'cdf': _doc_cdf, 'logcdf': _doc_logcdf, 'sf': _doc_sf, 'logsf': _doc_logsf, 'ppf': _doc_ppf, 'isf': _doc_isf, 'stats': _doc_stats, 'entropy': _doc_entropy, 'fit': _doc_fit, 'moment': _doc_moment, 'expect': _doc_expect, 'interval': _doc_interval, 'mean': _doc_mean, 'std': _doc_std, 'var': _doc_var, 'median': _doc_median, 'allmethods': _doc_allmethods, 'longsummary': _doc_default_longsummary, 'frozennote': _doc_default_frozen_note, 'example': _doc_default_example, 'default': _doc_default, 'before_notes': _doc_default_before_notes, 'after_notes': _doc_default_locscale }
# Reuse common content between continuous and discrete docs, change some # minor bits.
'ppf', 'isf', 'stats', 'entropy', 'expect', 'median', 'mean', 'var', 'std', 'interval']
'rv_continuous', 'rv_discrete')
Alternatively, the object may be called (as a function) to fix the shape and location parameters returning a "frozen" discrete RV object:
rv = %(name)s(%(shapes)s, loc=0) - Frozen RV object with the same methods but holding the given shape and location fixed. """
Examples -------- >>> from scipy.stats import %(name)s >>> import matplotlib.pyplot as plt >>> fig, ax = plt.subplots(1, 1)
Calculate a few first moments:
%(set_vals_stmt)s >>> mean, var, skew, kurt = %(name)s.stats(%(shapes)s, moments='mvsk')
Display the probability mass function (``pmf``):
>>> x = np.arange(%(name)s.ppf(0.01, %(shapes)s), ... %(name)s.ppf(0.99, %(shapes)s)) >>> ax.plot(x, %(name)s.pmf(x, %(shapes)s), 'bo', ms=8, label='%(name)s pmf') >>> ax.vlines(x, 0, %(name)s.pmf(x, %(shapes)s), colors='b', lw=5, alpha=0.5)
Alternatively, the distribution object can be called (as a function) to fix the shape and location. This returns a "frozen" RV object holding the given parameters fixed.
Freeze the distribution and display the frozen ``pmf``:
>>> rv = %(name)s(%(shapes)s) >>> ax.vlines(x, 0, rv.pmf(x), colors='k', linestyles='-', lw=1, ... label='frozen pmf') >>> ax.legend(loc='best', frameon=False) >>> plt.show()
Check accuracy of ``cdf`` and ``ppf``:
>>> prob = %(name)s.cdf(x, %(shapes)s) >>> np.allclose(x, %(name)s.ppf(prob, %(shapes)s)) True
Generate random numbers:
>>> r = %(name)s.rvs(%(shapes)s, size=1000) """
The probability mass function above is defined in the "standardized" form. To shift distribution use the ``loc`` parameter. Specifically, ``%(name)s.pmf(k, %(shapes)s, loc)`` is identically equivalent to ``%(name)s.pmf(k - loc, %(shapes)s)``. """
docdict_discrete['allmethods']])
docdict_discrete['allmethods'], docdict_discrete['frozennote'], docdict_discrete['example']])
# clean up all the separate docstring elements, we do not need them anymore # in Python 3, loop variables are not visible after the loop
if mu is None: mu = data.mean() return ((data - mu)**n).mean()
if (n == 0): return 1.0 elif (n == 1): if mu is None: val = moment_func(1, *args) else: val = mu elif (n == 2): if mu2 is None or mu is None: val = moment_func(2, *args) else: val = mu2 + mu*mu elif (n == 3): if g1 is None or mu2 is None or mu is None: val = moment_func(3, *args) else: mu3 = g1 * np.power(mu2, 1.5) # 3rd central moment val = mu3+3*mu*mu2+mu*mu*mu # 3rd non-central moment elif (n == 4): if g1 is None or g2 is None or mu2 is None or mu is None: val = moment_func(4, *args) else: mu4 = (g2+3.0)*(mu2**2.0) # 4th central moment mu3 = g1*np.power(mu2, 1.5) # 3rd central moment val = mu4+4*mu*mu3+6*mu*mu*mu2+mu*mu*mu*mu else: val = moment_func(n, *args)
return val
""" skew is third central moment / variance**(1.5) """ data = np.ravel(data) mu = data.mean() m2 = ((data - mu)**2).mean() m3 = ((data - mu)**3).mean() return m3 / np.power(m2, 1.5)
""" kurtosis is fourth central moment / variance**2 - 3 """ data = np.ravel(data) mu = data.mean() m2 = ((data - mu)**2).mean() m4 = ((data - mu)**4).mean() return m4 / m2**2 - 3
# Frozen RV class
self.args = args self.kwds = kwds
# create a new instance self.dist = dist.__class__(**dist._updated_ctor_param())
# a, b may be set in _argcheck, depending on *args, **kwds. Ouch. shapes, _, _ = self.dist._parse_args(*args, **kwds) self.dist._argcheck(*shapes) self.a, self.b = self.dist.a, self.dist.b
def random_state(self): return self.dist._random_state
def random_state(self, seed): self.dist._random_state = check_random_state(seed)
def pdf(self, x): # raises AttributeError in frozen discrete distribution return self.dist.pdf(x, *self.args, **self.kwds)
return self.dist.logpdf(x, *self.args, **self.kwds)
return self.dist.cdf(x, *self.args, **self.kwds)
return self.dist.logcdf(x, *self.args, **self.kwds)
return self.dist.ppf(q, *self.args, **self.kwds)
return self.dist.isf(q, *self.args, **self.kwds)
kwds = self.kwds.copy() kwds.update({'size': size, 'random_state': random_state}) return self.dist.rvs(*self.args, **kwds)
return self.dist.sf(x, *self.args, **self.kwds)
return self.dist.logsf(x, *self.args, **self.kwds)
kwds = self.kwds.copy() kwds.update({'moments': moments}) return self.dist.stats(*self.args, **kwds)
return self.dist.median(*self.args, **self.kwds)
return self.dist.mean(*self.args, **self.kwds)
return self.dist.var(*self.args, **self.kwds)
return self.dist.std(*self.args, **self.kwds)
return self.dist.moment(n, *self.args, **self.kwds)
return self.dist.entropy(*self.args, **self.kwds)
return self.dist.pmf(k, *self.args, **self.kwds)
return self.dist.logpmf(k, *self.args, **self.kwds)
return self.dist.interval(alpha, *self.args, **self.kwds)
# expect method only accepts shape parameters as positional args # hence convert self.args, self.kwds, also loc/scale # See the .expect method docstrings for the meaning of # other parameters. a, loc, scale = self.dist._parse_args(*self.args, **self.kwds) if isinstance(self.dist, rv_discrete): return self.dist.expect(func, a, loc, lb, ub, conditional, **kwds) else: return self.dist.expect(func, a, loc, scale, lb, ub, conditional, **kwds)
# This should be rewritten """Return the sequence of ravel(args[i]) where ravel(condition) is True in 1D.
Examples -------- >>> import numpy as np >>> rand = np.random.random_sample >>> A = rand((4, 5)) >>> B = 2 >>> C = rand((1, 5)) >>> cond = np.ones(A.shape) >>> [A1, B1, C1] = argsreduce(cond, A, B, C) >>> B1.shape (20,) >>> cond[2,:] = 0 >>> [A2, B2, C2] = argsreduce(cond, A, B, C) >>> B2.shape (15,)
""" newargs = [newargs, ]
def _parse_args(self, %(shape_arg_str)s %(locscale_in)s): return (%(shape_arg_str)s), %(locscale_out)s
def _parse_args_rvs(self, %(shape_arg_str)s %(locscale_in)s, size=None): return self._argcheck_rvs(%(shape_arg_str)s %(locscale_out)s, size=size)
def _parse_args_stats(self, %(shape_arg_str)s %(locscale_in)s, moments='mv'): return (%(shape_arg_str)s), %(locscale_out)s, moments """
# Both the continuous and discrete distributions depend on ncx2. # I think the function name ncx2 is an abbreviation for noncentral chi squared.
# We use (xs**2 + ns**2)/2 = (xs - ns)**2/2 + xs*ns, and include the factor # of exp(-xs*ns) into the ive function to improve numerical stability # at large values of xs. See also `rice.pdf`. df2 = df/2.0 - 1.0 xs, ns = np.sqrt(x), np.sqrt(nc) res = xlogy(df2/2.0, x/nc) - 0.5*(xs - ns)**2 res += np.log(ive(df2, xs*ns) / 2.0) return res
return np.exp(_ncx2_log_pdf(x, df, nc))
return chndtr(x, df, nc)
"""Class which encapsulates common functionality between rv_discrete and rv_continuous.
"""
# figure out if _stats signature has 'moments' keyword ('moments' in sign[0]))
def random_state(self): """ Get or set the RandomState object for generating random variates.
This can be either None or an existing RandomState object.
If None (or np.random), use the RandomState singleton used by np.random. If already a RandomState instance, use it. If an int, use a new RandomState instance seeded with seed.
""" return self._random_state
def random_state(self, seed): self._random_state = check_random_state(seed)
return self._updated_ctor_param(), self._random_state
ctor_param, r = state self.__init__(**ctor_param) self._random_state = r return self
self, meths_to_inspect, locscale_in, locscale_out): """Construct the parser for the shape arguments.
Generates the argument-parsing functions dynamically and attaches them to the instance. Is supposed to be called in __init__ of a class for each distribution.
If self.shapes is a non-empty string, interprets it as a comma-separated list of shape parameters.
Otherwise inspects the call signatures of `meths_to_inspect` and constructs the argument-parsing functions from these. In this case also sets `shapes` and `numargs`. """
# sanitize the user-supplied shapes if not isinstance(self.shapes, string_types): raise TypeError('shapes must be a string.')
shapes = self.shapes.replace(',', ' ').split()
for field in shapes: if keyword.iskeyword(field): raise SyntaxError('keywords cannot be used as shapes.') if not re.match('^[_a-zA-Z][_a-zA-Z0-9]*$', field): raise SyntaxError( 'shapes must be valid python identifiers') else: # find out the call signatures (_pdf, _cdf etc), deduce shape # arguments. Generic methods only have 'self, x', any further args # are shapes.
# *args or **kwargs are not allowed w/automatic shapes raise TypeError( '*args are not allowed w/out explicit shapes') raise TypeError( '**kwds are not allowed w/out explicit shapes') raise TypeError('defaults are not allowed for shapes')
# make sure the signatures are consistent raise TypeError('Shape arguments are inconsistent.') else:
# have the arguments, construct the method from template locscale_in=locscale_in, locscale_out=locscale_out, ) # NB: attach to the instance, not class instancemethod(ns[name], self, self.__class__) )
# allows more general subclassing with *args
"""Construct the instance docstring with string substitutions."""
else:
# remove shapes from call parameters if there are none "\n%(shapes)s : array_like\n shape parameters", "") # necessary because we use %(shapes)s in two forms (w w/o ", ")
# correct for empty shapes
docdict=None, discrete='continuous'): """Construct instance docstring from the default template.""" longname = 'A' extradoc = extradoc[2:] '\n\n%(before_notes)s\n', docheaders['notes'], extradoc, '\n%(example)s'])
"""Freeze the distribution for the given arguments.
Parameters ---------- arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution. Should include all the non-optional arguments, may include ``loc`` and ``scale``.
Returns ------- rv_frozen : rv_frozen instance The frozen distribution.
""" return rv_frozen(self, *args, **kwds)
return self.freeze(*args, **kwds)
# The actual calculation functions (no basic checking need be done) # If these are defined, the others won't be looked at. # Otherwise, the other set can be defined. return None, None, None, None
# Central moments # Silence floating point warnings from integration. olderr = np.seterr(all='ignore') vals = self.generic_moment(n, *args) np.seterr(**olderr) return vals
# Handle broadcasting and size validation of the rvs method. # Subclasses should not have to override this method. # The rule is that if `size` is not None, then `size` gives the # shape of the result (integer values of `size` are treated as # tuples with length 1; i.e. `size=3` is the same as `size=(3,)`.) # # `args` is expected to contain the shape parameters (if any), the # location and the scale in a flat tuple (e.g. if there are two # shape parameters `a` and `b`, `args` will be `(a, b, loc, scale)`). # The only keyword argument expected is 'size'. size = kwargs.get('size', None) all_bcast = np.broadcast_arrays(*args)
def squeeze_left(a): while a.ndim > 0 and a.shape[0] == 1: a = a[0] return a
# Eliminate trivial leading dimensions. In the convention # used by numpy's random variate generators, trivial leading # dimensions are effectively ignored. In other words, when `size` # is given, trivial leading dimensions of the broadcast parameters # in excess of the number of dimensions in size are ignored, e.g. # >>> np.random.normal([[1, 3, 5]], [[[[0.01]]]], size=3) # array([ 1.00104267, 3.00422496, 4.99799278]) # If `size` is not given, the exact broadcast shape is preserved: # >>> np.random.normal([[1, 3, 5]], [[[[0.01]]]]) # array([[[[ 1.00862899, 3.00061431, 4.99867122]]]]) # all_bcast = [squeeze_left(a) for a in all_bcast] bcast_shape = all_bcast[0].shape bcast_ndim = all_bcast[0].ndim
if size is None: size_ = bcast_shape else: size_ = tuple(np.atleast_1d(size))
# Check compatibility of size_ with the broadcast shape of all # the parameters. This check is intended to be consistent with # how the numpy random variate generators (e.g. np.random.normal, # np.random.beta) handle their arguments. The rule is that, if size # is given, it determines the shape of the output. Broadcasting # can't change the output size.
# This is the standard broadcasting convention of extending the # shape with fewer dimensions with enough dimensions of length 1 # so that the two shapes have the same number of dimensions. ndiff = bcast_ndim - len(size_) if ndiff < 0: bcast_shape = (1,)*(-ndiff) + bcast_shape elif ndiff > 0: size_ = (1,)*ndiff + size_
# This compatibility test is not standard. In "regular" broadcasting, # two shapes are compatible if for each dimension, the lengths are the # same or one of the lengths is 1. Here, the length of a dimension in # size_ must not be less than the corresponding length in bcast_shape. ok = all([bcdim == 1 or bcdim == szdim for (bcdim, szdim) in zip(bcast_shape, size_)]) if not ok: raise ValueError("size does not match the broadcast shape of " "the parameters.")
param_bcast = all_bcast[:-2] loc_bcast = all_bcast[-2] scale_bcast = all_bcast[-1]
return param_bcast, loc_bcast, scale_bcast, size_
## These are the methods you must define (standard form functions) ## NB: generic _pdf, _logpdf, _cdf are different for ## rv_continuous and rv_discrete hence are defined in there """Default check for correct values on args and keywords.
Returns condition array of 1's where arguments are correct and 0's where they are not.
"""
return (self.a <= x) & (x <= self.b)
# This method must handle self._size being a tuple, and it must # properly broadcast *args and self._size. self._size might be # an empty tuple, which means a scalar random variate is to be # generated.
## Use basic inverse cdf algorithm for RV generation as default. U = self._random_state.random_sample(self._size) Y = self._ppf(U, *args) return Y
return log(self._cdf(x, *args))
return 1.0-self._cdf(x, *args)
return log(self._sf(x, *args))
return self._ppfvec(q, *args)
return self._ppf(1.0-q, *args) # use correct _ppf for subclasses
# These are actually called, and should not be overwritten if you # want to keep error checking. """ Random variates of given type.
Parameters ---------- arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information). loc : array_like, optional Location parameter (default=0). scale : array_like, optional Scale parameter (default=1). size : int or tuple of ints, optional Defining number of random variates (default is 1). random_state : None or int or ``np.random.RandomState`` instance, optional If int or RandomState, use it for drawing the random variates. If None, rely on ``self.random_state``. Default is None.
Returns ------- rvs : ndarray or scalar Random variates of given `size`.
""" discrete = kwds.pop('discrete', None) rndm = kwds.pop('random_state', None) args, loc, scale, size = self._parse_args_rvs(*args, **kwds) cond = logical_and(self._argcheck(*args), (scale >= 0)) if not np.all(cond): raise ValueError("Domain error in arguments.")
if np.all(scale == 0): return loc*ones(size, 'd')
# extra gymnastics needed for a custom random_state if rndm is not None: random_state_saved = self._random_state self._random_state = check_random_state(rndm)
# `size` should just be an argument to _rvs(), but for, um, # historical reasons, it is made an attribute that is read # by _rvs(). self._size = size vals = self._rvs(*args)
vals = vals * scale + loc
# do not forget to restore the _random_state if rndm is not None: self._random_state = random_state_saved
# Cast to int if discrete if discrete: if size == (): vals = int(vals) else: vals = vals.astype(int)
return vals
""" Some statistics of the given RV.
Parameters ---------- arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional location parameter (default=0) scale : array_like, optional (continuous RVs only) scale parameter (default=1) moments : str, optional composed of letters ['mvsk'] defining which moments to compute: 'm' = mean, 'v' = variance, 's' = (Fisher's) skew, 'k' = (Fisher's) kurtosis. (default is 'mv')
Returns ------- stats : sequence of requested moments.
""" args, loc, scale, moments = self._parse_args_stats(*args, **kwds) # scale = 1 by construction for discrete RVs loc, scale = map(asarray, (loc, scale)) args = tuple(map(asarray, args)) cond = self._argcheck(*args) & (scale > 0) & (loc == loc) output = [] default = valarray(shape(cond), self.badvalue)
# Use only entries that are valid in calculation if np.any(cond): goodargs = argsreduce(cond, *(args+(scale, loc))) scale, loc, goodargs = goodargs[-2], goodargs[-1], goodargs[:-2]
if self._stats_has_moments: mu, mu2, g1, g2 = self._stats(*goodargs, **{'moments': moments}) else: mu, mu2, g1, g2 = self._stats(*goodargs) if g1 is None: mu3 = None else: if mu2 is None: mu2 = self._munp(2, *goodargs) if g2 is None: # (mu2**1.5) breaks down for nan and inf mu3 = g1 * np.power(mu2, 1.5)
if 'm' in moments: if mu is None: mu = self._munp(1, *goodargs) out0 = default.copy() place(out0, cond, mu * scale + loc) output.append(out0)
if 'v' in moments: if mu2 is None: mu2p = self._munp(2, *goodargs) if mu is None: mu = self._munp(1, *goodargs) mu2 = mu2p - mu * mu if np.isinf(mu): # if mean is inf then var is also inf mu2 = np.inf out0 = default.copy() place(out0, cond, mu2 * scale * scale) output.append(out0)
if 's' in moments: if g1 is None: mu3p = self._munp(3, *goodargs) if mu is None: mu = self._munp(1, *goodargs) if mu2 is None: mu2p = self._munp(2, *goodargs) mu2 = mu2p - mu * mu with np.errstate(invalid='ignore'): mu3 = mu3p - 3 * mu * mu2 - mu**3 g1 = mu3 / np.power(mu2, 1.5) out0 = default.copy() place(out0, cond, g1) output.append(out0)
if 'k' in moments: if g2 is None: mu4p = self._munp(4, *goodargs) if mu is None: mu = self._munp(1, *goodargs) if mu2 is None: mu2p = self._munp(2, *goodargs) mu2 = mu2p - mu * mu if mu3 is None: mu3p = self._munp(3, *goodargs) with np.errstate(invalid='ignore'): mu3 = mu3p - 3 * mu * mu2 - mu**3 with np.errstate(invalid='ignore'): mu4 = mu4p - 4 * mu * mu3 - 6 * mu * mu * mu2 - mu**4 g2 = mu4 / mu2**2.0 - 3.0 out0 = default.copy() place(out0, cond, g2) output.append(out0) else: # no valid args output = [] for _ in moments: out0 = default.copy() output.append(out0)
if len(output) == 1: return output[0] else: return tuple(output)
""" Differential entropy of the RV.
Parameters ---------- arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information). loc : array_like, optional Location parameter (default=0). scale : array_like, optional (continuous distributions only). Scale parameter (default=1).
Notes ----- Entropy is defined base `e`:
>>> drv = rv_discrete(values=((0, 1), (0.5, 0.5))) >>> np.allclose(drv.entropy(), np.log(2.0)) True
""" args, loc, scale = self._parse_args(*args, **kwds) # NB: for discrete distributions scale=1 by construction in _parse_args args = tuple(map(asarray, args)) cond0 = self._argcheck(*args) & (scale > 0) & (loc == loc) output = zeros(shape(cond0), 'd') place(output, (1-cond0), self.badvalue) goodargs = argsreduce(cond0, *args) place(output, cond0, self.vecentropy(*goodargs) + log(scale)) return output
""" n-th order non-central moment of distribution.
Parameters ---------- n : int, n >= 1 Order of moment. arg1, arg2, arg3,... : float The shape parameter(s) for the distribution (see docstring of the instance object for more information). loc : array_like, optional location parameter (default=0) scale : array_like, optional scale parameter (default=1)
""" args, loc, scale = self._parse_args(*args, **kwds) if not (self._argcheck(*args) and (scale > 0)): return nan if (floor(n) != n): raise ValueError("Moment must be an integer.") if (n < 0): raise ValueError("Moment must be positive.") mu, mu2, g1, g2 = None, None, None, None if (n > 0) and (n < 5): if self._stats_has_moments: mdict = {'moments': {1: 'm', 2: 'v', 3: 'vs', 4: 'vk'}[n]} else: mdict = {} mu, mu2, g1, g2 = self._stats(*args, **mdict) val = _moment_from_stats(n, mu, mu2, g1, g2, self._munp, args)
# Convert to transformed X = L + S*Y # E[X^n] = E[(L+S*Y)^n] = L^n sum(comb(n, k)*(S/L)^k E[Y^k], k=0...n) if loc == 0: return scale**n * val else: result = 0 fac = float(scale) / float(loc) for k in range(n): valk = _moment_from_stats(k, mu, mu2, g1, g2, self._munp, args) result += comb(n, k, exact=True)*(fac**k) * valk result += fac**n * val return result * loc**n
""" Median of the distribution.
Parameters ---------- arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional Location parameter, Default is 0. scale : array_like, optional Scale parameter, Default is 1.
Returns ------- median : float The median of the distribution.
See Also -------- stats.distributions.rv_discrete.ppf Inverse of the CDF
""" return self.ppf(0.5, *args, **kwds)
""" Mean of the distribution.
Parameters ---------- arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional location parameter (default=0) scale : array_like, optional scale parameter (default=1)
Returns ------- mean : float the mean of the distribution
""" kwds['moments'] = 'm' res = self.stats(*args, **kwds) if isinstance(res, ndarray) and res.ndim == 0: return res[()] return res
""" Variance of the distribution.
Parameters ---------- arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional location parameter (default=0) scale : array_like, optional scale parameter (default=1)
Returns ------- var : float the variance of the distribution
""" kwds['moments'] = 'v' res = self.stats(*args, **kwds) if isinstance(res, ndarray) and res.ndim == 0: return res[()] return res
""" Standard deviation of the distribution.
Parameters ---------- arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional location parameter (default=0) scale : array_like, optional scale parameter (default=1)
Returns ------- std : float standard deviation of the distribution
""" kwds['moments'] = 'v' res = sqrt(self.stats(*args, **kwds)) return res
""" Confidence interval with equal areas around the median.
Parameters ---------- alpha : array_like of float Probability that an rv will be drawn from the returned range. Each value should be in the range [0, 1]. arg1, arg2, ... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information). loc : array_like, optional location parameter, Default is 0. scale : array_like, optional scale parameter, Default is 1.
Returns ------- a, b : ndarray of float end-points of range that contain ``100 * alpha %`` of the rv's possible values.
""" alpha = asarray(alpha) if np.any((alpha > 1) | (alpha < 0)): raise ValueError("alpha must be between 0 and 1 inclusive") q1 = (1.0-alpha)/2 q2 = (1.0+alpha)/2 a = self.ppf(q1, *args, **kwds) b = self.ppf(q2, *args, **kwds) return a, b
## continuous random variables: implement maybe later ## ## hf --- Hazard Function (PDF / SF) ## chf --- Cumulative hazard function (-log(SF)) ## psf --- Probability sparsity function (reciprocal of the pdf) in ## units of percent-point-function (as a function of q). ## Also, the derivative of the percent-point function.
""" A generic continuous random variable class meant for subclassing.
`rv_continuous` is a base class to construct specific distribution classes and instances for continuous random variables. It cannot be used directly as a distribution.
Parameters ---------- momtype : int, optional The type of generic moment calculation to use: 0 for pdf, 1 (default) for ppf. a : float, optional Lower bound of the support of the distribution, default is minus infinity. b : float, optional Upper bound of the support of the distribution, default is plus infinity. xtol : float, optional The tolerance for fixed point calculation for generic ppf. badvalue : float, optional The value in a result arrays that indicates a value that for which some argument restriction is violated, default is np.nan. name : str, optional The name of the instance. This string is used to construct the default example for distributions. longname : str, optional This string is used as part of the first line of the docstring returned when a subclass has no docstring of its own. Note: `longname` exists for backwards compatibility, do not use for new subclasses. shapes : str, optional The shape of the distribution. For example ``"m, n"`` for a distribution that takes two integers as the two shape arguments for all its methods. If not provided, shape parameters will be inferred from the signature of the private methods, ``_pdf`` and ``_cdf`` of the instance. extradoc : str, optional, deprecated This string is used as the last part of the docstring returned when a subclass has no docstring of its own. Note: `extradoc` exists for backwards compatibility, do not use for new subclasses. seed : None or int or ``numpy.random.RandomState`` instance, optional This parameter defines the RandomState object to use for drawing random variates. If None (or np.random), the global np.random state is used. If integer, it is used to seed the local RandomState instance. Default is None.
Methods ------- rvs logpdf cdf logcdf sf logsf ppf isf moment stats entropy expect median mean std var interval __call__ fit fit_loc_scale nnlf
Notes ----- Public methods of an instance of a distribution class (e.g., ``pdf``, ``cdf``) check their arguments and pass valid arguments to private, computational methods (``_pdf``, ``_cdf``). For ``pdf(x)``, ``x`` is valid if it is within the support of a distribution, ``self.a <= x <= self.b``. Whether a shape parameter is valid is decided by an ``_argcheck`` method (which defaults to checking that its arguments are strictly positive.)
**Subclassing**
New random variables can be defined by subclassing the `rv_continuous` class and re-defining at least the ``_pdf`` or the ``_cdf`` method (normalized to location 0 and scale 1).
If positive argument checking is not correct for your RV then you will also need to re-define the ``_argcheck`` method.
Correct, but potentially slow defaults exist for the remaining methods but for speed and/or accuracy you can over-ride::
_logpdf, _cdf, _logcdf, _ppf, _rvs, _isf, _sf, _logsf
Rarely would you override ``_isf``, ``_sf`` or ``_logsf``, but you could.
**Methods that can be overwritten by subclasses** ::
_rvs _pdf _cdf _sf _ppf _isf _stats _munp _entropy _argcheck
There are additional (internal and private) generic methods that can be useful for cross-checking and for debugging, but might work in all cases when directly called.
A note on ``shapes``: subclasses need not specify them explicitly. In this case, `shapes` will be automatically deduced from the signatures of the overridden methods (`pdf`, `cdf` etc). If, for some reason, you prefer to avoid relying on introspection, you can specify ``shapes`` explicitly as an argument to the instance constructor.
**Frozen Distributions**
Normally, you must provide shape parameters (and, optionally, location and scale parameters to each call of a method of a distribution.
Alternatively, the object may be called (as a function) to fix the shape, location, and scale parameters returning a "frozen" continuous RV object:
rv = generic(<shape(s)>, loc=0, scale=1) frozen RV object with the same methods but holding the given shape, location, and scale fixed
**Statistics**
Statistics are computed using numerical integration by default. For speed you can redefine this using ``_stats``:
- take shape parameters and return mu, mu2, g1, g2 - If you can't compute one of these, return it as None - Can also be defined with a keyword argument ``moments``, which is a string composed of "m", "v", "s", and/or "k". Only the components appearing in string should be computed and returned in the order "m", "v", "s", or "k" with missing values returned as None.
Alternatively, you can override ``_munp``, which takes ``n`` and shape parameters and returns the n-th non-central moment of the distribution.
Examples -------- To create a new Gaussian distribution, we would do the following:
>>> from scipy.stats import rv_continuous >>> class gaussian_gen(rv_continuous): ... "Gaussian distribution" ... def _pdf(self, x): ... return np.exp(-x**2 / 2.) / np.sqrt(2.0 * np.pi) >>> gaussian = gaussian_gen(name='gaussian')
``scipy.stats`` distributions are *instances*, so here we subclass `rv_continuous` and create an instance. With this, we now have a fully functional distribution with all relevant methods automagically generated by the framework.
Note that above we defined a standard normal distribution, with zero mean and unit variance. Shifting and scaling of the distribution can be done by using ``loc`` and ``scale`` parameters: ``gaussian.pdf(x, loc, scale)`` essentially computes ``y = (x - loc) / scale`` and ``gaussian._pdf(y) / scale``.
""" badvalue=None, name=None, longname=None, shapes=None, extradoc=None, seed=None):
# save the ctor parameters, cf generic freeze momtype=momtype, a=a, b=b, xtol=xtol, badvalue=badvalue, name=name, longname=longname, shapes=shapes, extradoc=extradoc, seed=seed)
name = 'Distribution' locscale_in='loc=0, scale=1', locscale_out='loc, scale')
# nin correction
self.generic_moment = vectorize(self._mom0_sc, otypes='d') else: # Because of the *args argument of _mom0_sc, vectorize cannot count the # number of arguments correctly.
hstr = "An " else:
# Skip adding docstrings if interpreter is run with -OO extradoc=extradoc, docdict=docdict, discrete='continuous') else:
""" Return the current version of _ctor_param, possibly updated by user.
Used by freezing and pickling. Keep this in sync with the signature of __init__. """ dct = self._ctor_param.copy() dct['a'] = self.a dct['b'] = self.b dct['xtol'] = self.xtol dct['badvalue'] = self.badvalue dct['name'] = self.name dct['shapes'] = self.shapes dct['extradoc'] = self.extradoc return dct
return self.cdf(*(x, )+args)-q
left = right = None if self.a > -np.inf: left = self.a if self.b < np.inf: right = self.b
factor = 10. if not left: # i.e. self.a = -inf left = -1.*factor while self._ppf_to_solve(left, q, *args) > 0.: right = left left *= factor # left is now such that cdf(left) < q if not right: # i.e. self.b = inf right = factor while self._ppf_to_solve(right, q, *args) < 0.: left = right right *= factor # right is now such that cdf(right) > q
return optimize.brentq(self._ppf_to_solve, left, right, args=(q,)+args, xtol=self.xtol)
# moment from definition return x**m * self.pdf(x, *args)
return integrate.quad(self._mom_integ0, self.a, self.b, args=(m,)+args)[0]
# moment calculated using ppf return (self.ppf(q, *args))**m
return integrate.quad(self._mom_integ1, 0, 1, args=(m,)+args)[0]
return derivative(self._cdf, x, dx=1e-5, args=args, order=5)
## Could also define any of these return log(self._pdf(x, *args))
return integrate.quad(self._pdf, self.a, x, args=args)[0]
return self._cdfvec(x, *args)
## generic _argcheck, _logcdf, _sf, _logsf, _ppf, _isf, _rvs are defined ## in rv_generic
""" Probability density function at x of the given RV.
Parameters ---------- x : array_like quantiles arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional location parameter (default=0) scale : array_like, optional scale parameter (default=1)
Returns ------- pdf : ndarray Probability density function evaluated at x
""" args, loc, scale = self._parse_args(*args, **kwds) x, loc, scale = map(asarray, (x, loc, scale)) args = tuple(map(asarray, args)) dtyp = np.find_common_type([x.dtype, np.float64], []) x = np.asarray((x - loc)/scale, dtype=dtyp) cond0 = self._argcheck(*args) & (scale > 0) cond1 = self._support_mask(x) & (scale > 0) cond = cond0 & cond1 output = zeros(shape(cond), dtyp) putmask(output, (1-cond0)+np.isnan(x), self.badvalue) if np.any(cond): goodargs = argsreduce(cond, *((x,)+args+(scale,))) scale, goodargs = goodargs[-1], goodargs[:-1] place(output, cond, self._pdf(*goodargs) / scale) if output.ndim == 0: return output[()] return output
""" Log of the probability density function at x of the given RV.
This uses a more numerically accurate calculation if available.
Parameters ---------- x : array_like quantiles arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional location parameter (default=0) scale : array_like, optional scale parameter (default=1)
Returns ------- logpdf : array_like Log of the probability density function evaluated at x
""" args, loc, scale = self._parse_args(*args, **kwds) x, loc, scale = map(asarray, (x, loc, scale)) args = tuple(map(asarray, args)) dtyp = np.find_common_type([x.dtype, np.float64], []) x = np.asarray((x - loc)/scale, dtype=dtyp) cond0 = self._argcheck(*args) & (scale > 0) cond1 = self._support_mask(x) & (scale > 0) cond = cond0 & cond1 output = empty(shape(cond), dtyp) output.fill(NINF) putmask(output, (1-cond0)+np.isnan(x), self.badvalue) if np.any(cond): goodargs = argsreduce(cond, *((x,)+args+(scale,))) scale, goodargs = goodargs[-1], goodargs[:-1] place(output, cond, self._logpdf(*goodargs) - log(scale)) if output.ndim == 0: return output[()] return output
""" Cumulative distribution function of the given RV.
Parameters ---------- x : array_like quantiles arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional location parameter (default=0) scale : array_like, optional scale parameter (default=1)
Returns ------- cdf : ndarray Cumulative distribution function evaluated at `x`
""" args, loc, scale = self._parse_args(*args, **kwds) x, loc, scale = map(asarray, (x, loc, scale)) args = tuple(map(asarray, args)) dtyp = np.find_common_type([x.dtype, np.float64], []) x = np.asarray((x - loc)/scale, dtype=dtyp) cond0 = self._argcheck(*args) & (scale > 0) cond1 = self._open_support_mask(x) & (scale > 0) cond2 = (x >= self.b) & cond0 cond = cond0 & cond1 output = zeros(shape(cond), dtyp) place(output, (1-cond0)+np.isnan(x), self.badvalue) place(output, cond2, 1.0) if np.any(cond): # call only if at least 1 entry goodargs = argsreduce(cond, *((x,)+args)) place(output, cond, self._cdf(*goodargs)) if output.ndim == 0: return output[()] return output
""" Log of the cumulative distribution function at x of the given RV.
Parameters ---------- x : array_like quantiles arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional location parameter (default=0) scale : array_like, optional scale parameter (default=1)
Returns ------- logcdf : array_like Log of the cumulative distribution function evaluated at x
""" args, loc, scale = self._parse_args(*args, **kwds) x, loc, scale = map(asarray, (x, loc, scale)) args = tuple(map(asarray, args)) dtyp = np.find_common_type([x.dtype, np.float64], []) x = np.asarray((x - loc)/scale, dtype=dtyp) cond0 = self._argcheck(*args) & (scale > 0) cond1 = self._open_support_mask(x) & (scale > 0) cond2 = (x >= self.b) & cond0 cond = cond0 & cond1 output = empty(shape(cond), dtyp) output.fill(NINF) place(output, (1-cond0)*(cond1 == cond1)+np.isnan(x), self.badvalue) place(output, cond2, 0.0) if np.any(cond): # call only if at least 1 entry goodargs = argsreduce(cond, *((x,)+args)) place(output, cond, self._logcdf(*goodargs)) if output.ndim == 0: return output[()] return output
""" Survival function (1 - `cdf`) at x of the given RV.
Parameters ---------- x : array_like quantiles arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional location parameter (default=0) scale : array_like, optional scale parameter (default=1)
Returns ------- sf : array_like Survival function evaluated at x
""" return output
""" Log of the survival function of the given RV.
Returns the log of the "survival function," defined as (1 - `cdf`), evaluated at `x`.
Parameters ---------- x : array_like quantiles arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional location parameter (default=0) scale : array_like, optional scale parameter (default=1)
Returns ------- logsf : ndarray Log of the survival function evaluated at `x`.
""" args, loc, scale = self._parse_args(*args, **kwds) x, loc, scale = map(asarray, (x, loc, scale)) args = tuple(map(asarray, args)) dtyp = np.find_common_type([x.dtype, np.float64], []) x = np.asarray((x - loc)/scale, dtype=dtyp) cond0 = self._argcheck(*args) & (scale > 0) cond1 = self._open_support_mask(x) & (scale > 0) cond2 = cond0 & (x <= self.a) cond = cond0 & cond1 output = empty(shape(cond), dtyp) output.fill(NINF) place(output, (1-cond0)+np.isnan(x), self.badvalue) place(output, cond2, 0.0) if np.any(cond): goodargs = argsreduce(cond, *((x,)+args)) place(output, cond, self._logsf(*goodargs)) if output.ndim == 0: return output[()] return output
""" Percent point function (inverse of `cdf`) at q of the given RV.
Parameters ---------- q : array_like lower tail probability arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional location parameter (default=0) scale : array_like, optional scale parameter (default=1)
Returns ------- x : array_like quantile corresponding to the lower tail probability q.
""" args, loc, scale = self._parse_args(*args, **kwds) q, loc, scale = map(asarray, (q, loc, scale)) args = tuple(map(asarray, args)) cond0 = self._argcheck(*args) & (scale > 0) & (loc == loc) cond1 = (0 < q) & (q < 1) cond2 = cond0 & (q == 0) cond3 = cond0 & (q == 1) cond = cond0 & cond1 output = valarray(shape(cond), value=self.badvalue)
lower_bound = self.a * scale + loc upper_bound = self.b * scale + loc place(output, cond2, argsreduce(cond2, lower_bound)[0]) place(output, cond3, argsreduce(cond3, upper_bound)[0])
if np.any(cond): # call only if at least 1 entry goodargs = argsreduce(cond, *((q,)+args+(scale, loc))) scale, loc, goodargs = goodargs[-2], goodargs[-1], goodargs[:-2] place(output, cond, self._ppf(*goodargs) * scale + loc) if output.ndim == 0: return output[()] return output
""" Inverse survival function (inverse of `sf`) at q of the given RV.
Parameters ---------- q : array_like upper tail probability arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional location parameter (default=0) scale : array_like, optional scale parameter (default=1)
Returns ------- x : ndarray or scalar Quantile corresponding to the upper tail probability q.
""" args, loc, scale = self._parse_args(*args, **kwds) q, loc, scale = map(asarray, (q, loc, scale)) args = tuple(map(asarray, args)) cond0 = self._argcheck(*args) & (scale > 0) & (loc == loc) cond1 = (0 < q) & (q < 1) cond2 = cond0 & (q == 1) cond3 = cond0 & (q == 0) cond = cond0 & cond1 output = valarray(shape(cond), value=self.badvalue)
lower_bound = self.a * scale + loc upper_bound = self.b * scale + loc place(output, cond2, argsreduce(cond2, lower_bound)[0]) place(output, cond3, argsreduce(cond3, upper_bound)[0])
if np.any(cond): goodargs = argsreduce(cond, *((q,)+args+(scale, loc))) scale, loc, goodargs = goodargs[-2], goodargs[-1], goodargs[:-2] place(output, cond, self._isf(*goodargs) * scale + loc) if output.ndim == 0: return output[()] return output
return -np.sum(self._logpdf(x, *args), axis=0)
try: loc = theta[-2] scale = theta[-1] args = tuple(theta[:-2]) except IndexError: raise ValueError("Not enough input arguments.") return loc, scale, args
'''Return negative loglikelihood function.
Notes ----- This is ``-sum(log pdf(x, theta), axis=0)`` where `theta` are the parameters (including loc and scale). ''' loc, scale, args = self._unpack_loc_scale(theta) if not self._argcheck(*args) or scale <= 0: return inf x = asarray((x-loc) / scale) n_log_scale = len(x) * log(scale) if np.any(~self._support_mask(x)): return inf return self._nnlf(x, *args) + n_log_scale
cond0 = ~self._support_mask(x) n_bad = np.count_nonzero(cond0, axis=0) if n_bad > 0: x = argsreduce(~cond0, x)[0] logpdf = self._logpdf(x, *args) finite_logpdf = np.isfinite(logpdf) n_bad += np.sum(~finite_logpdf, axis=0) if n_bad > 0: penalty = n_bad * log(_XMAX) * 100 return -np.sum(logpdf[finite_logpdf], axis=0) + penalty return -np.sum(logpdf, axis=0)
''' Return penalized negative loglikelihood function, i.e., - sum (log pdf(x, theta), axis=0) + penalty where theta are the parameters (including loc and scale) ''' loc, scale, args = self._unpack_loc_scale(theta) if not self._argcheck(*args) or scale <= 0: return inf x = asarray((x-loc) / scale) n_log_scale = len(x) * log(scale) return self._nnlf_and_penalty(x, args) + n_log_scale
# return starting point for fit (shape arguments + loc + scale) if args is None: args = (1.0,)*self.numargs loc, scale = self._fit_loc_scale_support(data, *args) return args + (loc, scale)
# Return the (possibly reduced) function to optimize in order to find MLE # estimates for the .fit method # First of all, convert fshapes params to fnum: eg for stats.beta, # shapes='a, b'. To fix `a`, can specify either `f1` or `fa`. # Convert the latter into the former. if self.shapes: shapes = self.shapes.replace(',', ' ').split() for j, s in enumerate(shapes): val = kwds.pop('f' + s, None) or kwds.pop('fix_' + s, None) if val is not None: key = 'f%d' % j if key in kwds: raise ValueError("Duplicate entry for %s." % key) else: kwds[key] = val
args = list(args) Nargs = len(args) fixedn = [] names = ['f%d' % n for n in range(Nargs - 2)] + ['floc', 'fscale'] x0 = [] for n, key in enumerate(names): if key in kwds: fixedn.append(n) args[n] = kwds.pop(key) else: x0.append(args[n])
if len(fixedn) == 0: func = self._penalized_nnlf restore = None else: if len(fixedn) == Nargs: raise ValueError( "All parameters fixed. There is nothing to optimize.")
def restore(args, theta): # Replace with theta for all numbers not in fixedn # This allows the non-fixed values to vary, but # we still call self.nnlf with all parameters. i = 0 for n in range(Nargs): if n not in fixedn: args[n] = theta[i] i += 1 return args
def func(theta, x): newtheta = restore(args[:], theta) return self._penalized_nnlf(newtheta, x)
return x0, func, restore, args
""" Return MLEs for shape (if applicable), location, and scale parameters from data.
MLE stands for Maximum Likelihood Estimate. Starting estimates for the fit are given by input arguments; for any arguments not provided with starting estimates, ``self._fitstart(data)`` is called to generate such.
One can hold some parameters fixed to specific values by passing in keyword arguments ``f0``, ``f1``, ..., ``fn`` (for shape parameters) and ``floc`` and ``fscale`` (for location and scale parameters, respectively).
Parameters ---------- data : array_like Data to use in calculating the MLEs. args : floats, optional Starting value(s) for any shape-characterizing arguments (those not provided will be determined by a call to ``_fitstart(data)``). No default value. kwds : floats, optional Starting values for the location and scale parameters; no default. Special keyword arguments are recognized as holding certain parameters fixed:
- f0...fn : hold respective shape parameters fixed. Alternatively, shape parameters to fix can be specified by name. For example, if ``self.shapes == "a, b"``, ``fa``and ``fix_a`` are equivalent to ``f0``, and ``fb`` and ``fix_b`` are equivalent to ``f1``.
- floc : hold location parameter fixed to specified value.
- fscale : hold scale parameter fixed to specified value.
- optimizer : The optimizer to use. The optimizer must take ``func``, and starting position as the first two arguments, plus ``args`` (for extra arguments to pass to the function to be optimized) and ``disp=0`` to suppress output as keyword arguments.
Returns ------- mle_tuple : tuple of floats MLEs for any shape parameters (if applicable), followed by those for location and scale. For most random variables, shape statistics will be returned, but there are exceptions (e.g. ``norm``).
Notes ----- This fit is computed by maximizing a log-likelihood function, with penalty applied for samples outside of range of the distribution. The returned answer is not guaranteed to be the globally optimal MLE, it may only be locally optimal, or the optimization may fail altogether.
Examples --------
Generate some data to fit: draw random variates from the `beta` distribution
>>> from scipy.stats import beta >>> a, b = 1., 2. >>> x = beta.rvs(a, b, size=1000)
Now we can fit all four parameters (``a``, ``b``, ``loc`` and ``scale``):
>>> a1, b1, loc1, scale1 = beta.fit(x)
We can also use some prior knowledge about the dataset: let's keep ``loc`` and ``scale`` fixed:
>>> a1, b1, loc1, scale1 = beta.fit(x, floc=0, fscale=1) >>> loc1, scale1 (0, 1)
We can also keep shape parameters fixed by using ``f``-keywords. To keep the zero-th shape parameter ``a`` equal 1, use ``f0=1`` or, equivalently, ``fa=1``:
>>> a1, b1, loc1, scale1 = beta.fit(x, fa=1, floc=0, fscale=1) >>> a1 1
Not all distributions return estimates for the shape parameters. ``norm`` for example just returns estimates for location and scale:
>>> from scipy.stats import norm >>> x = norm.rvs(a, b, size=1000, random_state=123) >>> loc1, scale1 = norm.fit(x) >>> loc1, scale1 (0.92087172783841631, 2.0015750750324668) """ Narg = len(args) if Narg > self.numargs: raise TypeError("Too many input arguments.")
start = [None]*2 if (Narg < self.numargs) or not ('loc' in kwds and 'scale' in kwds): # get distribution specific starting locations start = self._fitstart(data) args += start[Narg:-2] loc = kwds.pop('loc', start[-2]) scale = kwds.pop('scale', start[-1]) args += (loc, scale) x0, func, restore, args = self._reduce_func(args, kwds)
optimizer = kwds.pop('optimizer', optimize.fmin) # convert string to function in scipy.optimize if not callable(optimizer) and isinstance(optimizer, string_types): if not optimizer.startswith('fmin_'): optimizer = "fmin_"+optimizer if optimizer == 'fmin_': optimizer = 'fmin' try: optimizer = getattr(optimize, optimizer) except AttributeError: raise ValueError("%s is not a valid optimizer" % optimizer)
# by now kwds must be empty, since everybody took what they needed if kwds: raise TypeError("Unknown arguments: %s." % kwds)
vals = optimizer(func, x0, args=(ravel(data),), disp=0) if restore is not None: vals = restore(args, vals) vals = tuple(vals) return vals
""" Estimate loc and scale parameters from data accounting for support.
Parameters ---------- data : array_like Data to fit. arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information).
Returns ------- Lhat : float Estimated location parameter for the data. Shat : float Estimated scale parameter for the data.
""" data = np.asarray(data)
# Estimate location and scale according to the method of moments. loc_hat, scale_hat = self.fit_loc_scale(data, *args)
# Compute the support according to the shape parameters. self._argcheck(*args) a, b = self.a, self.b support_width = b - a
# If the support is empty then return the moment-based estimates. if support_width <= 0: return loc_hat, scale_hat
# Compute the proposed support according to the loc and scale estimates. a_hat = loc_hat + a * scale_hat b_hat = loc_hat + b * scale_hat
# Use the moment-based estimates if they are compatible with the data. data_a = np.min(data) data_b = np.max(data) if a_hat < data_a and data_b < b_hat: return loc_hat, scale_hat
# Otherwise find other estimates that are compatible with the data. data_width = data_b - data_a rel_margin = 0.1 margin = data_width * rel_margin
# For a finite interval, both the location and scale # should have interesting values. if support_width < np.inf: loc_hat = (data_a - a) - margin scale_hat = (data_width + 2 * margin) / support_width return loc_hat, scale_hat
# For a one-sided interval, use only an interesting location parameter. if a > -np.inf: return (data_a - a) - margin, 1 elif b < np.inf: return (data_b - b) + margin, 1 else: raise RuntimeError
""" Estimate loc and scale parameters from data using 1st and 2nd moments.
Parameters ---------- data : array_like Data to fit. arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information).
Returns ------- Lhat : float Estimated location parameter for the data. Shat : float Estimated scale parameter for the data.
""" mu, mu2 = self.stats(*args, **{'moments': 'mv'}) tmp = asarray(data) muhat = tmp.mean() mu2hat = tmp.var() Shat = sqrt(mu2hat / mu2) Lhat = muhat - Shat*mu if not np.isfinite(Lhat): Lhat = 0 if not (np.isfinite(Shat) and (0 < Shat)): Shat = 1 return Lhat, Shat
def integ(x): val = self._pdf(x, *args) return entr(val)
# upper limit is often inf, so suppress warnings when integrating olderr = np.seterr(over='ignore') h = integrate.quad(integ, self.a, self.b)[0] np.seterr(**olderr)
if not np.isnan(h): return h else: # try with different limits if integration problems low, upp = self.ppf([1e-10, 1. - 1e-10], *args) if np.isinf(self.b): upper = upp else: upper = self.b if np.isinf(self.a): lower = low else: lower = self.a return integrate.quad(integ, lower, upper)[0]
conditional=False, **kwds): """Calculate expected value of a function with respect to the distribution.
The expected value of a function ``f(x)`` with respect to a distribution ``dist`` is defined as::
ubound E[x] = Integral(f(x) * dist.pdf(x)) lbound
Parameters ---------- func : callable, optional Function for which integral is calculated. Takes only one argument. The default is the identity mapping f(x) = x. args : tuple, optional Shape parameters of the distribution. loc : float, optional Location parameter (default=0). scale : float, optional Scale parameter (default=1). lb, ub : scalar, optional Lower and upper bound for integration. Default is set to the support of the distribution. conditional : bool, optional If True, the integral is corrected by the conditional probability of the integration interval. The return value is the expectation of the function, conditional on being in the given interval. Default is False.
Additional keyword arguments are passed to the integration routine.
Returns ------- expect : float The calculated expected value.
Notes ----- The integration behavior of this function is inherited from `integrate.quad`.
""" lockwds = {'loc': loc, 'scale': scale} self._argcheck(*args) if func is None: def fun(x, *args): return x * self.pdf(x, *args, **lockwds) else: def fun(x, *args): return func(x) * self.pdf(x, *args, **lockwds) if lb is None: lb = loc + self.a * scale if ub is None: ub = loc + self.b * scale if conditional: invfac = (self.sf(lb, *args, **lockwds) - self.sf(ub, *args, **lockwds)) else: invfac = 1.0 kwds['args'] = args # Silence floating point warnings from integration. olderr = np.seterr(all='ignore') vals = integrate.quad(fun, lb, ub, **kwds)[0] / invfac np.seterr(**olderr) return vals
# Helpers for the discrete distributions """Non-central moment of discrete distribution.""" def fun(x): return np.power(x, n) * self._pmf(x, *args) return _expect(fun, self.a, self.b, self.ppf(0.5, *args), self.inc)
b = self.b a = self.a if isinf(b): # Be sure ending point is > q b = int(max(100*q, 10)) while 1: if b >= self.b: qb = 1.0 break qb = self._cdf(b, *args) if (qb < q): b += 10 else: break else: qb = 1.0 if isinf(a): # be sure starting point < q a = int(min(-100*q, -10)) while 1: if a <= self.a: qb = 0.0 break qa = self._cdf(a, *args) if (qa > q): a -= 10 else: break else: qa = self._cdf(a, *args)
while 1: if (qa == q): return a if (qb == q): return b if b <= a+1: # testcase: return wrong number at lower index # python -c "from scipy.stats import zipf;print zipf.ppf(0.01, 2)" wrong # python -c "from scipy.stats import zipf;print zipf.ppf([0.01, 0.61, 0.77, 0.83], 2)" # python -c "from scipy.stats import logser;print logser.ppf([0.1, 0.66, 0.86, 0.93], 0.6)" if qa > q: return a else: return b c = int((a+b)/2.0) qc = self._cdf(c, *args) if (qc < q): if a != c: a = c else: raise RuntimeError('updating stopped, endless loop') qa = qc elif (qc > q): if b != c: b = c else: raise RuntimeError('updating stopped, endless loop') qb = qc else: return c
"""Calculate the entropy of a distribution for given probability values.
If only probabilities `pk` are given, the entropy is calculated as ``S = -sum(pk * log(pk), axis=0)``.
If `qk` is not None, then compute the Kullback-Leibler divergence ``S = sum(pk * log(pk / qk), axis=0)``.
This routine will normalize `pk` and `qk` if they don't sum to 1.
Parameters ---------- pk : sequence Defines the (discrete) distribution. ``pk[i]`` is the (possibly unnormalized) probability of event ``i``. qk : sequence, optional Sequence against which the relative entropy is computed. Should be in the same format as `pk`. base : float, optional The logarithmic base to use, defaults to ``e`` (natural logarithm).
Returns ------- S : float The calculated entropy.
""" pk = asarray(pk) pk = 1.0*pk / np.sum(pk, axis=0) if qk is None: vec = entr(pk) else: qk = asarray(qk) if len(qk) != len(pk): raise ValueError("qk and pk must have same length.") qk = 1.0*qk / np.sum(qk, axis=0) vec = rel_entr(pk, qk) S = np.sum(vec, axis=0) if base is not None: S /= log(base) return S
# Must over-ride one of _pmf or _cdf or pass in # x_k, p(x_k) lists in initialization
""" A generic discrete random variable class meant for subclassing.
`rv_discrete` is a base class to construct specific distribution classes and instances for discrete random variables. It can also be used to construct an arbitrary distribution defined by a list of support points and corresponding probabilities.
Parameters ---------- a : float, optional Lower bound of the support of the distribution, default: 0 b : float, optional Upper bound of the support of the distribution, default: plus infinity moment_tol : float, optional The tolerance for the generic calculation of moments. values : tuple of two array_like, optional ``(xk, pk)`` where ``xk`` are integers with non-zero probabilities ``pk`` with ``sum(pk) = 1``. inc : integer, optional Increment for the support of the distribution. Default is 1. (other values have not been tested) badvalue : float, optional The value in a result arrays that indicates a value that for which some argument restriction is violated, default is np.nan. name : str, optional The name of the instance. This string is used to construct the default example for distributions. longname : str, optional This string is used as part of the first line of the docstring returned when a subclass has no docstring of its own. Note: `longname` exists for backwards compatibility, do not use for new subclasses. shapes : str, optional The shape of the distribution. For example "m, n" for a distribution that takes two integers as the two shape arguments for all its methods If not provided, shape parameters will be inferred from the signatures of the private methods, ``_pmf`` and ``_cdf`` of the instance. extradoc : str, optional This string is used as the last part of the docstring returned when a subclass has no docstring of its own. Note: `extradoc` exists for backwards compatibility, do not use for new subclasses. seed : None or int or ``numpy.random.RandomState`` instance, optional This parameter defines the RandomState object to use for drawing random variates. If None, the global np.random state is used. If integer, it is used to seed the local RandomState instance. Default is None.
Methods ------- rvs pmf logpmf cdf logcdf sf logsf ppf isf moment stats entropy expect median mean std var interval __call__
Notes -----
This class is similar to `rv_continuous`. Whether a shape parameter is valid is decided by an ``_argcheck`` method (which defaults to checking that its arguments are strictly positive.) The main differences are:
- the support of the distribution is a set of integers - instead of the probability density function, ``pdf`` (and the corresponding private ``_pdf``), this class defines the *probability mass function*, `pmf` (and the corresponding private ``_pmf``.) - scale parameter is not defined.
To create a new discrete distribution, we would do the following:
>>> from scipy.stats import rv_discrete >>> class poisson_gen(rv_discrete): ... "Poisson distribution" ... def _pmf(self, k, mu): ... return exp(-mu) * mu**k / factorial(k)
and create an instance::
>>> poisson = poisson_gen(name="poisson")
Note that above we defined the Poisson distribution in the standard form. Shifting the distribution can be done by providing the ``loc`` parameter to the methods of the instance. For example, ``poisson.pmf(x, mu, loc)`` delegates the work to ``poisson._pmf(x-loc, mu)``.
**Discrete distributions from a list of probabilities**
Alternatively, you can construct an arbitrary discrete rv defined on a finite set of values ``xk`` with ``Prob{X=xk} = pk`` by using the ``values`` keyword argument to the `rv_discrete` constructor.
Examples --------
Custom made discrete distribution:
>>> from scipy import stats >>> xk = np.arange(7) >>> pk = (0.1, 0.2, 0.3, 0.1, 0.1, 0.0, 0.2) >>> custm = stats.rv_discrete(name='custm', values=(xk, pk)) >>> >>> import matplotlib.pyplot as plt >>> fig, ax = plt.subplots(1, 1) >>> ax.plot(xk, custm.pmf(xk), 'ro', ms=12, mec='r') >>> ax.vlines(xk, 0, custm.pmf(xk), colors='r', lw=4) >>> plt.show()
Random number generation:
>>> R = custm.rvs(size=100)
""" moment_tol=1e-8, values=None, inc=1, longname=None, shapes=None, extradoc=None, seed=None):
# dispatch to a subclass return super(rv_discrete, cls).__new__(rv_sample) else: # business as usual
moment_tol=1e-8, values=None, inc=1, longname=None, shapes=None, extradoc=None, seed=None):
# cf generic freeze a=a, b=b, name=name, badvalue=badvalue, moment_tol=moment_tol, values=values, inc=inc, longname=longname, shapes=shapes, extradoc=extradoc, seed=seed)
raise ValueError("rv_discrete.__init__(..., values != None, ...)")
locscale_in='loc=0', # scale=1 for discrete RVs locscale_out='loc, 1')
# nin correction needs to be after we know numargs # correct nin for generic moment vectorization self, rv_discrete)
# correct nin for ppf vectorization self, rv_discrete)
# now that self.numargs is defined, we can adjust nin
name = 'Distribution'
# generate docstring for subclass instances hstr = "An " else:
# Skip adding docstrings if interpreter is run with -OO self._construct_default_doc(longname=longname, extradoc=extradoc, docdict=docdict_discrete, discrete='discrete') else:
# discrete RV do not have the scale parameter, remove it '\n scale : array_like, ' 'optional\n scale parameter (default=1)', '')
""" Return the current version of _ctor_param, possibly updated by user.
Used by freezing and pickling. Keep this in sync with the signature of __init__. """ dct = self._ctor_param.copy() dct['a'] = self.a dct['b'] = self.b dct['badvalue'] = self.badvalue dct['moment_tol'] = self.moment_tol dct['inc'] = self.inc dct['name'] = self.name dct['shapes'] = self.shapes dct['extradoc'] = self.extradoc return dct
return floor(k) == k
return self._cdf(k, *args) - self._cdf(k-1, *args)
return log(self._pmf(k, *args))
m = arange(int(self.a), k+1) return np.sum(self._pmf(m, *args), axis=0)
k = floor(x) return self._cdfvec(k, *args)
# generic _logcdf, _sf, _logsf, _ppf, _isf, _rvs defined in rv_generic
""" Random variates of given type.
Parameters ---------- arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information). loc : array_like, optional Location parameter (default=0). size : int or tuple of ints, optional Defining number of random variates (Default is 1). Note that `size` has to be given as keyword, not as positional argument. random_state : None or int or ``np.random.RandomState`` instance, optional If int or RandomState, use it for drawing the random variates. If None, rely on ``self.random_state``. Default is None.
Returns ------- rvs : ndarray or scalar Random variates of given `size`.
""" kwargs['discrete'] = True return super(rv_discrete, self).rvs(*args, **kwargs)
""" Probability mass function at k of the given RV.
Parameters ---------- k : array_like Quantiles. arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) loc : array_like, optional Location parameter (default=0).
Returns ------- pmf : array_like Probability mass function evaluated at k
""" args, loc, _ = self._parse_args(*args, **kwds) k, loc = map(asarray, (k, loc)) args = tuple(map(asarray, args)) k = asarray((k-loc)) cond0 = self._argcheck(*args) cond1 = (k >= self.a) & (k <= self.b) & self._nonzero(k, *args) cond = cond0 & cond1 output = zeros(shape(cond), 'd') place(output, (1-cond0) + np.isnan(k), self.badvalue) if np.any(cond): goodargs = argsreduce(cond, *((k,)+args)) place(output, cond, np.clip(self._pmf(*goodargs), 0, 1)) if output.ndim == 0: return output[()] return output
""" Log of the probability mass function at k of the given RV.
Parameters ---------- k : array_like Quantiles. arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information). loc : array_like, optional Location parameter. Default is 0.
Returns ------- logpmf : array_like Log of the probability mass function evaluated at k.
""" args, loc, _ = self._parse_args(*args, **kwds) k, loc = map(asarray, (k, loc)) args = tuple(map(asarray, args)) k = asarray((k-loc)) cond0 = self._argcheck(*args) cond1 = (k >= self.a) & (k <= self.b) & self._nonzero(k, *args) cond = cond0 & cond1 output = empty(shape(cond), 'd') output.fill(NINF) place(output, (1-cond0) + np.isnan(k), self.badvalue) if np.any(cond): goodargs = argsreduce(cond, *((k,)+args)) place(output, cond, self._logpmf(*goodargs)) if output.ndim == 0: return output[()] return output
""" Cumulative distribution function of the given RV.
Parameters ---------- k : array_like, int Quantiles. arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information). loc : array_like, optional Location parameter (default=0).
Returns ------- cdf : ndarray Cumulative distribution function evaluated at `k`.
""" args, loc, _ = self._parse_args(*args, **kwds) k, loc = map(asarray, (k, loc)) args = tuple(map(asarray, args)) k = asarray((k-loc)) cond0 = self._argcheck(*args) cond1 = (k >= self.a) & (k < self.b) cond2 = (k >= self.b) cond = cond0 & cond1 output = zeros(shape(cond), 'd') place(output, (1-cond0) + np.isnan(k), self.badvalue) place(output, cond2*(cond0 == cond0), 1.0)
if np.any(cond): goodargs = argsreduce(cond, *((k,)+args)) place(output, cond, np.clip(self._cdf(*goodargs), 0, 1)) if output.ndim == 0: return output[()] return output
""" Log of the cumulative distribution function at k of the given RV.
Parameters ---------- k : array_like, int Quantiles. arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information). loc : array_like, optional Location parameter (default=0).
Returns ------- logcdf : array_like Log of the cumulative distribution function evaluated at k.
""" args, loc, _ = self._parse_args(*args, **kwds) k, loc = map(asarray, (k, loc)) args = tuple(map(asarray, args)) k = asarray((k-loc)) cond0 = self._argcheck(*args) cond1 = (k >= self.a) & (k < self.b) cond2 = (k >= self.b) cond = cond0 & cond1 output = empty(shape(cond), 'd') output.fill(NINF) place(output, (1-cond0) + np.isnan(k), self.badvalue) place(output, cond2*(cond0 == cond0), 0.0)
if np.any(cond): goodargs = argsreduce(cond, *((k,)+args)) place(output, cond, self._logcdf(*goodargs)) if output.ndim == 0: return output[()] return output
""" Survival function (1 - `cdf`) at k of the given RV.
Parameters ---------- k : array_like Quantiles. arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information). loc : array_like, optional Location parameter (default=0).
Returns ------- sf : array_like Survival function evaluated at k.
""" args, loc, _ = self._parse_args(*args, **kwds) k, loc = map(asarray, (k, loc)) args = tuple(map(asarray, args)) k = asarray(k-loc) cond0 = self._argcheck(*args) cond1 = (k >= self.a) & (k < self.b) cond2 = (k < self.a) & cond0 cond = cond0 & cond1 output = zeros(shape(cond), 'd') place(output, (1-cond0) + np.isnan(k), self.badvalue) place(output, cond2, 1.0) if np.any(cond): goodargs = argsreduce(cond, *((k,)+args)) place(output, cond, np.clip(self._sf(*goodargs), 0, 1)) if output.ndim == 0: return output[()] return output
""" Log of the survival function of the given RV.
Returns the log of the "survival function," defined as 1 - `cdf`, evaluated at `k`.
Parameters ---------- k : array_like Quantiles. arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information). loc : array_like, optional Location parameter (default=0).
Returns ------- logsf : ndarray Log of the survival function evaluated at `k`.
""" args, loc, _ = self._parse_args(*args, **kwds) k, loc = map(asarray, (k, loc)) args = tuple(map(asarray, args)) k = asarray(k-loc) cond0 = self._argcheck(*args) cond1 = (k >= self.a) & (k < self.b) cond2 = (k < self.a) & cond0 cond = cond0 & cond1 output = empty(shape(cond), 'd') output.fill(NINF) place(output, (1-cond0) + np.isnan(k), self.badvalue) place(output, cond2, 0.0) if np.any(cond): goodargs = argsreduce(cond, *((k,)+args)) place(output, cond, self._logsf(*goodargs)) if output.ndim == 0: return output[()] return output
""" Percent point function (inverse of `cdf`) at q of the given RV.
Parameters ---------- q : array_like Lower tail probability. arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information). loc : array_like, optional Location parameter (default=0).
Returns ------- k : array_like Quantile corresponding to the lower tail probability, q.
""" args, loc, _ = self._parse_args(*args, **kwds) q, loc = map(asarray, (q, loc)) args = tuple(map(asarray, args)) cond0 = self._argcheck(*args) & (loc == loc) cond1 = (q > 0) & (q < 1) cond2 = (q == 1) & cond0 cond = cond0 & cond1 output = valarray(shape(cond), value=self.badvalue, typecode='d') # output type 'd' to handle nin and inf place(output, (q == 0)*(cond == cond), self.a-1) place(output, cond2, self.b) if np.any(cond): goodargs = argsreduce(cond, *((q,)+args+(loc,))) loc, goodargs = goodargs[-1], goodargs[:-1] place(output, cond, self._ppf(*goodargs) + loc)
if output.ndim == 0: return output[()] return output
""" Inverse survival function (inverse of `sf`) at q of the given RV.
Parameters ---------- q : array_like Upper tail probability. arg1, arg2, arg3,... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information). loc : array_like, optional Location parameter (default=0).
Returns ------- k : ndarray or scalar Quantile corresponding to the upper tail probability, q.
""" args, loc, _ = self._parse_args(*args, **kwds) q, loc = map(asarray, (q, loc)) args = tuple(map(asarray, args)) cond0 = self._argcheck(*args) & (loc == loc) cond1 = (q > 0) & (q < 1) cond2 = (q == 1) & cond0 cond = cond0 & cond1
# same problem as with ppf; copied from ppf and changed output = valarray(shape(cond), value=self.badvalue, typecode='d') # output type 'd' to handle nin and inf place(output, (q == 0)*(cond == cond), self.b) place(output, cond2, self.a-1)
# call place only if at least 1 valid argument if np.any(cond): goodargs = argsreduce(cond, *((q,)+args+(loc,))) loc, goodargs = goodargs[-1], goodargs[:-1] # PB same as ticket 766 place(output, cond, self._isf(*goodargs) + loc)
if output.ndim == 0: return output[()] return output
if hasattr(self, 'pk'): return entropy(self.pk) else: return _expect(lambda x: entr(self.pmf(x, *args)), self.a, self.b, self.ppf(0.5, *args), self.inc)
conditional=False, maxcount=1000, tolerance=1e-10, chunksize=32): """ Calculate expected value of a function with respect to the distribution for discrete distribution.
Parameters ---------- func : callable, optional Function for which the expectation value is calculated. Takes only one argument. The default is the identity mapping f(k) = k. args : tuple, optional Shape parameters of the distribution. loc : float, optional Location parameter. Default is 0. lb, ub : int, optional Lower and upper bound for the summation, default is set to the support of the distribution, inclusive (``ul <= k <= ub``). conditional : bool, optional If true then the expectation is corrected by the conditional probability of the summation interval. The return value is the expectation of the function, `func`, conditional on being in the given interval (k such that ``ul <= k <= ub``). Default is False. maxcount : int, optional Maximal number of terms to evaluate (to avoid an endless loop for an infinite sum). Default is 1000. tolerance : float, optional Absolute tolerance for the summation. Default is 1e-10. chunksize : int, optional Iterate over the support of a distributions in chunks of this size. Default is 32.
Returns ------- expect : float Expected value.
Notes ----- For heavy-tailed distributions, the expected value may or may not exist, depending on the function, `func`. If it does exist, but the sum converges slowly, the accuracy of the result may be rather low. For instance, for ``zipf(4)``, accuracy for mean, variance in example is only 1e-5. increasing `maxcount` and/or `chunksize` may improve the result, but may also make zipf very slow.
The function is not vectorized.
""" if func is None: def fun(x): # loc and args from outer scope return (x+loc)*self._pmf(x, *args) else: def fun(x): # loc and args from outer scope return func(x+loc)*self._pmf(x, *args) # used pmf because _pmf does not check support in randint and there # might be problems(?) with correct self.a, self.b at this stage maybe # not anymore, seems to work now with _pmf
self._argcheck(*args) # (re)generate scalar self.a and self.b if lb is None: lb = self.a else: lb = lb - loc # convert bound for standardized distribution if ub is None: ub = self.b else: ub = ub - loc # convert bound for standardized distribution if conditional: invfac = self.sf(lb-1, *args) - self.sf(ub, *args) else: invfac = 1.0
# iterate over the support, starting from the median x0 = self.ppf(0.5, *args) res = _expect(fun, lb, ub, x0, self.inc, maxcount, tolerance, chunksize) return res / invfac
chunksize=32): """Helper for computing the expectation value of `fun`."""
# short-circuit if the support size is small enough if (ub - lb) <= chunksize: supp = np.arange(lb, ub+1, inc) vals = fun(supp) return np.sum(vals)
# otherwise, iterate starting from x0 if x0 < lb: x0 = lb if x0 > ub: x0 = ub
count, tot = 0, 0. # iterate over [x0, ub] inclusive for x in _iter_chunked(x0, ub+1, chunksize=chunksize, inc=inc): count += x.size delta = np.sum(fun(x)) tot += delta if abs(delta) < tolerance * x.size: break if count > maxcount: warnings.warn('expect(): sum did not converge', RuntimeWarning) return tot
# iterate over [lb, x0) for x in _iter_chunked(x0-1, lb-1, chunksize=chunksize, inc=-inc): count += x.size delta = np.sum(fun(x)) tot += delta if abs(delta) < tolerance * x.size: break if count > maxcount: warnings.warn('expect(): sum did not converge', RuntimeWarning) break
return tot
"""Iterate from x0 to x1 in chunks of chunksize and steps inc.
x0 must be finite, x1 need not be. In the latter case, the iterator is infinite. Handles both x0 < x1 and x0 > x1. In the latter case, iterates downwards (make sure to set inc < 0.)
>>> [x for x in _iter_chunked(2, 5, inc=2)] [array([2, 4])] >>> [x for x in _iter_chunked(2, 11, inc=2)] [array([2, 4, 6, 8]), array([10])] >>> [x for x in _iter_chunked(2, -5, inc=-2)] [array([ 2, 0, -2, -4])] >>> [x for x in _iter_chunked(2, -9, inc=-2)] [array([ 2, 0, -2, -4]), array([-6, -8])]
""" if inc == 0: raise ValueError('Cannot increment by zero.') if chunksize <= 0: raise ValueError('Chunk size must be positive; got %s.' % chunksize)
s = 1 if inc > 0 else -1 stepsize = abs(chunksize * inc)
x = x0 while (x - x1) * inc < 0: delta = min(stepsize, abs(x - x1)) step = delta * s supp = np.arange(x, x + step, inc) x += step yield supp
"""A 'sample' discrete distribution defined by the support and values.
The ctor ignores most of the arguments, only needs the `values` argument. """ moment_tol=1e-8, values=None, inc=1, longname=None, shapes=None, extradoc=None, seed=None):
super(rv_discrete, self).__init__(seed)
if values is None: raise ValueError("rv_sample.__init__(..., values=None,...)")
# cf generic freeze self._ctor_param = dict( a=a, b=b, name=name, badvalue=badvalue, moment_tol=moment_tol, values=values, inc=inc, longname=longname, shapes=shapes, extradoc=extradoc, seed=seed)
if badvalue is None: badvalue = nan self.badvalue = badvalue self.moment_tol = moment_tol self.inc = inc self.shapes = shapes self.vecentropy = self._entropy
xk, pk = values
if len(xk) != len(pk): raise ValueError("xk and pk need to have the same length.") if not np.allclose(np.sum(pk), 1): raise ValueError("The sum of provided pk is not 1.")
indx = np.argsort(np.ravel(xk)) self.xk = np.take(np.ravel(xk), indx, 0) self.pk = np.take(np.ravel(pk), indx, 0) self.a = self.xk[0] self.b = self.xk[-1] self.qvals = np.cumsum(self.pk, axis=0)
self.shapes = ' ' # bypass inspection self._construct_argparser(meths_to_inspect=[self._pmf], locscale_in='loc=0', # scale=1 for discrete RVs locscale_out='loc, 1')
self._construct_docstrings(name, longname, extradoc)
return np.select([x == k for k in self.xk], [np.broadcast_arrays(p, x)[0] for p in self.pk], 0)
xx, xxk = np.broadcast_arrays(x[:, None], self.xk) indx = np.argmax(xxk > xx, axis=-1) - 1 return self.qvals[indx]
qq, sqq = np.broadcast_arrays(q[..., None], self.qvals) indx = argmax(sqq >= qq, axis=-1) return self.xk[indx]
# Need to define it explicitly, otherwise .rvs() with size=None # fails due to explicit broadcasting in _ppf U = self._random_state.random_sample(self._size) if self._size is None: U = np.array(U, ndmin=1) Y = self._ppf(U)[0] else: Y = self._ppf(U) return Y
return entropy(self.pk)
n = asarray(n) return np.sum(self.xk**n[np.newaxis, ...] * self.pk, axis=0)
""" Collect names of statistical distributions and their generators.
Parameters ---------- namespace_pairs : sequence A snapshot of (name, value) pairs in the namespace of a module. rv_base_class : class The base class of random variable generator classes in a module.
Returns ------- distn_names : list of strings Names of the statistical distributions. distn_gen_names : list of strings Names of the generators of the statistical distributions. Note that these are not simply the names of the statistical distributions, with a _gen suffix added.
""" |