diff --git a/doc/src/modules/stats.rst b/doc/src/modules/stats.rst index 00352445ddb6..7ea6045312b4 100644 --- a/doc/src/modules/stats.rst +++ b/doc/src/modules/stats.rst @@ -25,9 +25,17 @@ Continuous Types .. 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 @@ -35,11 +43,14 @@ Continuous Types .. 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 diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py index c660e444b518..1d7aa0aa5d7a 100644 --- a/sympy/core/tests/test_args.py +++ b/sympy/core/tests/test_args.py @@ -696,6 +696,16 @@ def test_sympy__stats__crv_types__ChiPSpace(): assert _test_args(ChiPSpace('X', 1)) +def test_sympy__stats__crv_types__ChiNoncentralPSpace(): + from sympy.stats.crv_types import ChiNoncentralPSpace + assert _test_args(ChiNoncentralPSpace('X', 1,1)) + + +def test_sympy__stats__crv_types__ChiSquaredPSpace(): + from sympy.stats.crv_types import ChiSquaredPSpace + assert _test_args(ChiSquaredPSpace('X', 1)) + + def test_sympy__stats__crv_types__DagumPSpace(): from sympy.stats.crv_types import DagumPSpace assert _test_args(DagumPSpace('X', 1, 1, 1)) @@ -706,11 +716,35 @@ def test_sympy__stats__crv_types__ExponentialPSpace(): assert _test_args(ExponentialPSpace('X', 1)) +def test_sympy__stats__crv_types__FDistributionPSpace(): + from sympy.stats.crv_types import FDistributionPSpace + assert _test_args(FDistributionPSpace('X', 1, 1)) + + +def test_sympy__stats__crv_types__FisherZPSpace(): + from sympy.stats.crv_types import FisherZPSpace + assert _test_args(FisherZPSpace('X', 1, 1)) + + +def test_sympy__stats__crv_types__FrechetPSpace(): + from sympy.stats.crv_types import FrechetPSpace + assert _test_args(FrechetPSpace('X', 1, 1, 1)) + + +def test_sympy__stats__crv_types__GammaInversePSpace(): + from sympy.stats.crv_types import GammaInversePSpace + assert _test_args(GammaInversePSpace('X', 1, 1)) + + def test_sympy__stats__crv_types__GammaPSpace(): from sympy.stats.crv_types import GammaPSpace assert _test_args(GammaPSpace('X', 1, 1)) +def test_sympy__stats__crv_types__KumaraswamyPSpace(): + from sympy.stats.crv_types import KumaraswamyPSpace + assert _test_args(KumaraswamyPSpace('X', 1, 1)) + def test_sympy__stats__crv_types__LaplacePSpace(): from sympy.stats.crv_types import LaplacePSpace assert _test_args(LaplacePSpace('X', 0, 1)) @@ -746,6 +780,14 @@ def test_sympy__stats__crv_types__ParetoPSpace(): assert _test_args(ParetoPSpace('X', 1, 1)) +def test_sympy__stats__crv_types__QuadraticUPSpace(): + from sympy.stats.crv_types import QuadraticUPSpace + assert _test_args(QuadraticUPSpace('X', 1, 2)) + +def test_sympy__stats__crv_types__RaisedCosinePSpace(): + from sympy.stats.crv_types import RaisedCosinePSpace + assert _test_args(RaisedCosinePSpace('X', 1, 1)) + def test_sympy__stats__crv_types__RayleighPSpace(): from sympy.stats.crv_types import RayleighPSpace assert _test_args(RayleighPSpace('X', 1)) @@ -771,6 +813,11 @@ def test_sympy__stats__crv_types__UniformSumPSpace(): assert _test_args(UniformSumPSpace('X', 1)) +def test_sympy__stats__crv_types__VonMisesPSpace(): + from sympy.stats.crv_types import VonMisesPSpace + assert _test_args(VonMisesPSpace('X', 1, 1)) + + def test_sympy__stats__crv_types__WeibullPSpace(): from sympy.stats.crv_types import WeibullPSpace assert _test_args(WeibullPSpace('X', 1, 1)) diff --git a/sympy/stats/__init__.py b/sympy/stats/__init__.py index 8039c899a8a9..ca2bd5c9f34d 100644 --- a/sympy/stats/__init__.py +++ b/sympy/stats/__init__.py @@ -53,9 +53,12 @@ import crv_types from crv_types import ( - Arcsin, Benini, Beta, BetaPrime, Cauchy, Chi, ContinuousRV, Dagum, - Exponential, Gamma, Laplace, Logistic, LogNormal, Maxwell, Nakagami, - Normal, Pareto, Rayleigh, StudentT, Triangular, Uniform, UniformSum, - Weibull, WignerSemicircle, + 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 ) __all__.extend(crv_types.__all__) diff --git a/sympy/stats/crv_types.py b/sympy/stats/crv_types.py index 0c1dabf6612f..9a5e6511e5eb 100644 --- a/sympy/stats/crv_types.py +++ b/sympy/stats/crv_types.py @@ -9,9 +9,17 @@ BetaPrime Cauchy Chi +ChiNoncentral +ChiSquared Dagum +Erlang Exponential +FDistribution +FisherZ +Frechet Gamma +GammaInverse +Kumaraswamy Laplace Logistic LogNormal @@ -19,19 +27,23 @@ Nakagami Normal Pareto +QuadraticU +RaisedCosine Rayleigh StudentT Triangular Uniform UniformSum +VonMises Weibull WignerSemicircle """ from sympy import (exp, log, sqrt, pi, S, Dummy, Interval, S, sympify, gamma, Piecewise, And, Eq, binomial, factorial, Sum, floor, Abs, - Symbol, log) + Symbol, log, besseli) from sympy import beta as beta_fn +from sympy import cos, exp, besseli from crv import SingleContinuousPSpace from sympy.core.decorators import _sympifyit import random @@ -45,9 +57,17 @@ 'BetaPrime', 'Cauchy', 'Chi', +'ChiNoncentral', +'ChiSquared', 'Dagum', +'Erlang', 'Exponential', +'FDistribution', +'FisherZ', +'Frechet', 'Gamma', +'GammaInverse', +'Kumaraswamy', 'Laplace', 'Logistic', 'LogNormal', @@ -55,11 +75,14 @@ 'Nakagami', 'Normal', 'Pareto', +'QuadraticU', +'RaisedCosine', 'Rayleigh', 'StudentT', 'Triangular', 'Uniform', 'UniformSum', +'VonMises', 'Weibull', 'WignerSemicircle' ] @@ -487,6 +510,128 @@ def Chi(name, k): return ChiPSpace(name, k).value +#------------------------------------------------------------------------------- +# Non-central Chi distribution ------------------------------------------------- + + +class ChiNoncentralPSpace(SingleContinuousPSpace): + def __new__(cls, name, k, l): + k = sympify(k) + l = sympify(l) + x = Symbol(name) + pdf = exp(-(x**2+l**2)/2)*x**k*l / (l*x)**(k/2) * besseli(k/2-1, l*x) + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set = Interval(0, oo)) + return obj + + +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 :math:`x \geq 0`. + + Parameters + ========== + + k : `k` > 0 the number of degrees of freedom + l : Shift parameter + + Returns + ======= + + A RandomSymbol. + + Examples + ======== + + >>> from sympy.stats import ChiNoncentral, density, E, std + >>> from sympy import Symbol, simplify + + >>> k = Symbol("k", integer=True) + >>> l = Symbol("l") + + >>> X = ChiNoncentral("x", k, l) + + >>> density(X) + Lambda(_x, _x**k*l*(_x*l)**(-k/2)*exp(-_x**2/2 - l**2/2)*besseli(k/2 - 1, _x*l)) + + References + ========== + + [1] http://en.wikipedia.org/wiki/Noncentral_chi_distribution + """ + + return ChiNoncentralPSpace(name, k, l).value + +#------------------------------------------------------------------------------- +# Chi squared distribution ----------------------------------------------------- + + +class ChiSquaredPSpace(SingleContinuousPSpace): + def __new__(cls, name, k): + k = sympify(k) + x = Symbol(name) + pdf = 1/(2**(k/2)*gamma(k/2))*x**(k/2 - 1)*exp(-x/2) + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(0, oo)) + return obj + + +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 :math:`x \geq 0`. + + Parameters + ========== + + k : Integer, `k` > 0 the number of degrees of freedom + + Returns + ======= + + A RandomSymbol. + + Examples + ======== + + >>> from sympy.stats import ChiSquared, density, E, variance + >>> from sympy import Symbol, simplify, combsimp, expand_func + + >>> k = Symbol("k", integer=True, positive=True) + + >>> X = ChiSquared("x", k) + + >>> density(X) + Lambda(_x, 2**(-k/2)*_x**(k/2 - 1)*exp(-_x/2)/gamma(k/2)) + + >>> combsimp(E(X)) + k + + >>> simplify(expand_func(variance(X))) + 2*k + + References + ========== + + [1] http://en.wikipedia.org/wiki/Chi_squared_distribution + [2] http://mathworld.wolfram.com/Chi-SquaredDistribution.html + """ + + return ChiSquaredPSpace(name, k).value + #------------------------------------------------------------------------------- # Dagum distribution ----------------------------------------------------------- @@ -547,6 +692,72 @@ def Dagum(name, p, a, b): return DagumPSpace(name, p, a, b).value +#------------------------------------------------------------------------------- +# 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 :math:`x \in [0,\infty]`. + + Parameters + ========== + + k : Integer + l : Real number, :math:`\lambda` > 0 the rate + + Returns + ======= + + A RandomSymbol. + + Examples + ======== + + >>> from sympy.stats import Erlang, density, cdf, E, variance + >>> from sympy import Symbol, simplify, pprint + + >>> k = Symbol("k", integer=True, positive=True) + >>> l = Symbol("l", positive=True) + + >>> X = Erlang("x", k, l) + + >>> D = density(X) + >>> pprint(D, use_unicode=False) + / k - 1 k -x*l\ + | x *l *e | + Lambda|x, ---------------| + \ gamma(k) / + + >>> C = cdf(X, meijerg=True) + >>> pprint(C, use_unicode=False) + / / k*lowergamma(k, 0) k*lowergamma(k, z*l) \ + Lambda|z, |- ------------------ + -------------------- for z >= 0| + | < gamma(k + 1) gamma(k + 1) | + | | | + \ \ 0 otherwise / + + >>> simplify(E(X)) + k/l + + >>> simplify(variance(X)) + (gamma(k)*gamma(k + 2) - gamma(k + 1)**2)/(l**2*gamma(k)**2) + + References + ========== + + [1] http://en.wikipedia.org/wiki/Erlang_distribution + [2] http://mathworld.wolfram.com/ErlangDistribution.html + """ + + return GammaPSpace(name, k, 1/l).value + #------------------------------------------------------------------------------- # Exponential distribution ----------------------------------------------------- @@ -635,6 +846,202 @@ def Exponential(name, rate): return ExponentialPSpace(name, rate).value +#------------------------------------------------------------------------------- +# F distribution --------------------------------------------------------------- + +class FDistributionPSpace(SingleContinuousPSpace): + def __new__(cls, name, d1, d2): + d1 = sympify(d1) + d2 = sympify(d2) + x = Symbol(name) + pdf = (sqrt((d1*x)**d1*d2**d2 / (d1*x+d2)**(d1+d2)) + / (x * beta_fn(d1/2, d2/2))) + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(0, oo)) + return obj + +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 :math:`x > 0`. + + Parameters + ========== + + d1 : `d1` > 0 a parameter + d2 : `d2` > 0 a parameter + + Returns + ======= + + A RandomSymbol. + + Examples + ======== + + >>> from sympy.stats import FDistribution, density + >>> from sympy import Symbol, simplify, pprint + + >>> d1 = Symbol("d1", positive=True) + >>> d2 = Symbol("d2", positive=True) + + >>> X = FDistribution("x", d1, d2) + + >>> D = density(X) + >>> pprint(D, use_unicode=False) + / d2 \ + | -- ______________________________ | + | 2 / d1 -d1 - d2 /d1 d2\| + | d2 *\/ (x*d1) *(x*d1 + d2) *gamma|-- + --|| + | \2 2 /| + Lambda|x, -----------------------------------------------------| + | /d1\ /d2\ | + | x*gamma|--|*gamma|--| | + \ \2 / \2 / / + + References + ========== + + [1] http://en.wikipedia.org/wiki/F-distribution + [2] http://mathworld.wolfram.com/F-Distribution.html + """ + + return FDistributionPSpace(name, d1, d2).value + +#------------------------------------------------------------------------------- +# Fisher Z distribution -------------------------------------------------------- + +class FisherZPSpace(SingleContinuousPSpace): + def __new__(cls, name, d1, d2): + d1 = sympify(d1) + d2 = sympify(d2) + x = Symbol(name) + pdf = (2*d1**(d1/2)*d2**(d2/2) / beta_fn(d1/2, d2/2) * + exp(d1*x) / (d1*exp(2*x)+d2)**((d1+d2)/2)) + obj = SingleContinuousPSpace.__new__(cls, x, pdf) + return obj + +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}} + + Parameters + ========== + + d1 : `d1` > 0, degree of freedom + d2 : `d2` > 0, degree of freedom + + Returns + ======= + + A RandomSymbol. + + Examples + ======== + + >>> from sympy.stats import FisherZ, density + >>> from sympy import Symbol, simplify, pprint + + >>> d1 = Symbol("d1", positive=True) + >>> d2 = Symbol("d2", positive=True) + + >>> X = FisherZ("x", d1, d2) + + >>> D = density(X) + >>> pprint(D, use_unicode=False) + / d1 d2 \ + | d1 d2 - -- - -- | + | -- -- 2 2 | + | 2 2 / 2*x \ x*d1 /d1 d2\| + | 2*d1 *d2 *\d1*e + d2/ *e *gamma|-- + --|| + | \2 2 /| + Lambda|x, --------------------------------------------------------| + | /d1\ /d2\ | + | gamma|--|*gamma|--| | + \ \2 / \2 / / + + References + ========== + + [1] http://en.wikipedia.org/wiki/Fisher%27s_z-distribution + [2] http://mathworld.wolfram.com/Fishersz-Distribution.html + """ + + return FisherZPSpace(name, d1, d2).value + +#------------------------------------------------------------------------------- +# Frechet distribution --------------------------------------------------------- + +class FrechetPSpace(SingleContinuousPSpace): + def __new__(cls, name, a, s=0, m=0): + a = sympify(a) + s = sympify(s) + m = sympify(m) + x = Symbol(name) + pdf = a/s * ((x-m)/s)**(-1-a) * exp(-((x-m)/s)-a) + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set = Interval(m, oo)) + return obj + +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 :math:`x \geq m`. + + Parameters + ========== + + a : Real number, :math:`a \in \left(0, \infty\right)` the shape + s : Real number, :math:`s \in \left(0, \infty\right)` the scale + m : Real number, :math:`m \in \left(-\infty, \infty\right)` the minimum + + Returns + ======= + + A RandomSymbol. + + Examples + ======== + + >>> from sympy.stats import Frechet, density, E, std + >>> from sympy import Symbol, simplify + + >>> a = Symbol("a", positive=True) + >>> s = Symbol("s", positive=True) + >>> m = Symbol("m", real=True) + + >>> X = Frechet("x", a, s, m) + + >>> density(X) + Lambda(_x, a*((_x - m)/s)**(-a - 1)*exp(-a - (_x - m)/s)/s) + + References + ========== + + [1] http://en.wikipedia.org/wiki/Fr%C3%A9chet_distribution + """ + + return FrechetPSpace(name, a, s, m).value + #------------------------------------------------------------------------------- # Gamma distribution ----------------------------------------------------------- @@ -730,6 +1137,141 @@ def Gamma(name, k, theta): return GammaPSpace(name, k, theta).value +#------------------------------------------------------------------------------- +# Inverse Gamma distribution --------------------------------------------------- + +class GammaInversePSpace(SingleContinuousPSpace): + def __new__(cls, name, a, b): + a = sympify(a) + b = sympify(b) + _value_check(a > 0, "alpha must be positive") + _value_check(b > 0, "beta must be positive") + x = Symbol(name) + pdf = b**a/gamma(a) * x**(-a-1) * exp(-b/x) + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(0, oo)) + obj.a = a + obj.b = b + return obj + +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 :math:`x > 0`. + + Parameters + ========== + + a : Real number, `a` > 0 a shape + b : Real number, `b` > 0 a scale + + Returns + ======= + + A RandomSymbol. + + Examples + ======== + + >>> from sympy.stats import GammaInverse, density, cdf, E, variance + >>> from sympy import Symbol, pprint + + >>> a = Symbol("a", positive=True) + >>> b = Symbol("b", positive=True) + + >>> X = GammaInverse("x", a, b) + + >>> D = density(X) + >>> pprint(D, use_unicode=False) + / -b\ + | --| + | -a - 1 a x | + | x *b *e | + Lambda|x, --------------| + \ gamma(a) / + + References + ========== + + [1] http://en.wikipedia.org/wiki/Inverse-gamma_distribution + """ + + return GammaInversePSpace(name, a, b).value + +#------------------------------------------------------------------------------- +# Kumaraswamy distribution ----------------------------------------------------- + +class KumaraswamyPSpace(SingleContinuousPSpace): + def __new__(cls, name, a, b): + a, b = sympify(a), sympify(b) + + _value_check(a > 0, "a must be positive") + _value_check(b > 0, "b must be positive") + + x = Symbol(name) + pdf = a * b * x**(a-1) * (1-x**a)**(b-1) + + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(0, 1)) + obj.a = a + obj.b = b + return obj + +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 :math:`x \in [0,1]`. + + Parameters + ========== + + a : Real number, `a` > 0 a shape + b : Real number, `b` > 0 a shape + + Returns + ======= + + A RandomSymbol. + + Examples + ======== + + >>> from sympy.stats import Kumaraswamy, density, E, variance + >>> from sympy import Symbol, simplify, pprint + + >>> a = Symbol("a", positive=True) + >>> b = Symbol("b", positive=True) + + >>> X = Kumaraswamy("x", a, b) + + >>> D = density(X) + >>> pprint(D, use_unicode=False) + / b - 1\ + | a - 1 / a \ | + Lambda\x, x *a*b*\- x + 1/ / + + >>> simplify(E(X, meijerg=True)) + gamma(1 + 1/a)*gamma(b + 1)/gamma(b + 1 + 1/a) + + References + ========== + + [1] http://en.wikipedia.org/wiki/Kumaraswamy_distribution + """ + + return KumaraswamyPSpace(name, a, b).value + #------------------------------------------------------------------------------- # Laplace distribution --------------------------------------------------------- @@ -1220,6 +1762,148 @@ def Pareto(name, xm, alpha): return ParetoPSpace(name, xm, alpha).value +#------------------------------------------------------------------------------- +# QuadraticU distribution ------------------------------------------------------ + +class QuadraticUPSpace(SingleContinuousPSpace): + def __new__(cls, name, a, b): + a, b = sympify(a), sympify(b) + alpha = 12 / (b-a)**3 + beta = (a+b) / 2 + + x = Symbol(name) + pdf = Piecewise( + (alpha * (x-beta)**2, And(a<=x, x<=b)), + (S.Zero, True)) + + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(a, b)) + obj.a = a + obj.b = b + return obj + +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 :math:`x \in [a,b]`. + + Parameters + ========== + + a : Real number + b : Real number, :math:`a < b` + + Returns + ======= + + A RandomSymbol. + + Examples + ======== + + >>> from sympy.stats import QuadraticU, density, E, variance + >>> from sympy import Symbol, simplify, factor, pprint + + >>> a = Symbol("a", real=True) + >>> b = Symbol("b", real=True) + + >>> X = QuadraticU("x", a, b) + + >>> D = density(X) + >>> pprint(D, use_unicode=False) + / / 2 \ + | | / a b\ | + | |12*|x - - - -| | + | | \ 2 2/ | + Lambda|x, <--------------- for And(x <= b, a <= x)| + | | 3 | + | | (-a + b) | + | | | + \ \ 0 otherwise / + + References + ========== + + [1] http://en.wikipedia.org/wiki/U-quadratic_distribution + """ + + return QuadraticUPSpace(name, a, b).value + +#------------------------------------------------------------------------------- +# RaisedCosine distribution ---------------------------------------------------- + +class RaisedCosinePSpace(SingleContinuousPSpace): + def __new__(cls, name, mu, s): + mu, s = sympify(mu), sympify(s) + + _value_check(s > 0, "s must be positive") + + x = Symbol(name) + pdf = Piecewise( + ((1+cos(pi*(x-mu)/s)) / (2*s), And(mu-s<=x, x<=mu+s)), + (S.Zero, True)) + + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(mu-s, mu+s)) + obj.mu = mu + obj.s = s + return obj + +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 :math:`x \in [\mu-s,\mu+s]`. + + Parameters + ========== + + mu : Real number + s : Real number, `s` > 0 + + Returns + ======= + + A RandomSymbol. + + Examples + ======== + + >>> from sympy.stats import RaisedCosine, density, E, variance + >>> from sympy import Symbol, simplify, pprint + + >>> mu = Symbol("mu", real=True) + >>> s = Symbol("s", positive=True) + + >>> X = RaisedCosine("x", mu, s) + + >>> D = density(X) + >>> pprint(D, use_unicode=False) + / / /pi*(x - mu)\ \ + | |cos|-----------| + 1 | + | | \ s / | + Lambda|x, <-------------------- for And(x <= mu + s, mu - s <= x)| + | | 2*s | + | | | + \ \ 0 otherwise / + + References + ========== + + [1] http://en.wikipedia.org/wiki/Raised_cosine_distribution + """ + + return RaisedCosinePSpace(name, mu, s).value + #------------------------------------------------------------------------------- # Rayleigh distribution -------------------------------------------------------- @@ -1593,6 +2277,72 @@ def UniformSum(name, n): return UniformSumPSpace(name, n).value +#------------------------------------------------------------------------------- +# VonMises distribution -------------------------------------------------------- + +class VonMisesPSpace(SingleContinuousPSpace): + def __new__(cls, name, mu, k): + mu, k = sympify(mu), sympify(k) + + _value_check(k > 0, "k must be positive") + + x = Symbol(name) + pdf = exp(k*cos(x-mu)) / (2*pi*besseli(0, k)) + + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(0, 2*pi)) + obj.mu = mu + obj.k = k + return obj + +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 :math:`x \in [0,2\pi]`. + + Parameters + ========== + + mu : Real number, measure of location + k : Real number, measure of concentration + + Returns + ======= + + A RandomSymbol. + + Examples + ======== + + >>> from sympy.stats import VonMises, density, E, variance + >>> from sympy import Symbol, simplify, pprint + + >>> mu = Symbol("mu") + >>> k = Symbol("k", positive=True) + + >>> X = VonMises("x", mu, k) + + >>> D = density(X) + >>> pprint(D, use_unicode=False) + / k*cos(x - mu) \ + | e | + Lambda|x, ------------------| + \ 2*pi*besseli(0, k)/ + + References + ========== + + [1] http://en.wikipedia.org/wiki/Von_Mises_distribution + [2] http://mathworld.wolfram.com/vonMisesDistribution.html + """ + + return VonMisesPSpace(name, mu, k).value + #------------------------------------------------------------------------------- # Weibull distribution --------------------------------------------------------- diff --git a/sympy/stats/tests/test_continuous_rv.py b/sympy/stats/tests/test_continuous_rv.py index 18d5e24ecf3b..98b46115e7f8 100644 --- a/sympy/stats/tests/test_continuous_rv.py +++ b/sympy/stats/tests/test_continuous_rv.py @@ -1,14 +1,18 @@ from sympy.stats import (P, E, where, density, variance, covariance, skewness, given, pspace, cdf, ContinuousRV, sample, - Arcsin, Benini, Beta, BetaPrime, Cauchy, Chi, Dagum, - Exponential, Gamma, Laplace, Logistic, LogNormal, - Maxwell, Nakagami, Normal, Pareto, Rayleigh, - StudentT, Triangular, Uniform, UniformSum, Weibull, - WignerSemicircle) + Arcsin, Benini, Beta, BetaPrime, Cauchy, Chi, ChiSquared, + ChiNoncentral, 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 sympy import (Symbol, Dummy, Abs, exp, S, N, pi, simplify, Interval, erf, Eq, log, lowergamma, Sum, symbols, sqrt, And, gamma, beta, - Piecewise, Integral, sin, Lambda, factorial, binomial, - floor) + Eq, log, lowergamma, Sum, symbols, sqrt, And, gamma, beta, + Piecewise, Integral, sin, cos, besseli, Lambda, factorial, + binomial, floor) + from sympy.utilities.pytest import raises, XFAIL oo = S.Infinity @@ -199,6 +203,20 @@ def test_chi(): assert density(X) == (Lambda(_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) == (Lambda(_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) == (Lambda(_x, + 2**(-k/2)*_x**(k/2 - 1)*exp(-_x/2)/gamma(k/2))) def test_dagum(): p = Symbol("p", positive=True) @@ -209,6 +227,12 @@ def test_dagum(): assert density(X) == Lambda(_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) == Lambda(_x, _x**(k - 1)*l**k*exp(-_x*l)/gamma(k)) def test_exponential(): rate = Symbol('lambda', positive=True, real=True, bounded=True) @@ -223,6 +247,29 @@ def test_exponential(): assert where(X <= 1).set == Interval(0, 1) +def test_f_distribution(): + d1 = Symbol("d1", positive=True) + d2 = Symbol("d2", positive=True) + + X = FDistribution("x", d1, d2) + assert density(X) == (Lambda(_x, d2**(d2/2)*sqrt((_x*d1)**d1*(_x*d1 + d2)**(-d1 - d2))* + gamma(d1/2 + d2/2)/(_x*gamma(d1/2)*gamma(d2/2)))) + +def test_fisher_z(): + d1 = Symbol("d1", positive=True) + d2 = Symbol("d2", positive=True) + + X = FisherZ("x", d1, d2) + assert density(X) == (Lambda(_x, 2*d1**(d1/2)*d2**(d2/2)*(d1*exp(2*_x) + d2)**(-d1/2 - d2/2)* + exp(_x*d1)*gamma(d1/2 + d2/2)/(gamma(d1/2)*gamma(d2/2)))) + +def test_frechet(): + a = Symbol("a", positive=True) + s = Symbol("s", positive=True) + m = Symbol("m", real=True) + + X = Frechet("x", a, s, m) + assert density(X) == Lambda(_x, a*((_x - m)/s)**(-a - 1)*exp(-a - (_x - m)/s)/s) def test_gamma(): k = Symbol("k", positive=True) @@ -244,6 +291,19 @@ def test_gamma(): # 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) == Lambda(_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) == Lambda(_x, _x**(a - 1)*a*b*(-_x**a + 1)**(b - 1)) def test_laplace(): mu = Symbol("mu") @@ -252,7 +312,6 @@ def test_laplace(): X = Laplace('x', mu, b) assert density(X) == Lambda(_x, exp(-Abs(_x - mu)/b)/(2*b)) - def test_logistic(): mu = Symbol("mu", real=True) s = Symbol("s", positive=True) @@ -261,7 +320,6 @@ def test_logistic(): assert density(X) == Lambda(_x, exp((-_x + mu)/s)/(s*(exp((-_x + mu)/s) + 1)**2)) - def test_lognormal(): mean = Symbol('mu', real=True, bounded=True) std = Symbol('sigma', positive=True, real=True, bounded=True) @@ -288,7 +346,6 @@ def test_lognormal(): X = LogNormal('x', 0, 1) # Mean 0, standard deviation 1 assert density(X) == Lambda(_x, sqrt(2)*exp(-log(_x)**2/2)/(2*_x*sqrt(pi))) - def test_maxwell(): a = Symbol("a", positive=True) @@ -338,6 +395,15 @@ def test_pareto_numeric(): # Skewness tests too slow. Try shortcutting function? +def test_raised_cosine(): + mu = Symbol("mu", real=True) + s = Symbol("s", positive=True) + + X = RaisedCosine("x", mu, s) + assert density(X) == (Lambda(_x, Piecewise(((cos(pi*(_x - mu)/s) + 1)/(2*s), + And(_x <= mu + s, mu - s <= _x)), (0, True)))) + + def test_rayleigh(): sigma = Symbol("sigma", positive=True) @@ -370,6 +436,15 @@ def test_triangular(): (0, True))) +@XFAIL +def test_quadratic_u(): + a = Symbol("a", real=True) + b = Symbol("b", real=True) + + X = QuadraticU("x", a, b) + assert density(X) == (Lambda(_x, Piecewise((12*(_x - a/2 - b/2)**2/(-a + b)**3, + And(_x <= b, a <= _x)), (0, True)))) + def test_uniform(): l = Symbol('l', real=True, bounded=True) w = Symbol('w', positive=True, bounded=True) @@ -396,6 +471,14 @@ def test_uniformsum(): *binomial(n, _k), (_k, 0, floor(_x)))/factorial(n - 1))) +def test_von_mises(): + mu = Symbol("mu") + k = Symbol("k", positive=True) + + X = VonMises("x", mu, k) + assert density(X) == Lambda(_x, exp(k*cos(_x - mu))/(2*pi*besseli(0, k))) + + def test_weibull(): a, b = symbols('a b', positive=True) X = Weibull('x', a, b)