Gaussian processes

gpyreg.gaussian_process

This module contains the GP class, which is the main entry point for working with Gaussian processes. A class instance is initialized with specified dimension and covariance/mean/noise functions. Then, its various methods can be used to fit data, make predictions, plot the GP, etc.

Reproducible fitting

Since version 1.1.0, fitting accepts an rng argument. Given a GP instance gp and training arrays X and y, an integer seed initializes the stream used for both the space-filling design and hyperparameter sampling:

hyp, optimization, sampling = gp.fit(X, y, rng=1234)

Pass a numpy.random.Generator to share a stream across calls. See Random number generation for stream ownership and the default behavior.

Reusing prediction kernels

Since version 1.2.0, gpyreg.GP.predict() can return the latent kernel matrices between training and prediction inputs. Given a fitted gp and an (M, D) array X_star:

mu, variance = gp.predict(X_star)
mu, variance, kernels = gp.predict(
    X_star, return_cross_covariance=True
)

mu and variance have shape (M, 1) when predictions are averaged across hyperparameter samples, or (M, S) with separate_samples=True. The extra kernels tuple always contains one entry for each of the S hyperparameter samples, independently of separate_samples. Each entry is K(X, X_star) with shape (N, M), where N is the training-set size. A GP without training targets returns None for each sample.

These are unconditioned latent kernels. Observation noise is excluded even with add_noise=True. To obtain the posterior predictive covariance between prediction points, use gpyreg.GP.predict_full().

With return_lpd=True, the kernel tuple follows the log predictive density:

mu, variance, lpd, kernels = gp.predict(
    X_star, y_star=y_star, return_lpd=True,
    return_cross_covariance=True
)

Here y_star contains the observed values at the prediction points. Treat returned kernels as read-only. Retaining all matrices requires approximately 8 * N * M * S bytes for float64 data, in addition to other prediction allocations. GPyReg imposes no retention cap; request these matrices only when needed and release them after use. The option defaults to False.

GP

class gpyreg.GP(D: int, covariance: object, mean: object, noise: object, raise_on_cholesky_failure: bool = False)

A single Gaussian Process (GP).

Parameters:

Dint

The dimension of the Gaussian Process.

covarianceobject

The covariance function to use. This can be one of the objects from the following module: gpyreg.covariance_functions.

meanobject

The mean function to use. This can be one of the objects from the following module: gpyreg.mean_functions.

noiseobject

The noise function to use. This can be one of the objects from the following module: gpyreg.noise_functions.

raise_on_cholesky_failurebool, defaults to False

What a failed Cholesky factorization of the training covariance does, in the posteriors and in the objective of fit() (log_likelihood() and log_posterior()). By default the factorization is tried again with the noise variance multiplied tenfold, up to ten attempts, and a posterior keeps the multiplier of the attempt that succeeded (sn2_mult of Posterior), which the predictions apply to the noise, while get_hyperparameters() returns the noise as fitted. With True the first failure raises numpy.linalg.LinAlgError, as in MATLAB BADS (CholAttempts = 0), and no posterior carries a multiplier other than one. In fit(), the space-filling design ranks a starting point whose objective raises it last, as one of infinite value, and starts no optimization from it unless every point of the design raised, while a failure in the optimization, in the sampling or in the posteriors of the fitted hyperparameters raises from fit, which leaves the GP as it was. The GP keeps the choice as its attribute raise_on_cholesky_failure, which a copy or a pickle of the GP keeps; a GP unpickled from gpyreg 1.3.3 or earlier has no such attribute and inflates the noise.

Raises:

ValueError

Raised when the covariance function has periods (its attribute periods) and their number is not D.

bounds_to_dict(lower_bounds: ndarray, upper_bounds: ndarray)

Convert the given hyperparameter lower and upper bounds to a dict.

Parameters:

lower_boundsndarray, shape (hyp_N,)

The lower bounds.

upper_boundsndarray, shape (hyp_N,)

The upper bounds.

Returns:

bounds_dictdict

A dictionary of the current hyperparameter names and tuples of their lower and upper bounds.

clean()

