# Gaussian Process -Statistical rethinking

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 * M + self.slopes * 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)

``````

What does the input data look like here? Is it standardized?

Also, we have a built-in Linear mean function already, so no need to code your own.

Just for fun, try using half-normals instead of taking the absolute values of a normal (then add your offsets). Probably not the issue, but can’t hurt.

Yes the distance matrix Dmat is standardized where maximum value is 1 for this particular example.

Great thanks; this works fine.

Does it mean like
`1 + HalfNormal("etasq",0.25)` instead of `pm.abs(pm.Normal("etasq,1,0.25))`?

Also went into another issue, while exploring (link below), As I believe that using PYMC’s expEquad kernel function (below equation) is providing different estimates than when calculated manually! resulting in some loss in precession but I might have mis-specify something here…