Wishart fails with valid arguments

[fonnesbeck]

Models that use the Wishart as a prior currently fail with valid parameters. For example,

s = pm.Wishart('s', 4, 3, np.eye(3))

Fails with the following:

TypeError: Wrong number of dimensions: expected 0, got 2 with shape (3, 3)

I've tried a variety of 3x3 matrices, but all yield similar errors.

PS - I have a pull request that removed the redundant second parameter.

[jsalvatier]

I think the issue is that you didn't specify 'shape'. That's definitely
something we should be able to infer thought.

On Fri, May 2, 2014 at 12:27 PM, Chris Fonnesbeck
notifications@github.comwrote:

Models that use the Wishart as a prior currently fail with valid
parameters. For example,

s = pm.Wishart('s', 4, 3, np.eye(3))

Fails with the following:

TypeError: Wrong number of dimensions: expected 0, got 2 with shape (3, 3)

I've tried a variety of 3x3 matrices, but all yield similar errors.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/538
.

[fonnesbeck]

Yeah, I tried that, and got:

---------------------------------------------------------------------------
AsTensorError                             Traceback (most recent call last)
<ipython-input-73-6f02cca3962c> in <module>()
      1 with pm.Model():
----> 2     s = pm.Wishart('s', 4, 3, np.eye(3), shape=(3,3))

/Users/fonnescj/Code/pymc/pymc/distributions/distribution.pyc in __new__(cls, name, *args, **kwargs)
     17             data = kwargs.pop('observed', None)
     18             dist = cls.dist(*args, **kwargs)
---> 19             return model.Var(name, dist, data)
     20         elif name is None:
     21             return object.__new__(cls) #for pickle

/Users/fonnescj/Code/pymc/pymc/model.pyc in Var(self, name, dist, data)
    142         """
    143         if data is None:
--> 144             var = FreeRV(name=name, distribution=dist, model=self)
    145             self.free_RVs.append(var)
    146         else:

/Users/fonnescj/Code/pymc/pymc/model.pyc in __init__(self, type, owner, index, name, distribution, model)
    321             self.tag.test_value = np.ones(
    322                 distribution.shape, distribution.dtype) * distribution.default()
--> 323             self.logp_elemwiset = distribution.logp(self)
    324             self.model = model
    325 

/Users/fonnescj/Code/pymc/pymc/distributions/multivariate.pyc in logp(self, X)
    179              2) - n * log(IVI) - 2 * multigammaln(p, n / 2)) / 2,
    180 
--> 181             all(n > p - 1))

/Library/Python/2.7/site-packages/theano/tensor/basic.pyc in all(x, axis, keepdims)
   5078 
   5079 def all(x, axis=None, keepdims=False):
-> 5080     out = elemwise.All(axis)(x)
   5081 
   5082     if keepdims:

/Library/Python/2.7/site-packages/theano/gof/op.pyc in __call__(self, *inputs, **kwargs)
    409         """
    410         return_list = kwargs.pop('return_list', False)
--> 411         node = self.make_node(*inputs, **kwargs)
    412         if self.add_stack_trace_on_call:
    413             self.add_tag_trace(node)

/Library/Python/2.7/site-packages/theano/tensor/elemwise.pyc in make_node(self, input)
   1639 
   1640     def make_node(self, input):
-> 1641         input = as_tensor_variable(input)
   1642         if input.dtype not in ["int8", "uint8"]:
   1643             input = theano.tensor.neq(input, 0)