Clean auxiliary computational structures from the Gaussian Process, thus reducing memory usage. These can be reconstructed with a call to update() with compute_posterior=True.

Furthermore, the temporary_data attribute is being cleared.

fit(X: ndarray = None, y: ndarray = None, s2: ndarray = None, hyp0=None, options: dict = None, rng=None)

Train the hyperparameters of the Gaussian Process.

A fit that raises leaves the GP as it was before the call, with the training data, bounds, priors and posteriors that it held.

Parameters:

Xndarray, shape (N, D), optional

Training points that will replace the current training points of the GP. If not given the current training points are used.

yndarray, shape (N, 1), optional

Training targets that will replace the current training targets of the GP. If not given the current training targets are used.

s2ndarray, shape (N, 1), optional

Noise variance at training points that will replace the current noise variances of the GP. If not given the current noise variances are used.

optionsdict, optional

A dictionary of options for training. The counts among them, opts_N, init_N, n_samples, thin and burn, are whole numbers of an integer or a float type (2.0 is taken as 2), or 0-d arrays that hold one. The possible options are:

opts_Nint, defaults to 3

Number of hyperparameter optimization runs.

init_Nint, defaults to 1024

Initial design size for hyperparameter optimization.

df_baseint, defaults to 7

The degrees of freedom of a "student_t" or "smoothbox_student_t" prior whose df is NaN, filled into a copy of the priors for the duration of the fit, as gplite_train.m fills its local copy. The GP keeps the priors as they were set, and outside the fit a NaN df reads as "gaussian" or "smoothbox" (see set_priors()).

n_samplesint, defaults to 10

Number of hyperparameters to sample.

thinint, defaults to 5

Thinning parameter for slice sampling: one sample in thin is kept. A whole number greater than zero.

burnint or None, defaults to thin * n_samples

Burn parameter for slice sampling: the number of samples drawn and dropped before the first one kept. A whole number of at least zero, or None, which leaves it to SliceSampler.sample().

lower_boundsstr or ndarray, defaults to “current”

User-provided lower bounds. Any values which are nan will be filled with the recommended bounds. “recommended” means use all recommended bounds. “current” means use current bounds.

upper_boundsstr or ndarray, defaults to “current”

User-provided upper bounds. Any values which are nan will be filled with the recommended bounds. “recommended” means use all recommended bounds. “current” means use current bounds.

init_method{‘sobol’, ‘rand’}, defaults to ‘sobol’

Specify whether to use Sobol or random sequences for the initial space-filling design.

sampler_name{‘slicesample’}, defaults to ‘slicesample’

The name of the sampler to use. Currently only slice sampling is supported.

tol_optfloat, defaults to 1e-5

Optimization tolerance for stopping.

tol_opt_mcmcfloat, defaults to 1e-3

Preliminary optimization tolerance when doing MCMC.

widthsndarray, shape (hyp_n,), optional

Default widths to use for sampling. If not provided appropriate ones will be computed.

rngNone, numpy.random.Generator or seed, optional

Where the fit’s random draws come from (the space-filling initial design and the slice sampler). None (default) keeps NumPy’s global legacy stream, as before generators were supported, so np.random.seed still fixes a fit; a numpy.random.Generator is used as is (and shared with the caller); an integer or SeedSequence seeds a new generator. See gpyreg.rng.resolve_rng().

Returns:

hypndarray, shape (hyp_samples, hyp_N)

The fitted hyperparameters.

optimize_resultOptimizeResult

The optimization result represented as a OptimizeResult object. For more details see scipy.optimize.minimize().

sampling_resultdict

If sampling was performed this is a dictionary with info on the sampling run, and None otherwise.

Raises:

ValueError

Raised when the GP has no training data, neither given to the fit nor held from an earlier fit or update, before the fit changes anything.

ValueError

Raised by the check of the shapes of the training data given, before the fit changes anything: when X is neither a 1-D nor a 2-D array or does not have D columns, when y does not hold one value per row of X, or when an array s2 does not hold one variance per row of X.

TypeError

Raised by the same check, before the fit changes anything, when s2 is neither an array, a number nor None, such as a list.

