Vehicle-cpp/hkpc/software/sympy/stats/tests/test_symbolic_probability.py

from sympy.concrete.summations import Sum
from sympy.core.mul import Mul
from sympy.core.numbers import (oo, pi)
from sympy.core.relational import Eq
from sympy.core.symbol import (Dummy, symbols)
from sympy.functions.elementary.exponential import exp
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.functions.elementary.trigonometric import sin
from sympy.integrals.integrals import Integral
from sympy.core.expr import unchanged
from sympy.stats import (Normal, Poisson, variance, Covariance, Variance,
                         Probability, Expectation, Moment, CentralMoment)
from sympy.stats.rv import probability, expectation


def test_literal_probability():
    X = Normal('X', 2, 3)
    Y = Normal('Y', 3, 4)
    Z = Poisson('Z', 4)
    W = Poisson('W', 3)
    x = symbols('x', real=True)
    y, w, z = symbols('y, w, z')

    assert Probability(X > 0).evaluate_integral() == probability(X > 0)
    assert Probability(X > x).evaluate_integral() == probability(X > x)
    assert Probability(X > 0).rewrite(Integral).doit() == probability(X > 0)
    assert Probability(X > x).rewrite(Integral).doit() == probability(X > x)

    assert Expectation(X).evaluate_integral() == expectation(X)
    assert Expectation(X).rewrite(Integral).doit() == expectation(X)
    assert Expectation(X**2).evaluate_integral() == expectation(X**2)
    assert Expectation(x*X).args == (x*X,)
    assert Expectation(x*X).expand() == x*Expectation(X)
    assert Expectation(2*X + 3*Y + z*X*Y).expand() == 2*Expectation(X) + 3*Expectation(Y) + z*Expectation(X*Y)
    assert Expectation(2*X + 3*Y + z*X*Y).args == (2*X + 3*Y + z*X*Y,)
    assert Expectation(sin(X)) == Expectation(sin(X)).expand()
    assert Expectation(2*x*sin(X)*Y + y*X**2 + z*X*Y).expand() == 2*x*Expectation(sin(X)*Y) \
                            + y*Expectation(X**2) + z*Expectation(X*Y)
    assert Expectation(X + Y).expand() ==  Expectation(X) + Expectation(Y)
    assert Expectation((X + Y)*(X - Y)).expand() == Expectation(X**2) - Expectation(Y**2)
    assert Expectation((X + Y)*(X - Y)).expand().doit() == -12
    assert Expectation(X + Y, evaluate=True).doit() == 5
    assert Expectation(X + Expectation(Y)).doit() == 5
    assert Expectation(X + Expectation(Y)).doit(deep=False) == 2 + Expectation(Expectation(Y))
    assert Expectation(X + Expectation(Y + Expectation(2*X))).doit(deep=False) == 2 \
                                + Expectation(Expectation(Y + Expectation(2*X)))
    assert Expectation(X + Expectation(Y + Expectation(2*X))).doit() == 9
    assert Expectation(Expectation(2*X)).doit() == 4
    assert Expectation(Expectation(2*X)).doit(deep=False) == Expectation(2*X)
    assert Expectation(4*Expectation(2*X)).doit(deep=False) == 4*Expectation(2*X)
    assert Expectation((X + Y)**3).expand() == 3*Expectation(X*Y**2) +\
                3*Expectation(X**2*Y) + Expectation(X**3) + Expectation(Y**3)
    assert Expectation((X - Y)**3).expand() == 3*Expectation(X*Y**2) -\
                3*Expectation(X**2*Y) + Expectation(X**3) - Expectation(Y**3)
    assert Expectation((X - Y)**2).expand() == -2*Expectation(X*Y) +\
                Expectation(X**2) + Expectation(Y**2)

    assert Variance(w).args == (w,)
    assert Variance(w).expand() == 0
    assert Variance(X).evaluate_integral() == Variance(X).rewrite(Integral).doit() == variance(X)
    assert Variance(X + z).args == (X + z,)
    assert Variance(X + z).expand() == Variance(X)
    assert Variance(X*Y).args == (Mul(X, Y),)
    assert type(Variance(X*Y)) == Variance
    assert Variance(z*X).expand() == z**2*Variance(X)
    assert Variance(X + Y).expand() == Variance(X) + Variance(Y) + 2*Covariance(X, Y)
    assert Variance(X + Y + Z + W).expand() == (Variance(X) + Variance(Y) + Variance(Z) + Variance(W) +
                                       2 * Covariance(X, Y) + 2 * Covariance(X, Z) + 2 * Covariance(X, W) +
                                       2 * Covariance(Y, Z) + 2 * Covariance(Y, W) + 2 * Covariance(W, Z))
    assert Variance(X**2).evaluate_integral() == variance(X**2)
    assert unchanged(Variance, X**2)
    assert Variance(x*X**2).expand() == x**2*Variance(X**2)
    assert Variance(sin(X)).args == (sin(X),)
    assert Variance(sin(X)).expand() == Variance(sin(X))
    assert Variance(x*sin(X)).expand() == x**2*Variance(sin(X))

    assert Covariance(w, z).args == (w, z)
    assert Covariance(w, z).expand() == 0
    assert Covariance(X, w).expand() == 0
    assert Covariance(w, X).expand() == 0
    assert Covariance(X, Y).args == (X, Y)
    assert type(Covariance(X, Y)) == Covariance
    assert Covariance(z*X + 3, Y).expand() == z*Covariance(X, Y)
    assert Covariance(X, X).args == (X, X)
    assert Covariance(X, X).expand() == Variance(X)
    assert Covariance(z*X + 3, w*Y + 4).expand() == w*z*Covariance(X,Y)
    assert Covariance(X, Y) == Covariance(Y, X)
    assert Covariance(X + Y, Z + W).expand() == Covariance(W, X) + Covariance(W, Y) + Covariance(X, Z) + Covariance(Y, Z)
    assert Covariance(x*X + y*Y, z*Z + w*W).expand() == (x*w*Covariance(W, X) + w*y*Covariance(W, Y) +
                                                x*z*Covariance(X, Z) + y*z*Covariance(Y, Z))
    assert Covariance(x*X**2 + y*sin(Y), z*Y*Z**2 + w*W).expand() == (w*x*Covariance(W, X**2) + w*y*Covariance(sin(Y), W) +
                                                        x*z*Covariance(Y*Z**2, X**2) + y*z*Covariance(Y*Z**2, sin(Y)))
    assert Covariance(X, X**2).expand() == Covariance(X, X**2)
    assert Covariance(X, sin(X)).expand() == Covariance(sin(X), X)
    assert Covariance(X**2, sin(X)*Y).expand() == Covariance(sin(X)*Y, X**2)
    assert Covariance(w, X).evaluate_integral() == 0


