Bayesian model calibration with Gaussian Process

Hmm, I am not sure how this solves the problem of need for something like:

with pm.Model() as model:
theta = pm.MvNormal("theta", mu = mean, cov = covariance)

cov_func = pm.gp.cov.ExpQuad(len(mu) + X.shape[1], lengthscale)

gp = pm.gp.Marginal(cov_func = cov_func)

#Now I need to perfrom something that stacks all (theta,x). Lets call this X1
f1 = gp.marginal_likelihood('y', X = X1, y = y, noise = sigma)

I do need to perform inference about unknown RV \theta and get its posterior distribution.