Gaussian Process with 2 dimension input

I am new to pymc3 and gaussian process. I try to do a simple gaussian process with 2D inputs (i.e. X has 2 columns) and 1D output (i.e. y has 1 column). My code is like this:

import pymc3 as pm
import theano
import theano.tensor as tt
import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


nvar = 2
a0 = 2
a1 = 5

#define a function
def fxn(x0,x1,a0,a1):
return a0 * x0 + a1 * x1

#create the base function
m = 600
x0 = np.random.uniform(0,3,m)
x1 = np.random.uniform(0,5,m)
X = np.vstack([x0,x1])

y = (fxn(x0,x1,a0,a1).flatten()+np.random.normal(loc=0,scale=2,size=m))
y_true = fxn(x0,x1,a0,a1)

fig = plt.figure()
ax = fig.add_subplot(111, projection=‘3d’)
ax.scatter(x0, x1, y, marker=‘o’,color=‘k’,label=‘True Points’)
plt.plot(x0, x1, y_true,‘r.’,label=‘base function’)

y = y.reshape((m,1))

with pm.Model() as surrogate_model:
ls = pm.HalfNormal(‘l’,sd=2,shape=nvar)
eta = pm.HalfNormal(“eta”, sd=2, shape=1)
cov = eta**2.*,ls=ls)
sigma = pm.HalfNormal(‘sigma’,1.0)

gp =

y_ = gp.marginal_likelihood('y',X=X,y=y,noise=sigma)

trace = pm.sample(1000,tune=1000,njobs=2)

_ = pm.traceplot(trace,priors=[ls.distribution,eta.distribution,sigma.distribution])
summary = pm.summary(trace)

l_means = np.array([summary[‘mean’][key] for key in summary[‘mean’].keys() if ‘l_’ in key]).reshape(1,nvar)
eta_mean = np.array([summary[‘mean’][key] for key in summary[‘mean’].keys() if ‘eta’ in key]).reshape(1,1)
sigma_mean = [summary[‘mean’][key] for key in summary[‘mean’].keys() if ‘sigma’ in key]

mm = 600
xx0 = np.linspace(0,6,mm)
xx1 = np.linspace(0,10,mm)
XX = np.vstack([xx0,xx1])

with surrogate_model:
f_pred = gp.conditional(‘f_pred’,XX)

with surrogate_model:
pred_samples = pm.sample_ppc(trace,vars=[f_pred],samples=1000)

It gives me “incompatible dimensions” error when I tried to do “gp.conditional”, can someone help me with this?

Thanks in advance.

Please use this to format code and

It looks like this

import something_useful
import more

x = 4

with format as f:

There is an earlier discussion here: Cox Process in PyMC3

Did the error resolved…?..If yes can you share the code…Thanks in advance

Hi, you can check this notebook of Multidimensional GP with 2 predictors from Chris Fonnesbeck’s repos on this link

I copy the model here:

with pm.Model() as spatial_model:
    l = pm.HalfCauchy("l", beta=3, shape=(2,))
    sf2 = pm.HalfCauchy("sf2", beta=3)
    sn2 = pm.HalfCauchy("sn2", beta=3)

    K =, l) * sf2**2
    gp_spatial =, approx="FITC")
    obs = gp_spatial.marginal_likelihood("obs", X=X_obs, Xu=Xu, y=y_obs, noise=sn2)

    mp = pm.find_MAP()

As X_obs has 2 columns, we need to set 2 in, l).


1 Like

Hi, thanks for helping me out…Also, is there any reference to study the heatmap of output (mean) of gaussian process with 2d input?