Gaussian Process Class

class sam.GaussianProcess(x, y, yErr=None, kernel='squaredExp')

A class for doing Gaussian process computation and modeling.

Initializes the Gaussian process (GP) object with the observations and a kernel type.

Parameters:

• x – The locations of the observations to condition the GP on; should have shape [nSamples x nDimensions].
• y – The values of the input at the locations given by x. Should have shape [nSamples].
• kernel – The name of the kernel to use. Must be one of: Exponential, Squared Exponential, Matern(3/2), or Matern(5/2)

Returns:

The Gaussian processes object for further use.

draw(self, xTest, Size nDraws=1)

Draws a random vector from the Gaussian process at the specified test points.

Parameters:

• xTest – the points to make predictions at. Should be [nPoints x nDimensions].
• nDraws – the number of draws to produce. The first draw is much more computationally than subsequent draws.

Returns:

A matrix of values sampled from the Gaussian process. [nDraws x nPoints]

gradient(self, xTest)

Computes the expected gradient of the Gaussian process at a given point.

Each component i is dy/dx_i.

Parameters:xTest – the point to measure the gradient at. Should be [nDimensions].

Returns

The components of the gradient as a vector of length [nDimensions].

logLikelihood(self, __Pyx_memviewslice params=None) → double

Computes the log-likelihood of the Gaussian process for the given parameters, x, and y.

Parameters:

• params (optional) – the parameters of the kernel to use. Otherwise use whatever is
• in self.params. (currently) –

Returns:

The log-likelihood as a double.

optimizeParams(self, paramsGuess=None, tol=1e-3, logBounds=[(-10, 10), (-5, 5), (-10, 10)], fixWhiteNoise=None)

Attempts to locate the maximum likelihood parameters for the given x and y.

Parameters:

• params (optional) – the kernel parameters to start from. Otherwise use whatever is currently in self.params.
• tol (optional) – the error allowed by the optimizer before it is considered to be converged.
• logBounds (optional) – the bounds on the log of the parameters to be optimized. This is important for avoiding implausible parts of parameter space which would cause positive-definite errors in the Cholesky computation.
• whiteNoise (optional) – a value to fix the white noise at. If none, fits the white noise along with the other parameters.

Returns:

The optimized parameters (which are also written to self.params).

precompute(self, __Pyx_memviewslice params=None, bool force=False) → int

Conducts essential precomputation for the Gaussian Process.

Specifically, this constructs the covariance matrix, the Cholesky thereof (L), and alpha = L.T * L * y. For certain inputs, the numerics are not very stable, and so the Cholesky operation may incorrectly determine that the covariance matrix is not positive definite. This is best resolved by using a small amount of white noise (a diagonal addition to the covariance matrix); see parameter information.

Parameters:

• params (optional) – the parameters of the kernel to use. Otherwise use whatever is currently in self.params.
• force (default False) – Recompute even if cholesky and alpha are flagged as fresh. This shouldn’t be necessary under normal conditions, since it is automatically set to True if params is not None.

Returns:

0 if successful, otherwise -1.

predict(self, xTest, bool getCovariance=True)

Compute the Gaussian process’ prediction for a given set of points.

Parameters:

• xTest – the points to make predictions at. Should be [nPoints x nDimensions].
• getCovariance – Whether to get the predictive covariance as well as the mean.

Returns:

A vector of predicted means and (if getCovariance==True), the covariance matrix of the predictions. For reference, the marginal uncertainties are np.sqrt(np.diag(var)).

setY(self, newY)

Changes the measured y values in a way that minimizes the required recomputation.

Specifically, this will need to recompute alpha O(n^2) but not the Cholesky O(n^3)

Parameters:newY – the new values of y to use. Should be [nPoints]. Returns:None

Example Use

import numpy as np
from matplotlib import pyplot as plt
from sam import GaussianProcess

x = np.linspace(0,10,10)
y = np.sin(x)

f = GaussianProcess(x,y)
f.optimizeParams(5*ones(3))

xTest = np.linspace(0,10,100)
yTest, yErr = f.predict(xTest)
yErr = np.sqrt(np.diag(yErr))
plt.plot(xTest,yTest)
plt.fill_between(xTest,yTest-yErr,yTest+yErr,alpha=.5)
plt.plot(x,y,'.')