homGP

Methods

`copy`()	Make a copy of the model, which is useful in tandem with the update function
`dlogLikHom`(X0, Z0, Z, mult, theta, g[, ...])	Log likelihood gradient under homoskedastic noise
`get`(key)	General get item (retrives key from self.__dict__)
`logLikHom`(X0, Z0, Z, mult, theta, g[, ...])	Log Likelihood under homoskedastic noise
`mleHomGP`(X, Z[, lower, upper, known, ...])	Gaussian process modeling with homoskedastic noise.
`plot`([type])	Plot model behavior.
`predict`(x[, xprime, interval, ...])	Prediction under homoskedastic noise
`rebuild`([robust])	Rebuilds inverse covariance matrix in homGP object (usually after saving).
`strip`()	Removes larger model objects
`summary`()	Print a summary of the model object and the MLE routine
`update`(Xnew, Znew[, ginit, lower, upper, ...])	Update model object with new observations

copy()[source]

Make a copy of the model, which is useful in tandem with the update function

Parameters:: None
Returns:: newmodel
Return type:: (deep) copy of model

dlogLikHom(X0, Z0, Z, mult, theta, g, beta0=None, covtype='Gaussian', eps=np.float64(1.4901161193847656e-08), components=('theta', 'g'))[source]

Log likelihood gradient under homoskedastic noise

Parameters:

X0 (ndarray_like) – unique design matrix (of size nxd)
Z0 (ndarray_like) – averaged responses (size nx1)
Z (ndarray_like) – observation (output) vector (size N)
mult (ndarray_like) – number of replicates as each design location. Note that mult.sum()==N
theta (ndarray_like) – lengthscale
g (float) – noise variance
beta0 (float) – trend
covtype (str) – covariance kernel to use
eps (float) – amount of jitter on diagonal of covariance matrix
components (list) – directions for partial derivatives, one or both of ‘theta’ (lengthscales) or ‘g’ (for noise)

Returns:

out – gradient with respect to hyperparameters

Return type:

ndarray_like

get(key)[source]: General get item (retrives key from self.__dict__)

logLikHom(X0, Z0, Z, mult, theta, g, beta0=None, covtype='Gaussian', eps=np.float64(1.4901161193847656e-08))[source]

Log Likelihood under homoskedastic noise

Parameters:

X0 (ndarray_like) – unique design matrix (of size nxd)
Z0 (ndarray_like) – averaged responses (size nx1)
Z (ndarray_like) – observation (output) vector (size N)
mult (ndarray_like) – number of replicates as each design location. Note that mult.sum()==N
theta (ndarray_like) – lengthscale
g (float) – noise variance
beta0 (float) – trend
covtype (str) – covariance kernel to use
eps (float) – amount of jitter on diagonal of covariance matrix

Returns:

loglik – log-likelihood

Return type:

float

mleHomGP(X, Z, lower=None, upper=None, known={}, noiseControl={'g_bounds': (np.float64(1.4901161193847656e-08), 100.0)}, init={}, covtype='Gaussian', maxit=100, eps=np.float64(1.4901161193847656e-08), settings={'factr': 10000000.0, 'returnKi': True})[source]

Gaussian process modeling with homoskedastic noise.

You may also call this function as model.mle

Gaussian process regression under homoskedastic noise based on maximum likelihood estimation of the hyperparameters. This function is enhanced to deal with replicated observations.

Parameters:

X (ndarray_like) – matrix of all designs, one per row, or list with elements: - X0 matrix of unique design locations, one point per row - Z0 vector of averaged observations, of length len(X0) - mult number of replicates at designs in X0, of length len(X0)
Z (ndarray_like) – Z vector of all observations. If using a list with X, Z has to be ordered with respect to X0, and of length sum(mult)
lower (ndarray_like) – optional bounds for the theta parameter (see :func: covariance_functions.cov_gen for the exact parameterization). In the multivariate case, it is possible to give vectors for bounds (resp. scalars) for anisotropy (resp. isotropy)
upper (ndarray_like) – optional bounds for the theta parameter (see :func: covariance_functions.cov_gen for the exact parameterization). In the multivariate case, it is possible to give vectors for bounds (resp. scalars) for anisotropy (resp. isotropy)
noiseControl (dict) –
dict with element:
- g_bounds vector providing minimal and maximal noise to signal ratio (default to (sqrt(MACHINE_DOUBLE_EPS), 100)).
settings (dict) –
dict for options about the general modeling procedure, with elements:
- return_Ki boolean to include the inverse covariance matrix in the object for further use (e.g., prediction).
- factr (default to 1e7) and pgtol are available to be passed to options for L-BFGS-B in :func: scipy.optimize.minimize.
eps (float) – jitter used in the inversion of the covariance matrix for numerical stability
known (dict) – optional dict of known parameters (e.g. beta0, theta, g)
init (dict) –
optional lists of starting values for mle optimization:
- theta_init initial value of the theta parameters to be optimized over (default to 10% of the range determined with lower and upper)
- g_init vector of nugget parameter to be optimized over
covtype (str) – covariance kernel type, either 'Gaussian', 'Matern5_2' or 'Matern3_2', see :func: ~covariance_functions.cov_gen
maxit (int) – maximum number of iterations for L-BFGS-B of :func: scipy.optimize.minimize dedicated to maximum likelihood optimization

Return type:

None

Notes

