Gaussian Process used to model varying intercepts in GLM

I realize you’d want to put the covariance function inside the model by using pm.gp.cov. But when you also include the marginal_likelihood method it wants you to give it both x and y values, but my y values haven’t been calculated yet. Here is an example of my failed reparameterizing:

with pm.Model() as m_13_7_reparam:

etasq = pm.HalfCauchy(‘etasq’, 1)
rhosq = pm.HalfCauchy(‘rhosq’, 1)
Kij = etasq * pm.gp.cov.ExpQuad(2,ls=rhosq) # two dimensions

g = pm.gp.MarginalSparse('g',cov_func=Kij)

# priors
a = pm.Normal('a', 0, 10)
bp = pm.Normal('bp', 0, 1)
# linear model
lam = pm.math.exp(a + g[dk.society.values] + bp * dk.logpop)
# likelihood
obs = pm.Poisson('total_tools', lam, observed=dk.total_tools)
# marginal likelihood
g = gp.marginal_likelihood('g', X=Dmat_, y=dk.total_tools ,noise=.5)
# inference
trace_13_7_reparam = pm.sample(1000, tune=1000)

To be clear I want to model some output function say f* not only in terms of the metric of similarity x that goes into the kernel function k(x,x) but also in terms of a linear submodel at each point x. So its like a hiearchical linear regression but with shrinkage not solely due to discrete membership to a larger group X_{[{1...i]}} but also in a continuous sense by the metric of similarity modeled by k(x,x).