I have the problem of how to combine a bayesian linear regression and gaussian process. As an example, I have two features, a 2D spatial GP trying to fit a model as:
Created the following code in PyMC but still not sure about the implementation:
X = data[['x1', 'x2']].T.values
X_gp = data[['gp0', 'gp1']].values
y = data.y.values
with pm3.Model() as model_mlr:
## Bayesian linear regression
Ξ±_blr = pm3.Normal('Ξ±_blr', mu=0, sd=10)
Ξ²_blr = pm3.Normal('Ξ²_blr', mu=0, sd=1, shape=num_features)
Ο = pm3.HalfCauchy('Ο', 5)
ΞΌ_blr = Ξ±_blr + pm3.math.dot(Ξ²_blr, X)
## The spatial GP
Ξ·_spatial_trend = pm3.Uniform('Ξ·_spatial_trend', lower=0, upper=100)
β_spatial_trend = pm3.Uniform('β_spatial_trend', lower=0, upper=100, shape=gp_input_dim)
cov_spatial_trend = (
Ξ·_spatial_trend**2 * pm3.gp.cov.ExpQuad(input_dim=gp_input_dim, ls=β_spatial_trend)
)
gp_spatial_trend = pm3.gp.Latent(cov_func=cov_spatial_trend)
gp_spatial_func = gp_spatial_trend.prior('gp_spatial_func', X=X_gp)
# noise model
cov_noise = pm3.gp.cov.WhiteNoise(Ο)
y_hat = pm3.Normal(
'y_hat',
mu=ΞΌ_blr + gp_spatial_func,
sd=Ο,
observed=y)
start = None
step = pm3.Metropolis()
trace_mlr = pm3.sample(30000, start=start, step=step, burn=10000)
chain = trace_mlr[1000:]

Should I use
gp.Latent
orgp.Marginal
? 
How prediction works in this case?
when I usegp.Latent
, I get complain about not converging.
Here is a sample data
sample_data.csv (37.8 KB)