  Gaussian processes (GPs) are sophisticated distributions to model functional
data. Whilst theoretically appealing, they are computationally cumbersome
except for small datasets. We implement two methods for scaling GP inference in
Stan: First, a general sparse approximation using a directed acyclic dependency
graph; second, a fast, exact method for regularly spaced data modeled by GPs
with stationary kernels using the fast Fourier transform. Based on benchmark
experiments, we offer guidance for practitioners to decide between different
methods and parameterizations. We consider two real-world examples to
illustrate the package. The implementation follows Stan's design and exposes
performant inference through a familiar interface. Full posterior inference for
ten thousand data points is feasible on a laptop in less than 20 seconds.
Details on how to get started using the popular interfaces cmdstanpy for Python
and cmdstanr for R are provided.
