# Gaussian Process regression of a 3D scalar function

Hi all,

I am new to Pymc3 and Gaussian Processes so I apologize if this question seems trivial or I am missing something fundamental. I am trying to model a scalar function with multiple inputs. I have successfully done this using Pymc3 for a 2d toy function ( f(x,y) where x and y are coordinates ). However, I am unable to extend my model for a function with three inputs. I am trying to model now a toy function of the form f(x,y,z)=cos(x)+cos(y)+cos(z) + 2*(x2 + y2 + z**2). The code I am using is:

import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt

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

#Define Test Function

def true_funct(x,y,z):
val = np.cos(x) + np.cos(y) + np.cos(z) + 2*(x2 + y2 + z**2)

``````return val
``````

# Create Observations

n_points = 40
noise_true = 0.001
data_x = np.random.uniform(-1,1,n_points)
data_y = np.random.uniform(-1,1,n_points)
data_z = np.random.uniform(-1,1,n_points)

X_train = np.concatenate((data_x[:,None],data_y[:,None],data_z[:,None]),axis=1)
Y_true = true_funct(data_x,data_y,data_z)[:,None]
Y_train = Y_true + np.random.normal(0,noise_true, Y_true.shape)[:,None]

# sigma_noise = noise_true

with pm.Model() as model:

``````#List of Prior parameters
#-------------------------------

l_scale = pm.HalfCauchy('l_scale', 1, shape=3) #length scale parameter
sigma_f = pm.HalfCauchy('sigma_f', 1) #amplitude parameter
sigma_noise = pm.HalfCauchy('sigma_noise', .1)

#build the mean and cov matrix
cov_funct=cov_funct(X_train)
K = tt.fill_diagonal(cov_funct, sigma_f**2. + sigma_noise**2.)
mean_func = pm.gp.mean.Zero()

obs = pm.MvNormal('obs', mean_func(X_train), K, observed=Y_train)
trace = pm.sample(1000,chains=1,cores=1, step=pm.NUTS())
``````

pm.traceplot(trace, var_names=[‘sigma_f’,‘l_scale’]);

The issue I am having is that the length scales I get are extremely large and I pretty much get a constant response when sampling the GP at new locations. I have used successfully the above code for a function with: f(x,y) = sin(norm(x,y)) but I am having difficulty now with the three input case and I do not see why.

Thank you very much in advance for your time and apologies if this has been asked before.

I wanted to do something similar to you, so I used your example (reduced to 2D), but I can’t get either your example or my reduced version to work. I’m using theano 1.0.5 & pymc 3.9.3 on Linux (conda-forge install).

``````import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt

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

#Define Test Function

def true_funct(x,y):
val = np.cos(x) + np.cos(y) + 2*(x*x + y*y)
return val

# Create Observations

n_points = 40
noise_true = 0.001
data_x = np.random.uniform(-1,1,n_points)
data_y = np.random.uniform(-1,1,n_points)

X_train = np.concatenate((data_x[:,None],data_y[:,None]),axis=1)
Y_true = true_funct(data_x,data_y)[:,None]
Y_train = Y_true + np.random.normal(0,noise_true, Y_true.shape)[:,None]
sigma_noise = noise_true

with pm.Model() as model:
#List of Prior parameters
#-------------------------------

l_scale = pm.HalfCauchy('l_scale', 1, shape=2) #length scale parameter
sigma_f = pm.HalfCauchy('sigma_f', 1) #amplitude parameter
sigma_noise = pm.HalfCauchy('sigma_noise', .1)

# build the mean and cov matrix
cov_funct=cov_funct(X_train)
K = tt.fill_diagonal(cov_funct, sigma_f**2. + sigma_noise**2.)
mean_func = pm.gp.mean.Zero()

obs = pm.MvNormal('obs', mean_func(X_train), K, observed=Y_train)
trace = pm.sample(1000,chains=2,cores=2, step=pm.NUTS())

pm.traceplot(trace, var_names=["sigma_f","l_scale"]);
``````

This code basically stops with an error:

``````~/miniconda/envs/python3/lib/python3.7/site-packages/pymc3/step_methods/hmc/quadpotential.py in raise_ok()
262                               ' is zero.'.format(*name_slc[ii]))
--> 263             raise ValueError('\n'.join(errmsg))
264

ValueError: Mass matrix contains zeros on the diagonal.
The derivative of RV `sigma_f_log__`.ravel() is zero.

The above exception was the direct cause of the following exception:

RuntimeError                              Traceback (most recent call last)
<ipython-input-2-44d3d3d36dd9> in <module>
41
42     obs = pm.MvNormal('obs', mean_func(X_train), K, observed=Y_train)
---> 43     trace = pm.sample(1000,chains=2,cores=2, step=pm.NUTS())
44
45 pm.traceplot(trace, var_names=["sigma_f","l_scale"]);

~/miniconda/envs/python3/lib/python3.7/site-packages/pymc3/sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs)
543         _print_step_hierarchy(step)
544         try:
--> 545             trace = _mp_sample(**sample_args, **parallel_args)
546         except pickle.PickleError:
547             _log.warning("Could not pickle model, sampling singlethreaded.")

~/miniconda/envs/python3/lib/python3.7/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs)
1484         try:
1485             with sampler:
-> 1486                 for draw in sampler:
1487                     trace = traces[draw.chain - chain]
1488                     if trace.supports_sampler_stats and draw.stats is not None:

~/miniconda/envs/python3/lib/python3.7/site-packages/pymc3/parallel_sampling.py in __iter__(self)
490
491         while self._active: