Adding covariance matrix to Gaussian Processes in Pymc3

Questions
1

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             )
   2995         otype = TensorType(dtype=v.dtype, broadcastable=bcast)

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?

2

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