ValueError

Raised, before the fit changes anything, when the data it would hold, those given with those the GP holds for what is not given, would make the number of targets, or of the noise variances that the noise function reads (as gpyreg.noise_functions.GaussianNoise with user_provided_add does), differ from the number of inputs: when X is given without y, or without s2, and the targets, or those variances, that the GP holds, one per input it holds, are not one per row of X.

ValueError

Raised by get_recommended_bounds(), through which the fit fills the bounds that are not set, before the fit changes anything: when the option lower_bounds or upper_bounds is neither "recommended", None, "current" nor an array, when a lower bound given is above its upper bound, or when a column of the training inputs has no spread and a hyperparameter whose recommended bounds take their scale from it is not given a finite lower bound.

ValueError

Raised when the option opts_N, init_N or n_samples is not a whole number of at least zero, the option thin not a whole number of at least one, or the option burn neither None nor a whole number of at least zero: a fraction, a number below that, an infinity, NaN, a bool or a value that is not a number, before the fit changes anything.

ValueError

Raised when n_samples is positive and sampler_name is not 'slicesample', after the optimization.

ValueError

Raised by gpyreg.slice_sample.SliceSampler when n_samples is positive, after the optimization: when the widths of the sampler, from the option widths or computed by the fit, are not all positive and finite, or when the log posterior is not finite at the optimized hyperparameters, where the chain starts.

ValueError

Raised by update(), which computes the posterior of the fitted hyperparameters, when one of them is NaN.

LinAlgError

Raised when the Cholesky decomposition of the training covariance fails even after its noise is multiplied tenfold, up to ten times: by update() at the fitted hyperparameters, or by the objective at a starting point or during the optimization, which lets it propagate. Where the GP raises on a failed factorization (raise_on_cholesky_failure) the first failure raises, by update(), during the optimization or during the sampling, and a starting point of the space-filling design whose factorization fails is ranked last instead, and starts no optimization unless every point failed. Whether a covariance that holds NaN, as from a starting point that does, fails there or gives NaN depends on the LAPACK build.

get_bounds()

Gets a dictionary with the current lower and upper bounds of the hyperparameters.

Returns:

bounds_dictdict

A dictionary of the current hyperparameter names and their bounds.

get_hyperparameters(as_array: bool = False)

Get the current hyperparameters of the GP.

If hyperparameters have not been set yet, the result will be filled with NaN.

Parameters:

as_arraybool, defaults to False

Whether to return the hyperparameters as an array of shape (hyp_samples, hyp_N), or a list of dictionaries for each sample.

Returns:

hypobject

The hyperparameters in the form specified by as_array.

get_priors()

Return the current hyperparameter priors as a dict.

Returns:

hyper_priorsdict

A dictionary of the current hyperparameter names and their priors, in the form set_priors() takes, which set_priors writes back as the GP holds them, so that set_priors(get_priors()) changes nothing. A hyperparameter has None, as set_priors takes it for no prior, where all of its prior’s entries are NaN. Any other prior holds the arrays of its entries, NaN included. A block with a prior in some coordinate comes back under the name of the family it was set with, or of the matching Gaussian family where a Student’s t family was set with zero degrees of freedom throughout. A block with no prior in any coordinate, whose location and sigma are NaN throughout, comes back under the Gaussian or the Student’s t family that its degrees of freedom name: "gaussian" where it was set with "gaussian" or "smoothbox", or with "student_t" or "smoothbox_student_t" and zero degrees of freedom throughout; None where it was set with one of the latter two and NaN degrees of freedom throughout; and "student_t" where it was set with one of them and other degrees of freedom.

Raises:

ValueError

Raised when the priors that the GP holds are not in a form set_priors takes, as priors written into hyper_priors directly, or by a version of gpyreg that took them, may not be: a coordinate that has a location (mu, or the ends of a smooth box) and a sigma that is not finite and positive, or the reverse; an inverted smooth box; or a block that holds both a mu and the ends of a smooth box.

get_recommended_bounds(lower_bounds=None, upper_bounds=None)

Return the recommended hyperparameter lower and upper bounds as a dict.

