I am sorry, I reallized that I might not have been the most clear. I will try to elaborate. I am trying to model something like this
z_i = f(x_i, \theta) + e_i
y_j = f(x^*_j, \theta^*_j) +e_j
where
f \sim GP\{m(.,.), \Sigma[(.,.)(.,.)]\}
with
m(x,\theta) = x^T \theta
and
\Sigma[(x,\theta)(x',\theta')] = e^{-\sum (\frac{\theta_i - \theta'_i}{\gamma_i})^2 -\sum (\frac{x_j - x'_j}{\eta_i})^2 }.
And my observed responses are D = (z_1, \dots, z_n, y_1, \dots, y_m) with covariates \{(x_1, \theta), \dots, (x_n,\theta), (x^*_1, \theta^*_1), \dots, (x^*_m, \theta^*_m)\}. Where I have x_i, x^*_j, \theta^*_j as known fixed values (it is part of my data) and unknow calibration parameters \theta.
My goal is to do inference abut all \gamma_i, \eta_j, \theta and get their posterior distributions. So that is why I need to use unknown \theta as an input into a covariance function, and that is what I am struggling with.