Merge branch 'main' into statistical_rethinking

Author: fonnesbeck

examples/bart/bart_categorical_hawks.ipynb

@@ -14,7 +14,7 @@ myst:
     pip_dependencies: pymc-bart
 ---
 
-+++ {"editable": true, "slideshow": {"slide_type": ""}}
++++ {"slideshow": {"slide_type": ""}}
 
 (bart_categorical)=
 # Categorical regression
@@ -136,11 +136,11 @@ It may be that some of the input variables are not informative for classifying b
 
 ```{code-cell} ipython3
 ---
-editable: true
 slideshow:
   slide_type: ''
 ---
-pmb.plot_variable_importance(idata, μ, x_0, method="VI", random_seed=RANDOM_SEED);
+vi_results = pmb.compute_variable_importance(idata, μ, x_0, method="VI", random_seed=RANDOM_SEED)
+pmb.plot_variable_importance(vi_results);
 ```
 
 It can be observed that with the covariables `Hallux`, `Culmen`, and `Wing` we achieve the same R$^2$ value that we obtained with all the covariables, this is that the last two covariables contribute less than the other three to the classification. One thing we have to take into account in this is that the HDI is quite wide, which gives us less precision on the results, later we are going to see a way to reduce this.    
@@ -152,7 +152,7 @@ It can be observed that with the covariables `Hallux`, `Culmen`, and `Wing` we a
 Let's check the behavior of each covariable for each species with `pmb.plot_pdp()`, which shows the marginal effect a covariate has on the predicted variable, while we average over all the other covariates.  
 
 ```{code-cell} ipython3
-pmb.plot_pdp(μ, X=x_0, Y=y_0, grid=(5, 3), figsize=(6, 9));
+pmb.plot_pdp(μ, X=x_0, Y=y_0, grid=(5, 3), figsize=(12, 7));
 ```
 
 The pdp plot, together with the Variable Importance plot, confirms that `Tail` is the covariable with the smaller effect over the predicted variable. In the Variable Importance plot `Tail` is the last covariable to be added and does not improve the result, in the pdp plot `Tail` has the flattest response. For the rest of the covariables in this plot, it's hard to see which of them have more effect over the predicted variable, because they have great variability, showed in the HDI wide, same as before later we are going to see a way to reduce this variability. Finally, some variability depends on the amount of data for each species, which we can see  in the `counts` from one of the covariables using Pandas `.describe()` and grouping the data from "Species" with `.groupby("Species")`.  
@@ -215,11 +215,14 @@ with pm.Model(coords=coords) as model_t:
 Now we are going to reproduce the same analyses as before.  
 
 ```{code-cell} ipython3
-pmb.plot_variable_importance(idata_t, μ_t, x_0, method="VI", random_seed=RANDOM_SEED);
+vi_results = pmb.compute_variable_importance(
+    idata_t, μ_t, x_0, method="VI", random_seed=RANDOM_SEED
+)
+pmb.plot_variable_importance(vi_results);
 ```
 
 ```{code-cell} ipython3
-pmb.plot_pdp(μ_t, X=x_0, Y=y_0, grid=(5, 3), figsize=(6, 9));
+pmb.plot_pdp(μ_t, X=x_0, Y=y_0, grid=(5, 3), figsize=(12, 7));
 ```
 
 Comparing these two plots with the previous ones shows a marked reduction in the variance for each one. In the case of `pmb.plot_variable_importance()` there are smallers error bands with an R$^{2}$ value more close to 1. And for `pm.plot_pdp()` we can see thinner bands and a reduction in the limits on the y-axis, this is a representation of the reduction of the uncertainty due to adjusting the trees separately. Another benefit of this is that is more visible the behavior of each covariable for each one of the species.   
