def set(self):
     rvs = [i for i in random_symbols(self.args[1])]
     marginalise_out = [i for i in random_symbols(self.args[1]) \
      if i not in self.args[1]]
     for i in rvs:
         if i in marginalise_out:
             rvs.remove(i)
     return ProductSet((i.pspace.set for i in rvs))

 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.intersect(self.domain.set)
     return SingleContinuousDomain(rv.symbol, interval)

 def __new__(cls,dist, rvs):
     if not all([isinstance(rv, (Indexed, RandomSymbol))] for rv in rvs):
         raise ValueError(filldedent('''Marginal distribution can be
          intitialised only in terms of random variables or indexed random
          variables'''))
     rvs = Tuple.fromiter(rv for rv in rvs)
     if not isinstance(dist, JointDistribution) and len(random_symbols(dist)) == 0:
         return dist
     return Basic.__new__(cls, dist, rvs)

 def restricted_domain(self, condition):
     rvs = random_symbols(condition)
     assert all(r.symbol in self.symbols for r in rvs)
     if (len(rvs) > 1):
         raise NotImplementedError(filldedent('''Multivariate discrete
         random variables are not yet supported.'''))
     conditional_domain = reduce_rational_inequalities_wrap(condition,
         rvs[0])
     conditional_domain = conditional_domain.intersect(self.domain.set)
     return conditional_domain

 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 self.args[1]]
     syms = [i.pspace.symbol for i in self.args[1]]
     for i in expr.atoms(Indexed):
         if isinstance(i, Indexed) and isinstance(i.base, RandomSymbol)\
          and i not in rvs:
             marginalise_out.append(i)
     if isinstance(expr, CompoundDistribution):
         syms = Dummy('x', real=True)
         expr = expr.args[0].pdf(syms)
     elif isinstance(expr, JointDistribution):
         count = len(expr.domain.args)
         x = Dummy('x', real=True, finite=True)
         syms = [Indexed(x, i) for i in count]
         expr = expression.pdf(syms)
     return Lambda(syms, self.compute_pdf(expr, marginalise_out))(*x)

 def compute_expectation(self, expr, rvs=None, evaluate=False, **kwargs):
     syms = tuple(self.value[i] for i in range(self.component_count))
     rvs = rvs or syms
     if not any([i in rvs for i in syms]):
         return expr
     expr = expr*self.pdf
     for rv in rvs:
         if isinstance(rv, Indexed):
             expr = expr.xreplace({rv: Indexed(str(rv.base), rv.args[1])})
         elif isinstance(rv, RandomSymbol):
             expr = expr.xreplace({rv: rv.symbol})
     if self.value in random_symbols(expr):
         raise NotImplementedError(filldedent('''
         Expectations of expression with unindexed joint random symbols
         cannot be calculated yet.'''))
     limits = tuple((Indexed(str(rv.base),rv.args[1]),
         self.distribution.set.args[rv.args[1]]) for rv in syms)
     return Integral(expr, *limits)

def JointRV(symbol, pdf, _set=None):
    """
    Create a Joint Random Variable where each of its component is conitinuous,
    given the following:

    -- a symbol
    -- a PDF in terms of indexed symbols of the symbol given
     as the first argument

    NOTE: As of now, the set for each component for a `JointRV` is
    equal to the set of all integers, which can not be changed.

    Returns a RandomSymbol.

    Examples
    ========

    >>> from sympy import symbols, exp, pi, Indexed, S
    >>> from sympy.stats import density
    >>> from sympy.stats.joint_rv_types import JointRV

    >>> x1, x2 = (Indexed('x', i) for i in (1, 2))
    >>> pdf = exp(-x1**2/2 + x1 - x2**2/2 - S(1)/2)/(2*pi)

    >>> N1 = JointRV('x', pdf) #Multivariate Normal distribution
    >>> density(N1)(1, 2)
    exp(-2)/(2*pi)
    """
    #TODO: Add support for sets provided by the user
    symbol = sympify(symbol)
    syms = list(i for i in pdf.free_symbols if isinstance(i, Indexed)
        and i.base == IndexedBase(symbol))
    syms.sort(key = lambda index: index.args[1])
    _set = S.Reals**len(syms)
    pdf = Lambda(syms, pdf)
    dist = JointDistributionHandmade(pdf, _set)
    jrv = JointPSpace(symbol, dist).value
    rvs = random_symbols(pdf)
    if len(rvs) != 0:
        dist = MarginalDistribution(dist, (jrv,))
        return JointPSpace(symbol, dist).value
    return jrv

 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 where(self, condition):
     assert all(r.symbol in self.symbols for r in random_symbols(condition))
     return ConditionalFiniteDomain(self.domain, condition)

 def domain(self):
     rvs = random_symbols(self.distribution)
     if len(rvs) == 0:
         return SingleDomain(self.symbol, self.set)
     return ProductDomain(*[rv.pspace.domain for rv in rvs])

 def latent_distributions(self):
     return random_symbols(self.args[0])