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)