Repository: sympy/sympy

Issue:
[GSoC Discussion] Compound Distributions
The current compound distributions in the `stats` don't provide expected results and are not tested as:
```python
>>> N = Normal('N',1,3)
>>> M = Normal('M',N,2)
>>> E(M)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/smit/Smitlunagariya/sympy/sympy/stats/rv.py", line 763, in expectation
    result = pspace(expr).compute_expectation(expr, evaluate=evaluate, **kwargs)
  File "/home/smit/Smitlunagariya/sympy/sympy/stats/joint_rv.py", line 126, in compute_expectation
    cannot be calculated yet.'''))
NotImplementedError:
Expectations of expression with unindexed joint random symbols cannot
be calculated yet
```
Therefore, I plan to add compound distributions from scratch.
The API, for using them will be same of other random variables for example:
```python
>>> M = Normal('M', 1, 3)
>>> N = Normal('N', M, 2)
>>> density(N)(x)
...
>>> E(N)
...
>>> P(N>2)
...
```
The mathematical implementation according to the definition of compound distributions [here](https://en.wikipedia.org/wiki/Compound_probability_distribution#Definition) can become unreliable as of now because we are not always sure, that the `integrate` will evaluate the results correctly.
For example, the  implementation of Normal distribution compounded with another Normal Distribution:
```python
>>> from sympy import integrate, oo
>>> from sympy.stats import Normal, density
>>> from sympy.abc import x,y,z
>>> N = Normal('N', 1, 2)
>>> M = Normal('M', x, 1) # We will compound M with x and x will follow distribution of N
>>> expr1 = density(M)(y)
>>> expr1
sqrt(2)*exp(-(-x + y)**2/2)/(2*sqrt(pi))    # dist. of M parametrized by x
>>> expr2=density(N)(x) 
>>> expr2
sqrt(2)*exp(-(x - 1)**2/8)/(4*sqrt(pi))      # dist. of x
>>> expr3 = integrate(expr1*expr2, (x,-oo,oo)).simplify()   # the dist. of M free from x.
>>> expr3
sqrt(10)*exp(-y**2/10 + y/5 - 1/10)/(10*sqrt(pi))      # the dist. of M with mean 1, correctly reached
>>> integrate(expr3, (y,-oo,oo)).simplify()
-erfc(sqrt(10)/10)/2 - erf(sqrt(10)/10)/2 + 3/2    #  It is correct but should evaluate to one itself.
```
Related issues #17007, #16855, #11135
ping @czgdp1807 @Upabjojr @jmig5776 

Discussion:
sympy-bot: 

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v160). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).

:x: There was an issue with the release notes. **Please do not close this pull request;** instead edit the description after reading the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).

* No release notes were found. Please edit the PR description and write the release notes under `<!-- BEGIN RELEASE NOTES -->`.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it. <details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- Your title above should be a short description of what
    was changed. Do not include the issue number in the title. -->

    #### References to other Issues or PRs
    <!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
    format, e.g. "Fixes #1234" (see
    https://tinyurl.com/auto-closing for more information). Also, please
    write a comment on that issue linking back to this pull request once it is
    open. -->


    #### Brief description of what is fixed or changed
    Currently, I have implemented the initial handling of Compound Distributions, I will show some examples of how this is working. Also, it is not completed yet and so requires code cleanings and touch-ups.

    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below. See
    https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
    on how to write release notes. The bot will check your release notes
    automatically to see if they are formatted correctly. -->

    <!-- BEGIN RELEASE NOTES -->

    <!-- END RELEASE NOTES -->

</details><p>

sympy-bot: ### 🟠

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v160). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at https://github.com/sympy/sympy-bot/issues/75.

The following commits **add new files**:
* 38423d55a026b0983e19e2c31d8f54474345d534:
  - `sympy/stats/compound_rv.py`

If these files were added/deleted on purpose, you can ignore this message.

czgdp1807: Please add tests otherwise it's difficult to figure out what's going on.
Smit-create: Some of the examples where this works perfectly:
1. Compounding a Gaussian distribution with mean distributed according to another Gaussian distribution yields (again) a Gaussian distribution.

```python
>>> from sympy.stats import *
>>> from sympy import symbols
>>> x,y,z=symbols('x y z', real=True)
>>> X = Normal('X',3,4)
>>> Y = Normal('Y',X, 5)
>>> E(Y)
3
>>> density(Y)(x).simplify()
sqrt(82)*exp(-x**2/82 + 3*x/41 - 9/82)/(82*sqrt(pi))
>>> variance(Y).simplify()== 4**2 +  5**2 # https://math.stackexchange.com/questions/1484451/density-of-a-normal-rv-whose-mean-is-drawn-from-a-normal-distribution-compound
True
```

