Merge pull request #19724 from Smit-create/gsoc_symbolic_moment

[GSoC] Added symbolic classes for Moment and Centralmoment

sympy/core/tests/test_args.py

@@ -1204,6 +1204,20 @@ def test_sympy__stats__symbolic_probability__Variance():
     assert _test_args(Variance(X))
 
 
+def test_sympy__stats__symbolic_probability__Moment():
+    from sympy.stats.symbolic_probability import Moment
+    from sympy.stats import Normal
+    X = Normal('X', 0, 1)
+    assert _test_args(Moment(X, 3, 2, X > 3))
+
+
+def test_sympy__stats__symbolic_probability__CentralMoment():
+    from sympy.stats.symbolic_probability import CentralMoment
+    from sympy.stats import Normal
+    X = Normal('X', 0, 1)
+    assert _test_args(CentralMoment(X, 2, X > 1))
+
+
 def test_sympy__stats__frv_types__DiscreteUniformDistribution():
     from sympy.stats.frv_types import DiscreteUniformDistribution
     from sympy.core.containers import Tuple

sympy/stats/__init__.py

@@ -142,7 +142,7 @@
     'joint_eigen_distribution', 'JointEigenDistribution',
     'level_spacing_distribution',
 
-    'Probability', 'Expectation', 'Variance', 'Covariance',
+    'Probability', 'Expectation', 'Variance', 'Covariance', 'Moment', 'CentralMoment',
 
     'ExpectationMatrix', 'VarianceMatrix', 'CrossCovarianceMatrix'
 
@@ -189,7 +189,7 @@
         JointEigenDistribution, level_spacing_distribution)
 
 from .symbolic_probability import (Probability, Expectation, Variance,
-        Covariance)
+        Covariance, Moment, CentralMoment)
 
 from .symbolic_multivariate_probability import (ExpectationMatrix, VarianceMatrix,
         CrossCovarianceMatrix)

sympy/stats/rv_interface.py

@@ -1,6 +1,7 @@
 from __future__ import print_function, division
 from sympy.sets import FiniteSet
