Customized covariance for a subset of training data

Consider the following model:

y_i = f(x_i) + \delta(x_i) +\epsilon_i \quad \,, i=1,2,\cdots,n
y^*_i = f(x^*_i) +\epsilon^*_j \quad \,, j=n+1,n+2,\cdots,n+m

where f(x) and \delta(x) are two GPs.
We can define the combined data as Z=[y_1,...y_n,y^*_{n+1},\cdots,y^*_{n+m}] corresponding to X=[x_1,\cdots,x_{n+m}].
To define the model in PyMC3, I have done the following:

eta = pm.HalfCauchy("eta", beta=5)
lam = pm.Gamma("lam", alpha=1, beta=5)    
cov_ = eta ** 2.*, lam)
eta_del = pm.HalfCauchy("eta_del", beta=1)                              
lam_del = pm.Gamma("lam_del", alpha=1, beta=1, shape=1)     
cov_del_ = eta_del**2*ExpQuad_custom(1, lam_del, nTot=n+m, nCut=n)
noise_y = pm.HalfCauchy("sigma_y", beta=5)
noise_ys = pm.HalfCauchy("sigma_ys", beta=5)
y_=gp.marginal_likelihood("y_", X=X, y=Z, noise=noise_)

where X,Z have been defined using the hint in

The ‘CustomWhiteNoise’ is taken from

And I have tried to implement ExpQuad_custom as:

class ExpQuad_custom(Stationary_custom):
    def full(self, X, Xs=None):        
        X, Xs = self._slice(X, Xs)                
        cov_=tt.exp(-0.5 * self.square_dist(X, Xs))               
        if Xs is None: 
        return cov_

My Questions are:

  1. Have I correctly implemented the model in PyMC3? And is there any other way?
  2. I am not sure about the customized covariance function for delta(x), do you think it is correct?
  3. The implemented ExpQuad_custom works fine in the training but I have problem in prediction. Should the \delta(x) covariance be constructed between the first n training data and the test sample? If yes, is there any way to do it in PyMC3?