Gaussian Process -Statistical rethinking

Hi all!
I am trying to recreate the Mcelreath’s lecture problem (snapshot below) :point_down: 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(
    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 *, ls_inv=2*rhosq)
    # specify the GP:
    gp =, cov_func=SIGMA)
    sigma = pm.Exponential("sigma", 1.0)
    B = gp.marginal_likelihood("B", X=Dmat, y=B_, sigma=sigma)

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.

Hi @fonnesbeck

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…