-from sympy import sqrt, log, exp, FallingFactorial, Rational, Eq, Dummy, piecewise_fold, solveset
+from sympy import (sqrt, log, exp, FallingFactorial, Rational, Eq, Dummy,
+                piecewise_fold, solveset, Integral)
 from .rv import (probability, expectation, density, where, given, pspace, cdf, PSpace,
                  characteristic_function, sample, sample_iter, random_symbols, independent, dependent,
                  sampling_density, moment_generating_function, quantile, is_random,
@@ -33,7 +34,10 @@ def moment(X, n, c=0, condition=None, **kwargs):
     >>> moment(X, 1) == E(X)
     True
     """
-    return expectation((X - c)**n, condition, **kwargs)
+    from sympy.stats.symbolic_probability import Moment
+    if kwargs.pop('evaluate', True):
+        return Moment(X, n, c, condition).doit()
+    return Moment(X, n, c, condition).rewrite(Integral)
 
 
 def variance(X, condition=None, **kwargs):
@@ -211,8 +215,10 @@ def cmoment(X, n, condition=None, **kwargs):
     >>> cmoment(X, 2) == variance(X)
     True
     """
-    mu = expectation(X, condition, **kwargs)
-    return moment(X, n, mu, condition, **kwargs)
+    from sympy.stats.symbolic_probability import CentralMoment
+    if kwargs.pop('evaluate', True):
+        return CentralMoment(X, n, condition).doit()
+    return CentralMoment(X, n, condition).rewrite(Integral)
 
 
 def smoment(X, n, condition=None, **kwargs):

sympy/stats/symbolic_probability.py

@@ -561,3 +561,111 @@ def _eval_rewrite_as_Integral(self, arg1, arg2, condition=None, **kwargs):
 
     def evaluate_integral(self):
         return self.rewrite(Integral).doit()
+
+
+class Moment(Expr):
+    """
+    Symbolic class for Moment
+
+    Examples
+    ========
+
+    >>> from sympy import Symbol, Integral
+    >>> from sympy.stats import Normal, Expectation, Probability, Moment
+    >>> mu = Symbol('mu', real=True)
+    >>> sigma = Symbol('sigma', real=True, positive=True)
+    >>> X = Normal('X', mu, sigma)
+    >>> M = Moment(X, 3, 1)
+
+    To evaluate the result of Moment use `doit`:
+    >>> M.doit()
+    mu**3 - 3*mu**2 + 3*mu*sigma**2 + 3*mu - 3*sigma**2 - 1
+
+    Rewrite the Moment expression in terms of Expectation:
+    >>> M.rewrite(Expectation)
+    Expectation((X - 1)**3)
+
+    Rewrite the Moment expression in terms of Probability:
+    >>> M.rewrite(Probability)
+    Integral((x - 1)**3*Probability(Eq(X, x)), (x, -oo, oo))
+
+    Rewrite the Moment expression in terms of Integral:
+    >>> M.rewrite(Integral)
+    Integral(sqrt(2)*(X - 1)**3*exp(-(X - mu)**2/(2*sigma**2))/(2*sqrt(pi)*sigma), (X, -oo, oo))
+
+    """
+    def __new__(cls, X, n, c=0, condition=None, **kwargs):
+        X = _sympify(X)
+        n = _sympify(n)
+        c = _sympify(c)
+        if condition is not None:
+            condition = _sympify(condition)
+        return Expr.__new__(cls, X, n, c, condition)
+
+    def doit(self, **hints):
+        if not is_random(self.args[0]):
+            return self.args[0]
+        return self.rewrite(Expectation).doit(**hints)
+
+    def _eval_rewrite_as_Expectation(self, X, n, c=0, condition=None, **kwargs):
+        return Expectation((X - c)**n, condition)
+
+    def _eval_rewrite_as_Probability(self, X, n, c=0, condition=None, **kwargs):
+        return self.rewrite(Expectation).rewrite(Probability)
+
+    def _eval_rewrite_as_Integral(self, X, n, c=0, condition=None, **kwargs):
+        return self.rewrite(Expectation).rewrite(Integral)
+
+
+class CentralMoment(Expr):
+    """
+    Symbolic class Central Moment
+
+    Examples
+    ========
+
+    >>> from sympy import Symbol, Integral
+    >>> from sympy.stats import Normal, Expectation, Probability, CentralMoment
+    >>> mu = Symbol('mu', real=True)
+    >>> sigma = Symbol('sigma', real=True, positive=True)
+    >>> X = Normal('X', mu, sigma)
+    >>> CM = CentralMoment(X, 4)
+
+    To evaluate the result of CentralMoment use `doit`:
+    >>> CM.doit().simplify()
+    3*sigma**4
+
+    Rewrite the CentralMoment expression in terms of Expectation:
+    >>> CM.rewrite(Expectation)
+    Expectation((X - Expectation(X))**4)
+
+    Rewrite the CentralMoment expression in terms of Probability:
+    >>> CM.rewrite(Probability)
+    Integral((x - Integral(x*Probability(True), (x, -oo, oo)))**4*Probability(Eq(X, x)), (x, -oo, oo))
+
+    Rewrite the CentralMoment expression in terms of Integral:
+    >>> CM.rewrite(Integral)
+    Integral(sqrt(2)*(X - Integral(sqrt(2)*X*exp(-(X - mu)**2/(2*sigma**2))/(2*sqrt(pi)*sigma), (X, -oo, oo)))**4*exp(-(X - mu)**2/(2*sigma**2))/(2*sqrt(pi)*sigma), (X, -oo, oo))
+
+    """
+    def __new__(cls, X, n, condition=None, **kwargs):
+        X = _sympify(X)
+        n = _sympify(n)
+        if condition is not None:
+            condition = _sympify(condition)
+        return Expr.__new__(cls, X, n, condition)
+
+    def doit(self, **hints):
+        if not is_random(self.args[0]):
+            return self.args[0]
+        return self.rewrite(Expectation).doit(**hints)
+
+    def _eval_rewrite_as_Expectation(self, X, n, condition=None, **kwargs):
+        mu = Expectation(X, condition, **kwargs)
+        return Moment(X, n, mu, condition, **kwargs).rewrite(Expectation)
+
+    def _eval_rewrite_as_Probability(self, X, n, condition=None, **kwargs):
+        return self.rewrite(Expectation).rewrite(Probability)
+
+    def _eval_rewrite_as_Integral(self, X, n, condition=None, **kwargs):
+        return self.rewrite(Expectation).rewrite(Integral)

sympy/stats/tests/test_symbolic_probability.py

@@ -1,7 +1,7 @@
-from sympy import symbols, Mul, sin, Integral, oo, Eq, Sum
+from sympy import symbols, Mul, sin, Integral, oo, Eq, Sum, sqrt, exp, pi, Dummy
 from sympy.core.expr import unchanged
 from sympy.stats import (Normal, Poisson, variance, Covariance, Variance,
-                         Probability, Expectation)
+                         Probability, Expectation, Moment, CentralMoment)
 from sympy.stats.rv import probability, expectation
 
 
