fixed posterior covariance calculation of Code 13.15, closes #53

[chdamianos]

Closes issue pymc-devs/pymc#53 by modifying the calculation of L_chol in section Code 13.15 from L_chol[np.triu_indices(2, 0)] = Chol_cov to L_chol[np.tril_indices(2, 0)] = Chol_cov. This modification corrects L_chol to a lower triangular matrix of the Cholesky decomposition and not an upper triangular matrix.

[review-notebook-app]

Check out this pull request on ReviewNB: https://app.reviewnb.com/pymc-devs/resources/pull/54

You'll be able to see visual diffs and write comments on notebook cells. Powered by ReviewNB.

[AlexAndorra]

Thank you, that's well spotted!
And I agree with you (#53 ): maybe we should seize on this PR to change the way we compute back the covariance matrix from the trace.
Your proposal with pm.expand_packed_triangular is elegant and PyMC3onic.
Your second proposal (doing as in the book) is more pedagogical though, as it shows the details of the computation.
So I don't have a set preference here on what to do. What do you think @aloctavodia ?

[chdamianos]

@AlexAndorra I agree that pm.expand_packed_triangular is easier to read, also taking into account that this was the way the model was defined. On the other hand, as you mentioned, in it's current implementation it makes the reader think a bit more on how the Cholesky decomposition is used to calculate the covariance. Also using np.tril_indices is more consistent with the rest of code. I'm leaning towards np.tril_indices for consistency.

[AlexAndorra]

Yeah I get your point. However, if we're not going with pm.expand_packed_triangular I think doing the computation like in the book would be better than sticking with the current implentation: I find it clearer to the user, more pedagogical and consistent with the spirit of this whole repo (staying as close to the book as possible).

[chdamianos]

OK that makes sense. I've changed the calculation to make use of pm.expand_packed_triangular. I also supressed FutureWarning because I thought it was making the notebook a bit difficult to read.

As a first time reader I think it would benefit me to have seen one or two examples of using the pm.LKJCorr prior, instead of pm.LKJCholeskyCov for the covariance, which would be closer to the code examples in the book. After I finish the book I might have a go adding an example with pm.LKJCorr.

[AlexAndorra]

Thank you, LGTM now!

I don't know exactly but I think pm.LKJCholeskyCovwas prefered to pm.LKJCorr for computation efficiency reasons. That being said, I think you're right, one or two examples of how to use the correlation prior could be useful to readers. I can help you with that later if you need :)

And about the FutureWarning, do you know where it comes from? I have the same annoying message but I'd like to get rid of it entirely instead of just hiding it.

[chdamianos]

I'm under the same impression about using pm.LKJCholeskyCov. Numerically it make sense to be more efficient to work with an upper triangular matrix (L) than the whole covariance matrix. Thanks for offering help! If I try to get some code examples working with pm.LKJCorr I'll be very happy to collaborate.

The warnings are a mystery for me at the moment but I think are the same as in this Theano issue discussed here Theano/Theano#6695

[lucianopaz]

I just wanted to chime in about the future warnings. I managed to track them down in this issue, but I didn't see an immediate fix at the time.

[AlexAndorra]

Oh ok, thanks Luciano, that's good to know! So looks like the only solution for the moment is indeed to hide them in the notebooks.

[AlexAndorra]

FYI @chdamianos pm.LKJCorr seems faulty and may be deprecated - while pm.LKJCholeskyCov is right

[chdamianos]

thanks @AlexAndorra for the heads up! Yes it looks like pm.LKJCholeskyCov is the safest prior for covariance matrices.

[AlexAndorra]

Does anyone with write access see why there is a conflict with this PR?
@chdamianos are you sure you based this on the latest master?
As wrote above, the change in the posterior covariance calculation looks good to me :)

[chdamianos]

@AlexAndorra thanks for the heads up. It looks like the notebook has been updated since my pull request. I will pull latest master and update

[chdamianos]

@AlexAndorra I rebased to latest master. It should be ok now.

[junpenglao]

Could you also please rerun the notebook from the latest PyMC3 release?

[chdamianos]

@junpenglao good catch, it's now run with the latest release

[junpenglao]

Thanks!