The global covariance matrix of the model is parameterized as nu_hat * (C + g * diag(1/mult)) = nu_hat * K, with C the correlation matrix between unique designs, depending on the family of kernel used (see :func: ~hetgpy.covariance_functions.cov_gen for available choices) and values of lengthscale parameters. nu_hat is the plugin estimator of the variance of the process.

It is generally recommended to use :func: ~find_reps.find_reps to pre-process the data, to rescale the inputs to the unit cube and to normalize the outputs.

Returns:

theta: unless given, maximum likelihood estimate (mle) of the lengthscale parameter(s),
nu_hat: plugin estimator of the variance,
g: unless given, mle of the nugget of the noise/log-noise process,
trendtype: either "SK" if beta0 is provided, else "OK",
beta0: constant trend of the mean process, plugin-estimator unless given,
ll: log-likelihood value, (ll_non_pen) is the value without the penalty,
nit_opt, msg: counts and message returned by :func:scipy.optimize.minimize
used_args: list with arguments provided in the call to the function,
- Ki, Kgi: inverse of the covariance matrices of the mean and noise processes (not scaled by nu_hat and nu_hat_var),
- X0, Z0, Z, eps, logN, covtype: values given in input,
time: time to train the model, in seconds.

See also ~hetgpy.homGP.homGP.predict for predictions, ~hetgpy.homGP.update for updating an existing model. summary and plot functions are available as well. ~homTP.mleHomTP provide a Student-t equivalent.

Return type:

self, with the following attributes

Parameters:

X (ArrayLike)
Z (ArrayLike)
lower (TypeAliasForwardRef('ArrayLike') | None)
upper (TypeAliasForwardRef('ArrayLike') | None)
known (dict)
noiseControl (dict)
init (dict)
covtype (str)
maxit (int)
eps (float)
settings (dict)

Examples

>>> from hetgpy import homGP
>>> from hetgpy.example_data import mcycle
>>> m = mcycle()
>>> model = homGP()
>>> model.mle(m['times'],m['accel'],lower=[1.0],upper=[10.0],covtype="Matern5_2")

References

M. Binois, Robert B. Gramacy, M. Ludkovski (2018), Practical heteroskedastic Gaussian process modeling for large simulation experiments, Journal of Computational and Graphical Statistics, 27(4), 808–821. Preprint available on arXiv:1611.05902.

plot(type='diagnostics')[source]

Plot model behavior. See hetgpy.plot for more details on types of plots

Parameters:: type (str) – type of plot: currently supported are “diagnostics” and “iterates”
Returns:: fig, ax
Return type:: figure and axis object

predict(x, xprime=None, interval=None, interval_lower=None, interval_upper=None, **kw)[source]

Prediction under homoskedastic noise

Parameters:

x (ndarray_like) – matrix of designs locations to predict at (one point per row)
xprime (ndarray_like) – optional second matrix of predictive locations to obtain the predictive covariance matrix between x and xprime
interval (str) – one of ‘confidence’ or ‘predictive’ which is a convenience method to return confidence/predictive intervals corresponding to interval_lower and interval_upper
interval_lower (float) – lower of confidence/predictive interval
interval_upper (float) – upper of confidence/predictive interval
nugs_only (bool (default False)) – if True, only return noise variance prediction
kwargs (dict) – optional additional elements (only used for nugs_only)

Returns:

mean: kriging mean;
sd2: kriging variance (filtered, e.g. without the nugget values)
nugs: noise variance prediction
cov: (returned if xprime is given) predictive covariance matrix between x and xprime
confidence_interval: prediction with kriging variance only
predictive_interval: prediction with kriging and noise variance

Return type:

dict with elements

rebuild(robust=False)[source]

Rebuilds inverse covariance matrix in homGP object (usually after saving). Works in tandem with strip

Parameters:: robust (bool) – use np.linalg.pinv for covariance matrix inversion. Otherwise use cholesky
Return type:: self with rebuilt inverse covariance matrix

strip()[source]

Removes larger model objects

Can be rebuilt with rebuild

summary()[source]: Print a summary of the model object and the MLE routine

update(Xnew, Znew, ginit=0.01, lower=None, upper=None, noiseControl=None, settings=None, known={}, maxit=100)[source]

Update model object with new observations

Parameters:

Xnew (ndarray_like) – new inputs (one observation per row)
Znew (ndarray_like) – new observations
lower (ndarray_like) – lower bound for lengthscales. If not provided extracted from self
upper (ndarray_like) – upper bound for lengthscales. If not provided extracted from self
noiseControl (dict) – noise bounds. If not provided extracted from self
settings (dict) – options for optimization. If not provided extracted from self
known (dict) – known hyperparameters to fix at values
maxit (int) – maximum number of iterations for hyperparameter optimization
ginit (float)

Return type:

self, potentially with updated hyperparameters

Notes

If hyperparameters do not need to be updated, maxit can be set to 0. In this case it is possible to pass NAs in Znew, then the model can still be used to provide updated variance predictions

Example

>>> import numpy as np, from hetgpy import homGP
>>> X = np.array([1, 2, 3, 4, 12]).reshape(-1,1)
>>> Y = np.array([0, -1.75, -2, -0.5, 5])
>>> model = homGP().mle(X,Y,known={'theta':np.array([10.0]),'g':1e-8})
>>> Xnew  = np.array([2.5]).reshape(-1,1)
>>> Znew  = model.predict(Xnew)['mean']
>>> model.update(Xnew,Znew)