pymc3 multiple gaussian process regression

Question

I am trying to run a gaussian process regression with two features by extending the first example in https://pymc-devs.github.io/pymc3/notebooks/GP-introduction.html

n = 20
X = np.array([list(a) for a in zip(np.sort(3*np.random.rand(n)), np.sort(3*np.random.rand(n)))])
y = np.random.normal(size=n)

with pm.Model() as model:
    # priors on the covariance function hyperparameters
    l = np.array([pm.Uniform('l1', 0, 10), pm.Uniform('l2', 0, 10)])

    # uninformative prior on the function variance
    log_s2_f = pm.Uniform('log_s2_f', lower=-10, upper=5)
    s2_f = pm.Deterministic('s2_f', tt.exp(log_s2_f))

    # uninformative prior on the noise variance
    log_s2_n = pm.Uniform('log_s2_n', lower=-10, upper=5)
    s2_n = pm.Deterministic('s2_n', tt.exp(log_s2_n))

    # covariance functions for the function f and the noise
    f_cov = s2_f * pm.gp.cov.ExpQuad(input_dim=2, lengthscales=l)

    y_obs = pm.gp.GP('y_obs', cov_func=f_cov, sigma=s2_n, observed={'X':X, 'Y':y})

Here the inputs of X and y are for testing the shape of the inputs. When I run the code I get a theano AsTensorError error which is traced to this in pymc3

/usr/local/lib/python2.7/site-packages/pymc3/gp/cov.pyc in square_dist(self, X, Z)
    124 
    125     def square_dist(self, X, Z):
--> 126         X = tt.mul(X, 1.0 / self.lengthscales)
    127         Xs = tt.sum(tt.square(X), 1)
    128         if Z is None:

Is it possible to run multiple gaussian regression in pymc3? If so I am sure I have messed up with the dimensions somewhere.

Answer 1

I found the solution to my question in the following blog which apparently is the point of reference for pymc3.

https://discourse.pymc.io/t/multidimensional-input-using-gaussian-process/128,

Instead of defining the covariance priors as an array of distributions, they are defined as multinomial distributions with the corresponding number of components. By changing the following line of the code above everything works as expected

with pm.Model() as model:
    # priors on the covariance function hyperparameters
    l = pm.Gamma('l', 1, 1, shape=2)