Parameters:

lower_boundsndarray, optional

If present, override the recommended lower bounds with these values. Any nan values will be replaced by the corresponding recommended bounds. Defaults to all nan values.

upper_boundsndarray, optional

If present, override the recommended upper bounds with these values. Any nan values will be replaced by the corresponding recommended bounds. Defaults to all nan values.

Returns:

bounds_dictdict

A dictionary of the hyperparameter names and tuples of lower upper bounds.

Raises:

ValueError

Raise when GP does not have X or y set yet, or when provided bounds are not one of “recommended”/None, “current”, or array_like.

ValueError

Raised when a lower bound given is above the upper bound given for the same hyperparameter. A pair that comes out inverted once the recommendations fill its NaN entries is collapsed onto its lower bound instead.

ValueError

Raised when a column of the training inputs has no spread, as every column of a single training point has none, and a hyperparameter whose recommended bounds take their scale from the spread of that column (a length scale of the kernel, or the scale of gpyreg.mean_functions.NegativeQuadratic) is not given a finite lower bound: its recommended bounds are both -inf, and an upper bound left unset collapses onto the lower bound given.

hyperparameters_from_dict(hyp_dict_list)

Convert a list of hyperparameter dictionaries to a hyperparameter array.

Parameters:

hyp_dict_listobject

A list of hyperparameter dictionaries with hyperparameter names and values. One can also pass just one dictionary instead of a list with one element.

Returns:

hyp_arrndarray, shape (hyp_samples, hyp_N)

The hyperparameter array where hyp_samples is the length of the list hyp_dict_list.

hyperparameters_to_dict(hyp_arr: ndarray)

Convert a hyperparameter array to a list which contains a dictionary with hyperparameter names and values for each hyperparameter sample.

Parameters:

hyp_arrndarray

An array of shape (hyp_samples, hyp_N) or shape (hyp_N,), which is interpreted as shape (1, hyp_N), containing hyperparameters.

Returns:

hyp_dictobject

A list which contains a dictonary with hyperparameter names and values for each sample.

Raises:

ValueError

Raised when the input hyperparameter array has the wrong shape.

log_likelihood(hyp: object, compute_grad: bool = False)

Compute the (positive) log marginal likelihood of the GP for given hyperparameters.

Parameters:

hypobject

Either an 1D array or a dictionary of hyperparameters.

compute_gradbool, defaults to False

Whether to compute the gradient with respect to hyperparameters.

Returns:

lZfloat

The positive log marginal likelihood.

dlZndarray, shape (hyp_N,), optional

The gradient with respect to hyperparameters.

Raises:

LinAlgError

Raised when the Cholesky decomposition failed multiple times even by adding numerical stability values to the matrix, or once where the GP raises on a failed factorization (raise_on_cholesky_failure).

log_posterior(hyp: object, compute_grad: bool = False)

Compute the (positive) log marginal likelihood of the GP with added log prior for given hyperparameters (that is, the unnormalized log posterior).

Each hyperparameter’s prior is renormalized over its bounds, so the value carries a constant that a hand-computed sum of the log marginal likelihood and the prior densities does not. The constant does not depend on the hyperparameters, and is zero for a hyperparameter with no bounds.

Parameters:

hypobject

Either an 1D array or a dictionary of hyperparameters.

compute_gradbool, defaults to False

Whether to compute the gradient with respect to hyperparameters.

Returns:

lZ_plus_posteriorfloat

The positive log marginal likelihood with added log prior.

dlZ_plus_d_posteriorndarray, shape (hyp_N,), optional

The gradient with respect to hyperparameters.

Raises:

LinAlgError

Raised when the Cholesky decomposition failed multiple times even by adding numerical stability values to the matrix, or once where the GP raises on a failed factorization (raise_on_cholesky_failure).

plot(x0: ndarray = None, lb: ndarray = None, ub: ndarray = None, delta_y: float = None, max_min_flag: bool = True)

Plot the Gaussian Process profile centered around a given point.

