Hi everyone!

I’m a bit new to PyMC (and to Gaussian processes in general), so sorry in advance for some imprecision.

I have a function that depends on three parameters f(x, \alpha, \beta), that will be part of my model. However, I only have access to its values on a discrete set of (x, \alpha, \beta) points. My idea is to perform an interpolation via Gaussian processes. In this way, whenever I need the value of f I would make a sampling of the Gaussian process.

My question is then: after having done the “interpolation”, how can I sample it within a different Bayesian model? Something like:

```
with pm.Model() as model:
x = pm.Data('x', X, mutable=True)
# Prior
alpha = pm.Uniform('alpha', 0, 20)
beta = pm.Uniform('beta', 0, 20)
sigma = pm.HalfNormal('sigma', 10)
# Part from GP (what I don't know how to do...)
prms_new = at.stack((x, alpha, beta)).T
with gp_model:
f_pred = gp.conditional('f_pred', prms_new)
y = pm.Deterministic('y', f_pred)
# Likelihood
pm.Normal('obs', mu=y, sigma=sigma, observed=Y_obs)
trace = pm.sample()
```

I’ve come across similar questions in other threads, but if this question was answered before, I’m afraid I didn’t understand it…

Thanks in advance!