fehiepsi / rethinking-numpyro Goto Github PK

Chapter 4. Geocentric Models
4.4. Linear prediction
Code 4.48 (and Code 4.49 to display it)

Whenever N is equal or bigger than 54, sampled posterior returns rank-1 matrix (aka, all the rows are the same). Sounds like a numpyro package issue, but wanted to report it here since I cannot get to the bottom of it. Could you help?

In notebooks/05_the_many_variables_and_the_spurious_waffles.ipynb, there is code section Code 5.23 which I think should not be there. Am I missing something? If you want, I can make a PR that removes this and fix the following code block numbers

It will be nice if you had some quick tests to run that will validate that the examples run correctly with updates to NumPyro.

@fehiepsi Would you want any help in converting these notes into colab notebooks?

5.17 is repeated, the first occurrence should be 5.16.

Hi there,

I was trying to replicate the following Code 4.54 and I get the following:

Where the expected result is:

Please find a complete and inclusive minimum viable example below, which produces the same unexpected output:

In Chapter 2, code section 2.6 (2nd Edition), McElreath reports quadratic approximation values of

p 0.67   0.16 0.42  0.92

On my machine (numpyro 0.11.0, jax 0.4.10) I get

 p      0.63      0.13      0.43      0.83

Are these numerical differences to be expected? I do not have any experience with numpyro and not much with quadratic approximation, but I'm surprised at the divergence for this 1D problem. McElreath's solution is also more accurate to the true proportion.

I think it will be useful to pin the numpyro version in requirements.txt. As we update numpyro, we can check compatibility and bump up the number in this file if all the examples work correctly.

In code 11.6, if you use bw=0.3, it doesn't show the U-shaped curve of the prior predictive.
It is better to write just az.plot_kde(p).

• Use svi.run instead of init+update+lax.scan+get_params.

The code indexes from file 05_the_many_variables_and_the_spurious_waffles.ipynb, from 5.23 to 5.54, are wrong by +1.

The pull request would correct it.

Resolutions of the exercises could help those early learners following the book and this repo.

From the perspective of someone learning bayesian statistics and numpyro simultaneously, reading the book and following the examples in this repo is a bit challenging.

Mainly because the documentation of numpyro methods used in the code examples utilizes terms from bayesian statistics yet to be explained by the book. It seems unlikely that the first three chapter exercises could be done without much effort from the code examples without previous bayesian statistics knowledge, just based on the code examples. At least, that has been my personal experience following the book and the repo.

My intention is compile the solutions of all exercises of each chapter in numpyro in GitHub to help other readers.

What do you think about it? Would you already have the solutions in numpyro?

TLDR: I think a piece of code in the book is hacky in one place and works correctly by sheer coincidence; perhaps worth correcting for it here because it took quite a while for me to debug this (and probably took/will take for someone else) as i was following the example notebooks.

Proposed solution:
Turn this:

kl_data = dict(
    N=kl_dyads.shape[0],
    N_households=kl_dyads.hidB.max(),
    did=kl_dyads.did.values - 1,
    hidA=kl_dyads.hidA.values - 1,
    hidB=kl_dyads.hidB.values - 1,
    giftsAB=kl_dyads.giftsAB.values,
    giftsBA=kl_dyads.giftsBA.values,
)

Into this:

kl_data = dict(
    N=kl_dyads.shape[0],
    # N_households=kl_dyads.hidB.max(),
    did=kl_dyads.did.values - 1,
    hidA=kl_dyads.hidA.values - 1,
    hidB=kl_dyads.hidB.values - 1,
    giftsAB=kl_dyads.giftsAB.values,
    giftsBA=kl_dyads.giftsBA.values,
)

kl_data["N_households"] = len(set(kl_data["hidA"]) | set(kl_data["hidB"]))

Details:
While reproducing models from 2022 lectures and the book with pandas + numpyro + plotly stack and found a peculiar thing:

in 14.31 code N_households is computed as kl_dyads.hidB.max(), same as in the book. It's implied that hidB and hidA in original data are 1-indexed, so -1 is performed at kl_dyads init. But N_households is computed without -1.
And one must pay attention to it if one computes it not at init but somewhere downstream instead (like I insolently did). The model results (coefficients) with will look quite the same but the sampling will be bad (very low n_eff), as if not reparameterized. A recipe for a very fun time debugging correlated varying effects :)
The book's code, I think, is a hack since proper logic would imply instead computing a number of unique entries in hidB and hidA (see solution).

In Code4.36 of notebooks/04_geocentric_models.ipynb,

mu, sigma = jnp.mean(samples_flat, axis=1), jnp.cov(samples_flat) # axis =1 is missing