This is a use case where PyMC3 should really work well! The GP functions can definitely handle this. You can find a somewhat similar example here, where a flat prior is placed over the inducing point locations, which are input to a covariances function.
Is there a different theta for each x_i? or is it the same theta? (is it a scalar or a column vector?) If it’s a scalar, maybe try something like (you’ll have to double check these shapes)
theta = pm.Normal("theta", mu=0.0, sd=1.0)
theta_col = tt.alloc(theta, len(x1))
input_1 = tt.concatenate((x1[:,None], theta_col[:,None], axis=1)
Also, if by chance what you’re working on is open source / online id love to see a link!