I am trying to work with the multivariate normal distribution in symbolic form using SymPy. However, it doesn't seem capable of realizing that the matrix expression I want to exponentiate evaluates to a scalar. Here's the code:
from sympy import *
# Sample size, number of covariates
n, k = symbols('n k')
# Data
X = MatrixSymbol('X', n, k)
y = MatrixSymbol('y', n, 1)
# Parameters
beta = MatrixSymbol('beta', k, 1)
Sigma = MatrixSymbol('Sigma', n, n)
# Exponent expression, OK:
(-Rational(1, 2)*(y - X*beta).T*Sigma*(y - X*beta))
# Trying to use as exponent, not OK:
exp(-Rational(1, 2)*(y - X*beta).T*Sigma*(y - X*beta))
The error message that follows:
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
<ipython-input-76-62c854b75962> in <module>()
----> 1 exp((-Rational(1, 2)*(y - X*beta).T*Sigma*(y - X*beta)))
D:\Programming\Anaconda3\lib\site-packages\sympy\core\function.py in __new__(cls, *args, **options)
425
426 evaluate = options.get('evaluate', global_evaluate[0])
--> 427 result = super(Function, cls).__new__(cls, *args, **options)
428 if not evaluate or not isinstance(result, cls):
429 return result
D:\Programming\Anaconda3\lib\site-packages\sympy\core\function.py in __new__(cls, *args, **options)
248
249 if evaluate:
--> 250 evaluated = cls.eval(*args)
251 if evaluated is not None:
252 return evaluated
D:\Programming\Anaconda3\lib\site-packages\sympy\functions\elementary\exponential.py in eval(cls, arg)
300
301 elif arg.is_Matrix:
--> 302 return arg.exp()
303
304 @property
AttributeError: 'MatMul' object has no attribute 'exp'
The exponential function does not support matrix expressions; it expects an explicit matrix, meaning a matrix of definite size with definite entries (which can be symbolic). To obtain an explicit matrix from a matrix expression, use as_explicit method.
exp((-Rational(1, 2)*(y - X*beta).T*Sigma*(y - X*beta)).as_explicit())
returns a 1 by 1 matrix, namely
Matrix([[exp(Sum(-y[_k, 0]*Sum(Sigma[_k, _k]*y[_k, 0], (_k, 0, n - 1)) + y[_k, 0]*Sum(Sigma[_k, _k]*Sum(X[_k, _k]*beta[_k, 0], (_k, 0, k - 1)), (_k, 0, n - 1)) + Sum(Sigma[_k, _k]*y[_k, 0], (_k, 0, n - 1))*Sum(X[_k, _k]*beta[_k, 0], (_k, 0, k - 1)) - Sum(Sigma[_k, _k]*Sum(X[_k, _k]*beta[_k, 0], (_k, 0, k - 1)), (_k, 0, n - 1))*Sum(X[_k, _k]*beta[_k, 0], (_k, 0, k - 1)), (_k, 0, n - 1))/2)]])
Or you can do
exp((-Rational(1, 2)*(y - X*beta).T*Sigma*(y - X*beta)).as_explicit()[0, 0])
dropping from a 1 by 1 matrix to a scalar quantity.