# Using a Gaussian process at intermediate level in hierarchy

Hello !

I am trying to use a Gaussian process at an intermediate level of a hierarchical model in PyMC3, e.g.

f(x)\sim\mathcal{GP}(\mu(x),k(x, x')),
y\sim \mathcal{N}(e^{f(x)}, \sigma),

where say \mu, \sigma and k are all known.

I am interested in computing the (conditional) posterior predictive p(y^*|y, x^*) at new locations x^*.
How can I do this nicely with PyMC3 ?

You can use Marginal Gaussian Process Model from pymc3.gp submodule.

import theano
import theano.tensor as tt
import pymc3 as pm

X = np.linspace(0, 1, num=5)[:, None] # Your dataset
y = np.exp(X.squeeze()) + np.random.randn(5) # Your observed data
Xnew = np.linspace(0, 1, num=10)[:, None] # Your new data points

with pm.Model() as model:
# Put priors on the lenght scale and amplitude
ℓ = pm.Gamma("ℓ", alpha=1, beta=2)
η = pm.HalfCauchy("η", beta=3)
# Create a kernel to use in GP Model
cov = η**2 * pm.gp.cov.Matern52(1, ℓ)
gp = pm.gp.Marginal(cov_func=cov)
# Get the marginal likelihood P(f | X, y)
f = gp.marginal_likelihood('f', X=X, y=y, noise=1e-2, is_observed=False)
# Put a prior on Noise
σ = pm.HalfCauchy("σ", beta=5)
# Putting a distributions on y.
ypred = pm.Normal('ypred', mu=tt.exp(f), sigma=σ, observed=y)
# fstar is the conditional P(f* | f, X, y)
# Then, we put a distribution on the conditional
fstar = gp.conditional('fstar', Xnew=Xnew, given={'X': X, 'y': ypred, 'f': f})
ystar = pm.Normal('ystar', mu=tt.exp(fstar+1e-4), sigma=σ, shape=Xnew.shape)

# Sample!

You can use any inference technique you like. Like pm.find_MAP, pm.sample, or variational inference. I don’t think using pm.sample_posterior_predictive (same as sample_ppc) without first calling a inference method makes sense. You can use pm.sample_posterior_predictive once you have a trace or MAP estimates.