Using PyMC5, I know:

- How to perform Gaussian Process regression.
- How to construct priors and likelihood, then make sampling to get posterior.

But I’m not sure, in PyMC5, how can I use the GP regression model to construct the likelihood. To illustrate, I shall use the example in Introductory Overview of PyMC — PyMC dev documentation.

```
with pm.Model() as inference_model:
alpha = pm.Normal("alpha", mu=0, sigma=10)
beta = pm.Normal("beta", mu=0, sigma=10, shape=2)
sigma = pm.HalfNormal("sigma", sigma=1)
mu = alpha + beta[0] * X1 + beta[1] * X2 # REPLACE WITH GP MODEL???
Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=Y)
idata = pm.sample()
```

In this example, `mu`

is calculated explicitly from an analytic expression. But what if we must use the Gaussian Process model to calculate `mu`

? \mu = \mathcal{GP}(X)

The GP may be modeled as following:

```
with pm.Model() as model_gp:
ls = pm.Gamma('ls', 2, 1, shape=3)
cov = pm.gp.cov.ExpQuad(3, ls=ls)
gp = pm.gp.Marginal(cov_func=cov)
f_pred = gp.marginal_likelihood('f_pred', X=X_train, y=y_train, sigma=sigma)
gp_data = pm.sample()
```

I had posted another relevant question at `predictt` in PyMC5?, explaining that I was able to do it in PyMC3, but not in PyMC5. But there was no reply. So I try to ask it again in a more general way here. I apologize if this is undesirable.