# Using a Gaussian process at intermediate level in hierarchy

**URL:** <https://discourse.pymc.io/t/using-a-gaussian-process-at-intermediate-level-in-hierarchy/5102>\
**Category:** Questions\
**Created:** [May 16, 2020, 1:44am UTC](https://discourse.pymc.io/t/using-a-gaussian-process-at-intermediate-level-in-hierarchy/5102 "2020-05-16T01:44:04Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![lneufcourt](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/lneufcourt/32/2743_2.png) [@lneufcourt](https://discourse.pymc.io/u/lneufcourt)\
**Post date:** [May 16, 2020, 1:44am UTC](https://discourse.pymc.io/t/using-a-gaussian-process-at-intermediate-level-in-hierarchy/5102/1 "2020-05-16T01:44:04Z")

</div>

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 ?

---

<div class="post-metadata">

**Author:** ![tirthasheshpatel](https://avatars.discourse-cdn.com/v4/letter/t/ac91a4/32.png) [@tirthasheshpatel](https://discourse.pymc.io/u/tirthasheshpatel)\
**Post date:** [May 16, 2020, 4:56am UTC](https://discourse.pymc.io/t/using-a-gaussian-process-at-intermediate-level-in-hierarchy/5102/2 "2020-05-16T04:56:38Z")

</div>

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

```python
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
    # of ExpQuad Kernel
    ℓ = 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[0])

    # Sample!
    trace = pm.sample(1000, init="advi", chains=1)

```

You can see other kernel functions at [Mean and Covariance functions](https://docs.pymc.io/notebooks/GP-MeansAndCovs.html) docs. If you have a large dataset, you can use MarginalSparse Model. [Docs of MarginalSparse GP](https://docs.pymc.io/notebooks/GP-SparseApprox.html).

---

<div class="post-metadata">

**Author:** ![lneufcourt](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/lneufcourt/32/2743_2.png) [@lneufcourt](https://discourse.pymc.io/u/lneufcourt)\
**Post date:** [May 16, 2020, 8:39pm UTC](https://discourse.pymc.io/t/using-a-gaussian-process-at-intermediate-level-in-hierarchy/5102/3 "2020-05-16T20:39:29Z")

</div>

Thanks a lot @tirthasheshpatel. I think that this is just what I was looking for.

What would be the difference if I would use sample\_ppc to sample from the conditional, instead of sample ?

---

<div class="post-metadata">

**Author:** ![tirthasheshpatel](https://avatars.discourse-cdn.com/v4/letter/t/ac91a4/32.png) [@tirthasheshpatel](https://discourse.pymc.io/u/tirthasheshpatel)\
**Post date:** [May 17, 2020, 10:27am UTC](https://discourse.pymc.io/t/using-a-gaussian-process-at-intermediate-level-in-hierarchy/5102/4 "2020-05-17T10:27:00Z")

</div>

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.