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))

This results in the fitted function as the image.

enter image description here

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?

  • 2
    $\begingroup$ 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). $\endgroup$
    – adityar
    Oct 18, 2018 at 15:10
  • $\begingroup$ @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? $\endgroup$
    – Luca
    Oct 23, 2018 at 7:49
  • 1
    $\begingroup$ 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. $\endgroup$
    – adityar
    Oct 23, 2018 at 8:37


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.