gpr.py
# Copyright 2016 James Hensman, Valentine Svensson, alexggmatthews, fujiisoup
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

from typing import Optional
import tensorflow as tf

from .. import likelihoods
from ..conditionals import base_conditional
from ..kernels import Kernel
from ..logdensities import multivariate_normal
from ..mean_functions import MeanFunction
from .model import GPModel

class GPR(GPModel):
"""
Gaussian Process Regression.

This is a vanilla implementation of GP regression with a Gaussian
likelihood.  Multiple columns of Y are treated independently.

The log likelihood i this models is sometimes referred to as the 'marginal log likelihood', and is given by

.. math::

\\log p(\\mathbf y \\,|\\, \\mathbf f) = \\mathcal N\\left(\\mathbf y\,|\, 0, \\mathbf K + \\sigma_n \\mathbf I\\right)
"""

def __init__(self,
X: tf.Tensor,
Y: tf.Tensor,
kernel: Kernel,
mean_function: Optional[MeanFunction] = None):
"""
X is a data matrix, size [N, D]
Y is a data matrix, size [N, R]
# kernel, mean_function are appropriate GPflow objects
"""
likelihood = likelihoods.Gaussian()

self.X = X
self.Y = Y
super().__init__(kernel, likelihood, mean_function)

def log_likelihood(self):
"""
Construct a tensorflow function to compute the likelihood.

\log p(Y | theta).

"""
K = self.kernel(self.X)
S = tf.eye(self.X.shape,
dtype=self.X.dtype) * self.likelihood.variance
L = tf.linalg.cholesky(K + S)
m = self.mean_function(self.X)
logpdf = multivariate_normal(
self.Y, m,
L)  # (R,) log-likelihoods for each independent dimension of Y
return tf.reduce_sum(logpdf)

def predict_f(self, Xnew, full_cov=False, full_output_cov=False):
"""
Xnew is a data matrix, point at which we want to predict

This method computes

p(F* | Y )

where F* are points on the GP at Xnew, Y are noisy observations at X.
"""
X = self.X
y = self.Y - self.mean_function(X)
Kmn = self.kernel(X, Xnew)
S = tf.eye(X.shape, dtype=X.dtype) * self.likelihood.variance
Kmm = self.kernel(X)
Knn = self.kernel(Xnew, full=full_cov)
f_mean, f_var = base_conditional(
Kmn, Kmm + S, Knn, y, full_cov=full_cov,
white=False)  # [N, P], [N, P] or [P, N, N]
return f_mean + self.mean_function(Xnew), f_var