Add Wishart with Bartlett decomposition. #701
Conversation
|
I think this looks pretty good, but I think long term it might be nicer to maybe have something like that does this? Maybe that doesn't really work here. Also an example of usage would be nice. |
|
Right, in this case it would be pretty difficult to do that. What do you think of adding new RVs under the hood? Can't really think of a better way unfortunately. I'll add an example. |
|
Depends on #731 |
twiecki
changed the title
|
Hey with the new automatic transforms I think this should be possible. |
|
I'm not sure but I think this is not possible here. It's not really a transform of the Wishart but rather a reformulation. The Wishart is actually never used explicitly. Or do you have another idea? |
|
Hmm, yes I see that now. We could leave it as is and I think that would be fine. We can also create a transform that does the bartlett decomposition (and the reverse), X = L * A * A.T *L.T. That's a transform from unconstrained spaces to the positive definite space, right? So that should work as a transform for using the Wishart. |
|
But there's no function I know of to go from X to L and A. |
|
Ah, interesting. Right now we require it, but I wonder if that's really necessary. I think probably the only thing we use it for is to move the default point, but we could probably not require that. The other thing that would be needed is the log determinant of the jacobian, which might be tricky, but I don't know. |
|
Well we definitely need the direction from Wishart -> transformed Wishart, and that's the missing step. Also, the transformation only works by adding two other RVs so I don't see a way to fit this into the current transform framework (although it would be great if it were possible). |
|
I don't see why we need the direction Wishart -> transformed Wishart. It mostly seems like we need the direction inputs -> wishart -> rest of the model, which is transformed Wishart -> Wishart -> rest of the model. Except for getting a good starting value from the wishart distribution, but again, I think we can probably be clever about that. |
|
I might be confused about transforms, help me clear that up: But I'm not sure if we even have that. Bartlett starts with separate RVs that it combines that just happen to be distributed the same as a Wishart, but there is no function I know of to go from one to the other. |
|
Maybe I'm missing something. If it's easy you should be able to provide me with that function and I'd be happy to include it, but I really don't know how to do that. |
def transform(x):
diag_idx = np.diag_indices_from(S)
tril_idx = np.tril_indices_from(S, k=-1)
n_diag = len(diag_idx[0])
n_tril = len(tril_idx[0])
c = x[:n_diag]
z = x[n_diag:]
# Construct A matrix
A = zeros(S.shape, dtype=np.float32)
A = set_subtensor(A[diag_idx], c)
A = set_subtensor(A[tril_idx], z)
# L * A * A.T * L.T ~ Wishart(L*L.T, nu)
return dot(dot(dot(L, A), A.T), L.T)Of course there are a few problems. One, you have to make x the size of n_diag before its passed in, which we currently do by doing a backwards transform, but that's not an option here. You also clearly need to pass in L somehow, which is not currently supported. You also need to compute the jacobian determinant, which might be difficult. I haven't attempted that. These things seem fixable. Is there something I'm missing? |
|
Another problem is that it's not enforcing x and S to have the right distributions: |
|
Or do you assume that those would be specified in the model itself? |
|
If I understand things correctly (and I may not), we're just using the pdf of the Wishart instead of the Chi and Normal densities, so we don't need to do anything extra (except the jacobian det). |
|
Well, the key is that we don't need to deal with the wishart which is impossible to sample from. Instead we are sampling from a chi and normal and combine them in a way so that the result is distributed like a Wishart (but the Wishart pdf never enters into the equation here). Seems like you are talking about a different transform -- what is the input x for your transform? A wishart pdf? |
|
x is a pure theano variable with no likelihood attached. The issue is that we can't sample from the Wishart because its too constrained. All the values we propose are illegal. So we reformulate in a lower dimensional space (with roughly half the elements), x is in that lower dimensional space. Then we transform back into the original space to compute the density (which now needs a log jacobian determinant adjustment). Basically the thing I'm proposing is to compute the density in the original wishart space rather than the new Chi/Normal space. It may turn out that computing the density is much more efficient in the Chi/Normal space, in which case I think your original proposal is the roughly right (I actually have an alternative, but we'll leave that for later). |
|
Ah, I had assumed you wanted to do the Bartlett decomposition but you're proposing something else entirely. Instead, you just sample the upper triangular and then create a positive symmetric matrix from there, correct? I think we do something similar for the LKJ prior. |
|
Yeah, exactly. On Mon, Jun 8, 2015 at 1:07 AM, Thomas Wiecki [email protected]
|
|
That might help with the symmetry but I don't think it enforces positive definitness, no? |
|
I think it should since in the bartlett decomposition the off diagonal
|
|
Is there a status update on this, or should we close it? |
|
I do plan to get back to this. It works fine except for the travis issue. |
|
Hey @twiecki do you need some help getting this into Master? Anything I can do in my period of unemployment? |
|
@springcoil Thanks for offering. The last issue are the failing unittests on travis I wasn't able to pin down yet (and can't reproduce locally). |
|
Just rebased this PR off of master. |
|
|
Does this work now with the recent updates of Theano etc? I haven't had time to try it. |
|
I haven't tried. It should. |
|
Ok I pushed this through - after doing some rebasing - it should work. |
DOC Add encoding information. Add WishartBartlett example. ENH Add run function to wishart example. BUG Wishart example did not use proper namespace. ENH Fix imports. Add verbosity arguments. ENH: Fixing up imports ENH: Fixing up multivariate imports
|
Since the tests pass I'll merge this. |
|
Awesome! |
Closes #538.
This adds the Wishart distribution that can be sampled from because it's using matrix decomposition to sample from an unconstrained space.
However, it's a bit odd as it adds RVs itself to the model. I think it's pretty convenient though. @jsalvatier @fonnesbeck thoughts on adding this?