CGP (version 2.1-1)

CGP: Fit 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

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 (\(\lambda\) 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 (\(\theta\) in the paper)

alpha_l

Optional. The lower bound for all correlation parameters in the local GP (\(\alpha\) in the paper). It is also the upper bound for all correlation parameters in the global GP (the \(\theta\)). 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 \(\kappa\), where we define \(\alpha_j=\theta_j+\kappa\) for \(j=1,\ldots,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)))

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 \((\lambda)\)

theta

Vector of estimated correlation parameters \((\theta)\) in the global GP

alpha

Vector of estimated correlation parameters \((\alpha)\) 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 \((\mu)\) for global trend

tau2

Estimated variance \((\tau^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}+\lambda\mathbf{\Sigma}^{1/2}\mathbf{L}\mathbf{\Sigma}^{1/2})^{-1}\)

temp_matrix

Matrix of \((\mathbf{G}+\lambda\mathbf{\Sigma}^{1/2}\mathbf{L}\mathbf{\Sigma}^{1/2})^{-1} (\mathbf{y}- \hat{\mu}\mathbf{1})\)

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 \(\alpha_j=\theta_j+\kappa\) for \(j=1,\ldots,p\), fitting the CGP model only requires estimating \(p+3\) unknown parameters, which is comparable to fitting a stationary GP model (\(p\) unknown parameters).

References

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

See Also

predict.CGP, print.CGP, summary.CGP

Examples

# NOT RUN {
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
# }