2. Compounding a Poisson distribution with rate parameter distributed according to a gamma distribution yields a negative binomial distribution.

```python
>>> k,t = symbols('k t', positive=True, real=True)
>>> G = Gamma('G',k,t)
>>> D = Poisson('P', G)
>>> density(D)(x).simplify()
Piecewise((t**x*(t + 1)**(-k - x)*gamma(k + x)/(gamma(k)*gamma(x + 1)), k > -x), (t**(-k)*Integral(exp(-G - G/t)*G**(k + x - 1), (G, 0, oo))/(gamma(k)*gamma(x + 1)), True))
>>> E(D).simplify() # https://en.wikipedia.org/wiki/Negative_binomial_distribution
k*t
```

Some of the examples where this fails due to much complex integral:
```python
>>> R = Rayleigh('R', 4)
>>> Y = Normal('Y', 2, R)
>>> density(Y)(z).simplify()
/               1   z                                        
|               - - -                                        
|               2   4                                        
|              e                                           pi
|              ------                for 2*|arg(z - 2)| <= --
|                8                                         2 
|                                                            
|       oo                                                   
|        /                                                   
|       |                                                    
|       |    /          2    4\                              
<       |   -\16*(z - 2)  + R /                              
|       |   --------------------                             
|       |              2                                     
|  ___  |          32*R                                      
|\/ 2 * |  e                     dR                          
|       |                                                    
|      /                                                     
|      0                                                     
|----------------------------------         otherwise        
|                 ____                                       
|            32*\/ pi                                        
\                                                            
>>> E(Y, evaluate=False)
 oo                                                                           
  /                                                                           
 |                                                                            
 |  /  /    ____     /Y   1\     ____     /Y   1\\                            
 |  |Y*|- \/ pi *sinh|- - -| + \/ pi *cosh|- - -||                            
 |  |  \             \4   2/              \4   2//                        pi  
 |  |---------------------------------------------  for 2*|arg(Y - 2)| <= --  
 |  |                       ____                                          2   
 |  |                   8*\/ pi                                               
 |  |                                                                         
 |  |         oo                                                              
 |  |          /                                                              
 |  |         |                                                               
 |  |         |                 2                                             
 |  <         |         -(Y - 2)      2                                      dY
 |  |         |         ----------  -R                                        
 |  |         |               2     ----                                      
 |  |         |    ___     2*R       32                                       
 |  |         |  \/ 2 *e          *e                                          
 |  |      Y* |  ----------------------- dR                otherwise          
 |  |         |              ____                                             
 |  |         |         32*\/ pi                                              
 |  |         |                                                               
 |  |        /                                                                
 |  |        0                                                                
 |  \                                                                         
 |                                                                            
/                                                                             
```