The plot is a D-by-D panel matrix, in which panels on the diagonal show the profile of the Gaussian Process prediction (mean and +/- 1 SD) by varying one dimension at a time, whereas off-diagonal panels show 2-D contour plots of the GP mean and standard deviation (respectively, above and below diagonal). In each panel, black lines indicate the location of the reference point.

Parameters:

x0ndarray, shape (D,), optional

The reference point.

lbndarray, shape (D,), optional

Lower bounds for the plotting.

ubndarray, shape (D,), optional

Upper bounds for the plotting.

delta_yfloat, optional

Range of the plot such that the plotted predictive GP mean approximately brackets [y0-delta_y, y0+delta_y] where y0 is the predictive GP mean at x0. If lower or upper bounds are given this will do nothing.

max_min_flagbool, defaults to True

If set to False then the minimum, and if set to True then the maximum of the GP training input is used as the reference point.

predict(x_star: ndarray, y_star: ndarray = None, s2_star: ndarray = None, add_noise: bool = False, separate_samples: bool = False, return_lpd: bool = False, *, return_cross_covariance: bool = False)

Compute the GP posterior mean and variance at given points.

Parameters:

x_starndarray, shape (M, D)

The points we want to predict the values at.

y_starndarray, shape (M, 1), optional

True values at the points.

s2_starndarray, shape (M, 1), optional

Noise variance at the points.

add_noisebool, defaults to False

Whether to add noise to the prediction results.

separate_samplesbool, defaults to False

Whether to return the results separately for each hyperparameter sample or averaged.

return_lpdbool, defaults to False

Whether to return the log predictive density at the input points. The density always carries the observation noise, whichever variance add_noise selects for s2. With separate_samples False it is the log density of a Gaussian whose mean is the mean of the per-sample means and whose variance is the mean of the per-sample predictive variances plus the sample variance (ddof=1) of the per-sample means, as gplite_pred.m pools them, and not the average of the per-sample log densities.

return_cross_covariancebool, defaults to False

Whether to append the latent training-to-prediction kernel matrices to the return values. The matrices are kept separate for each hyperparameter sample even when separate_samples is False.

Returns:

mundarray

Posterior mean at the requested points. If we requested separate samples the shape is (M, sample_N) while otherwise it is (M, 1).

s2ndarray

Variance at each point: the latent posterior variance, or, with add_noise, that variance plus the observation noise. Pooled over several hyperparameter samples, it is the mean of the per-sample variances plus the sample variance (ddof=1) of the per-sample means. If we requested separate samples the shape is (M, sample_N) while otherwise it is (M, 1).

lpdndarray, optional

Log predictive density at each point. Returned when return_lpd is True and shaped like mu.

cross_covariancetuple of ndarray or None, optional

Returned when return_cross_covariance is True. Entry s is the latent, unconditioned kernel matrix K(X, x_star) with shape (N, M) for hyperparameter sample s. Prior-only GPs, which have no training targets, return None for every sample. The matrices are call-local and should be treated as read-only; custom covariance results may be copied to give each tuple entry stable values. Observation-noise inputs and add_noise do not affect these latent kernel matrices.

Raises:

ValueError

Raised by the check of the shapes of the data given: when x_star is neither a 1-D nor a 2-D array or does not have D columns, when y_star does not hold one value per row of x_star, or when an array s2_star does not hold one variance per row of x_star.

TypeError

Raised when s2_star is neither an array, a number nor None, such as a list.

ValueError

Raised when return_lpd is True and y_star is None.

LinAlgError

Raised when a posterior holds the negative inverse of the training covariance without its Cholesky factor, as one pickled by gpyreg 1.3.1 or earlier does where the smallest noise variance at the training inputs is below 1e-6, and the factorization that computes the factor again fails even after its noise is multiplied tenfold, up to ten times, or fails at all where the GP raises on a failed factorization (raise_on_cholesky_failure).

predict_full(x_star: ndarray, y_star: ndarray = None, s2_star: ndarray = None, add_noise: bool = False)

Compute the GP posterior mean and full covariance matrix for each hyperparameter sample.

Parameters:

x_starndarray, shape (M, D)

The points we want to predict the values at.

y_starndarray, shape (M, 1), optional

