如何使用 sympy 创建多元正态密度?
How do I create a multivariate normal density with sympy?
我在使用 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
问题很简单:计算密度的公式不是支撑矩阵,看看:
https://github.com/sympy/sympy/blob/sympy-0.7.6.1/sympy/stats/crv_types.py#L1641
在这个表达式中,(x-self.mean)得到平方(即2的次方),但非方阵的平方不是定义。
简而言之,似乎不支持多元正态分布,但您可以通过定义新分布来尝试解决方法:
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 ]
我在使用 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
问题很简单:计算密度的公式不是支撑矩阵,看看:
https://github.com/sympy/sympy/blob/sympy-0.7.6.1/sympy/stats/crv_types.py#L1641
在这个表达式中,(x-self.mean)得到平方(即2的次方),但非方阵的平方不是定义。
简而言之,似乎不支持多元正态分布,但您可以通过定义新分布来尝试解决方法:
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 ]