Patch:
diff --git a/sympy/stats/compound_rv.py b/sympy/stats/compound_rv.py
new file mode 100644
index 0000000000..e0dba16c2e
--- /dev/null
+++ b/sympy/stats/compound_rv.py
@@ -0,0 +1,203 @@
+from sympy import Basic, integrate, Sum, Dummy, Lambda
+from sympy.stats.rv import (NamedArgsMixin, random_symbols, _symbol_converter,
+                        PSpace, RandomSymbol)
+from sympy.stats.crv import ContinuousDistribution, SingleContinuousPSpace
+from sympy.stats.drv import DiscreteDistribution, SingleDiscretePSpace
+from sympy.stats.frv import SingleFiniteDistribution, SingleFinitePSpace
+from sympy.stats.crv_types import ContinuousDistributionHandmade
+from sympy.stats.drv_types import DiscreteDistributionHandmade
+from sympy.stats.frv_types import FiniteDistributionHandmade
+
+
+class CompoundPSpace(PSpace):
+    """
+    A temporary Probability Space for the Compound Distribution. After
+    Marginalization, this returns the corresponding Probability Space of the
+    parent distribution.
+    """
+
+    def __new__(cls, s, distribution):
+        s = _symbol_converter(s)
+        if isinstance(distribution, ContinuousDistribution):
+            return SingleContinuousPSpace(s, distribution)
+        if isinstance(distribution, DiscreteDistribution):
+            return SingleDiscretePSpace(s, distribution)
+        if isinstance(distribution, SingleFiniteDistribution):
+            return SingleFinitePSpace(s, distribution)
+        if not isinstance(distribution, CompoundDistribution):
+            raise ValueError("%s should be an isinstance of "
+                        "CompoundDistribution"%(distribution))
+        return Basic.__new__(cls, s, distribution)
+
+    @property
+    def value(self):
+        return RandomSymbol(self.symbol, self)
+
+    @property
+    def symbol(self):
+        return self.args[0]
+
+    @property
+    def distribution(self):
+        return self.args[1]
+
+    @property
+    def pdf(self):
+        return self.distribution.pdf(self.symbol)
+
+    @property
+    def set(self):
+        return self.distribution.set
+
+    @property
+    def domain(self):
+        return self._get_newpspace().domain
+
+    def _get_newpspace(self):
+        x = Dummy('x')
+        parent_dist = self.distribution.args[0]
+        new_pspace = self._transform_pspace(self.symbol, parent_dist,
+                            Lambda(x, self.distribution.pdf(x)))
+        if new_pspace is not None:
+            return new_pspace
+        message = ("Compound Distribution for %s is not implemeted yet" % str(parent_dist))
+        raise NotImplementedError(message)
+
+    def _transform_pspace(self, sym, dist, pdf):
+        """
+        This function returns the new pspace of the distribution using handmade
+        Distributions and their corresponding pspace.
+        """
+        pdf = Lambda(sym, pdf(sym))
+        _set = dist.set
+        if isinstance(dist, ContinuousDistribution):
+            return SingleContinuousPSpace(sym, ContinuousDistributionHandmade(pdf, _set))
+        elif isinstance(dist, DiscreteDistribution):
+            return SingleDiscretePSpace(sym, DiscreteDistributionHandmade(pdf, _set))
+        elif isinstance(dist, SingleFiniteDistribution):
+            dens = dict((k, pdf(k)) for k in _set)
+            return SingleFinitePSpace(sym, FiniteDistributionHandmade(dens))
+
+    def compute_density(self, expr, **kwargs):
+        new_pspace = self._get_newpspace()
+        expr = expr.subs({self.value: new_pspace.value})
+        return new_pspace.compute_density(expr, **kwargs)
+
+    def compute_cdf(self, expr, **kwargs):
+        new_pspace = self._get_newpspace()
+        expr = expr.subs({self.value: new_pspace.value})
+        return new_pspace.compute_cdf(expr, **kwargs)
+
+    def compute_expectation(self, expr, rvs=None, evaluate=False, **kwargs):
+        new_pspace = self._get_newpspace()
+        expr = expr.subs({self.value: new_pspace.value})
+        if rvs:
+            rvs = rvs.subs({self.value: new_pspace.value})
+        if isinstance(new_pspace, SingleFinitePSpace):
+            return new_pspace.compute_expectation(expr, rvs, **kwargs)
+        return new_pspace.compute_expectation(expr, rvs, evaluate, **kwargs)
+
+    def probability(self, condition, **kwargs):
+        new_pspace = self._get_newpspace()
+        condition = condition.subs({self.value: new_pspace.value})
+        return new_pspace.probability(condition)
+
+    def conditional_space(self, condition, **kwargs):
+        new_pspace = self._get_newpspace()
+        condition = condition.subs({self.value: new_pspace.value})
+        return new_pspace.conditional_space(condition)
+
+
+class CompoundDistribution(Basic, NamedArgsMixin):
+    """
+    Class for Compound Distributions.
+
+    Parameters
+    ==========
+
+    dist : Distribution
+        Distribution must contain a random parameter
+
+    Examples
+    ========
+
+    >>> from sympy.stats.compound_rv import CompoundDistribution
+    >>> from sympy.stats.crv_types import NormalDistribution
+    >>> from sympy.stats import Normal
+    >>> from sympy.abc import x
+    >>> X = Normal('X', 2, 4)
+    >>> N = NormalDistribution(X, 4)
+    >>> C = CompoundDistribution(N)
+    >>> C.set
+    Interval(-oo, oo)
+    >>> C.pdf(x).simplify()
+    exp(-x**2/64 + x/16 - 1/16)/(8*sqrt(pi))
+
+    References
+    ==========
+
+    .. [1] https://en.wikipedia.org/wiki/Compound_probability_distribution
+
+    """
+
+    is_Finite = None
+    is_Continuous = None
+    is_Discrete = None
+
+    def __new__(cls, dist):
+        if isinstance(dist, ContinuousDistribution):
+            cls.is_Continuous = True
+        elif isinstance(dist, DiscreteDistribution):
+            cls.is_Discrete = True
+        elif isinstance(dist, SingleFiniteDistribution):
+            cls.is_Finite = True
+        else:
+            message = "Compound Distribution for %s is not implemeted yet" % str(dist)
+            raise NotImplementedError(message)
+        if not cls._compound_check(dist):
+            return dist
+        return Basic.__new__(cls, dist)
+
+    @property
+    def set(self):
+        return self.args[0].set
+
+    def pdf(self, x):
+        dist = self.args[0]
+        randoms = []
+        for arg in dist.args:
+            randoms.extend(random_symbols(arg))
+        if len(randoms) > 1:
+            raise NotImplementedError("Compound Distributions for more than"
+                " one random argument is not implemeted yet.")
+        rand_sym = randoms[0]
+        if isinstance(dist, SingleFiniteDistribution):
+            y = Dummy('y', integer=True, negative=False)
+            _pdf = dist.pmf(y)
+        else:
+            y = Dummy('y')
+            _pdf = dist.pdf(y)
+        if isinstance(rand_sym.pspace.distribution, SingleFiniteDistribution):
+            rand_dens = rand_sym.pspace.distribution.pmf(rand_sym)
+        else:
+            rand_dens = rand_sym.pspace.distribution.pdf(rand_sym)
+        rand_sym_dom = rand_sym.pspace.domain.set
+        if rand_sym.pspace.is_Discrete or rand_sym.pspace.is_Finite:
+            _pdf = Sum(_pdf*rand_dens, (rand_sym, rand_sym_dom._inf,
+                    rand_sym_dom._sup)).doit()
+        else:
+            _pdf = integrate(_pdf*rand_dens, (rand_sym, rand_sym_dom._inf,
+                    rand_sym_dom._sup))
+        return Lambda(y, _pdf)(x)
+
+    @classmethod
+    def _compound_check(self, dist):
+        """
+        Checks if the given distribution contains random parameters.
+        """
+        randoms = []
+        for arg in dist.args:
+            randoms.extend(random_symbols(arg))
+        if len(randoms) == 0:
+            return False
+        return True
diff --git a/sympy/stats/crv_types.py b/sympy/stats/crv_types.py
index dcfd78259f..958f43dd5f 100644
--- a/sympy/stats/crv_types.py
+++ b/sympy/stats/crv_types.py
@@ -65,7 +65,6 @@
                    hyper, uppergamma, sinh, Ne, expint, Rational, integrate)
 from sympy.matrices import MatrixBase, MatrixExpr
 from sympy.stats.crv import SingleContinuousPSpace, SingleContinuousDistribution
