Как создать многомерную нормальную плотность с sympy?

Question

У меня возникла проблема создания многомерной нормальной плотности с sympy 0.7.6.1.

Вот мой код.

from sympy import * 
from sympy.stats import * 

mu = Matrix([5, 13]) 
Sigma = Matrix([[2, 0], [0, 2]]) 
X = Normal('X', mu, Sigma) 
y = MatrixSymbol('y', 2, 1) 
density(X)(y) 

Последняя строка дает мне эту ошибку:

Power of non-square matrix Matrix([ 
[ -5], 
[-13]]) + y 

Answer 1

Проблема проста: формула для расчета плотности не один поддерживает матрицы, имеют вид:

https://github.com/sympy/sympy/blob/sympy-0.7.6.1/sympy/stats/crv_types.py#L1641

В этом выражении (x-self.mean) получает квадрат (т.е. поднят до степени 2), но квадрат o f неквадратичная матрица не определена.

Короче говоря, это выглядит как многомерные нормальные распределения не поддерживаются, но вы можете попробовать обходной путь, определив новое распределение:

from sympy.stats.crv_types import rv, SingleContinuousDistribution, _value_check 

class MultivariateNormalDistribution(SingleContinuousDistribution): 
     _argnames = ('mean', 'std') 
     @staticmethod 
     def check(mean, std): 
       _value_check(std > 0, "Standard deviation must be positive") 
     def pdf(self, x): 
       return exp(-S.Half * (x - self.mean).T * (self.std.inv()) * (x - self.mean))/(sqrt(2*pi)**(self.std.shape[0])*self.std.det()) 
     def sample(self): 
       pass 
       # define sampling function here 

def MultivariateNormal(name, mean, std): 
     return rv(name, MultivariateNormalDistribution, (mean, std)) 

К сожалению, ваш пример еще не работает, из-за отсутствия признаков в модуле матрицы (то есть, не экспоненцирование выражений с MatrixSymbol не поддерживаются, пока), но вы можете получить плотность точек:

In[12]: X = MultivariateNormal('X', mu, Sigma) 

In [13]: density(X)(Matrix([0, 0])) 
Out[13]: 
[ -97/2] 
[e  ] 
[------] 
[ 8*pi ] 

Или с символами в матрице:

In [14]: x1, x2 = symbols('x1, x2') 

In [15]: density(X)(Matrix([x1, x2])) 
Out[15]: 
[  2   2    ] 
[ x1 5*x1 x2 13*x2 97] 
[ - --- + ---- - --- + ----- - --] 
[ 4  2  4  2  2 ] 
[e        ] 
[--------------------------------] 
[    8*pi    ]