Repository: sympy/sympy

Issue:
MultivariateNormal density calculation fails when matrix params are used
The following code fails:

```python
import sympy as sy
import sympy.stats

S = sy.MatrixSymbol('S', 3, 3)
m = sy.MatrixSymbol('m', 3, 1)

dist = sy.stats.MultivariateNormal('X', m, S)
x = sy.MatrixSymbol('x', 3, 1)

sy.stats.density(dist)(x)
```
with error:
```
~/miniconda3/lib/python3.7/site-packages/sympy/stats/joint_rv_types.py in pdf(self, *args)
    151         return  S.One/sqrt((2*pi)**(k)*det(sigma))*exp(
    152             Rational(-1, 2)*x.transpose()*(sigma.inv()*\
--> 153                 x))[0]
    154 
    155     def _marginal_distribution(self, indices, sym):

TypeError: 'exp' object is not subscriptable
```

Is this intended behaviour? 

Discussion:

Patch:
diff --git a/sympy/stats/joint_rv_types.py b/sympy/stats/joint_rv_types.py
index 1f5c06d7c6..51ccb7089b 100644
--- a/sympy/stats/joint_rv_types.py
+++ b/sympy/stats/joint_rv_types.py
@@ -4,6 +4,7 @@
                    Intersection, Matrix, symbols, Product, IndexedBase)
 from sympy.matrices import ImmutableMatrix, MatrixSymbol
 from sympy.matrices.expressions.determinant import det
+from sympy.matrices.expressions.matexpr import MatrixElement
 from sympy.stats.joint_rv import JointDistribution, JointPSpace, MarginalDistribution
 from sympy.stats.rv import _value_check, random_symbols
 
@@ -146,11 +147,14 @@ def check(mu, sigma):
     def pdf(self, *args):
         mu, sigma = self.mu, self.sigma
         k = mu.shape[0]
-        args = ImmutableMatrix(args)
+        if len(args) == 1 and args[0].is_Matrix:
+            args = args[0]
+        else:
+            args = ImmutableMatrix(args)
         x = args - mu
-        return  S.One/sqrt((2*pi)**(k)*det(sigma))*exp(
-            Rational(-1, 2)*x.transpose()*(sigma.inv()*\
-                x))[0]
+        density = S.One/sqrt((2*pi)**(k)*det(sigma))*exp(
+            Rational(-1, 2)*x.transpose()*(sigma.inv()*x))
+        return MatrixElement(density, 0, 0)
 
     def _marginal_distribution(self, indices, sym):
         sym = ImmutableMatrix([Indexed(sym, i) for i in indices])
@@ -188,7 +192,7 @@ def MultivariateNormal(name, mu, sigma):
     ========
 
     >>> from sympy.stats import MultivariateNormal, density, marginal_distribution
-    >>> from sympy import symbols
+    >>> from sympy import symbols, MatrixSymbol
     >>> X = MultivariateNormal('X', [3, 4], [[2, 1], [1, 2]])
     >>> y, z = symbols('y z')
     >>> density(X)(y, z)
@@ -200,6 +204,21 @@ def MultivariateNormal(name, mu, sigma):
     >>> marginal_distribution(X, X[0])(y)
     exp(-(y - 3)**2/4)/(2*sqrt(pi))
 
+    The example below shows that it is also possible to use
+    symbolic parameters to define the MultivariateNormal class.
+
+    >>> n = symbols('n', natural=True)
+    >>> Sg = MatrixSymbol('Sg', n, n)
+    >>> mu = MatrixSymbol('mu', n, 1)
+    >>> obs = MatrixSymbol('obs', n, 1)
+    >>> X = MultivariateNormal('X', mu, Sg)
+
+    The density of a multivariate normal can be
+    calculated using a matrix argument, as shown below.
+
+    >>> density(X)(obs)
+    (exp(((1/2)*mu.T - (1/2)*obs.T)*Sg**(-1)*(-mu + obs))/sqrt((2*pi)**n*Determinant(Sg)))[0, 0]
+
     References
     ==========
 


Test patch:
diff --git a/sympy/stats/tests/test_joint_rv.py b/sympy/stats/tests/test_joint_rv.py
index fed1eb1c61..13da27b000 100644
--- a/sympy/stats/tests/test_joint_rv.py
+++ b/sympy/stats/tests/test_joint_rv.py
@@ -1,9 +1,11 @@
 from sympy import (symbols, pi, oo, S, exp, sqrt, besselk, Indexed, Sum, simplify,
                    Rational, factorial, gamma, Piecewise, Eq, Product, Interval,
-                   IndexedBase, RisingFactorial, polar_lift, ProductSet, Range)
+                   IndexedBase, RisingFactorial, polar_lift, ProductSet, Range, eye,
+                   Determinant)
 from sympy.core.numbers import comp
 from sympy.integrals.integrals import integrate
 from sympy.matrices import Matrix, MatrixSymbol
+from sympy.matrices.expressions.matexpr import MatrixElement
 from sympy.stats import density, median, marginal_distribution, Normal, Laplace, E, sample
 from sympy.stats.joint_rv_types import (JointRV, MultivariateNormalDistribution,
                 JointDistributionHandmade, MultivariateT, NormalGamma,
@@ -46,6 +48,38 @@ def test_Normal():
     # assert E(X) == mu
     # assert variance(X) == sigma
 
+    # test symbolic multivariate normal densities
+    n = 3
+
+    Sg = MatrixSymbol('Sg', n, n)
+    mu = MatrixSymbol('mu', n, 1)
+    obs = MatrixSymbol('obs', n, 1)
+
+    X = MultivariateNormal('X', mu, Sg)
+    density_X = density(X)
+
+    eval_a = density_X(obs).subs({Sg: eye(3),
+        mu: Matrix([0, 0, 0]), obs: Matrix([0, 0, 0])}).doit()
+    eval_b = density_X(0, 0, 0).subs({Sg: eye(3), mu: Matrix([0, 0, 0])}).doit()
+
+    assert eval_a == sqrt(2)/(4*pi**Rational(3/2))
+    assert eval_b == sqrt(2)/(4*pi**Rational(3/2))
+
+    n = symbols('n', natural=True)
+
+    Sg = MatrixSymbol('Sg', n, n)
+    mu = MatrixSymbol('mu', n, 1)
+    obs = MatrixSymbol('obs', n, 1)
+
+    X = MultivariateNormal('X', mu, Sg)
+    density_X_at_obs = density(X)(obs)
+
+    expected_density = MatrixElement(
+        exp((S(1)/2) * (mu.T - obs.T) * Sg**(-1) * (-mu + obs)) / \
+        sqrt((2*pi)**n * Determinant(Sg)), 0, 0)
+
+    assert density_X_at_obs == expected_density
+
 def test_MultivariateTDist():
     t1 = MultivariateT('T', [0, 0], [[1, 0], [0, 1]], 2)
     assert(density(t1))(1, 1) == 1/(8*pi)