Example for Hilbert space approximation of Gaussian processes

[omarfsosa]

This example replicates a few of the models in the excellent case
study by Aki Vehtari [1] (originally written using R and Stan).
The case study uses approximate Gaussian processes [2] to model the
relative number of births per day in the US from 1969 to 1988.
The Hilbert space approximation is way faster than the exact Gaussian
processes because it circumvents the need for inverting the
covariance matrix.

The original case study presented by Aki also emphasizes the iterative
process of building a Bayesian model, which is excellent as a pedagogical
resource. Here, however, I replicate only 4 out of all the models available in [1].
There are a few minor differences in the mathematical details of my model GP4,
which I had to make in order for the chains to mix properly. I have clearly
commented on the places where our models are different.

References:
1. Gelman, Vehtari, Simpson, et al (2020), "Bayesian workflow book - Birthdays" <https://avehtari.github.io/casestudies/Birthdays/birthdays.html>.
2. Riutort-Mayol G, Bürkner PC, Andersen MR, et al (2020),
"Practical Hilbert space approximate bayesian gaussian processes for probabilistic programming"

[fehiepsi]

you can also use a wrapper numpyro.enable_x64()

[omarfsosa]

Nice, thanks for this tip. I've updated the PR.

[fehiepsi]

This is superb, @omarfsosa! I really like your code structure and I also think that you've incorporated the best practices of using (Num)Pyro for modeling. Please let me know if you need me to take a closer look at something.

Could you do a few steps:

• run make license
• if you prefer using I in the docstring, do you want to add your name too? You can put "Author" line at any place in the docstring.
• add your script to rtd index

[omarfsosa]

This is superb, @omarfsosa! I really like your code structure and I also think that you've incorporated the best practices of using (Num)Pyro for modeling. Please let me know if you need me to take a closer look at something.

Could you do a few steps:

• run make license
• if you prefer using I in the docstring, do you want to add your name too? You can put "Author" line at any place in the docstring.
• add your script to rtd index

Thanks for the feedback and suggestions @fehiepsi ! I've made some changes according to your recommendations. If you have a better idea for how to get load the data I'd love to hear it. At the moment I'm pulling the csv directly from the repo of the original case study, which feels a bit fragile maybe?

[fehiepsi]

I guess it is fine to use the link in your script. Currently, we have some load_data utility in datasets that you can mimic.

feels a bit fragile

How about adding a test? :)

Btw, you'll need to run make format to pass CI.

[fehiepsi]

@omarfsosa One of the test is failing

Traceback (most recent call last):
  File "/home/runner/work/numpyro/numpyro/examples/hsgp.py", line 717, in <module>
    main(args)
  File "/home/runner/work/numpyro/numpyro/examples/hsgp.py", line 686, in main
    data = model.get_data()
  File "/home/runner/work/numpyro/numpyro/examples/hsgp.py", line 596, in get_data
    **self._get_floating_days(data),
  File "/home/runner/work/numpyro/numpyro/examples/hsgp.py", line 541, in _get_floating_days
    & data["date"].dt.day.ge(25),
AttributeError: 'DatetimeProperties' object has no attribute 'day_of_week'
test/test_examples.py::test_cpu[hsgp.py --model tywd --num-samples 10 --num-warmup 10 --num-chains 2] Running:
python examples/hsgp.py --model tywd --num-samples 10 --num-warmup 10 --num-chains 2
FAILED

[omarfsosa]

@omarfsosa One of the test is failing

Traceback (most recent call last):
  File "/home/runner/work/numpyro/numpyro/examples/hsgp.py", line 717, in <module>
    main(args)
  File "/home/runner/work/numpyro/numpyro/examples/hsgp.py", line 686, in main
    data = model.get_data()
  File "/home/runner/work/numpyro/numpyro/examples/hsgp.py", line 596, in get_data
    **self._get_floating_days(data),
  File "/home/runner/work/numpyro/numpyro/examples/hsgp.py", line 541, in _get_floating_days
    & data["date"].dt.day.ge(25),
AttributeError: 'DatetimeProperties' object has no attribute 'day_of_week'
test/test_examples.py::test_cpu[hsgp.py --model tywd --num-samples 10 --num-warmup 10 --num-chains 2] Running:
python examples/hsgp.py --model tywd --num-samples 10 --num-warmup 10 --num-chains 2
FAILED

Thanks @fehiepsi . My tests were all passing locally but, for some reason, they ran using a later version of pandas (1.3 instead of 1.1.5). I've changed the offending lines to use weekday instead of day_of_week and this fixed it for me. I've updated the PR.

[fehiepsi]

Thanks, @omarfsosa!