Update case study: LKJ

[chiral-carbon]

Addresses issue #53 and aims to advance it to best practices.

I did not see much to update here apart from using numpy.random.default_rng() so let me know if any improvements could be made in the code or plotting.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

chiral-carbon commented on 2021-08-09T17:11:22Z
----------------------------------------------------------------

get too many warnings here but not sure how to address these.

OriolAbril commented on 2021-08-10T07:28:20Z
----------------------------------------------------------------

see which variables are the culprits of the warnings. I think only the diagonal of chol_corr which is constant should raise these warnings, and it should, as a constant the kde doesn't work. I think that with chol_corr being a symmetric matrix of size 2 with ones on the diagonal, only one of the off diagonal elements should be plotted.

Once using coords via idata_kwargs this can be done here and should at least reduce the warnings. It would also be good to link to https://arviz-devs.github.io/arviz/user_guide/label_guide.html#custom-labellers because for >2d matrices pointwise selection is needed.

chiral-carbon commented on 2021-08-10T22:43:38Z
----------------------------------------------------------------

sure. will try that

[OriolAbril]

see which variables are the culprits of the warnings. I think only the diagonal of chol_corr which is constant should raise these warnings, and it should, as a constant the kde doesn't work. I think that with chol_corr being a symmetric matrix of size 2 with ones on the diagonal, only one of the off diagonal elements should be plotted.

Once using coords via idata_kwargs this can be done here and should at least reduce the warnings. It would also be good to link to https://arviz-devs.github.io/arviz/user_guide/label_guide.html#custom-labellers because for >2d matrices pointwise selection is needed.

---

View entire conversation on ReviewNB

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-08-10T07:31:45Z
----------------------------------------------------------------

Line #2.        trace = pm.sample(random_seed=RANDOM_SEED, init="adapt_diag", return_inferencedata=True)

It would be good to show how one can use named dims and coords for cholesky matrices. There is one example in model covariation_intercept_slope at https://docs.pymc.io/notebooks/multilevel_modeling.html.

Line #3.    az.summary(trace, var_names=["~chol"], round_to=2)

Here I am not sure this ~chol is doing anything. Is there a variable named exactly chol? Same in plot_trace below.

chiral-carbon commented on 2021-08-12T16:22:25Z
----------------------------------------------------------------

yup, there is a variable named chol. if you look where model is defined 2 cells before:

chol, corr, stds = pm.LKJCholeskyCov(
       "chol", n=2, eta=2.0, sd_dist=pm.Exponential.dist(1.0), compute_corr=True
   )

OriolAbril commented on 2021-08-12T16:36:25Z
----------------------------------------------------------------

You have to look this in the InferenceData, not in the model. I am sure this line creates the chol_corr and chol_stds variables, not completely sure it creates a chol alone one, hence my asking.

chiral-carbon commented on 2021-08-12T19:15:46Z
----------------------------------------------------------------

oh yes, that it does. I have printed the entire data here.

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-08-10T07:31:46Z
----------------------------------------------------------------

As personal preference, here I would use non filled ellipses so it's easier to see the differences. I would also use some other colors because those two are hard to distinguish for people with tritanomaly.

chiral-carbon commented on 2021-08-12T19:17:38Z
----------------------------------------------------------------

I removed fill color and wanted to use colors that were contrasting enough but red and yellow, for example, kind of blended together. used black and yellow in the end, these are also not very clear but better than before. if you have a better suggestion though then let me know.

[OriolAbril]

I did not see much to update here apart from using numpy.random.default_rng() so let me know if any improvements could be made in the code or plotting.

also fix the comment on the issue description about using .values instead of .data

[chiral-carbon]

sure. will try that

---

View entire conversation on ReviewNB

[chiral-carbon]

yup, there is a variable named chol. if you look where model is defined 2 cells before:

chol, corr, stds = pm.LKJCholeskyCov(
       "chol", n=2, eta=2.0, sd_dist=pm.Exponential.dist(1.0), compute_corr=True
   )

---

View entire conversation on ReviewNB

[OriolAbril]

You have to look this in the InferenceData, not in the model. I am sure this line creates the chol_corr and chol_stds variables, not completely sure it creates a chol alone one, hence my asking.

---

View entire conversation on ReviewNB

[chiral-carbon]

oh yes, that it does. I have printed the entire data here.

---

View entire conversation on ReviewNB

[chiral-carbon]

I removed fill color and wanted to use colors that were contrasting enough but red and yellow, for example, kind of blended together. used black and yellow in the end, these are also not very clear but better than before. if you have a better suggestion though then let me know.

---

View entire conversation on ReviewNB

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

MarcoGorelli commented on 2021-08-13T10:35:03Z
----------------------------------------------------------------

 no divergences and good r-hats