def test_probability_rewrite():
    X = Normal('X', 2, 3)
    Y = Normal('Y', 3, 4)
    Z = Poisson('Z', 4)
    W = Poisson('W', 3)
    x, y, w, z = symbols('x, y, w, z')

    assert Variance(w).rewrite(Expectation) == 0
    assert Variance(X).rewrite(Expectation) == Expectation(X ** 2) - Expectation(X) ** 2
    assert Variance(X, condition=Y).rewrite(Expectation) == Expectation(X ** 2, Y) - Expectation(X, Y) ** 2
    assert Variance(X, Y) != Expectation(X**2) - Expectation(X)**2
    assert Variance(X + z).rewrite(Expectation) == Expectation((X + z) ** 2) - Expectation(X + z) ** 2
    assert Variance(X * Y).rewrite(Expectation) == Expectation(X ** 2 * Y ** 2) - Expectation(X * Y) ** 2

    assert Covariance(w, X).rewrite(Expectation) == -w*Expectation(X) + Expectation(w*X)
    assert Covariance(X, Y).rewrite(Expectation) == Expectation(X*Y) - Expectation(X)*Expectation(Y)
    assert Covariance(X, Y, condition=W).rewrite(Expectation) == Expectation(X * Y, W) - Expectation(X, W) * Expectation(Y, W)

    w, x, z = symbols("W, x, z")
    px = Probability(Eq(X, x))
    pz = Probability(Eq(Z, z))

    assert Expectation(X).rewrite(Probability) == Integral(x*px, (x, -oo, oo))
    assert Expectation(Z).rewrite(Probability) == Sum(z*pz, (z, 0, oo))
    assert Variance(X).rewrite(Probability) == Integral(x**2*px, (x, -oo, oo)) - Integral(x*px, (x, -oo, oo))**2
    assert Variance(Z).rewrite(Probability) == Sum(z**2*pz, (z, 0, oo)) - Sum(z*pz, (z, 0, oo))**2
    assert Covariance(w, X).rewrite(Probability) == \
           -w*Integral(x*Probability(Eq(X, x)), (x, -oo, oo)) + Integral(w*x*Probability(Eq(X, x)), (x, -oo, oo))

    # To test rewrite as sum function
    assert Variance(X).rewrite(Sum) == Variance(X).rewrite(Integral)
    assert Expectation(X).rewrite(Sum) == Expectation(X).rewrite(Integral)

    assert Covariance(w, X).rewrite(Sum) == 0

    assert Covariance(w, X).rewrite(Integral) == 0

    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', 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**5


def test_symbolic_CentralMoment():
    mu = symbols('mu', real=True)
    sigma = symbols('sigma', 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**5