From 1531cb72f50b2807fac1e6a3efaa84f00ef9bf08 Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Thu, 5 Jan 2023 10:04:11 +0300 Subject: [PATCH] Drop the stats module (separate package?) --- diofant/printing/latex.py | 3 - diofant/printing/pretty/pretty.py | 5 - diofant/printing/str.py | 3 - diofant/stats/__init__.py | 69 - diofant/stats/crv.py | 401 ---- diofant/stats/crv_types.py | 2438 --------------------- diofant/stats/drv.py | 80 - diofant/stats/drv_types.py | 134 -- diofant/stats/frv.py | 349 --- diofant/stats/frv_types.py | 330 --- diofant/stats/rv.py | 1032 --------- diofant/stats/rv_interface.py | 216 -- diofant/tests/core/test_args.py | 363 +-- diofant/tests/printing/test_latex.py | 15 - diofant/tests/printing/test_pretty.py | 14 - diofant/tests/printing/test_str.py | 23 +- diofant/tests/stats/__init__.py | 1 - diofant/tests/stats/test_continuous_rv.py | 665 ------ diofant/tests/stats/test_discrete_rv.py | 38 - diofant/tests/stats/test_finite_rv.py | 302 --- diofant/tests/stats/test_mix.py | 16 - diofant/tests/stats/test_rv.py | 264 --- docs/modules/index.rst | 1 - docs/modules/stats.rst | 182 -- docs/release/notes-0.10.rst | 2 +- docs/release/notes-0.14.rst | 1 + 26 files changed, 8 insertions(+), 6939 deletions(-) delete mode 100644 diofant/stats/__init__.py delete mode 100644 diofant/stats/crv.py delete mode 100644 diofant/stats/crv_types.py delete mode 100644 diofant/stats/drv.py delete mode 100644 diofant/stats/drv_types.py delete mode 100644 diofant/stats/frv.py delete mode 100644 diofant/stats/frv_types.py delete mode 100644 diofant/stats/rv.py delete mode 100644 diofant/stats/rv_interface.py delete mode 100644 diofant/tests/stats/__init__.py delete mode 100644 diofant/tests/stats/test_continuous_rv.py delete mode 100644 diofant/tests/stats/test_discrete_rv.py delete mode 100644 diofant/tests/stats/test_finite_rv.py delete mode 100644 diofant/tests/stats/test_mix.py delete mode 100644 diofant/tests/stats/test_rv.py delete mode 100644 docs/modules/stats.rst diff --git a/diofant/printing/latex.py b/diofant/printing/latex.py index 05adad8699..037e6f8c2a 100644 --- a/diofant/printing/latex.py +++ b/diofant/printing/latex.py @@ -1471,9 +1471,6 @@ def _print_ProductSet(self, p): else: return r' \times '.join(self._print(set) for set in p.sets) - def _print_RandomDomain(self, d): - return 'Domain: ' + self._print(d.as_boolean()) - def _print_FiniteSet(self, s): items = sorted(s.args, key=default_sort_key) return self._print_set(items) diff --git a/diofant/printing/pretty/pretty.py b/diofant/printing/pretty/pretty.py index df7a0e1eb7..6d2cddb808 100644 --- a/diofant/printing/pretty/pretty.py +++ b/diofant/printing/pretty/pretty.py @@ -1604,11 +1604,6 @@ def _print_KroneckerDelta(self, e): bot = stringPict(*a.right(' '*b.width())) return prettyForm(binding=prettyForm.POW, *bot.below(top)) - def _print_RandomDomain(self, d): - pform = self._print('Domain: ') - pform = prettyForm(*pform.right(self._print(d.as_boolean()))) - return pform - def _print_Tr(self, p): # TODO: Handle indices pform = self._print(p.args[0]) diff --git a/diofant/printing/str.py b/diofant/printing/str.py index 3bac8cf8d9..f35bc95a39 100644 --- a/diofant/printing/str.py +++ b/diofant/printing/str.py @@ -113,9 +113,6 @@ def _print_dict(self, d): def _print_Dict(self, expr): return self._print_dict(expr) - def _print_RandomDomain(self, d): - return 'Domain: ' + self._print(d.as_boolean()) - def _print_Dummy(self, expr): return '_' + expr.name diff --git a/diofant/stats/__init__.py b/diofant/stats/__init__.py deleted file mode 100644 index 97f17b3ef6..0000000000 --- a/diofant/stats/__init__.py +++ /dev/null @@ -1,69 +0,0 @@ -""" -Diofant statistics module - -Introduces a random variable type into the Diofant language. - -Random variables may be declared using prebuilt functions such as -Normal, Exponential, Coin, Die, etc... or built with functions like FiniteRV. - -Queries on random expressions can be made using the functions - -========================= ============================= - Expression Meaning -------------------------- ----------------------------- - ``P(condition)`` Probability - ``E(expression)`` Expected value - ``variance(expression)`` Variance - ``density(expression)`` Probability Density Function - ``sample(expression)`` Produce a realization - ``where(condition)`` Where the condition is true -========================= ============================= - -Examples -======== - ->>> from diofant.stats import Die, E, Normal, P, variance ->>> X, Y = Die('X', 6), Die('Y', 6) # Define two six sided dice ->>> Z = Normal('Z', 0, 1) # Declare a Normal random variable with mean 0, std 1 ->>> P(X > 3) # Probability X is greater than 3 -1/2 ->>> E(X + Y) # Expectation of the sum of two dice -7 ->>> variance(X + Y) # Variance of the sum of two dice -35/6 ->>> simplify(P(Z > 1)) # Probability of Z being greater than 1 --erf(sqrt(2)/2)/2 + 1/2 - -""" - -from .crv_types import (Arcsin, Benini, Beta, BetaPrime, Cauchy, Chi, - ChiNoncentral, ChiSquared, ContinuousRV, Dagum, Erlang, - Exponential, FDistribution, FisherZ, Frechet, Gamma, - GammaInverse, Kumaraswamy, Laplace, Logistic, - LogNormal, Maxwell, Nakagami, Normal, Pareto, - QuadraticU, RaisedCosine, Rayleigh, StudentT, - Triangular, Uniform, UniformSum, VonMises, Weibull, - WignerSemicircle) -from .drv_types import Geometric, Poisson -from .frv_types import (Bernoulli, Binomial, Coin, Die, DiscreteUniform, - FiniteRV, Hypergeometric, Rademacher) -from .rv_interface import (E, P, cdf, cmoment, correlation, covariance, - density, dependent, given, independent, moment, - pspace, random_symbols, sample, sample_iter, - sampling_density, skewness, smoment, std, variance, - where) - - -__all__ = ('Arcsin', 'Benini', 'Beta', 'BetaPrime', 'Cauchy', 'Chi', - 'ChiNoncentral', 'ChiSquared', 'ContinuousRV', 'Dagum', 'Erlang', - 'Exponential', 'FDistribution', 'FisherZ', 'Frechet', 'Gamma', - 'GammaInverse', 'Kumaraswamy', 'Laplace', 'Logistic', 'LogNormal', - 'Maxwell', 'Nakagami', 'Normal', 'Pareto', 'QuadraticU', - 'RaisedCosine', 'Rayleigh', 'StudentT', 'Triangular', 'Uniform', - 'UniformSum', 'VonMises', 'Weibull', 'WignerSemicircle', - 'Geometric', 'Poisson', 'Bernoulli', 'Binomial', 'Coin', 'Die', - 'DiscreteUniform', 'FiniteRV', 'Hypergeometric', 'Rademacher', - 'E', 'P', 'cdf', 'cmoment', 'correlation', 'covariance', 'density', - 'dependent', 'given', 'independent', 'moment', 'pspace', - 'random_symbols', 'sample', 'sample_iter', 'sampling_density', - 'skewness', 'smoment', 'std', 'variance', 'where') diff --git a/diofant/stats/crv.py b/diofant/stats/crv.py deleted file mode 100644 index 86fcca767d..0000000000 --- a/diofant/stats/crv.py +++ /dev/null @@ -1,401 +0,0 @@ -""" -Continuous Random Variables Module - -See Also -======== -diofant.stats.crv_types -diofant.stats.rv -diofant.stats.frv - -""" - -import random - -from ..core import (Dummy, Expr, Integer, Lambda, Mul, S, cacheit, oo, symbols, - sympify) -from ..functions import DiracDelta, Piecewise -from ..integrals import Integral, integrate -from ..logic import And -from ..polys.polyerrors import PolynomialError -from ..sets import Interval -from ..solvers import solve -from ..solvers.inequalities import reduce_rational_inequalities -from .rv import (ConditionalDomain, NamedArgsMixin, ProductDomain, - ProductPSpace, PSpace, RandomDomain, SingleDomain, - SinglePSpace, random_symbols) - - -class ContinuousDomain(RandomDomain): - """ - A domain with continuous support - - Represented using symbols and Intervals. - - """ - - is_Continuous = True - - def as_boolean(self): - raise NotImplementedError('Not Implemented for generic Domains') - - -class SingleContinuousDomain(ContinuousDomain, SingleDomain): - """ - A univariate domain with continuous support - - Represented using a single symbol and interval. - - """ - - def integrate(self, expr, variables=None, **kwargs): - # assumes only intervals - ivl = (self.set.left, self.set.right) - return Integral(expr, (self.symbol, *ivl), **kwargs) - - def as_boolean(self): - return self.set.as_relational(self.symbol) - - -class ProductContinuousDomain(ProductDomain, ContinuousDomain): - """A collection of independent domains with continuous support.""" - - def integrate(self, expr, variables=None, **kwargs): - for domain in self.domains: - domain_vars = frozenset(variables) & frozenset(domain.symbols) - if domain_vars: - expr = domain.integrate(expr, domain_vars, **kwargs) - return expr - - def as_boolean(self): - return And(*[domain.as_boolean() for domain in self.domains]) - - -class ConditionalContinuousDomain(ContinuousDomain, ConditionalDomain): - """ - A domain with continuous support that has been further restricted by a - condition such as x > 3 - - """ - - def integrate(self, expr, variables=None, **kwargs): - if variables is None: - variables = self.symbols - if not variables: - return expr - # Extract the full integral - fullintgrl = self.fulldomain.integrate(expr, variables) - # separate into integrand and limits - integrand, limits = fullintgrl.function, list(fullintgrl.limits) - - conditions = [self.condition] - while conditions: - cond = conditions.pop() - if cond.is_Boolean: - if isinstance(cond, And): - conditions.extend(cond.args) - else: - raise NotImplementedError('Or not implemented here') - elif cond.is_Relational: - if cond.is_Equality: - # Add the appropriate Delta to the integrand - integrand *= DiracDelta(cond.lhs - cond.rhs) - else: - symbols = cond.free_symbols & set(self.symbols) - if len(symbols) != 1: # Can't handle x > y - raise NotImplementedError( - 'Multivariate Inequalities not yet implemented') - # Can handle x > 0 - symbol = symbols.pop() - # Find the limit with x, such as (x, -oo, oo) - for i, limit in enumerate(limits): - assert limit[0] == symbol - # Make condition into an Interval like [0, oo] - cintvl = reduce_rational_inequalities_wrap(cond, symbol) - # Make limit into an Interval like [-oo, oo] - lintvl = Interval(limit[1], limit[2]) - # Intersect them to get [0, oo] - intvl = cintvl.intersection(lintvl) - # Put back into limits list - limits[i] = (symbol, intvl.left, intvl.right) - else: - raise TypeError( - f'Condition {cond} is not a relational or Boolean') - - return Integral(integrand, *limits, **kwargs) - - def as_boolean(self): - return And(self.fulldomain.as_boolean(), self.condition) - - @property - def set(self): - if len(self.symbols) == 1: - return (self.fulldomain.set & reduce_rational_inequalities_wrap( - self.condition, tuple(self.symbols)[0])) - else: - raise NotImplementedError( - 'Set of Conditional Domain not Implemented') - - -class ContinuousDistribution(Expr): - """Base class for continuous distributions.""" - - def __call__(self, *args): - return self.pdf(*args) - - -class SingleContinuousDistribution(ContinuousDistribution, NamedArgsMixin): - """Continuous distribution of a single variable. - - Serves as superclass for Normal/Exponential/UniformDistribution etc.... - - Represented by parameters for each of the specific classes. E.g - NormalDistribution is represented by a mean and standard deviation. - - Provides methods for pdf, cdf, and sampling - - See Also: - diofant.stats.crv_types.* - - """ - - set = Interval(-oo, oo, True, True) - - def __new__(cls, *args): - args = list(map(sympify, args)) - return Expr.__new__(cls, *args) - - @staticmethod - def check(*args): - pass - - def sample(self): - """A random realization from the distribution.""" - icdf = self._inverse_cdf_expression() - return icdf(random.uniform(0, 1)) - - @cacheit - def _inverse_cdf_expression(self): - """Inverse of the CDF. - - Used by sample - - """ - x, z = symbols('x, z', extended_real=True, positive=True, cls=Dummy) - # Invert CDF - inverse_cdf = solve(self.cdf(x) - z, x) - if not inverse_cdf or len(inverse_cdf) != 1: - raise NotImplementedError('Could not invert CDF') - - return Lambda(z, inverse_cdf[0][x]) - - @cacheit - def compute_cdf(self, **kwargs): - """Compute the CDF from the PDF - - Returns a Lambda - - """ - x, z = symbols('x, z', real=True, cls=Dummy) - left_bound = self.set.start - - # CDF is integral of PDF from left bound to z - pdf = self.pdf(x) - cdf = integrate(pdf, (x, left_bound, z), **kwargs) - # CDF Ensure that CDF left of left_bound is zero - cdf = Piecewise((cdf, z >= left_bound), (0, True)) - return Lambda(z, cdf) - - def cdf(self, x, **kwargs): - """Cumulative density function.""" - return self.compute_cdf(**kwargs)(x) - - def expectation(self, expr, var, evaluate=True, **kwargs): - """Expectation of expression over distribution.""" - integral = Integral(expr * self.pdf(var), (var, self.set), **kwargs) - return integral.doit() if evaluate else integral - - -class ContinuousDistributionHandmade(SingleContinuousDistribution): - """Continuous distribution with custom pdf and support.""" - - _argnames = 'pdf', - - @property - def set(self): - return self.args[1] - - def __new__(cls, pdf, set=Interval(-oo, oo, True, True)): - return Expr.__new__(cls, pdf, set) - - -class ContinuousPSpace(PSpace): - """Continuous Probability Space. - - Represents the likelihood of an event space defined over a continuum. - - Represented with a ContinuousDomain and a PDF (Lambda-Like) - - """ - - is_Continuous = True - is_real = True - - @property - def domain(self): - return self.args[0] - - @property - def density(self): - return self.args[1] - - @property - def pdf(self): - return self.density(*self.domain.symbols) - - def integrate(self, expr, **kwargs): - rvs = self.values - expr = expr.xreplace({rv: rv.symbol for rv in rvs}) - domain_symbols = frozenset(rv.symbol for rv in rvs) - return self.domain.integrate(self.pdf * expr, - domain_symbols, **kwargs) - - def compute_density(self, expr, **kwargs): - # Common case Density(X) where X in self.values - if expr in self.values: - # Marginalize all other random symbols out of the density - randomsymbols = tuple(set(self.values) - frozenset([expr])) - symbols = tuple(rs.symbol for rs in randomsymbols) - pdf = self.domain.integrate(self.pdf, symbols, **kwargs) - return Lambda(expr.symbol, pdf) - - z = Dummy('z', real=True) - return Lambda(z, self.integrate(DiracDelta(expr - z), **kwargs)) - - @cacheit - def compute_cdf(self, expr, **kwargs): - if not self.domain.set.is_Interval: - raise ValueError( - 'CDF not well defined on multivariate expressions') - - d = self.compute_density(expr, **kwargs) - x, z = symbols('x, z', real=True, cls=Dummy) - left_bound = self.domain.set.start - - # CDF is integral of PDF from left bound to z - cdf = integrate(d(x), (x, left_bound, z), **kwargs) - # CDF Ensure that CDF left of left_bound is zero - cdf = Piecewise((cdf, z >= left_bound), (0, True)) - return Lambda(z, cdf) - - def probability(self, condition, **kwargs): - z = Dummy('z', real=True) - # Univariate case can be handled by where - domain = self.where(condition) - rv = [rv for rv in self.values if rv.symbol == domain.symbol][0] - # Integrate out all other random variables - pdf = self.compute_density(rv, **kwargs) - if domain.set is S.EmptySet: - return Integer(0) - # Integrate out the last variable over the special domain - return Integral(pdf(z), (z, domain.set), **kwargs) - - def where(self, condition): - rvs = frozenset(random_symbols(condition)) - if not (len(rvs) == 1 and rvs.issubset(self.values)): - raise NotImplementedError( - 'Multiple continuous random variables not supported') - rv = tuple(rvs)[0] - interval = reduce_rational_inequalities_wrap(condition, rv) - interval = interval.intersection(self.domain.set) - return SingleContinuousDomain(rv.symbol, interval) - - def conditional_space(self, condition, **kwargs): - condition = condition.xreplace({rv: rv.symbol for rv in self.values}) - - domain = ConditionalContinuousDomain(self.domain, condition) - pdf = self.pdf / domain.integrate(self.pdf, **kwargs) - density = Lambda(domain.symbols, pdf) - - return ContinuousPSpace(domain, density) - - -class SingleContinuousPSpace(ContinuousPSpace, SinglePSpace): - """ - A continuous probability space over a single univariate variable - - These consist of a Symbol and a SingleContinuousDistribution - - This class is normally accessed through the various random variable - functions, Normal, Exponential, Uniform, etc.... - - """ - - @property - def set(self): - return self.distribution.set - - @property - def domain(self): - return SingleContinuousDomain(sympify(self.symbol), self.set) - - def sample(self): - """ - Internal sample method - - Returns dictionary mapping RandomSymbol to realization value. - - """ - return {self.value: self.distribution.sample()} - - def integrate(self, expr, rvs=None, **kwargs): - rvs = rvs or (self.value,) - expr = expr.xreplace({rv: rv.symbol for rv in rvs}) - x = self.value.symbol - return self.distribution.expectation(expr, x, evaluate=False, **kwargs) - - def compute_cdf(self, expr, **kwargs): - if expr == self.value: - return self.distribution.compute_cdf(**kwargs) - else: - raise NotImplementedError - - def compute_density(self, expr, **kwargs): - # https://en.wikipedia.org/wiki/Random_variable#Functions_of_random_variables - if expr == self.value: - return self.density - y = Dummy('y') - gs = solve(expr - y, self.value) - fx = self.compute_density(self.value) - fy = sum(fx(g[self.value]) * abs(g[self.value].diff(y)) for g in gs) - return Lambda(y, fy) - - -class ProductContinuousPSpace(ProductPSpace, ContinuousPSpace): - """A collection of independent continuous probability spaces.""" - - @property - def pdf(self): - p = Mul(*[space.pdf for space in self.spaces]) - return p.subs({rv: rv.symbol for rv in self.values}) - - -def _reduce_inequalities(conditions, var, **kwargs): - try: - return reduce_rational_inequalities(conditions, var, **kwargs) - except PolynomialError as exc: - raise ValueError('Reduction of condition ' - f'failed {conditions[0]}\n') from exc - - -def reduce_rational_inequalities_wrap(condition, var): - if condition.is_Relational: - return _reduce_inequalities([[condition]], var, relational=False) - elif condition.__class__ is And: - intervals = [_reduce_inequalities([[arg]], var, relational=False) - for arg in condition.args] - I = intervals[0] - for i in intervals: - I = I.intersection(i) - return I - else: - raise NotImplementedError diff --git a/diofant/stats/crv_types.py b/diofant/stats/crv_types.py deleted file mode 100644 index b1fdd16a65..0000000000 --- a/diofant/stats/crv_types.py +++ /dev/null @@ -1,2438 +0,0 @@ -""" -Continuous Random Variables - Prebuilt variables - -Contains -======== -Arcsin -Benini -Beta -BetaPrime -Cauchy -Chi -ChiNoncentral -ChiSquared -Dagum -Erlang -Exponential -FDistribution -FisherZ -Frechet -Gamma -GammaInverse -Kumaraswamy -Laplace -Logistic -LogNormal -Maxwell -Nakagami -Normal -Pareto -QuadraticU -RaisedCosine -Rayleigh -StudentT -Triangular -Uniform -UniformSum -VonMises -Weibull -WignerSemicircle -""" - -import random - -from ..concrete import Sum -from ..core import Dummy, Eq, Expr, Lambda, Rational, oo, pi, sympify -from ..functions import Abs, Piecewise, besseli -from ..functions import beta as beta_fn -from ..functions import binomial, cos, exp, factorial, floor, gamma, log, sqrt -from ..logic import And -from ..sets import Interval -from .crv import (ContinuousDistributionHandmade, SingleContinuousDistribution, - SingleContinuousPSpace) -from .rv import _value_check - - -__all__ = ('ContinuousRV', - 'Arcsin', - 'Benini', - 'Beta', - 'BetaPrime', - 'Cauchy', - 'Chi', - 'ChiNoncentral', - 'ChiSquared', - 'Dagum', - 'Erlang', - 'Exponential', - 'FDistribution', - 'FisherZ', - 'Frechet', - 'Gamma', - 'GammaInverse', - 'Kumaraswamy', - 'Laplace', - 'Logistic', - 'LogNormal', - 'Maxwell', - 'Nakagami', - 'Normal', - 'Pareto', - 'QuadraticU', - 'RaisedCosine', - 'Rayleigh', - 'StudentT', - 'Triangular', - 'Uniform', - 'UniformSum', - 'VonMises', - 'Weibull', - 'WignerSemicircle') - - -def ContinuousRV(symbol, density, set=Interval(-oo, oo, True, True)): - """ - Create a Continuous Random Variable given the following: - - -- a symbol - -- a probability density function - -- set on which the pdf is valid (defaults to entire real line) - - Returns a RandomSymbol. - - Many common continuous random variable types are already implemented. - This function should be necessary only very rarely. - - Examples - ======== - - >>> from diofant.stats import E, P - - >>> pdf = sqrt(2)*exp(-x**2/2)/(2*sqrt(pi)) # Normal distribution - >>> X = ContinuousRV(x, pdf) - - >>> E(X) - 0 - >>> P(X > 0) - 1/2 - - """ - pdf = Lambda(symbol, density) - dist = ContinuousDistributionHandmade(pdf, set) - return SingleContinuousPSpace(symbol, dist).value - - -def rv(symbol, cls, args): - args = list(map(sympify, args)) - dist = cls(*args) - dist.check(*args) - return SingleContinuousPSpace(symbol, dist).value - -######################################## -# Continuous Probability Distributions # -######################################## - -# ------------------------------------------------------------------------------ -# Arcsin distribution ---------------------------------------------------------- - - -class ArcsinDistribution(SingleContinuousDistribution): - _argnames = ('a', 'b') - - def pdf(self, x): - return 1/(pi*sqrt((x - self.a)*(self.b - x))) - - -def Arcsin(name, a=0, b=1): - r""" - Create a Continuous Random Variable with an arcsin distribution. - - The density of the arcsin distribution is given by - - .. math:: - f(x) := \frac{1}{\pi\sqrt{(x-a)(b-x)}} - - with `x \in [a,b]`. It must hold that `-\infty < a < b < \infty`. - - Parameters - ========== - - a : Real number, the left interval boundary - b : Real number, the right interval boundary - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> a, b = symbols('a b', real=True) - - >>> X = Arcsin('x', a, b) - >>> density(X)(z) - 1/(pi*sqrt((-a + z)*(b - z))) - - References - ========== - - * https://en.wikipedia.org/wiki/Arcsine_distribution - - """ - return rv(name, ArcsinDistribution, (a, b)) - -# ------------------------------------------------------------------------------ -# Benini distribution ---------------------------------------------------------- - - -class BeniniDistribution(SingleContinuousDistribution): - _argnames = ('alpha', 'beta', 'sigma') - - @property - def set(self): - return Interval(self.sigma, oo, False, True) - - def pdf(self, x): - alpha, beta, sigma = self.alpha, self.beta, self.sigma - return (exp(-alpha*log(x/sigma) - beta*log(x/sigma)**2) - * (alpha/x + 2*beta*log(x/sigma)/x)) - - -def Benini(name, alpha, beta, sigma): - r""" - Create a Continuous Random Variable with a Benini distribution. - - The density of the Benini distribution is given by - - .. math:: - f(x) := e^{-\alpha\log{\frac{x}{\sigma}} - -\beta\log^2\left[{\frac{x}{\sigma}}\right]} - \left(\frac{\alpha}{x}+\frac{2\beta\log{\frac{x}{\sigma}}}{x}\right) - - This is a heavy-tailed distribution and is also known as the log-Rayleigh - distribution. - - Parameters - ========== - - alpha : Real number, `\alpha > 0`, a shape - beta : Real number, `\beta > 0`, a shape - sigma : Real number, `\sigma > 0`, a scale - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> alpha = Symbol('alpha', positive=True) - >>> beta = Symbol('beta', positive=True) - >>> sigma = Symbol('sigma', positive=True) - - >>> X = Benini('x', alpha, beta, sigma) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - / z \ 2/ z \ / / z \\ - - alpha*log|-----| - beta*log |-----| | 2*beta*log|-----|| - \sigma/ \sigma/ |alpha \sigma/| - E *|----- + -----------------| - \ z z / - - References - ========== - - * https://en.wikipedia.org/wiki/Benini_distribution - * https://reference.wolfram.com/legacy/v8/ref/BeniniDistribution.html - - """ - return rv(name, BeniniDistribution, (alpha, beta, sigma)) - -# ------------------------------------------------------------------------------ -# Beta distribution ------------------------------------------------------------ - - -class BetaDistribution(SingleContinuousDistribution): - _argnames = ('alpha', 'beta') - - set = Interval(0, 1) - - @staticmethod - def check(alpha, beta): - _value_check(alpha > 0, 'Alpha must be positive') - _value_check(beta > 0, 'Beta must be positive') - - def pdf(self, x): - alpha, beta = self.alpha, self.beta - return x**(alpha - 1) * (1 - x)**(beta - 1) / beta_fn(alpha, beta) - - def sample(self): - return random.betavariate(self.alpha, self.beta) - - -def Beta(name, alpha, beta): - r""" - Create a Continuous Random Variable with a Beta distribution. - - The density of the Beta distribution is given by - - .. math:: - f(x) := \frac{x^{\alpha-1}(1-x)^{\beta-1}} {\mathrm{B}(\alpha,\beta)} - - with `x \in [0,1]`. - - Parameters - ========== - - alpha : Real number, `\alpha > 0`, a shape - beta : Real number, `\beta > 0`, a shape - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import E, density - - >>> alpha = Symbol('alpha', positive=True) - >>> beta = Symbol('beta', positive=True) - - >>> X = Beta('x', alpha, beta) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - alpha - 1 beta - 1 - z *(-z + 1) - --------------------------- - beta(alpha, beta) - - >>> expand_func(simplify(E(X, meijerg=True))) - alpha/(alpha + beta) - - References - ========== - - * https://en.wikipedia.org/wiki/Beta_distribution - * https://mathworld.wolfram.com/BetaDistribution.html - - """ - return rv(name, BetaDistribution, (alpha, beta)) - -# ------------------------------------------------------------------------------ -# Beta prime distribution ------------------------------------------------------ - - -class BetaPrimeDistribution(SingleContinuousDistribution): - _argnames = ('alpha', 'beta') - - set = Interval(0, oo, False, True) - - def pdf(self, x): - alpha, beta = self.alpha, self.beta - return x**(alpha - 1)*(1 + x)**(-alpha - beta)/beta_fn(alpha, beta) - - -def BetaPrime(name, alpha, beta): - r""" - Create a continuous random variable with a Beta prime distribution. - - The density of the Beta prime distribution is given by - - .. math:: - f(x) := \frac{x^{\alpha-1} (1+x)^{-\alpha -\beta}}{B(\alpha,\beta)} - - with `x > 0`. - - Parameters - ========== - - alpha : Real number, `\alpha > 0`, a shape - beta : Real number, `\beta > 0`, a shape - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> alpha = Symbol('alpha', positive=True) - >>> beta = Symbol('beta', positive=True) - - >>> X = BetaPrime('x', alpha, beta) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - alpha - 1 -alpha - beta - z *(z + 1) - ------------------------------- - beta(alpha, beta) - - References - ========== - - * https://en.wikipedia.org/wiki/Beta_prime_distribution - * https://mathworld.wolfram.com/BetaPrimeDistribution.html - - """ - return rv(name, BetaPrimeDistribution, (alpha, beta)) - -# ------------------------------------------------------------------------------ -# Cauchy distribution ---------------------------------------------------------- - - -class CauchyDistribution(SingleContinuousDistribution): - _argnames = ('x0', 'gamma') - - def pdf(self, x): - return 1/(pi*self.gamma*(1 + ((x - self.x0)/self.gamma)**2)) - - -def Cauchy(name, x0, gamma): - r""" - Create a continuous random variable with a Cauchy distribution. - - The density of the Cauchy distribution is given by - - .. math:: - f(x) := \frac{1}{\pi} \arctan\left(\frac{x-x_0}{\gamma}\right) - +\frac{1}{2} - - Parameters - ========== - - x0 : Real number, the location - gamma : Real number, `\gamma > 0`, the scale - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> x0 = Symbol('x0') - >>> gamma = Symbol('gamma', positive=True) - - >>> X = Cauchy('x', x0, gamma) - - >>> density(X)(z) - 1/(pi*gamma*(1 + (-x0 + z)**2/gamma**2)) - - References - ========== - - * https://en.wikipedia.org/wiki/Cauchy_distribution - * https://mathworld.wolfram.com/CauchyDistribution.html - - """ - return rv(name, CauchyDistribution, (x0, gamma)) - -# ------------------------------------------------------------------------------ -# Chi distribution ------------------------------------------------------------- - - -class ChiDistribution(SingleContinuousDistribution): - _argnames = 'k', - - set = Interval(0, oo, False, True) - - def pdf(self, x): - return 2**(1 - self.k/2)*x**(self.k - 1)*exp(-x**2/2)/gamma(self.k/2) - - -def Chi(name, k): - r""" - Create a continuous random variable with a Chi distribution. - - The density of the Chi distribution is given by - - .. math:: - f(x) := \frac{2^{1-k/2}x^{k-1}e^{-x^2/2}}{\Gamma(k/2)} - - with `x \geq 0`. - - Parameters - ========== - - k : A positive Integer, `k > 0`, the number of degrees of freedom - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import Chi, density - - >>> X = Chi('x', k) - - >>> density(X)(z) - 2**(-k/2 + 1)*E**(-z**2/2)*z**(k - 1)/gamma(k/2) - - References - ========== - - * https://en.wikipedia.org/wiki/Chi_distribution - * https://mathworld.wolfram.com/ChiDistribution.html - - """ - return rv(name, ChiDistribution, (k,)) - -# ------------------------------------------------------------------------------ -# Non-central Chi distribution ------------------------------------------------- - - -class ChiNoncentralDistribution(SingleContinuousDistribution): - _argnames = ('k', 'l') - - set = Interval(0, oo, False, True) - - def pdf(self, x): - k, l = self.k, self.l - return exp(-(x**2+l**2)/2)*x**k*l / (l*x)**(k/2) * besseli(k/2-1, l*x) - - -def ChiNoncentral(name, k, l): - r""" - Create a continuous random variable with a non-central Chi distribution. - - The density of the non-central Chi distribution is given by - - .. math:: - f(x) := \frac{e^{-(x^2+\lambda^2)/2} x^k\lambda} - {(\lambda x)^{k/2}} I_{k/2-1}(\lambda x) - - with `x \geq 0`. Here, `I_\nu (x)` is the - :ref:`modified Bessel function of the first kind `. - - Parameters - ========== - - k : A positive Integer, `k > 0`, the number of degrees of freedom - l : Shift parameter - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> l = Symbol('l') - - >>> X = ChiNoncentral('x', k, l) - - >>> density(X)(z) - E**(-l**2/2 - z**2/2)*l*z**k*(l*z)**(-k/2)*besseli(k/2 - 1, l*z) - - References - ========== - - * https://en.wikipedia.org/wiki/Noncentral_chi_distribution - - """ - return rv(name, ChiNoncentralDistribution, (k, l)) - -# ------------------------------------------------------------------------------ -# Chi squared distribution ----------------------------------------------------- - - -class ChiSquaredDistribution(SingleContinuousDistribution): - _argnames = 'k', - - set = Interval(0, oo, False, True) - - def pdf(self, x): - k = self.k - return 1/(2**(k/2)*gamma(k/2))*x**(k/2 - 1)*exp(-x/2) - - -def ChiSquared(name, k): - r""" - Create a continuous random variable with a Chi-squared distribution. - - The density of the Chi-squared distribution is given by - - .. math:: - f(x) := \frac{1}{2^{\frac{k}{2}}\Gamma\left(\frac{k}{2}\right)} - x^{\frac{k}{2}-1} e^{-\frac{x}{2}} - - with `x \geq 0`. - - Parameters - ========== - - k : A positive Integer, `k > 0`, the number of degrees of freedom - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import E, density, variance - - >>> k = Symbol('k', integer=True, positive=True) - - >>> X = ChiSquared('x', k) - - >>> density(X)(z) - 2**(-k/2)*E**(-z/2)*z**(k/2 - 1)/gamma(k/2) - - >>> combsimp(E(X)) - k - - >>> simplify(expand_func(variance(X))) - 2*k - - References - ========== - - * https://en.wikipedia.org/wiki/Chi_squared_distribution - * https://mathworld.wolfram.com/Chi-SquaredDistribution.html - - """ - return rv(name, ChiSquaredDistribution, (k, )) - -# ------------------------------------------------------------------------------ -# Dagum distribution ----------------------------------------------------------- - - -class DagumDistribution(SingleContinuousDistribution): - _argnames = ('p', 'a', 'b') - - def pdf(self, x): - p, a, b = self.p, self.a, self.b - return a*p/x*((x/b)**(a*p)/(((x/b)**a + 1)**(p + 1))) - - -def Dagum(name, p, a, b): - r""" - Create a continuous random variable with a Dagum distribution. - - The density of the Dagum distribution is given by - - .. math:: - f(x) := \frac{a p}{x} \left( \frac{\left(\tfrac{x}{b}\right)^{a p}} - {\left(\left(\tfrac{x}{b}\right)^a + 1 \right)^{p+1}} \right) - - with `x > 0`. - - Parameters - ========== - - p : Real number, `p > 0`, a shape - a : Real number, `a > 0`, a shape - b : Real number, `b > 0`, a scale - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> p = Symbol('p', positive=True) - >>> b = Symbol('b', positive=True) - >>> a = Symbol('a', positive=True) - - >>> X = Dagum('x', p, a, b) - - >>> density(X)(z) - a*p*(z/b)**(a*p)*((z/b)**a + 1)**(-p - 1)/z - - References - ========== - - * https://en.wikipedia.org/wiki/Dagum_distribution - - """ - return rv(name, DagumDistribution, (p, a, b)) - -# ------------------------------------------------------------------------------ -# Erlang distribution ---------------------------------------------------------- - - -def Erlang(name, k, l): - r""" - Create a continuous random variable with an Erlang distribution. - - The density of the Erlang distribution is given by - - .. math:: - f(x) := \frac{\lambda^k x^{k-1} e^{-\lambda x}}{(k-1)!} - - with `x \in [0,\infty]`. - - Parameters - ========== - - k : Integer - l : Real number, `\lambda > 0`, the rate - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import E, cdf, density, variance - - >>> k = Symbol('k', integer=True, positive=True) - >>> l = Symbol('l', positive=True) - - >>> X = Erlang('x', k, l) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - -l*z k k - 1 - E *l *z - --------------- - gamma(k) - - >>> C = cdf(X, meijerg=True)(z) - >>> pprint(C, use_unicode=False) - / k*lowergamma(k, 0) k*lowergamma(k, l*z) - |- ------------------ + -------------------- for z >= 0 - < gamma(k + 1) gamma(k + 1) - | - \ 0 otherwise - - >>> simplify(E(X)) - k/l - - >>> simplify(variance(X)) - k/l**2 - - References - ========== - - * https://en.wikipedia.org/wiki/Erlang_distribution - * https://mathworld.wolfram.com/ErlangDistribution.html - - """ - return rv(name, GammaDistribution, (k, 1/l)) - -# ------------------------------------------------------------------------------ -# Exponential distribution ----------------------------------------------------- - - -class ExponentialDistribution(SingleContinuousDistribution): - _argnames = 'rate', - - set = Interval(0, oo, False, True) - - @staticmethod - def check(rate): - _value_check(rate > 0, 'Rate must be positive.') - - def pdf(self, x): - return self.rate * exp(-self.rate*x) - - def sample(self): - return random.expovariate(self.rate) - - -def Exponential(name, rate): - r""" - Create a continuous random variable with an Exponential distribution. - - The density of the exponential distribution is given by - - .. math:: - f(x) := \lambda \exp(-\lambda x) - - with `x > 0`. Note that the expected value is `1/\lambda`. - - Parameters - ========== - - rate : A positive Real number, `\lambda > 0`, the rate (or inverse scale/inverse mean) - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import E, cdf, density, skewness, std, variance - - >>> l = Symbol('lambda', positive=True) - - >>> X = Exponential('x', l) - - >>> density(X)(z) - E**(-lambda*z)*lambda - - >>> cdf(X)(z) - Piecewise((1 - E**(-lambda*z), z >= 0), (0, true)) - - >>> E(X) - 1/lambda - - >>> variance(X) - lambda**(-2) - - >>> skewness(X) - 2 - - >>> X = Exponential('x', 10) - - >>> density(X)(z) - 10*E**(-10*z) - - >>> E(X) - 1/10 - - >>> std(X) - 1/10 - - References - ========== - - * https://en.wikipedia.org/wiki/Exponential_distribution - * https://mathworld.wolfram.com/ExponentialDistribution.html - - """ - return rv(name, ExponentialDistribution, (rate, )) - -# ------------------------------------------------------------------------------ -# F distribution --------------------------------------------------------------- - - -class FDistributionDistribution(SingleContinuousDistribution): - _argnames = ('d1', 'd2') - - set = Interval(0, oo, False, True) - - def pdf(self, x): - d1, d2 = self.d1, self.d2 - return (sqrt((d1*x)**d1*d2**d2 / (d1*x+d2)**(d1+d2)) - / (x * beta_fn(d1/2, d2/2))) - - -def FDistribution(name, d1, d2): - r""" - Create a continuous random variable with a F distribution. - - The density of the F distribution is given by - - .. math:: - f(x) := \frac{\sqrt{\frac{(d_1 x)^{d_1} d_2^{d_2}} - {(d_1 x + d_2)^{d_1 + d_2}}}} - {x \mathrm{B} \left(\frac{d_1}{2}, \frac{d_2}{2}\right)} - - with `x > 0`. - - .. TODO - What do these parameters mean? - - Parameters - ========== - - d1 : `d_1 > 0` a parameter - d2 : `d_2 > 0` a parameter - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> d1 = Symbol('d1', positive=True) - >>> d2 = Symbol('d2', positive=True) - - >>> X = FDistribution('x', d1, d2) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - d2 - -- ______________________________ - 2 / d1 -d1 - d2 - d2 *\/ (d1*z) *(d1*z + d2) - -------------------------------------- - /d1 d2\ - z*beta|--, --| - \2 2 / - - References - ========== - - * https://en.wikipedia.org/wiki/F-distribution - * https://mathworld.wolfram.com/F-Distribution.html - - """ - return rv(name, FDistributionDistribution, (d1, d2)) - -# ------------------------------------------------------------------------------ -# Fisher Z distribution -------------------------------------------------------- - - -class FisherZDistribution(SingleContinuousDistribution): - _argnames = ('d1', 'd2') - - def pdf(self, x): - d1, d2 = self.d1, self.d2 - return (2*d1**(d1/2)*d2**(d2/2) / beta_fn(d1/2, d2/2) * - exp(d1*x) / (d1*exp(2*x)+d2)**((d1+d2)/2)) - - -def FisherZ(name, d1, d2): - r""" - Create a Continuous Random Variable with an Fisher's Z distribution. - - The density of the Fisher's Z distribution is given by - - .. math:: - f(x) := \frac{2d_1^{d_1/2} d_2^{d_2/2}} {\mathrm{B}(d_1/2, d_2/2)} - \frac{e^{d_1z}}{\left(d_1e^{2z}+d_2\right)^{\left(d_1+d_2\right)/2}} - - - .. TODO - What is the difference between these degrees of freedom? - - Parameters - ========== - - d1 : `d_1 > 0`, degree of freedom - d2 : `d_2 > 0`, degree of freedom - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> d1 = Symbol('d1', positive=True) - >>> d2 = Symbol('d2', positive=True) - - >>> X = FisherZ('x', d1, d2) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - d1 d2 - d1 d2 - -- - -- - -- -- 2 2 - d1*z 2 2 / 2*z \ - 2*E *d1 *d2 *\E *d1 + d2/ - ----------------------------------------- - /d1 d2\ - beta|--, --| - \2 2 / - - References - ========== - - * https://en.wikipedia.org/wiki/Fisher%27s_z-distribution - * https://mathworld.wolfram.com/Fishersz-Distribution.html - - """ - return rv(name, FisherZDistribution, (d1, d2)) - -# ------------------------------------------------------------------------------ -# Frechet distribution --------------------------------------------------------- - - -class FrechetDistribution(SingleContinuousDistribution): - _argnames = ('a', 's', 'm') - - set = Interval(0, oo, False, True) - - def __new__(cls, a, s=1, m=0): - a, s, m = list(map(sympify, (a, s, m))) - return Expr.__new__(cls, a, s, m) - - def pdf(self, x): - a, s, m = self.a, self.s, self.m - return a/s * ((x-m)/s)**(-1-a) * exp(-((x-m)/s)**(-a)) - - -def Frechet(name, a, s=1, m=0): - r""" - Create a continuous random variable with a Frechet distribution. - - The density of the Frechet distribution is given by - - .. math:: - f(x) := \frac{\alpha}{s} \left(\frac{x-m}{s}\right)^{-1-\alpha} - e^{-(\frac{x-m}{s})^{-\alpha}} - - with `x \geq m`. - - Parameters - ========== - - a : Real number, `a \in \left(0, \infty\right)` the shape - s : Real number, `s \in \left(0, \infty\right)` the scale - m : Real number, `m \in \left(-\infty, \infty\right)` the minimum - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> a, s = symbols('a s', positive=True) - >>> m = Symbol('m', real=True) - - >>> X = Frechet('x', a, s, m) - - >>> density(X)(z) - E**(-((-m + z)/s)**(-a))*a*((-m + z)/s)**(-a - 1)/s - - References - ========== - - * https://en.wikipedia.org/wiki/Fr%C3%A9chet_distribution - - """ - return rv(name, FrechetDistribution, (a, s, m)) - -# ------------------------------------------------------------------------------ -# Gamma distribution ----------------------------------------------------------- - - -class GammaDistribution(SingleContinuousDistribution): - _argnames = ('k', 'theta') - - set = Interval(0, oo, False, True) - - @staticmethod - def check(k, theta): - _value_check(k > 0, 'k must be positive') - _value_check(theta > 0, 'Theta must be positive') - - def pdf(self, x): - k, theta = self.k, self.theta - return x**(k - 1) * exp(-x/theta) / (gamma(k)*theta**k) - - def sample(self): - return random.gammavariate(self.k, self.theta) - - -def Gamma(name, k, theta): - r""" - Create a continuous random variable with a Gamma distribution. - - The density of the Gamma distribution is given by - - .. math:: - f(x) := \frac{1}{\Gamma(k) \theta^k} x^{k - 1} e^{-\frac{x}{\theta}} - - with `x \in [0,1]`. - - Parameters - ========== - - k : Real number, `k > 0`, a shape - theta : Real number, `\theta > 0`, a scale - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import E, cdf, density, variance - - >>> k = Symbol('k', positive=True) - >>> theta = Symbol('theta', positive=True) - - >>> X = Gamma('x', k, theta) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - -z - ----- - theta -k k - 1 - E *theta *z - --------------------- - gamma(k) - - >>> C = cdf(X, meijerg=True)(z) - >>> pprint(C, use_unicode=False) - / / z \ - | k*lowergamma|k, -----| - | k*lowergamma(k, 0) \ theta/ - <- ------------------ + ---------------------- for z >= 0 - | gamma(k + 1) gamma(k + 1) - | - \ 0 otherwise - - >>> E(X) - theta*gamma(k + 1)/gamma(k) - - >>> V = simplify(variance(X)) - >>> pprint(V, use_unicode=False) - 2 - k*theta - - References - ========== - - * https://en.wikipedia.org/wiki/Gamma_distribution - * https://mathworld.wolfram.com/GammaDistribution.html - - """ - return rv(name, GammaDistribution, (k, theta)) - -# ------------------------------------------------------------------------------ -# Inverse Gamma distribution --------------------------------------------------- - - -class GammaInverseDistribution(SingleContinuousDistribution): - _argnames = ('a', 'b') - - set = Interval(0, oo, False, True) - - @staticmethod - def check(a, b): - _value_check(a > 0, 'alpha must be positive') - _value_check(b > 0, 'beta must be positive') - - def pdf(self, x): - a, b = self.a, self.b - return b**a/gamma(a) * x**(-a-1) * exp(-b/x) - - -def GammaInverse(name, a, b): - r""" - Create a continuous random variable with an inverse Gamma distribution. - - The density of the inverse Gamma distribution is given by - - .. math:: - f(x) := \frac{\beta^\alpha}{\Gamma(\alpha)} x^{-\alpha - 1} - \exp\left(\frac{-\beta}{x}\right) - - with `x > 0`. - - Parameters - ========== - - a : Real number, `a > 0` a shape - b : Real number, `b > 0` a scale - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> a = Symbol('a', positive=True) - >>> b = Symbol('b', positive=True) - - >>> X = GammaInverse('x', a, b) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - -b - --- - z a -a - 1 - E *b *z - --------------- - gamma(a) - - References - ========== - - * https://en.wikipedia.org/wiki/Inverse-gamma_distribution - - """ - return rv(name, GammaInverseDistribution, (a, b)) - -# ------------------------------------------------------------------------------ -# Kumaraswamy distribution ----------------------------------------------------- - - -class KumaraswamyDistribution(SingleContinuousDistribution): - _argnames = ('a', 'b') - - set = Interval(0, oo, False, True) - - @staticmethod - def check(a, b): - _value_check(a > 0, 'a must be positive') - _value_check(b > 0, 'b must be positive') - - def pdf(self, x): - a, b = self.a, self.b - return a * b * x**(a-1) * (1-x**a)**(b-1) - - -def Kumaraswamy(name, a, b): - r""" - Create a Continuous Random Variable with a Kumaraswamy distribution. - - The density of the Kumaraswamy distribution is given by - - .. math:: - f(x) := a b x^{a-1} (1-x^a)^{b-1} - - with `x \in [0,1]`. - - Parameters - ========== - - a : Real number, `a > 0` a shape - b : Real number, `b > 0` a shape - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> a = Symbol('a', positive=True) - >>> b = Symbol('b', positive=True) - - >>> X = Kumaraswamy('x', a, b) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - b - 1 - a - 1 / a \ - a*b*z *\- z + 1/ - - - References - ========== - - * https://en.wikipedia.org/wiki/Kumaraswamy_distribution - - """ - return rv(name, KumaraswamyDistribution, (a, b)) - -# ------------------------------------------------------------------------------ -# Laplace distribution --------------------------------------------------------- - - -class LaplaceDistribution(SingleContinuousDistribution): - _argnames = ('mu', 'b') - - def pdf(self, x): - mu, b = self.mu, self.b - return 1/(2*b)*exp(-Abs(x - mu)/b) - - -def Laplace(name, mu, b): - r""" - Create a continuous random variable with a Laplace distribution. - - The density of the Laplace distribution is given by - - .. math:: - f(x) := \frac{1}{2 b} \exp \left(-\frac{|x-\mu|}b \right) - - Parameters - ========== - - mu : Real number, the location (mean) - b : Real number, `b > 0`, a scale - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> mu = Symbol('mu') - >>> b = Symbol('b', positive=True) - - >>> X = Laplace('x', mu, b) - - >>> density(X)(z) - E**(-Abs(mu - z)/b)/(2*b) - - References - ========== - - * https://en.wikipedia.org/wiki/Laplace_distribution - * https://mathworld.wolfram.com/LaplaceDistribution.html - - """ - return rv(name, LaplaceDistribution, (mu, b)) - -# ------------------------------------------------------------------------------ -# Logistic distribution -------------------------------------------------------- - - -class LogisticDistribution(SingleContinuousDistribution): - _argnames = ('mu', 's') - - def pdf(self, x): - mu, s = self.mu, self.s - return exp(-(x - mu)/s)/(s*(1 + exp(-(x - mu)/s))**2) - - -def Logistic(name, mu, s): - r""" - Create a continuous random variable with a logistic distribution. - - The density of the logistic distribution is given by - - .. math:: - f(x) := \frac{e^{-(x-\mu)/s}} {s\left(1+e^{-(x-\mu)/s}\right)^2} - - Parameters - ========== - - mu : Real number, the location (mean) - s : Real number, `s > 0` a scale - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> mu = Symbol('mu', real=True) - >>> s = Symbol('s', positive=True) - - >>> X = Logistic('x', mu, s) - - >>> density(X)(z) - E**((mu - z)/s)/(s*(E**((mu - z)/s) + 1)**2) - - References - ========== - - * https://en.wikipedia.org/wiki/Logistic_distribution - * https://mathworld.wolfram.com/LogisticDistribution.html - - """ - return rv(name, LogisticDistribution, (mu, s)) - -# ------------------------------------------------------------------------------ -# Log Normal distribution ------------------------------------------------------ - - -class LogNormalDistribution(SingleContinuousDistribution): - _argnames = ('mean', 'std') - - set = Interval(0, oo, False, True) - - def pdf(self, x): - mean, std = self.mean, self.std - return exp(-(log(x) - mean)**2 / (2*std**2)) / (x*sqrt(2*pi)*std) - - def sample(self): - return random.lognormvariate(self.mean, self.std) - - -def LogNormal(name, mean, std): - r""" - Create a continuous random variable with a log-normal distribution. - - The density of the log-normal distribution is given by - - .. math:: - f(x) := \frac{1}{x\sqrt{2\pi\sigma^2}} - e^{-\frac{\left(\ln x-\mu\right)^2}{2\sigma^2}} - - with `x \geq 0`. - - Parameters - ========== - - mu : Real number, the log-scale - sigma : Real number, `\sigma^2 > 0` a shape - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> mu = Symbol('mu', real=True) - >>> sigma = Symbol('sigma', positive=True) - - >>> X = LogNormal('x', mu, sigma) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - 2 - -(-mu + log(z)) - ----------------- - 2 - ___ 2*sigma - \/ 2 *E - ------------------------ - ____ - 2*\/ pi *sigma*z - - >>> X = LogNormal('x', 0, 1) # Mean 0, standard deviation 1 - - >>> density(X)(z) - sqrt(2)*E**(-log(z)**2/2)/(2*sqrt(pi)*z) - - References - ========== - - * https://en.wikipedia.org/wiki/Lognormal - * https://mathworld.wolfram.com/LogNormalDistribution.html - - """ - return rv(name, LogNormalDistribution, (mean, std)) - -# ------------------------------------------------------------------------------ -# Maxwell distribution --------------------------------------------------------- - - -class MaxwellDistribution(SingleContinuousDistribution): - _argnames = 'a', - - set = Interval(0, oo, False, True) - - def pdf(self, x): - a = self.a - return sqrt(2/pi)*x**2*exp(-x**2/(2*a**2))/a**3 - - -def Maxwell(name, a): - r""" - Create a continuous random variable with a Maxwell distribution. - - The density of the Maxwell distribution is given by - - .. math:: - f(x) := \sqrt{\frac{2}{\pi}} \frac{x^2 e^{-x^2/(2a^2)}}{a^3} - - with `x \geq 0`. - - .. TODO - what does the parameter mean? - - Parameters - ========== - - a : Real number, `a > 0` - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import E, density, variance - - >>> a = Symbol('a', positive=True) - - >>> X = Maxwell('x', a) - - >>> density(X)(z) - sqrt(2)*E**(-z**2/(2*a**2))*z**2/(sqrt(pi)*a**3) - - >>> E(X) - 2*sqrt(2)*a/sqrt(pi) - - >>> simplify(variance(X)) - a**2*(-8 + 3*pi)/pi - - References - ========== - - * https://en.wikipedia.org/wiki/Maxwell_distribution - * https://mathworld.wolfram.com/MaxwellDistribution.html - - """ - return rv(name, MaxwellDistribution, (a, )) - -# ------------------------------------------------------------------------------ -# Nakagami distribution -------------------------------------------------------- - - -class NakagamiDistribution(SingleContinuousDistribution): - _argnames = ('mu', 'omega') - - set = Interval(0, oo, False, True) - - def pdf(self, x): - mu, omega = self.mu, self.omega - return 2*mu**mu/(gamma(mu)*omega**mu)*x**(2*mu - 1)*exp(-mu/omega*x**2) - - -def Nakagami(name, mu, omega): - r""" - Create a continuous random variable with a Nakagami distribution. - - The density of the Nakagami distribution is given by - - .. math:: - f(x) := \frac{2\mu^\mu}{\Gamma(\mu)\omega^\mu} x^{2\mu-1} - \exp\left(-\frac{\mu}{\omega}x^2 \right) - - with `x > 0`. - - Parameters - ========== - - mu : Real number, `\mu \geq \frac{1}{2}` a shape - omega : Real number, `\omega > 0`, the spread - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import E, density, variance - - >>> mu = Symbol('mu', positive=True) - >>> omega = Symbol('omega', positive=True) - - >>> X = Nakagami('x', mu, omega) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - 2 - -mu*z - ------- - omega mu -mu 2*mu - 1 - 2*E *mu *omega *z - ---------------------------------- - gamma(mu) - - >>> simplify(E(X, meijerg=True)) - sqrt(mu)*sqrt(omega)*gamma(mu + 1/2)/gamma(mu + 1) - - >>> V = simplify(variance(X, meijerg=True)) - >>> pprint(V, use_unicode=False) - 2 - omega*gamma (mu + 1/2) - omega - ----------------------- - gamma(mu)*gamma(mu + 1) - - References - ========== - - * https://en.wikipedia.org/wiki/Nakagami_distribution - - """ - return rv(name, NakagamiDistribution, (mu, omega)) - -# ------------------------------------------------------------------------------ -# Normal distribution ---------------------------------------------------------- - - -class NormalDistribution(SingleContinuousDistribution): - _argnames = ('mean', 'std') - - @staticmethod - def check(mean, std): - _value_check(std > 0, 'Standard deviation must be positive') - - def pdf(self, x): - return exp(-(x - self.mean)**2 / (2*self.std**2)) / (sqrt(2*pi)*self.std) - - def sample(self): - return random.normalvariate(self.mean, self.std) - - -def Normal(name, mean, std): - r""" - Create a continuous random variable with a Normal distribution. - - The density of the Normal distribution is given by - - .. math:: - f(x) := \frac{1}{\sigma\sqrt{2\pi}} e^{ -\frac{(x-\mu)^2}{2\sigma^2} } - - Parameters - ========== - - mu : Real number, the mean - sigma : Real number, `\sigma^2 > 0` the variance - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import E, cdf, density, skewness, std - - >>> mu = Symbol('mu') - >>> sigma = Symbol('sigma', positive=True) - - >>> X = Normal('x', mu, sigma) - - >>> density(X)(z) - sqrt(2)*E**(-(-mu + z)**2/(2*sigma**2))/(2*sqrt(pi)*sigma) - - >>> C = simplify(expand(cdf(X)))(z) # it needs a little more help... - >>> pprint(C, use_unicode=False) - / ___ \ - |\/ 2 *(mu - z)| - erf|--------------| - \ 2*sigma / 1 - - ------------------- + - - 2 2 - - >>> simplify(skewness(X)) - 0 - - >>> X = Normal('x', 0, 1) # Mean 0, standard deviation 1 - >>> density(X)(z) - sqrt(2)*E**(-z**2/2)/(2*sqrt(pi)) - - >>> E(2*X + 1) - 1 - - >>> simplify(std(2*X + 1)) - 2 - - References - ========== - - * https://en.wikipedia.org/wiki/Normal_distribution - * https://mathworld.wolfram.com/NormalDistributionFunction.html - - """ - return rv(name, NormalDistribution, (mean, std)) - -# ------------------------------------------------------------------------------ -# Pareto distribution ---------------------------------------------------------- - - -class ParetoDistribution(SingleContinuousDistribution): - _argnames = ('xm', 'alpha') - - @property - def set(self): - return Interval(self.xm, oo, False, True) - - @staticmethod - def check(xm, alpha): - _value_check(xm > 0, 'Xm must be positive') - _value_check(alpha > 0, 'Alpha must be positive') - - def pdf(self, x): - xm, alpha = self.xm, self.alpha - return alpha * xm**alpha / x**(alpha + 1) - - def sample(self): - return random.paretovariate(self.alpha) - - -def Pareto(name, xm, alpha): - r""" - Create a continuous random variable with the Pareto distribution. - - The density of the Pareto distribution is given by - - .. math:: - f(x) := \frac{\alpha\,x_m^\alpha}{x^{\alpha+1}} - - with `x \in [x_m,\infty]`. - - Parameters - ========== - - xm : Real number, `x_m > 0`, a scale - alpha : Real number, `\alpha > 0`, a shape - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> xm = Symbol('xm', positive=True) - >>> beta = Symbol('beta', positive=True) - - >>> X = Pareto('x', xm, beta) - - >>> density(X)(z) - beta*xm**beta*z**(-beta - 1) - - References - ========== - - * https://en.wikipedia.org/wiki/Pareto_distribution - * https://mathworld.wolfram.com/ParetoDistribution.html - - """ - return rv(name, ParetoDistribution, (xm, alpha)) - -# ------------------------------------------------------------------------------ -# QuadraticU distribution ------------------------------------------------------ - - -class QuadraticUDistribution(SingleContinuousDistribution): - _argnames = ('a', 'b') - - @property - def set(self): - return Interval(self.a, self.b) - - def pdf(self, x): - a, b = self.a, self.b - alpha = 12/(b - a)**3 - beta = (a + b)/2 - return Piecewise((alpha*(x - beta)**2, And(a <= x, x <= b)), - (0, True)) - - -def QuadraticU(name, a, b): - r""" - Create a Continuous Random Variable with a U-quadratic distribution. - - The density of the U-quadratic distribution is given by - - .. math:: - f(x) := \alpha (x-\beta)^2 - - with `x \in [a,b]`. - - Parameters - ========== - - a : Real number - b : Real number, `a < b` - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> a, b = symbols('a b', real=True) - - >>> X = QuadraticU('x', a, b) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - / 2 - | / a b \ - |12*|- - - - + z| - | \ 2 2 / - <----------------- for And(a <= z, z <= b) - | 3 - | (-a + b) - | - \ 0 otherwise - - References - ========== - - * https://en.wikipedia.org/wiki/U-quadratic_distribution - - """ - return rv(name, QuadraticUDistribution, (a, b)) - -# ------------------------------------------------------------------------------ -# RaisedCosine distribution ---------------------------------------------------- - - -class RaisedCosineDistribution(SingleContinuousDistribution): - _argnames = ('mu', 's') - - @property - def set(self): - return Interval(self.mu - self.s, self.mu + self.s) - - @staticmethod - def check(mu, s): - _value_check(s > 0, 's must be positive') - - def pdf(self, x): - mu, s = self.mu, self.s - return Piecewise( - ((1 + cos(pi*(x - mu)/s))/(2*s), And(mu - s <= x, x <= mu + s)), - (0, True)) - - -def RaisedCosine(name, mu, s): - r""" - Create a Continuous Random Variable with a raised cosine distribution. - - The density of the raised cosine distribution is given by - - .. math:: - f(x) := \frac{1}{2s}\left(1+\cos\left(\frac{x-\mu}{s}\pi\right)\right) - - with `x \in [\mu-s,\mu+s]`. - - Parameters - ========== - - mu : Real number - s : Real number, `s > 0` - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> mu = Symbol('mu', real=True) - >>> s = Symbol('s', positive=True) - - >>> X = RaisedCosine('x', mu, s) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - / /pi*(-mu + z)\ - |cos|------------| + 1 - | \ s / - <--------------------- for And(z <= mu + s, mu - s <= z) - | 2*s - | - \ 0 otherwise - - References - ========== - - * https://en.wikipedia.org/wiki/Raised_cosine_distribution - - """ - return rv(name, RaisedCosineDistribution, (mu, s)) - -# ------------------------------------------------------------------------------ -# Rayleigh distribution -------------------------------------------------------- - - -class RayleighDistribution(SingleContinuousDistribution): - _argnames = 'sigma', - - set = Interval(0, oo, False, True) - - def pdf(self, x): - sigma = self.sigma - return x/sigma**2*exp(-x**2/(2*sigma**2)) - - -def Rayleigh(name, sigma): - r""" - Create a continuous random variable with a Rayleigh distribution. - - The density of the Rayleigh distribution is given by - - .. math :: - f(x) := \frac{x}{\sigma^2} e^{-x^2/2\sigma^2} - - with `x > 0`. - - Parameters - ========== - - sigma : Real number, `\sigma > 0` - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import E, density, variance - - >>> sigma = Symbol('sigma', positive=True) - - >>> X = Rayleigh('x', sigma) - - >>> density(X)(z) - E**(-z**2/(2*sigma**2))*z/sigma**2 - - >>> E(X) - sqrt(2)*sqrt(pi)*sigma/2 - - >>> variance(X) - -pi*sigma**2/2 + 2*sigma**2 - - References - ========== - - * https://en.wikipedia.org/wiki/Rayleigh_distribution - * https://mathworld.wolfram.com/RayleighDistribution.html - - """ - return rv(name, RayleighDistribution, (sigma, )) - -# ------------------------------------------------------------------------------ -# StudentT distribution -------------------------------------------------------- - - -class StudentTDistribution(SingleContinuousDistribution): - _argnames = 'nu', - - def pdf(self, x): - nu = self.nu - return 1/(sqrt(nu)*beta_fn(Rational(1, 2), nu/2))*(1 + x**2/nu)**(-(nu + 1)/2) - - -def StudentT(name, nu): - r""" - Create a continuous random variable with a student's t distribution. - - The density of the student's t distribution is given by - - .. math:: - f(x) := \frac{\Gamma \left(\frac{\nu+1}{2} \right)} - {\sqrt{\nu\pi}\Gamma \left(\frac{\nu}{2} \right)} - \left(1+\frac{x^2}{\nu} \right)^{-\frac{\nu+1}{2}} - - Parameters - ========== - - nu : Real number, `\nu > 0`, the degrees of freedom - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> nu = Symbol('nu', positive=True) - - >>> X = StudentT('x', nu) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - nu 1 - - -- - - - 2 2 - / 2\ - | z | - |1 + --| - \ nu/ - -------------------- - ____ / nu\ - \/ nu *beta|1/2, --| - \ 2 / - - References - ========== - - * https://en.wikipedia.org/wiki/Student_t-distribution - * https://mathworld.wolfram.com/Studentst-Distribution.html - - """ - return rv(name, StudentTDistribution, (nu, )) - -# ------------------------------------------------------------------------------ -# Triangular distribution ------------------------------------------------------ - - -class TriangularDistribution(SingleContinuousDistribution): - _argnames = ('a', 'b', 'c') - - def pdf(self, x): - a, b, c = self.a, self.b, self.c - return Piecewise( - (2*(x - a)/((b - a)*(c - a)), And(a <= x, x < c)), - (2/(b - a), Eq(x, c)), - (2*(b - x)/((b - a)*(b - c)), And(c < x, x <= b)), - (0, True)) - - -def Triangular(name, a, b, c): - r""" - Create a continuous random variable with a triangular distribution. - - The density of the triangular distribution is given by - - .. math:: - f(x) := \begin{cases} - 0 & \mathrm{for\ } x < a, \\ - \frac{2(x-a)}{(b-a)(c-a)} & \mathrm{for\ } a \le x < c, \\ - \frac{2}{b-a} & \mathrm{for\ } x = c, \\ - \frac{2(b-x)}{(b-a)(b-c)} & \mathrm{for\ } c < x \le b, \\ - 0 & \mathrm{for\ } b < x. - \end{cases} - - Parameters - ========== - - a : Real number, `a \in \left(-\infty, \infty\right)` - b : Real number, `a < b` - c : Real number, `a \leq c \leq b` - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> X = Triangular('x', a, b, c) - - >>> pprint(density(X)(z), use_unicode=False) - / -2*a + 2*z - |----------------- for And(a <= z, z < c) - |(-a + b)*(-a + c) - | - | 2 - | ------ for z = c - < -a + b - | - | 2*b - 2*z - |---------------- for And(z <= b, c < z) - |(-a + b)*(b - c) - | - \ 0 otherwise - - References - ========== - - * https://en.wikipedia.org/wiki/Triangular_distribution - * https://mathworld.wolfram.com/TriangularDistribution.html - - """ - return rv(name, TriangularDistribution, (a, b, c)) - -# ------------------------------------------------------------------------------ -# Uniform distribution --------------------------------------------------------- - - -class UniformDistribution(SingleContinuousDistribution): - _argnames = ('left', 'right') - - def pdf(self, x): - left, right = self.left, self.right - return Piecewise( - (1/(right - left), And(left <= x, x <= right)), - (0, True)) - - def compute_cdf(self, **kwargs): - from ..functions import Min - z = Dummy('z', real=True) - result = SingleContinuousDistribution.compute_cdf(self, **kwargs)(z) - reps = { - Min(z, self.right): z, - Min(z, self.left, self.right): self.left, - Min(z, self.left): self.left} - result = result.subs(reps) - return Lambda(z, result) - - def expectation(self, expr, var, **kwargs): - from ..functions import Max, Min - kwargs['evaluate'] = True - result = SingleContinuousDistribution.expectation(self, expr, var, **kwargs) - result = result.subs({Max(self.left, self.right): self.right, - Min(self.left, self.right): self.left}) - return result - - def sample(self): - return random.uniform(self.left, self.right) - - -def Uniform(name, left, right): - r""" - Create a continuous random variable with a uniform distribution. - - The density of the uniform distribution is given by - - .. math:: - f(x) := \begin{cases} - \frac{1}{b - a} & \text{for } x \in [a,b] \\ - 0 & \text{otherwise} - \end{cases} - - with `x \in [a,b]`. - - Parameters - ========== - - a : Real number, `-\infty < a` the left boundary - b : Real number, `a < b < \infty` the right boundary - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import E, cdf, density, variance - - >>> a = Symbol('a', negative=True) - >>> b = Symbol('b', positive=True) - - >>> X = Uniform('x', a, b) - - >>> density(X)(z) - Piecewise((1/(-a + b), (a <= z) & (z <= b)), (0, true)) - - >>> cdf(X)(z) - -a/(-a + b) + z/(-a + b) - - >>> simplify(E(X)) - a/2 + b/2 - - >>> simplify(variance(X)) - a**2/12 - a*b/6 + b**2/12 - - References - ========== - - * https://en.wikipedia.org/wiki/Uniform_distribution_%28continuous%29 - * https://mathworld.wolfram.com/UniformDistribution.html - - """ - return rv(name, UniformDistribution, (left, right)) - -# ------------------------------------------------------------------------------ -# UniformSum distribution ------------------------------------------------------ - - -class UniformSumDistribution(SingleContinuousDistribution): - _argnames = 'n', - - @property - def n(self): - return self.args[self._argnames.index('n')] - - @property - def set(self): - return Interval(0, self.n) - - def pdf(self, x): - n = self.n - k = Dummy('k') - return 1/factorial( - n - 1)*Sum((-1)**k*binomial(n, k)*(x - k)**(n - 1), (k, 0, floor(x))) - - -def UniformSum(name, n): - r""" - Create a continuous random variable with an Irwin-Hall distribution. - - The probability distribution function depends on a single parameter - `n` which is an integer. - - The density of the Irwin-Hall distribution is given by - - .. math :: - f(x) := \frac{1}{(n-1)!}\sum_{k=0}^{\lfloor x\rfloor}(-1)^k - \binom{n}{k}(x-k)^{n-1} - - Parameters - ========== - - n : A positive Integer, `n > 0` - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> X = UniformSum('x', n) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - floor(z) - ___ - \ ` - \ k n - 1 /n\ - ) (-1) *(z - k) *| | - / \k/ - /__, - k = 0 - ------------------------------- - (n - 1)! - - References - ========== - - * https://en.wikipedia.org/wiki/Uniform_sum_distribution - * https://mathworld.wolfram.com/UniformSumDistribution.html - - """ - return rv(name, UniformSumDistribution, (n, )) - -# ------------------------------------------------------------------------------ -# VonMises distribution -------------------------------------------------------- - - -class VonMisesDistribution(SingleContinuousDistribution): - _argnames = ('mu', 'k') - - set = Interval(0, 2*pi) - - @staticmethod - def check(mu, k): - _value_check(k > 0, 'k must be positive') - - def pdf(self, x): - mu, k = self.mu, self.k - return exp(k*cos(x-mu)) / (2*pi*besseli(0, k)) - - -def VonMises(name, mu, k): - r""" - Create a Continuous Random Variable with a von Mises distribution. - - The density of the von Mises distribution is given by - - .. math:: - f(x) := \frac{e^{\kappa\cos(x-\mu)}}{2\pi I_0(\kappa)} - - with `x \in [0,2\pi]`. - - Parameters - ========== - - mu : Real number, measure of location - k : Real number, measure of concentration - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> mu = Symbol('mu') - >>> k = Symbol('k', positive=True) - - >>> X = VonMises('x', mu, k) - - >>> D = density(X)(z) - >>> pprint(D, use_unicode=False) - k*cos(mu - z) - E - ------------------ - 2*pi*besseli(0, k) - - References - ========== - - * https://en.wikipedia.org/wiki/Von_Mises_distribution - * https://mathworld.wolfram.com/vonMisesDistribution.html - - """ - return rv(name, VonMisesDistribution, (mu, k)) - -# ------------------------------------------------------------------------------ -# Weibull distribution --------------------------------------------------------- - - -class WeibullDistribution(SingleContinuousDistribution): - _argnames = ('alpha', 'beta') - - set = Interval(0, oo, False, True) - - @staticmethod - def check(alpha, beta): - _value_check(alpha > 0, 'Alpha must be positive') - _value_check(beta > 0, 'Beta must be positive') - - def pdf(self, x): - alpha, beta = self.alpha, self.beta - return beta * (x/alpha)**(beta - 1) * exp(-(x/alpha)**beta) / alpha - - def sample(self): - return random.weibullvariate(self.alpha, self.beta) - - -def Weibull(name, alpha, beta): - r""" - Create a continuous random variable with a Weibull distribution. - - The density of the Weibull distribution is given by - - .. math:: - f(x) := \begin{cases} - \frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1} - e^{-(x/\lambda)^{k}} & x\geq0\\ - 0 & x<0 - \end{cases} - - Parameters - ========== - - lambda : Real number, `\lambda > 0` a scale - k : Real number, `k > 0` a shape - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import E, density, variance - - >>> l = Symbol('lambda', positive=True) - >>> k = Symbol('k', positive=True, real=True) - - >>> X = Weibull('x', l, k) - - >>> density(X)(z) - E**(-(z/lambda)**k)*k*(z/lambda)**(k - 1)/lambda - - >>> simplify(E(X)) - lambda*gamma(1 + 1/k) - - >>> simplify(variance(X)) - lambda**2*(-gamma(1 + 1/k)**2 + gamma(1 + 2/k)) - - References - ========== - - * https://en.wikipedia.org/wiki/Weibull_distribution - * https://mathworld.wolfram.com/WeibullDistribution.html - - """ - return rv(name, WeibullDistribution, (alpha, beta)) - -# ------------------------------------------------------------------------------ -# Wigner semicircle distribution ----------------------------------------------- - - -class WignerSemicircleDistribution(SingleContinuousDistribution): - _argnames = 'R', - - @property - def set(self): - return Interval(-self.R, self.R) - - def pdf(self, x): - R = self.R - return 2/(pi*R**2)*sqrt(R**2 - x**2) - - -def WignerSemicircle(name, R): - r""" - Create a continuous random variable with a Wigner semicircle distribution. - - The density of the Wigner semicircle distribution is given by - - .. math:: - f(x) := \frac2{\pi R^2}\,\sqrt{R^2-x^2} - - with `x \in [-R,R]`. - - Parameters - ========== - - R : Real number, `R > 0`, the radius - - Returns - ======= - - A `RandomSymbol`. - - Examples - ======== - - >>> from diofant.stats import E, density - - >>> R = Symbol('R', positive=True) - - >>> X = WignerSemicircle('x', R) - - >>> density(X)(z) - 2*sqrt(R**2 - z**2)/(pi*R**2) - - >>> E(X) - 0 - - References - ========== - - * https://en.wikipedia.org/wiki/Wigner_semicircle_distribution - * https://mathworld.wolfram.com/WignersSemicircleLaw.html - - """ - return rv(name, WignerSemicircleDistribution, (R,)) diff --git a/diofant/stats/drv.py b/diofant/stats/drv.py deleted file mode 100644 index 4268a8df32..0000000000 --- a/diofant/stats/drv.py +++ /dev/null @@ -1,80 +0,0 @@ -from ..concrete import Sum, summation -from ..core import Dummy, Expr, Lambda, cacheit, symbols, sympify -from ..functions import Piecewise -from ..sets import Integers -from .rv import NamedArgsMixin, SingleDomain, SinglePSpace - - -class SingleDiscreteDistribution(Expr, NamedArgsMixin): - """Discrete distribution of a single variable. - - Serves as superclass for PoissonDistribution etc.... - - Provides methods for pdf, cdf, and sampling - - See Also: - diofant.stats.crv_types.* - - """ - - set = Integers - - def __new__(cls, *args): - args = list(map(sympify, args)) - return Expr.__new__(cls, *args) - - @cacheit - def compute_cdf(self, **kwargs): - """Compute the CDF from the PDF. - - Returns a Lambda - - """ - x, z = symbols('x, z', integer=True, finite=True, cls=Dummy) - left_bound = self.set.inf - - # CDF is integral of PDF from left bound to z - pdf = self.pdf(x) - cdf = summation(pdf, (x, left_bound, z), **kwargs) - # CDF Ensure that CDF left of left_bound is zero - cdf = Piecewise((cdf, z >= left_bound), (0, True)) - return Lambda(z, cdf) - - def cdf(self, x, **kwargs): - """Cumulative density function.""" - return self.compute_cdf(**kwargs)(x) - - def expectation(self, expr, var, evaluate=True, **kwargs): - """Expectation of expression over distribution.""" - # TODO: support discrete sets with non integer stepsizes - if evaluate: - return summation(expr * self.pdf(var), - (var, self.set.inf, self.set.sup), **kwargs) - else: - return Sum(expr * self.pdf(var), - (var, self.set.inf, self.set.sup), **kwargs) - - def __call__(self, *args): - return self.pdf(*args) - - -class SingleDiscreteDomain(SingleDomain): - """Base class for a discrete domain.""" - - -class SingleDiscretePSpace(SinglePSpace): - """Discrete probability space over a single univariate variable.""" - - @property - def set(self): - return self.distribution.set - - @property - def domain(self): - return SingleDiscreteDomain(self.symbol, self.set) - - def integrate(self, expr, rvs=None, **kwargs): - rvs = rvs or (self.value,) - expr = expr.xreplace({rv: rv.symbol for rv in rvs}) - x = self.value.symbol - return self.distribution.expectation(expr, x, evaluate=False, **kwargs) diff --git a/diofant/stats/drv_types.py b/diofant/stats/drv_types.py deleted file mode 100644 index 196ee54e5c..0000000000 --- a/diofant/stats/drv_types.py +++ /dev/null @@ -1,134 +0,0 @@ -from ..core import sympify -from ..functions import exp, factorial -from ..sets import Naturals, Naturals0 -from .drv import SingleDiscreteDistribution, SingleDiscretePSpace -from .rv import _value_check - - -__all__ = 'Geometric', 'Poisson' - - -def rv(symbol, cls, *args): - args = list(map(sympify, args)) - dist = cls(*args) - dist.check(*args) - return SingleDiscretePSpace(symbol, dist).value - - -class PoissonDistribution(SingleDiscreteDistribution): - _argnames = 'lamda', - - set = Naturals0 - - @staticmethod - def check(lamda): - _value_check(lamda > 0, 'Lambda must be positive') - - def pdf(self, k): - return self.lamda**k / factorial(k) * exp(-self.lamda) - - -def Poisson(name, lamda): - r""" - Create a discrete random variable with a Poisson distribution. - - The density of the Poisson distribution is given by - - .. math:: - f(k) := \frac{\lambda^{k} e^{- \lambda}}{k!} - - Parameters - ========== - - lamda: Positive number, a rate - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import E, density, variance - - >>> rate = Symbol('lambda', positive=True) - - >>> X = Poisson('x', rate) - - >>> density(X)(z) - E**(-lambda)*lambda**z/factorial(z) - - >>> E(X) - lambda - - >>> simplify(variance(X)) - lambda - - References - ========== - - [1] https://en.wikipedia.org/wiki/Poisson_distribution - [2] https://mathworld.wolfram.com/PoissonDistribution.html - - """ - return rv(name, PoissonDistribution, lamda) - - -class GeometricDistribution(SingleDiscreteDistribution): - _argnames = 'p', - set = Naturals - - @staticmethod - def check(p): - _value_check(0 < p <= 1, 'p must be between 0 and 1') - - def pdf(self, k): - return (1 - self.p)**(k - 1) * self.p - - -def Geometric(name, p): - r""" - Create a discrete random variable with a Geometric distribution. - - The density of the Geometric distribution is given by - - .. math:: - f(k) := p (1 - p)^{k - 1} - - Parameters - ========== - - p: A probability between 0 and 1 - - Returns - ======= - - A RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import E, density, variance - - >>> p = Rational(1, 5) - - >>> X = Geometric('x', p) - - >>> density(X)(z) - (4/5)**(z - 1)/5 - - >>> E(X) - 5 - - >>> variance(X) - 20 - - References - ========== - - [1] https://en.wikipedia.org/wiki/Geometric_distribution - [2] https://mathworld.wolfram.com/GeometricDistribution.html - - """ - return rv(name, GeometricDistribution, p) diff --git a/diofant/stats/frv.py b/diofant/stats/frv.py deleted file mode 100644 index 65b928f854..0000000000 --- a/diofant/stats/frv.py +++ /dev/null @@ -1,349 +0,0 @@ -""" -Finite Discrete Random Variables Module - -See Also -======== -diofant.stats.frv_types -diofant.stats.rv -diofant.stats.crv - -""" - -import random -from itertools import product - -from ..core import Dict, Eq, Expr, Lambda, Mul, Symbol, Tuple, cacheit, sympify -from ..functions import Piecewise -from ..logic import And, Or -from ..sets import FiniteSet -from .rv import (ConditionalDomain, NamedArgsMixin, ProductDomain, - ProductPSpace, PSpace, RandomDomain, SinglePSpace, - random_symbols, rv_subs, sumsets) - - -class FiniteDensity(dict): - """Finite probabibility density.""" - - def __call__(self, item): - item = sympify(item) - if item in self: - return self[item] - else: - return 0 - - @property - def dict(self): - return dict(self) - - -class FiniteDomain(RandomDomain): - """ - A domain with discrete finite support - - Represented using a FiniteSet. - - """ - - is_Finite = True - - @property - def dict(self): - return FiniteSet(*[Dict(dict(el)) for el in self.elements]) - - def as_boolean(self): - return Or(*[And(*[Eq(sym, val) for sym, val in item]) for item in self]) - - def __iter__(self): - raise NotImplementedError - - -class SingleFiniteDomain(FiniteDomain): - """ - A FiniteDomain over a single symbol/set - - Example: The possibilities of a *single* die roll. - - """ - - def __new__(cls, symbol, set): - if not isinstance(set, FiniteSet): - set = FiniteSet(*set) - return Expr.__new__(cls, symbol, set) - - @property - def symbol(self): - return self.args[0] - - @property - def symbols(self): - return FiniteSet(self.symbol) - - @property - def set(self): - return self.args[1] - - @property - def elements(self): - return FiniteSet(*[frozenset(((self.symbol, elem), )) for elem in self.set]) - - def __iter__(self): - return (frozenset(((self.symbol, elem),)) for elem in self.set) - - def __contains__(self, other): - sym, val = tuple(other)[0] - return sym == self.symbol and val in self.set - - -class ProductFiniteDomain(ProductDomain, FiniteDomain): - """ - A Finite domain consisting of several other FiniteDomains - - Example: The possibilities of the rolls of three independent dice - - """ - - def __iter__(self): - proditer = product(*self.domains) - return (sumsets(items) for items in proditer) - - @property - def elements(self): - return FiniteSet(*self) - - -class ConditionalFiniteDomain(ConditionalDomain, ProductFiniteDomain): - """ - A FiniteDomain that has been restricted by a condition - - Example: The possibilities of a die roll under the condition that the - roll is even. - - """ - - def __new__(cls, domain, condition): - if condition is True: - return domain - cond = rv_subs(condition) - # Check that we aren't passed a condition like die1 == z - # where 'z' is a symbol that we don't know about - # We will never be able to test this equality through iteration - if not cond.free_symbols.issubset(domain.free_symbols): - raise ValueError(f'Condition "{condition}" contains foreign symbols' - f' \n{tuple(cond.free_symbols - domain.free_symbols)}.\n' - 'Will be unable to iterate using this condition') - - return Expr.__new__(cls, domain, cond) - - def _test(self, elem): - val = self.condition.xreplace(dict(elem)) - if val in [True, False]: - return val - elif val.is_Equality: - return val.lhs == val.rhs - raise ValueError(f'Undeciable if {val!s}') - - def __contains__(self, other): - return other in self.fulldomain and self._test(other) - - def __iter__(self): - return (elem for elem in self.fulldomain if self._test(elem)) - - @property - def set(self): - if self.fulldomain.__class__ is SingleFiniteDomain: - return FiniteSet(*[elem for elem in self.fulldomain.set - if frozenset(((self.fulldomain.symbol, elem),)) in self]) - else: - raise NotImplementedError( - 'Not implemented on multi-dimensional conditional domain') - - def as_boolean(self): - return FiniteDomain.as_boolean(self) - - -class SingleFiniteDistribution(Expr, NamedArgsMixin): - """Base class for finite distributions.""" - - def __new__(cls, *args): - args = list(map(sympify, args)) - return Expr.__new__(cls, *args) - - @property # type: ignore[misc] - @cacheit - def dict(self): - return {k: self.pdf(k) for k in self.set} - - @property - def pdf(self): - x = Symbol('x') - return Lambda(x, Piecewise(*( - [(v, Eq(k, x)) for k, v in self.dict.items()] + [(0, True)]))) - - @property - def set(self): - return list(self.dict) - - values = property(lambda self: self.dict.values) - items = property(lambda self: self.dict.items) - __iter__ = property(lambda self: self.dict.__iter__) - __getitem__ = property(lambda self: self.dict.__getitem__) - - __call__ = pdf - - def __contains__(self, other): - return other in self.set - - -##################### -# Probability Space # -##################### - - -class FinitePSpace(PSpace): - """ - A Finite Probability Space - - Represents the probabilities of a finite number of events. - - """ - - is_Finite = True - - @property - def domain(self): - return self.args[0] - - @property - def density(self): - return self.args[0] - - def __new__(cls, domain, density): - density = {sympify(key): sympify(val) - for key, val in density.items()} - public_density = Dict(density) - - obj = PSpace.__new__(cls, domain, public_density) - obj._density = density - return obj - - def prob_of(self, elem): - return self._density.get(elem, 0) - - def where(self, condition): - return ConditionalFiniteDomain(self.domain, condition) - - def compute_density(self, expr): - expr = expr.xreplace({rs: rs.symbol for rs in self.values}) - d = FiniteDensity() - for elem in self.domain: - val = expr.xreplace(dict(elem)) - prob = self.prob_of(elem) - d[val] = d.get(val, 0) + prob - return d - - @cacheit - def compute_cdf(self, expr): - d = self.compute_density(expr) - cum_prob = 0 - cdf = [] - for key in sorted(d): - prob = d[key] - cum_prob += prob - cdf.append((key, cum_prob)) - - return dict(cdf) - - @cacheit - def sorted_cdf(self, expr, python_float=False): - cdf = self.compute_cdf(expr) - items = list(cdf.items()) - sorted_items = sorted(items, key=lambda val_cumprob: val_cumprob[1]) - if python_float: - sorted_items = [(v, float(cum_prob)) - for v, cum_prob in sorted_items] - return sorted_items - - def integrate(self, expr, rvs=None): - rvs = rvs or self.values - expr = expr.xreplace({rs: rs.symbol for rs in rvs}) - return sum(expr.xreplace(dict(elem)) * self.prob_of(elem) - for elem in self.domain) - - def probability(self, condition): - cond_symbols = frozenset(rs.symbol for rs in random_symbols(condition)) - assert cond_symbols.issubset(self.symbols) - return sum(self.prob_of(elem) for elem in self.where(condition)) - - def conditional_space(self, condition): - domain = self.where(condition) - prob = self.probability(condition) - density = {key: val / prob - for key, val in self._density.items() if key in domain} - return FinitePSpace(domain, density) - - def sample(self): - """ - Internal sample method - - Returns dictionary mapping RandomSymbol to realization value. - - """ - expr = Tuple(*self.values) - cdf = self.sorted_cdf(expr, python_float=True) - - x = random.uniform(0, 1) - # Find first occurrence with cumulative probability less than x - # This should be replaced with binary search - for value, cum_prob in cdf: # pragma: no branch - if x < cum_prob: - # return dictionary mapping RandomSymbols to values - return dict(zip(expr, value)) - - -class SingleFinitePSpace(SinglePSpace, FinitePSpace): - """ - A single finite probability space - - Represents the probabilities of a set of random events that can be - attributed to a single variable/symbol. - - This class is implemented by many of the standard FiniteRV types such as - Die, Bernoulli, Coin, etc.... - - """ - - @property - def domain(self): - return SingleFiniteDomain(self.symbol, self.distribution.set) - - @property # type: ignore[misc] - @cacheit - def _density(self): - return {frozenset(((self.symbol, val),)): prob - for val, prob in self.distribution.dict.items()} - - -class ProductFinitePSpace(ProductPSpace, FinitePSpace): - """A collection of several independent finite probability spaces.""" - - @property - def domain(self): - return ProductFiniteDomain(*[space.domain for space in self.spaces]) - - @property # type: ignore[misc] - @cacheit - def _density(self): - proditer = product(*[iter(space._density.items()) - for space in self.spaces]) - d = {} - for items in proditer: - elems, probs = list(zip(*items)) - elem = sumsets(elems) - prob = Mul(*probs) - d[elem] = d.get(elem, 0) + prob - return Dict(d) - - @property # type: ignore[misc] - @cacheit - def density(self): - return self._density diff --git a/diofant/stats/frv_types.py b/diofant/stats/frv_types.py deleted file mode 100644 index 63c3d8f537..0000000000 --- a/diofant/stats/frv_types.py +++ /dev/null @@ -1,330 +0,0 @@ -""" -Finite Discrete Random Variables - Prebuilt variable types - -Contains -======== -FiniteRV -DiscreteUniform -Die -Bernoulli -Coin -Binomial -Hypergeometric -""" - -from ..concrete import Sum -from ..core import Dict, Dummy, Expr, Integer, Rational, cacheit, sympify -from ..core.compatibility import as_int -from ..core.logic import fuzzy_and, fuzzy_not -from ..functions import KroneckerDelta, binomial -from .frv import SingleFiniteDistribution, SingleFinitePSpace - - -__all__ = ('FiniteRV', 'DiscreteUniform', 'Die', 'Bernoulli', 'Coin', - 'Binomial', 'Hypergeometric') - - -def rv(name, cls, *args): - density = cls(*args) - return SingleFinitePSpace(name, density).value - - -class FiniteDistributionHandmade(SingleFiniteDistribution): - @property - def dict(self): - return self.args[0] - - def __new__(cls, density): - density = Dict(density) - return Expr.__new__(cls, density) - - -def FiniteRV(name, density): - """ - Create a Finite Random Variable given a dict representing the density. - - Returns a RandomSymbol. - - >>> from diofant.stats import E, P - - >>> density = {0: .1, 1: .2, 2: .3, 3: .4} - >>> X = FiniteRV('X', density) - - >>> E(X) - 2.00000000000000 - >>> P(X >= 2) - 0.700000000000000 - - """ - return rv(name, FiniteDistributionHandmade, density) - - -class DiscreteUniformDistribution(SingleFiniteDistribution): - @property - def p(self): - return Rational(1, len(self.args)) - - @property # type: ignore[misc] - @cacheit - def dict(self): - return {k: self.p for k in self.set} - - @property - def set(self): - return self.args - - def pdf(self, x): # pylint: disable=invalid-overridden-method - if x in self.args: - return self.p - else: - return Integer(0) - - -def DiscreteUniform(name, items): - """ - Create a Finite Random Variable representing a uniform distribution over - the input set. - - Returns a RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> X = DiscreteUniform('X', (a, b, c)) # equally likely over a, b, c - >>> density(X).dict - {a: 1/3, b: 1/3, c: 1/3} - - >>> Y = DiscreteUniform('Y', list(range(5))) # distribution over a range - >>> density(Y).dict - {0: 1/5, 1: 1/5, 2: 1/5, 3: 1/5, 4: 1/5} - - """ - return rv(name, DiscreteUniformDistribution, *items) - - -class DieDistribution(SingleFiniteDistribution): - _argnames = 'sides', - - def __new__(cls, sides): - sides_sym = sympify(sides) - if fuzzy_not(fuzzy_and((sides_sym.is_integer, sides_sym.is_positive))): - raise ValueError("'sides' must be a positive integer.") - return super().__new__(cls, sides) - - @property # type: ignore[misc] - @cacheit - def dict(self): - as_int(self.sides) - return super().dict - - @property - def set(self): - return list(map(Integer, range(1, self.sides + 1))) - - def pdf(self, x): # pylint: disable=invalid-overridden-method - x = sympify(x) - if x.is_number: - if x.is_Integer and 1 <= x <= self.sides: - return Rational(1, self.sides) - return Integer(0) - if x.is_Symbol: - i = Dummy('i', integer=True, positive=True) - return Sum(KroneckerDelta(x, i)/self.sides, (i, 1, self.sides)) - raise ValueError("'x' expected as an argument of type 'number' or 'symbol', " - f'not {type(x)}') - - -def Die(name, sides=6): - """ - Create a Finite Random Variable representing a fair die. - - Returns a RandomSymbol. - - >>> from diofant.stats import density - - >>> D6 = Die('D6', 6) # Six sided Die - >>> density(D6).dict - {1: 1/6, 2: 1/6, 3: 1/6, 4: 1/6, 5: 1/6, 6: 1/6} - - >>> D4 = Die('D4', 4) # Four sided Die - >>> density(D4).dict - {1: 1/4, 2: 1/4, 3: 1/4, 4: 1/4} - - """ - return rv(name, DieDistribution, sides) - - -class BernoulliDistribution(SingleFiniteDistribution): - _argnames = ('p', 'succ', 'fail') - - @property # type: ignore[misc] - @cacheit - def dict(self): - return {self.succ: self.p, self.fail: 1 - self.p} - - -def Bernoulli(name, p, succ=1, fail=0): - """ - Create a Finite Random Variable representing a Bernoulli process. - - Returns a RandomSymbol - - >>> from diofant.stats import density - - >>> X = Bernoulli('X', Rational(3, 4)) # 1-0 Bernoulli variable, probability = 3/4 - >>> density(X).dict - {0: 1/4, 1: 3/4} - - >>> X = Bernoulli('X', Rational(1, 2), 'Heads', 'Tails') # A fair coin toss - >>> density(X).dict - {Heads: 1/2, Tails: 1/2} - - """ - return rv(name, BernoulliDistribution, p, succ, fail) - - -def Coin(name, p=Rational(1, 2)): - """ - Create a Finite Random Variable representing a Coin toss. - - Probability p is the chance of getting "Heads." Half by default - - Returns a RandomSymbol. - - >>> from diofant.stats import density - - >>> H, T = Symbol('H'), Symbol('T') - - >>> C = Coin('C') # A fair coin toss - >>> density(C).dict - {H: 1/2, T: 1/2} - - >>> C2 = Coin('C2', Rational(3, 5)) # An unfair coin - >>> density(C2).dict - {H: 3/5, T: 2/5} - - """ - return rv(name, BernoulliDistribution, p, 'H', 'T') - - -class BinomialDistribution(SingleFiniteDistribution): - _argnames = ('n', 'p', 'succ', 'fail') - - def __new__(cls, *args): - n = args[BinomialDistribution._argnames.index('n')] - p = args[BinomialDistribution._argnames.index('p')] - n_sym = sympify(n) - p_sym = sympify(p) - - if fuzzy_not(fuzzy_and((n_sym.is_integer, n_sym.is_nonnegative))): - raise ValueError(f"'n' must be positive integer. n = {n!s}.") - if fuzzy_not(fuzzy_and((p_sym.is_nonnegative, (p_sym - 1).is_nonpositive))): - raise ValueError(f"'p' must be: 0 <= p <= 1 . p = {p!s}") - return super().__new__(cls, *args) - - @property - def n(self): - return self.args[self._argnames.index('n')] - - @property # type: ignore[misc] - @cacheit - def dict(self): - n, p, succ, fail = self.n, self.p, self.succ, self.fail - n = as_int(n) - return {k*succ + (n - k)*fail: - binomial(n, k) * p**k * (1 - p)**(n - k) for k in range(n + 1)} - - -def Binomial(name, n, p, succ=1, fail=0): - """ - Create a Finite Random Variable representing a binomial distribution. - - Returns a RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> X = Binomial('X', 4, Rational(1, 2)) # Four "coin flips" - >>> density(X).dict - {0: 1/16, 1: 1/4, 2: 3/8, 3: 1/4, 4: 1/16} - - """ - return rv(name, BinomialDistribution, n, p, succ, fail) - - -class HypergeometricDistribution(SingleFiniteDistribution): - _argnames = ('N', 'm', 'n') - - @property # type: ignore[misc] - @cacheit - def dict(self): - N, m, n = self.N, self.m, self.n - N, m, n = list(map(sympify, (N, m, n))) - density = {sympify(k): - Rational(binomial(m, k) * binomial(N - m, n - k), - binomial(N, n)) - for k in range(max(0, n + m - N), min(m, n) + 1)} - return density - - @property - def n(self): - return self.args[self._argnames.index('n')] - - -def Hypergeometric(name, N, m, n): - """ - Create a Finite Random Variable representing a hypergeometric distribution. - - Returns a RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> X = Hypergeometric('X', 10, 5, 3) # 10 marbles, 5 white (success), 3 draws - >>> density(X).dict - {0: 1/12, 1: 5/12, 2: 5/12, 3: 1/12} - - """ - return rv(name, HypergeometricDistribution, N, m, n) - - -class RademacherDistribution(SingleFiniteDistribution): - @property # type: ignore[misc] - @cacheit - def dict(self): - return {-1: Rational(1, 2), 1: Rational(1, 2)} - - -def Rademacher(name): - """ - Create a Finite Random Variable representing a Rademacher distribution. - - Return a RandomSymbol. - - Examples - ======== - - >>> from diofant.stats import density - - >>> X = Rademacher('X') - >>> density(X).dict - {-1: 1/2, 1: 1/2} - - See Also - ======== - - diofant.stats.Bernoulli - - References - ========== - - * https://en.wikipedia.org/wiki/Rademacher_distribution - - """ - return rv(name, RademacherDistribution) diff --git a/diofant/stats/rv.py b/diofant/stats/rv.py deleted file mode 100644 index 6e8d8a393c..0000000000 --- a/diofant/stats/rv.py +++ /dev/null @@ -1,1032 +0,0 @@ -""" -Main Random Variables Module. - -Defines abstract random variable type. -Contains interfaces for probability space object (PSpace) as well as standard -operators, P, E, sample, density, where - -See Also -======== -diofant.stats.crv -diofant.stats.frv -diofant.stats.rv_interface - -""" - -from __future__ import annotations - -from ..abc import x -from ..core import (Add, Eq, Equality, Expr, Integer, Lambda, Symbol, Tuple, - oo, sympify) -from ..core.logic import fuzzy_or -from ..core.relational import Relational -from ..functions import DiracDelta -from ..logic.boolalg import Boolean, false, true -from ..sets import FiniteSet, ProductSet -from ..solvers import solve -from ..utilities import lambdify - - -class RandomDomain(Expr): - """ - Represents a set of variables and the values which they can take - - See Also - ======== - diofant.stats.crv.ContinuousDomain - diofant.stats.frv.FiniteDomain - - """ - - is_ProductDomain = False - is_Finite = False - is_Continuous = False - - def __new__(cls, symbols, *args): - symbols = FiniteSet(*symbols) - return Expr.__new__(cls, symbols, *args) - - @property - def set(self): - return self.args[1] - - def __contains__(self, other): - raise NotImplementedError - - def integrate(self, expr): - raise NotImplementedError - - -class SingleDomain(RandomDomain): - """ - A single variable and its domain - - See Also - ======== - - diofant.stats.crv.SingleContinuousDomain - diofant.stats.frv.SingleFiniteDomain - - """ - - def __new__(cls, symbol, set): - assert symbol.is_Symbol - return Expr.__new__(cls, symbol, set) - - @property - def symbol(self): - return self.args[0] - - @property - def symbols(self): - return FiniteSet(self.symbol) - - -class ConditionalDomain(RandomDomain): - """ - A RandomDomain with an attached condition - - See Also - ======== - - diofant.stats.crv.ConditionalContinuousDomain - diofant.stats.frv.ConditionalFiniteDomain - - """ - - def __new__(cls, fulldomain, condition): - condition = condition.xreplace({rs: rs.symbol - for rs in random_symbols(condition)}) - return Expr.__new__(cls, fulldomain, condition) - - @property - def symbols(self): - return self.fulldomain.symbols - - @property - def fulldomain(self): - return self.args[0] - - @property - def condition(self): - return self.args[1] - - @property - def set(self): - raise NotImplementedError('Set of Conditional Domain not Implemented') - - -class PSpace(Expr): - """ - A Probability Space - - Probability Spaces encode processes that equal different values - probabilistically. These underly Random Symbols which occur in Diofant - expressions and contain the mechanics to evaluate statistical statements. - - See Also - ======== - diofant.stats.crv.ContinuousPSpace - diofant.stats.frv.FinitePSpace - - """ - - is_Finite: bool | None = None - is_Continuous: bool | None = None - - @property - def values(self): - return frozenset(RandomSymbol(self, sym) for sym in self.domain.symbols) - - @property - def symbols(self): - return self.domain.symbols - - def where(self, condition): - raise NotImplementedError - - def compute_density(self, expr): - raise NotImplementedError - - def sample(self): - raise NotImplementedError - - def probability(self, condition): - raise NotImplementedError - - def integrate(self, expr): - raise NotImplementedError - - -class SinglePSpace(PSpace): - """ - Represents the probabilities of a set of random events that can be - attributed to a single variable/symbol. - - """ - - def __new__(cls, s, distribution): - if isinstance(s, str): - s = Symbol(s) - if not isinstance(s, Symbol): - raise TypeError('s should have been string or Symbol') - return Expr.__new__(cls, s, distribution) - - @property - def value(self): - return RandomSymbol(self, self.symbol) - - @property - def symbol(self): - return self.args[0] - - @property - def distribution(self): - return self.args[1] - - -class RandomSymbol(Expr): - """ - Random Symbols represent ProbabilitySpaces in Diofant Expressions - In principle they can take on any value that their symbol can take on - within the associated PSpace with probability determined by the PSpace - Density. - - Random Symbols contain pspace and symbol properties. - The pspace property points to the represented Probability Space - The symbol is a standard Diofant Symbol that is used in that probability space - for example in defining a density. - - You can form normal Diofant expressions using RandomSymbols and operate on - those expressions with the Functions - - E - Expectation of a random expression - P - Probability of a condition - density - Probability Density of an expression - given - A new random expression (with new random symbols) given a condition - - An object of the RandomSymbol type should almost never be created by the - user. They tend to be created instead by the PSpace class's value method. - Traditionally a user doesn't even do this but instead calls one of the - convenience functions Normal, Exponential, Coin, Die, FiniteRV, etc.... - - """ - - def __new__(cls, pspace, symbol): - if not isinstance(symbol, Symbol): - raise TypeError('symbol should be of type Symbol') - if not isinstance(pspace, PSpace): - raise TypeError('pspace variable should be of type PSpace') - return Expr.__new__(cls, pspace, symbol) - - is_finite = True - is_Symbol = True - is_Atom = True - - _diff_wrt = True - - pspace = property(lambda self: self.args[0]) - symbol = property(lambda self: self.args[1]) - name = property(lambda self: self.symbol.name) - - def _eval_is_positive(self): - return self.symbol.is_positive - - def _eval_is_integer(self): - return self.symbol.is_integer - - def _eval_is_extended_real(self): - return fuzzy_or([self.symbol.is_extended_real, - self.pspace.is_extended_real]) - - def _eval_is_commutative(self): - return self.symbol.is_commutative - - def _hashable_content(self): - return self.pspace, self.symbol - - @property - def free_symbols(self): - return {self} - - -class ProductPSpace(PSpace): - """ - A probability space resulting from the merger of two independent probability - spaces. - - Often created using the function, pspace - - """ - - def __new__(cls, *spaces): - rs_space_dict = {} - for space in spaces: - for value in space.values: - rs_space_dict[value] = space - - symbols = FiniteSet(*[val.symbol for val in rs_space_dict]) - - # Overlapping symbols - if len(symbols) < sum(len(space.symbols) for space in spaces): - raise ValueError('Overlapping Random Variables') - - new_cls = cls - if all(space.is_Finite for space in spaces): - from .frv import ProductFinitePSpace - new_cls = ProductFinitePSpace - if all(space.is_Continuous for space in spaces): - from .crv import ProductContinuousPSpace - new_cls = ProductContinuousPSpace - - obj = Expr.__new__(new_cls, *FiniteSet(*spaces)) - - return obj - - @property - def rs_space_dict(self): - d = {} - for space in self.spaces: - for value in space.values: - d[value] = space - return d - - @property - def symbols(self): - return FiniteSet(*[val.symbol for val in self.rs_space_dict]) - - @property - def spaces(self): - return FiniteSet(*self.args) - - @property - def values(self): - return sumsets(space.values for space in self.spaces) - - def integrate(self, expr, rvs=None, **kwargs): - rvs = rvs or self.values - rvs = frozenset(rvs) - for space in self.spaces: - expr = space.integrate(expr, rvs & space.values, **kwargs) - return expr - - @property - def domain(self): - return ProductDomain(*[space.domain for space in self.spaces]) - - @property - def density(self): - raise NotImplementedError('Density not available for ProductSpaces') - - def sample(self): - return {k: v for space in self.spaces - for k, v in space.sample().items()} - - -class ProductDomain(RandomDomain): - """ - A domain resulting from the merger of two independent domains - - See Also - ======== - - diofant.stats.crv.ProductContinuousDomain - diofant.stats.frv.ProductFiniteDomain - - """ - - is_ProductDomain = True - - def __new__(cls, *domains): - # Flatten any product of products - domains2 = [] - for domain in domains: - if not domain.is_ProductDomain: - domains2.append(domain) - else: - domains2.extend(domain.domains) - domains2 = FiniteSet(*domains2) - - new_cls = cls - if all(domain.is_Finite for domain in domains2): - from .frv import ProductFiniteDomain - new_cls = ProductFiniteDomain - if all(domain.is_Continuous for domain in domains2): - from .crv import ProductContinuousDomain - new_cls = ProductContinuousDomain - - return Expr.__new__(new_cls, *domains2) - - @property - def symbols(self): - return FiniteSet(*[sym for domain in self.domains - for sym in domain.symbols]) - - @property - def domains(self): - return self.args - - @property - def set(self): - return ProductSet(domain.set for domain in self.domains) - - def __contains__(self, other): - # Split event into each subdomain - for domain in self.domains: - # Collect the parts of this event which associate to this domain - elem = frozenset(item for item in other - if domain.symbols.contains(item[0]) == true) - # Test this sub-event - if elem not in domain: - return False - # All subevents passed - return True - - -def random_symbols(expr): - """Returns all RandomSymbols within a Diofant Expression.""" - try: - return list(expr.atoms(RandomSymbol)) - except AttributeError: - return [] - - -def pspace(expr): - """ - Returns the underlying Probability Space of a random expression. - - For internal use. - - Examples - ======== - - >>> from diofant.stats import Normal - >>> X = Normal('X', 0, 1) - >>> pspace(2*X + 1) == X.pspace - True - - """ - expr = sympify(expr) - rvs = random_symbols(expr) - if not rvs: - raise ValueError(f'Expression containing Random Variable expected, not {expr}') - # If only one space present - if all(rv.pspace == rvs[0].pspace for rv in rvs): - return rvs[0].pspace - # Otherwise make a product space - return ProductPSpace(*[rv.pspace for rv in rvs]) - - -def sumsets(sets): - """Union of sets.""" - return frozenset().union(*sets) - - -def rs_swap(a, b): - """ - Build a dictionary to swap RandomSymbols based on their underlying symbol. - - i.e. - if ``X = ('x', pspace1)`` - and ``Y = ('x', pspace2)`` - then ``X`` and ``Y`` match and the key, value pair - ``{X:Y}`` will appear in the result - - Inputs: collections a and b of random variables which share common symbols - Output: dict mapping RVs in a to RVs in b - - """ - d = {} - for rsa in a: - d[rsa] = [rsb for rsb in b if rsa.symbol == rsb.symbol][0] - return d - - -def given(expr, condition=None, **kwargs): - r"""Conditional Random Expression. - - From a random expression and a condition on that expression creates a new - probability space from the condition and returns the same expression on that - conditional probability space. - - Examples - ======== - - >>> from diofant.stats import Die - >>> X = Die('X', 6) - >>> Y = given(X, X > 3) - >>> density(Y).dict - {4: 1/3, 5: 1/3, 6: 1/3} - - Following convention, if the condition is a random symbol then that symbol - is considered fixed. - - >>> from diofant.stats import Normal - - >>> X = Normal('X', 0, 1) - >>> Y = Normal('Y', 0, 1) - >>> pprint(density(X + Y, Y)(z), use_unicode=False) - 2 - -(-Y + z) - ----------- - ___ 2 - \/ 2 *E - ------------------ - ____ - 2*\/ pi - - """ - if not random_symbols(condition) or pspace_independent(expr, condition): - return expr - - if isinstance(condition, RandomSymbol): - condition = Eq(condition, condition.symbol) - - condsymbols = random_symbols(condition) - if (isinstance(condition, Equality) and len(condsymbols) == 1 and - not isinstance(pspace(expr).domain, ConditionalDomain)): - rv = tuple(condsymbols)[0] - results = solve(condition, rv) - return sum(expr.subs(res) for res in results) - - # Get full probability space of both the expression and the condition - fullspace = pspace(Tuple(expr, condition)) - # Build new space given the condition - space = fullspace.conditional_space(condition, **kwargs) - # Dictionary to swap out RandomSymbols in expr with new RandomSymbols - # That point to the new conditional space - swapdict = rs_swap(fullspace.values, space.values) - # Swap random variables in the expression - expr = expr.xreplace(swapdict) - return expr - - -def expectation(expr, condition=None, numsamples=None, evaluate=True, **kwargs): - """ - Returns the expected value of a random expression - - Parameters - ========== - - expr : Expr containing RandomSymbols - The expression of which you want to compute the expectation value - given : Expr containing RandomSymbols - A conditional expression. E(X, X>0) is expectation of X given X > 0 - numsamples : int - Enables sampling and approximates the expectation with this many samples - evalf : Bool (defaults to True) - If sampling return a number rather than a complex expression - evaluate : Bool (defaults to True) - In case of continuous systems return unevaluated integral - - Examples - ======== - - >>> from diofant.stats import Die, E - >>> X = Die('X', 6) - >>> E(X) - 7/2 - >>> E(2*X + 1) - 8 - - >>> E(X, X > 3) # Expectation of X given that it is above 3 - 5 - - """ - if not random_symbols(expr): # expr isn't random? - return expr - if numsamples: # Computing by monte carlo sampling? - return sampling_E(expr, condition, numsamples=numsamples, **kwargs) - - # Create new expr and recompute E - if condition is not None: # If there is a condition - return expectation(given(expr, condition), evaluate=evaluate) - - # A few known statements for efficiency - - if expr.is_Add: # We know that E is Linear - return Add(*[expectation(arg, evaluate=evaluate) - for arg in expr.args]) - - # Otherwise case is simple, pass work off to the ProbabilitySpace - result = pspace(expr).integrate(expr) - if evaluate and hasattr(result, 'doit'): - return result.doit(**kwargs) - else: - return result - - -def probability(condition, given_condition=None, numsamples=None, - evaluate=True, **kwargs): - """ - Probability that a condition is true, optionally given a second condition - - Parameters - ========== - - condition : Combination of Relationals containing RandomSymbols - The condition of which you want to compute the probability - given_condition : Combination of Relationals containing RandomSymbols - A conditional expression. P(X > 1, X > 0) is expectation of X > 1 - given X > 0 - numsamples : int - Enables sampling and approximates the probability with this many samples - evaluate : Bool (defaults to True) - In case of continuous systems return unevaluated integral - - Examples - ======== - - >>> from diofant.stats import Die, P - >>> X, Y = Die('X', 6), Die('Y', 6) - >>> P(X > 3) - 1/2 - >>> P(Eq(X, 5), X > 2) # Probability that X == 5 given that X > 2 - 1/4 - >>> P(X > Y) - 5/12 - - """ - condition = sympify(condition) - given_condition = sympify(given_condition) - - if given_condition is not None and \ - not isinstance(given_condition, (Relational, Boolean)): - raise ValueError(f'{given_condition} is not a relational or combination of relationals') - if given_condition == false: - return Integer(0) - if not isinstance(condition, (Relational, Boolean)): - raise ValueError(f'{condition} is not a relational or combination of relationals') - if condition == true: - return Integer(1) - if condition == false: - return Integer(0) - - if numsamples: - return sampling_P(condition, given_condition, numsamples=numsamples, - **kwargs) - if given_condition is not None: # If there is a condition - # Recompute on new conditional expr - return probability(given(condition, given_condition, **kwargs), **kwargs) - - # Otherwise pass work off to the ProbabilitySpace - result = pspace(condition).probability(condition, **kwargs) - if evaluate and hasattr(result, 'doit'): - return result.doit() - else: - return result - - -class Density(Expr): - """Probability density.""" - - expr = property(lambda self: self.args[0]) - - @property - def condition(self): - if len(self.args) > 1: - return self.args[1] - - def doit(self, **kwargs): - evaluate = kwargs.pop('evaluate', True) - - expr, condition = self.expr, self.condition - if condition is not None: - # Recompute on new conditional expr - expr = given(expr, condition, **kwargs) - if not random_symbols(expr): - return Lambda(x, DiracDelta(x - expr)) - if (isinstance(expr, RandomSymbol) and - hasattr(expr.pspace, 'distribution') and - isinstance(pspace(expr), SinglePSpace)): - return expr.pspace.distribution - result = pspace(expr).compute_density(expr, **kwargs) - - if evaluate and hasattr(result, 'doit'): - return result.doit() - else: - return result - - -def density(expr, condition=None, evaluate=True, numsamples=None, **kwargs): - """ - Probability density of a random expression, optionally given a second - condition. - - This density will take on different forms for different types of - probability spaces. Discrete variables produce Dicts. Continuous - variables produce Lambdas. - - Parameters - ========== - - expr : Expr containing RandomSymbols - The expression of which you want to compute the density value - condition : Relational containing RandomSymbols - A conditional expression. density(X > 1, X > 0) is density of X > 1 - given X > 0 - numsamples : int - Enables sampling and approximates the density with this many samples - - Examples - ======== - - >>> from diofant.stats import Die, Normal - - >>> D = Die('D', 6) - >>> X = Normal(x, 0, 1) - - >>> density(D).dict - {1: 1/6, 2: 1/6, 3: 1/6, 4: 1/6, 5: 1/6, 6: 1/6} - >>> density(2*D).dict - {2: 1/6, 4: 1/6, 6: 1/6, 8: 1/6, 10: 1/6, 12: 1/6} - >>> density(X)(x) - sqrt(2)*E**(-x**2/2)/(2*sqrt(pi)) - - """ - if numsamples: - return sampling_density(expr, condition, numsamples=numsamples, - **kwargs) - - kwargs['evaluate'] = evaluate - - return Density(expr, condition).doit(**kwargs) - - -def cdf(expr, condition=None, evaluate=True, **kwargs): - """ - Cumulative Distribution Function of a random expression. - - optionally given a second condition - - This density will take on different forms for different types of - probability spaces. - Discrete variables produce Dicts. - Continuous variables produce Lambdas. - - Examples - ======== - - >>> from diofant.stats import Die, Normal - - >>> D = Die('D', 6) - >>> X = Normal('X', 0, 1) - - >>> density(D).dict - {1: 1/6, 2: 1/6, 3: 1/6, 4: 1/6, 5: 1/6, 6: 1/6} - >>> cdf(D) - {1: 1/6, 2: 1/3, 3: 1/2, 4: 2/3, 5: 5/6, 6: 1} - >>> cdf(3*D, D > 2) - {9: 1/4, 12: 1/2, 15: 3/4, 18: 1} - - >>> cdf(X) - Lambda(_z, erf(sqrt(2)*_z/2)/2 + 1/2) - - """ - if condition is not None: # If there is a condition - # Recompute on new conditional expr - return cdf(given(expr, condition, **kwargs), **kwargs) - - # Otherwise pass work off to the ProbabilitySpace - result = pspace(expr).compute_cdf(expr, **kwargs) - - if evaluate and hasattr(result, 'doit'): - return result.doit() - else: - return result - - -def where(condition, given_condition=None, **kwargs): - """ - Returns the domain where a condition is True. - - Examples - ======== - - >>> from diofant.stats import Die, Normal - - >>> D1, D2 = Die('a', 6), Die('b', 6) - >>> a, b = D1.symbol, D2.symbol - >>> X = Normal('x', 0, 1) - - >>> where(X**2 < 1) - Domain: (-1 < x) & (x < 1) - - >>> where(X**2 < 1).set - (-1, 1) - - >>> where(And(D1 <= D2, D2 < 3)) - Domain: (Eq(a, 1) & Eq(b, 1)) | (Eq(a, 1) & Eq(b, 2)) | (Eq(a, 2) & Eq(b, 2)) - - """ - if given_condition is not None: # If there is a condition - # Recompute on new conditional expr - return where(given(condition, given_condition, **kwargs), **kwargs) - - # Otherwise pass work off to the ProbabilitySpace - return pspace(condition).where(condition, **kwargs) - - -def sample(expr, condition=None, **kwargs): - """ - A realization of the random expression - - Examples - ======== - - >>> from diofant.stats import Die - >>> X, Y, Z = Die('X', 6), Die('Y', 6), Die('Z', 6) - - >>> die_roll = sample(X + Y + Z) # A random realization of three dice - - """ - return next(sample_iter(expr, condition, numsamples=1)) - - -def sample_iter(expr, condition=None, numsamples=oo, **kwargs): - """ - Returns an iterator of realizations from the expression given a condition - - expr: Random expression to be realized - condition: A conditional expression (optional) - numsamples: Length of the iterator (defaults to infinity) - - Examples - ======== - - >>> from diofant.stats import Normal - >>> X = Normal('X', 0, 1) - >>> expr = X*X + 3 - >>> iterator = sample_iter(expr, numsamples=3) - >>> list(iterator) # doctest: +SKIP - [12, 4, 7] - - See Also - ======== - - diofant.stats.sample - diofant.stats.rv.sampling_P - diofant.stats.rv.sampling_E - - """ - if condition is not None: - ps = pspace(Tuple(expr, condition)) - else: - ps = pspace(expr) - - rvs = list(ps.values) - fn = lambdify(rvs, expr, **kwargs) - if condition is not None: - given_fn = lambdify(rvs, condition, **kwargs) - - # Check that lambdify can handle the expression - # Some operations like Sum can prove difficult - try: - d = ps.sample() # a dictionary that maps RVs to values - args = [d[rv] for rv in rvs] - fn(*args) - if condition is not None: - given_fn(*args) - except (TypeError, ValueError) as exc: - raise TypeError('Expr/condition too complex for lambdify') from exc - - def return_generator(): - count = 0 - while count < numsamples: - d = ps.sample() # a dictionary that maps RVs to values - args = [d[rv] for rv in rvs] - - if condition is not None: # Check that these values satisfy the condition - gd = given_fn(*args) - if gd not in (True, False): - raise ValueError( - 'Conditions must not contain free symbols') - if not gd: # If the values don't satisfy then try again - continue - - yield fn(*args) - count += 1 - return return_generator() - - -def sampling_P(condition, given_condition=None, numsamples=1, - evalf=True, **kwargs): - """ - Sampling version of P - - See Also - ======== - diofant.stats.P - diofant.stats.rv.sampling_E - diofant.stats.rv.sampling_density - - """ - count_true = 0 - count_false = 0 - - samples = sample_iter(condition, given_condition, - numsamples=numsamples, **kwargs) - - for x in samples: - if x: - count_true += 1 - else: - count_false += 1 - - result = Integer(count_true) / numsamples - return result.evalf() - - -def sampling_E(expr, given_condition=None, numsamples=1, - evalf=True, **kwargs): - """ - Sampling version of E - - See Also - ======== - diofant.stats.P - diofant.stats.rv.sampling_P - diofant.stats.rv.sampling_density - - """ - samples = sample_iter(expr, given_condition, - numsamples=numsamples, **kwargs) - - result = Add(*list(samples)) / numsamples - return result.evalf(strict=False) - - -def sampling_density(expr, given_condition=None, numsamples=1, **kwargs): - """ - Sampling version of density - - See Also - ======== - diofant.stats.density - diofant.stats.rv.sampling_P - diofant.stats.rv.sampling_E - - """ - results = {} - for result in sample_iter(expr, given_condition, - numsamples=numsamples, **kwargs): - results[result] = results.get(result, 0) + 1 - return results - - -def dependent(a, b): - """ - Dependence of two random expressions - - Two expressions are independent if knowledge of one does not change - computations on the other. - - Examples - ======== - - >>> from diofant.stats import Normal - - >>> X, Y = Normal('X', 0, 1), Normal('Y', 0, 1) - >>> dependent(X, Y) - False - >>> dependent(2*X + Y, -Y) - True - >>> X, Y = given(Tuple(X, Y), Eq(X + Y, 3)) - >>> dependent(X, Y) - True - - See Also - ======== - diofant.stats.rv.independent - - """ - if pspace_independent(a, b): - return False - - z = Symbol('z', extended_real=True) - # Dependent if density is unchanged when one is given information about - # the other - return (density(a, Eq(b, z)) != density(a) or - density(b, Eq(a, z)) != density(b)) - - -def independent(a, b): - """ - Independence of two random expressions - - Two expressions are independent if knowledge of one does not change - computations on the other. - - Examples - ======== - - >>> from diofant.stats import Normal - - >>> X, Y = Normal('X', 0, 1), Normal('Y', 0, 1) - >>> independent(X, Y) - True - >>> independent(2*X + Y, -Y) - False - >>> X, Y = given(Tuple(X, Y), Eq(X + Y, 3)) - >>> independent(X, Y) - False - - See Also - ======== - diofant.stats.rv.dependent - - """ - return not dependent(a, b) - - -def pspace_independent(a, b): - """ - Tests for independence between a and b by checking if their PSpaces have - overlapping symbols. This is a sufficient but not necessary condition for - independence and is intended to be used internally. - - Notes - ===== - - pspace_independent(a, b) implies independent(a, b) - independent(a, b) does not imply pspace_independent(a, b) - - """ - a_symbols = set(pspace(b).symbols) - b_symbols = set(pspace(a).symbols) - - if len(a_symbols.intersection(b_symbols)) == 0: - return True - - -def rv_subs(expr): - """Given a random expression replace all random variables with their symbols.""" - symbols = random_symbols(expr) - if not symbols: - return expr - swapdict = {rv: rv.symbol for rv in symbols} - return expr.xreplace(swapdict) - - -class NamedArgsMixin: - """Helper class for named arguments.""" - - _argnames: tuple[str, ...] = () - - def __getattr__(self, attr): - try: - return self.args[self._argnames.index(attr)] - except ValueError as exc: - raise AttributeError(f"'{type(self).__name__}' object " - f"has not attribute '{attr}'") from exc - - -def _value_check(condition, message): - """ - Check a condition on input value. - - Raises ValueError with message if condition is not True - - """ - if condition != true: - raise ValueError(message) diff --git a/diofant/stats/rv_interface.py b/diofant/stats/rv_interface.py deleted file mode 100644 index 147f17fcc9..0000000000 --- a/diofant/stats/rv_interface.py +++ /dev/null @@ -1,216 +0,0 @@ -from ..functions import sqrt -from .rv import (cdf, density, dependent, expectation, given, independent, - probability, pspace, random_symbols, sample, sample_iter, - sampling_density, where) - - -__all__ = ('P', 'E', 'density', 'where', 'given', 'sample', 'cdf', 'pspace', - 'sample_iter', 'variance', 'std', 'skewness', 'covariance', - 'dependent', 'independent', 'random_symbols', 'correlation', - 'moment', 'cmoment', 'sampling_density') - - -def moment(X, n, c=0, condition=None, **kwargs): - """ - Return the nth moment of a random expression about c i.e. E((X-c)**n) - Default value of c is 0. - - Examples - ======== - - >>> from diofant.stats import Die, E - >>> X = Die('X', 6) - >>> moment(X, 1, 6) - -5/2 - >>> moment(X, 2) - 91/6 - >>> moment(X, 1) == E(X) - True - - """ - return expectation((X - c)**n, condition, **kwargs) - - -def variance(X, condition=None, **kwargs): - """ - Variance of a random expression - - Expectation of (X-E(X))**2 - - Examples - ======== - - >>> from diofant.stats import Bernoulli, Die - - >>> X = Die('X', 6) - >>> p = Symbol('p') - >>> B = Bernoulli('B', p, 1, 0) - - >>> variance(2*X) - 35/3 - - >>> simplify(variance(B)) - p*(-p + 1) - - """ - return cmoment(X, 2, condition, **kwargs) - - -def standard_deviation(X, condition=None, **kwargs): - """ - Standard Deviation of a random expression - - Square root of the Expectation of (X-E(X))**2 - - Examples - ======== - - >>> from diofant.stats import Bernoulli - - >>> p = Symbol('p') - >>> B = Bernoulli('B', p, 1, 0) - - >>> simplify(std(B)) - sqrt(p*(-p + 1)) - - """ - return sqrt(variance(X, condition, **kwargs)) - - -std = standard_deviation - - -def covariance(X, Y, condition=None, **kwargs): - """ - Covariance of two random expressions - - The expectation that the two variables will rise and fall together - - Covariance(X,Y) = E( (X-E(X)) * (Y-E(Y)) ) - - Examples - ======== - - >>> from diofant.stats import Exponential - - >>> rate = Symbol('lambda', positive=True, real=True) - >>> X = Exponential('X', rate) - >>> Y = Exponential('Y', rate) - - >>> covariance(X, X) - lambda**(-2) - >>> covariance(X, Y) - 0 - >>> covariance(X, Y + rate*X) - 1/lambda - - """ - return expectation( - (X - expectation(X, condition, **kwargs)) * - (Y - expectation(Y, condition, **kwargs)), - condition, **kwargs) - - -def correlation(X, Y, condition=None, **kwargs): - """ - Correlation of two random expressions, also known as correlation - coefficient or Pearson's correlation - - The normalized expectation that the two variables will rise - and fall together - - Correlation(X,Y) = E( (X-E(X)) * (Y-E(Y)) / (sigma(X) * sigma(Y)) ) - - Examples - ======== - - >>> from diofant.stats import Exponential - - >>> rate = Symbol('lambda', positive=True, real=True) - >>> X = Exponential('X', rate) - >>> Y = Exponential('Y', rate) - - >>> correlation(X, X) - 1 - >>> correlation(X, Y) - 0 - >>> correlation(X, Y + rate*X) - 1/sqrt(1 + lambda**(-2)) - - """ - return covariance(X, Y, condition, **kwargs)/(std(X, condition, **kwargs) - * std(Y, condition, **kwargs)) - - -def cmoment(X, n, condition=None, **kwargs): - """ - Return the nth central moment of a random expression about its mean - i.e. E((X - E(X))**n) - - Examples - ======== - - >>> from diofant.stats import Die - >>> X = Die('X', 6) - >>> cmoment(X, 3) - 0 - >>> cmoment(X, 2) - 35/12 - >>> cmoment(X, 2) == variance(X) - True - - """ - mu = expectation(X, condition, **kwargs) - return moment(X, n, mu, condition, **kwargs) - - -def smoment(X, n, condition=None, **kwargs): - """ - Return the nth Standardized moment of a random expression i.e. - E( ((X - mu)/sigma(X))**n ) - - Examples - ======== - - >>> from diofant.stats import Exponential - >>> rate = Symbol('lambda', positive=True, real=True) - >>> Y = Exponential('Y', rate) - >>> smoment(Y, 4) - 9 - >>> smoment(Y, 4) == smoment(3*Y, 4) - True - >>> smoment(Y, 3) == skewness(Y) - True - - """ - sigma = std(X, condition, **kwargs) - return (1/sigma)**n*cmoment(X, n, condition, **kwargs) - - -def skewness(X, condition=None, **kwargs): - """ - Measure of the asymmetry of the probability distribution - - Positive skew indicates that most of the values lie to the right of - the mean - - skewness(X) = E( ((X - E(X))/sigma)**3 ) - - Examples - ======== - - >>> from diofant.stats import Exponential, Normal - >>> X = Normal('X', 0, 1) - >>> skewness(X) - 0 - >>> rate = Symbol('lambda', positive=True, real=True) - >>> Y = Exponential('Y', rate) - >>> skewness(Y) - 2 - - """ - return smoment(X, 3, condition, **kwargs) - - -P = probability -E = expectation diff --git a/diofant/tests/core/test_args.py b/diofant/tests/core/test_args.py index f5ab7127f6..58227d24a9 100644 --- a/diofant/tests/core/test_args.py +++ b/diofant/tests/core/test_args.py @@ -23,7 +23,7 @@ RootSum, S, Set, StrictGreaterThan, StrictLessThan, Subs, Subset, Sum, Symbol, SymmetricDifference, Trace, Transpose, Tuple, Unequality, Union, Wild, WildFunction, - Xor, ZeroMatrix, divisor_sigma, false, mobius, oo, sin, + Xor, ZeroMatrix, divisor_sigma, false, mobius, sin, symbols, totient, true) from diofant.abc import a, b, w, x, y, z from diofant.concrete.expr_with_intlimits import ExprWithIntLimits @@ -90,52 +90,6 @@ Naturals0, Rationals, Reals) from diofant.sets.sets import EmptySet, UniversalSet from diofant.simplify.hyperexpand import G_Function, Hyper_Function -from diofant.stats.crv import (ConditionalContinuousDomain, - ContinuousDistributionHandmade, - ContinuousDomain, ContinuousPSpace, - ProductContinuousDomain, - ProductContinuousPSpace, SingleContinuousDomain, - SingleContinuousPSpace) -from diofant.stats.crv_types import (ArcsinDistribution, BeniniDistribution, - BetaDistribution, BetaPrimeDistribution, - CauchyDistribution, ChiDistribution, - ChiNoncentralDistribution, - ChiSquaredDistribution, DagumDistribution, - ExponentialDistribution, - FDistributionDistribution, - FisherZDistribution, FrechetDistribution, - GammaDistribution, - GammaInverseDistribution, - KumaraswamyDistribution, - LaplaceDistribution, LogisticDistribution, - LogNormalDistribution, - MaxwellDistribution, NakagamiDistribution, - Normal, NormalDistribution, - ParetoDistribution, - QuadraticUDistribution, - RaisedCosineDistribution, - RayleighDistribution, - StudentTDistribution, - TriangularDistribution, - UniformDistribution, - UniformSumDistribution, - VonMisesDistribution, WeibullDistribution, - WignerSemicircleDistribution) -from diofant.stats.drv import SingleDiscreteDomain, SingleDiscretePSpace -from diofant.stats.drv_types import GeometricDistribution, PoissonDistribution -from diofant.stats.frv import (ConditionalFiniteDomain, FiniteDomain, - FinitePSpace, ProductFiniteDomain, - ProductFinitePSpace, SingleFiniteDomain, - SingleFinitePSpace) -from diofant.stats.frv_types import (BernoulliDistribution, - BinomialDistribution, DieDistribution, - DiscreteUniformDistribution, - FiniteDistributionHandmade, - HypergeometricDistribution, - RademacherDistribution) -from diofant.stats.rv import (ConditionalDomain, Density, ProductDomain, - ProductPSpace, PSpace, RandomDomain, - RandomSymbol, SingleDomain) from diofant.tensor import ImmutableDenseNDimArray, ImmutableSparseNDimArray from diofant.tensor.tensor import (TensAdd, TensorHead, TensorIndex, TensorIndexType, TensorSymmetry, TensorType, @@ -568,321 +522,6 @@ def test_diofant__sets__contains__Contains(): assert _test_args(Contains(x, Range(0, 10, 2))) -# STATS - - -nd = NormalDistribution(0, 1) -die = DieDistribution(6) - - -def test_diofant__stats__crv__ContinuousDomain(): - assert _test_args(ContinuousDomain({x}, Interval(-oo, oo))) - - -def test_diofant__stats__crv__SingleContinuousDomain(): - assert _test_args(SingleContinuousDomain(x, Interval(-oo, oo))) - - -def test_diofant__stats__crv__ProductContinuousDomain(): - D = SingleContinuousDomain(x, Interval(-oo, oo)) - E = SingleContinuousDomain(y, Interval(0, oo)) - assert _test_args(ProductContinuousDomain(D, E)) - - -def test_diofant__stats__crv__ConditionalContinuousDomain(): - D = SingleContinuousDomain(x, Interval(-oo, oo)) - assert _test_args(ConditionalContinuousDomain(D, x > 0)) - - -def test_diofant__stats__crv__ContinuousPSpace(): - D = SingleContinuousDomain(x, Interval(-oo, oo)) - assert _test_args(ContinuousPSpace(D, nd)) - - -def test_diofant__stats__crv__SingleContinuousPSpace(): - assert _test_args(SingleContinuousPSpace(x, nd)) - - -def test_diofant__stats__crv__ProductContinuousPSpace(): - A = SingleContinuousPSpace(x, nd) - B = SingleContinuousPSpace(y, nd) - assert _test_args(ProductContinuousPSpace(A, B)) - - -def test_diofant__stats__crv__SingleContinuousDistribution(): - pass - - -def test_diofant__stats__drv__SingleDiscreteDomain(): - assert _test_args(SingleDiscreteDomain(x, S.Naturals)) - - -def test_diofant__stats__drv__SingleDiscretePSpace(): - assert _test_args(SingleDiscretePSpace(x, PoissonDistribution(1))) - - -def test_diofant__stats__drv__SingleDiscreteDistribution(): - pass - - -def test_diofant__stats__rv__RandomDomain(): - assert _test_args(RandomDomain(FiniteSet(x), FiniteSet(1, 2, 3))) - - -def test_diofant__stats__rv__SingleDomain(): - assert _test_args(SingleDomain(x, FiniteSet(1, 2, 3))) - - -def test_diofant__stats__rv__ConditionalDomain(): - D = RandomDomain(FiniteSet(x), FiniteSet(1, 2)) - assert _test_args(ConditionalDomain(D, x > 1)) - - -def test_diofant__stats__rv__PSpace(): - D = RandomDomain(FiniteSet(x), FiniteSet(1, 2, 3, 4, 5, 6)) - assert _test_args(PSpace(D, die)) - - -def test_diofant__stats__rv__SinglePSpace(): - pass - - -def test_diofant__stats__rv__RandomSymbol(): - A = SingleContinuousPSpace(x, nd) - assert _test_args(RandomSymbol(A, x)) - - -def test_diofant__stats__rv__ProductPSpace(): - A = SingleContinuousPSpace(x, nd) - B = SingleContinuousPSpace(y, nd) - assert _test_args(ProductPSpace(A, B)) - - -def test_diofant__stats__rv__ProductDomain(): - D = SingleDomain(x, Interval(-oo, oo)) - E = SingleDomain(y, Interval(0, oo)) - assert _test_args(ProductDomain(D, E)) - - -def test_diofant__stats__frv_types__DiscreteUniformDistribution(): - assert _test_args(DiscreteUniformDistribution(Tuple(*range(6)))) - - -def test_diofant__stats__frv_types__DieDistribution(): - assert _test_args(DieDistribution(6)) - - -def test_diofant__stats__frv_types__BernoulliDistribution(): - assert _test_args(BernoulliDistribution(Rational(1, 2), 0, 1)) - - -def test_diofant__stats__frv_types__BinomialDistribution(): - assert _test_args(BinomialDistribution(5, Rational(1, 2), 1, 0)) - - -def test_diofant__stats__frv_types__HypergeometricDistribution(): - assert _test_args(HypergeometricDistribution(10, 5, 3)) - - -def test_diofant__stats__frv_types__RademacherDistribution(): - assert _test_args(RademacherDistribution()) - - -def test_diofant__stats__frv__FiniteDomain(): - assert _test_args(FiniteDomain({(x, 1), (x, 2)})) # x can be 1 or 2 - - -def test_diofant__stats__frv__SingleFiniteDomain(): - assert _test_args(SingleFiniteDomain(x, {1, 2})) # x can be 1 or 2 - - -def test_diofant__stats__frv__ProductFiniteDomain(): - xd = SingleFiniteDomain(x, {1, 2}) - yd = SingleFiniteDomain(y, {1, 2}) - assert _test_args(ProductFiniteDomain(xd, yd)) - - -def test_diofant__stats__frv__ConditionalFiniteDomain(): - xd = SingleFiniteDomain(x, {1, 2}) - assert _test_args(ConditionalFiniteDomain(xd, x > 1)) - - -def test_diofant__stats__frv__FinitePSpace(): - xd = SingleFiniteDomain(x, {1, 2}) - assert _test_args(FinitePSpace(xd, {(x, 1): Rational(1, 2), (x, 2): Rational(1, 2)})) - - -def test_diofant__stats__frv__SingleFinitePSpace(): - assert _test_args(SingleFinitePSpace(Symbol('x'), die)) - - -def test_diofant__stats__frv__ProductFinitePSpace(): - xp = SingleFinitePSpace(Symbol('x'), die) - yp = SingleFinitePSpace(Symbol('y'), die) - assert _test_args(ProductFinitePSpace(xp, yp)) - - -def test_diofant__stats__frv__SingleFiniteDistribution(): - pass - - -def test_diofant__stats__crv__ContinuousDistribution(): - pass - - -def test_diofant__stats__frv_types__FiniteDistributionHandmade(): - assert _test_args(FiniteDistributionHandmade({1: 1})) - - -def test_diofant__stats__crv__ContinuousDistributionHandmade(): - assert _test_args(ContinuousDistributionHandmade(Symbol('x'), - Interval(0, 2))) - - -def test_diofant__stats__rv__Density(): - assert _test_args(Density(Normal('x', 0, 1))) - - -def test_diofant__stats__crv_types__ArcsinDistribution(): - assert _test_args(ArcsinDistribution(0, 1)) - - -def test_diofant__stats__crv_types__BeniniDistribution(): - assert _test_args(BeniniDistribution(1, 1, 1)) - - -def test_diofant__stats__crv_types__BetaDistribution(): - assert _test_args(BetaDistribution(1, 1)) - - -def test_diofant__stats__crv_types__BetaPrimeDistribution(): - assert _test_args(BetaPrimeDistribution(1, 1)) - - -def test_diofant__stats__crv_types__CauchyDistribution(): - assert _test_args(CauchyDistribution(0, 1)) - - -def test_diofant__stats__crv_types__ChiDistribution(): - assert _test_args(ChiDistribution(1)) - - -def test_diofant__stats__crv_types__ChiNoncentralDistribution(): - assert _test_args(ChiNoncentralDistribution(1, 1)) - - -def test_diofant__stats__crv_types__ChiSquaredDistribution(): - assert _test_args(ChiSquaredDistribution(1)) - - -def test_diofant__stats__crv_types__DagumDistribution(): - assert _test_args(DagumDistribution(1, 1, 1)) - - -def test_diofant__stats__crv_types__ExponentialDistribution(): - assert _test_args(ExponentialDistribution(1)) - - -def test_diofant__stats__crv_types__FDistributionDistribution(): - assert _test_args(FDistributionDistribution(1, 1)) - - -def test_diofant__stats__crv_types__FisherZDistribution(): - assert _test_args(FisherZDistribution(1, 1)) - - -def test_diofant__stats__crv_types__FrechetDistribution(): - assert _test_args(FrechetDistribution(1, 1, 1)) - - -def test_diofant__stats__crv_types__GammaInverseDistribution(): - assert _test_args(GammaInverseDistribution(1, 1)) - - -def test_diofant__stats__crv_types__GammaDistribution(): - assert _test_args(GammaDistribution(1, 1)) - - -def test_diofant__stats__crv_types__KumaraswamyDistribution(): - assert _test_args(KumaraswamyDistribution(1, 1)) - - -def test_diofant__stats__crv_types__LaplaceDistribution(): - assert _test_args(LaplaceDistribution(0, 1)) - - -def test_diofant__stats__crv_types__LogisticDistribution(): - assert _test_args(LogisticDistribution(0, 1)) - - -def test_diofant__stats__crv_types__LogNormalDistribution(): - assert _test_args(LogNormalDistribution(0, 1)) - - -def test_diofant__stats__crv_types__MaxwellDistribution(): - assert _test_args(MaxwellDistribution(1)) - - -def test_diofant__stats__crv_types__NakagamiDistribution(): - assert _test_args(NakagamiDistribution(1, 1)) - - -def test_diofant__stats__crv_types__NormalDistribution(): - assert _test_args(NormalDistribution(0, 1)) - - -def test_diofant__stats__crv_types__ParetoDistribution(): - assert _test_args(ParetoDistribution(1, 1)) - - -def test_diofant__stats__crv_types__QuadraticUDistribution(): - assert _test_args(QuadraticUDistribution(1, 2)) - - -def test_diofant__stats__crv_types__RaisedCosineDistribution(): - assert _test_args(RaisedCosineDistribution(1, 1)) - - -def test_diofant__stats__crv_types__RayleighDistribution(): - assert _test_args(RayleighDistribution(1)) - - -def test_diofant__stats__crv_types__StudentTDistribution(): - assert _test_args(StudentTDistribution(1)) - - -def test_diofant__stats__crv_types__TriangularDistribution(): - assert _test_args(TriangularDistribution(-1, 0, 1)) - - -def test_diofant__stats__crv_types__UniformDistribution(): - assert _test_args(UniformDistribution(0, 1)) - - -def test_diofant__stats__crv_types__UniformSumDistribution(): - assert _test_args(UniformSumDistribution(1)) - - -def test_diofant__stats__crv_types__VonMisesDistribution(): - assert _test_args(VonMisesDistribution(1, 1)) - - -def test_diofant__stats__crv_types__WeibullDistribution(): - assert _test_args(WeibullDistribution(1, 1)) - - -def test_diofant__stats__crv_types__WignerSemicircleDistribution(): - assert _test_args(WignerSemicircleDistribution(1)) - - -def test_diofant__stats__drv_types__PoissonDistribution(): - assert _test_args(PoissonDistribution(1)) - - -def test_diofant__stats__drv_types__GeometricDistribution(): - assert _test_args(GeometricDistribution(.5)) - - def test_diofant__core__symbol__BaseSymbol(): pass diff --git a/diofant/tests/printing/test_latex.py b/diofant/tests/printing/test_latex.py index 2b691e48d8..debccb2a90 100644 --- a/diofant/tests/printing/test_latex.py +++ b/diofant/tests/printing/test_latex.py @@ -37,7 +37,6 @@ from diofant.printing.latex import (LatexPrinter, greek_letters_set, latex, other_symbols, tex_greek_dictionary, translate) -from diofant.stats import Die, Exponential, Normal, pspace, where from diofant.tensor import (ImmutableDenseNDimArray, ImmutableSparseNDimArray, MutableDenseNDimArray, MutableSparseNDimArray, tensorproduct) @@ -1104,20 +1103,6 @@ def test_latex_MatrixSlice(): r'X\left[5, :5:2\right]' -def test_latex_RandomDomain(): - X = Normal('x1', 0, 1) - assert latex(where(X > 0)) == r'Domain: 0 < x_{1} \wedge x_{1} < \infty' - - D = Die('d1', 6) - assert latex(where(D > 4)) == r'Domain: d_{1} = 5 \vee d_{1} = 6' - - A = Exponential('a', 1) - B = Exponential('b', 1) - assert latex( - pspace(Tuple(A, B)).domain) == \ - r'Domain: 0 \leq a \wedge 0 \leq b \wedge a < \infty \wedge b < \infty' - - def test_PrettyPoly(): R = QQ.inject(x, y) F = R.field diff --git a/diofant/tests/printing/test_pretty.py b/diofant/tests/printing/test_pretty.py index f3ed6560d2..5cfec725ea 100644 --- a/diofant/tests/printing/test_pretty.py +++ b/diofant/tests/printing/test_pretty.py @@ -26,7 +26,6 @@ from diofant.core.trace import Tr from diofant.printing.pretty import pretty as xpretty from diofant.printing.pretty.pretty_symbology import U, xobj -from diofant.stats import Die, Exponential, Normal, pspace, where from diofant.tensor import (ImmutableDenseNDimArray, ImmutableSparseNDimArray, MutableDenseNDimArray, MutableSparseNDimArray, tensorproduct) @@ -4941,19 +4940,6 @@ def test_elliptic_functions(): assert upretty(expr) == ucode_str -def test_RandomDomain(): - X = Normal('x1', 0, 1) - assert upretty(where(X > 0)) == 'Domain: 0 < x₁ ∧ x₁ < ∞' - - D = Die('d1', 6) - assert upretty(where(D > 4)) == 'Domain: d₁ = 5 ∨ d₁ = 6' - - A = Exponential('a', 1) - B = Exponential('b', 1) - assert upretty(pspace(Tuple(A, B)).domain) == \ - 'Domain: 0 ≤ a ∧ 0 ≤ b ∧ a < ∞ ∧ b < ∞' - - def test_PrettyPoly(): R = QQ.inject(x, y) F = R.field diff --git a/diofant/tests/printing/test_str.py b/diofant/tests/printing/test_str.py index 67247e98ab..aabefecd96 100644 --- a/diofant/tests/printing/test_str.py +++ b/diofant/tests/printing/test_str.py @@ -8,16 +8,15 @@ I, Integer, Integral, Interval, Lambda, Limit, Matrix, MatrixSymbol, Mul, Ne, O, Poly, Pow, Rational, Reals, Rel, RootOf, RootSum, S, SparseMatrix, StrPrinter, Sum, Symbol, - SymmetricDifference, Tuple, Wild, WildFunction, Xor, - ZeroMatrix, cbrt, cos, exp, factor, factorial, factorial2, - false, field, grlex, groebner, nan, oo, pi, ring, root, - sin, sqrt, sstr, sstrrepr, subfactorial, summation, - symbols, true, zeta, zoo) + SymmetricDifference, Wild, WildFunction, Xor, ZeroMatrix, + cbrt, cos, exp, factor, factorial, factorial2, false, + field, grlex, groebner, nan, oo, pi, ring, root, sin, + sqrt, sstr, sstrrepr, subfactorial, summation, symbols, + true, zeta, zoo) from diofant.abc import w, x, y, z from diofant.combinatorics import AbelianGroup, Cycle, Permutation from diofant.core.trace import Tr from diofant.geometry import Circle, Point -from diofant.stats import Die, Exponential, Normal, pspace, where from diofant.tensor.array import ImmutableDenseNDimArray @@ -677,18 +676,6 @@ def test_settings(): pytest.raises(TypeError, lambda: sstr(Integer(4), method='garbage')) -def test_RandomDomain(): - X = Normal('x1', 0, 1) - assert str(where(X > 0)) == 'Domain: (0 < x1) & (x1 < oo)' - - D = Die('d1', 6) - assert str(where(D > 4)) == 'Domain: Eq(d1, 5) | Eq(d1, 6)' - - A = Exponential('a', 1) - B = Exponential('b', 1) - assert str(pspace(Tuple(A, B)).domain) == 'Domain: (0 <= a) & (0 <= b) & (a < oo) & (b < oo)' - - def test_FiniteSet(): assert str(FiniteSet(*range(1, 51))) == '{1, 2, 3, ..., 48, 49, 50}' assert str(FiniteSet(*range(1, 6))) == '{1, 2, 3, 4, 5}' diff --git a/diofant/tests/stats/__init__.py b/diofant/tests/stats/__init__.py deleted file mode 100644 index c2eaad7d01..0000000000 --- a/diofant/tests/stats/__init__.py +++ /dev/null @@ -1 +0,0 @@ -"""Tests for stats package.""" diff --git a/diofant/tests/stats/test_continuous_rv.py b/diofant/tests/stats/test_continuous_rv.py deleted file mode 100644 index 13bf9be958..0000000000 --- a/diofant/tests/stats/test_continuous_rv.py +++ /dev/null @@ -1,665 +0,0 @@ -import pytest - -from diofant import (And, Eq, Integer, Integral, Interval, N, Piecewise, - Rational, Sum, Symbol, besseli, beta, binomial, cos, erf, - exp, expand_func, factorial, floor, gamma, log, - lowergamma, oo, pi, simplify, sin, sqrt, symbols) -from diofant.abc import k, x, y, z -from diofant.stats import (Arcsin, Benini, Beta, BetaPrime, Cauchy, Chi, - ChiNoncentral, ChiSquared, ContinuousRV, Dagum, E, - Erlang, Exponential, FDistribution, FisherZ, - Frechet, Gamma, GammaInverse, Kumaraswamy, Laplace, - Logistic, LogNormal, Maxwell, Nakagami, Normal, P, - Pareto, QuadraticU, RaisedCosine, Rayleigh, - StudentT, Triangular, Uniform, UniformSum, VonMises, - Weibull, WignerSemicircle, cdf, cmoment, - correlation, covariance, density, given, moment, - pspace, sample, skewness, smoment, variance, where) -from diofant.stats.crv_types import NormalDistribution -from diofant.stats.rv import ProductPSpace - - -__all__ = () - - -def test_single_normal(): - mu = Symbol('mu', real=True) - sigma = Symbol('sigma', real=True, positive=True) - X = Normal('x', 0, 1) - Y = X*sigma + mu - - assert simplify(E(Y)) == mu - assert simplify(variance(Y)) == sigma**2 - pdf = density(Y) - x = Symbol('x') - assert (pdf(x) == - sqrt(2)*exp(-(x - mu)**2/(2*sigma**2))/(2*sqrt(pi)*sigma)) - - assert P(X**2 < 1) == erf(sqrt(2)/2) - - assert E(X, Eq(X, mu)) == mu - - -def test_conditional_1d(): - X = Normal('x', 0, 1) - Y = given(X, X >= 0) - - assert density(Y)(x) == 2 * density(X)(x) - - assert Y.pspace.domain.set == Interval(0, oo, right_open=True) - assert E(Y) == sqrt(2/pi) - assert E(X**2) == E(Y**2) - - -def test_ContinuousDomain(): - X = Normal('x', 0, 1) - assert where(X**2 <= 1).set == Interval(-1, 1) - assert where(X**2 <= 1).symbol == X.symbol - assert where(And(X**2 <= 1, X >= 0)).set == Interval(0, 1) - pytest.raises(ValueError, lambda: where(sin(X) > 1)) - - Y = given(X, X >= 0) - - assert Y.pspace.domain.set == Interval(0, oo, False, True) - - -@pytest.mark.slow -def test_multiple_normal(): - X, Y = Normal('x', 0, 1), Normal('y', 0, 1) - - assert E(X + Y) == 0 - assert variance(X + Y) == 2 - assert variance(X + X) == 4 - assert covariance(X, Y) == 0 - assert covariance(2*X + Y, -X) == -2*variance(X) - assert skewness(X) == 0 - assert skewness(X + Y) == 0 - assert correlation(X, Y) == 0 - assert correlation(X, X + Y) == correlation(X, X - Y) - assert moment(X, 2) == 1 - assert cmoment(X, 3) == 0 - assert moment(X + Y, 4) == 12 - assert cmoment(X, 2) == variance(X) - assert smoment(X*X, 2) == 1 - assert smoment(X + Y, 3) == skewness(X + Y) - assert E(X, Eq(X + Y, 0)) == 0 - assert variance(X, Eq(X + Y, 0)) == Rational(1, 2) - - -@pytest.mark.slow -def test_symbolic(): - mu1, mu2 = symbols('mu1 mu2', real=True) - s1, s2 = symbols('sigma1 sigma2', real=True, positive=True) - rate = Symbol('lambda', real=True, positive=True) - X = Normal('x', mu1, s1) - Y = Normal('y', mu2, s2) - Z = Exponential('z', rate) - a, b = symbols('a b', real=True) - - assert E(X) == mu1 - assert E(X + Y) == mu1 + mu2 - assert E(a*X + b) == a*E(X) + b - assert variance(X) == s1**2 - assert simplify(variance(X + a*Y + b)) == variance(X) + a**2*variance(Y) - - assert E(Z) == 1/rate - assert E(a*Z + b) == a*E(Z) + b - assert E(X + a*Z + b) == mu1 + a/rate + b - - -def test_cdf(): - X = Normal('x', 0, 1) - - d = cdf(X) - assert P(X < 1) == d(1) - assert d(0) == Rational(1, 2) - - d = cdf(X, X > 0) # given X>0 - assert d(0) == 0 - - Y = Exponential('y', 10) - d = cdf(Y) - assert d(-5) == 0 - assert P(Y > 3) == 1 - d(3) - - pytest.raises(ValueError, lambda: cdf(X + Y)) - - Z = Exponential('z', 1) - f = cdf(Z) - z = Symbol('z') - assert f(z) == Piecewise((1 - exp(-z), z >= 0), (0, True)) - - U = Uniform('x', 3, 5) - u = cdf(U) - assert u(z) == z/2 - Rational(3, 2) - - -def test_sample(): - z = Symbol('z') - Z = ContinuousRV(z, exp(-z), set=Interval(0, oo)) - assert sample(Z) in Z.pspace.domain.set - sym, val = list(Z.pspace.sample().items())[0] - assert sym == Z - assert val in Interval(0, oo) - - -def test_ContinuousRV(): - x = Symbol('x') - pdf = sqrt(2)*exp(-x**2/2)/(2*sqrt(pi)) # Normal distribution - # X and Y should be equivalent - X = ContinuousRV(x, pdf) - Y = Normal('y', 0, 1) - - assert variance(X) == variance(Y) - assert P(X > 0) == P(Y > 0) - - -def test_arcsin(): - a = Symbol('a', extended_real=True) - b = Symbol('b', extended_real=True) - - X = Arcsin('x', a, b) - assert density(X)(x) == 1/(pi*sqrt((-x + b)*(x - a))) - - -def test_benini(): - alpha = Symbol('alpha', positive=True) - b = Symbol('beta', positive=True) - sigma = Symbol('sigma', positive=True) - - X = Benini('x', alpha, b, sigma) - assert density(X)(x) == ((alpha/x + 2*b*log(x/sigma)/x) - * exp(-alpha*log(x/sigma) - b*log(x/sigma)**2)) - assert X.pspace.domain.set == Interval(sigma, oo, False, True) - - -def test_beta(): - a, b = symbols('alpha beta', positive=True) - - B = Beta('x', a, b) - - assert pspace(B).domain.set == Interval(0, 1) - - dens = density(B) - x = Symbol('x') - assert dens(x) == x**(a - 1)*(1 - x)**(b - 1) / beta(a, b) - - # This is too slow - # assert E(B) == a / (a + b) - # assert variance(B) == (a*b) / ((a+b)**2 * (a+b+1)) - - # Full symbolic solution is too much, test with numeric version - a, b = Integer(1), Integer(2) - B = Beta('x', a, b) - assert expand_func(E(B)) == a/(a + b) - assert expand_func(variance(B)) == (a*b)/(a + b)**2/(a + b + 1) - - -def test_betaprime(): - alpha = Symbol('alpha', positive=True) - betap = Symbol('beta', positive=True) - - X = BetaPrime('x', alpha, betap) - assert density(X)(x) == x**(alpha - 1)*(x + 1)**(-alpha - betap)/beta(alpha, betap) - - -def test_cauchy(): - x0 = Symbol('x0') - gamma = Symbol('gamma', positive=True) - - X = Cauchy('x', x0, gamma) - assert density(X)(x) == 1/(pi*gamma*(1 + (x - x0)**2/gamma**2)) - - -def test_chi(): - k = Symbol('k', integer=True) - - X = Chi('x', k) - assert density(X)(x) == 2**(-k/2 + 1)*x**(k - 1)*exp(-x**2/2)/gamma(k/2) - - -def test_chi_noncentral(): - k = Symbol('k', integer=True) - l = Symbol('l') - - X = ChiNoncentral('x', k, l) - assert density(X)(x) == (x**k*l*(x*l)**(-k/2) * - exp(-x**2/2 - l**2/2)*besseli(k/2 - 1, x*l)) - - -def test_chi_squared(): - k = Symbol('k', integer=True) - - X = ChiSquared('x', k) - assert density(X)(x) == 2**(-k/2)*x**(k/2 - 1)*exp(-x/2)/gamma(k/2) - - -def test_dagum(): - p = Symbol('p', positive=True) - b = Symbol('b', positive=True) - a = Symbol('a', positive=True) - - X = Dagum('x', p, a, b) - assert density(X)(x) == a*p*(x/b)**(a*p)*((x/b)**a + 1)**(-p - 1)/x - - -def test_erlang(): - k = Symbol('k', integer=True, positive=True) - l = Symbol('l', positive=True) - - X = Erlang('x', k, l) - assert density(X)(x) == x**(k - 1)*l**k*exp(-x*l)/gamma(k) - - -def test_exponential(): - rate = Symbol('lambda', positive=True, real=True) - X = Exponential('x', rate) - - assert E(X) == 1/rate - assert variance(X) == 1/rate**2 - assert skewness(X) == 2 - assert skewness(X) == smoment(X, 3) - assert smoment(2*X, 4) == smoment(X, 4) - assert moment(X, 3) == 3*2*1/rate**3 - assert P(X > 0) == Integer(1) - assert P(X > 1) == exp(-rate) - assert P(X > 10) == exp(-10*rate) - - assert where(X <= 1).set == Interval(0, 1) - - # issue sympy/sympy#10003 - X = Exponential('x', 3) - G = Gamma('g', 1, 2) - assert P(X < -1) == 0 - assert P(G < -1) == 0 - - -def test_f_distribution(): - d1 = Symbol('d1', positive=True) - d2 = Symbol('d2', positive=True) - - X = FDistribution('x', d1, d2) - assert density(X)(x) == (d2**(d2/2)*sqrt((d1*x)**d1*(d1*x + d2)**(-d1 - d2)) - / (x*beta(d1/2, d2/2))) - - -def test_fisher_z(): - d1 = Symbol('d1', positive=True) - d2 = Symbol('d2', positive=True) - - X = FisherZ('x', d1, d2) - assert density(X)(x) == (2*d1**(d1/2)*d2**(d2/2) * - (d1*exp(2*x) + d2)**(-d1/2 - d2/2) * - exp(d1*x)/beta(d1/2, d2/2)) - - -def test_frechet(): - a = Symbol('a', positive=True) - s = Symbol('s', positive=True) - m = Symbol('m', extended_real=True) - - X = Frechet('x', a, s=s, m=m) - assert density(X)(x) == a*((x - m)/s)**(-a - 1)*exp(-((x - m)/s)**(-a))/s - - -def test_gamma(): - k = Symbol('k', positive=True) - theta = Symbol('theta', positive=True) - - X = Gamma('x', k, theta) - assert density(X)(x) == x**(k - 1)*theta**(-k)*exp(-x/theta)/gamma(k) - assert cdf(X, meijerg=True)(z) == Piecewise( - (-k*lowergamma(k, 0)/gamma(k + 1) + - k*lowergamma(k, z/theta)/gamma(k + 1), z >= 0), - (0, True)) - # assert simplify(variance(X)) == k*theta**2 # handled numerically below - assert E(X) == moment(X, 1) - - k, theta = symbols('k theta', real=True, positive=True) - X = Gamma('x', k, theta) - assert simplify(E(X)) == k*theta - # can't get things to simplify on this one so we use subs - assert variance(X).subs({k: 5}) == (k*theta**2).subs({k: 5}) - # The following is too slow - # assert simplify(skewness(X)).subs({k: 5}) == (2/sqrt(k)).subs({k: 5}) - - -def test_gamma_inverse(): - a = Symbol('a', positive=True) - b = Symbol('b', positive=True) - - X = GammaInverse('x', a, b) - assert density(X)(x) == x**(-a - 1)*b**a*exp(-b/x)/gamma(a) - - -def test_kumaraswamy(): - a = Symbol('a', positive=True) - b = Symbol('b', positive=True) - - X = Kumaraswamy('x', a, b) - assert density(X)(x) == x**(a - 1)*a*b*(-x**a + 1)**(b - 1) - - -def test_laplace(): - mu = Symbol('mu') - b = Symbol('b', positive=True) - - X = Laplace('x', mu, b) - assert density(X)(x) == exp(-abs(x - mu)/b)/(2*b) - - -def test_logistic(): - mu = Symbol('mu', extended_real=True) - s = Symbol('s', positive=True) - - X = Logistic('x', mu, s) - assert density(X)(x) == exp((-x + mu)/s)/(s*(exp((-x + mu)/s) + 1)**2) - - -def test_lognormal(): - # Right now, only density function and sampling works - # Test sampling: Only e^mean in sample std of 0 - for i in range(3): - X = LogNormal('x', i, 0) - assert sample(X) == N(exp(i)) - # The diofant integrator can't do this too well - # assert E(X) == - - mu = Symbol('mu', extended_real=True) - sigma = Symbol('sigma', positive=True) - - X = LogNormal('x', mu, sigma) - assert density(X)(x) == (sqrt(2)*exp(-(-mu + log(x))**2 - / (2*sigma**2))/(2*x*sqrt(pi)*sigma)) - - X = LogNormal('x', 0, 1) # Mean 0, standard deviation 1 - assert density(X)(x) == sqrt(2)*exp(-log(x)**2/2)/(2*x*sqrt(pi)) - - -@pytest.mark.xfail -def test_lognormal_xfail(): - mean = Symbol('mu', real=True) - std = Symbol('sigma', positive=True, real=True) - X = LogNormal('x', mean, std) - assert E(X) == exp(mean+std**2/2) - assert variance(X) == (exp(std**2)-1) * exp(2*mean + std**2) - - -def test_maxwell(): - a = Symbol('a', positive=True) - - X = Maxwell('x', a) - - assert density(X)(x) == (sqrt(2)*x**2*exp(-x**2/(2*a**2)) / - (sqrt(pi)*a**3)) - assert E(X) == 2*sqrt(2)*a/sqrt(pi) - assert simplify(variance(X)) == a**2*(-8 + 3*pi)/pi - - -def test_nakagami(): - mu = Symbol('mu', positive=True) - omega = Symbol('omega', positive=True) - - X = Nakagami('x', mu, omega) - assert density(X)(x) == (2*x**(2*mu - 1)*mu**mu*omega**(-mu) - * exp(-x**2*mu/omega)/gamma(mu)) - assert simplify(E(X, meijerg=True)) == (sqrt(mu)*sqrt(omega) - * gamma(mu + Rational(1, 2))/gamma(mu + 1)) - assert simplify(variance(X, meijerg=True)) == ( - omega - omega*gamma(mu + Rational(1, 2))**2/(gamma(mu)*gamma(mu + 1))) - - -def test_pareto(): - xm, beta = symbols('xm beta', positive=True, finite=True) - alpha = beta + 5 - X = Pareto('x', xm, alpha) - - dens = density(X) - x = Symbol('x') - assert dens(x) == x**(-(alpha + 1))*alpha*xm**alpha - - # These fail because Diofant can not deduce that 1/xm != 0 - # assert simplify(E(X)) == alpha*xm/(alpha-1) - # assert simplify(variance(X)) == xm**2*alpha / ((alpha-1)**2*(alpha-2)) - - -def test_pareto_numeric(): - xm, beta = Integer(3), Integer(2) - alpha = beta + 5 - X = Pareto('x', xm, alpha) - - assert E(X) == alpha*xm/(alpha - 1) - assert variance(X) == xm**2*alpha/((alpha - 1)**2*(alpha - 2)) - # Skewness tests too slow. Try shortcutting function? - - -def test_raised_cosine(): - mu = Symbol('mu', extended_real=True) - s = Symbol('s', positive=True) - - X = RaisedCosine('x', mu, s) - assert density(X)(x) == (Piecewise(((cos(pi*(x - mu)/s) + 1)/(2*s), - And(x <= mu + s, mu - s <= x)), (0, True))) - assert X.pspace.domain.set == Interval(mu - s, mu + s) - - -def test_rayleigh(): - sigma = Symbol('sigma', positive=True) - - X = Rayleigh('x', sigma) - assert density(X)(x) == x*exp(-x**2/(2*sigma**2))/sigma**2 - assert E(X) == sqrt(2)*sqrt(pi)*sigma/2 - assert variance(X) == -pi*sigma**2/2 + 2*sigma**2 - - -def test_studentt(): - nu = Symbol('nu', positive=True) - - X = StudentT('x', nu) - assert density(X)(x) == (1 + x**2/nu)**(-nu/2 - 1/2)/(sqrt(nu)*beta(1/2, nu/2)) - - -def test_triangular(): - a = Symbol('a') - b = Symbol('b') - c = Symbol('c') - - X = Triangular('x', a, b, c) - assert density(X)(x) == Piecewise( - ((2*x - 2*a)/((-a + b)*(-a + c)), And(a <= x, x < c)), - (2/(-a + b), Eq(x, c)), - ((-2*x + 2*b)/((-a + b)*(b - c)), And(x <= b, c < x)), - (0, True)) - - -def test_quadratic_u(): - a = Symbol('a', extended_real=True) - b = Symbol('b', extended_real=True) - - X = QuadraticU('x', a, b) - assert density(X)(x) == (Piecewise((12*(x - a/2 - b/2)**2/(-a + b)**3, - And(x <= b, a <= x)), (0, True))) - assert X.pspace.domain.set == Interval(a, b) - - -def test_uniform(): - l = Symbol('l', real=True) - w = Symbol('w', positive=True, finite=True) - X = Uniform('x', l, l + w) - - assert simplify(E(X)) == l + w/2 - assert simplify(variance(X)) == w**2/12 - - # With numbers all is well - X = Uniform('x', 3, 5) - assert P(X < 3) == 0 - assert P(X > 5) == 0 - assert P(X < 4) == P(X > 4) == Rational(1, 2) - - -def test_uniform_P(): - # This stopped working because SingleContinuousPSpace.compute_density no - # longer calls integrate on a DiracDelta but rather just solves directly. - # integrate used to call UniformDistribution.expectation which special-cased - # subsed out the Min and Max terms that Uniform produces. - # - # I decided to regress on this class for general cleanliness (and I suspect - # speed) of the algorithm. - l = Symbol('l', real=True) - w = Symbol('w', positive=True, finite=True) - X = Uniform('x', l, l + w) - assert P(X < l) == 0 - assert P(X > l + w) == 0 - - -def test_uniformsum(): - n = Symbol('n', integer=True) - - X = UniformSum('x', n) - assert X.pspace.domain.set == Interval(0, n) - - d = density(X)(z) - assert d == Sum((-1)**k*(z - k)**(n - 1)*binomial(n, k), (k, 0, floor(z)))/factorial(n - 1) - - -def test_von_mises(): - mu = Symbol('mu') - k = Symbol('k', positive=True) - - X = VonMises('x', mu, k) - assert density(X)(x) == exp(k*cos(x - mu))/(2*pi*besseli(0, k)) - - -def test_weibull(): - a, b = symbols('a b', positive=True, real=True) - X = Weibull('x', a, b) - - assert simplify(E(X)) == simplify(a * gamma(1 + 1/b)) - assert simplify(variance(X)) == simplify(a**2 * gamma(1 + 2/b) - E(X)**2) - # Skewness tests too slow. Try shortcutting function? - - -def test_weibull_numeric(): - # Test for integers and rationals - a = 1 - bvals = [Rational(1, 2), 1, Rational(3, 2), Integer(5)] - for b in bvals: - X = Weibull('x', a, b) - assert simplify(E(X)) == simplify(a * gamma(1 + 1/b)) - assert simplify(variance(X)) == simplify( - a**2 * gamma(1 + 2/b) - E(X)**2) - # Not testing Skew... it's slow with int/frac values > 3/2 - - -def test_wignersemicircle(): - R = Symbol('R', positive=True) - - X = WignerSemicircle('x', R) - assert density(X)(x) == 2*sqrt(-x**2 + R**2)/(pi*R**2) - assert E(X) == 0 - - -def test_prefab_sampling(): - N = Normal('X', 0, 1) - L = LogNormal('L', 0, 1) - E = Exponential('Ex', 1) - P = Pareto('P', 1, 3) - W = Weibull('W', 1, 1) - U = Uniform('U', 0, 1) - B = Beta('B', 2, 5) - G = Gamma('G', 1, 3) - - variables = [N, L, E, P, W, U, B, G] - niter = 10 - for var in variables: - for _ in range(niter): - assert sample(var) in var.pspace.domain.set - - -def test_input_value_assertions(): - a, b = symbols('a b') - p, q = symbols('p q', positive=True) - - pytest.raises(ValueError, lambda: Normal('x', 3, 0)) - pytest.raises(ValueError, lambda: Normal('x', a, b)) - Normal('X', a, p) # No error raised - pytest.raises(ValueError, lambda: Exponential('x', a)) - Exponential('Ex', p) # No error raised - for fn in [Pareto, Weibull, Beta, Gamma]: - pytest.raises(ValueError, lambda: fn('x', a, p)) - pytest.raises(ValueError, lambda: fn('x', p, a)) - fn('x', p, q) # No error raised - - -def test_unevaluated(): - X = Normal('x', 0, 1) - assert E(X, evaluate=False) == ( - Integral(sqrt(2)*x*exp(-x**2/2)/(2*sqrt(pi)), (x, -oo, oo))) - - assert E(X + 1, evaluate=False) == ( - Integral(sqrt(2)*x*exp(-x**2/2)/(2*sqrt(pi)), (x, -oo, oo)) + 1) - - -@pytest.mark.xfail -def test_unevaluated_fail(): - X = Normal('x', 0, 1) - assert P(X > 0, evaluate=False) == ( - Integral(sqrt(2)*exp(-x**2/2)/(2*sqrt(pi)), (x, 0, oo))) - - assert P(X > 0, X**2 < 1, evaluate=False) == ( - Integral(sqrt(2)*exp(-x**2/2)/(2*sqrt(pi) * - Integral(sqrt(2)*exp(-x**2/2)/(2*sqrt(pi)), - (x, -1, 1))), (x, 0, 1))) - - -def test_probability_unevaluated(): - T = Normal('T', 30, 3) - assert type(P(T > 33, evaluate=False)) == Integral - - -def test_density_unevaluated(): - X = Normal('X', 0, 1) - Y = Normal('Y', 0, 2) - assert isinstance(density(X+Y, evaluate=False)(z), Integral) - - -def test_NormalDistribution(): - nd = NormalDistribution(0, 1) - x = Symbol('x') - assert nd.cdf(x) == erf(sqrt(2)*x/2)/2 + Rational(1, 2) - assert isinstance(nd.sample(), float) or nd.sample().is_Number - assert nd.expectation(1, x) == 1 - assert nd.expectation(x, x) == 0 - assert nd.expectation(x**2, x) == 1 - - -def test_random_parameters(): - mu = Normal('mu', 2, 3) - meas = Normal('T', mu, 1) - assert density(meas, evaluate=False)(z) - assert isinstance(pspace(meas), ProductPSpace) - # assert density(meas, evaluate=False)(z) == Integral(mu.pspace.pdf * - # meas.pspace.pdf, (mu.symbol, -oo, oo)).subs({meas.symbol: z}) - - -def test_random_parameters_given(): - mu = Normal('mu', 2, 3) - meas = Normal('T', mu, 1) - assert given(meas, Eq(mu, 5)) == Normal('T', 5, 1) - - -def test_conjugate_priors(): - mu = Normal('mu', 2, 3) - x = Normal('x', mu, 1) - assert isinstance(simplify(density(mu, Eq(x, y), evaluate=False)(z)), - Integral) - - -def test_difficult_univariate(): - # Since using solve in place of deltaintegrate we're able to perform - # substantially more complex density computations on single continuous - # random variables. - x = Normal('x', 0, 1) - assert density(x**3) - assert density(exp(x**2)) - assert density(log(x)) diff --git a/diofant/tests/stats/test_discrete_rv.py b/diofant/tests/stats/test_discrete_rv.py deleted file mode 100644 index 9ac2ce19d3..0000000000 --- a/diofant/tests/stats/test_discrete_rv.py +++ /dev/null @@ -1,38 +0,0 @@ -from diofant import Rational, Sum, exp, factorial -from diofant.abc import x, z -from diofant.stats import E, density, variance -from diofant.stats.drv_types import (Geometric, GeometricDistribution, Poisson, - PoissonDistribution) - - -__all__ = () - - -def test_PoissonDistribution(): - l = 3 - p = PoissonDistribution(l) - assert abs(p.cdf(10).evalf() - 1) < .001 - assert p.expectation(x, x) == l - assert p.expectation(x**2, x) - p.expectation(x, x)**2 == l - - -def test_Poisson(): - l = 3 - x = Poisson('x', l) - assert E(x) == l - assert variance(x) == l - assert density(x) == PoissonDistribution(l) - assert isinstance(E(x, evaluate=False), Sum) - assert isinstance(E(2*x, evaluate=False), Sum) - assert density(x)(z) == 3**z/(exp(3)*factorial(z)) - - -def test_GeometricDistribution(): - p = Rational(1, 5) - d = GeometricDistribution(p) - assert d.expectation(x, x) == 1/p - assert d.expectation(x**2, x) - d.expectation(x, x)**2 == (1-p)/p**2 - assert abs(d.cdf(20000).evalf() - 1) < .001 - - X = Geometric('x', p) - assert density(X)(z) == Rational(4, 5)**(z - 1)/5 diff --git a/diofant/tests/stats/test_finite_rv.py b/diofant/tests/stats/test_finite_rv.py deleted file mode 100644 index 2e6a0c3775..0000000000 --- a/diofant/tests/stats/test_finite_rv.py +++ /dev/null @@ -1,302 +0,0 @@ -import pytest - -from diofant import (And, Dict, Eq, FiniteSet, Integer, Matrix, Or, Rational, - Symbol, Tuple, binomial, cos, pi, simplify, sqrt, symbols, - sympify) -from diofant.abc import p, x -from diofant.stats import (Bernoulli, Binomial, Coin, Die, DiscreteUniform, E, - FiniteRV, Hypergeometric, P, Rademacher, cdf, - cmoment, correlation, covariance, density, moment, - pspace, sample, skewness, smoment, variance, where) -from diofant.stats.frv_types import DieDistribution - - -__all__ = () - - -def BayesTest(A, B): - assert P(A, B) == P(And(A, B)) / P(B) - assert P(A, B) == P(B, A) * P(A) / P(B) - - -def test_discreteuniform(): - # Symbolic - a, b, c = symbols('a b c') - X = DiscreteUniform('X', [a, b, c]) - assert X.pspace.distribution.pdf(a) == Rational(1, 3) - assert X.pspace.distribution.pdf(p) == 0 - - assert E(X) == (a + b + c)/3 - assert simplify(variance(X) - - ((a**2 + b**2 + c**2)/3 - (a/3 + b/3 + c/3)**2)) == 0 - assert P(Eq(X, a)) == P(Eq(X, b)) == P(Eq(X, c)) == Rational(1, 3) - - Y = DiscreteUniform('Y', range(-5, 5)) - - # Numeric - assert E(Y) == -Rational(1, 2) - assert variance(Y) == Rational(33, 4) - - for x in range(-5, 5): - assert P(Eq(Y, x)) == Rational(1, 10) - assert P(Y <= x) == Rational(x + 6, 10) - assert P(Y >= x) == Rational(5 - x, 10) - - assert dict(density(Die('D', 6)).items()) == \ - dict(density(DiscreteUniform('U', range(1, 7))).items()) - - -def test_dice(): - # TODO: Make iid method! - X, Y, Z = Die('X', 6), Die('Y', 6), Die('Z', 6) - a, b = symbols('a b') - - assert E(X) == 3 + Rational(1, 2) - assert variance(X) == Rational(35, 12) - assert E(X + Y) == 7 - assert E(X + X) == 7 - assert E(a*X + b) == a*E(X) + b - assert variance(X + Y) == variance(X) + variance(Y) == cmoment(X + Y, 2) - assert variance(X + X) == 4 * variance(X) == cmoment(X + X, 2) - assert cmoment(X, 0) == 1 - assert cmoment(4*X, 3) == 64*cmoment(X, 3) - assert covariance(X, Y) == 0 - assert covariance(X, X + Y) == variance(X) - assert density(Eq(cos(X*pi), 1))[True] == Rational(1, 2) - assert correlation(X, Y) == 0 - assert correlation(X, Y) == correlation(Y, X) - assert smoment(X + Y, 3) == skewness(X + Y) - assert smoment(X, 0) == 1 - assert P(X > 3) == Rational(1, 2) - assert P(2*X > 6) == Rational(1, 2) - assert P(X > Y) == Rational(5, 12) - assert P(Eq(X, Y)) == P(Eq(X, 1)) - - assert E(X, X > 3) == 5 == moment(X, 1, 0, X > 3) - assert E(X, Y > 3) == E(X) == moment(X, 1, 0, Y > 3) - assert E(X + Y, Eq(X, Y)) == E(2*X) - assert moment(X, 0) == 1 - assert moment(5*X, 2) == 25*moment(X, 2) - - assert P(X > 3, X > 3) == 1 - assert P(X > Y, Eq(Y, 6)) == 0 - assert P(Eq(X + Y, 12)) == Rational(1, 36) - assert P(Eq(X + Y, 12), Eq(X, 6)) == Rational(1, 6) - - assert density(X + Y) == density(Y + Z) != density(X + X) - d = density(2*X + Y**Z) - assert d[22] == Rational(1, 108) - assert d[4100] == Rational(1, 216) - assert 3130 not in d - - assert pspace(X).domain.as_boolean() == Or( - *[Eq(X.symbol, i) for i in [1, 2, 3, 4, 5, 6]]) - - assert where(X > 3).set == FiniteSet(4, 5, 6) - - X = Die('X', 2) - x = X.symbol - - assert X.pspace.compute_cdf(X) == {1: Rational(1, 2), 2: 1} - assert X.pspace.sorted_cdf(X) == [(1, Rational(1, 2)), (2, 1)] - - assert X.pspace.compute_density(X)(1) == Rational(1, 2) - assert X.pspace.compute_density(X)(0) == 0 - assert X.pspace.compute_density(X)(8) == 0 - - assert X.pspace.density == x - - -def test_given(): - X = Die('X', 6) - assert density(X, X > 5) == {6: 1} - assert where(X > 2, X > 5).as_boolean() == Eq(X.symbol, 6) - assert sample(X, X > 5) == 6 - - -def test_domains(): - X, Y = Die('x', 6), Die('y', 6) - x, y = X.symbol, Y.symbol - # Domains - d = where(X > Y) - assert d.condition == (x > y) - d = where(And(X > Y, Y > 3)) - assert d.as_boolean() == Or(And(Eq(x, 5), Eq(y, 4)), And(Eq(x, 6), - Eq(y, 5)), And(Eq(x, 6), Eq(y, 4))) - assert len(d.elements) == 3 - - assert len(pspace(X + Y).domain.elements) == 36 - - Z = Die('x', 4) - - pytest.raises(ValueError, lambda: P(X > Z)) # Two domains with same internal symbol - - assert pspace(X + Y).domain.set == FiniteSet(1, 2, 3, 4, 5, 6)**2 - - assert where(X > 3).set == FiniteSet(4, 5, 6) - assert X.pspace.domain.dict == FiniteSet( - *[Dict({X.symbol: i}) for i in range(1, 7)]) - - assert where(X > Y).dict == FiniteSet(*[Dict({X.symbol: i, Y.symbol: j}) - for i in range(1, 7) for j in range(1, 7) if i > j]) - - X, Y = Die('x', 7), Die('y', 3) - x, y = X.symbol, Y.symbol - pytest.raises(ValueError, lambda: X.pspace.where(Y < 3)) - cset = X.pspace.where(X < 3) - assert ((x, 1),) in cset - assert ((x, 3),) not in cset - - assert X.pspace.where(True) == X.pspace.domain - - -def test_dice_bayes(): - X, Y, Z = Die('X', 6), Die('Y', 6), Die('Z', 6) - - BayesTest(X > 3, X + Y < 5) - BayesTest(Eq(X - Y, Z), Z > Y) - BayesTest(X > 3, X > 2) - - -def test_die_args(): - pytest.raises(ValueError, lambda: Die('X', -1)) # issue sympy/sympy#8105: negative sides. - pytest.raises(ValueError, lambda: Die('X', 0)) - pytest.raises(ValueError, lambda: Die('X', 1.5)) # issue sympy/sympy#8103: non integer sides. - - k = Symbol('k') - sym_die = Die('X', k) - pytest.raises(ValueError, lambda: density(sym_die).dict) - - -def test_bernoulli(): - p, a, b = symbols('p a b') - X = Bernoulli('B', p, a, b) - - assert E(X) == a*p + b*(-p + 1) - assert density(X)[a] == p - assert density(X)[b] == 1 - p - - X = Bernoulli('B', p, 1, 0) - - assert E(X) == p - assert simplify(variance(X)) == p*(1 - p) - assert E(a*X + b) == a*E(X) + b - assert simplify(variance(a*X + b)) == simplify(a**2 * variance(X)) - - -def test_cdf(): - D = Die('D', 6) - o = Integer(1) - - assert cdf( - D) == sympify({1: o/6, 2: o/3, 3: o/2, 4: 2*o/3, 5: 5*o/6, 6: o}) - - -def test_coins(): - C, D = Coin('C'), Coin('D') - H, T = symbols('H, T') - assert P(Eq(C, D)) == Rational(1, 2) - assert density(Tuple(C, D)) == {(H, H): Rational(1, 4), (H, T): Rational(1, 4), - (T, H): Rational(1, 4), (T, T): Rational(1, 4)} - assert dict(density(C).items()) == {H: Rational(1, 2), T: Rational(1, 2)} - - F = Coin('F', Rational(1, 10)) - assert P(Eq(F, H)) == Rational(1, 10) - - d = pspace(C).domain - - assert d.as_boolean() == Or(Eq(C.symbol, H), Eq(C.symbol, T)) - - pytest.raises(ValueError, lambda: P(C > D)) # Can't intelligently compare H to T - - -def test_binomial_verify_parameters(): - pytest.raises(ValueError, lambda: Binomial('b', .2, .5)) - pytest.raises(ValueError, lambda: Binomial('b', 3, 1.5)) - - -def test_binomial_numeric(): - nvals = range(5) - pvals = [0, Rational(1, 4), Rational(1, 2), Rational(3, 4), 1] - - for n in nvals: - for p in pvals: - X = Binomial('X', n, p) - assert E(X) == n*p - assert variance(X) == n*p*(1 - p) - if n > 0 and 0 < p < 1: - assert skewness(X) == (1 - 2*p)/sqrt(n*p*(1 - p)) - for k in range(n + 1): - assert P(Eq(X, k)) == binomial(n, k)*p**k*(1 - p)**(n - k) - - -@pytest.mark.slow -def test_binomial_symbolic(): - n = 10 # Because we're using for loops, can't do symbolic n - p = symbols('p', positive=True) - X = Binomial('X', n, p) - assert simplify(E(X)) == n*p == simplify(moment(X, 1)) - assert simplify(variance(X)) == n*p*(1 - p) == simplify(cmoment(X, 2)) - assert simplify(skewness(X) - (1-2*p)/sqrt(n*p*(1-p))) == 0 - - # Test ability to change success/failure winnings - H, T = symbols('H T') - Y = Binomial('Y', n, p, succ=H, fail=T) - assert simplify(E(Y) - (n*(H*p + T*(1 - p)))) == 0 - - -def test_hypergeometric_numeric(): - for N in range(1, 5): - for m in range(N + 1): - for n in range(1, N + 1): - X = Hypergeometric('X', N, m, n) - N, m, n = map(sympify, (N, m, n)) - assert sum(density(X).values()) == 1 - assert E(X) == n * m / N - if N > 1: - assert variance(X) == n*(m/N)*(N - m)/N*(N - n)/(N - 1) - # Only test for skewness when defined - if N > 2 and 0 < m < N and n < N: - assert skewness(X) == simplify((N - 2*m)*sqrt(N - 1)*(N - 2*n) - / (sqrt(n*m*(N - m)*(N - n))*(N - 2))) - - -def test_rademacher(): - X = Rademacher('X') - - assert E(X) == 0 - assert variance(X) == 1 - assert density(X)[-1] == Rational(1, 2) - assert density(X)[1] == Rational(1, 2) - - -def test_FiniteRV(): - F = FiniteRV('F', {1: Rational(1, 2), 2: Rational(1, 4), 3: Rational(1, 4)}) - - assert dict(density(F).items()) == {1: Rational(1, 2), 2: Rational(1, 4), 3: Rational(1, 4)} - assert P(F >= 2) == Rational(1, 2) - - assert pspace(F).domain.as_boolean() == Or( - *[Eq(F.symbol, i) for i in [1, 2, 3]]) - - -def test_density_call(): - x = Bernoulli('x', p) - d = density(x) - assert d(0) == 1 - p - assert d(5) == 0 - - assert 0 in d - assert 5 not in d - assert d(0) == d[0] - - -def test_DieDistribution(): - X = DieDistribution(6) - assert X.pdf(Rational(1, 2)) == 0 - assert X.pdf(x).subs({x: 1}).doit() == Rational(1, 6) - assert X.pdf(x).subs({x: 7}).doit() == 0 - assert X.pdf(x).subs({x: -1}).doit() == 0 - assert X.pdf(x).subs({x: Rational(1, 3)}).doit() == 0 - pytest.raises(TypeError, lambda: X.pdf(x).subs({x: Matrix([0, 0])})) - pytest.raises(ValueError, lambda: X.pdf(x**2 - 1)) diff --git a/diofant/tests/stats/test_mix.py b/diofant/tests/stats/test_mix.py deleted file mode 100644 index 12d711e583..0000000000 --- a/diofant/tests/stats/test_mix.py +++ /dev/null @@ -1,16 +0,0 @@ -from diofant import Eq, Symbol -from diofant.stats import Beta, Poisson -from diofant.stats.drv_types import PoissonDistribution -from diofant.stats.rv import ProductPSpace, density, pspace - - -__all__ = () - - -def test_density(): - x = Symbol('x') - l = Symbol('l', positive=True) - rate = Beta(l, 2, 3) - X = Poisson(x, rate) - assert isinstance(pspace(X), ProductPSpace) - assert density(X, Eq(rate, rate.symbol)) == PoissonDistribution(l) diff --git a/diofant/tests/stats/test_rv.py b/diofant/tests/stats/test_rv.py deleted file mode 100644 index ae70bf8eb1..0000000000 --- a/diofant/tests/stats/test_rv.py +++ /dev/null @@ -1,264 +0,0 @@ -import pytest - -from diofant import (And, Basic, Dict, DiracDelta, Eq, Interval, Rational, Sum, - Symbol, Tuple, cos, integrate, oo, sin, symbols) -from diofant.abc import x, z -from diofant.stats import (Die, E, Exponential, Normal, P, density, dependent, - given, independent, pspace, random_symbols, sample, - variance, where) -from diofant.stats.rv import (Density, NamedArgsMixin, ProductDomain, - ProductPSpace, RandomSymbol, rs_swap, - sample_iter) - - -__all__ = () - - -def test_where(): - X, Y = Die('X'), Die('Y') - Z = Normal('Z', 0, 1) - - assert where(Z**2 <= 1).set == Interval(-1, 1) - assert where( - Z**2 <= 1).as_boolean() == Interval(-1, 1).as_relational(Z.symbol) - assert where(And(X > Y, Y > 4)).as_boolean() == And( - Eq(X.symbol, 6), Eq(Y.symbol, 5)) - - assert len(where(X < 3).set) == 2 - assert 1 in where(X < 3).set - - X, Y = Normal('X', 0, 1), Normal('Y', 0, 1) - assert where(And(X**2 <= 1, X >= 0)).set == Interval(0, 1) - XX = given(X, And(X**2 <= 1, X >= 0)) - assert XX.pspace.domain.set == Interval(0, 1) - assert XX.pspace.domain.as_boolean() == \ - And(0 <= X.symbol, X.symbol**2 <= 1, -oo < X.symbol, X.symbol < oo) - - with pytest.raises(TypeError): - XX = given(X, X + 3) - - -def test_random_symbols(): - X, Y = Normal('X', 0, 1), Normal('Y', 0, 1) - - assert set(random_symbols(2*X + 1)) == {X} - assert set(random_symbols(2*X + Y)) == {X, Y} - assert set(random_symbols(2*X + Y.symbol)) == {X} - assert set(random_symbols(2)) == set() - - assert X.is_commutative - - -def test_pspace(): - X, Y = Normal('X', 0, 1), Normal('Y', 0, 1) - - pytest.raises(ValueError, lambda: pspace(5 + 3)) - pytest.raises(ValueError, lambda: pspace(x < 1)) - assert pspace(X) == X.pspace - assert pspace(2*X + 1) == X.pspace - assert pspace(2*X + Y) == ProductPSpace(Y.pspace, X.pspace) - - -def test_rs_swap(): - X = Normal('x', 0, 1) - Y = Exponential('y', 1) - - XX = Normal('x', 0, 2) - YY = Normal('y', 0, 3) - - expr = 2*X + Y - assert expr.subs(rs_swap((X, Y), (YY, XX))) == 2*XX + YY - - -def test_RandomSymbol(): - pytest.raises(TypeError, lambda: RandomSymbol(0, 1)) - pytest.raises(TypeError, lambda: RandomSymbol(0, Symbol('x'))) - - X = Normal('x', 0, 1) - Y = Normal('x', 0, 2) - assert X.symbol == Y.symbol - assert X != Y - - assert X.name == X.symbol.name - - Normal('lambda', 0, 1) # make sure we can use protected terms - X = Normal('Lambda', 0, 1) # make sure we can use Diofant terms - - pytest.raises(TypeError, lambda: Normal(1, 0, 1)) - - -def test_RandomSymbol_diff(): - X = Normal('x', 0, 1) - assert (2*X).diff(X) - - -def test_overlap(): - X = Normal('x', 0, 1) - Y = Normal('x', 0, 2) - - pytest.raises(ValueError, lambda: P(X > Y)) - - -def test_ProductPSpace(): - X = Normal('X', 0, 1) - Y = Normal('Y', 0, 1) - px = X.pspace - py = Y.pspace - assert pspace(X + Y) == ProductPSpace(px, py) - assert pspace(X + Y) == ProductPSpace(py, px) - - X = Die('X', 2) - Y = Die('Y', 2) - - assert (pspace(X + Y).density == - Dict((frozenset({('X', 1), ('Y', 2)}), Rational(1, 4)), - (frozenset({('X', 1), ('Y', 1)}), Rational(1, 4)), - (frozenset({('X', 2), ('Y', 1)}), Rational(1, 4)), - (frozenset({('X', 2), ('Y', 2)}), Rational(1, 4)))) - d = pspace(X + Y).domain - assert ((X.symbol, 1), (Y.symbol, 2)) in d - assert ((X.symbol, 0), (Y.symbol, 2)) not in d - - Z = Die('Z', 2) - d1 = pspace(X + Y).domain - assert ProductDomain(d1, Z.pspace.domain) == pspace(X + Y + Z).domain - - -def test_E(): - assert E(5) == 5 - - -def test_Sample(): - X = Die('X', 6) - Y = Normal('Y', 0, 1) - - assert sample(X) in [1, 2, 3, 4, 5, 6] - assert sample(X + Y).is_Float - - assert P(X + Y > 0, Y < 0, numsamples=10).is_number - assert P(X > 10, numsamples=10).is_number - assert E(X + Y, numsamples=10).is_number - assert variance(X + Y, numsamples=10).is_number - - pytest.raises(TypeError, lambda: P(Y > z, numsamples=5)) - - assert P(sin(Y) <= 1, numsamples=10, modules=['math']) == 1 - assert P(sin(Y) <= 1, cos(Y) < 1, numsamples=10, modules=['math']) == 1 - - assert all(i in range(1, 7) for i in density(X, numsamples=10)) - assert all(i in range(4, 7) for i in density(X, X > 3, numsamples=10)) - - # Make sure this doesn't raise an error - Y = Normal('Y', 0, 1) - E(Sum(1/z**Y, (z, 1, oo)), Y > 2, numsamples=3, modules='mpmath') - - -def test_given(): - X = Normal('X', 0, 1) - Y = Normal('Y', 0, 1) - A = given(X, True) - B = given(X, Y > 2) - - assert X == A == B - - -def test_dependence(): - X, Y = Die('X'), Die('Y') - assert independent(X, 2*Y) - assert not dependent(X, 2*Y) - - X, Y = Normal('X', 0, 1), Normal('Y', 0, 1) - assert independent(X, Y) - assert dependent(X, 2*X) - - # Create a dependency - XX, YY = given(Tuple(X, Y), Eq(X + Y, 3)) - assert dependent(XX, YY) - - -@pytest.mark.xfail -def test_dependent_finite(): - X, Y = Die('X'), Die('Y') - # Dependence testing requires symbolic conditions which currently break - # finite random variables - assert dependent(X, Y + X) - - XX, YY = given(Tuple(X, Y), X + Y > 5) # Create a dependency - assert dependent(XX, YY) - - -def test_normality(): - X, Y = Normal('X', 0, 1), Normal('Y', 0, 1) - x, z = symbols('x, z', real=True) - dens = density(X - Y, Eq(X + Y, z)) - - assert integrate(dens(x), (x, -oo, oo)) == 1 - - -def test_Density(): - X = Die('X', 6) - d = Density(X) - assert d.doit() == density(X) - r = Rational(1, 6) - assert density(2*X).dict == {2: r, 4: r, 6: r, 8: r, 10: r, 12: r} - - X = Normal('X', 0, 1) - Y = Normal('Y', 0, 1) - assert str(density(X - Y, Y)(z)) == 'sqrt(2)*E**(-(Y + z)**2/2)/(2*sqrt(pi))' - - -def test_NamedArgsMixin(): - class Foo(Basic, NamedArgsMixin): - _argnames = 'foo', 'bar' - - a = Foo(1, 2) - - assert a.foo == 1 - assert a.bar == 2 - - pytest.raises(AttributeError, lambda: a.baz) - - class Bar(Basic, NamedArgsMixin): - pass - - pytest.raises(AttributeError, lambda: Bar(1, 2).foo) - - -def test_density_constant(): - assert density(3)(2) == 0 - assert density(3)(3) == DiracDelta(0) - - -def test_real(): - x = Normal('x', 0, 1) - assert x.is_extended_real - - -def test_exponential(): - # issue sympy/sympy#10052: - X = Exponential('X', 3) - assert P(X < oo) == 1 - assert P(X > oo) == 0 - assert P(X < 2, X > oo) == 0 - assert P(X < oo, X > oo) == 0 - assert P(X < oo, X > 2) == 1 - assert P(X < 3, X == 2) == 0 - pytest.raises(ValueError, lambda: P(1)) - pytest.raises(ValueError, lambda: P(X < 1, 2)) - - -def test_sample_iter(): - X = Normal('X', 0, 1) - expr = X*X + 3 - assert all(_ > 3 for _ in sample_iter(expr, numsamples=3)) - assert all(_ > 5 for _ in sample_iter(expr, numsamples=3, condition=X > 2)) - pytest.raises(ValueError, lambda: tuple(sample_iter(expr, numsamples=3, - condition=X))) - - -def test_sympyissue_8129(): - # pylint: disable=comparison-with-itself - X = Exponential('X', 4) - assert P(X >= X) == 1 - assert P(X > X) == 0 - assert P(X > X+1) == 0 diff --git a/docs/modules/index.rst b/docs/modules/index.rst index 75ee29b4e7..72fad71e7d 100644 --- a/docs/modules/index.rst +++ b/docs/modules/index.rst @@ -37,7 +37,6 @@ automatically generated using Diofant's docstrings. plotting sets simplify/simplify - stats solvers/index tensor/index utilities/index diff --git a/docs/modules/stats.rst b/docs/modules/stats.rst deleted file mode 100644 index 22803f23da..0000000000 --- a/docs/modules/stats.rst +++ /dev/null @@ -1,182 +0,0 @@ -Stats -===== - -.. automodule:: diofant.stats - -Random Variable Types -^^^^^^^^^^^^^^^^^^^^^ - -Finite Types ---------------- -.. autofunction:: DiscreteUniform -.. autofunction:: Die -.. autofunction:: Bernoulli -.. autofunction:: Coin -.. autofunction:: Binomial -.. autofunction:: Hypergeometric -.. autofunction:: FiniteRV - -Discrete Types ------------------ -.. autofunction:: Geometric -.. autofunction:: Poisson - -Continuous Types -------------------- - -.. autofunction:: Arcsin -.. autofunction:: Benini -.. autofunction:: Beta -.. autofunction:: BetaPrime -.. autofunction:: Cauchy -.. autofunction:: Chi -.. autofunction:: ChiNoncentral -.. autofunction:: ChiSquared -.. autofunction:: Dagum -.. autofunction:: Erlang -.. autofunction:: Exponential -.. autofunction:: FDistribution -.. autofunction:: FisherZ -.. autofunction:: Frechet -.. autofunction:: Gamma -.. autofunction:: GammaInverse -.. autofunction:: Kumaraswamy -.. autofunction:: Laplace -.. autofunction:: Logistic -.. autofunction:: LogNormal -.. autofunction:: Maxwell -.. autofunction:: Nakagami -.. autofunction:: Normal -.. autofunction:: Pareto -.. autofunction:: QuadraticU -.. autofunction:: RaisedCosine -.. autofunction:: Rayleigh -.. autofunction:: StudentT -.. autofunction:: Triangular -.. autofunction:: Uniform -.. autofunction:: UniformSum -.. autofunction:: VonMises -.. autofunction:: Weibull -.. autofunction:: WignerSemicircle -.. autofunction:: ContinuousRV - -Interface -^^^^^^^^^ - -.. autofunction:: P -.. autofunction:: E -.. autofunction:: density -.. autofunction:: given -.. autofunction:: where -.. autofunction:: variance -.. autofunction:: std -.. autofunction:: sample -.. autofunction:: sample_iter - -Mechanics -^^^^^^^^^ -.. module:: diofant.stats.rv - -Diofant Stats employs a relatively complex class hierarchy. - -``RandomDomain``\s are a mapping of variables to possible values. For example we -might say that the symbol ``Symbol('x')`` can take on the values -`\{1,2,3,4,5,6\}`. - -.. class:: RandomDomain - -A ``PSpace``, or Probability Space, combines a ``RandomDomain`` with a density to -provide probabilistic information. For example the above domain could be -enhanced by a finite density ``{1:1/6, 2:1/6, 3:1/6, 4:1/6, 5:1/6, 6:1/6}`` to -fully define the roll of a fair die named ``x``. - -.. class:: PSpace - -A RandomSymbol represents the PSpace's symbol 'x' inside of Diofant expressions. - -.. class:: RandomSymbol - -The RandomDomain and PSpace classes are almost never directly instantiated. -Instead they are subclassed for a variety of situations. - -RandomDomains and PSpaces must be sufficiently general to represent domains and -spaces of several variables with arbitrarily complex densities. This generality -is often unnecessary. Instead we often build SingleDomains and SinglePSpaces to -represent single, univariate events and processes such as a single die or a -single normal variable. - -.. class:: SinglePSpace -.. class:: SingleDomain - - -Another common case is to collect together a set of such univariate random -variables. A collection of independent SinglePSpaces or SingleDomains can be -brought together to form a ProductDomain or ProductPSpace. These objects would -be useful in representing three dice rolled together for example. - -.. class:: ProductDomain - -.. class:: ProductPSpace - -The Conditional adjective is added whenever we add a global condition to a -RandomDomain or PSpace. A common example would be three independent dice where -we know their sum to be greater than 12. - -.. class:: ConditionalDomain - -We specialize further into Finite and Continuous versions of these classes to -represent finite (such as dice) and continuous (such as normals) random -variables. - -.. module:: diofant.stats.frv -.. class:: FiniteDomain -.. class:: FinitePSpace - -.. module:: diofant.stats.crv -.. class:: ContinuousDomain -.. class:: ContinuousPSpace - -Additionally there are a few specialized classes that implement certain common -random variable types. There is for example a DiePSpace that implements -SingleFinitePSpace and a NormalPSpace that implements SingleContinuousPSpace. - -.. module:: diofant.stats.frv_types -.. class:: DiePSpace - -.. module:: diofant.stats.crv_types -.. class:: NormalPSpace - -RandomVariables can be extracted from these objects using the PSpace.values -method. - -As previously mentioned Diofant Stats employs a relatively complex class -structure. Inheritance is widely used in the implementation of end-level -classes. This tactic was chosen to balance between the need to allow Diofant to -represent arbitrarily defined random variables and optimizing for common cases. -This complicates the code but is structured to only be important to those -working on extending Diofant Stats to other random variable types. - -Users will not use this class structure. Instead these mechanics are exposed -through variable creation functions Die, Coin, FiniteRV, Normal, Exponential, -etc.... These build the appropriate SinglePSpaces and return the corresponding -RandomVariable. Conditional and Product spaces are formed in the natural -construction of Diofant expressions and the use of interface functions E, Given, -Density, etc.... - - -.. function:: diofant.stats.Die -.. function:: diofant.stats.Normal - -There are some additional functions that may be useful. They are largely used -internally. - - -.. autofunction:: diofant.stats.rv.random_symbols -.. autofunction:: diofant.stats.rv.pspace -.. autofunction:: diofant.stats.rv.rs_swap - -.. autofunction:: diofant.stats.rv.sampling_P -.. autofunction:: diofant.stats.rv.sampling_E -.. autofunction:: diofant.stats.rv.sampling_density -.. autofunction:: diofant.stats.rv.independent -.. autofunction:: diofant.stats.rv.dependent diff --git a/docs/release/notes-0.10.rst b/docs/release/notes-0.10.rst index d03a37e5e4..22dde57e63 100644 --- a/docs/release/notes-0.10.rst +++ b/docs/release/notes-0.10.rst @@ -81,7 +81,7 @@ Minor changes * Added ``FALLBACK_GCD_ZZ_METHOD`` configuration option to specify GCD algorithm for polynomials with integer coefficients if heuristic GCD was off or just unlucky, see :pull:`721`. * Added ``GCD_AA_METHOD`` configuration option to specify GCD algorithm for polynomials with algebraic coefficients, see :pull:`721`. * :meth:`~diofant.polys.polytools.Poly.sqf_part`, :meth:`~diofant.polys.polytools.Poly.sqf_norm`, :meth:`~diofant.polys.polytools.Poly.sqf_list` methods and :attr:`~diofant.polys.polytools.Poly.is_squarefree` property use notion of being square-free w.r.t. to all polynomial variables, see :pull:`726`. -* 100% test coverage for :mod:`~diofant.core`, :mod:`~diofant.polys` and :mod:`~diofant.stats` modules. Overall test coverage is around 97%. +* 100% test coverage for :mod:`~diofant.core`, :mod:`~diofant.polys` and ``stats`` modules. Overall test coverage is around 97%. Developer changes ================= diff --git a/docs/release/notes-0.14.rst b/docs/release/notes-0.14.rst index 7a22ec6d06..9cf2de86a1 100644 --- a/docs/release/notes-0.14.rst +++ b/docs/release/notes-0.14.rst @@ -26,6 +26,7 @@ Compatibility breaks * Removed ``diofant.calculus.euler`` and ``diofant.calculus.finite_diff`` modules, see :pull:`1271`. * Removed ``diofant.vector`` module, see :pull:`1274`. * Removed ``diofant.diffgeom`` module, see :pull:`1281`. +* Removed ``diofant.stats`` module, see :pull:`1276`. Minor changes =============