The GaussianProcess Class¶
Implementation of a Gaussian Process Emulator.
This class provides a representation of a Gaussian Process Emulator. It contains methods for fitting the GP to a given set of hyperparameters, computing the negative log marginal likelihood plus prior (so negative log posterior) and its derivatives, and making predictions on unseen data. Note that routines to estimate hyperparameters are not included in the class definition, and are instead provided externally to facilitate implementation of high performance versions.
The required arguments to initialize a GP is a set of training data consisting of the inputs and targets for each of those inputs. These must be numpy arrays whose first axis is of the same length. Targets must be a 1D array, while inputs can be 2D or 1D. If inputs is 1D, then it is assumed that the length of the second axis is unity.
Optional arguments are the particular mean function to use
(default is zero mean), the covariance kernel to use (default is
the squared exponential covariance), a list of prior distributions
for each hyperparameter (default is no prior information on any
hyperparameters) and the method for handling the nugget parameter.
The nugget is additional “noise” that is added to the diagonal of
the covariance kernel as a variance. This nugget can represent
uncertainty in the target values themselves, or simply be used to
stabilize the numerical inversion of the covariance matrix. The
nugget can be fixed (a non-negative float), can be found
adaptively (by passing the string "adaptive" to make the noise
only as large as necessary to successfully invert the covariance
matrix), can be fit as a hyperparameter (by passing the string
"fit"), or pivoting can be used to ignore any collinear matrix
rows and ensure a zero nugget is used (by passing the string
"pivot").
The internal emulator structure involves arrays for the inputs,
targets, and hyperparameters. Other useful information are the
number of training examples n, the number of input parameters
D, and the number of hyperparameters n_params. These
parameters can be obtained externally by accessing these
attributes.
Example:
>>> import numpy as np
>>> from mogp_emulator import GaussianProcess
>>> x = np.array([[1., 2., 3.], [4., 5., 6.]])
>>> y = np.array([4., 6.])
>>> gp = GaussianProcess(x, y)
>>> print(gp)
Gaussian Process with 2 training examples and 3 input variables
>>> gp.n
2
>>> gp.D
3
>>> gp.n_params
5
>>> gp.fit(np.zeros(gp.n_params))
>>> x_predict = np.array([[2., 3., 4.], [7., 8., 9.]])
>>> gp.predict(x_predict)
(array([4.74687618, 6.84934016]), array([0.01639298, 1.05374973]),
array([[8.91363045e-05, 7.18827798e-01, 3.74439445e-16],
[4.64005897e-06, 3.74191346e-02, 1.94917337e-17]]))
-
class
mogp_emulator.GaussianProcess.GaussianProcess(inputs, targets, mean=None, kernel=<mogp_emulator.Kernel.SquaredExponential object>, priors=None, nugget='adaptive', inputdict={}, use_patsy=True)¶ Implementation of a Gaussian Process Emulator.
This class provides a representation of a Gaussian Process Emulator. It contains methods for fitting the GP to a given set of hyperparameters, computing the negative log marginal likelihood plus prior (so negative log posterior) and its derivatives, and making predictions on unseen data. Note that routines to estimate hyperparameters are not included in the class definition, and are instead provided externally to facilitate implementation of high performance versions.
The required arguments to initialize a GP is a set of training data consisting of the inputs and targets for each of those inputs. These must be numpy arrays whose first axis is of the same length. Targets must be a 1D array, while inputs can be 2D or 1D. If inputs is 1D, then it is assumed that the length of the second axis is unity.
Optional arguments are the particular mean function to use (default is zero mean), the covariance kernel to use (default is the squared exponential covariance), a list of prior distributions for each hyperparameter (default is no prior information on any hyperparameters) and the method for handling the nugget parameter. The nugget is additional “noise” that is added to the diagonal of the covariance kernel as a variance. This nugget can represent uncertainty in the target values themselves, or simply be used to stabilize the numerical inversion of the covariance matrix. The nugget can be fixed (a non-negative float), can be found adaptively (by passing the string
"adaptive"to make the noise only as large as necessary to successfully invert the covariance matrix), can be fit as a hyperparameter (by passing the string"fit"), or pivoting can be used to ignore any collinear matrix rows and ensure a zero nugget is used (by passing the string"pivot").The internal emulator structure involves arrays for the inputs, targets, and hyperparameters. Other useful information are the number of training examples
n, the number of input parametersD, and the number of hyperparametersn_params. These parameters can be obtained externally by accessing these attributes.Example:
>>> import numpy as np >>> from mogp_emulator import GaussianProcess >>> x = np.array([[1., 2., 3.], [4., 5., 6.]]) >>> y = np.array([4., 6.]) >>> gp = GaussianProcess(x, y) >>> print(gp) Gaussian Process with 2 training examples and 3 input variables >>> gp.n 2 >>> gp.D 3 >>> gp.n_params 5 >>> gp.fit(np.zeros(gp.n_params)) >>> x_predict = np.array([[2., 3., 4.], [7., 8., 9.]]) >>> gp.predict(x_predict) (array([4.74687618, 6.84934016]), array([0.01639298, 1.05374973]), array([[8.91363045e-05, 7.18827798e-01, 3.74439445e-16], [4.64005897e-06, 3.74191346e-02, 1.94917337e-17]]))
-
__init__(inputs, targets, mean=None, kernel=<mogp_emulator.Kernel.SquaredExponential object>, priors=None, nugget='adaptive', inputdict={}, use_patsy=True)¶ Create a new GaussianProcess Emulator
Creates a new GaussianProcess Emulator from either the input data and targets to be fit and optionally a mean function, covariance kernel, and nugget parameter/method.
Required arguments are numpy arrays
inputsandtargets, described below. Additional arguments aremeanto specify the mean function (default isNonefor zero mean),kernelto specify the covariance kernel (default is the squared exponential kernel),priorsto indicate prior distributions on the hyperparameters, andnuggetto specify how the nugget parameter is handled (see below; default is to fit the nugget adaptively).inputsis a 2D array-like object holding the input data, whose shape isnbyD, wherenis the number of training examples to be fit andDis the number of input variables to each simulation. Ifinputsis a 1D array, it is assumed thatD = 1.targetsis the target data to be fit by the emulator, also held in an array-like object. This must be a 1D array of lengthn.priormust be a list of lengthn_paramswhose elements are eitherPrior-derived objects orNone. Each element is used as the prior for the corresponding parameter (withNoneindicating an uninformative prior). Passing the empty list orNoneas this argument (in its entirety) may be used as an abbreviation for a list ofn_paramswhere all list elements areNone.nuggetcontrols how additional noise is added to the emulator targets when fitting. This can be specified in several ways. If a string is provided, it can take the values of"adaptive"or"fit", which indicate that the nugget will be chosen in the fitting process. The nugget can also be chosen as part of the fitting process, with three options for nugget handling:"adaptive"means that the nugget will be made only as large as necessary to invert the covariance matrix,"fit"means that the nugget will be treated as a hyperparameter to be optimized, and"pivot"indicates that pivoting will be used to ignore any collinear rows and ensure use of a nugget of zero. Alternatively, a non-negative float can be used to specify a fixed noise level. If no value is specified for thenuggetparameter,"adaptive"is the default.Parameters: - inputs – Numpy array holding emulator input
parameters. Must be 1D with length
nor 2D with shapenbyD, wherenis the number of training examples andDis the number of input parameters for each output. :type inputs: ndarray - targets (ndarray) – Numpy array holding emulator targets. Must be
1D with length
n - mean (None or MeanFunction) – Mean function to be used (optional, default is
Nonefor a zero mean) - kernel (Kernel or str) – Covariance kernel to be used (optional, default
is Squared Exponential) Can provide either a
Kernelobject or a string matching the kernel type to be used. - priors (list or None) – List of priors to be used. Must be None
(default) or an empty list (indicates
uninformative priors) or list of length
n_params. Any parameter for which you wish to specify an uninformative prior, passNone. Number of parameters is the number of parameters in the mean function plusD + 2(one correlation length per input plus a covariance scale and a nugget). If the nugget will not be fit, the list can have lengthn_params - 1. - nugget (float or str) – Noise to be added to the diagonal, specified as
a string or a float. A non-negative float
specifies the noise level explicitly, while a
string indicates that the nugget will be found
via fitting, either as
"adaptive","fit", or"pivot"(see above for a description). Default is"adaptive".
Returns: New
GaussianProcessinstanceReturn type: - inputs – Numpy array holding emulator input
parameters. Must be 1D with length
-
D¶ Returns number of inputs for the emulator
Returns: Number of inputs for the emulator object Return type: int
-
fit(theta)¶ Fits the emulator and sets the parameters
Pre-calculates the matrices needed to compute the log-likelihood and its derivatives and make subsequent predictions. This is called any time the hyperparameter values are changed in order to ensure that all the information is needed to evaluate the log-likelihood and its derivatives, which are needed when fitting the optimal hyperparameters.
The method computes the mean function and covariance matrix and inverts the covariance matrix using the method specified by the value of
nugget_type. The factorized matrix and the product of the inverse with the difference between the targets and the mean are cached for later use, and the negative marginal log-likelihood is also cached. This method has no return value, but it does modify the state of the object.Parameters: theta (ndarray) – Values of the hyperparameters to use in fitting. Must be a numpy array with length n_paramsReturns: None
-
get_K_matrix()¶ Returns current value of the covariance matrix as a numpy array. Does not include the nugget parameter, as this is dependent on how the nugget is fit.
-
inputs¶ Returns inputs for the emulator as a numpy array
Returns: Emulator inputs, 2D array with shape (n, D)Return type: ndarray
-
logpost_deriv(theta)¶ Calculate the partial derivatives of the negative log-posterior
Calculate the partial derivatives of the negative log-posterior with respect to the hyperparameters. Note that this function is normally used only when fitting the hyperparameters, and it is not needed to make predictions.
During normal use, the
logpost_derivmethod is called after evaluating thelogposteriormethod. The implementation takes advantage of this by reusing cached results, as the factorized covariance matrix is expensive to compute and is used by thelogposterior,logpost_deriv, andlogpost_hessianmethods. If the function is evaluated with a different set of parameters than was previously used to set the log-posterior, the method callsfit(and subsequently resets the cached information).Parameters: theta (ndarray) – Value of the hyperparameters. Must be array-like with shape (n_params,)Returns: partial derivatives of the negative log-posterior with respect to the hyperparameters (array with shape (n_params,))Return type: ndarray
-
logpost_hessian(theta)¶ Calculate the Hessian of the negative log-posterior
Calculate the Hessian of the negative log-posterior with respect to the hyperparameters. Note that this function is normally used only when fitting the hyperparameters, and it is not needed to make predictions. It is also used to estimate an appropriate step size when fitting hyperparameters using the lognormal approximation or MCMC sampling.
When used in an optimization routine, the
logpost_hessianmethod is called after evaluating thelogposteriormethod. The implementation takes advantage of this by storing the inverse of the covariance matrix, which is expensive to compute and is used by thelogposteriorandlogpost_derivmethods as well. If the function is evaluated with a different set of parameters than was previously used to set the log-posterior, the method callsfitto compute the needed information and changes the cached values.Parameters: theta (ndarray) – Value of the hyperparameters. Must be array-like with shape (n_params,)Returns: Hessian of the negative log-posterior (array with shape (n_params, n_params))Return type: ndarray
-
logposterior(theta)¶ Calculate the negative log-posterior at a particular value of the hyperparameters
Calculate the negative log-posterior for the given set of parameters. Calling this method sets the parameter values and computes the needed inverse matrices in order to evaluate the log-posterior and its derivatives. In addition to returning the log-posterior value, it stores the current value of the hyperparameters and log-posterior in attributes of the object.
Parameters: theta (ndarray) – Value of the hyperparameters. Must be array-like with shape (n_params,)Returns: negative log-posterior Return type: float
-
n¶ Returns number of training examples for the emulator
Returns: Number of training examples for the emulator object Return type: int
-
n_params¶ Returns number of hyperparameters
Returns the number of hyperparameters for the emulator. The number depends on the choice of mean function, covariance function, and nugget strategy, and possibly the number of inputs for certain choices of the mean function.
Returns: Number of hyperparameters Return type: int
-
nugget¶ Returns emulator nugget parameter
Returns current value of the nugget parameter. If the nugget is to be selected adaptively, by pivoting, or by fitting the emulator and the nugget has not been fit, returns
None.Returns: Current nugget value, either a float or NoneReturn type: float or None The
nuggetparameter controls how noise is added to the covariance matrix in order to stabilize the inversion or smooth the emulator predictions. Ifnuggetis a non-negative float, then that particular value is used for the nugget. Note that setting this parameter to be zero enforces that the emulator strictly interpolates between points. Alternatively, a string can be provided. A value of"fit"means that the nugget is treated as a hyperparameter, and is the last entry in thethetaarray. Alternatively, ifnuggetis set to be"adaptive", the fitting routine will adaptively make the noise parameter as large as is needed to ensure that the emulator can be fit. Finally, pivoting can be selected by setting to"pivot", which will ignore any collinear rows in the covariance matrix.Internally, this modifies both the way the nugget is chosen (which can be determined via the
nugget_typeproperty) and the value itself (thenuggetproperty)Parameters: nugget (float or str) – Noise to be added to the diagonal, specified as a string or a float. A non-negative float specifies the noise level explicitly, while a string indicates that the nugget will be found via fitting, either as "adaptive","pivot", or"fit"(see above for a description).Returns: None Return type: None
-
nugget_type¶ Returns method used to select nugget parameter
Returns a string indicating how the nugget parameter is treated, either
"adaptive","pivot","fit", or"fixed". This is automatically set when changing thenuggetproperty.Returns: Current nugget fitting method Return type: str
-
predict(testing, unc=True, deriv=True, include_nugget=True)¶ Make a prediction for a set of input vectors for a single set of hyperparameters
Makes predictions for the emulator on a given set of input vectors. The input vectors must be passed as a
(n_predict, D),(n_predict,)or(D,)shaped array-like object, wheren_predictis the number of different prediction points under consideration andDis the number of inputs to the emulator. If the prediction inputs array is 1D andD == 1for the GP instance, then the 1D array must have shape(n_predict,). Otherwise, if the array is 1D it must have shape(D,), and the method assumesn_predict == 1. The prediction is returned as an(n_predict, )shaped numpy array as the first return value from the method.Optionally, the emulator can also calculate the variances in the predictions and the derivatives with respect to each input parameter. If the uncertainties are computed, they are returned as the second output from the method as an
(n_predict,)shaped numpy array. If the derivatives are computed, they are returned as the third output from the method as an(n_predict, D)shaped numpy array.The final input to the method determines if the predictive variance should include the nugget or not. For situations where the nugget represents observational error and predictions are estimating the true underlying function, this should be set to
False. However, most other cases should probably include the nugget, as the emulator is using it to represent some of the uncertainty in the underlying simulator, so the default value isTrue.Parameters: - testing (ndarray) – Array-like object holding the points where
predictions will be made. Must have shape
(n_predict, D)or(D,)(for a single prediction) - unc (bool) – (optional) Flag indicating if the uncertainties
are to be computed. If
Falsethe method returnsNonein place of the uncertainty array. Default value isTrue. - deriv (bool) – (optional) Flag indicating if the derivatives
are to be computed. If
Falsethe method returnsNonein place of the derivative array. Default value isTrue. - include_nugget (bool) – (optional) Flag indicating if the
nugget should be included in the
predictive variance. Only relevant if
unc = True. Default isTrue.
Returns: Tuple of numpy arrays holding the predictions, uncertainties, and derivatives, respectively. Predictions and uncertainties have shape
(n_predict,)while the derivatives have shape(n_predict, D). If theuncorderivflags are set toFalse, then those arrays are replaced byNone.Return type: tuple
- testing (ndarray) – Array-like object holding the points where
predictions will be made. Must have shape
-
priors¶ The current list priors used in computing the log posterior
To set the priors, must be a list or
None. Entries can beNoneor a subclass ofPrior.Noneindicates weak prior information. An empty list orNonemeans all uninformative priors. Otherwise list should have the same length as the number of hyperparameters, or alternatively can be one shorter than the number of hyperparameters ifnugget_typeis"adaptive","pivot"or"fixed"meaning that the nugget hyperparameter is not fit but is instead fixed or found adaptively. If the nugget hyperparameter is not fit, the prior for the nugget will automatically be set toNoneeven if a distribution is provided.
-
targets¶ Returns targets for the emulator as a numpy array
Returns: Emulator targets, 1D array with shape (n,)Return type: ndarray
-
theta¶ Returns emulator hyperparameters
Returns current hyperparameters for the emulator as a numpy array if they have been fit. If no parameters have been fit, returns
None. Note that the number of parameters depends on the mean function, so the length of this array will vary across instances.Returns: Current parameter values (numpy array of length n_params), orNoneif the parameters have not been fit.Return type: ndarray or None When set, pre-calculates the matrices needed to compute the log-likelihood and its derivatives and make subsequent predictions. This is called any time the hyperparameter values are changed in order to ensure that all the information is needed to evaluate the log-likelihood and its derivatives, which are needed when fitting the optimal hyperparameters.
The method computes the mean function and covariance matrix and inverts the covariance matrix using the method specified by the value of
nugget. The factorized matrix and the product of the inverse with the difference between the targets and the mean are cached for later use, and the negative marginal log-likelihood is also cached. This method has no return value, but it does modify the state of the object.Parameters: theta (ndarray) – Values of the hyperparameters to use in fitting. Must be a numpy array with length n_params
-
The PredictResult Class¶
-
class
mogp_emulator.GaussianProcess.PredictResult¶ Prediction results object
Dictionary-like object containing mean, uncertainty (variance), and derivatives with respect to the inputs of an emulator prediction. Values can be accessed like a dictionary with keys
'mean','unc', and'deriv'(or indices 0, 1, and 2 for the mean, uncertainty, and derivative for backwards compatability), or using attributes (p.meanifpis an instance ofPredictResult). Also supports iteration and unpacking with the ordering(mean, unc, deriv)to be consistent with indexing behavior.Code is mostly based on scipy’s
OptimizeResultclass, with some additional code to support iteration and integer indexing.Variables: - mean – Predicted mean for each input point. Numpy array with
shape
(n_predict,) - unc – Predicted variance for each input point. Numpy array
with shape
(n_predict,) - deriv – Predicted derivative with respect to the inputs for
each input point. Numpy array with shape
(n_predict, D)
- mean – Predicted mean for each input point. Numpy array with
shape