Symmetry, non-negative-definiteness not enforced for Wishart sampling

[jangevaare]

Each sample generated from the Wishart distribution should be symmetric and non-negative definite, but that doesn't seem to be enforced in pymc.

e.g.

import pymc as pm
import numpy as np

cov=np.array([[2,1],[1,3]])
mean=([2,7])
tau=np.linalg.inv(cov)
N=1
z_data=np.ndarray.flatten(np.random.multivariate_normal(mean, cov, N))

with pm.Model() as model:
    wishart_test=pm.Wishart('wishart_test', n=2, p=2, V=np.eye(2), shape=(2,2))
    z=pm.MvNormal('z', mu=mean, tau=wishart_test, shape=2, observed=z_data)

with model:
    start = {'wishart_test':tau}

with model:
    step = pm.Metropolis()    

with model: 
    trace = pm.sample(3000, step, start)

trace[wishart_test][500]

np.linalg.eig(trace[wishart_test][500])[0]

For this one test, the 500th iteration was neither symmetric nor non-negative definite.

Out[1]: 
array([[-119.98956617,  -93.5080308 ],
       [-137.06886717, -119.48174073]])

Out[2]: 
array([-232.94830435,   -6.52300255])

[jsalvatier]

You're right, that needs to be added.

On Wed, Nov 20, 2013 at 10:08 AM, Justin Angevaare <[email protected]

wrote:

Each sample generated from the Wishart distribution should be symmetric
and non-negative definite, but that doesn't seem to be enforced in pymc.

e.g.

import pymc as pm
import numpy as np

cov=np.array([[2,1],[1,3]])
mean=([2,7])
tau=np.linalg.inv(cov)
N=1
z_data=np.ndarray.flatten(np.random.multivariate_normal(mean, cov, N))

with pm.Model() as model:
wishart_test=pm.Wishart('wishart_test', n=2, p=2, V=np.eye(2), shape=(2,2))
z=pm.MvNormal('z', mu=mean, tau=wishart_test, shape=2, observed=z_data)

with model:
start = {'wishart_test':tau}

with model:
step = pm.Metropolis()

with model:
trace = pm.sample(3000, step, start)

trace[wishart_test][500]

np.linalg.eig(trace[wishart_test][500])[0]

For this one test, the 500th iteration was neither symmetric nor
non-negative definite.

Out[1]:
array([[-119.98956617, -93.5080308 ],
[-137.06886717, -119.48174073]])

Out[2]:
array([-232.94830435, -6.52300255])

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

[jangevaare]

Is there a more efficient way of doing this other than adding conditions in bound(logp, *conditions) for the Wishart?

[fonnesbeck]

The most efficient way is to ensure that only non-negative definite matrices are proposed in the first place. You can do this by constructing matrices in your model with some constraints. Symmetry is easy, for example.

[jangevaare]

I see. So getting Bartlett Decomposition (http://en.wikipedia.org/wiki/Wishart_distribution#Bartlett_decomposition) to work would probably be the best way? This is what I was going after when I ran into #379

[twiecki]

I added the bounds now but it renders the Wishart unusable: 642f639. Still going to close this as we'll better track this issue over here: #538