# Adding positivity constraint to GP

I am fitting a simple one-dimensional GP model to my observations. My observations are simply the mean number of people/minute who visited a store between different time intervals.

So, my data looks as follows:

import numpy as np

# Observations
y = np.asarray([3.7, 7.2, 30.1, 20.3, 10.2, 18.1, 15.2, 30.5, 25.0, 11.1])
# Independent variable
X = np.asarray([2, 4, 7, 9, 11, 13, 15, 17, 19, 23])[:, np.newaxis]


I then fit the scikit-learn GP model to the data.

kernel = 1.0 * RBF(length_scale=100.0, length_scale_bounds=(1e-2, 1e3)) \
+ WhiteKernel(noise_level=1, noise_level_bounds=(1e-10, 1e+1))
gp = GaussianProcessRegressor(kernel=kernel,
alpha=0.0).fit(X, y)
X_ = np.linspace(0, 24, 500)
y_mean, y_cov = gp.predict(X_[:, np.newaxis], return_cov=True)
plt.plot(X_, y_mean, 'k', lw=3, zorder=9)

plt.scatter(X[:, 0], y, c='r', s=50, zorder=10, edgecolors=(0, 0, 0))
plt.tight_layout()
plt.show()


This results in the fitted function as the image.

Now, my output variable should always be >=0. I am guessing there is nothing in this model that guarantees that and I was wondering if it is possible to add that constraint somehow to the GP fitting procedure. Although, I am guessing that because of the fact that Gaussians are being used it is not possible to add this sort of constraint?

• You could log-transform your y-variable. If you need a MAP solution, you could either use monte-carlo simulations to find it, or, noting that the log is a monotonic transformation, you could find the analytical expression of the transformed variable's density and optimize it (see this for reference). Oct 18, 2018 at 15:10
• @InfProbSciX Thank you for the comment. I did log-transform the variable and fitted the GP to it. I then applied the exponential to the mean and the drawn functions from it to transform it back to the original function. Is that a valid approach?
– Luca
Oct 23, 2018 at 7:49
• To find the MAP, you need to do a bit of work; the MAP path of the not-logged function, is not equal to the $exp$ of the MAP of the GP fitted on logged values (see the Jensen inequality for more info). If you want simulations of the GP, then applying the $exp$ on draws from the GP fitted on logged values is a fair way to go. To obtain the MAP however, you'd need to define a function like $gp.log\_marginal\_likelihood() + ln J$ where the extra term is the log jacobian of the transform. You need to find the path that maximises this function. Oct 23, 2018 at 8:37