from abc import abstractmethod
import numpy as np
[docs]
class GaussianProcessRegressorBase:
"""Base class (interface) for Gaussian process regressors."""
@abstractmethod
def __call__(self, x, return_std=False):
"""Compute posterior means (and, optional, standard deviations) of the Gaussian process
at locations `x`.
Parameters
----------
x : array_like
2-D array with locations.
return_std : boolean, optional
If True, also return standard deviations.
Returns
-------
means : ndarray of float
Array with posterior means at locations `x`.
stds : ndarray of float, optional
Array with posterior standard deviations at locations `x`.
"""
raise NotImplementedError
[docs]
@abstractmethod
def fit(self, x, y, noise):
"""Fit the Gaussian process to data `y` observed at training locations `x` with
noise standard deviations `noise`.
Parameters
----------
x : array_like
2-D array with training locations.
y : array_like
1-D array with training data observed at `x`. The length has to be the same as with `x`.
noise : float or array_like
Standard deviations of observational noise. If array_like, the length has to be the
same as with `y`.
"""
raise NotImplementedError
@property
@abstractmethod
def x_train(self):
"""Training locations used.
Returns
-------
x_train : ndarray of float
"""
raise NotImplementedError
@property
@abstractmethod
def y_train(self):
"""Training data used.
Returns
-------
y_train : ndarray of float
"""
raise NotImplementedError
@property
@abstractmethod
def noise(self):
"""Standard deviations of observational noise.
Returns
-------
noise : ndarray of float
"""
raise NotImplementedError
@property
@abstractmethod
def kernel(self):
"""Kernel optimized using training data.
Returns
-------
kernel : instance of a kernel class (depending on the underlying GPR implementation)
"""
raise NotImplementedError
@property
@abstractmethod
def kernel_hyperparameters(self):
"""Hyperparameters of `kernel`.
Returns
-------
kernel_hyperparameters : ndarray of float
"""
raise NotImplementedError
@property
@abstractmethod
def kernel_length_scales(self):
"""Length scales of `kernel`.
Returns
-------
kernel_length_scales : ndarray of float
"""
raise NotImplementedError
[docs]
@abstractmethod
def add_data(self, x_new, y_new, noise):
"""Add new data `y_new` measured at training locations `x_new` to the Gaussian process
and update the posterior distribution. `noise` is again the standard deviation of the noise.
Parameters
----------
x_new : array_like
2-D array with new training locations.
y_new : array_like
1-D array with new training data observed at `x_new`. The length has to be the
same as with `x_new`.
"""
raise NotImplementedError
[docs]
@abstractmethod
def stop_optimization(self):
"""Stop the optimization of kernel hyperparameters. The current `kernel` is now fixed."""
raise NotImplementedError
[docs]
class GaussianProcessRegressorImplementation(GaussianProcessRegressorBase):
"""Base class (interface) for implementations of Gaussian process regressors."""
pass
[docs]
class GaussianProcessRegressor(GaussianProcessRegressorBase):
"""Class for wrapping an implementation of a Gaussian process regressor."""
def __init__(self, gpr_impl: GaussianProcessRegressorImplementation):
self.gpr_impl = gpr_impl
def __call__(self, x, return_std=False):
return self.gpr_impl.__call__(x, return_std)
[docs]
def fit(self, x, y, noise):
assert len(x) == len(y)
self.gpr_impl.fit(x, y, noise)
@property
def x_train(self):
return self.gpr_impl.x_train
@property
def y_train(self):
return self.gpr_impl.y_train
@property
def noise(self):
return self.gpr_impl.noise
@property
def kernel(self):
return self.gpr_impl.kernel
@property
def kernel_hyperparameters(self):
return self.gpr_impl.kernel_hyperparameters
@property
def kernel_length_scales(self):
return self.gpr_impl.kernel_length_scales
[docs]
def add_data(self, x_new, y_new, noise):
assert len(x_new) == len(y_new)
self.gpr_impl.add_data(x_new, y_new, noise)
[docs]
def stop_optimization(self):
self.gpr_impl.stop_optimization()
[docs]
class LogGaussianProcessRegressor(GaussianProcessRegressorBase):
"""Class for log-transforming and wrapping an implementation of a Gaussian process regressor."""
def __init__(self, gpr_impl: GaussianProcessRegressorImplementation):
self.gpr = GaussianProcessRegressor(gpr_impl)
def __call__(self, x, return_std=False):
mean_gpr, std_gpr = self.gpr(x, return_std=True)
# transform mean to log-normal mean
pred = np.exp(mean_gpr + 0.5 * std_gpr**2)
if not return_std:
return pred
# transform standard deviation to log-normal standard deviation
pred_std = np.sqrt((np.exp(std_gpr**2) - 1) * np.exp(2 * mean_gpr + std_gpr**2))
return pred, pred_std
[docs]
def fit(self, x, y, noise):
assert len(x) == len(y)
x = np.asarray(x)
y = np.asarray(y)
if noise is None or np.isscalar(noise) or len(noise) != len(y):
raise Exception("noise for each data point as an array expected")
noise = np.asarray(noise)
# transform noise to context of `gpr` and let it fit log-transformed data
self.gpr.fit(x, np.log(y), noise_from_log(noise, y))
@property
def x_train(self):
return self.gpr.x_train
@property
def y_train(self):
return np.exp(self.gpr.y_train)
@property
def noise(self):
return noise_to_log(self.gpr.noise, self.y_train)
@property
def kernel(self):
return self.gpr.kernel
@property
def kernel_hyperparameters(self):
return self.gpr.kernel_hyperparameters
@property
def kernel_length_scales(self):
return self.gpr.kernel_length_scales
[docs]
def add_data(self, x_new, y_new, noise):
assert len(x_new) == len(y_new) == len(noise)
x_new = np.asarray(x_new)
y_new = np.asarray(y_new)
noise = np.asarray(noise)
self.gpr.add_data(x_new, np.log(y_new), noise_from_log(noise, y_new))
[docs]
def stop_optimization(self):
self.gpr.stop_optimization()
def noise_from_log(noise, y):
"""Transform noise related to a Gaussian process to the noise related to a log-Gaussian process."""
assert len(noise) == len(y)
return np.log(0.5 * (np.sqrt(4 * noise / y**2 + 1) + 1))
def noise_to_log(noise, y):
"""Transform noise related to a log-Gaussian process to the noise related to a Gaussian process."""
assert len(noise) == len(y)
return y**2 * (np.exp(noise) - 1) * np.exp(noise)
