BUG Fix bounds of Wishart and add warning.

The Wishart has never worked correctly and produced
divergent samples. One of the reasons was that we did
not enforce the matrix to be symmetric nor positive
semi-definite.

Now, however, it is impossible to randomly sample a
valid matrix so we raise a warning that the Wishart
is basically unusable currently and to instead
use the LKJ prior which has many desriable properties.

For more discussions, see #538.

Author: twiecki

pymc3/distributions/multivariate.py

@@ -1,9 +1,11 @@
+import warnings
+
 from .dist_math import *
 
 import numpy as np
 
-from theano.tensor.nlinalg import det, matrix_inverse, trace
-from theano.tensor import dot, cast, eye, diag, eq, le, ge, all
+from theano.tensor.nlinalg import det, matrix_inverse, trace, eigh
+from theano.tensor import dot, cast, eye, diag, eq, le, ge, gt, all
 from theano.printing import Print
 
 __all__ = ['MvNormal', 'Dirichlet', 'Multinomial', 'Wishart', 'LKJCorr']
@@ -160,6 +162,7 @@ class Wishart(Continuous):
     """
     def __init__(self, n, V, *args, **kwargs):
         super(Wishart, self).__init__(*args, **kwargs)
+        warnings.warn('The Wishart distribution can currently not be used for MCMC sampling. The probability of sampling a symmetric matrix is basically zero. Instead, please use the LKJCorr prior. For more information on the issues surrounding the Wishart see here: https://github.com/pymc-devs/pymc3/issues/538.', UserWarning)
         self.n = n
         self.p = p = V.shape[0]
         self.V = V
@@ -179,7 +182,10 @@ def logp(self, X):
         return bound(
             ((n - p - 1) * log(IXI) - trace(matrix_inverse(V).dot(X)) -
                 n * p * log(2) - n * log(IVI) - 2 * multigammaln(n / 2., p)) / 2,
-             n > (p - 1))
+            gt(n, (p - 1)),
+            all(gt(eigh(X)[0], 0)),
+            eq(X, X.T)
+        )
 
 
 class LKJCorr(Continuous):
@@ -249,7 +255,7 @@ def logp(self, x):
 
         X = x[self.tri_index]
         X = t.fill_diagonal(X, 1)
-        
+
         result = self._normalizing_constant(n, p)
         result += (n - 1.) * log(det(X))
         return bound(result,