Maybe I’m misunderstanding (I haven’t read the Kennedy paper) but your \theta values are your parameter values in your more complicated function correct? If you never vary x_i but treat it as known, why do you need it as an input to your gaussian process? Isn’t that information already included in your y vector? i.e.
y = np.array([some_function(x,thetas) for thetas in theta_vals]) #(A n x len(x) matrix)
So now your gp is really only a function of your \theta values which would then be your X?
Please disregard if I’m still not getting it.