@@ -127,3 +127,41 @@ def test_probability_rewrite():
 
     assert Variance(X, condition=Y).rewrite(Probability) == Integral(x**2*Probability(Eq(X, x), Y), (x, -oo, oo)) - \
                                                             Integral(x*Probability(Eq(X, x), Y), (x, -oo, oo))**2
+
+
+def test_symbolic_Moment():
+    mu = symbols('mu', real=True)
+    sigma = symbols('sigma', real=True, positive=True)
+    x = symbols('x')
+    X = Normal('X', mu, sigma)
+    M = Moment(X, 4, 2)
+    assert M.rewrite(Expectation) == Expectation((X - 2)**4)
+    assert M.rewrite(Probability) == Integral((x - 2)**4*Probability(Eq(X, x)),
+                                    (x, -oo, oo))
+    k = Dummy('k')
+    expri = Integral(sqrt(2)*(k - 2)**4*exp(-(k - \
+                mu)**2/(2*sigma**2))/(2*sqrt(pi)*sigma), (k, -oo, oo))
+    assert M.rewrite(Integral).dummy_eq(expri)
+    assert M.doit() == (mu**4 - 8*mu**3 + 6*mu**2*sigma**2 + \
+                24*mu**2 - 24*mu*sigma**2 - 32*mu + 3*sigma**4 + 24*sigma**2 + 16)
+    M = Moment(2, 5)
+    assert M.doit() == 2
+
+
+def test_symbolic_CentralMoment():
+    mu = symbols('mu', real=True)
+    sigma = symbols('sigma', real=True, positive=True)
+    x = symbols('x')
+    X = Normal('X', mu, sigma)
+    CM = CentralMoment(X, 6)
+    assert CM.rewrite(Expectation) == Expectation((X - Expectation(X))**6)
+    assert CM.rewrite(Probability) == Integral((x - Integral(x*Probability(True),
+                    (x, -oo, oo)))**6*Probability(Eq(X, x)), (x, -oo, oo))
+    k = Dummy('k')
+    expri = Integral(sqrt(2)*(k - Integral(sqrt(2)*k*exp(-(k - \
+        mu)**2/(2*sigma**2))/(2*sqrt(pi)*sigma), (k, -oo, oo)))**6*exp(-(k - \
+        mu)**2/(2*sigma**2))/(2*sqrt(pi)*sigma), (k, -oo, oo))
+    assert CM.rewrite(Integral).dummy_eq(expri)
+    assert CM.doit().simplify() == 15*sigma**6
+    CM = Moment(5, 5)
+    assert CM.doit() == 5