Gaussian Process with 2 dimension input

Questions
1

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

np.random.seed(123)

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ā€™)
plt.legend()
plt.show()

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.*pm.gp.cov.ExpQuad(nvar,ls=ls)
sigma = pm.HalfNormal(ā€˜sigmaā€™,1.0)

gp = pm.gp.Marginal(cov_func=cov)

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.

2

Please use this to format code https://help.github.com/articles/basic-writing-and-formatting-syntax/#quoting-code and https://help.github.com/articles/creating-and-highlighting-code-blocks/

It looks like this

import something_useful
import more

x = 4

with format as f:
3

There is an earlier discussion here: Cox Process in PyMC3

4

Did the error resolvedā€¦?..If yes can you share the codeā€¦Thanks in advance

5

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 = pm.gp.cov.ExpQuad(2, l) * sf2**2
    
    gp_spatial = pm.gp.MarginalSparse(cov_func=K, 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 pm.gp.cov.ExpQuad(2, l).

Regards

1 Like
6

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