187 lines
8.6 KiB
Python
187 lines
8.6 KiB
Python
|
from sympy.concrete.products import Product
|
||
|
|
from sympy.core.numbers import pi
|
||
|
|
from sympy.core.singleton import S
|
||
|
|
from sympy.core.symbol import (Dummy, symbols)
|
||
|
|
from sympy.functions.elementary.exponential import exp
|
||
|
|
from sympy.functions.elementary.miscellaneous import sqrt
|
||
|
|
from sympy.functions.special.gamma_functions import gamma
|
||
|
|
from sympy.matrices import Determinant, Matrix, Trace, MatrixSymbol, MatrixSet
|
||
|
|
from sympy.stats import density, sample
|
||
|
|
from sympy.stats.matrix_distributions import (MatrixGammaDistribution,
|
||
|
|
MatrixGamma, MatrixPSpace, Wishart, MatrixNormal, MatrixStudentT)
|
||
|
|
from sympy.testing.pytest import raises, skip
|
||
|
|
from sympy.external import import_module
|
||
|
|
|
||
|
|
|
||
|
|
def test_MatrixPSpace():
|
||
|
|
M = MatrixGammaDistribution(1, 2, [[2, 1], [1, 2]])
|
||
|
|
MP = MatrixPSpace('M', M, 2, 2)
|
||
|
|
assert MP.distribution == M
|
||
|
|
raises(ValueError, lambda: MatrixPSpace('M', M, 1.2, 2))
|
||
|
|
|
||
|
|
def test_MatrixGamma():
|
||
|
|
M = MatrixGamma('M', 1, 2, [[1, 0], [0, 1]])
|
||
|
|
assert M.pspace.distribution.set == MatrixSet(2, 2, S.Reals)
|
||
|
|
assert isinstance(density(M), MatrixGammaDistribution)
|
||
|
|
X = MatrixSymbol('X', 2, 2)
|
||
|
|
num = exp(Trace(Matrix([[-S(1)/2, 0], [0, -S(1)/2]])*X))
|
||
|
|
assert density(M)(X).doit() == num/(4*pi*sqrt(Determinant(X)))
|
||
|
|
assert density(M)([[2, 1], [1, 2]]).doit() == sqrt(3)*exp(-2)/(12*pi)
|
||
|
|
X = MatrixSymbol('X', 1, 2)
|
||
|
|
Y = MatrixSymbol('Y', 1, 2)
|
||
|
|
assert density(M)([X, Y]).doit() == exp(-X[0, 0]/2 - Y[0, 1]/2)/(4*pi*sqrt(
|
||
|
|
X[0, 0]*Y[0, 1] - X[0, 1]*Y[0, 0]))
|
||
|
|
# symbolic
|
||
|
|
a, b = symbols('a b', positive=True)
|
||
|
|
d = symbols('d', positive=True, integer=True)
|
||
|
|
Y = MatrixSymbol('Y', d, d)
|
||
|
|
Z = MatrixSymbol('Z', 2, 2)
|
||
|
|
SM = MatrixSymbol('SM', d, d)
|
||
|
|
M2 = MatrixGamma('M2', a, b, SM)
|
||
|
|
M3 = MatrixGamma('M3', 2, 3, [[2, 1], [1, 2]])
|
||
|
|
k = Dummy('k')
|
||
|
|
exprd = pi**(-d*(d - 1)/4)*b**(-a*d)*exp(Trace((-1/b)*SM**(-1)*Y)
|
||
|
|
)*Determinant(SM)**(-a)*Determinant(Y)**(a - d/2 - S(1)/2)/Product(
|
||
|
|
gamma(-k/2 + a + S(1)/2), (k, 1, d))
|
||
|
|
assert density(M2)(Y).dummy_eq(exprd)
|
||
|
|
raises(NotImplementedError, lambda: density(M3 + M)(Z))
|
||
|
|
raises(ValueError, lambda: density(M)(1))
|
||
|
|
raises(ValueError, lambda: MatrixGamma('M', -1, 2, [[1, 0], [0, 1]]))
|
||
|
|
raises(ValueError, lambda: MatrixGamma('M', -1, -2, [[1, 0], [0, 1]]))
|
||
|
|
raises(ValueError, lambda: MatrixGamma('M', -1, 2, [[1, 0], [2, 1]]))
|
||
|
|
raises(ValueError, lambda: MatrixGamma('M', -1, 2, [[1, 0], [0]]))
|
||
|
|
|
||
|
|
def test_Wishart():
|
||
|
|
W = Wishart('W', 5, [[1, 0], [0, 1]])
|
||
|
|
assert W.pspace.distribution.set == MatrixSet(2, 2, S.Reals)
|
||
|
|
X = MatrixSymbol('X', 2, 2)
|
||
|
|
term1 = exp(Trace(Matrix([[-S(1)/2, 0], [0, -S(1)/2]])*X))
|
||
|
|
assert density(W)(X).doit() == term1 * Determinant(X)/(24*pi)
|
||
|
|
assert density(W)([[2, 1], [1, 2]]).doit() == exp(-2)/(8*pi)
|
||
|
|
n = symbols('n', positive=True)
|
||
|
|
d = symbols('d', positive=True, integer=True)
|
||
|
|
Y = MatrixSymbol('Y', d, d)
|
||
|
|
SM = MatrixSymbol('SM', d, d)
|
||
|
|
W = Wishart('W', n, SM)
|
||
|
|
k = Dummy('k')
|
||
|
|
exprd = 2**(-d*n/2)*pi**(-d*(d - 1)/4)*exp(Trace(-(S(1)/2)*SM**(-1)*Y)
|
||
|
|
)*Determinant(SM)**(-n/2)*Determinant(Y)**(
|
||
|
|
-d/2 + n/2 - S(1)/2)/Product(gamma(-k/2 + n/2 + S(1)/2), (k, 1, d))
|
||
|
|
assert density(W)(Y).dummy_eq(exprd)
|
||
|
|
raises(ValueError, lambda: density(W)(1))
|
||
|
|
raises(ValueError, lambda: Wishart('W', -1, [[1, 0], [0, 1]]))
|
||
|
|
raises(ValueError, lambda: Wishart('W', -1, [[1, 0], [2, 1]]))
|
||
|
|
raises(ValueError, lambda: Wishart('W', 2, [[1, 0], [0]]))
|
||
|
|
|
||
|
|
def test_MatrixNormal():
|
||
|
|
M = MatrixNormal('M', [[5, 6]], [4], [[2, 1], [1, 2]])
|
||
|
|
assert M.pspace.distribution.set == MatrixSet(1, 2, S.Reals)
|
||
|
|
X = MatrixSymbol('X', 1, 2)
|
||
|
|
term1 = exp(-Trace(Matrix([[ S(2)/3, -S(1)/3], [-S(1)/3, S(2)/3]])*(
|
||
|
|
Matrix([[-5], [-6]]) + X.T)*Matrix([[S(1)/4]])*(Matrix([[-5, -6]]) + X))/2)
|
||
|
|
assert density(M)(X).doit() == (sqrt(3)) * term1/(24*pi)
|
||
|
|
assert density(M)([[7, 8]]).doit() == sqrt(3)*exp(-S(1)/3)/(24*pi)
|
||
|
|
d, n = symbols('d n', positive=True, integer=True)
|
||
|
|
SM2 = MatrixSymbol('SM2', d, d)
|
||
|
|
SM1 = MatrixSymbol('SM1', n, n)
|
||
|
|
LM = MatrixSymbol('LM', n, d)
|
||
|
|
Y = MatrixSymbol('Y', n, d)
|
||
|
|
M = MatrixNormal('M', LM, SM1, SM2)
|
||
|
|
exprd = (2*pi)**(-d*n/2)*exp(-Trace(SM2**(-1)*(-LM.T + Y.T)*SM1**(-1)*(-LM + Y)
|
||
|
|
)/2)*Determinant(SM1)**(-d/2)*Determinant(SM2)**(-n/2)
|
||
|
|
assert density(M)(Y).doit() == exprd
|
||
|
|
raises(ValueError, lambda: density(M)(1))
|
||
|
|
raises(ValueError, lambda: MatrixNormal('M', [1, 2], [[1, 0], [0, 1]], [[1, 0], [2, 1]]))
|
||
|
|
raises(ValueError, lambda: MatrixNormal('M', [1, 2], [[1, 0], [2, 1]], [[1, 0], [0, 1]]))
|
||
|
|
raises(ValueError, lambda: MatrixNormal('M', [1, 2], [[1, 0], [0, 1]], [[1, 0], [0, 1]]))
|
||
|
|
raises(ValueError, lambda: MatrixNormal('M', [1, 2], [[1, 0], [2]], [[1, 0], [0, 1]]))
|
||
|
|
raises(ValueError, lambda: MatrixNormal('M', [1, 2], [[1, 0], [2, 1]], [[1, 0], [0]]))
|
||
|
|
raises(ValueError, lambda: MatrixNormal('M', [[1, 2]], [[1, 0], [0, 1]], [[1, 0]]))
|
||
|
|
raises(ValueError, lambda: MatrixNormal('M', [[1, 2]], [1], [[1, 0]]))
|
||
|
|
|
||
|
|
def test_MatrixStudentT():
|
||
|
|
M = MatrixStudentT('M', 2, [[5, 6]], [[2, 1], [1, 2]], [4])
|
||
|
|
assert M.pspace.distribution.set == MatrixSet(1, 2, S.Reals)
|
||
|
|
X = MatrixSymbol('X', 1, 2)
|
||
|
|
D = pi ** (-1.0) * Determinant(Matrix([[4]])) ** (-1.0) * Determinant(Matrix([[2, 1], [1, 2]])) \
|
||
|
|
** (-0.5) / Determinant(Matrix([[S(1) / 4]]) * (Matrix([[-5, -6]]) + X)
|
||
|
|
* Matrix([[S(2) / 3, -S(1) / 3], [-S(1) / 3, S(2) / 3]]) * (
|
||
|
|
Matrix([[-5], [-6]]) + X.T) + Matrix([[1]])) ** 2
|
||
|
|
assert density(M)(X) == D
|
||
|
|
|
||
|
|
v = symbols('v', positive=True)
|
||
|
|
n, p = 1, 2
|
||
|
|
Omega = MatrixSymbol('Omega', p, p)
|
||
|
|
Sigma = MatrixSymbol('Sigma', n, n)
|
||
|
|
Location = MatrixSymbol('Location', n, p)
|
||
|
|
Y = MatrixSymbol('Y', n, p)
|
||
|
|
M = MatrixStudentT('M', v, Location, Omega, Sigma)
|
||
|
|
|
||
|
|
exprd = gamma(v/2 + 1)*Determinant(Matrix([[1]]) + Sigma**(-1)*(-Location + Y)*Omega**(-1)*(-Location.T + Y.T))**(-v/2 - 1) / \
|
||
|
|
(pi*gamma(v/2)*sqrt(Determinant(Omega))*Determinant(Sigma))
|
||
|
|
|
||
|
|
assert density(M)(Y) == exprd
|
||
|
|
raises(ValueError, lambda: density(M)(1))
|
||
|
|
raises(ValueError, lambda: MatrixStudentT('M', 1, [1, 2], [[1, 0], [0, 1]], [[1, 0], [2, 1]]))
|
||
|
|
raises(ValueError, lambda: MatrixStudentT('M', 1, [1, 2], [[1, 0], [2, 1]], [[1, 0], [0, 1]]))
|
||
|
|
raises(ValueError, lambda: MatrixStudentT('M', 1, [1, 2], [[1, 0], [0, 1]], [[1, 0], [0, 1]]))
|
||
|
|
raises(ValueError, lambda: MatrixStudentT('M', 1, [1, 2], [[1, 0], [2]], [[1, 0], [0, 1]]))
|
||
|
|
raises(ValueError, lambda: MatrixStudentT('M', 1, [1, 2], [[1, 0], [2, 1]], [[1], [2]]))
|
||
|
|
raises(ValueError, lambda: MatrixStudentT('M', 1, [[1, 2]], [[1, 0], [0, 1]], [[1, 0]]))
|
||
|
|
raises(ValueError, lambda: MatrixStudentT('M', 1, [[1, 2]], [1], [[1, 0]]))
|
||
|
|
raises(ValueError, lambda: MatrixStudentT('M', -1, [1, 2], [[1, 0], [0, 1]], [4]))
|
||
|
|
|
||
|
|
def test_sample_scipy():
|
||
|
|
distribs_scipy = [
|
||
|
|
MatrixNormal('M', [[5, 6]], [4], [[2, 1], [1, 2]]),
|
||
|
|
Wishart('W', 5, [[1, 0], [0, 1]])
|
||
|
|
]
|
||
|
|
|
||
|
|
size = 5
|
||
|
|
scipy = import_module('scipy')
|
||
|
|
if not scipy:
|
||
|
|
skip('Scipy not installed. Abort tests for _sample_scipy.')
|
||
|
|
else:
|
||
|
|
for X in distribs_scipy:
|
||
|
|
samps = sample(X, size=size)
|
||
|
|
for sam in samps:
|
||
|
|
assert Matrix(sam) in X.pspace.distribution.set
|
||
|
|
M = MatrixGamma('M', 1, 2, [[1, 0], [0, 1]])
|
||
|
|
raises(NotImplementedError, lambda: sample(M, size=3))
|
||
|
|
|
||
|
|
def test_sample_pymc():
|
||
|
|
distribs_pymc = [
|
||
|
|
MatrixNormal('M', [[5, 6], [3, 4]], [[1, 0], [0, 1]], [[2, 1], [1, 2]]),
|
||
|
|
Wishart('W', 7, [[2, 1], [1, 2]])
|
||
|
|
]
|
||
|
|
size = 3
|
||
|
|
pymc = import_module('pymc')
|
||
|
|
if not pymc:
|
||
|
|
skip('PyMC is not installed. Abort tests for _sample_pymc.')
|
||
|
|
else:
|
||
|
|
for X in distribs_pymc:
|
||
|
|
samps = sample(X, size=size, library='pymc')
|
||
|
|
for sam in samps:
|
||
|
|
assert Matrix(sam) in X.pspace.distribution.set
|
||
|
|
M = MatrixGamma('M', 1, 2, [[1, 0], [0, 1]])
|
||
|
|
raises(NotImplementedError, lambda: sample(M, size=3))
|
||
|
|
|
||
|
|
def test_sample_seed():
|
||
|
|
X = MatrixNormal('M', [[5, 6], [3, 4]], [[1, 0], [0, 1]], [[2, 1], [1, 2]])
|
||
|
|
|
||
|
|
libraries = ['scipy', 'numpy', 'pymc']
|
||
|
|
for lib in libraries:
|
||
|
|
try:
|
||
|
|
imported_lib = import_module(lib)
|
||
|
|
if imported_lib:
|
||
|
|
s0, s1, s2 = [], [], []
|
||
|
|
s0 = sample(X, size=10, library=lib, seed=0)
|
||
|
|
s1 = sample(X, size=10, library=lib, seed=0)
|
||
|
|
s2 = sample(X, size=10, library=lib, seed=1)
|
||
|
|
for i in range(10):
|
||
|
|
assert (s0[i] == s1[i]).all()
|
||
|
|
assert (s1[i] != s2[i]).all()
|
||
|
|
|
||
|
|
except NotImplementedError:
|
||
|
|
continue
|