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

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:
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:

- Have I correctly implemented the model in PyMC3? And is there any other way?
- I am not sure about the customized covariance function for delta(x), do you think it is correct?
- 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?