Customized covariance for a subset of training data

[saleh]

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.*pm.gp.cov.ExpQuad(1, 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)
noise_=CustomWhiteNoise(noise_y,noise_ys,n,m)        
gp_f=pm.gp.Marginal(cov_func=cov_)        
gp_del=pm.gp.Marginal(cov_func=cov_del_)        
gp=gp_f+gp_del
y_=gp.marginal_likelihood("y_", X=X, y=Z, noise=noise_)

where X,Z have been defined using the hint in

Bayesian model calibration with Gaussian Process

I am attempting to implement bayesian model calibration under the classical Kennedy-O’Hagan framework using PyMC3. Briefly, the setup is as follows: I have n observations z_i of the form z_i = f(x_i, \theta) + g(x_i) + e_i, where x_i are known imputs and \theta are unknown calibration parameters and e_i are iid error terms. We also have m model evaluations y_j of the form y_j = f(x^*_j, \theta_j), where both x^*_j (different than x_i above) and \theta_j are known. Therefore, the data consi…

The ‘CustomWhiteNoise’ is taken from
https://docs.pymc.io/notebooks/GP-MaunaLoa2.html

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: 
           N_=self.nTot-self.nCut
           a12=tt.tile(np.array([0]),(self.nCut,N_),ndim=2)
           a22=tt.tile(np.array([0]),(N_,self.nTot),ndim=2)        
           cov_=tt.concatenate([cov_,a12],axis=1)
           cov_=tt.concatenate([cov_,a22],axis=0)                
        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?