@@ -254,7 +257,8 @@ all
 ```
 
 ## Authors
-- Authored by [Pablo Garay](https://github.com/PabloGGaray) and [Osvaldo Martin](https://aloctavodia.github.io/) in May, 2024  
+- Authored by [Pablo Garay](https://github.com/PabloGGaray) and [Osvaldo Martin](https://aloctavodia.github.io/) in May, 2024
+- Updated by Osvaldo Martin in Dec, 2024
 
 +++
 

examples/bart/bart_categorical_hawks.myst.md

@@ -147,6 +147,7 @@ The fit looks good! In fact, we see that the mean and variance increase as a fun
 - Authored by [Juan Orduz](https://juanitorduz.github.io/) in Feb, 2023 
 - Rerun by Osvaldo Martin in Mar, 2023
 - Rerun by Osvaldo Martin in Nov, 2023
+- Rerun by Osvaldo Martin in Dec, 2024
 
 +++
 

examples/bart/bart_heteroscedasticity.ipynb

@@ -199,14 +199,15 @@ Finally, like with other regression methods, we should be careful that the effec
 
 ### Variable importance
 
-As we saw in the previous section a partial dependence plot can visualize and give us an idea of how much each covariable contributes to the predicted outcome. Moreover, PyMC-BART provides a novel method to assess the importance of each variable in the model. You can see an example in the following figure. 
+As we saw in the previous section a partial dependence plot can visualize give us an idea of how much each covariable contributes to the predicted outcome. Moreover, PyMC-BART provides a novel method to assess the importance of each variable in the model. You can see an example in the following figure. 
 
 On the x-axis we have the number of covariables and on the y-axis R² (the the square of the Pearson correlation coefficient) between the predictions made for the full model (all variables included) and the restricted models, those with only a subset of the variables. 
 
 In this example, the most important variable is `hour`, then `temperature`, `humidity`, and finally `workingday`.  Notice that the first value of R², is the value of a model that only includes the variable `hour`, the second R² is for a model with two variables, `hour` and  `temperature`, and so on. Besides this ranking, we can see that even a model with a single component, `hour`, is very close to the full model. Even more, the model with two components `hour`, and `temperature` is on average indistinguishable from the full model. The error bars represent the 94 \% HDI from the posterior predictive distribution. This means that we should expect a model with only `hour` and `temperature` to have a similar predictice performance than a model with the four variables, `hour`,  `temperature`, `humidity`, and `workingday`.
 
 ```{code-cell} ipython3
-pmb.plot_variable_importance(idata_bikes, μ, X);
+vi_results = pmb.compute_variable_importance(idata_bikes, μ, X)
+pmb.plot_variable_importance(vi_results);
 ```
 
 `plot_variable_importance` is fast because it makes two assumptions:
@@ -405,6 +406,7 @@ This plot helps us understand the reason behind the bad performance on the test
 * Juan Orduz added out-of-sample section in Jan, 2023
 * Updated by Osvaldo Martin in Mar, 2023
 * Updated by Osvaldo Martin in Nov, 2023
+* Updated by Osvaldo Martin in Dec, 2024
 
 +++
 

examples/bart/bart_heteroscedasticity.myst.md

@@ -149,6 +149,7 @@ We can see that when we use a Normal likelihood, and from that fit we compute th
 * Authored by Osvaldo Martin in Jan, 2023
 * Rerun by Osvaldo Martin in Mar, 2023
 * Rerun by Osvaldo Martin in Nov, 2023
+* Rerun by Osvaldo Martin in Dec, 2024
 
 +++
 

examples/bart/bart_introduction.ipynb

@@ -10,9 +10,16 @@ kernelspec:
   name: python3
 ---
 
+(GP-TProcess)=
 # Student-t Process
 
-PyMC3 also includes T-process priors.  They are a generalization of a Gaussian process prior to the multivariate Student's T distribution.  The usage is identical to that of `gp.Latent`, except they require a degrees of freedom parameter when they are specified in the model.  For more information, see chapter 9 of [Rasmussen+Williams](http://www.gaussianprocess.org/gpml/), and [Shah et al.](https://arxiv.org/abs/1402.4306).
+:::{post} August 2017
+:tags: t-process, gaussian process, bayesian non-parametrics
+:category: intermediate
+:author: Bill Engels
+:::
+
+PyMC also includes T-process priors.  They are a generalization of a Gaussian process prior to the multivariate Student's T distribution.  The usage is identical to that of `gp.Latent`, except they require a degrees of freedom parameter when they are specified in the model.  For more information, see chapter 9 of [Rasmussen+Williams](http://www.gaussianprocess.org/gpml/), and [Shah et al.](https://arxiv.org/abs/1402.4306).
 
 Note that T processes aren't additive in the same way as GPs, so addition of `TP` objects are not supported.
 
@@ -23,10 +30,13 @@ Note that T processes aren't additive in the same way as GPs, so addition of `TP
 The following code draws samples from a T process prior with 3 degrees of freedom and a Gaussian process, both with the same covariance matrix.
 
 ```{code-cell} ipython3
+import arviz as az
 import matplotlib.pyplot as plt
 import numpy as np
-import pymc3 as pm
-import theano.tensor as tt
+import pymc as pm
+import pytensor.tensor as pt
+
+from pymc.gp.util import plot_gp_dist
 
 %matplotlib inline
 ```
