Jonathan Lalou's blog

Bill Engels and Chris Fonnesbeck, both brilliant software developers from PyMC Labs, delivered an insightful 90-minute tutorial at PyData Global 2024 titled “Making Gaussian Processes Useful.” Aimed at demystifying Gaussian processes (GPs) for practicing data scientists, their session bridged the gap between theoretical complexity and practical application. Using baseball analytics as a motivating example, Chris introduced Bayesian modeling and GPs, while Bill provided hands-on strategies for overcoming computational and identifiability challenges. This post explores their comprehensive approach, offering actionable insights for leveraging GPs in real-world scenarios.

Bayesian Inference and Probabilistic Programming

Chris kicked off the tutorial by grounding the audience in Bayesian inference, often implemented through probabilistic programming. He described it as writing software with partially random outputs, enabled by languages like PyMC that provide primitives for random variables. Unlike deterministic programming, probabilistic programming allows modeling distributions over variables, including functions via GPs. Chris explained that Bayesian inference involves specifying a joint probability model for data and parameters, using Bayes’ formula to derive the posterior distribution. This posterior reflects what we learn about unknown parameters after observing data, with the likelihood and priors as key components. The computational challenge lies in the normalizing constant, a multidimensional integral that probabilistic programming libraries handle numerically, freeing data scientists to focus on model specification.

Hierarchical Modeling with Baseball Data

To illustrate Bayesian modeling, Chris used the example of estimating home run probabilities for baseball players. He introduced a simple unpooled model where each player’s home run rate is modeled with a beta prior and a binomial likelihood, reflecting hits over plate appearances. Using PyMC, this model is straightforward to implement, with each line of code corresponding to a mathematical component. However, Chris highlighted its limitations: players with few at-bats yield highly uncertain estimates, leaning heavily on the flat prior. This led to the introduction of hierarchical modeling, or partial pooling, where individual home run rates are drawn from a population distribution with hyperparameters (mean and standard deviation). This approach shrinks extreme estimates, producing more realistic rates, as seen when comparing unpooled estimates (with outliers up to 80%) to pooled ones (clustered below 10%, aligning with real-world data like Barry Bonds’ 15% peak).

Gaussian Processes as a Hierarchical Extension

Chris transitioned to GPs, framing them as a generalization of hierarchical models for continuous predictors, such as player age affecting home run rates. Unlike categorical groups, GPs model relationships where similarity decreases with distance (e.g., younger players’ performance is more similar). A GP is a distribution over functions, parameterized by a mean function (often zero) and a covariance function, which defines how outputs covary based on input proximity. Chris emphasized two key properties of multivariate Gaussians—easy marginalization and conditioning—that make GPs computationally tractable despite their infinite dimensionality. By evaluating a covariance function at specific inputs, a GP yields a finite multivariate normal, enabling flexible, nonlinear modeling without explicitly parameterizing the function’s form.

Computational Challenges and the HSGP Approximation

One of the biggest hurdles with GPs is their computational cost, particularly for latent GPs used with non-Gaussian data like binomial home run counts. Chris explained that the posterior covariance function requires inverting a matrix, which scales cubically with the number of data points (e.g., thousands of players). This makes exact GPs infeasible for large datasets. To address this, he introduced the Hilbert Space Gaussian Process (HSGP) approximation, which reduces cubic compute time to linear by approximating the GP with a finite set of basis functions. These functions depend on the data, while coefficients rely on hyperparameters like length scale and amplitude. Chris demonstrated implementing an HSGP in PyMC to model age effects, specifying 100 basis functions and a boundary three times the data range, resulting in a model that ran in minutes rather than years.

Practical Debugging with GPs

Bill took over to provide practical tips for fitting GPs, emphasizing their sensitivity to priors and the need for debugging. He revisited the baseball example, modeling batting averages with a hierarchical model before introducing a GP to account for age effects. Bill showed that a standard hierarchical model treats players as exchangeable, pooling information equally across all players. A GP, however, allows local pooling, where players of similar ages inform each other more strongly. He introduced the exponentiated quadratic covariance function, which uses a length scale to define “closeness” in age and a scale parameter for effect size. Bill highlighted common pitfalls, such as small length scales reducing a GP to a standard hierarchical model or large length scales causing identifiability issues with intercepts, and provided solutions like informative priors (e.g., inverse gamma, log-normal) to constrain length scales to realistic ranges.

Advanced GP Modeling for Slugging Percentage

Bill concluded with a sophisticated model for slugging percentage, a metric reflecting hitting power, using 10 years of baseball data. The model included player, park, and season effects, with an HSGP to capture age effects. He initially used an exponentiated quadratic covariance function but encountered sampling issues (divergences), a common problem with GPs. Bill fixed this by switching to a Matern 5/2 covariance function, which assumes less smoothness and better suits real-world data, and adopting a centered parameterization for stronger age effects. These changes reduced divergences to near zero, producing a reliable model. The resulting age curve peaked at 26, aligning with baseball wisdom, and showed a decline for older players, demonstrating the GP’s ability to capture nonlinear trends.

Key Takeaways and Resources

Bill and Chris emphasized that GPs extend hierarchical models by enabling local pooling over continuous variables, but their computational and identifiability challenges require careful handling. Informative priors, appropriate covariance functions (e.g., Matern over exponential quadratic), and approximations like HSGP are critical for practical use. They encouraged using PyMC for its high-level interface and the Nutpie sampler for efficiency, while noting alternatives like GPFlow for specialized needs. Their GitHub repository, linked below, includes slides and notebooks for further exploration, making this tutorial a valuable resource for data scientists aiming to apply GPs effectively.

Links:

• Bill Engels’ LinkedIn
• Chris Fonnesbeck’s LinkedIn
• PyMC Labs
• PyMC
• GitHub Repository
• PyData Global 2024
• NumFOCUS