One of the r-hats is NaN - do you know why?

OriolAbril commented on 2021-08-13T12:37:31Z
----------------------------------------------------------------

the diagonal of the correlation matrix is not a random variable, it's always 1. As rhat compares the distributions from different chains to see if they are compatible with being the same with terms like within chain variance and between chain variance, both 0 for constants, and some are on the denominator rhat should be NaN as it doesn't really make sense to use the diagnostic on constants. What I don't understand is why chol_corr[1,1] is not NaN too

MarcoGorelli commented on 2021-08-15T16:40:36Z
----------------------------------------------------------------

What I don't understand is why chol_corr[1,1] is not NaN too

I'm seeing the same in NumPyro - looks like, for that element, there's some small numerical errors which make it not exactly 1 everywhere:

ipdb>  x.shape
(4, 1000, 1, 2, 2)
ipdb>  x[:, 0, :, 0, 0]
array([[1.],

[1.],

[1.],

[1.]], dtype=float32)
ipdb>  x[:, 0, :, 1, 1]
array([[1.0000001 ],

[1.        ],

[0.99999994],

[0.99999994]], dtype=float32)

Perhaps a quick explanation would make the " no divergences and good r-hats" clearer? something like:

 no divergences and good r-hats (except for the diagonal elements of the correlation matrix - however, these are not a concern, because, they should be equal 1 for each sample for each chain and the variance of a constant value isn't defined. If one of the diagonal elements has r_hat defined, it's likely due to tiny numerical errors)

chiral-carbon commented on 2021-08-17T13:09:51Z
----------------------------------------------------------------

@MarcoGorelli should I add this explanation in notebook?

MarcoGorelli commented on 2021-08-17T13:32:59Z
----------------------------------------------------------------

Hey @chiral-carbon - yeah, I think so - just here in the markdown cell would be fine, no need to re-run the whole notebook

[OriolAbril]

the diagonal of the correlation matrix is not a random variable, it's always 1. As rhat compares the distributions from different chains to see if they are compatible with being the same with terms like within chain variance and between chain variance, both 0 for constants, and some are on the denominator rhat should be NaN as it doesn't really make sense to use the diagnostic on constants. What I don't understand is why chol_corr[1,1] is not NaN too

---

View entire conversation on ReviewNB

[MarcoGorelli]

What I don't understand is why chol_corr[1,1] is not NaN too

I'm seeing the same in NumPyro - looks like, for that element, there's some small numerical errors which make it not exactly 1 everywhere:

ipdb>  x.shape
(4, 1000, 1, 2, 2)
ipdb>  x[:, 0, :, 0, 0]
array([[1.],

[1.],

[1.],

[1.]], dtype=float32)
ipdb>  x[:, 0, :, 1, 1]
array([[1.0000001 ],

[1.        ],

[0.99999994],

[0.99999994]], dtype=float32)

Perhaps a quick explanation would make the " no divergences and good r-hats" clearer? something like:

 no divergences and good r-hats (except for the diagonal elements of the correlation matrix - however, these are not a concern, because, they should be equal 1 for each sample for each chain and the variance of a constant value isn't defined. If one of the diagonal elements has r_hat defined, it's likely due to tiny numerical errors)

---

View entire conversation on ReviewNB

[chiral-carbon]

@MarcoGorelli should I add this explanation in notebook?

---

View entire conversation on ReviewNB

[MarcoGorelli]

Hey @chiral-carbon - yeah, I think so - just here in the markdown cell would be fine, no need to re-run the whole notebook

---

View entire conversation on ReviewNB

[MarcoGorelli]

looks good to me (@OriolAbril not sure how actionable the comment about all the warnings is at the moment or whether that's OK to leave as a follow-up, will leave this with you 😄 )

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-08-24T17:12:38Z
----------------------------------------------------------------

Line #29.    rect = plt.Rectangle((0, 0), 1, 1, fc="C0", alpha=0.5)

This should be updated to use Line instead. The Rectangle class in this line is only being used to generate the "example" in the legend. now with the non-filled ellipse, lines match better the plot. I would also not use any transparency now with the non-filled ellipses, maybe also make the lines a bit thicker and/or make one dashed

[OriolAbril]

We'll go over the two minor nits still missing in my opinion and merge in the next few days. Trying to make the final plot a bit more clear and see if we can remove some (or even all) warnings in the plot_trace figure

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-08-26T15:26:49Z
----------------------------------------------------------------

just seeing this now, the mean should use named dimension names. same in the cell below. It should also use .values after instead of .data

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-08-26T15:26:50Z
----------------------------------------------------------------

Line #14.    e.set_alpha(0.5)

I would not make the lines transparent here.

And I would make the dashed line be above the continuous line. This can be modified with the zorder

Line #16.    e.set_zorder(10)