@@ -39,33 +49,34 @@ n = 100  # The number of data points
 X = np.linspace(0, 10, n)[:, None]  # The inputs to the GP, they must be arranged as a column vector
 
 # Define the true covariance function and its parameters
-ℓ_true = 1.0
-η_true = 3.0
-cov_func = η_true**2 * pm.gp.cov.Matern52(1, ℓ_true)
+ell_true = 1.0
+eta_true = 3.0
+cov_func = eta_true**2 * pm.gp.cov.Matern52(1, ell_true)
 
 # A mean function that is zero everywhere
 mean_func = pm.gp.mean.Zero()
 
 # The latent function values are one sample from a multivariate normal
 # Note that we have to call `eval()` because PyMC3 built on top of Theano
-tp_samples = pm.MvStudentT.dist(mu=mean_func(X).eval(), cov=cov_func(X).eval(), nu=3).random(size=8)
+tp_samples = pm.draw(pm.MvStudentT.dist(mu=mean_func(X).eval(), scale=cov_func(X).eval(), nu=3), 8)
 
 ## Plot samples from TP prior
 fig = plt.figure(figsize=(12, 5))
-ax = fig.gca()
-ax.plot(X.flatten(), tp_samples.T, lw=3, alpha=0.6)
-ax.set_xlabel("X")
-ax.set_ylabel("y")
-ax.set_title("Samples from TP with DoF=3")
+ax0 = fig.gca()
+ax0.plot(X.flatten(), tp_samples.T, lw=3, alpha=0.6)
+ax0.set_xlabel("X")
+ax0.set_ylabel("y")
+ax0.set_title("Samples from TP with DoF=3")
 
 
-gp_samples = pm.MvNormal.dist(mu=mean_func(X).eval(), cov=cov_func(X).eval()).random(size=8)
+gp_samples = pm.draw(pm.MvNormal.dist(mu=mean_func(X).eval(), cov=cov_func(X).eval()), 8)
 fig = plt.figure(figsize=(12, 5))
