I was checking sklearn's implementation of log marginal likelihood of a Gaussian Process (GP). The implementation is based on Algorithm 2.1 in Rasmussen's Gaussian Processes for Machine Learning which I also attached a snapshot of it for convenience:
However, I constantly came across some cases wherethe log likelihood computed by this formula is positive. One specific example is the following example code:
import numpy as np
from scipy.linalg import cholesky, cho_solve, solve_triangular
from sklearn.gaussian_process.kernels import Matern
kernel=Matern(nu=2.5)
x = np.array([1, 2, 3]).reshape(-1,1)
y = np.array([0.1, 0.2, 0.3])
noise=0
amp2=0.05
K = kernel(x)
cov = amp2 * (K + 0*np.eye(x.shape[0]))
cov[np.diag_indices_from(cov)] += noise
L = cholesky(cov, lower=True)
alpha = cho_solve((L, True), y)
logprob = -0.5 * np.dot(y, alpha) - np.log(np.diag(L)).sum() - \
x.shape[0] / 2. * np.log(2 * np.pi)
print(logprob) # Result: 1.1359631938135135
I believe that the log marginal likelihood of a GP log Pr(y|x, M)
should always be non-positive. Why the code above produces a positive log marginal likelihood?
Thank you!