-from sympy.stats.joint_rv import JointPSpace, CompoundDistribution
 from sympy.stats.rv import _value_check, is_random
 
 oo = S.Infinity
@@ -134,7 +133,8 @@ def rv(symbol, cls, args):
     dist.check(*args)
     pspace = SingleContinuousPSpace(symbol, dist)
     if any(is_random(arg) for arg in args):
-        pspace = JointPSpace(symbol, CompoundDistribution(dist))
+        from sympy.stats.compound_rv import CompoundPSpace, CompoundDistribution
+        pspace = CompoundPSpace(symbol, CompoundDistribution(dist))
     return pspace.value
 
 
diff --git a/sympy/stats/drv_types.py b/sympy/stats/drv_types.py
index f32d9afd82..c0fb01946f 100644
--- a/sympy/stats/drv_types.py
+++ b/sympy/stats/drv_types.py
@@ -19,7 +19,6 @@
                    hyper, binomial, Piecewise, floor, besseli, sqrt, Sum, Dummy,
                    Lambda)
 from sympy.stats.drv import SingleDiscreteDistribution, SingleDiscretePSpace
-from sympy.stats.joint_rv import JointPSpace, CompoundDistribution
 from sympy.stats.rv import _value_check, is_random
 
 
@@ -40,7 +39,8 @@ def rv(symbol, cls, *args):
     dist.check(*args)
     pspace = SingleDiscretePSpace(symbol, dist)
     if any(is_random(arg) for arg in args):
-        pspace = JointPSpace(symbol, CompoundDistribution(dist))
+        from sympy.stats.compound_rv import CompoundPSpace, CompoundDistribution
+        pspace = CompoundPSpace(symbol, CompoundDistribution(dist))
     return pspace.value
 
 
diff --git a/sympy/stats/frv_types.py b/sympy/stats/frv_types.py
index 6933367ff2..ef89682832 100644
--- a/sympy/stats/frv_types.py
+++ b/sympy/stats/frv_types.py
@@ -40,7 +40,11 @@ def rv(name, cls, *args):
     args = list(map(sympify, args))
     dist = cls(*args)
     dist.check(*args)