True values at the points.

s2_starndarray, shape (M, 1), optional

Noise variance at the points.

add_noisebool, defaults to False

Whether to add the observation noise, which enters on the diagonal, to the returned covariance.

Returns:

mundarray, shape (M, sample_N)

Posterior mean at the requested points for each hyperparameter sample.

covndarray, shape (M, M, sample_N)

Covariance matrix for each hyperparameter sample. Its diagonal is not clamped at zero, unlike the variances predict() returns, so on a nearly singular posterior an entry can come out slightly negative.

Raises:

ValueError

Raised by the check of the shapes of the data given: when x_star is neither a 1-D nor a 2-D array or does not have D columns, when y_star does not hold one value per row of x_star, or when an array s2_star does not hold one variance per row of x_star.

TypeError

Raised when s2_star is neither an array, a number nor None, such as a list.

LinAlgError

Raised when a posterior holds the negative inverse of the training covariance without its Cholesky factor, as one pickled by gpyreg 1.3.1 or earlier does where the smallest noise variance at the training inputs is below 1e-6, and the factorization that computes the factor again fails even after its noise is multiplied tenfold, up to ten times, or fails at all where the GP raises on a failed factorization (raise_on_cholesky_failure).

quad(mu, sigma, compute_var: bool = False, separate_samples: bool = False)

Bayesian quadrature for a Gaussian Process.

Compute the integral of a function represented by a Gaussian Process with respect to a given Gaussian measure.

Parameters:

muarray_like

Either a array of shape (N, D) with each row containing the mean of a single Gaussian measure, or a single floating point number which is interpreted as an array of shape (1, D). A one-dimensional array of length D is one measure of D dimensions, as gplite_quad.m’s size(mu, 1) reads a row vector.

sigmaarray_like

Either a array of shape (N, D) with each row containing the standard deviation of a single Gaussian measure, or a single floating point number which is interpreted as an array of shape (1, D), or an array of shape (N, 1) with one standard deviation per measure, the same in every dimension. A one-dimensional array of length D is one measure, as for mu.

compute_varbool, defaults to False

Whether to compute variance for each integral.

separate_samplesbool, defaults to False

Whether to return the results separately for each hyperparameter sample, or averaged.

Returns:

Fndarray

The computed integrals in an array with shape (N, 1) if samples are averaged and shape (N, hyp_samples) if requested separately.

F_varndarray, optional

The computed variances of the integrals in an array with shape (N, 1) if samples are averaged and shape (N, hyp_samples) if requested separately.

Raises:

ValueError

Raised when the method is called and the covariance of the GP is not squared exponential or has periods, or the mean function is none of the zero, constant and negative quadratic means.

ValueError

Raised when the GP has no training data or no posterior factors, or when mu does not have one column per input dimension, or sigma neither one nor one per dimension.

LinAlgError

Raised when compute_var is True and a posterior holds the negative inverse of the training covariance without its Cholesky factor, as one pickled by gpyreg 1.3.1 or earlier does where the smallest noise variance at the training inputs is below 1e-6, and the factorization that computes the factor again fails even after its noise is multiplied tenfold, up to ten times, or fails at all where the GP raises on a failed factorization (raise_on_cholesky_failure).

random_function(X_star: ndarray, add_noise: bool = False, rng=None)

Draw a random function from the Gaussian Process.

Parameters:

X_starndarray, shape (M, D)

The points at which to evaluate the drawn function.

add_noisebool, defaults to False

Whether to add noise to the values of the drawn function.

rngNone, numpy.random.Generator or seed, optional

Where the draws come from (the hyperparameter sample and the function values). None (default) keeps NumPy’s global legacy stream, as before generators were supported. See gpyreg.rng.resolve_rng().

Returns:

f_starndarray, shape (M, 1)

The values of the drawn function at the requested points.

Raises:

LinAlgError

Raised when the covariance of the draw, as computed, has a negative eigenvalue beyond the rounding of the prior variance at X_star.

LinAlgError

