 # 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)

~\Anaconda3\envs\pmf\lib\site-packages\pymc3\gp\gp.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

~\Anaconda3\envs\pmf\lib\site-packages\theano\tensor\basic.py in make_node(self, value,     *shape)

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

TypeError: ('The Alloc value to use has more dimensions than the specified dimensions', 2, 1)
``````

The only way I can avoid the error is by placing a diagonal matrix for Sigma, for example:

``````Sigma = [1,2,3,4]
``````

Any suggestions?

You’ll have to make a thin wrapper subclassing the base `Covariance` class here.

``````import theano.tensor as tt
import pymc3 as pm

class ConstantMatrix(pm.gp.cov.Covariance):
def __init__(self, cov):
self.cov = tt.as_tensor_variable(cov)

def full(self, X, Xs):
# covariance matrix known, not explicitly function of X
return self.cov

def diag(self, X):
return tt.diag(self.cov)
``````

Then wrap your `Sigma` with it:

``````Sigma = ConstantMatrix(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
``````

You should be good to go then for inference and prediction.

2 Likes