My data:

I have a matrix of `N`

predictors, with `t`

observations for the model training:

`X (t x N)`

I have a target vector (the thing I would like to predict):

`y (t,)`

I have a vector of test cases for each predictor, in order to perform the forecast:

`X_new (1 x N)`

I would like to perform a multivariate Gaussian Process regression with hyperpriors `p(θ)`

, whereby optimum hyperparameters `(θ = ℓ, a, σ)`

(see code below) are first determined by Metropolis MCMC and **then** the optimum regression weight parameters `β`

are determined (`y = βX + ε`

) in order to perform a forecast `y_new = βX_new`

. However I’m not quite sure if I am implementing this correctly in pyMC3…

So my understanding is that parameters `μ`

and `Σ`

(see code below) correspond to the mean and covariance of the model prior `p(β)`

. I’m wondering, does my code below just take **one** random sample from the defined distributions of each of the hyperparameters and then perform MCMC sampling on the model posterior `p(β | y, X, θ)`

to get the regression weights? Or does it perform MCMC sampling on the posterior of the hyperparameters `p(θ | y, X)`

**and** the model posterior `p(β | y, X, θ)`

? If it’s the former, how can I change my code to implement hyperpriors and optimise them with Metropolis MCMC, and finally perform a prediction of `y`

?

Side note - the prior on my hyperparameters would be weakly informative. All I know is that they should all be > 0.

Thanks.

```
with pm.Model() as model:
ℓ = pm.Gamma('ℓ', alpha=1, beta=10, shape=X.shape[1])
a = pm.HalfCauchy('a', beta=10)
σ = pm.HalfCauchy("σ", beta=10)
μ = pm.gp.mean.Zero()
Σ = a * pm.gp.cov.ExpQuad(input_dim=2,active_dims=np.arange(X.shape[1]),ls=ℓ)
gp = pm.gp.Marginal(mean_func=μ, cov_func=Σ)
y_ = gp.marginal_likelihood('y_',X=X, y=y, noise=σ)
MAP = pm.find_MAP()
method = pm.Metropolis()
trace = pm.sample(1000, method, MAP, random_seed=1, progressbar=True)
func = gp.conditional('func', Xnew=X_new)
fmean, fvar = gp.predict(Xnew=X_new, point=MAP)
```