Raised when the posterior drawn from holds the negative inverse of the training covariance without its Cholesky factor, as one pickled by gpyreg 1.3.1 or earlier does where the smallest noise variance at the training inputs is below 1e-6, and the factorization that computes the factor again fails even after its noise is multiplied tenfold, up to ten times, or fails at all where the GP raises on a failed factorization (raise_on_cholesky_failure).

set_bounds(bounds: dict = None)

Set the hyperparameter lower and upper bounds.

Parameters:

boundsdict, optional

A dictionary of GP hyperparameter names and tuples of their lower and upper bounds. All hyperparameters need to appear in the dictionary. Use the value None to set the bounds of a specific hyperparameter to the default recommended values. If bounds=None, all hyperparameter bounds will be set to their default recommended values.

Raises:

ValueError

Raised when bounds is missing the entry of an expected hyperparameter.

ValueError

Raised when bounds is given, but a specified hyperparameter is unknown.

ValueError

Raised when a lower bound is above the upper bound of the same hyperparameter. Equal bounds fix a hyperparameter.

set_hyperparameters(hyp_new: object, compute_posterior: bool = True)

Set new hyperparameters for the Gaussian Process.

The hyperparameters go to update(), and a call that raises leaves the GP as it was before the call, with the hyperparameters and posteriors that it held.

Parameters:

hyp_newobject

The new hyperparameters. This can be an array of shape (hyp_samples, hyp_N) where hyp_N is the number of hyperparameters, and hyp_samples is the amount of hyperparameter samples, a single dictionary with hyperparameter names and values, or a list of dictionaries. Passing a single dictionary or a list with one dictionary is equivalent.

compute_posteriorbool, defaults to True

Whether to compute the posterior for the new hyperparameters.

Raises:

ValueError

Raised when hyp_new is an array of the wrong shape.

ValueError

Raised by update(), before it changes anything, when compute_posterior is True, the GP holds training inputs and targets, and a hyperparameter is NaN (not set).

LinAlgError

Raised by update() when the Cholesky decomposition of the training covariance fails even after its noise is multiplied tenfold, up to ten times, or fails at all where the GP raises on a failed factorization (raise_on_cholesky_failure); the GP is left as it was.

set_priors(priors: dict = None)

Set the hyperparameter priors.

Parameters:

priorsdict, optional

A dictionary of GP hyperparameter names and tuples of their priors. All hyperparameters need to appear in the dictionary. Use the value None to set no priors for a hyperparameter. If priors=None, all hyperparameter priors are removed. Within a block of several hyperparameters, a coordinate whose location (mu, or a and b for the smooth-box families) and sigma are both NaN has no prior; every other coordinate needs a finite location, with a <= b, and a finite, positive sigma. A block with no prior in any coordinate holds none, and a GP whose blocks all hold none has no priors. A smooth box with a == b is the Gaussian (or Student’s t) centred at a. Degrees of freedom df that are zero, infinite or NaN make a "student_t" prior "gaussian", as gplite_hypprior.m reads them, and a "smoothbox_student_t" prior "smoothbox", gpyreg’s own reading (gplite has no smooth-box priors). For the duration of fit(), a NaN df takes the value of its option df_base instead.

Raises:

ValueError

Raised when priors is given, but missing the entry of an expected hyperparameter.

ValueError

Raised when priors is given, but a specified hyperparameter is unknown.

ValueError

Raised when a coordinate that has a prior is given a sigma that is not finite and positive, or a location that is not finite, or a smooth box whose a is above its b.

update(X_new: ndarray = None, y_new: ndarray = None, s2_new: ndarray = None, hyp: ndarray = None, compute_posterior: bool = True)

Add new data to the Gaussian Process.

An update that raises leaves the GP as it was before the call, with the training data and posteriors that it held.

Parameters:

X_newndarray, shape (N, D), optional

New training inputs that will be added to the old training inputs.

y_newnarray, shape (N, 1) optional

New training targets that will be added to the old training targets.

s2_newndarray, shape (N, 1), optional

New input-dependent noise that will be added to the old training inputs. When only some training points come with a user-provided variance, the others are assigned zero.

hypndarray, shape (hyp_samples, hyp_N), optional

