mre-additive-gp

import numpy as np
import pymc3 as pm

n = 4

X_2D_train = np.random.randn(n,2)
X_2D_star  = np.random.randn(n-1, 2)
y          = np.random.randn(n,1)

with pm.Model() as model:
    sigma1 = pm.HalfCauchy("η1", beta=5)
    sigma2 = pm.HalfCauchy("η2", beta=5)
    ls1 = pm.HalfCauchy('ls1', beta=10)
    ls2 = pm.HalfCauchy('ls2', beta=10)
    cov1 = sigma1 * pm.gp.cov.ExpQuad(2, ls=ls1, active_dims=[0])
    cov2 = sigma2 * pm.gp.cov.ExpQuad(2, ls=ls2, active_dims=[1]) 
    gp1 = pm.gp.Marginal(cov1)
    gp2 = pm.gp.Marginal(cov2)
    gp = gp1 + gp2
    sigma = pm.HalfCauchy("sigma", beta=3)
    f = gp.marginal_likelihood("f", X_2D_train, y, sigma)
    f_star = gp.conditional("f_star", X_2D_star)
    trace = pm.sample()