-ax = fig.gca()
-ax.plot(X.flatten(), gp_samples.T, lw=3, alpha=0.6)
-ax.set_xlabel("X")
-ax.set_ylabel("y")
-ax.set_title("Samples from GP");
+ax1 = fig.gca()
+ax1.plot(X.flatten(), gp_samples.T, lw=3, alpha=0.6)
+ax1.set_xlabel("X")
+ax1.set_ylabel("y")
+ax1.set_ylim(ax0.get_ylim())
+ax1.set_title("Samples from GP");
 ```
 
 ## Poisson data generated by a T process
@@ -79,16 +90,16 @@ n = 150  # The number of data points
 X = np.linspace(0, 10, n)[:, None]  # The inputs to the GP, they must be arranged as a column vector
 
 # Define the true covariance function and its parameters
-ℓ_true = 1.0
-η_true = 3.0
-cov_func = η_true**2 * pm.gp.cov.ExpQuad(1, ℓ_true)
+ell_true = 1.0
+eta_true = 3.0
+cov_func = eta_true**2 * pm.gp.cov.ExpQuad(1, ell_true)
 
 # A mean function that is zero everywhere
 mean_func = pm.gp.mean.Zero()
 
 # The latent function values are one sample from a multivariate normal
 # Note that we have to call `eval()` because PyMC3 built on top of Theano
-f_true = pm.MvStudentT.dist(mu=mean_func(X).eval(), cov=cov_func(X).eval(), nu=3).random(size=1)
+f_true = pm.draw(pm.MvStudentT.dist(mu=mean_func(X).eval(), scale=cov_func(X).eval(), nu=3), 1)
 y = np.random.poisson(f_true**2)
 
 fig = plt.figure(figsize=(12, 5))
@@ -102,23 +113,22 @@ plt.legend();
 
 ```{code-cell} ipython3
 with pm.Model() as model:
-    ℓ = pm.Gamma("ℓ", alpha=2, beta=2)
-    η = pm.HalfCauchy("η", beta=3)
-    cov = η**2 * pm.gp.cov.ExpQuad(1, ℓ)
+    ell = pm.Gamma("ell", alpha=2, beta=2)
+    eta = pm.HalfCauchy("eta", beta=3)
+    cov = eta**2 * pm.gp.cov.ExpQuad(1, ell)
 
     # informative prior on degrees of freedom < 5
-    ν = pm.Gamma("ν", alpha=2, beta=1)
-    tp = pm.gp.TP(cov_func=cov, nu=ν)
+    nu = pm.Gamma("nu", alpha=2, beta=1)
+    tp = pm.gp.TP(scale_func=cov, nu=nu)
     f = tp.prior("f", X=X)
 
-    # adding a small constant seems to help with numerical stability here
-    y_ = pm.Poisson("y", mu=tt.square(f) + 1e-6, observed=y)
+    pm.Poisson("y", mu=pt.square(f), observed=y)
 
-    tr = pm.sample(1000)
+    tr = pm.sample(target_accept=0.9, nuts_sampler="nutpie", chains=2)
 ```
 
 ```{code-cell} ipython3
-pm.traceplot(tr, var_names=["ℓ", "ν", "η"]);
+az.plot_trace(tr, var_names=["ell", "nu", "eta"]);
 ```
 
 ```{code-cell} ipython3
@@ -131,24 +141,37 @@ with model:
 
 # Sample from the GP conditional distribution
 with model:
-    pred_samples = pm.sample_posterior_predictive(tr, vars=[f_pred], samples=1000)
+    pm.sample_posterior_predictive(tr, var_names=["f_pred"], extend_inferencedata=True)
 ```
 
 ```{code-cell} ipython3
 fig = plt.figure(figsize=(12, 5))
 ax = fig.gca()
-from pymc3.gp.util import plot_gp_dist
 
-plot_gp_dist(ax, np.square(pred_samples["f_pred"]), X_new)
+f_pred_samples = np.square(
+    az.extract(tr.posterior_predictive).astype(np.float32)["f_pred"].values
+).T
+plot_gp_dist(ax, f_pred_samples, X_new)
 plt.plot(X, np.square(f_true), "dodgerblue", lw=3, label="True f")
 plt.plot(X, y, "ok", ms=3, alpha=0.5, label="Observed data")
 plt.xlabel("X")
 plt.ylabel("True f(x)")
-plt.ylim([-2, 20])
 plt.title("Conditional distribution of f_*, given f")
 plt.legend();
 ```
 
+## Authors
+
+* Authored by Bill Engels
+* Updated by Chris Fonnesbeck to use PyMC v5
+
++++
+
+## References
+:::{bibliography}
+:filter: docname in docnames
+:::
+
 ```{code-cell} ipython3
 %load_ext watermark
 %watermark -n -u -v -iv -w