New hyperparameters that will replace the old ones, one row per hyperparameter sample.

compute_posteriorbool, defaults to True

Whether to compute the new posterior or not.

Raises:

ValueError

Raised, before the update changes anything, when the update would make the number of targets, or of noise variances, that the GP holds differ from its number of inputs; a GP may hold inputs without targets, and data without noise variances. So y_new given without X_new, one value per input held, is taken where the GP holds no targets, as s2_new is where it holds no variances, and either is refused where the GP holds what it gives or holds no inputs; X_new without y_new is refused where the GP holds targets, and X_new with y_new where it holds inputs without targets.

ValueError

Raised by the check of the shapes of the training data given, before the update changes anything: when X_new is neither a 1-D nor a 2-D array or does not have D columns, when y_new does not hold one value per row of X_new, or when an array s2_new does not hold one variance per row of X_new; without X_new, the rows are the inputs that the GP holds.

TypeError

Raised when s2_new is neither an array, a number nor None, such as a list, before the update changes anything.

LinAlgError

Raised when the Cholesky decomposition of the training covariance fails even after its noise is multiplied tenfold, up to ten times, or fails at all where the GP raises on a failed factorization (raise_on_cholesky_failure), in a posterior computed in full or in the full recomputation to which a single-point update falls back; the GP is left as it was.

ValueError

Raised when hyp is not a 2D array with one column per hyperparameter of the GP.

ValueError

Raised, before the update changes anything, when compute_posterior is True, the GP would hold training inputs and targets, every posterior is computed in full, and a hyperparameter is NaN (not set), as it is on a GP whose hyperparameters were never given.

Posterior

class gpyreg.gaussian_process.Posterior(hyp, alpha, sW, L, sn2_mult, Lchol, sl=None, L_factor=None)

GP posterior for one hyperparameter vector.

Stores the coefficients and matrix factor that prediction, sampling and quadrature need for a fixed set of hyperparameters, so that repeated calls reuse them instead of factorizing the training covariance again. GP.posteriors holds one instance per hyperparameter sample. Below, K is the training covariance, sn2 the vector of training noise variances, m the mean function evaluated at the training inputs and sl = min(sn2) * sn2_mult at the time the posterior was computed; sl scales the Cholesky factor when L_chol is True. Every attribute except hyp is None when the posterior has not been computed, after GP.clean() or GP.update() with compute_posterior=False.

Attributes:

hypndarray, shape (hyp_N,)

The hyperparameters this posterior was computed for.

alphandarray, shape (N, 1)

inv(K + sn2_mult * diag(sn2)) @ (y - m), the weights of the training points in the posterior mean.

sWndarray, shape (N, 1)

1 / sqrt(sl) in every entry; the square root of the noise precision used to scale the factorization.

Lndarray, shape (N, N)

If L_chol is True, the upper triangular Cholesky factor of (K + sn2_mult * diag(sn2)) / sl. Otherwise -inv(K + sn2_mult * diag(sn2)), the low-noise representation, as gplite_core.m has it where the noise is small.

sn2_multint

Multiplier applied to the noise variances, increased in powers of ten until the Cholesky decomposition succeeds.

L_cholbool

Whether L is a Cholesky factor (min(sn2) >= 1e-6) or a negative inverse.

slfloat

min(sn2) * sn2_mult when the posterior was computed; the scale of L when L_chol is True. Rank-one updates extend the factorization with this scale, so it can differ from the minimum of the current training noise after observations have been appended.

L_factorndarray, shape (N, N) or None

If L_chol is False, the upper triangular Cholesky factor of K + sn2_mult * diag(sn2), the matrix whose negative inverse L is. GP.predict(), GP.predict_full(), GP.quad() and GP.random_function() form their covariances from it, and a single-point GP.update() takes from it the predictive variance of the new point and extends it: formed from L, these would carry the rounding of the inverse, which grows as the noise shrinks. None if L_chol is True, where L is that factor, scaled. A posterior pickled by gpyreg 1.3.1 or earlier has no such attribute: the predictions, the quadrature and the draws compute the factor again at each call, and a single-point update recomputes the posterior in full.