-    return SingleFinitePSpace(name, dist).value
+    pspace = SingleFinitePSpace(name, dist)
+    if any(is_random(arg) for arg in args):
+        from sympy.stats.compound_rv import CompoundPSpace, CompoundDistribution
+        pspace = CompoundPSpace(name, CompoundDistribution(dist))
+    return pspace.value
 
 class FiniteDistributionHandmade(SingleFiniteDistribution):
 
@@ -261,8 +265,12 @@ def set(self):
         return set([self.succ, self.fail])
 
     def pmf(self, x):
-        return Piecewise((self.p, x == self.succ),
-                         (1 - self.p, x == self.fail),
+        if isinstance(self.succ, Symbol) and isinstance(self.fail, Symbol):
+            return Piecewise((self.p, x == self.succ),
+                             (1 - self.p, x == self.fail),
+                             (S.Zero, True))
+        return Piecewise((self.p, Eq(x, self.succ)),
+                         (1 - self.p, Eq(x, self.fail)),
                          (S.Zero, True))
 
 
diff --git a/sympy/stats/joint_rv.py b/sympy/stats/joint_rv.py
index 30a880f5fd..e6b6ac88f9 100644
--- a/sympy/stats/joint_rv.py
+++ b/sympy/stats/joint_rv.py
@@ -19,10 +19,8 @@
 from sympy.core.containers import Tuple
 from sympy.integrals.integrals import Integral, integrate
 from sympy.matrices import ImmutableMatrix
-from sympy.stats.crv import (ContinuousDistribution,
-                             SingleContinuousDistribution, SingleContinuousPSpace)
-from sympy.stats.drv import (DiscreteDistribution,
-                             SingleDiscreteDistribution, SingleDiscretePSpace)
+from sympy.stats.crv import SingleContinuousDistribution, SingleContinuousPSpace
+from sympy.stats.drv import SingleDiscreteDistribution, SingleDiscretePSpace
 from sympy.stats.rv import (ProductPSpace, NamedArgsMixin,
                             ProductDomain, RandomSymbol, random_symbols, SingleDomain)
 from sympy.utilities.misc import filldedent
@@ -194,42 +192,6 @@ def __getitem__(self, key):
             return Indexed(self, key)
 
 
-class CompoundDistribution(Basic, NamedArgsMixin):
-    """
-    Represents a compound probability distribution.
-
-    Constructed using a single probability distribution with a parameter
-    distributed according to some given distribution.
-    """
-    def __new__(cls, dist):
-        if not isinstance(dist, (ContinuousDistribution, DiscreteDistribution)):
-            raise ValueError(filldedent('''CompoundDistribution can only be
-             initialized from ContinuousDistribution or DiscreteDistribution
-             '''))
-        _args = dist.args
-        if not any([isinstance(i, RandomSymbol) for i in _args]):
-            return dist
-        return Basic.__new__(cls, dist)
-
-    @property
-    def latent_distributions(self):
-        return random_symbols(self.args[0])
-
-    def pdf(self, *x):
-        dist = self.args[0]
-        z = Dummy('z')
-        if isinstance(dist, ContinuousDistribution):
-            rv = SingleContinuousPSpace(z, dist).value
-        elif isinstance(dist, DiscreteDistribution):
-            rv = SingleDiscretePSpace(z, dist).value
-        return MarginalDistribution(self, (rv,)).pdf(*x)
-
-    def set(self):
-        return self.args[0].set
-
-    def __call__(self, *args):
-        return self.pdf(*args)
-
 
 class MarginalDistribution(Basic):
     """
@@ -268,10 +230,7 @@ def symbols(self):
     def pdf(self, *x):
         expr, rvs = self.args[0], self.args[1]
         marginalise_out = [i for i in random_symbols(expr) if i not in rvs]
-        if isinstance(expr, CompoundDistribution):
-            syms = Dummy('x', real=True)
-            expr = expr.args[0].pdf(syms)
-        elif isinstance(expr, JointDistribution):
+        if isinstance(expr, JointDistribution):
             count = len(expr.domain.args)
             x = Dummy('x', real=True, finite=True)
             syms = tuple(Indexed(x, i) for i in count)
diff --git a/sympy/stats/rv.py b/sympy/stats/rv.py
index 4e9fc1e5c0..88e46954cd 100644
--- a/sympy/stats/rv.py
+++ b/sympy/stats/rv.py
@@ -360,7 +360,8 @@ def __new__(cls, *spaces):
         symbols = FiniteSet(*[val.symbol for val in rs_space_dict.keys()])
 
         # Overlapping symbols
-        from sympy.stats.joint_rv import MarginalDistribution, CompoundDistribution
+        from sympy.stats.joint_rv import MarginalDistribution
+        from sympy.stats.compound_rv import CompoundDistribution
         if len(symbols) < sum(len(space.symbols) for space in spaces if not
          isinstance(space.distribution, (
             CompoundDistribution, MarginalDistribution))):
@@ -604,6 +605,10 @@ def pspace(expr):
     # If only one space present
     if all(rv.pspace == rvs[0].pspace for rv in rvs):
         return rvs[0].pspace
+    from sympy.stats.compound_rv import CompoundPSpace
+    for rv in rvs:
+        if isinstance(rv.pspace, CompoundPSpace):
+            return rv.pspace
     # Otherwise make a product space
     return IndependentProductPSpace(*[rv.pspace for rv in rvs])
 


Test patch:
diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py
index c9cac546ef..26b7ec759e 100644
--- a/sympy/core/tests/test_args.py
+++ b/sympy/core/tests/test_args.py
@@ -1090,12 +1090,21 @@ def test_sympy__stats__joint_rv__MarginalDistribution():
     assert _test_args(MarginalDistribution(r, (r,)))
 
 
-def test_sympy__stats__joint_rv__CompoundDistribution():
-    from sympy.stats.joint_rv import CompoundDistribution
-    from sympy.stats.drv_types import PoissonDistribution
-    r = PoissonDistribution(x)
+def test_sympy__stats__compound_rv__CompoundDistribution():
+    from sympy.stats.compound_rv import CompoundDistribution
+    from sympy.stats.drv_types import PoissonDistribution, Poisson
+    r = Poisson('r', 10)
     assert _test_args(CompoundDistribution(PoissonDistribution(r)))
 
+
+def test_sympy__stats__compound_rv__CompoundPSpace():
+    from sympy.stats.compound_rv import CompoundPSpace, CompoundDistribution
+    from sympy.stats.drv_types import PoissonDistribution, Poisson
+    r = Poisson('r', 5)
+    C = CompoundDistribution(PoissonDistribution(r))
+    assert _test_args(CompoundPSpace('C', C))
+
+
 @SKIP("abstract class")
 def test_sympy__stats__drv__SingleDiscreteDistribution():
     pass
diff --git a/sympy/stats/tests/test_compound_rv.py b/sympy/stats/tests/test_compound_rv.py
new file mode 100644
index 0000000000..5b32abeb1e
--- /dev/null
+++ b/sympy/stats/tests/test_compound_rv.py
@@ -0,0 +1,121 @@
+from sympy import (symbols, S, erf, sqrt, pi, exp, gamma, Interval, oo, beta,
+                    Eq, Piecewise, Integral, Abs, arg, Dummy)
+from sympy.stats import (Normal, P, E, density, Gamma, Poisson, Rayleigh,
+                        variance, Bernoulli, Beta, Uniform, cdf)
+from sympy.stats.compound_rv import CompoundDistribution, CompoundPSpace
+from sympy.stats.crv_types import NormalDistribution
+from sympy.stats.drv_types import PoissonDistribution
+from sympy.stats.frv_types import BernoulliDistribution
+from sympy.testing.pytest import raises
+from sympy.stats.joint_rv_types import MultivariateNormalDistribution
+
+
+x = symbols('x')
+def test_normal_CompoundDist():
+    X = Normal('X', 1, 2)
+    Y = Normal('X', X, 4)
+    assert density(Y)(x).simplify() == sqrt(10)*exp(-x**2/40 + x/20 - S(1)/40)/(20*sqrt(pi))
+    assert E(Y) == 1 # it is always equal to mean of X
+    assert P(Y > 1) == S(1)/2 # as 1 is the mean
+    assert P(Y > 5).simplify() ==  S(1)/2 - erf(sqrt(10)/5)/2
+    assert variance(Y) == variance(X) + 4**2 # 2**2 + 4**2
+    # https://math.stackexchange.com/questions/1484451/
+    # (Contains proof of E and variance computation)
+
+
+def test_poisson_CompoundDist():
+    k, t, y = symbols('k t y', positive=True, real=True)
+    G = Gamma('G', k, t)
+    D = Poisson('P', G)
+    assert density(D)(y).simplify() == t**y*(t + 1)**(-k - y)*gamma(k + y)/(gamma(k)*gamma(y + 1))
+    # https://en.wikipedia.org/wiki/Negative_binomial_distribution#Gamma%E2%80%93Poisson_mixture
+    assert E(D).simplify() == k*t # mean of NegativeBinomialDistribution
+
+
+def test_bernoulli_CompoundDist():
+    X = Beta('X', 1, 2)
+    Y = Bernoulli('Y', X)
+    assert density(Y).dict == {0: S(2)/3, 1: S(1)/3}
+    assert E(Y) == P(Eq(Y, 1)) == S(1)/3
+    assert variance(Y) == S(2)/9
+    assert cdf(Y) == {0: S(2)/3, 1: 1}
+
+    # test issue 8128
+    a = Bernoulli('a', S(1)/2)
+    b = Bernoulli('b', a)
+    assert density(b).dict == {0: S(1)/2, 1: S(1)/2}
+    assert P(b > 0.5) == S(1)/2
+
+    X = Uniform('X', 0, 1)
+    Y = Bernoulli('Y', X)
+    assert E(Y) == S(1)/2
+    assert P(Eq(Y, 1)) == E(Y)
+
+
+def test_unevaluated_CompoundDist():
+    # these tests need to be removed once they work with evaluation as they are currently not
+    # evaluated completely in sympy.
+    R = Rayleigh('R', 4)
+    X = Normal('X', 3, R)
+    _k = Dummy('k')
+    exprd = Piecewise((exp(S(3)/4 - x/4)/8, 2*Abs(arg(x - 3)) <= pi/2),
+    (sqrt(2)*Integral(exp(-(_k**4 + 16*(x - 3)**2)/(32*_k**2)),
+    (_k, 0, oo))/(32*sqrt(pi)), True))
+    assert (density(X)(x).simplify()).dummy_eq(exprd.simplify())
+
+    expre = Integral(Piecewise((_k*exp(S(3)/4 - _k/4)/8, 2*Abs(arg(_k - 3)) <= pi/2),
+    (sqrt(2)*_k*Integral(exp(-(_k**4 + 16*(_k - 3)**2)/(32*_k**2)),
+    (_k, 0, oo))/(32*sqrt(pi)), True)), (_k, -oo, oo))
+    assert (E(X, evaluate=False).simplify()).dummy_eq(expre.simplify())
+
+
+def test_Compound_Distribution():
+    X = Normal('X', 2, 4)
+    N = NormalDistribution(X, 4)
+    C = CompoundDistribution(N)
+    assert C.is_Continuous
+    assert C.set == Interval(-oo, oo)
+    assert C.pdf(x).simplify() == exp(-x**2/64 + x/16 - S(1)/16)/(8*sqrt(pi))
+
+    assert not isinstance(CompoundDistribution(NormalDistribution(2, 3)),
+                            CompoundDistribution)
+    M = MultivariateNormalDistribution([1, 2], [[2, 1], [1, 2]])
+    raises(NotImplementedError, lambda: CompoundDistribution(M))
+
+    X = Beta('X', 2, 4)
+    B = BernoulliDistribution(X, 1, 0)
+    C = CompoundDistribution(B)
+    assert C.is_Finite
+    assert C.set == {0, 1}
+    y = symbols('y', negative=False, integer=True)
+    assert C.pdf(y) == Piecewise((S(1)/(30*beta(2, 4)), Eq(y, 0)),
+                (S(1)/(60*beta(2, 4)), Eq(y, 1)), (0, True))
+
+    k, t, z = symbols('k t z', positive=True, real=True)
+    G = Gamma('G', k, t)
+    X = PoissonDistribution(G)
+    C = CompoundDistribution(X)
+    assert C.is_Discrete
+    assert C.set == S.Naturals0
+    assert C.pdf(z).simplify() == t**z*(t + 1)**(-k - z)*gamma(k + z)/(gamma(k)*gamma(z + 1))
+
+
+def test_compound_pspace():
+    X = Normal('X', 2, 4)
+    Y = Normal('Y', 3, 6)
+    assert not isinstance(Y.pspace, CompoundPSpace)
+    N = NormalDistribution(1, 2)
+    D = PoissonDistribution(3)
+    B = BernoulliDistribution(0.2, 1, 0)
+    pspace1 = CompoundPSpace('N', N)
+    pspace2 = CompoundPSpace('D', D)
+    pspace3 = CompoundPSpace('B', B)
+    assert not isinstance(pspace1, CompoundPSpace)
+    assert not isinstance(pspace2, CompoundPSpace)
+    assert not isinstance(pspace3, CompoundPSpace)
+    M = MultivariateNormalDistribution([1, 2], [[2, 1], [1, 2]])
+    raises(ValueError, lambda: CompoundPSpace('M', M))
+    Y = Normal('Y', X, 6)
+    assert isinstance(Y.pspace, CompoundPSpace)
+    assert Y.pspace.distribution == CompoundDistribution(NormalDistribution(X, 6))
+    assert Y.pspace.domain.set == Interval(-oo, oo)
diff --git a/sympy/stats/tests/test_continuous_rv.py b/sympy/stats/tests/test_continuous_rv.py
index 0f08e7b7eb..1e249c68ea 100644
--- a/sympy/stats/tests/test_continuous_rv.py
+++ b/sympy/stats/tests/test_continuous_rv.py
@@ -24,7 +24,7 @@
 from sympy.stats.crv_types import NormalDistribution, ExponentialDistribution, ContinuousDistributionHandmade
 from sympy.stats.joint_rv_types import MultivariateLaplaceDistribution, MultivariateNormalDistribution
 from sympy.stats.crv import SingleContinuousPSpace, SingleContinuousDomain
-from sympy.stats.joint_rv import JointPSpace
+from sympy.stats.compound_rv import CompoundPSpace
 from sympy.testing.pytest import raises, XFAIL, slow, skip, ignore_warnings
 from sympy.testing.randtest import verify_numerically as tn
 
@@ -1395,11 +1395,10 @@ def test_random_parameters():
     mu = Normal('mu', 2, 3)
     meas = Normal('T', mu, 1)
     assert density(meas, evaluate=False)(z)
-    assert isinstance(pspace(meas), JointPSpace)
+    assert isinstance(pspace(meas), CompoundPSpace)
     X = Normal('x', [1, 2], [[1, 0], [0, 1]])
     assert isinstance(pspace(X).distribution, MultivariateNormalDistribution)
-    #assert density(meas, evaluate=False)(z) == Integral(mu.pspace.pdf *
-    #        meas.pspace.pdf, (mu.symbol, -oo, oo)).subs(meas.symbol, z)
+    assert density(meas)(z).simplify() == sqrt(5)*exp(-z**2/20 + z/5 - S(1)/5)/(10*sqrt(pi))
 
 
 def test_random_parameters_given():
diff --git a/sympy/stats/tests/test_mix.py b/sympy/stats/tests/test_mix.py
index 85e41651ee..9cbd336646 100644
--- a/sympy/stats/tests/test_mix.py
+++ b/sympy/stats/tests/test_mix.py
@@ -6,7 +6,8 @@
                         Multinomial, MultivariateBeta)
 from sympy.stats.crv_types import Normal
 from sympy.stats.drv_types import PoissonDistribution
-from sympy.stats.joint_rv import JointPSpace, CompoundDistribution, MarginalDistribution
+from sympy.stats.compound_rv import CompoundPSpace, CompoundDistribution
+from sympy.stats.joint_rv import MarginalDistribution
 from sympy.stats.rv import pspace, density
 
 def test_density():
@@ -14,13 +15,14 @@ def test_density():
     l = Symbol('l', positive=True)
     rate = Beta(l, 2, 3)
     X = Poisson(x, rate)
-    assert isinstance(pspace(X), JointPSpace)
+    assert isinstance(pspace(X), CompoundPSpace)
     assert density(X, Eq(rate, rate.symbol)) == PoissonDistribution(l)
     N1 = Normal('N1', 0, 1)
     N2 = Normal('N2', N1, 2)
     assert density(N2)(0).doit() == sqrt(10)/(10*sqrt(pi))
     assert simplify(density(N2, Eq(N1, 1))(x)) == \
         sqrt(2)*exp(-(x - 1)**2/8)/(4*sqrt(pi))
+    assert simplify(density(N2)(x)) == sqrt(10)*exp(-x**2/10)/(10*sqrt(pi))
 
 def test_MarginalDistribution():
     a1, p1, p2 = symbols('a1 p1 p2', positive=True)
@@ -48,7 +50,7 @@ def test_MarginalDistribution():
 def test_compound_distribution():
     Y = Poisson('Y', 1)
     Z = Poisson('Z', Y)
-    assert isinstance(pspace(Z), JointPSpace)
+    assert isinstance(pspace(Z), CompoundPSpace)
     assert isinstance(pspace(Z).distribution, CompoundDistribution)
     assert Z.pspace.distribution.pdf(1).doit() == exp(-2)*exp(exp(-1))