Bayesian model calibration with Gaussian Process

I have a different x_j 's that give y_j. Sorry for the mistake, I corrected it in the innitial question. I have still hard time to see how to exactly define it within the pm.Model() block. The gaussian processes also depend on some parameters \gamma which are in the case of f shared across z_j and y_j. I guess my question is, if I do something like this:

'gp1 = pm.gp.Marginal(mean_func1, cov_func1)
gp2 = pm.gp.Marginal(mean_func2, cov_func2)
gp3 = gp1 + gp2

f_1 = gp1.marginal_likelihood(“f_1”, X_1, y_1, noise)
f_2 = gp3.marginal_likelihood(“f_2”, X_2, y_2, noise)

will the inference be done properly?

In the case of \theta, it is a vector parameter that is not the parameter of the covariance functions but it plays a role of x in the covariance matrix. i.e. I have a input (\theta, x_i) for each z_i (same \theta for each)