/Library/Python/2.7/site-packages/theano/tensor/basic.pyc in as_tensor_variable(x, name, ndim)
    181     if isinstance(x, bool):
    182         raise AsTensorError(
--> 183             "Cannot cast True or False as a tensor variable. Please use 1 or "
    184             "0. This error might be caused by using the == operator on "
    185             "Variables. v == w does not do what you think it does, "

AsTensorError: Cannot cast True or False as a tensor variable. Please use 1 or 0. This error might be caused by using the == operator on Variables. v == w does not do what you think it does, use theano.tensor.eq(v, w) instead.

It does not seem to like the boundary check.

[jsalvatier]

Hmm that seems bad. Do we have tests for it?
On May 2, 2014 6:19 PM, "Chris Fonnesbeck" notifications@github.com wrote:

Yeah, I tried that, and got:

---

AsTensorError Traceback (most recent call last)
in ()
1 with pm.Model():
----> 2 s = pm.Wishart('s', 4, 3, np.eye(3), shape=(3,3))

/Users/fonnescj/Code/pymc/pymc/distributions/distribution.pyc in new(cls, name, _args, *_kwargs)
17 data = kwargs.pop('observed', None)
18 dist = cls.dist(_args, *_kwargs)
---> 19 return model.Var(name, dist, data)
20 elif name is None:
21 return object.new(cls) #for pickle

/Users/fonnescj/Code/pymc/pymc/model.pyc in Var(self, name, dist, data)
142 """
143 if data is None:
--> 144 var = FreeRV(name=name, distribution=dist, model=self)
145 self.free_RVs.append(var)
146 else:

/Users/fonnescj/Code/pymc/pymc/model.pyc in init(self, type, owner, index, name, distribution, model)
321 self.tag.test_value = np.ones(
322 distribution.shape, distribution.dtype) * distribution.default()
--> 323 self.logp_elemwiset = distribution.logp(self)
324 self.model = model
325

/Users/fonnescj/Code/pymc/pymc/distributions/multivariate.pyc in logp(self, X)
179 2) - n * log(IVI) - 2 * multigammaln(p, n / 2)) / 2,
180
--> 181 all(n > p - 1))

/Library/Python/2.7/site-packages/theano/tensor/basic.pyc in all(x, axis, keepdims)
5078
5079 def all(x, axis=None, keepdims=False):
-> 5080 out = elemwise.All(axis)(x)
5081
5082 if keepdims:

/Library/Python/2.7/site-packages/theano/gof/op.pyc in call(self, _inputs, *_kwargs)
409 """
410 return_list = kwargs.pop('return_list', False)
--> 411 node = self.make_node(_inputs, *_kwargs)
412 if self.add_stack_trace_on_call:
413 self.add_tag_trace(node)

/Library/Python/2.7/site-packages/theano/tensor/elemwise.pyc in make_node(self, input)
1639
1640 def make_node(self, input):
-> 1641 input = as_tensor_variable(input)
1642 if input.dtype not in ["int8", "uint8"]:
1643 input = theano.tensor.neq(input, 0)

/Library/Python/2.7/site-packages/theano/tensor/basic.pyc in as_tensor_variable(x, name, ndim)
181 if isinstance(x, bool):
182 raise AsTensorError(
--> 183 "Cannot cast True or False as a tensor variable. Please use 1 or "
184 "0. This error might be caused by using the == operator on "
185 "Variables. v == w does not do what you think it does, "

AsTensorError: Cannot cast True or False as a tensor variable. Please use 1 or 0. This error might be caused by using the == operator on Variables. v == w does not do what you think it does, use theano.tensor.eq(v, w) instead.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/538#issuecomment-42092620
.

[fonnesbeck]

We have one, which it passes. Does my example fail for you?

[jsalvatier]

Yup, it fails for me too.

On Fri, May 2, 2014 at 7:56 PM, Chris Fonnesbeck
notifications@github.comwrote:

We have one, which it passes. Does my example fail for you?

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/538#issuecomment-42094468
.

[twiecki]

#547 has a bunch of changes to Wishart. Might that fix it? @kadeng?

[kadeng]

@twiecki: You should definitely take a look at the changes I made there. I think that the current Wishart Implementation does not really work. I found discrepancies between the likelihood formula in the code and the formulas I found in literature (for example, in Gelman's Bayesian Data Analysis and in Murphy's Machine Learning, a probabilistic perspective). What's problematic there is that different authors seem to use different parametrizations and seem to name the parameters inconsistently, so in the end I was pretty confused.

But it's not like I have spent a lot of time to validate that my fix in turn is correct and works, so please don't rely on it. Use it as a suggestion :) The following paper is also an interesting read when it comes to actually using the Wishart or inverse Wishart as prior for the precision or covariance matrix for MVNs: http://ba.stat.cmu.edu/journal/2013/vol08/issue02/huang.pdf

[wgreenr]

Recently have tried to sample from Wishart prior using:

import pymc3 as pm
import numpy as np

prec_prior = np.array([[ 25.3968254,   -1.58730159],
                       [ -1.58730159,   6.34920635]])
with pm.Model() as model:
    prec = pm.Wishart('prec', 100.0, prec_prior / 100.0, shape=(2, 2))
    step = pm.Metropolis()
    trace = pm.sample(10000, step)

And instead I obtained clearly diverging sampling process, with some entries of sampled matrices moving towards negative infinity, and not in any way resembling the prior precision matrix:

Then I decided to check whether my notions about Wishart distribution is fundamentally flawed, and using the code from http://www.mit.edu/~mattjj/released-code/hsmm/stats_util.py,
I obtained samples which look as one might expect, given a quite narrow prior centered at precision matrix prec_prior:

    sample_wishart(prec_prior / 100.0, 100.0)
    array([[ 26.37928925,  -3.88792305],
           [ -3.88792305,   7.6372514 ]])

The question is: does this behavior of MCMC sampler is something expected, or is this a bug, and how to properly sample from Wishart?

[kiudee]

Yes, this is a bug and the same behavior I observed with the Wishart distribution (see issue #672).
If you can find the source of the error in the current code, we would greatly appreciate it.

I recommend taking a look at the LKJ distribution, which is a prior distribution for the correlation matrix instead of the covariance/precision matrix.
This has the advantage that you can specify separate priors for the standard deviations and the correlations.

[wgreenr]

Will run some tests, Wishart logp seems to look good

[twiecki]

Seems like we're not enforcing the tau > 0 constraint.

[kiudee]

Indeed, this is the case. (See below)

I added the constraint and tried it with the given example:

import numpy as np
import pymc3 as pm
import scipy.optimize as opt
prec_prior = np.array([[ 25.3968254,   -1.58730159],
                       [ -1.58730159,   6.34920635]])
with pm.Model() as model:
    prec = pm.Wishart('prec', 100.0, prec_prior / 100.0, shape=(2, 2))

    start = pm.find_MAP(fmin=opt.fmin_powell)
    step = pm.Metropolis()
    trace = pm.sample(10000, step, start)

Here is the resulting trace (without any thinning):

PS: I’m sorry, it is too late here. I just realized that this is not true. The matrices need to be positive definite.

[twiecki]

Oh, right. I experimented a bit with all(gt(eig(X), 0)) as a condition but didn't get it to work yet.