# CGP: Fit composite Gaussian process models In CGP: Composite Gaussian Process Models

## Description

Estimate parameters in the composite Gaussian process (CGP) model using maximum likelihood methods. Calculate the root mean squared (leave-one-out) cross validation error for diagnosis, and export intermediate values to facilitate predict.CGP function.

## Usage

 1 2 CGP(X, yobs, nugget_l = 0.001, num_starts = 5, theta_l = NULL, alpha_l = NULL, kappa_u = NULL) 

## Arguments

 X The design matrix yobs The vector of response values, corresponding to the rows of X nugget_l Optional, default is “0.001”. The lower bound for the nugget value (λ in the paper) num_starts Optional, default is “5”. Number of random starts for optimizing the likelihood function theta_l Optional, default is “0.0001”. The lower bound for all correlation parameters in the global GP (θ in the paper) alpha_l Optional. The lower bound for all correlation parameters in the local GP (α in the paper). It is also the upper bound for all correlation parameters in the global GP (the θ). Default is log(100)*mean(1/dist(Stand_X)^2), where Stand_X<-apply(X,2,function(x) (x-min(x))/max(x-min(x))). Please refer to Ba and Joseph (2012) for details kappa_u Optional. The upper bound for κ, where we define α_j=θ_j+κ for j=1,…,p. Default value is log(10^6)*mean(1/dist(Stand_X)^2), \ where Stand_X<-apply(X,2,function(x) (x-min(x))/max(x-min(x)))

## Details

This function fits a composite Gaussian process (CGP) model based on the given design matrix X and the observed responses yobs. The fitted model consists of a smooth GP to caputre the global trend and a local GP which is augmented with a flexible variance model to capture the change of local volatilities. For p input variables, such two GPs involve 2p+2 unknown parameters in total. As demonstrated in Ba and Joseph (2012), by assuming α_j=θ_j+κ for j=1,…,p, fitting the CGP model only requires estimating p+3 unknown parameters, which is comparable to fitting a stationary GP model (p unknown parameters).

## Value

This function fits the CGP model and returns an object of class “CGP”. Function predict.CGP can be further used for making new predictions and function summary.CGP can be used to print a summary of the “CGP” object.

An object of class “CGP” is a list containing at least the following components:

 lambda Estimated nugget value (λ) theta Vector of estimated correlation parameters (θ) in the global GP alpha Vector of estimated correlation parameters (α) in the local GP bandwidth Estimated bandwidth parameter (b) in the variance model rmscv Root mean squared (leave-one-out) cross validation error Yp_jackknife Vector of Jackknife (leave-one-out) predicted values mu Estimated mean value (μ) for global trend tau2 Estimated variance (τ^2) for global trend beststart Best starting value found for optimizing the log-likelihood objval Optimal objective value for the negative log-likelihood (up to a constant) var_names Vector of input variable names Sig_matrix Diagonal matrix containing standardized local variances at each of the design points sf Standardization constant for rescaling the local variance model res2 Vector of squared residuals from the estimated global trend invQ Matrix of (\mathbf{G}+λ\mathbf{Σ}^{1/2}\mathbf{L}\mathbf{Σ}^{1/2})^{-1} temp_matrix Matrix of (\mathbf{G}+λ\mathbf{Σ}^{1/2}\mathbf{L}\mathbf{Σ}^{1/2})^{-1} (\mathbf{y}- \hat{μ}\mathbf{1})

## Author(s)

Shan Ba <[email protected]> and V. Roshan Joseph <[email protected]>

## References

Ba, S. and V. Roshan Joseph (2012) “Composite Gaussian Process Models for Emulating Expensive Functions”. Annals of Applied Statistics, 6, 1838-1860.

predict.CGP, print.CGP, summary.CGP
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 x1<-c(0,.02,.075,.08,.14,.15,.155,.156,.18,.22,.29,.32,.36, .37,.42,.5,.57,.63,.72,.785,.8,.84,.925,1) x2<-c(.29,.02,.12,.58,.38,.87,.01,.12,.22,.08,.34,.185,.64, .02,.93,.15,.42,.71,1,0,.21,.5,.785,.21) X<-cbind(x1,x2) yobs<-sin(1/((x1*0.7+0.3)*(x2*0.7+0.3))) ## Not run: #Fit the CGP model #Increase the lower bound for nugget to 0.01 (Optional) mod<-CGP(X,yobs,nugget_l=0.01) summary(mod) mod$objval #-27.4537 mod$lambda #0.6210284 mod$theta #6.065497 8.093402 mod$alpha #143.1770 145.2049 mod$bandwidth #1 mod$rmscv #0.5714969 ## End(Not run)