Hi all!

I’m trying to use a Gaussian Process as the input for the covariance matrix of an MvNormal likelihood:

```
with pm.Model() as m:
etasq = 1. + pm.HalfNormal('etasq', 0.25)
rhosq = 3. + pm.HalfNormal('rhosq', 0.25)
cov = etasq * pm.gp.cov.ExpQuad(input_dim=1, ls=rhosq)
gp = pm.gp.Latent(cov_func=cov)
SIGMA = gp.prior('SIGMA', X=Dist_mat ** 2)
a = pm.Normal("a", 0., 1.)
b = pm.Normal("b", 0., 0.5, shape=2)
mu = a + b[0] * M + b[1] * G
B = pm.MvNormal("B", mu, cov=SIGMA, observed=B_)
trace = pm.sample()
```

But running this prints a `ValueError: cov must be two dimensional`

, because, indeed, `SIGMA.tag.test_value.shape`

returns `(151,)`

.

I’m not very familiar yet with the GP module, so I’m not sure what to make of this Does anyone see a way forward?

Thansk a lot in advance

PS: I wanted to make the post as short as possible, but I can share the data if needed