Adding covariance matrix to Gaussian Processes in Pymc3

I am trying to add a covariance matrix inside a Gaussian Processes using pymc3, this is my code:

import matplotlib.pyplot as plt
import numpy as np
import pymc3 as pm
%matplotlib inline

#Data
X = np.array([1,2,3,4])
Y= np.array([1,4,9,16])
Sigma = np.matrix([[1,0.1, 0,0], [0.1, 2,0.2,0.4], [0, 0.2,3,0.4],[0, 0.4,0.4,4]]) #Covariance  Matrix
le= len(X)

#Reorder the data stored in X
Xg = X.reshape((le,1))

with pm.Model() as model:
ℓ = pm.Uniform("ℓ", lower=0,upper=2.5)
η =  pm.Uniform("η", lower=0,upper=0.5)
cov = η ** 2 * pm.gp.cov.Matern52(1, ℓ)
gp = pm.gp.Marginal(cov_func=cov)
y_ = gp.marginal_likelihood("y", X=Xg, y=Y, noise=Sigma)
mp = pm.find_MAP()

However I get the following error:

TypeError                                 Traceback (most recent call last)
<ipython-input-69-4911417cefce> in <module>
---> 13     y_ = gp.marginal_likelihood("y", X=Xg, y=Y, noise=Sigma)

~Anaconda3envspmflibsite-packagespymc3gpgp.py in marginal_likelihood(self, name, X, y,     noise, is_observed, **kwargs)
441         if not isinstance(noise, Covariance):
442             noise = pm.gp.cov.WhiteNoise(noise)
--> 443         mu, cov = self._build_marginal_likelihood(X, noise)
444         self.X = X
445         self.y = y

~Anaconda3envspmflibsite-packagestheanotensorbasic.py in make_node(self, value,     *shape)

2991                 " than the specified dimensions",
2992                 v.ndim,
-> 2993                 len(sh),
2994             )