Hi there,

I want to model an underlying function f:

f ~ MvN(mu, Sigma), where Sigma is fixed.

and I assume that my observations y are i.i.d with Gaussian Noise:

y_i ~ N(f_i, eps)

I tried:

```
with pm.Model() as model:
epsilon = pm.HalfNormal('epsilon', sd=1)
fs = pm.MvNormal('fs', mu=tt.zeros(len(ys)), cov=Sigma)
ys = pm.MvNormal('ys', mu=fs, cov=epsilon*tt.eye(len(fs)), observed=ys)
# or: ys = pm.Normal('ys', mu=fs, sd=epsilon, shape=len(_ys), observed=ys)
print(pm.find_MAP())
```

which yields the error â€śInvalid dimension for value: 0â€ť. ys has shape (n,)

I tried different things and cannot get it to work. Iâ€™d be very thankful for some help