Vanilla implementation of Gaussian Process in PyMC3?

Questions
11
  • The PyMC GP API doesn’t take distance matrices directly. To reproduce the Rethinking example, I would recommend bypassing the GP API. If you check the source code for pm.gp.Latent, you’ll see it’s really not doing much. Feeding the distance matrix into pm.gp.Latent.prior as X is incorrect, its not supposed to be some distances of distances thing, it’s just over distance! The PyMC GP API is set up the way it is to make predicting easy.
  • The notation for a GP is usually f(x) \sim N(0, K(x, x')). So the GP can be evaluated over any x, but the covariance function K(x, x’) is a over all pair of x.
  • When you have a Gaussian likelihood and your model is:
\begin{align} f_i &\sim \text{MvN}(0, k(x_i, x')) \\ y_i &= N(f_i, \sigma^2) \end{align}

You can marginalize out f so you don’t have to sample it. \int p(y | f) p(f ) df. The result is

\mathbf{y} \sim \text{MvN}(0, \mathbf{K} + \sigma^2 \mathbf{I})

So these two models are equivalent (you should just say you have a GP with a normal likelihood either way), but the latter is a trick you can apply so you don’t need to sample the GP itself. If you do this, you can also use optimization to find good values for the covariance function hyperparameters. If you put priors on them, you’d call it MAP estimation. If you don’t, you’d call it type 2 maximum likelihood, or “maximizing the marginal likelihood”.

1 Like