Hi all!
I am trying to recreate the Mcelreath’s lecture problem (snapshot below) 
on Gaussian Process in PYMC.
I need some assistance in understanding the proper method of implementing the above model using GP classes in PYMC? I have referenced the PYMC lecture translation (link) as a guide, however it utilizes different priors. Despite this, I was able to run the model to some extent, but it is producing different parameter estimates and is running slower than the hard-coded model(i.e without using GP clases). I am wondering if you could help me identify where I may have gone wrong in my implementation.
Below is my implementation…
class LinearMean(pm.gp.mean.Mean):
def __init__(self, intercept, slopes):
self.intercept = intercept
self.slopes = slopes
def __call__(self, X):
"""In the trace summary, b0 will be bM and b1 will be bG"""
return self.intercept + self.slopes[0] * M + self.slopes[1] * G
with pm.Model() as mB_OU3:
a = pm.Normal("a", 0.0, 1.0)
b = pm.Normal("b", 0.0, 0.5, shape=2)
# specify the mean function:
mu = LinearMean(intercept=a, slopes=b)
eta_raw = pm.Normal('eta_raw', 1., 0.25)
rho_raw = pm.Normal('rho_raw', 3., 0.25)
etasq = pm.Deterministic("etasq",at.abs(eta_raw))
rhosq = pm.Deterministic("rhosq",at.abs(rho_raw))
SIGMA = etasq * pm.gp.cov.Exponential(input_dim=1, ls_inv=2*rhosq)
# specify the GP:
gp = pm.gp.Marginal(mean_func=mu, cov_func=SIGMA)
sigma = pm.Exponential("sigma", 1.0)
B = gp.marginal_likelihood("B", X=Dmat, y=B_, sigma=sigma)
idata_mB_OU3=pm.sample(cores=2,target_accept=0.98)