Movatterモバイル変換

Type:	Package
Title:	Gaussian Process Laboratory
Version:	0.5.8
Date:	2024-11-19
Description:	Gaussian process regression with an emphasis on kernels. Quantitative and qualitative inputs are accepted. Some pre-defined kernels are available, such as radial or tensor-sum for quantitative inputs, and compound symmetry, low rank, group kernel for qualitative inputs. The user can define new kernels and composite kernels through a formula mechanism. Useful methods include parameter estimation by maximum likelihood, simulation, prediction and leave-one-out validation.
License:	GPL-3
Depends:	Rcpp (≥ 0.10.5), methods, testthat, nloptr, lattice
Suggests:	DiceKriging, DiceDesign, inline, foreach, knitr, ggplot2,reshape2, corrplot
Imports:	MASS, numDeriv, stats4, doParallel, doFuture, utils
LinkingTo:	Rcpp
RoxygenNote:	6.0.1
Collate:	'CovFormulas.R' 'allGenerics.R' 'checkGrad.R' 'covComp.R''covMan.R' 'covQual.R' 'q1CompSymm.R' 'q1Symm.R' 'q1LowRank.R''covQualNested.R' 'covQualOrd.R' 'covRadial.R' 'covTS.R''covTP.R' 'covANOVA.R' 'covZZAll.R' 'gp.R' 'kFuns.R''kernelNorm.R' 'kernels1d_Call.R' 'logLikFuns.R' 'methodGLS.R''methodMLE.R' 'miscUtils.R' 'prinKrige.R' 'q1Diag.R''simulate_gp.R' 'warpFuns.R'
NeedsCompilation:	yes
Packaged:	2024-11-19 13:49:14 UTC; yves
Author:	Yves Deville [aut], David Ginsbourger [aut], Olivier Roustant [aut, cre], Nicolas Durrande [ctb]
Maintainer:	Olivier Roustant <roustant@insa-toulouse.fr>
Repository:	CRAN
Date/Publication:	2024-11-19 14:40:03 UTC

Gaussian Process Laboratory

Description

Laboratory Package for Gaussian Process interpolation, regression andsimulation, with an emphasis on user-defined covariance kernels.

Details

Warning

As a lab,kergp may strongly evolve in its future life. Usersinterested in stable software for the Analysis of Computer Experimentsare encouraged to use other packages such asDiceKriginginstead.

Note

This package was developed within the frame of the ReDice Consortium,gathering industrial partners (CEA, EDF, IFPEN, IRSN, Renault) andacademic partners (Mines Saint-Étienne, INRIA, and the University ofBern) around advanced methods for Computer Experiments.

Author(s)

Yves Deville (Alpestat), David Ginsbourger (University of Bern),Olivier Roustant (INSA Toulouse), with contributions fromNicolas Durrande (Mines Saint-Étienne).

Maintainer: Olivier Roustant, <roustant@insa-toulouse.fr>

References

Nicolas Durrande, David Ginsbourger, Olivier Roustant (2012)."Additive covariance kernels for high-dimensional gaussian process modeling".Annales de la Faculté des Sciences de Toulouse, 21 (3):481-499,doi:10.5802/afst.1342.

Nicolas Durrande, David Ginsbourger, Olivier Roustant, Laurent Carraro(2013)."ANOVA kernels and RKHS of zero mean functions for model-based sensitivity analysis".Journal of Multivariate Analysis, 115, 57-67,doi:10.1016/j.jmva.2012.08.016.

David Ginsbourger, Xavier Bay, Olivier Roustant, Laurent Carraro(2012)."Argumentwise invariant kernels for the approximation of invariant functions".Annales de la Faculté des Sciences de Toulouse, 21 (3):501-527,doi:10.5802/afst.1343.

David Ginsbourger, Nicolas Durrande, Olivier Roustant (2013)."Kernels and designs for modelling invariant functions: From group invariance to additivity"mODa 10 - Advances in Model-Oriented Design andAnalysis. Contributions to Statistics, 107-115,doi:10.1007/978-3-319-00218-7_13.

Olivier Roustant, David Ginsbourger, Yves Deville (2012)."DiceKriging, DiceOptim: Two R Packages for the Analysis ofComputer Experiments by Kriging-Based Metamodeling and Optimization".Journal of Statistical Software, 51(1), 1-55,doi:10.18637/jss.v051.i01.

Examples

## ------------------------------------------------------------------## Gaussian process modelling of function with invariance properties, ## by using an argumentwise invariant kernel## ------------------------------------------------------------------## -- define manually an argumentwise invariant kernel --kernFun <- function(x1, x2, par) {  h <- (abs(x1) - abs(x2)) / par[1]  S <- sum(h^2)  d2 <- exp(-S)  K <- par[2] * d2  d1 <- 2 * K * S / par[1]     attr(K, "gradient") <- c(theta = d1,  sigma2 = d2)  return(K)}## ---------------------------------------------------------------## quicker: with Rcpp; see also an example  with package inline## in "gp" doc. file. Note that the Rcpp "sugar" fucntions are## vectorized, so no for loops is required.## ---------------------------------------------------------------## Not run:     cppFunction('        NumericVector cppKernFun(NumericVector x1, NumericVector x2,                                  NumericVector par){        int n1 = x1.size();        double S, d1, d2;         NumericVector K(1), h(n1);        h = (abs(x1) - abs(x2)) / par[0];  // sugar function "abs"        S = sum(h * h);                    // sugar "*" and "sum"         d2 = exp(-S);        K[0] = par[1] * d2;        d1 = 2 * K[0] * S / par[0];           K.attr("gradient") = NumericVector::create(Named("theta", d1),                                                   Named("sigma2", d2));        return K;     }')## End(Not run)## ---------------------------------------------------------------## Below: with the R-based code for the kernel namely 'kernFun'.## You can also replace 'kernFun' by 'cppKernFun' for speed.## ---------------------------------------------------------------covSymGauss <- covMan(kernel = kernFun,                      hasGrad = TRUE,                      label = "argumentwise invariant",                      d = 2,                      parLower = c(theta = 0.0, sigma2 = 0.0),                      parUpper = c(theta = Inf, sigma2 = Inf),                      parNames = c("theta", "sigma2"),                      par = c(theta = 0.5, sigma2 = 2))covSymGauss## -- simulate a path from the corresponding GP --nGrid <- 24; n <- nGrid^2; d <- 2xGrid <- seq(from = -1, to = 1, length.out = nGrid)Xgrid <- expand.grid(x1 = xGrid, x2 = xGrid)Kmat <- covMat(object = covSymGauss, X = Xgrid,               compGrad = FALSE, index = 1L)library(MASS)set.seed(1)ygrid <- mvrnorm(mu = rep(0, n), Sigma = Kmat)## -- extract a design and the corr. response from the grid --nDesign <- 25tab <- subset(cbind(Xgrid, ygrid), x1 > 0 & x2 > 0)rowIndex <- seq(1, nrow(tab), length = nDesign)X <- tab[rowIndex, 1:2]y <- tab[rowIndex, 3]opar <- par(mfrow = c(1, 3))contour(x = xGrid, y = xGrid,        z = matrix(ygrid, nrow = nGrid, ncol = nGrid),         nlevels = 15)abline(h = 0, v = 0, col = "SpringGreen3")points(x2 ~ x1, data = X, type = "p", pch = 21,       col = "orangered", bg = "yellow", cex = 0.8)title("GRF Simulation")## -- Fit the Gaussian process model (trend + covariance parameters) -- covSymGausssymgp <- gp(formula = y ~ 1, data = data.frame(y, X),            inputs = names(X),            cov = covSymGauss,            parCovIni = c(0.1, 2),            varNoiseIni = 1.0e-8,            varNoiseLower = 0.9e-8, varNoiseUpper = 1.1e-8)# mind that the noise is not a symmetric kernel# so varNoiseUpper should be chosen as small as possible.summary(symgp)## -- predict and compare --predSymgp <- predict(object = symgp, newdata = Xgrid, type = "UK")contour(x = xGrid, y = xGrid,        z = matrix(predSymgp$mean, nrow = nGrid, ncol = nGrid),        nlevels = 15)abline(h = 0, v = 0, col = "SpringGreen3")points(x2 ~ x1, data = X, type = "p", pch = 21,       col = "orangered", bg = "yellow", cex = 0.8)title("Kriging mean")contour(x = xGrid, y = xGrid,        z = matrix(predSymgp$sd, nrow = nGrid, ncol = nGrid),        nlevels = 15)abline(h = 0, v = 0, col = "SpringGreen3")points(x2 ~ x1, data = X, type = "p", pch = 21,       col = "orangered", bg = "yellow", cex = 0.8)title("Kriging s.d.")par(opar)

Coerce acovTP Object into a List

Description

Coerce acovTP object representing a Tensor-Productcovariance kernel on thed-dimensional Euclidean spaceinto a list containingd one-dimensional kernels.

Usage

## S4 method for signature 'covTP'as.list(x)

Arguments

Value

A list with lengthd ord + 1 whered is the"dimension" slotx@d of the objectx. The firstdelements of the list are one-dimensionalcorrelation kernelobjects with class"covTP". Whenx is acovariance kernel (as opposed to acorrelation kernel),the list contains one more element which gives the variance.

Caution

Whenx is not a correlation kernel the(d + 1)-th element of the returned list may be different infuture versions: it may be a constant covariance kernel.

See Also

covTP andcovTP-class.

Examples

set.seed(123)d <- 6myCov1 <- covTP(d = d, cov = "corr")coef(myCov1) <- as.vector(simulPar(myCov1, nsim = 1))as.list(myCov1)## more examples and check the value of a 'covMat'L <- list()myCov <- list()myCov[[1]] <- covTP(d = d, cov = "corr")coef(myCov[[1]]) <- as.vector(simulPar(myCov[[1]], nsim = 1))L[[1]] <- as.list(myCov[[1]])myCov[[2]] <- covTP(k1Fun1 = k1Fun1PowExp, d = d, cov = "corr")coef(myCov[[2]]) <- as.vector(simulPar(myCov[[2]], nsim = 1))L[[2]] <- as.list(myCov[[2]])myCov[[3]] <- covTP(k1Fun1 = k1Fun1PowExp, d = d, iso1 = 0L, cov = "corr")coef(myCov[[3]]) <- as.vector(simulPar(myCov[[3]], nsim = 1))L[[3]] <- as.list(myCov[[3]])n <- 10X <- matrix(runif(n * d), nrow = n,            dimnames = list(NULL, paste("x", 1:d, sep = "")))for (iTest in 1:3) {   C <- covMat(L[[iTest]][[1]], X[ , 1, drop = FALSE])   for (j in 2:d) {      C <- C * covMat(L[[iTest]][[j]], X[ , j, drop = FALSE])   }   CTest <- covMat(myCov[[iTest]], X)   print(max(abs(abs(C - CTest))))}

Check the Gradient Provided in acovMan Object

Description

Check the gradient provided in acovMan object.

Usage

checkGrad(object, sym = TRUE,          x1 = NULL, n1 = 10,          x2 = NULL, n2 = NULL,          XLower = NULL, XUpper = NULL,          plot = TRUE)

Arguments

Details

Each of the two matricesx1 andx2 withn1 andn2 rows can be given or instead be drawn at random. The matrixof kernel values with dimensionc(n1, n2) is computed, togetherwith its gradient with dimensionc(n1, n2, npar) wherenpar is the number of parameters of the kernel. A numericaldifferentiation w.r.t. the kernel parameters is performed for thekernel value atx1 andx2, and the result is compared tothat provided by the kernel function (the function described in theslot named"kernel" ofobject). Note that the value ofthe parameter vector is the value provided bycoef(object) andit can be changed by using the replacement method`coef<-` ifneeded.

Value

A list of results related to the Jacobians

Caution

For now the function only works whenobject has class"covMan".

Note

As a rule of thumb, a gradient coded without error gives a value oftest less than1e-4, and usually the value is muchsmaller than that.

Author(s)

Yves Deville

Check Length and Names of a Vector of Values for Parameters orBounds

Description

Check length/names for a vector of values for parameters orbounds.

Usage

checkPar(value, parN, parNames, default)

Arguments

Value

A numeric vector.

Examples

checkPar(value = c(1, 2), parN = 2L, parNames = c("theta", "sigma2"),         default = 1.0)checkPar(value = NULL, parN = 2L, parNames = c("theta", "sigma2"),         default = 1.0)checkPar(value = c("sigma2" = 100, "theta" = 1),         parN = 2L, parNames = c("theta", "sigma2"),         default = 1.0)

Generic function: Check the Compatibility of a Design Matrix with aGiven Covariance Object

Description

Generic function to check the compatibility of a design matrix with acovariance object.

Usage

   checkX(object, X, ...)

Arguments

Value

A matrix with columns taken fromX and with column namesidentical toinputNames(object).

See Also

TheinputNames method.

Check the Compatibility of a Design with a GivenCovariance Object

Description

Check the compatibility of a design matrix with a covariance object.

Usage

   ## S4 method for signature 'covAll'checkX(object, X, strict = FALSE, ...)

Arguments

Details

The matrixX must have the number of columns expected from thecovariance kernel object description, and it must have named columnsconforming to the kernel input names as returned by theinputNames method. If the two sets of names are identical butthe names are in a different order, the columns are permuted in orderto be in the same order as the input names. If the names sets differ,an error occurs.

Value

A matrix with columns names identical to the input names attached withthe kernel object, i.e.inputNames(object). The columns arecopies of those found under the same names inX, but are put inthe order ofinputNames(object). When an input name does notexist incolnames(X) an error occurs.

See Also

TheinputNames method.

Generic Function: Replacement of Coefficient Values

Description

Generic function for the replacement of coefficient values.

Usage

`coef<-`(object, ..., value)

Arguments

Value

The modified object.

Extract Coefficients of a Covariance Kernel Object as Vector, List orMatrix

Description

Extract some of or all the coefficients of a covariance kernel object asvector, list or matrix.

Usage

## S4 method for signature 'covMan'coef(object)## S4 method for signature 'covTS'coef(object, type = "all", as = "vector")

Arguments

Value

A numeric vector of coefficients or a structure as specified byas containing the coefficients selected bytype.

See Also

Thecoef<- replacement method which takes a vector ofreplacement values.

Examples

d <- 3myCov1 <- covTS(d = d, kernel = "k1Exp", dep = c(range = "input"),                value = c(range = 1.1))myCov1## versatile 'coef' methodcoef(myCov1)coef(myCov1, as = "matrix")coef(myCov1, as = "list")coef(myCov1, as = "matrix", type = "range")coef(myCov1) <- c(0.2, 0.3, 0.4, 4, 16, 25)coef(myCov1, as = "matrix")

Extract or Set Lower/Upper Bounds on Coefficients

Description

Extract or set lower/upper bounds on coefficients for covariancekernel objects.

Usage

coefLower(object, ...)coefUpper(object, ...)

Arguments

Value

The lower or upper bounds on the covariance kernel parameters.

Modified Helmert Contrast Matrix

Description

Modified Helmert contrast (or coding) matrix.

Usage

contr.helmod(n)

Arguments

Details

The returned matrix is a scaled version ofcontr.helmert(A).

Value

An orthogonal matrix withn rows andn - 1 columns. Thecolumns form a basis of the subspace orthogonal to a vector ofn ones.

Examples

A <- contr.helmod(6)crossprod(A)

Correlation Matrix for the Compound Symmetry Structure

Description

Compute the correlation matrix for a the compound symmetry structure.

Usage

corLevCompSymm(par, nlevels, levels, lowerSQRT = FALSE, compGrad = TRUE,  cov = FALSE, impl = c("C", "R"))

Arguments

Value

A correlation matrix (or its Cholesky root) with theoptionalgradient attribute.

Note

WhenlowerSQRT isFALSE, the implementationused is always in R because no gain would then result from animplementation in C.

Author(s)

Yves Deville

Examples

checkGrad <- TRUElowerSQRT <- FALSEnlevels <- 12set.seed(1234)par <- runif(1L, min = 0, max = pi)##============================================================================## Compare R and C implementations for 'lowerSQRT = TRUE'##============================================================================tR <- system.time(TR <- corLevCompSymm(nlevels = nlevels, par = par,                                       lowerSQRT = lowerSQRT, impl = "R"))tC <- system.time(T <- corLevCompSymm(nlevels = nlevels, par = par,                                      lowerSQRT = lowerSQRT))tC2 <- system.time(T2 <- corLevCompSymm(nlevels = nlevels, par = par,                                        lowerSQRT = lowerSQRT, compGrad = FALSE))## timerbind(R = tR, C = tC, C2 = tC2)## resultsmax(abs(T - TR))max(abs(T2 - TR))##===========================================================================## Compare the gradients##===========================================================================if (checkGrad) {    library(numDeriv)    ##=======================    ## lower SQRT case only    ##========================    JR <- jacobian(fun = corLevCompSymm, x = par, nlevels = nlevels,                   lowerSQRT = lowerSQRT, impl = "R", method = "complex")    J <- attr(T, "gradient")    ## redim and compare.    dim(JR) <- dim(J)    max(abs(J - JR))    nG <- length(JR)    plot(1:nG, as.vector(JR), type = "p", pch = 21, col = "SpringGreen3",         cex = 0.8, ylim = range(J, JR),         main = paste("gradient check, lowerSQRT =", lowerSQRT))    points(x = 1:nG, y = as.vector(J), pch = 16, cex = 0.6, col = "orangered")}

Correlation or Covariance Matrix for a Diagonal Structure

Description

Compute the correlation or covariance matrix for a diagonalstructure.

Usage

corLevDiag(par, nlevels, levels, lowerSQRT = FALSE, compGrad = TRUE,  cov = 0)

Arguments

Value

A correlation matrix (or its Cholesky root) with theoptionalgradient attribute.

Examples

set.seed(123)checkGrad <- TRUEnlevels <- 12sigma2 <- rexp(n = nlevels)T0 <- corLevDiag(nlevels = nlevels, par = sigma2, cov = 2)L0 <- corLevDiag(nlevels = nlevels, par = sigma2, cov = 2,                 lowerSQRT = TRUE)

Correlation Matrix for a Low-Rank Structure

Description

Compute the correlation matrix for a low-rank structure.

Usage

corLevLowRank(par, nlevels, rank, levels,              lowerSQRT = FALSE, compGrad = TRUE,              cov = 0, impl = c("C", "R"))

Arguments

Details

The correlation matrix with sizem is the general symmetriccorrelation matrix with rank\leq r wherer isgiven, as described by Rapisarda et al. It depends on(r - 1) \times (m - r / 2) / 2 parameters\theta_{ij} where the indicesi andjare such that1 \leq j < i fori \leq r or such that1 \leq j < r forr < i \leq n. The parameters\theta_{ij} are anglesand are to be taken to be in[0, 2\pi) ifj = 1 and in[0, \pi) otherwise.

Value

A correlation matrix (or its root) with the optionalgradientattribute.

Note

This function is essentially for internal use and the correspondingcorrelation or covariance kernels are created ascovQualobjects by using theq1LowRank creator.

Here the parameters\theta_{ij} are usedinrow order rather than in the column order. This order simplifies thecomputation of the gradient.

References

Francesco Rapisarda, Damanio Brigo, Fabio Mercurio (2007)."Parameterizing Correlations a Geometric Interpretation".IMA Journal of Management Mathematics,18(1):55-73.

Igor Grubišić, Raoul Pietersz(2007). "Efficient Rank Reduction of Correlation Matrices".LinearAlgebra and its Applications,422: 629-653.

See Also

Theq1LowRank creator of a corresponding kernel objectwith class"covQual", and the similarcorLevSymmfunction for the full-rank case.

Correlation Matrix for a General Symmetric Correlation Structure

Description

Compute the correlation matrix for a general symmetric correlationstructure.

Usage

corLevSymm(par, nlevels, levels, lowerSQRT = FALSE, compGrad = TRUE,           cov = 0, impl = c("C", "R"))

Arguments

Details

The correlation matrix with dimensionn is thegeneralsymmetric correlation matrix as described by Pinheiro and Bates andimplemented in thenlme package. It depends onn \times (n - 1) / 2 parameters\theta_{ij} wherethe indicesi andj are such that1 \leq j < i \leq n. The parameters\theta_{ij} areangles and are to be taken to be in[0, \pi) for aone-to-one parameterisation.

Value

A correlation matrix (or its Cholesky root) with the optionalgradient attribute.

Note

This function is essentially for internal use and the correspondingcorrelation or covariance kernels are created ascovQualobjects by using theq1Symm creator.

The parameters\theta_{ij} are usedinrow order rather than in the column order as in the reference or in thenlme package. This order simplifies the computation of thegradients.

References

Jose C. Pinheiro and Douglas M. Bates(1996). "Unconstrained Parameterizations for Variance-Covariance matrices".Statistics and Computing, 6(3) pp. 289-296.

Jose C. Pinheiro and Douglas M. Bates (2000)Mixed-EffectsModels in S and S-PLUS, Springer.

See Also

ThecorSymm correlation structure in thenlmepackage.

Examples

checkGrad <- TRUEnlevels <- 12npar <- nlevels * (nlevels - 1) / 2par <- runif(npar, min = 0, max = pi)##============================================================================## Compare R and C implementations for 'lowerSQRT = TRUE'##============================================================================tR <- system.time(TR <- corLevSymm(nlevels = nlevels,                                   par = par, lowerSQRT = TRUE, impl = "R"))tC <- system.time(T <- corLevSymm(nlevels = nlevels, par = par,                                  lowerSQRT = TRUE))tC2 <- system.time(T2 <- corLevSymm(nlevels = nlevels, par = par,                                    lowerSQRT = TRUE, compGrad = FALSE))## timerbind(R = tR, C = tC, C2 = tC2)## resultsmax(abs(T - TR))max(abs(T2 - TR))##============================================================================## Compare R and C implementations for 'lowerSQRT = FALSE'##============================================================================tR <- system.time(TRF <- corLevSymm(nlevels = nlevels, par = par,                                    lowerSQRT = FALSE, impl = "R"))tC <- system.time(TCF <- corLevSymm(nlevels = nlevels, par = par,                                    compGrad = FALSE, lowerSQRT = FALSE))tC2 <- system.time(TCF2 <- corLevSymm(nlevels = nlevels, par = par,                                      compGrad = TRUE, lowerSQRT = FALSE))rbind(R = tR, C = tC, C2 = tC2)max(abs(TCF - TRF))max(abs(TCF2 - TRF))##===========================================================================## Compare the gradients##===========================================================================if (checkGrad) {    library(numDeriv)    ##==================    ## lower SQRT case    ##==================    JR <- jacobian(fun = corLevSymm, x = par, nlevels = nlevels,                   lowerSQRT = TRUE, method = "complex", impl = "R")    J <- attr(T, "gradient")    ## redim and compare.    dim(JR) <- dim(J)    max(abs(J - JR))    nG <- length(JR)    plot(1:nG, as.vector(JR), type = "p", pch = 21, col = "SpringGreen3",         cex = 0.8, ylim = range(J, JR),         main = "gradient check, lowerSQRT = TRUE")    points(x = 1:nG, y = as.vector(J), pch = 16, cex = 0.6, col = "orangered")    ##==================    ## Symmetric case    ##==================    JR <- jacobian(fun = corLevSymm, x = par, nlevels = nlevels,                   lowerSQRT = FALSE, impl = "R", method = "complex")    J <- attr(TCF2, "gradient")    ## redim and compare.    dim(JR) <- dim(J)    max(abs(J - JR))    nG <- length(JR)    plot(1:nG, as.vector(JR), type = "p", pch = 21, col = "SpringGreen3",         cex = 0.8,         ylim = range(J, JR),         main = "gradient check, lowerSQRT = FALSE")    points(x = 1:nG, y = as.vector(J), pch = 16, cex = 0.6, col = "orangered")}

Creator for the Class"covANOVA"

Description

Creator for the class"covANOVA".

Usage

covANOVA(k1Fun1 = k1Fun1Gauss,      cov = c("corr", "homo"),      iso = 0, iso1 = 1L,      hasGrad = TRUE,      inputs = NULL,      d = NULL,      parNames,      par = NULL, parLower = NULL, parUpper = NULL,      label = "ANOVA kernel",      ...)

Arguments

Details

A ANOVA kernel on thed-dimensional Euclidean spacetakes the form

K(\mathbf{x},\,\mathbf{x}') = \delta^2 \prod_{\ell = 1}^d (1 + \tau_\ell^2 \kappa(r_\ell))

where\kappa(r) is a suitable correlation kernel for a one-dimensional input, andr_\ell is given byr_\ell := [x_\ell - x'_\ell] / \theta_\ell for\ell = 1 tod.

In this default form, the ANOVA kernel depends on2d + 1parameters: theranges\theta_\ell >0, thevariance ratios\tau_\ell^2, and the variance parameter\delta^2.

Anisotropic form uses the same range\thetafor all inputs, i.e. sets\theta_\ell = \theta for all\ell. This is obtained byusingiso = TRUE.

Acorrelation version uses\delta^2 = 1. This is obtained by usingcov = "corr". Mind that it does not correspond to a correlation kernel. Indeed, in general, the variance is equal to

K(\mathbf{x},\,\mathbf{x}) = \delta^2 \prod_{\ell = 1}^d (1 + \tau_\ell^2).

Finally, the correlation kernel\kappa(r) can depend ona "shape" parameter, e.g. have the form\kappa(r;\,\alpha). The extra shape parameter\alpha will be considered then as a parameter of theresulting ANOVA kernel, making it possible to estimate itby ML along with the range(s) and the variance.

Value

An object with class"covANOVA".

Examples

## Not run: if (require(DiceKriging)) {    ## a 16-points factorial design and the corresponding response    d <- 2; n <- 16; x <- seq(from = 0.0, to = 1.0, length.out = 4)    X <- expand.grid(x1 = x, x2 = x)    y <- apply(X, 1, DiceKriging::branin)    ## kriging model with matern5_2 covariance structure, constant    ## trend. A crucial point is to set the upper bounds!    mycov <- covANOVA(k1Fun1 = k1Fun1Matern5_2, d = 2, cov = "homo")    coefUpper(mycov) <- c(2.0, 2.0, 5.0, 5.0, 1e10)    mygp <- gp(y ~ 1, data = data.frame(X, y),               cov = mycov, multistart = 100, noise = TRUE)    nGrid <- 50; xGrid <- seq(from = 0, to = 1, length.out = nGrid)    XGrid <- expand.grid(x1 = xGrid, x2 = xGrid)    yGrid <- apply(XGrid, 1, DiceKriging::branin)    pgp <- predict(mygp, XGrid)$mean    mykm <- km(design = X, response = y)    pkm <- predict(mykm, XGrid, "UK")$mean    c("km" = sqrt(mean((yGrid - pkm)^2)),      "gp" = sqrt(mean((yGrid - pgp)^2)))    }## End(Not run)

Class"covANOVA"

Description

S4 class representing a Tensor Product (ANOVA) covariance kernel.

Objects from the Class

Objects can be created by calls of the formnew("covANOVA", ...)or by using thecovANOVA function.

Slots

k1Fun1:

Object of class"function" A function of a scalar numericvariable.

k1Fun1Char:

Object of class"character" describing the function in theslotk1Fun1.

hasGrad:

Object of class"logical". Tells if the value returned bythe functionkern1Fun has an attribute named"der"giving the derivative(s).

cov:

Object of class"integer". The value1L correspondsto a general covariance kernel. The value of0L sets the variance parameter to1, which doesnot correspond to a correlation kernel. See Section 'details' ofcovANOVA.

iso:

Object of class"integer". The value1L correspondsto an isotropic covariance, with all the inputs sharing the samerange value.

iso1:

Object of class"integer" used only when the function inthe slotk1Fun1 depends on parameters i.e. has more thanone formal argument. NOT IMPLEMENTED YET.

label:

Object of class"character". Short description of theobject.

d:

Object of class"integer". Dimension, i.e. number ofinputs.

inputNames:

Object of class"optCharacter". Names of the inputs.

parLower:

Object of class"numeric". Numeric values for the lowerbounds on the parameters. Can be-Inf.

parUpper:

Object of class"numeric". Numeric values for the upperbounds on the parameters. Can beInf.

par:

Object of class"numeric". Numeric values for theparameters. Can beNA.

kern1ParN1:

Object of class"integer". The number of parameters ink1Fun1 (such as a shape).

parN1:

Object of class"integer". Number of parameters of thefunctionkern1Fun (such as a shape).

parN:

Object of class"integer". Number of parameters for theobject. The include:direct parameters in the functionkern1Fun, ranges, and variance.

kern1ParNames:

Object of class"character". Names of thedirectparameters.

kernParNames:

Object of class"character". Names of the parameters.

Extends

Class"covAll", directly.

Methods

coef

signature(object = "covANOVA"): Get the vector of values forthe parameters.

coef<-

signature(object = "covANOVA", value = "numeric"): Set thevector of values for the parameters.

coefLower

signature(object = "covANOVA"): Get the vector oflower bounds on the parameters.

coefLower<-

signature(object = "covANOVA"): Set the vector of lowerbounds on the parameters.

coefUpper

signature(object = "covANOVA"): Get the vector of upperbounds on the parameters.

coefUpper<-

signature(object = "covANOVA"): Set the vector of upperbounds on the parameters.

covMat

signature(object = "covANOVA"): Compute the covariancematrix for given sites.

npar

signature(object = "covANOVA"): Get the number ofparameters.

scores

signature(object = "covANOVA"): Compute the scoresi.e. the derivatives w.r.t. the parameters of the contribution ofthe covariance in the log-likelihood of agp.

show

signature(object = "covANOVA"): Print or show the object.

varVec

signature(object = "covANOVA"): Compute the variancevector for given sites.

See Also

covTP.

Examples

showClass("covANOVA")

Virtual Class"covAll"

Description

Virtual class"covAll", union of classes including"covTS","covMan".

Methods

checkX

signature(object = "covAll", X = "matrix"): checks thecompatibility of a design with a given covariance object.

checkX

signature(object = "covAll", X = "data.frame"): checks thecompatibility of a design with a given covariance object.

inputNames

signature(object = "covAll"): returns the charactervector of input names.

hasGrad

signature(object = "covAll"): returns the logical slot hasGrad.

simulPar

signature(object = "covTS"): simulates random values for the parameters.

Examples

showClass("covAll")

Creator for the Class"covComp" for Composite CovarianceKernels

Description

Creator for the class "covComp" for Composite Covariance kernels.

Usage

covComp(formula, where = .GlobalEnv, topParLower = NULL,  topParUpper = NULL, trace = 0, ...)

Arguments

Details

A covariance object is built usingformula which involveskernel objects inheriting from the class"covAll" andpossibly of other scalar numeric parameters calledtopparameters. The formula can be thought of as involving thecovariance matrices rather than the kernel objects, each kernelobject sayobj being replaced bycovMat(obj, X) forsome design matrix or data frameX. Indeed, the sum or theproduct of two kernel objects lead to a covariance which is simplythe sum or product of the kernel covariances. The top parametersare considered as parameters of the covariance structure, as wellas the parameters of the covariance objects used in theformula. Their value at the creation time will be used and thuswill serve as initial value in estimation.

Value

An object with S4 class"covComp".

Caution

The class definition and its creator are toregarded as a DRAFT, many changes being necessary until a stableimplementation will be reached. The functions relating to thisclass are not for final users of GP models, but rather to thoseinterested in the conception and specification in view of a futurerelease of thekergp package.

Examples

## =========================================================================## build some kernels (with their inputNames) in the global environment## =========================================================================myCovExp3 <- kMatern(d = 3, nu = "1/2")inputNames(myCovExp3) <- c("x", "y", "z")myCovGauss2 <- kGauss(d = 2)inputNames(myCovGauss2) <- c("temp1", "temp2")k <- kMatern(d = 1)inputNames(k) <- "x"ell <- kMatern(d = 1)inputNames(ell) <- "y"tau2 <- 100sigma2 <- 4myCovComp <- covComp(formula = ~ tau2 * myCovGauss2() * myCovExp3() + sigma2 * k())myCovComp1 <- covComp(formula = ~ myCovGauss2() * myCovExp3() + k())inputNames(myCovComp)coef(myCovComp)n <- 5set.seed(1234)X <- data.frame(x = runif(n), y = runif(n), z = runif(n),                temp1 = runif(n), temp2 = runif(n))C <- covMat(myCovComp, X = X)Cg <- covMat(myCovComp, X = X, compGrad = TRUE)## Simulation: purely formal example, not meaningful.Y <- simulate(myCovComp, X = X, nsim = 100)

Class"covComp"

Description

Class"covComp" representing a composite kernel combiningseveral kernels objects inheriting from the class"covAll".

Objects from the Class

Objects can be created by calls of the formnew("covComp", ...) or by usingcovComp.

Slots

def:

Object of class"expression" defining theThis is a parsed and cleaned version of the valueof theformula formal incovComp.

covAlls:

Object of class"list" containing the kernelobjects used by the formula. The coefficients of thesekernels can be changed.

hasGrad:

Object of class"logical": can we differentiatethe kernel w.r.t. all its parameters?

label:

Object of class"character" A label attached to the kernelto describe it.

d:

Object of class"integer": dimension (or number of inputs).

parN:

Object of class"integer": number of parameters.

parNames:

Object of class"character": vector of parameter names. Itslength is in slotparN.

inputNames:

Object of class"character": names of the inputs used bythe kernel.

topParN:

Object of class"integer": number oftop parameters.

topParNames:

Object of class"character". Names of the top parameters.

topPar:

Object of class"numeric". Values of the top parameters.

topParLower:

Object of class"numeric". Lower bounds for the topparameters.

topParUpper:

Object of class"numeric". Upper bounds for the topparameters.

parsedFormula:

Object of class"list". Ugly draft for some slots to beadded in the next versions.

Extends

Class"covAll", directly.

Methods

as.list

signature(object = "covComp"): coerceobject into alist of covariance kernels, each inheriting from the virual class"covAll". This is useful e.g., to extract the coefficientsor to plot a covariance component.

checkX

signature(object = "covComp", X = "data.frame"): check thatthe inputs exist with suitable column names and suitablefactorcontent. The levels should match the prescribed levels. Returns amatrix with the input columns in the order prescribed byobject.

coef, coef<-

signature(object = "covComp"): extract or replace thevector of coefficients.

coefLower, coefUpper

signature(object = "covComp"): extract the vector of Loweror Upper bounds on the coefficients.

scores

signature(object = "covComp"): return the vector ofscores, i.e. the derivative of the log-likelihood w.r.t. theparameter vector at the current parameter values.

See Also

ThecovComp creator.

Examples

showClass("covComp")

Creator Function forcovMan Objects

Description

Creator function forcovMan objects representing a covariancekernel entered manually.

Usage

   covMan(kernel, hasGrad = FALSE, acceptMatrix = FALSE,          inputs = paste("x", 1:d, sep = ""),           d = length(inputs), parNames,          par = NULL, parLower = NULL, parUpper = NULL,          label = "covMan", ...)

Arguments

Details

The formals and the returned value of thekernel functionmust be in accordance with the value ofacceptMatrix.

• WhenacceptMatrix isFALSE, the formal argumentsx1 andx2 ofkernel are numeric vectors withlengthd. The returned result is a numeric vector of length1. The attribute named"gradient" of the returnedvalue (if provided in accordance with the value ofhasGrad)must then be a numeric vector with length equal to the number ofcovariance parameters. It must contain the derivative of thekernel valueK(\mathbf{x}_1,\,\mathbf{x}_2;\,\boldsymbol{\theta}) with respect to the parameter vector\boldsymbol{\theta}.

• WhenacceptMatrix isTRUE, the formalsx1 andx2 are matrices withd columns and withn_1andn_2 rows. The result is then a covariance matrixwithn_1 rows andn_2 columns. The gradientattribute (if provided in accordance with the value ofhasGrad) must be a list with length equal to the number ofcovariance parameters. The list element\ell must containa numeric matrix with dimension(n_1, n_2)which is the derivative of the covariance matrix w.r.t. thecovariance parameter\theta_\ell.

Note

The kernel function must be symmetric with respect to itsfirst two arguments, and it must be positive definite,which is not checked. If the function returns an objectwith a"gradient" attribute buthasGrad wasset toFALSE, the gradient willnot be usedin optimization.

The name of the class was motivated by earlier stages in thedevelopment.

Examples

myCovMan <-       covMan(         kernel = function(x1, x2, par) {          htilde <- (x1 - x2) / par[1]         SS2 <- sum(htilde^2)         d2 <- exp(-SS2)         kern <- par[2] * d2         d1 <- 2 * kern * SS2 / par[1]                     attr(kern, "gradient") <- c(theta = d1,  sigma2 = d2)         return(kern)      },      hasGrad = TRUE,          d = 1,      label = "myGauss",      parLower = c(theta = 0.0, sigma2 = 0.0),      parUpper = c(theta = Inf, sigma2 = Inf),      parNames = c("theta", "sigma2"),      par = c(NA, NA)      )      # Let us now code the same kernel in CkernCode <- "       SEXP kern, dkern;       int nprotect = 0, d;       double SS2 = 0.0, d2, z, *rkern, *rdkern;       d = LENGTH(x1);       PROTECT(kern = allocVector(REALSXP, 1)); nprotect++;       PROTECT(dkern = allocVector(REALSXP, 2)); nprotect++;       rkern = REAL(kern);       rdkern = REAL(dkern);       for (int i = 0; i < d; i++) {         z = ( REAL(x1)[i] - REAL(x2)[i] ) / REAL(par)[0];         SS2 += z * z;        }       d2 = exp(-SS2);       rkern[0] = REAL(par)[1] * d2;       rdkern[1] =  d2;        rdkern[0] =  2 * rkern[0] * SS2 / REAL(par)[0];       SET_ATTR(kern, install(\"gradient\"), dkern);       UNPROTECT(nprotect);       return kern;     "     myCovMan  ## "inline" the C function into an R function: much more efficient! ## Not run: require(inline)kernC <- cfunction(sig = signature(x1 = "numeric", x2 = "numeric",                                   par = "numeric"),                    body = kernCode)myCovMan <- covMan(kernel = kernC, hasGrad = TRUE, d = 1,                   parNames = c("theta", "sigma2"))myCovMan## End(Not run)## A kernel admitting matricial argumentsmyCov <- covMan(        kernel = function(x1, x2, par) {       # x1 : matrix of size n1xd      # x2 : matrix of size n2xd                  d <- ncol(x1)                  SS2 <- 0      for (j in 1:d){        Aj <- outer(x1[, j], x2[, j], "-")        Aj2 <- Aj^2        SS2 <- SS2 + Aj2 / par[j]^2      }      D2 <- exp(-SS2)      kern <- par[d + 1] * D2    },    acceptMatrix = TRUE,    d = 2,    label = "myGauss",    parLower = c(theta1 = 0.0, theta2 = 0.0, sigma2 = 0.0),    parUpper = c(theta1 = Inf, theta2 = Inf, sigma2 = Inf),    parNames = c("theta1", "theta2", "sigma2"),    par = c(NA, NA, NA))      coef(myCov) <- c(0.5, 1, 4)show(myCov)## computing the covariance kernel between two pointsX <- matrix(c(0, 0), ncol = 2)Xnew <- matrix(c(0.5, 1), ncol = 2)colnames(X) <- colnames(Xnew) <- inputNames(myCov)covMat(myCov, X)            ## same pointscovMat(myCov, X, Xnew)      ## two different points# computing covariances between sets of given locationsX <- matrix(c(0, 0.5, 0.7, 0, 0.5, 1), ncol = 2)t <- seq(0, 1, length.out = 3)Xnew <- as.matrix(expand.grid(t, t))colnames(X) <- colnames(Xnew) <- inputNames(myCov)covMat(myCov, X)         ## covariance matrixcovMat(myCov, X, Xnew)   ## covariances between design and new data

Class"covMan"

Description

S4 class representing a covariance kernel defined manually by a(semi-)positive definite function.

Objects from the Class

Objects can be created by callingnew("covMan", ...)or by using thecovMan function.

Slots

kernel:

object of class"function" defining the kernel (supposed tobe (semi-)positive definite).

hasGrad:

logical indicating whetherkernel returns the gradient(w.r.t. the vector of parameters) as"gradient" attributeof the result.

acceptMatrix:

logical indicating whetherkernel admits matrixarguments. Default isFALSE.

label:

object of class character, typically one or two words, used todescribe the kernel.

d:

object of class"integer", the spatial dimension or numberof inputs of the covariance.

inputNames:

object of class"character", vector of input names. Lengthd.

parLower:

,

parUpper:

object of class"numeric", vector of (possibly infinite)lower/upper bounds on parameters.

par:

object of class"numeric", numeric vector of parametervalues.

parN:

object of class"integer", total number of parameters.

kernParNames:

object of class"character", name of the kernelparameters.

Methods

coef<-

signature(object = "covMan"): replace the whole vector ofcoefficients, as required during ML estimation.

coefLower<-

signature(object = "covMan"): replacement method for lowerbounds on covMan coefficients.

coefLower

signature(object = "covMan"): extracts the numeric values ofthe lower bounds.

coef

signature(object = "covMan"): extracts the numeric valuesof the covariance parameters.

coefUpper<-

signature(object = "covMan"): replacement method for upperbounds on covMan coefficients.

coefUpper

signature(object = "covMan"): ...

covMat

signature(object = "covMan"): builds the covariance matrixor the cross covariance matrix between two sets of locations for acovMan object.

scores

signature(object = "covMan"): computes the scores(derivatives of the log-likelihood w.r.t. the covarianceparameters.

show

signature(object = "covMan"): prints in a custom format.

Note

While thecoef<- replacement method is typically intended forinternal use during likelihood maximization, thecoefLower<-andcoefUpper<- replacement methods can be used when someinformation is available on the possible values of the parameters.

Author(s)

Y. Deville, O. Roustant, D. Ginsbourger and N. Durrande.

See Also

ThecovMan function providing a creator.

Examples

showClass("covMan")

Generic Function: Covariance or Cross-CovarianceMatrix Between two Sets of Locations

Description

Generic function returning a covariance or a cross-covariance matrixbetween two sets of locations.

Usage

covMat(object, X, Xnew, ...)

Arguments

Value

A rectangular matrix withnrow(X) rows andnrow(Xnew)columns containing the covariancesK(\mathbf{x}_1, \mathbf{x}_2) for all the couples of sites\mathbf{x}_1 and\mathbf{x}_2.

Covariance Matrix for a Covariance Kernel Object

Description

Covariance matrix for a covariance kernel object.

Usage

## S4 method for signature 'covMan'covMat(object, X, Xnew, compGrad = hasGrad(object),        checkNames = NULL, index = 1L, ...)## S4 method for signature 'covTS'covMat(object, X, Xnew, compGrad = FALSE,        checkNames = TRUE, index = 1L, ...)

Arguments

Details

The covariance matrix is computed in a C program using the.Call interface. The R kernel function is evaluated within theC code usingeval.

Value

An_1 \times n_2 matrix with general elementC_{ij} := K(\mathbf{x}_{1,i},\,\mathbf{x}_{2,j};\,\boldsymbol{\theta}) whereK(\mathbf{x}_1,\,\mathbf{x}_2;\,\boldsymbol{\theta}) is the covariance kernel function.

Note

The value of the parameter\boldsymbol{\theta} can beextracted from the object with thecoef method.

Author(s)

Y. Deville, O. Roustant, D. Ginsbourger, N. Durrande.

See Also

coef method

Examples

myCov <- covTS(inputs = c("Temp", "Humid", "Press"),               kernel = "k1PowExp",               dep = c(range = "cst", shape = "cst"),               value = c(shape = 1.8, range = 1.1))n <- 100; X <- matrix(runif(n*3), nrow = n, ncol = 3)try(C1 <- covMat(myCov, X)) ## bad colnamescolnames(X) <- inputNames(myCov)C2 <- covMat(myCov, X)Xnew <- matrix(runif(n * 3), nrow = n, ncol = 3)colnames(Xnew) <- inputNames(myCov)C2 <- covMat(myCov, X, Xnew)## check with the same matrix in 'X' and 'Xnew'CMM <- covMat(myCov, X, X)CM <- covMat(myCov, X)max(abs(CM - CMM))

Warping-Based Covariance for an Ordinal Input

Description

Creator function for the classcovOrd-class

Usage

covOrd(ordered,        k1Fun1 = k1Fun1Matern5_2,        warpFun = c("norm", "unorm", "power", "spline1", "spline2"),        cov = c("corr", "homo"),        hasGrad = TRUE, inputs = "u",        par = NULL, parLower = NULL, parUpper = NULL,        warpKnots = NULL, nWarpKnots = NULL,       label = "Ordinal kernel",        intAsChar = TRUE,        ...)

Arguments

Details

Covariance kernel for qualitative ordered inputs obtained by warping.

Letu be an ordered factor with levelsu_1, \dots, u_L. Letk_1 be a 1-dimensional stationary kernel (with no or fixed parameters),F a warping function i.e. an increasing function on the interval[0,1] and\theta a scale parameter. Thenk is defined by:

k(u_i, u_j) = k_1([F(z_i) - F(z_{j})]/\theta)

wherez_1, \dots, z_L form a regular sequence from0 to1 (included). At this stage, the possible choices are:

• A distribution function (cdf) truncated to[0,1], among the Power and Normal cdfs.

• For the Normal distribution, an unnormalized version, corresponding to the restriction of the cdf on[0,1], is also implemented (warp = "unorm").

• An increasing spline of degree 1 (piecewise linear function) or 2. In this case,F is unnormalized. For degree 2, the implementation depends on scaling functions from DiceKriging package, which must be loaded here.

Notice that for unnormalizedF, we set\theta to 1, in order to avoid overparameterization.

Value

An object of classcovOrd-class, inheriting fromcovQual-class.

See Also

covOrd-class,warpFun

Examples

u <- ordered(1:6, labels = letters[1:6])myCov <- covOrd(ordered = u, cov = "homo", intAsChar = FALSE)myCovcoef(myCov) <- c(mean = 0.5, sd = 1, theta = 3, sigma2 = 2)myCovcheckX(myCov, X = data.frame(u = c(1L, 3L)))covMat(myCov, X = data.frame(u = c(1L, 3L)))myCov2 <- covOrd(ordered = u, k1Fun1 = k1Fun1Cos, warpFun = "power")coef(myCov2) <- c(pow = 1, theta = 1) myCov2plot(myCov2, type = "cor", method = "ellipse")plot(myCov2, type = "warp", col = "blue", lwd = 2)myCov3 <- covOrd(ordered = u, k1Fun1 = k1Fun1Cos, warpFun = "spline1")coef(myCov3) <- c(rep(0.5, 2), 2, rep(0.5, 2))myCov3plot(myCov3, type = "cor", method = "ellipse")plot(myCov3, type = "warp", col = "blue", lwd = 2)str(warpPower)  # details on the list describing the Power cdfstr(warpNorm)   # details on the list describing the Normal cdf

Class"covOrd"

Description

Covariance kernel for qualitative ordered inputs obtained by warping.

Letu be an ordered factor with levelsu_1, \dots, u_L. Letk_1 be a 1-dimensional stationary kernel (with no or fixed parameters),F a warping function i.e. an increasing function on the interval[0,1] and\theta a scale parameter. Thenk is defined by:

k(u_i, u_j) = k_1([F(z_i) - F(z_{j})]/\theta)

wherez_1, \dots, z_L form a regular sequence from0 to1 (included). Notice that an example of warping is a distribution function (cdf) restricted to[0,1].

Objects from the Class

Objects can be created by calls of the formnew("covOrd", ...).

Slots

covLevels:

Same as forcovQual-class.

covLevMat:

Same as forcovQual-class.

hasGrad:

Same as forcovQual-class.

acceptLowerSQRT:

Same as forcovQual-class.

label:

Same as forcovQual-class.

d:

Same as forcovQual-class. Here equal to 1.

inputNames:

Same as forcovQual-class.

nlevels:

Same as forcovQual-class.

levels:

Same as forcovQual-class.

parLower:

Same as forcovQual-class.

parUpper:

Same as forcovQual-class.

par:

Same as forcovQual-class.

parN:

Same as forcovQual-class.

kernParNames:

Same as forcovQual-class.

k1Fun1:

A function representing a 1-dimensional stationary kernel function, with no or fixed parameters.

warpFun:

A cumulative density function representing a warping.

cov:

Object of class"integer". The value0L correspondsto a correlation kernel while1L is for a covariancekernel.

parNk1:

Object of class"integer". Number of parameters ofk1Fun1. Equal to0 at this stage.

parNwarp:

Object of class"integer". Number of parameters ofwarpFun.

k1ParNames:

Object of class"character". Parameter names ofk1Fun1.

warpParNames:

Object of class"character". Parameter names ofwarpFun.

warpKnots:

Object of class"numeric". Parameters ofwarpFun.

ordered:

Object of class"logical".TRUE for an ordinal input.

intAsChar:

Object of class"logical". IfTRUE (default), an integer-valued input will be coerced into a character.Otherwise, it will be coerced into a factor.

Methods

checkX

signature(object = "covOrd", X = "data.frame"): check thatthe inputs exist with suitable column names and suitablefactorcontent. The levels should match the prescribed levels. Returns amatrix with the input columns in the order prescribed byobject.

signature(object = "covOrd", X = "matrix"): check that theinputs exist with suitable column names and suitablenumericcontent for coercion into a factor with the prescribed levels.Returns a data frame with the input columns in the orderprescribed byobject.

coef<-

signature(object = "covOrd"): replace the whole vector ofcoefficients, as required during ML estimation.

coefLower<-

signature(object = "covOrd"): replacement method for lowerbounds on covOrd coefficients.

coefLower

signature(object = "covOrd"): extracts the numeric values ofthe lower bounds.

coef

signature(object = "covOrd"): extracts the numeric valuesof the covariance parameters.

coefUpper<-

signature(object = "covOrd"): replacement method for upperbounds oncovOrd coefficients.

coefUpper

signature(object = "covOrd"): ...

covMat

signature(object = "covOrd"): build the covariance matrixor the cross covariance matrix between two sets of locations for acovOrd object.

npar

signature(object = "covOrd"): returns the number ofparameters.

scores

signature(object = "covOrd"): return the vector ofscores, i.e. the derivative of the log-likelihood w.r.t. theparameter vector at the current parameter values.

simulate

signature(object = "covOrd"): simulatensim pathsfrom a Gaussian Process having the covariance structure. The pathsare indexed by the finite set of levels of factor inputs, and theyare returned as columns of a matrix.

varVec

signature(object = "covOrd"): build the variance vectorcorresponding to a set locations for acovOrd object.

Note

This class is to be regarded as experimental. The slot names or listmay be changed in the future. The methodsnpar,inputNames or`inputNames<-` should provide a morerobust access to some slot values.

See Also

SeecovMan for a comparable structure dedicatedto kernels with continuous inputs.

Examples

showClass("covOrd")

Class"covQual"

Description

Covariance kernel for qualitative inputs.

Objects from the Class

Objects can be created by calls of the formnew("covQual", ...).

Slots

covLevels:

Object of class"function". This function hasarguments'par' and optional argumentslowerSQRT andcompGrad. It returns the covariance matrix for an input corresponding to all the levels.

covLevMat:

Object of class"matrix". This is the result returned by the functioncovLevels (former slot) withlowerSQRT = FALSEandgradient = FALSE.

hasGrad:

Object of class"logical". WhenTRUE, thecovariance matrix returned by the function in slotcovLevels must compute the gradients. The returned covariance matrix must have a"gradient" attribute;this must be an array with dimensionc(m, m, np) wherem stands for the number of levels andnp is thenumber of parameters.

acceptLowerSQRT:

Object of class"logical". WhenTRUE, the functionin slotcovLevels must have a formallowerSQRTwhich can receive a logical value. When the value isTRUEthe Cholesky (lower) root of the covariance is returned instead ofthe covariance.

label:

Object of class"character". A description of the kernelwhich will remained attached with it.

d:

Object of class"integer". The dimension or number of(qualitative) inputs of the kernel.

inputNames:

Object of class"character". The names of the (qualitative)inputs. These will be matched against the columns of a data framewhen the kernel will be evaluated.

nlevels:

Object of class"integer". A vector with lengthd giving the number of levels for each of thedinputs.

levels:

Object of class"list". A list of lengthdcontaining thed character vectors of levels forthed (qualitative) inputs.

parLower:

Object of class"numeric". Vector ofparN lowervalues for the parameters of the structure. The value-Inf can be used when needed.

parUpper:

Object of class"numeric". Vector ofparN uppervalues for the parameters of the structure. The valueInf can be used when needed.

par:

Object of class"numeric". Vector ofparNcurrent values for the structure.

parN:

Object of class"integer". Number of parameters for thestructure, as returned by thenpar method.

kernParNames:

Object of class"character". Vector of lengthparNgiving the names of the parameters. E.g."range","var","sigma2" are popular names.

ordered:

Vector of class"logical" indicating whether the factorsare ordered or not.

intAsChar:

Object of class"logical" indicating how to cope with aninteger input. WhenintAsChar isTRUE the input iscoerced into a character; the values taken by this charactervector should then match the levels in thecovQual objectas given bylevels(object)[[1]]. If insteadintAsChar isFALSE, the integer values are assumedto correspond to the levels of thecovQual object in thesame order.

Methods

checkX

signature(object = "covQual", X = "data.frame"): check thatthe inputs exist with suitable column names and suitablefactorcontent. The levels should match the prescribed levels. Returns amatrix with the input columns in the order prescribed byobject.

signature(object = "covQual", X = "matrix"): check that theinputs exist with suitable column names and suitablenumericcontent for coercion into a factor with the prescribed levels.Returns a data frame with the input columns in the orderprescribed byobject.

coef<-

signature(object = "covQual"): replace the whole vector ofcoefficients, as required during ML estimation.

coefLower<-

signature(object = "covQual"): replacement method for lowerbounds on covQual coefficients.

coefLower

signature(object = "covQual"): extracts the numeric values ofthe lower bounds.

coef

signature(object = "covQual"): extracts the numeric valuesof the covariance parameters.

coefUpper<-

signature(object = "covQual"): replacement method for upperbounds oncovQual coefficients.

coefUpper

signature(object = "covQual"): ...

covMat

signature(object = "covQual"): build the covariance matrixor the cross covariance matrix between two sets of locations for acovQual object.

npar

signature(object = "covQual"): returns the number ofparameters.

plot

signature(x = "covQual"): seeplot,covQual-method.

scores

signature(object = "covQual"): return the vector ofscores, i.e. the derivative of the log-likelihood w.r.t. theparameter vector at the current parameter values.

simulate

signature(object = "covQual"): simulatensim pathsfrom a Gaussian Process having the covariance structure. The pathsare indexed by the finite set of levels of factor inputs, and theyare returned as columns of a matrix.

varVec

signature(object = "covQual"): build the variance vectorcorresponding to a set locations for acovQual object.

Note

This class is to be regarded as experimental. The slot names or listmay be changed in the future. The methodsnpar,inputNames or`inputNames<-` should provide a morerobust access to some slot values.

See Also

SeecovMan for a comparable structure dedicatedto kernels with continuous inputs.

Examples

showClass("covQual")

Nested Qualitative Covariance

Description

Nested Qualitative Covariance

Usage

covQualNested(input = "x",              groupList = NULL,              group = NULL,              nestedLevels = NULL,              between = "Symm",              within = "Diag",              covBet = c("corr", "homo", "hete"),               covWith = c("corr", "homo", "hete"),              compGrad = TRUE,              contrasts = contr.helmod,              intAsChar = TRUE)

Arguments

Value

An object with class"covQualNested".

Caution

WhencovBet andcovWith are zero, the resulting matrixis not a correlation matrix, due to the mode ofconstruction. The "between" covariance matrix is a correlation butdiagonal blocks are added to the extended matrix obtained by re-sizingthe "between" covariance into an \times n matrix.

Note

For now the replacement method such as'coef<-' are inheritedfrom the classcovQuall. Consequently when these methods areused they do not update the covariance structure in thebetweenslot nor those in thewithin (list) slot.

This covariance kernel involvestwo categorical (i.e. factor)inputs, but these are nested. It could be aliased in the future asq1Nested orq2Nested.

Examples

### Ex 1. See the vignette "groupKernel" for an example ###       inspired from computer experiments.### Ex 2. Below an example in data analysis.country <- c("B", "B", "B", "F", "F" ,"F", "D", "D", "D")cities <- c("AntWerp", "Ghent" , "Charleroi", "Paris", "Marseille",            "Lyon", "Berlin", "Hamburg", "Munchen")myGroupList <- list(B = cities[1:3],                    F = cities[4:6],                    D = cities[7:9])## create a nested covariance. # first way, with argument 'groupList':nest1 <- covQualNested(input = "ccities",                       groupList = myGroupList,                       between = "Symm", within = "Diag",                       compGrad = TRUE,                       covBet = "corr", covWith = "corr")# second way, with arguments 'group' and 'nestedLevels'nest2 <- covQualNested(input = "ccities",                       group = country, nestedLevels = cities,                       between = "Symm", within = "Diag",                       compGrad = TRUE,                       covBet = "corr", covWith = "corr")## 'show' and 'plot' method as automatically invocatednest2plot(nest2, type = "cor")## check that the covariance matrices match for nest1 and nest2max(abs(covMat(nest1) - covMat(nest2)))## When the groups are not given in order, an error occurs!countryBad <- c("B", "B", "F", "F", "F", "D", "D", "D", "B")cities <- c("AntWerp", "Ghent", "Paris", "Marseille", "Lyon",            "Berlin", "Hamburg", "Munchen", "Charleroi")nestBad <- try(covQualNested(input = "ccities",                             group = countryBad, nestedLevels = cities,                             between = "Symm", within = "Diag",                             compGrad = TRUE,                             covBet = "corr", covWith = "corr"))

Class"covQualNested"

Description

Correlation or covariance structure for qualitative inputs(i.e. factors) obtained by nesting.

Objects from the Class

Objects can be created by calls of the formnew("covQualNested", ...).

Slots

covLevels:

Object of class"function" computing the covariance matrixfor the set of all levels.

covLevMat:

Object of class"matrix". The matrix returned by thefunction in slotcovLevels. Since this matrix is oftenneeded, it can be stored rather than recomputed.

hasGrad:

Object of class"logical". IfTRUE, the analyticalgradient can be computed.

acceptLowerSQRT:

Object of class"logical". IfTRUE, the lower squareroot of the matrix can be returned

label:

Object of class"character". A label to describe thekernel.

d:

Object of class"integer". The number of inputs.

inputNames:

Object of class"character" Names of the inputs.

nlevels:

Object of class"integer" with lengthd give thenumber of levels for the factors.

levels:

Object of class"list" with lengthd. Gives thelevels for the inputs.

parLower:

Object of class"numeric". Lower bounds on the (hyper)parameters.

parUpper:

Object of class"numeric". Upper bounds on the (hyper)parameters.

par:

Object of class"numeric". Value of the (hyper) parameters.

parN:

Object of class"integer". Number of (hyper) parameters.

kernParNames:

Object of class"character". Name of the parameters.

group:

Object of class"integer". Group numbers: one for each finallevel.

groupLevels:

Object of class"character". Vector of labels for thegroups.

between:

Object of class"covQual". A covariance or correlationstructure that can be used between groups.

within:

Object of class"list". A list of covariance or correlationstructures that are used within the groups. Each item has class"covQual".

parNCum:

Object of class"integer". Cumulated number ofparameters. Used for technical computations.

contrasts:

Object of class"function". A contrast function likecontr.helmod. This function must return a contrast matrixwith columns having unit norm.

Extends

Class"covQual", directly.Class"covAll", by class "covQual", distance 2.

Methods

No methods defined with class "covQualNested" in the signature.

Examples

showClass("covQualNested")

Creator for the Class"covRadial"

Description

Creator for the class"covRadial", which describesradialkernels.

Usage

   covRadial(k1Fun1 = k1Fun1Gauss,             cov = c("corr", "homo"),              iso = 0, hasGrad = TRUE,             inputs = NULL, d = NULL,             parNames, par = NULL,             parLower = NULL, parUpper = NULL,             label = "Radial kernel",             ...)

Arguments

Details

A radial kernel on thed-dimensional Euclidean spacetakes the form

K(\mathbf{x},\,\mathbf{x}') = \sigma^2 k_1(r)

wherek_1(r) is a suitable correlation kernel for aone-dimensional input, andr is given by

r = \left\{\sum_{\ell = 1}^d [x_\ell - x'_\ell]^2 / \theta_\ell^2 \right\}^{1/2}.

In this default form, the radial kernel depends ond + 1 parameters:theranges\theta_\ell >0 and thevariance\sigma^2.

Anisotropic form uses the same range\theta forall inputs, i.e. sets\theta_\ell = \theta for all\ell. This is obtainedby usingiso = TRUE.

Acorrelation version uses\sigma^2 = 1. Thisis obtained by usingcov = "corr".

Finally, the correlation kernelk_1(r) can depend on a"shape" parameter, e.g. have the formk_1(r;\,\alpha). The extra shape parameter\alpha will beconsidered then as a parameter of the resulting radial kernel, makingit possible to estimate it by ML along with the range(s) and thevariance.

Value

An object with class"covRadial".

Note

Whenk1Fun1 has more than one formal argument, its argumentswith position> 1 are assumed to be "shape" parameters of themodel. Examples are functions with formalsfunction(x, shape = 1.0) orfunction(x, alpha = 2.0, beta = 3.0), corresponding tovector of parameter namesc("shape") andc("alpha", "beta"). Using more than one shape parameter has not been testedyet.

Remind that using a one-dimensional correlation kernelk_1(r) heredoes not warrant that a positivesemi-definite kernel will result forany dimensiond. This question relates to Schoenberg's theorem and the conceptof completely monotone functions.

References

Gregory Fassauher and Michael McCourt (2016)Kernel-basedApproximation Methods using MATLAB. World Scientific.

See Also

k1Fun1Exp,k1Fun1Matern3_2,k1Fun1Matern5_2 ork1Fun1Gauss forexamples of functions that can be used as values for thek1Fun1formal.

Examples

set.seed(123)d <- 2; ng <- 20xg <- seq(from = 0, to = 1, length.out = ng)X <- as.matrix(expand.grid(x1 = xg, x2 = xg))## ============================================================================## A radial kernel using the power-exponential one-dimensional## function## ============================================================================d <- 2myCovRadial <- covRadial(k1Fun1 = k1Fun1PowExp, d = 2, cov = "homo", iso = 1)coef(myCovRadial)inputNames(myCovRadial) <- colnames(X)coef(myCovRadial) <- c(alpha = 1.8, theta = 2.0, sigma2 = 4.0)y <- simulate(myCovRadial, X = X, nsim = 1)persp(x = xg, y = xg, z = matrix(y, nrow = ng))## ============================================================================## Define the inverse multiquadric kernel function. We return the first two## derivatives and the gradient as attributes of the result.## ============================================================================myk1Fun <- function(x, beta = 2) {    prov <- 1 + x * x    res <- prov^(-beta)    der <- matrix(NA, nrow = length(x), ncol = 2)    der[ , 1] <- - beta * 2 * x * res / prov    der[ , 2] <- -2 * beta * (1 - (1 + 2 * beta) * x * x) * res / prov / prov    grad <- -log(prov) * res    attr(res, "gradient") <- grad    attr(res, "der") <- der    res}myCovRadial1 <- covRadial(k1Fun1 = myk1Fun, d = 2, cov = "homo", iso = 1)coef(myCovRadial1)inputNames(myCovRadial1) <- colnames(X)coef(myCovRadial1) <- c(beta = 0.2, theta = 0.4, sigma2 = 4.0)y1 <- simulate(myCovRadial1, X = X, nsim = 1)persp(x = xg, y = xg, z = matrix(y1, nrow = ng))

Class"covRadial"

Description

Class of radial covariance kernels.

Objects from the Class

Objects can be created by calls of the formcovRadial(...) ofnew("covRadial", ...).

Slots

k1Fun1:

Object of class"function" A function of a scalar numericvariable. Note that using a one-dimensional kernel heredoesnot warrant that a positive semi-definite kernel results for anydimensiond.

k1Fun1Char:

Object of class"character" describing the function in theslotk1Fun1.

hasGrad:

Object of class"logical". Tells if the value returned bythe functionkern1Fun has an attribute named"der"giving the derivative(s).

cov:

Object of class"integer". The value0L correspondsto a correlation kernel while1L is for a covariancekernel.

iso:

Object of class"integer". The value1L correspondsto an isotropic covariance, with all the inputs sharing the samerange value.

label:

Object of class"character". Short description of theobject.

d:

Object of class"integer". Dimension, i.e. number ofinputs.

inputNames:

Object of class"optCharacter". Names of the inputs.

parLower:

Object of class"numeric". Numeric values for the lowerbounds on the parameters. Can be-Inf.

parUpper:

Object of class"numeric". Numeric values for the upperbounds on the parameters. Can beInf.

par:

Object of class"numeric". Numeric values for theparameters. Can beNA.

parN1:

Object of class"integer". Number of parameters of thefunctionkern1Fun (such as a shape).

parN:

Object of class"integer". Number of parameters for theobject. The include:direct parameters in the functionkern1Fun, ranges, and variance.

kern1ParNames:

Object of class"character". Names of thedirectparameters.

kernParNames:

Object of class"character". Names of the parameters.

Extends

Class"covAll", directly.

Methods

coef<-

signature(object = "covRadial", value = "numeric"): Set thevector of values for the parameters.

coefLower<-

signature(object = "covRadial"): Set the vector of lowerbounds on the parameters.

coefLower

signature(object = "covRadial"): Get the vector oflower bounds on the parameters.

coef

signature(object = "covRadial"): Get the vector of valuesfor the parameters.

coefUpper<-

signature(object = "covRadial"): Set the vector of upperbounds on the parameters.

coefUpper

signature(object = "covRadial"): Get the vector of upperbounds on the parameters.

covMat

signature(object = "covRadial"): Compute the covariancematrix for given sites.

npar

signature(object = "covRadial"): Get the number ofparameters.

scores

signature(object = "covRadial"): Compute the scoresi.e. the derivatives w.r.t. the parameters of the contribution ofthe covariance in the log-likelihood of agp.

show

signature(object = "covRadial"): Print or show the object.

varVec

signature(object = "covRadial"): Compute the variancevector for given sites.

See Also

The creator functioncovRadial, where some details aregiven on the form of kernel.covMan andcovMan for a comparable but more general class.

Examples

showClass("covRadial")

Creator for the Class"covTP"

Description

Creator for the class"covTP".

Usage

covTP(k1Fun1 = k1Fun1Gauss,      cov = c("corr", "homo"),      iso = 0, iso1 = 1L,      hasGrad = TRUE,      inputs = NULL,      d = NULL,      parNames,      par = NULL, parLower = NULL, parUpper = NULL,      label = "Tensor product kernel",      ...)

Arguments

Details

A tensor-product kernel on thed-dimensional Euclidean spacetakes the form

K(\mathbf{x},\,\mathbf{x}') = \sigma^2 \prod_{\ell = 1}^d \kappa(r_\ell)

where\kappa(r) is a suitable correlationkernel for a one-dimensional input, andr_\ell is given byr_\ell := [x_\ell - x'_\ell] / \theta_\ell for\ell = 1 tod.

In this default form, the tensor-product kernel depends ond + 1parameters: theranges\theta_\ell >0 andthe variance\sigma^2.

Anisotropic form uses the same range\thetafor all inputs, i.e. sets\theta_\ell = \theta for all\ell. This is obtained byusingiso = TRUE.

Acorrelation version uses\sigma^2 = 1. This is obtained by usingcov = "corr".

Finally, the correlation kernel\kappa(r) can depend ona "shape" parameter, e.g. have the form\kappa(r;\,\alpha). The extra shape parameter\alpha will be considered then as a parameter of theresulting tensor-product kernel, making it possible to estimate itby ML along with the range(s) and the variance.

Value

An object with class"covTP".

Examples

## Not run: if (require(DiceKriging)) {    ## a 16-points factorial design and the corresponding response    d <- 2; n <- 16; x <- seq(from = 0.0, to = 1.0, length.out = 4)    X <- expand.grid(x1 = x, x2 = x)    y <- apply(X, 1, DiceKriging::branin)    ## kriging model with matern5_2 covariance structure, constant    ## trend. A crucial point is to set the upper bounds!    mycov <- covTP(k1Fun1 = k1Fun1Matern5_2, d = 2, cov = "homo")    coefUpper(mycov) <- c(2.0, 2.0, 1e10)    mygp <- gp(y ~ 1, data = data.frame(X, y),               cov = mycov, multistart = 100, noise = FALSE)    nGrid <- 50; xGrid <- seq(from = 0, to = 1, length.out = nGrid)    XGrid <- expand.grid(x1 = xGrid, x2 = xGrid)    yGrid <- apply(XGrid, 1, DiceKriging::branin)    pgp <- predict(mygp, XGrid)$mean    mykm <- km(design = X, response = y)    pkm <- predict(mykm, XGrid, "UK")$mean    c("km" = sqrt(mean((yGrid - pkm)^2)),      "gp" = sqrt(mean((yGrid - pgp)^2)))    }## End(Not run)

Class"covTP"

Description

S4 class representing a Tensor Product (TP) covariance kernel.

Objects from the Class

Objects can be created by calls of the formnew("covTP", ...)or by using thecovTP function.

Slots

k1Fun1:

Object of class"function" A function of a scalar numericvariable.

k1Fun1Char:

Object of class"character" describing the function in theslotk1Fun1.

hasGrad:

Object of class"logical". Tells if the value returned bythe functionkern1Fun has an attribute named"der"giving the derivative(s).

cov:

Object of class"integer". The value0L correspondsto a correlation kernel while1L is for a covariancekernel.

iso:

Object of class"integer". The value1L correspondsto an isotropic covariance, with all the inputs sharing the samerange value.

iso1:

Object of class"integer" used only when the function inthe slotk1Fun1 depends on parameters i.e. has more thanone formal argument. NOT IMPLEMENTED YET.

label:

Object of class"character". Short description of theobject.

d:

Object of class"integer". Dimension, i.e. number ofinputs.

inputNames:

Object of class"optCharacter". Names of the inputs.

parLower:

Object of class"numeric". Numeric values for the lowerbounds on the parameters. Can be-Inf.

parUpper:

Object of class"numeric". Numeric values for the upperbounds on the parameters. Can beInf.

par:

Object of class"numeric". Numeric values for theparameters. Can beNA.

kern1ParN1:

Object of class"integer". The number of parameters ink1Fun1 (such as a shape).

parN1:

Object of class"integer". Number of parameters of thefunctionkern1Fun (such as a shape).

parN:

Object of class"integer". Number of parameters for theobject. The include:direct parameters in the functionkern1Fun, ranges, and variance.

kern1ParNames:

Object of class"character". Names of thedirectparameters.

kernParNames:

Object of class"character". Names of the parameters.

Extends

Class"covAll", directly.

Methods

coef

signature(object = "covTP"): Get the vector of values forthe parameters.

coef<-

signature(object = "covTP", value = "numeric"): Set thevector of values for the parameters.

coefLower

signature(object = "covTP"): Get the vector oflower bounds on the parameters.

coefLower<-

signature(object = "covTP"): Set the vector of lowerbounds on the parameters.

coefUpper

signature(object = "covTP"): Get the vector of upperbounds on the parameters.

coefUpper<-

signature(object = "covTP"): Set the vector of upperbounds on the parameters.

covMat

signature(object = "covTP"): Compute the covariancematrix for given sites.

npar

signature(object = "covTP"): Get the number ofparameters.

scores

signature(object = "covTP"): Compute the scoresi.e. the derivatives w.r.t. the parameters of the contribution ofthe covariance in the log-likelihood of agp.

show

signature(object = "covTP"): Print or show the object.

varVec

signature(object = "covTP"): Compute the variancevector for given sites.

See Also

covRadial which is a similar covariance class andcovTP which is intended to be the standard creatorfunction for this class.

Examples

showClass("covTP")

Creator Function forcovTS Objects

Description

Creator function forcovTS objects representing a Tensor Sumcovariance kernel.

Usage

covTS(inputs = paste("x", 1:d, sep = ""),      d = length(inputs), kernel = "k1Matern5_2",      dep = NULL, value = NULL, var = 1, ...)

Arguments

Details

AcovTS object represents ad-dimensional kernel objectK of the form

K(\mathbf{x}, \mathbf{x}'; \boldsymbol{\theta}) = \sum_{i=1}^d k(x_i, x_i';\boldsymbol{\theta}_{\mathbf{s}_i})

wherek isthe covariance kernel for a Gaussian ProcessY_x indexed by a scalarx. Thed numbersx_i stand for the components of thed-dimensional location vector\mathbf{x}. The lengthp of all the vectors\mathbf{s}_i is the number ofparameters of the one-dimensional kernelk, i.e. 2 or 3 forclassical covariance kernels.

The package comes with the following covariance kernels which canbe given askernel argument.

Note that the exponential kernel ofk1Exp is identical to theMatérn kernel for\nu = 1/2, and that the three Matérns kernelsprovided here for\nu = 1/2,\nu = 3/2 and\nu = 5/2are special cases of Continuous AutoRegressive (CAR) processcovariances, with respective order1,2 and3.

Value

An object with S4 class"covTS".

Caution

The1d kernelk as given inkernel is alwaysassumed to have a variance parameter with namevar. Thisassumption may be relaxed in future versions.

Note

Most arguments receive default values or are recycled if necessary.

Author(s)

Y. Deville, O. Roustant D. Ginsbourger

References

N. Durrande, D. Ginsbourger, and O. Roustant (2012) Additive"Covariance kernels for high-dimensional Gaussian Process modeling",Annales de la Faculté des Sciences de Toulouse 21(3),pp. 481–499.

Examples

myCov1 <- covTS(kernel = "k1Exp", inputs = c("v1", "v2", "v3"),                dep = c(range = "input"))coef(myCov1) <- c(range = c(0.3, 0.7, 0.9), sigma2 = c(2, 2, 8))myCov1coef(myCov1)coef(myCov1, as = "matrix")coef(myCov1, as = "list")coef(myCov1, as = "matrix", type = "range")# with a common range parametermyCov2 <- covTS(kernel = "k1Exp", inputs = c("v1", "v2", "v3"),                 dep = c(range = "cst"), value = c(range = 0.7),                var = c(2, 2, 8))myCov2myCov3 <- covTS(d = 3, kernel = "k1PowExp",                dep = c(range = "cst", shape = "cst"),                value = c(shape = 1.8, range = 1.1),                var = c(2, 2, 8))myCov3

Class"covTS"

Description

S4 class representing a Tensor Sum (TS) covariance kernel.

Objects from the Class

Objects can be created by call of the formnew("covTS", ...) orby using thecovTS function.

Slots

d:

Object of class"integer", the spatial dimension or numberof inputs of the covariance.

inputNames:

Object of class"character", vector of input names. Lengthd.

kernel:

Object of class"covMan" representing a 1d kernel.

kernParNames:

Object of class"character", name of the kernel (among theallowed ones).

kernParCodes:

Object of class"integer", an integer code stating thedependence of the parameter to the input.

par:

Object of class"numeric", numeric vector of parametervalues.

parN:

Object of class"integer", total number of parameters.

parInput:

Object of class"integer", the number of the inputs foreach parameter. Same length aspar, values between1andd.

parLower:

,

parUpper:

Object of class"numeric" numeric, vector of (possiblyinfinite) lower/upper bounds on parameters.

parBlock:

Object of class"integer"

Methods

coef

signature(object = "covTS"): extracts the numeric values ofthe covariance parameters.

coef<-

signature(object = "covTS"): replaces the whole vector ofcoefficients, as required during ML estimation.

coefLower

signature(object = "covTS"): extracts the numeric values ofthe lower bounds.

coefLower<-

signature(object = "covTS"): replacement method for lowerbounds on covTS coefficients.

coefUpper

signature(object = "covTS"): ...

coefUpper<-

signature(object = "covTS"): replacement method for upperbounds on covTS coefficients.

covMat

signature(object = "covTS"): builds the covariance matrix,or the cross covariance matrix between two sets of locations for acovTS object.

kernelName

signature(object = "covTS"): return the character value of the kernel name.

parMap

signature(object = "covTS"): an integer matrix used to mapthecovTS parameters on the inputs and kernel parametersduring the computations.

scores

signature(object = "covTS"): computes the scores.

show

signature(object = "covTS"): prints in a custom format.

simulPar

signature(object = "covTS"): simulates random values forthe covariance parameters.

Note

The names of the methods strive to respect acamelCase namingconvention.

While thecoef<- replacement method is typically intendedfor internal use during likelihood maximization, thecoefLower<-andcoefUpper<- replacement methods can be used when somerough information exists on the possible values of the parameters.

Author(s)

Y. Deville, O. Roustant, D. Ginsbourger.

See Also

ThecovTS function providing a creator.

Examples

showClass("covTS")

Generic Function: Generalized Least Squares Estimationwith a Given Covariance Kernel

Description

Generic function computing a Generalized Least Squares estimationwith a given covariance kernel.

Usage

gls(object, ...)

Arguments

Value

A list with several elements corresponding to the estimation results.

Generalized Least Squares Estimation with a Given Covariance Kernel

Description

Generalized Least Squares (GLS) estimation for a linear model witha covariance given by the covariance kernel object. The method givesauxiliary variables as needed in many algebraic computations.

Usage

## S4 method for signature 'covAll'gls(object,    y, X, F = NULL, varNoise = NULL,     beta = NULL, checkNames = TRUE,    ...)

Arguments

Details

There are two options: for unknown trend, this is the usual GLSestimation with given covariance kernel; for a known trend, it returnsthe corresponding auxiliary variables (seevalue below).

Value

A list with several elements.

All objects having lengthp or having one of their dimensionequal top will beNULL whenF isNULL,meaning thatp = 0.

Author(s)

Y. Deville, O. Roustant

References

Kenneth Lange (2010),Numerical Analysis for Statisticians 2nd ed.pp. 102-103. Springer-Verlag,

Examples

## a possible 'covTS'myCov <- covTS(inputs = c("Temp", "Humid"),               kernel = "k1Matern5_2",               dep = c(range = "input"),               value = c(range = 0.4))d <- myCov@d; n <- 100;X <- matrix(runif(n*d), nrow = n, ncol = d)colnames(X) <- inputNames(myCov)## generate the 'GP part'  C <- covMat(myCov, X = X)L <- t(chol(C))zeta <- L %*% rnorm(n)## trend matrix 'F' for Ordinary KrigingF <- matrix(1, nrow = n, ncol = 1)varNoise <- 0.5epsilon <- rnorm(n, sd = sqrt(varNoise))beta <- 10y <- F %*% beta + zeta + epsilonfit <- gls(myCov, X = X, y = y, F = F, varNoise = varNoise)

Gaussian Process Model

Description

Gaussian Process model.

Usage

gp(formula, data, inputs = inputNames(cov), cov, estim = TRUE, ...)

Arguments

Value

A list object which is given the S3 class"gp". The list contentis very likely to change, and should be used through methods.

Note

Whenestim isTRUE, the covariance object incovis expected to provide a gradient when used to compute a covariancematrix, since the default value ofcompGrad itTRUE, seemle,covAll-method.

Author(s)

Y. Deville, D. Ginsbourger, O. Roustant

See Also

mle,covAll-method for a detailed example ofmaximum-likelihood estimation.

Examples

## ==================================================================## Example 1. Data sampled from a GP model with a known covTS object## ==================================================================set.seed(1234)myCov <- covTS(inputs = c("Temp", "Humid"),               kernel = "k1Matern5_2",               dep = c(range = "input"),               value = c(range = 0.4))## change coefficients (variances)coef(myCov) <- c(0.5, 0.8, 2, 16)d <- myCov@d; n <- 20## design matrixX <- matrix(runif(n*d), nrow = n, ncol = d)colnames(X) <- inputNames(myCov)## generate the GP realizationmyGp <- gp(formula = y ~ 1, data = data.frame(y = rep(0, n), X),             cov = myCov, estim = FALSE,            beta = 10, varNoise = 0.05)y <- simulate(myGp, cond = FALSE)$sim## parIni: add noise to true parametersparCovIni <- coef(myCov)parCovIni[] <- 0.9 * parCovIni[] +  0.1 * runif(length(parCovIni))coefLower(myCov) <- rep(1e-2, 4)coefUpper(myCov) <- c(5, 5, 20, 20)est <- gp(y ~ 1, data = data.frame(y = y, X),          cov = myCov,           noise = TRUE,          varNoiseLower = 1e-2,          varNoiseIni = 1.0,          parCovIni = parCovIni) summary(est)coef(est)## =======================================================================## Example 2. Predicting an additive function with an additive GP model## =======================================================================## Not run:         addfun6d <- function(x){       res <- x[1]^3 + cos(pi * x[2]) + abs(x[3]) * sin(x[3]^2) +           3 * x[4]^3 + 3 * cos(pi * x[5]) + 3 * abs(x[6]) * sin(x[6]^2)    }    ## 'Fit' is for the learning set, 'Val' for the validation set    set.seed(123)    nFit <- 50       nVal <- 200    d <- 6     inputs <- paste("x", 1L:d, sep = "")    ## create design matrices with DiceDesign package     require(DiceDesign)    require(DiceKriging)    set.seed(0)    dataFitIni <- DiceDesign::lhsDesign(nFit, d)$design     dataValIni <- DiceDesign::lhsDesign(nVal, d)$design     dataFit <- DiceDesign::maximinSA_LHS(dataFitIni)$design    dataVal <- DiceDesign::maximinSA_LHS(dataValIni)$design    colnames(dataFit) <- colnames(dataVal) <- inputs    testfun <- addfun6d    dataFit <- data.frame(dataFit, y = apply(dataFit, 1, testfun))    dataVal <- data.frame(dataVal, y = apply(dataVal, 1, testfun))    ## Creation of "CovTS" object with one range by input    myCov <- covTS(inputs = inputs, d = d, kernel = "k1Matern3_2",                    dep = c(range = "input"))    ## Creation of a gp object    fitgp <- gp(formula = y ~ 1, data = dataFit,                 cov = myCov, noise = TRUE,                 parCovIni = rep(1, 2*d),                parCovLower = c(rep(1e-4, 2*d)),                parCovUpper = c(rep(5, d), rep(10,d)))     predTS <- predict(fitgp, newdata = as.matrix(dataVal[ , inputs]), type = "UK")$mean    ## Classical tensor product kernel as a reference for comparison    fitRef <- DiceKriging::km(formula = ~1,                              design = dataFit[ , inputs],                              response = dataFit$y,  covtype="matern3_2")    predRef <- predict(fitRef,                       newdata = as.matrix(dataVal[ , inputs]),                       type = "UK")$mean    ## Compare TS and Ref    RMSE <- data.frame(TS = sqrt(mean((dataVal$y - predTS)^2)),                       Ref = sqrt(mean((dataVal$y - predRef)^2)),                       row.names = "RMSE")    print(RMSE)    Comp <- data.frame(y = dataVal$y, predTS, predRef)    plot(predRef ~ y, data = Comp, col = "black", pch = 4,         xlab = "True", ylab = "Predicted",         main = paste("Prediction on a validation set (nFit = ",                      nFit, ", nVal = ", nVal, ").", sep = ""))    points(predTS ~ y, data = Comp, col = "red", pch = 20)    abline(a = 0, b = 1, col = "blue", lty = "dotted")    legend("bottomright", pch = c(4, 20), col = c("black", "red"),           legend = c("Ref", "Tensor Sum"))## End(Not run)##=======================================================================## Example 3: a 'covMan' kernel with 3 implementations##=======================================================================d <- 4## -- Define a 4-dimensional covariance structure with a kernel in RmyGaussFunR <- function(x1, x2, par) {     h <- (x1 - x2) / par[1]    SS2 <- sum(h^2)    d2 <- exp(-SS2)    kern <- par[2] * d2    d1 <- 2 * kern * SS2 / par[1]                attr(kern, "gradient") <- c(theta = d1,  sigma2 = d2)    return(kern)}myGaussR <- covMan(kernel = myGaussFunR,                   hasGrad = TRUE,                   d = d,                   parLower = c(theta = 0.0, sigma2 = 0.0),                   parUpper = c(theta = Inf, sigma2 = Inf),                   parNames = c("theta", "sigma2"),                   label = "Gaussian kernel: R implementation")## -- The same, still in R, but with a kernel admitting matrices as argumentsmyGaussFunRVec <- function(x1, x2, par) {     # x1, x2 : matrices with same number of columns 'd' (dimension)    n <- nrow(x1)    d <- ncol(x1)         SS2 <- 0      for (j in 1:d){        Aj <- outer(x1[ , j], x2[ , j], "-")        Hj2 <- (Aj / par[1])^2        SS2 <- SS2 + Hj2    }    D2 <- exp(-SS2)    kern <- par[2] * D2    D1 <- 2 * kern * SS2 / par[1]     attr(kern, "gradient") <- list(theta = D1,  sigma2 = D2)      return(kern)}myGaussRVec <- covMan(    kernel = myGaussFunRVec,    hasGrad = TRUE,    acceptMatrix = TRUE,    d = d,    parLower = c(theta = 0.0, sigma2 = 0.0),    parUpper = c(theta = Inf, sigma2 = Inf),    parNames = c("theta", "sigma2"),    label = "Gaussian kernel: vectorised R implementation")## -- The same, with inlined C code## (see also another example with Rcpp by typing: ?kergp).## Not run: if (require(inline)) {    kernCode <- "       SEXP kern, dkern;       int nprotect = 0, d;       double SS2 = 0.0, d2, z, *rkern, *rdkern;       d = LENGTH(x1);       PROTECT(kern = Rf_allocVector(REALSXP, 1)); nprotect++;       PROTECT(dkern = Rf_allocVector(REALSXP, 2)); nprotect++;       rkern = REAL(kern);       rdkern = REAL(dkern);       for (int i = 0; i < d; i++) {          z = ( REAL(x1)[i] - REAL(x2)[i] ) / REAL(par)[0];          SS2 += z * z;        }       d2 = exp(-SS2);       rkern[0] = REAL(par)[1] * d2;       rdkern[1] =  d2;        rdkern[0] =  2 * rkern[0] * SS2 / REAL(par)[0];       SET_ATTR(kern, Rf_install(\"gradient\"), dkern);       UNPROTECT(nprotect);       return kern;   "    myGaussFunC <- cfunction(sig = signature(x1 = "numeric", x2 = "numeric",                                          par = "numeric"),                             body = kernCode)    myGaussC <- covMan(kernel = myGaussFunC,                       hasGrad = TRUE,                       d = d,                       parLower = c(theta = 0.0, sigma2 = 0.0),                       parUpper = c(theta = Inf, sigma2 = Inf),                       parNames = c("theta", "sigma2"),                       label = "Gaussian kernel: C/inline implementation")}## End(Not run)## == Simulate data for covMan and trend ==n <- 100; p <- d + 1X <- matrix(runif(n * d), nrow = n)colnames(X) <- inputNames(myGaussRVec)design <- data.frame(X)coef(myGaussRVec) <- myPar <- c(theta = 0.5, sigma2 = 2)myGp <- gp(formula = y ~ 1, data = data.frame(y = rep(0, n), design),             cov = myGaussRVec, estim = FALSE,            beta = 0, varNoise = 1e-8)y <- simulate(myGp, cond = FALSE)$simF <- matrix(runif(n * p), nrow = n, ncol = p)beta <- (1:p) / py <- tcrossprod(F, t(beta)) + y## == ML estimation. ==tRVec <- system.time(    resRVec <- gp(formula = y ~ ., data = data.frame(y = y, design),                  cov = myGaussRVec,                  compGrad = TRUE,                   parCovIni = c(0.5, 0.5), varNoiseLower = 1e-4,                  parCovLower = c(1e-5, 1e-5), parCovUpper = c(Inf, Inf)))summary(resRVec)coef(resRVec)pRVec <- predict(resRVec, newdata = design, type = "UK")    tAll <- tRVeccoefAll <- coef(resRVec)## compare time required by the 3 implementations## Not run:     tR <- system.time(        resR <- gp(formula = y ~ ., data = data.frame(y = y, design),                   cov = myGaussR,                   compGrad = TRUE,                    parCovIni = c(0.5, 0.5), varNoiseLower = 1e-4,                   parCovLower = c(1e-5, 1e-5), parCovUpper = c(Inf, Inf))    )    tAll <- rbind(tRVec = tAll, tR)    coefAll <- rbind(coefAll, coef(resR))    if (require(inline)) {        tC <- system.time(            resC <- gp(formula = y ~ ., data = data.frame(y = y, design),                       cov = myGaussC,                       compGrad = TRUE,                        parCovIni = c(0.5, 0.5), varNoiseLower = 1e-4,                       parCovLower = c(1e-5, 1e-5), parCovUpper = c(Inf, Inf))        )        tAll <- rbind(tAll, tC)        coefAll <- rbind(coefAll, coef(resC))    }## End(Not run)tAll## rows must be identical coefAll

Generic Function: Extract slot hasGrad of a Covariance Kernel

Description

Generic function returning the slot hasGrad of a Covariance Kernel.

Usage

hasGrad(object, ...)## S4 method for signature 'covAll'hasGrad(object, ...)

Arguments

Value

A logical indicating whether the gradient is supplied inobject (as indicated in the slot 'hasGrad').

Diagnostics for a Gaussian Process Model, Based on Leave-One-Out

Description

Cross Validation by leave-one-out for agp object.

Usage

## S3 method for class 'gp'influence(model, type = "UK", trend.reestim = TRUE, ...)

Arguments

Details

Leave-one-out (LOO) consists in computing the prediction at a designpoint when the corresponding observation is removed from the learningset (and this, for all design points). A quick version of LOO based onDubrule's formula is also implemented; It is limited to 2 cases:

• (type == "SK") & !trend.reestim and

• (type == "UK") & trend.reestim.

Value

A list composed of the following elements, wheren is the totalnumber of observations.

Warning

Only trend parameters are re-estimated when removing oneobservation. When the numbern of observations is small, there-estimated values can be far away from those obtained with theentire learning set.

Author(s)

O. Roustant, D. Ginsbourger.

References

F. Bachoc (2013), "Cross Validation and Maximum Likelihood estimations ofhyper-parameters of Gaussian processes with modelmisspecification".Computational Statistics and Data Analysis,66, 55-69doi:10.1016/j.csda.2013.03.016

N.A.C. Cressie (1993),Statistics for spatial data. Wiley seriesin probability and mathematical statistics.

O. Dubrule (1983), "Cross validation of Kriging in a uniqueneighborhood".Mathematical Geology,15, 687-699.

J.D. Martin and T.W. Simpson (2005), "Use of kriging models toapproximate deterministic computer models".AIAA Journal,43 no. 4, 853-863.

M. Schonlau (1997),Computer experiments and global optimization.Ph.D. thesis, University of Waterloo.

See Also

predict.gp,plot.gp

Generic Function: Names of the Inputs of a Covariance Kernel

Description

Generic function returning or setting the names of the inputs attachedwith a covariance kernel.

Usage

inputNames(object, ...)## S4 replacement method for signature 'covAll'inputNames(object, ...) <- value

Arguments

Value

A character vector withdistinct input names that are usede.g. in prediction.

Note

The input names are usually checked to control that they match thecolnames of a spatial design matrix. They play an important role sincein general the inputs are found among other columns of a data frame,and their order is not fixed. For instance in a data frame used asnewdata formal in thepredict method, the inputs aregenerally found at positions which differ from those in the data frameused at the creation of the object.

See Also

checkX

Predefined covMan Objects for 1D Kernels

Description

Predefined kernel Objects ascovMan objects.

Usage

k1Expk1Matern3_2k1Matern5_2k1Gauss

Format

Objects with class"covMan".

Details

These objects are provided mainly as examples. They are usedcovTS.

Examples

set.seed(1234)x <- sort(runif(40))X <- cbind(x = x)yExp <- simulate(k1Exp, nsim = 20, X = X)matplot(X, yExp, type = "l", col = "SpringGreen", ylab = "")yGauss <- simulate(k1Gauss, nsim = 20, X = X)matlines(X, yGauss, col = "orangered")title("Simulated paths from 'k1Exp' (green) and 'k1Gauss' (orange)")## ============================================================================## You can build a similar object using a creator of## 'covMan'. Although the objects 'k1Gauss' and 'myk1Gauss' differ,## they describe the same mathematical object.## ============================================================================myk1Gauss <- kGauss(d = 1)

One-Dimensional Classical Covariance Kernel Functions

Description

One-dimensional classical covariance kernel Functions.

Usage

k1FunExp(x1, x2, par)k1FunGauss(x1, x2, par)k1FunPowExp(x1, x2, par)k1FunMatern3_2(x1, x2, par)k1FunMatern5_2(x1, x2, par)k1Fun1Cos(x)k1Fun1Exp(x)k1Fun1Gauss(x)k1Fun1PowExp(x, alpha = 1.5)k1Fun1Matern3_2(x)k1Fun1Matern5_2(x)

Arguments

Details

These kernel functions are described in the Roustant et al (2012),table 1 p. 8. More details are given in chap. 4 of Rasmussen et al(2006).

Value

A matrix with a"gradient" attribute. This matrix hasn_1rows andn_2 columns wheren_1 andn_2 are thelength ofx1 andx2. Ifx1 andx2 havelength 1, the attribute is a vector of the same lengthp aspar and gives the derivative of the kernel with respect to theparameters in the same order. Ifx1 orx2 have length> 1, the attribute is an array with dimension(n_1, n_2, p).

Note

The kernel functions are coded in C through the.Call interfaceand are mainly intended for internal use. They are used by thecovTS class.

Be aware that very few checks are done (length of objects, order ofthe parameters, ...).

Author(s)

Oivier Roustant, David Ginsbourger, Yves Deville

References

C.E. Rasmussen and C.K.I. Williams (2006),Gaussian Processesfor Machine Learning, the MIT Press,doi:10.7551/mitpress/3206.001.0001

O. Roustant, D. Ginsbourger, Y. Deville (2012)."DiceKriging, DiceOptim: Two R Packages for the Analysis ofComputer Experiments by Kriging-Based Metamodeling and Optimization."Journal of Statistical Software, 51(1), 1-55.doi:10.18637/jss.v051.i01

Examples

## show the functionsn <- 300x0 <- 0x <- seq(from = 0, to = 3, length.out = n)kExpVal <- k1FunExp(x0, x, par = c(range = 1, var = 2))kGaussVal <- k1FunGauss(x0, x, par = c(range = 1, var = 2))kPowExpVal <- k1FunPowExp(x0, x, par = c(range = 1, shape = 1.5, var = 2))kMatern3_2Val <- k1FunMatern3_2(x0, x, par = c(range = 1, var = 2))kMatern5_2Val <- k1FunMatern5_2(x0, x, par = c(range = 1, var = 2))kerns <- cbind(as.vector(kExpVal), as.vector(kGaussVal), as.vector(kPowExpVal),               as.vector(kMatern3_2Val), as.vector(kMatern5_2Val))matplot(x, kerns, type = "l", main = "five 'kergp' 1d-kernels", lwd = 2)## extract gradienthead(attr(kPowExpVal, "gradient"))

Gauss (Squared-Exponential) Kernel

Description

Gauss (or squared exponential) covariance function.

Usage

kGauss(d)

Arguments

Value

An object of class"covMan" with default parameters: 1 forranges and variance values.

References

C.E. Rasmussen and C.K.I. Williams (2006),Gaussian Processesfor Machine Learning, the MIT Press,doi:10.7551/mitpress/3206.001.0001

Examples

kGauss()  # default: d = 1, nu = 5/2myGauss <- kGauss(d = 2)coef(myGauss) <- c(range = c(2, 5), sigma2 = 0.1)myGauss

Matérn Kernels

Description

Matérn kernels, obtained by plugging the Euclidian norm into a1-dimensional Matérn function.

Usage

   kMatern(d, nu = "5/2")

Arguments

Value

An object of class"covMan" with default parameters: 1 forranges and variance values.

Note

Notice that these kernels are NOT obtained by tensor product.

References

C.E. Rasmussen and C.K.I. Williams (2006),Gaussian Processesfor Machine Learning, the MIT Press,doi:10.7551/mitpress/3206.001.0001

Examples

kMatern()  # default: d = 1, nu = 5/2kMatern(d = 2)myMatern <- kMatern(d = 5, nu = "3/2")coef(myMatern) <- c(range = 1:5, sigma2 = 0.1)myMaterntry(kMatern(nu = 2))  # error

Name of the One-Dimensional Kernel in a Composite Kernel Object

Description

Name of the 1d kernel in a composite kernel object.

Usage

kernelName(object, ...)

Arguments

Value

A character string giving the kernel name.

Generic Function: Maximum Likelihood Estimation of a Gaussian ProcessModel

Description

Generic function for the Maximum Likelihood estimation of a Gaussian Processmodel.

Usage

mle(object, ...)

Arguments

Value

An estimated model, typically a list.

See Also

mle-methods for examples.

Maximum Likelihood Estimation of Gaussian Process Model Parameters

Description

Maximum Likelihood estimation of Gaussian Process models whichcovariance structure is described in a covariance kernel object.

Usage

## S4 method for signature 'covAll'mle(object,     y, X, F = NULL, beta = NULL,    parCovIni = coef(object),    parCovLower = coefLower(object),     parCovUpper = coefUpper(object),    noise = TRUE, varNoiseIni = var(y) / 10,    varNoiseLower = 0, varNoiseUpper = Inf,    compGrad = hasGrad(object),    doOptim = TRUE,    optimFun = c("nloptr::nloptr", "stats::optim"),    optimMethod = ifelse(compGrad, "NLOPT_LD_LBFGS", "NLOPT_LN_COBYLA"),    optimCode = NULL,    multistart = 1,    parTrack = FALSE, trace  = 0, checkNames = TRUE,    ...)

Arguments

Details

The choice of optimization method is as follows.

• WhenoptimFun isnloptr:nloptr, it is assumedthat we are minimizing the negative log-likelihood- \log L. Note that both predefined methods"NLOPT_LD_LBFGS" and"NLOPT_LN_COBYLA" can cope with anon-finite value of the objective, except for the initial value ofthe parameter. Non-finite values of- \log L areoften met in practice during optimization steps. The method"NLOPT_LD_LBFGS" used whencompGrad isTRUErequires that the gradient is provided by/with the covarianceobject. You can try other values ofoptimMethod correspondingto the possible choice of the"algorithm" element in theopts argument ofnloptr:nloptr. It may be useful togive other options in order to control the optimization and itsstopping rule.

• WhenoptimFun is"stats:optim", it is assumedthat we are maximizing the log-likelihood\log L. Wesuggest to use one of the methods"L-BFGS-B" or"BFGS". Notice thatcontrol can be provided in..., butcontrol$fnscale is forced to be- 1,corresponding to maximization. Note that"L-BFGS-B" uses boxconstraints, but the optimization stops as soon as thelog-likelihood is non-finite orNA. The method"BFGS"does not use constraints but allows the log-likelihood to benon-finite orNA. Both methods can be used without gradientor with gradient ifobject allows this.

The vectorsparCovIni,parCovLower,parCovUppermust have elements corresponding to those of the vector of kernelparameters given bycoef(object). These vectors should havesuitably named elements.

Value

A list with elements hopefully having understandable names.

Note

The checks concerning the parameter names, dimensions of providedobjects, ... are not fully implemented yet.

Using theoptimCode possibility requires a bit of programmingeffort, although a typical code only contains a few lines.

Author(s)

Y. Deville, O. Roustant

See Also

gp for various examples,optimMethodsto see the possible values of the argumentoptimMethod.

Examples

set.seed(29770)##=======================================================================## Example. A 4-dimensional "covMan" kernel##=======================================================================d <- 4myCovMan <-       covMan(         kernel = function(x1, x2, par) {          htilde <- (x1 - x2) / par[1]         SS2 <- sum(htilde^2)         d2 <- exp(-SS2)         kern <- par[2] * d2         d1 <- 2 * kern * SS2 / par[1]                     attr(kern, "gradient") <- c(theta = d1,  sigma2 = d2)         return(kern)      },      label = "myGauss",      hasGrad = TRUE,      d = 4,          parLower = c(theta = 0.0, sigma2 = 0.0),      parUpper = c(theta = +Inf, sigma2 = +Inf),      parNames = c("theta", "sigma2"),      par = c(NA, NA)      )kernCode <- "       SEXP kern, dkern;       int nprotect = 0, d;       double SS2 = 0.0, d2, z, *rkern, *rdkern;       d = LENGTH(x1);       PROTECT(kern = allocVector(REALSXP, 1)); nprotect++;       PROTECT(dkern = allocVector(REALSXP, 2)); nprotect++;       rkern = REAL(kern);       rdkern = REAL(dkern);       for (int i = 0; i < d; i++) {         z = ( REAL(x1)[i] - REAL(x2)[i] ) / REAL(par)[0];         SS2 += z * z;        }       d2 = exp(-SS2);       rkern[0] = REAL(par)[1] * d2;       rdkern[1] =  d2;        rdkern[0] =  2 * rkern[0] * SS2 / REAL(par)[0];       SET_ATTR(kern, install(\"gradient\"), dkern);       UNPROTECT(nprotect);       return kern;     "## inline the C function into an R function: MUCH MORE EFFICIENT!!!## Not run: require(inline)kernC <- cfunction(sig = signature(x1 = "numeric", x2 = "numeric",                                   par = "numeric"),                    body = kernCode)myCovMan <- covMan(kernel = kernC, hasGrad = TRUE, label = "myGauss", d = 4,                   parNames = c("theta", "sigma2"),                   parLower = c("theta" = 0.0, "sigma2" = 0.0),                   parUpper = c("theta" = Inf, "sigma2" = Inf))## End(Not run)##=======================================================================## Example (continued). Simulate data for covMan and trend##=======================================================================n <- 100; X <- matrix(runif(n * d), nrow = n)colnames(X) <- inputNames(myCovMan)coef(myCovMan) <- myPar <- c(theta = 0.5, sigma2 = 2)C <- covMat(object = myCovMan, X = X,            compGrad = FALSE,  index = 1L)library(MASS)set.seed(456)y <- mvrnorm(mu = rep(0, n), Sigma = C)p <- rpois(1, lambda = 4)if (p > 0) {  F <- matrix(runif(n * p), nrow = n, ncol = p)  beta <- rnorm(p)  y <- F %*% beta + y} else F <- NULLpar <- parCovIni <- c("theta" = 0.6, "sigma2" = 4)##=======================================================================## Example (continued). ML estimation. Note the 'partrack' argument##=======================================================================           est <- mle(object = myCovMan,           parCovIni = parCovIni,           y = y, X = X, F = F,           parCovLower = c(0.05, 0.05), parCovUpper = c(10, 100),           parTrack = TRUE, noise = FALSE, checkNames = FALSE)est$opt$value## change the (constrained) optimization  method## Not run: est1 <- mle(object = myCovMan,            parCovIni = parCovIni,            optimFun = "stats::optim",            optimMethod = "L-BFGS-B",            y = y, X = X, F = F,            parCovLower = c(0.05, 0.05), parCovUpper = c(10, 100),            parTrack = TRUE, noise = FALSE, checkNames = FALSE)est1$opt$value## End(Not run)##=======================================================================## Example (continued). Grid for graphical analysis##=======================================================================## Not run:     theta.grid <- seq(from = 0.1, to = 0.7, by = 0.2)    sigma2.grid <- seq(from = 0.3, to = 6, by = 0.4)    par.grid <- expand.grid(theta = theta.grid, sigma2 = sigma2.grid)    ll <- apply(as.matrix(par.grid), 1, est$logLikFun)    llmat <- matrix(ll, nrow = length(theta.grid),                    ncol = length(sigma2.grid))## End(Not run)                ##=======================================================================## Example (continued). Explore the surface ?##=======================================================================## Not run:    require(rgl)   persp3d(x = theta.grid, y = sigma2.grid, z = ll,           xlab = "theta", ylab = "sigma2", zlab = "logLik",           col = "SpringGreen3", alpha = 0.6)## End(Not run)##=======================================================================## Example (continued). Draw a contour plot for the log-lik ##                        and show iterates##=======================================================================## Not run:     contour(x = theta.grid, y = sigma2.grid, z = llmat,            col = "SpringGreen3", xlab = "theta", ylab = "sigma2",            main = "log-likelihood contours and iterates",            xlim = range(theta.grid, est$parTracked[ , 1], na.rm = TRUE),            ylim = range(sigma2.grid, est$parTracked[ , 2], na.rm = TRUE))    abline(v = est$coef.kernel[1], h = est$coef.kernel[2], lty = "dotted")    niter <- nrow(est$parTracked)    points(est$parTracked[1:niter-1, ],           col = "orangered", bg = "yellow", pch = 21, lwd = 2, type = "o")    points(est$parTracked[niter, , drop = FALSE],           col = "blue", bg = "blue", pch = 21, lwd = 2, type = "o", cex = 1.5)    ann <- seq(from = 1, to = niter, by = 5)    text(x = est$parTracked[ann, 1], y = est$parTracked[ann, 2],         labels = ann - 1L, pos = 4, cex = 0.8, col = "orangered")    points(x = myPar["theta"], y = myPar["sigma2"],           bg = "Chartreuse3", col = "ForestGreen",           pch = 22, lwd = 2, cex = 1.4)    legend("topright", legend = c("optim", "optim (last)", "true"),           pch = c(21, 21, 22), lwd = c(2, 2, 2), lty = c(1, 1, NA),           col = c("orangered", "blue", "ForestGreen"),           pt.bg = c("yellow", "blue", "Chartreuse3")) ## End(Not run)

Generic function: Number of Free Parameters in a Covariance Kernel

Description

Generic function returning the number of free parameters in acovariance kernel.

Usage

npar(object, ...)

Arguments

Value

The number of parameters.

Number of Parameters for a Covariance Kernel Object

Description

Number of parameters for a covariance kernel object.

Usage

## S4 method for signature 'covMan'npar(object, ...)## S4 method for signature 'covTS'npar(object, ...)

Arguments

Value

The number of parameters.

See Also

coef method

Optimization Methods (or Algorithms) for themleMethod

Description

Optimization methods (or algorithms) for themle method.

Usage

optimMethods(optimMethod = NULL,             optimFun = c("both", "nloptr::nloptr", "stats::optim"))

Arguments

Value

A data frame with four character columns:optimFun,optimMethod,globLoc andderNo. The columnglobLoc indicate whether the method is global ("G") orlocal ("L"). The columnderNo indicates whether themethod uses derivatives (D) or not ("N") orpossibly uses it ("P"). Only methods corresponding theoptimFun = "stats::optim" can have the value"P" forderNo. The data frame can be zero-row ifoptimMethod isgiven and no method match.

Caution

The optimization method given in the argumentoptimMethod of themle method should be compliantwith thecompGrad argument. Only a small number ofpossibilities have been tested, including the default values.

References

SeeThe NLopt website.

See Also

mle-methods,optim,nloptr.

Examples

optimMethods()optimMethods(optimMethod = "cobyla")optimMethods(optimMethod = "nelder")optimMethods(optimMethod = "BFGS")optimMethods("CMAES")

Generic Function: Map the Parameters of a Composite CovarianceKernel

Description

Map the parameter of a composite covariance kernel on the inputs andthe parameters of the 1d kernel.

Usage

parMap(object, ...)

Arguments

Value

A matrix with one row by input and one column for each of theparameters of the 1d kernel.

Map the Parameters of a Structure on the Inputs and Kernel Parameters

Description

Map the parameters of a structure on the inputs and kernel parameters.

Usage

## S4 method for signature 'covTS'parMap(object, ...)

Arguments

Value

A matrix with integer values. The rows correspond to the inputs of theobject and the columns to the1d kernel parameters.The matrix element is the number of the corresponding officialcoefficient. The same parameter of the structure can be used forseveral inputs but not (yet) for several kernel parameters. Soeach row must have different integer elements, while the sameelement can be repeated within a column.

Note

This function is for internal use only.

Examples

myCov <- covTS(d = 3, kernel = "k1Gauss",               dep = c(range = "input"), value = c(range = 1.1))parMap(myCov)

Vector of Names for the General 'Symm' Parameterisation

Description

Vector of names for the general 'Symm' parameterisation.

Usage

parNamesSymm(nlev)

Arguments

Value

Character vector of names.

Examples

parNamesSymm(nlev = 4)

Parse a Formula or Expression Describing a CompositeCovariance Kernel

Description

Parse a formula (or expression) describing a composite covariancekernel.

Usage

parseCovFormula(formula, where = .GlobalEnv, trace = 0)

Arguments

Details

The formula involves existing covariance kernel objects and mustdefine a valid kernel composition rule. For instance the sum andthe product of kernels, the convex combination of kernels areclassically used. The kernels objects are used in the formula withparentheses as is they where functions calls with no formalarguments e.g.obj( ). Non-kernel objects used in theformula must be numeric scalar parameters and are calledtopparameters. The covariance objects must exist in the environmentdefined bywhere because their slots will be used toidentify the inputs and the parameters of the composite kerneldefined by the formula.

Value

A list with the results of parsing. Although the resultscontent is easy to understand, the function is not intended to beused by the final user, and the results may change in futureversions.

Caution

Only relatively simple formulas are correctlyparsed. So use only formulas having a structure similar to one ofthose given in the examples. In case of problems, error messagesare likely to be difficult to understand.

Note

The parsing separates covariance objects from topparameters. It retrieves information about the kernel inputs andparameters from the slots. Obviously, any change in thecovariances objects after the parsing (e.g. change in theparameters names or values) will not be reported in the results ofthe parsing, so kernel any needed customization must be done priorto the parsing.

Author(s)

Yves Deville

Examples

## =========================================================================## build some kernels (with their inputNames) in the global environment## =========================================================================myCovExp3 <- kMatern(d = 3, nu = "1/2")inputNames(myCovExp3) <- c("x", "y", "z")myCovGauss2 <- kGauss(d = 2)inputNames(myCovGauss2) <- c("temp1", "temp2")k <- kMatern(d = 1)inputNames(k) <- "x"ell <- kMatern(d = 1)inputNames(ell) <- "y"## =========================================================================## Parse a formula. This formula is stupid because 'myCovGauss2'## and 'myCovExp3' should be CORRELATION kernels and not## covariance kernels to produce an identifiable model.## =========================================================================cov <- ~ tau2 * myCovGauss2() * myCovExp3() + sigma2 * k()pf <- parseCovFormula(cov, trace = 1)## =========================================================================## Parse a formula with ANOVA composition## =========================================================================cov1 <- ~ tau2 * myCovGauss2() * myCovExp3() + sigma2 * (1 + k()) * (1 + ell())pf1 <- parseCovFormula(cov1, trace = 1)

Plot for a qualitative input

Description

Plots of the covariance matrix or the correlation matrix of a qualitative input. For an ordinal factor, the warping function can also be plotted.

Usage

## S4 method for signature 'covQual'plot(x, y, type = c("cov", "cor", "warping"), ...)

Arguments

Details

Covariance / correlation plots are done with packagecorrplot if loaded, orlattice else.

See Also

covOrd.

Examples

u <- ordered(1:6, levels = letters[1:6])myCov2 <- covOrd(ordered = u, k1Fun1 = k1Fun1Cos, warpFun = "norm")coef(myCov2) <- c(mean = 0.5, sd = 0.05, theta = 0.1) plot(myCov2, type = "cor", method = "ellipse")plot(myCov2, type = "warp", col = "blue", lwd = 2)

Diagnostic Plot for the Validation of agp Object

Description

Three plots are currently available, based on theinfluenceresults: one plot of fitted values against response values, one plotof standardized residuals, and one qqplot of standardized residuals.

Usage

## S3 method for class 'gp'plot(x, y, kriging.type = "UK",    trend.reestim = TRUE, which = 1:3, ...)

Arguments

Details

The standardized residuals are defined by[y(\mathbf{x}_i) - \widehat{y}_{-i}(\mathbf{x}_i)] / \widehat{\sigma}_{-i}(\mathbf{x}_i), wherey(\mathbf{x}_i) is the response at thelocation\mathbf{x}_i,\widehat{y}_{-i}(\mathbf{x}_i) is the fittedvalue when thei-th observation is omitted (seeinfluence.gp), and\widehat{\sigma}_{-i}(\mathbf{x}_i) is thecorresponding kriging standard deviation.

Value

A list composed of the following elements wheren is the totalnumber of observations.

Warning

Only trend parameters are re-estimated when removing oneobservation. When the numbern of observations is small,re-estimated values can substantially differ from those obtained withthe whole learning set.

References

F. Bachoc (2013), "Cross Validation and Maximum Likelihood estimations ofhyper-parameters of Gaussian processes with modelmisspecification".Computational Statistics and Data Analysis,66, 55-69.

N.A.C. Cressie (1993),Statistics for spatial data. Wiley seriesin probability and mathematical statistics.

O. Dubrule (1983), "Cross validation of Kriging in a uniqueneighborhood".Mathematical Geology,15, 687-699.

J.D. Martin and T.W. Simpson (2005), "Use of kriging models toapproximate deterministic computer models".AIAA Journal,43 no. 4, 853-863.

M. Schonlau (1997),Computer experiments and global optimization.Ph.D. thesis, University of Waterloo.

See Also

predict.gp andinfluence.gp, thepredict andinfluence methods for"gp".

Plot Simulations from agp Object

Description

Function to plot simulations from agp object.

Usage

## S3 method for class 'simulate.gp'plot(x, y,        col = list(sim = "SpringGreen3", trend = "orangered"),        show = c(sim = TRUE, trend = TRUE, y = TRUE),        ...)

Arguments

Value

Nothing.

Note

For now, this function can be used only when the number of inputs isone.

See Also

simulate.gp.

Prediction Method for the"gp" S3 Class

Description

Prediction method for the"gp" S3 class.

Usage

## S3 method for class 'gp'predict(object, newdata,        type = ifelse(object$trendKnown, "SK", "UK"),         seCompute = TRUE, covCompute = FALSE,        lightReturn = FALSE, biasCorrect = FALSE,        forceInterp,                ...)

Arguments

Details

The estimated (UK) variance and covariances are NOT multiplied byn/(n-p) by default (n andp denoting the number ofrows and columns of the trend matrix\mathbf{F}). Recall thatthis correction would contribute to limit bias: it would totallyremove it if the correlation parameters were known (which is not thecase here). However, this correction is often ignored in the contextof computer experiments, especially in adaptive strategies. It can beactivated by turningbiasCorrect toTRUE, whentype = "UK"

Value

A list with the following elements.

Author(s)

O. Roustant, D. Ginsbourger, Y. Deville

See Also

gp for the creation/estimation of a model. Seegls-methods for the signification of the auxiliary variables.

Principal Kriging Functions

Description

Principal Kriging Functions.

Usage

prinKrige(object)

Arguments

Details

The Principal Kriging Functions (PKF) are the eigenvectors of asymmetric positive matrix\mathbf{B} named theBendingEnergy Matrix which is met when combining a linear trend and acovariance kernel as done ingp. This matrix hasdimensionn \times n and rankn - p. The PKF aregiven in theascending order of the eigenvaluese_i

e_1 = e_2 = \dots = e_p = 0 < e_{p+1} \leq e_{p+2} \leq \dots \leq e_n.

Thep first PKF generate the same space as do thep columns of the trend matrix\mathbf{F}, say\textrm{colspan}(\mathbf{F}). The followingn-p PKFs generate a supplementary of the subspace\textrm{colspan}(\mathbf{F}), and they have a decreasinginfluence on the response. So thep +1-th PKF can give a hint ona possible deterministic trend functions that could be added to thep existing ones.

The matrix\mathbf{B} is such that\mathbf{B} \mathbf{F} = \mathbf{0}, so the columns of\mathbf{F} can bethought of as the eigenvectors that are associated with the zeroeigenvaluese_1,\dots,e_p.

Value

A list

Note

When an eigenvaluee_i is such thate_{i-1} < e_i < e_{i+1} (which can happen only fori > p), the corresponding PKF is unique up to a change of sign. However arun ofr > 1 identical eigenvalues is associated with ar-dimensional eigenspace and the corresponding PKFs have nomeaning when they are considered individually.

References

Sahu S.K. and Mardia K.V. (2003). A Bayesian kriged Kalmanmodel for short-term forecasting of air pollution levels.Appl. Statist. 54 (1), pp. 223-244.

Examples

library(kergp)set.seed(314159)n <- 100x <- sort(runif(n))y <- 2 + 4 * x  + 2 * x^2 + 3 * sin(6 * pi * x ) + 1.0 * rnorm(n)nNew <- 60; xNew <- sort(runif(nNew))df <- data.frame(x = x, y = y)##-------------------------------------------------------------------------## use a Matern 3/2 covariance and a mispecified trend. We should guess## that it lacks a mainily linear and slightly quadratic part.##-------------------------------------------------------------------------myKern <- k1Matern3_2inputNames(myKern) <- "x"mygp <- gp(formula = y ~ sin(6 * pi * x),           data = df,            parCovLower = c(0.01, 0.01), parCovUpper = c(10, 100),           cov = myKern, estim = TRUE, noise = TRUE)PK <- prinKrige(mygp)## the third PKF suggests a possible linear trend term, and the## fourth may suggest a possible quadratic linear trendmatplot(x, PK$vectors[ , 1:4], type = "l", lwd = 2)

Qualitative Correlation or Covariance Kernel with one Input andCompound Symmetric Correlation

Description

Qualitative correlation or covariance kernel with one input andcompound symmetric correlation.

Usage

q1CompSymm(factor, input = "x", cov = c("corr", "homo"), intAsChar = TRUE)

Arguments

Value

An object with class"covQual" withd = 1qualitative input.

Note

Correlation kernels are needed in tensor products becausethe tensor product of two covariance kernels each with unknownvariance would not be identifiable.

See Also

ThecorLevCompSymm function used to computethe correlation matrix and its gradients w.r.t. the correlationparameters.

Examples

School <- factor(1L:3L, labels = c("Bad", "Mean" , "Good"))myCor <- q1CompSymm(School, input = "School")coef(myCor) <- 0.26plot(myCor, type = "cor")## Use a data.frame with a factorset.seed(246)newSchool <- factor(sample(1L:3L, size = 20, replace = TRUE),                    labels = c("Bad", "Mean" , "Good"))C1 <- covMat(myCor, X = data.frame(School = newSchool),             compGrad = FALSE, lowerSQRT = FALSE)

Qualitative Correlation or Covariance Kernel with one Input andDiagonal Structure

Description

Qualitative correlation or covariance kernel with one input and diagonalstructure.

Usage

q1Diag(factor, input = "x", cov = c("corr", "homo", "hete"), intAsChar = TRUE)

Arguments

Value

An object with class"covQual" withd = 1 qualitativeinput.

Note

The correlation version obtained withcov = "corr" has noparameters.

See Also

q1Symm,q1CompSymm are other covariancestructures for one qualitative input.

Examples

School <- factor(1L:3L, labels = c("Bad", "Mean" , "Good"))## correlation: no parameter!myCor <- q1Diag(School, input = "School")## covariance myCov <- q1Diag(School, input = "School", cov = "hete")coef(myCov) <- c(1.1, 2.2, 3.3)

Qualitative Correlation or Covariance Kernel with one Input andLow-Rank Correlation

Description

Qualitative correlation or covariance kernel with one input andlow-rank correlation.

Usage

q1LowRank(factor, rank = 2L, input = "x",          cov = c("corr", "homo", "hete"), intAsChar = TRUE)

Arguments

Details

The correlation structure involves(r - 1)(m - r /2) parameters.The parameterization of Rapisarda et al is used: the correlationparameters are angles\theta_{i,j} corresponding to1 < i \leq r and1 \leq j < ior tor < i \leq m and1 \leq j < r. Thecorrelation matrix\mathbf{C} for the levels, with sizem, factors as\mathbf{C} = \mathbf{L}\mathbf{L}^\topwhere\mathbf{L} is a lower-triangularmatrix with dimensionm \times r with all its rowshaving unit Euclidean norm. Note that the diagonal elements of\mathbf{L} can be negative and correspondingly the angles\theta_{i,1} are taken in the interval[0, 2\pi) for1 < i \leq r. Thematrix\mathbf{L} is not unique. As explained in Grubišić andPietersz, the parameterization is surjective: any correlation withrank\leq r is obtained by choosing a suitable vectorof parameters, but this vector is not unique.

Correlation kernels are needed in tensor products because the tensorproduct of two covariance kernels each with unknown variance would notbe identifiable.

Value

An object with class"covQual" withd = 1 qualitativeinput.

References

Francesco Rapisarda, Damanio Brigo, Fabio Mercurio (2007)."Parameterizing Correlations a Geometric Interpretation".IMA Journal of Management Mathematics,18(1):55-73.

Igor Grubišić, Raoul Pietersz(2007). "Efficient Rank Reduction of Correlation Matrices".LinearAlgebra and its Applications,422: 629-653.

See Also

Theq1Symm function to create a kernel object for thefull-rank case andcorLevLowRank for the correlationfunction.

Examples

myFact <- factor(letters[1:8])myCov <- q1LowRank(factor = myFact, rank = 3)## corrplotplot(myCov)## find the rank using a pivoted Choleskychol(covMat(myCov), pivot = TRUE)

Qualitative Correlation or Covariance Kernel with one Inputand General Symmetric Correlation

Description

Qualitative correlation or covariance kernel with one input andgeneral symmetric correlation.

Usage

  q1Symm(factor, input = "x", cov = c("corr", "homo", "hete"), intAsChar = TRUE)

Arguments

Value

An object with class"covQual" withd = 1 qualitativeinput.

Note

Correlation kernels are needed in tensor products because the tensorproduct of two covariance kernels each with unknown variance would notbe identifiable.

See Also

ThecorLevSymm function used to compute thecorrelation matrix and its gradients w.r.t. the correlationparameters.

Examples

School <- factor(1L:3L, labels = c("Bad", "Mean" , "Good"))myCor <- q1Symm(School, input = "School")coef(myCor) <- c(theta_2_1 = pi / 3, theta_3_1 = pi / 4, theta_3_2 = pi / 8)plot(myCor, type = "cor")## Use a data.frame with a factorset.seed(246)newSchool <- factor(sample(1L:3L, size = 20, replace = TRUE),                    labels = c("Bad", "Mean" , "Good"))C1 <- covMat(myCor, X = data.frame(School = newSchool),             compGrad = FALSE, lowerSQRT = FALSE)

Generic Function: Scores for a Covariance Kernel Object

Description

Generic function returning the scores for a covariance kernel object.

Usage

scores(object, ...)

Arguments

Details

Compute the derivatives\partial_{\theta_k}\ellfor the (possibly concentrated) log-likelihood\ell := \log L of a covariance object with parameter vector\boldsymbol{\theta}. The score for\theta_k is obtained as a matrix scalar product

\partial_{\theta_k} \ell = \textrm{trace}(\mathbf{W} \mathbf{D})

where\mathbf{D} := \partial_{\theta_k} \mathbf{C}and where\mathbf{W} is the matrix \mathbf{W} := \mathbf{e}\mathbf{e}^\top - \mathbf{C}^{-1}.The vector\mathbf{e} is the vector of residualsand the matrix\mathbf{C}is the covariance computed for the design\mathbf{X}.

Value

A numeric vector of lengthnpar(object) containing the scores.

Note

The scores can be efficiently computed when the matrix\mathbf{W} has already been pre-computed.

Extracts the Slots of a Structure

Description

Extract the slot of a structure.

Usage

shapeSlot(object, slotName = "par", type = "all", as = "vector")

Arguments

Value

A vector, list or matrix containing the extraction.

Note

This function is for internal use only.

Generic function: Draw Random Values for the Parameters of aCovariance Kernel

Description

Generic function to draw random values for the parameters of acovariance kernel object.

Usage

  simulPar(object, nsim = 1L, seed = NULL, ...)

Arguments

Details

Draw random values for the parameters of a covariance kernel objectusing the informationscoefLower andcoefUpper.

Value

A matrix withnsim rows andnpar(object) columns.

Draw Random Values for the Parameters of a Covariance Kernel

Description

Draw random values for the parameters of a covariance kernel

object.

Usage

## S4 method for signature 'covAll'simulPar(object, nsim = 1L, seed = NULL)

Arguments

Details

Draw random values for the parameters of a covariance kernelobject using the informationscoefLower andcoefUpper.

Value

A matrix withnsim rows andnpar(object) columns.

Simulation of acovAll Object

Description

Simulation of acovAll object.

Usage

## S4 method for signature 'covAll'simulate(object, nsim = 1, seed = NULL,         X, mu = NULL, method = "mvrnorm", checkNames = TRUE,         ...)

Arguments

Value

A numeric matrix withnrow(X) rows andnsim columns.Each column is the vector of the simulated path at the simulationlocations.

Note

The simulation is unconditional.

See Also

Themvrnorm function.

Examples

## -- as in example(kergp) define an argumentwise invariant kernel --kernFun <- function(x1, x2, par) {  h <- (abs(x1) - abs(x2)) / par[1]  S <- sum(h^2)  d2 <- exp(-S)  K <- par[2] * d2  d1 <- 2 * K * S / par[1]     attr(K, "gradient") <- c(theta = d1,  sigma2 = d2)  return(K)}covSymGauss <- covMan(kernel = kernFun,                      hasGrad = TRUE,                      label = "argumentwise invariant",                      d = 2,                      parNames = c("theta", "sigma2"),                      par = c(theta = 0.5, sigma2 = 2))## -- simulate a path from the corresponding GP --nGrid <- 24; n <- nGrid^2; d <- 2xGrid <- seq(from = -1, to = 1, length.out = nGrid)Xgrid <- expand.grid(x1 = xGrid, x2 = xGrid)ySim <- simulate(covSymGauss, X = Xgrid)contour(x = xGrid, y = xGrid,        z = matrix(ySim, nrow = nGrid, ncol = nGrid),         nlevels = 15)

Simulation of Paths from agp Object

Description

Simulation of paths from agp object.

Usage

## S3 method for class 'gp'simulate(object, nsim = 1L, seed = NULL,         newdata = NULL,         cond = TRUE,         trendKnown = FALSE,         newVarNoise = NULL,         nuggetSim = 1e-8,         checkNames = TRUE,         output = c("list", "matrix"),         label = "y", unit = "",         ...)

Arguments

Value

A matrix with the simulated paths as its columns or a more completelist with more results. This list which is given the S3 class"simulate.gp" has the following elements.

Note

WhenbetaKnown isFALSE, thetrend and thesmooth GP parts of a simulation are usually correlated, andtheir sum will show less dispersion than each of the twocomponents. The covariance of the vector\widehat{\boldsymbol{\beta}} can be regarded as theposterior distribution corresponding to a non-informative prior, thedistribution from which a new path is drawn being the predictivedistribution.

Author(s)

Yves Deville

Examples

set.seed(314159)n <- 40x <- sort(runif(n))y <- 2 + 4 * x  + 2 * x^2 + 3 * sin(6 * pi * x ) + 1.0 * rnorm(n)df <- data.frame(x = x, y = y)##-------------------------------------------------------------------------## use a Matern 3/2 covariance. With model #2, the trend is mispecified,## so the smooth GP part of the model captures a part of the trend.##-------------------------------------------------------------------------myKern <- k1Matern3_2inputNames(myKern) <- "x"mygp <- list()mygp[[1]] <- gp(formula = y ~ x + I(x^2) + sin(6 * pi * x), data = df,                 parCovLower = c(0.01, 0.01), parCovUpper = c(10, 100),                cov = myKern, estim = TRUE, noise = TRUE)mygp[[2]] <- gp(formula = y ~ sin(6 * pi * x), data = df,                 parCovLower = c(0.01, 0.01), parCovUpper = c(10, 100),                cov = myKern, estim = TRUE, noise = TRUE)##-------------------------------------------------------------------------## New data##-------------------------------------------------------------------------nNew <- 150xNew <- seq(from = -0.2, to= 1.2, length.out = nNew)dfNew <- data.frame(x = xNew)opar <- par(mfrow = c(2L, 2L))nsim <- 40for (i in 1:2) {    ##--------------------------------------------------------------------    ## beta known or not, conditional    ##--------------------------------------------------------------------    simTU <- simulate(object = mygp[[i]], newdata = dfNew,  nsim = nsim,                      trendKnown = FALSE)    plot(simTU, main = "trend unknown, conditional")    simTK <- simulate(object = mygp[[i]], newdata = dfNew, nsim = nsim,                      trendKnown = TRUE)    plot(simTK, main = "trend known, conditional")    ##--------------------------------------------------------------------    ## The same but UNconditional    ##--------------------------------------------------------------------    simTU <- simulate(object = mygp[[i]], newdata = dfNew,  nsim = nsim,                     trendKnown = FALSE, cond = FALSE)    plot(simTU, main = "trend unknown, unconditional")    simTK <- simulate(object = mygp[[i]], newdata = dfNew, nsim = nsim,                      trendKnown = TRUE, cond = FALSE)    plot(simTK, main = "trend known, unconditional")}par(opar)

Vector of Indices Useful for Symmetric or Anti-Symmetric Matrices.

Description

Vector of indices useful for symmetric or anti-symmetric matrices

Usage

symIndices(n, diag = FALSE)

Arguments

Details

This function is intended to provide computations which are faster thanlower.tri andupper.tri.

Value

A list containing the following integer vectors, each with length(n - 1) n / 2.

Examples

n <- rpois(1, lambda = 10)L <- symIndices(n)X <- matrix(1L:(n * n), nrow = n)max(abs(X[lower.tri(X, diag = FALSE)] - L$kL))max(abs(t(X)[lower.tri(X, diag = FALSE)] - L$kU))cbind(row = L$i, col = L$j)

Make Translucent colors

Description

Make translucent colors.

Usage

  translude(colors, alpha = 0.6)

Arguments

Value

A vector of translucent (or semi-transparent) colors.

Generic Function: Variance of Gaussian Process at Specific Locations

Description

Generic function returning a variance vector

Usage

varVec(object, X, ...)

Arguments

Value

A numeric vector with lengthnrow(X) containing the variancesK(\mathbf{x}, \mathbf{x}) for all the sites\mathbf{x}.

Covariance Matrix for a Covariance Kernel Object

Description

Covariance matrix for a covariance kernel object.

Usage

## S4 method for signature 'covMan'varVec(object, X, compGrad = FALSE,        checkNames = NULL, index = -1L, ...)## S4 method for signature 'covTS'varVec(object, X, compGrad = FALSE,        checkNames = TRUE, index = -1L, ...)

Arguments

Details

The variance vector is computed in a C program using the.Callinterface. The R kernel function is evaluated within the C code usingeval.

Value

A vector of lengthnrow(X) with general elementV_{i} := K(\mathbf{x}_{i},\,\mathbf{x}_{i};\,\boldsymbol{\theta}) whereK(\mathbf{x}_1,\,\mathbf{x}_2;\,\boldsymbol{\theta}) is the covariance kernel function.

Note

The value of the parameter\boldsymbol{\theta} can beextracted from the object with thecoef method.

See Also

coef method

Examples

myCov <- covTS(inputs = c("Temp", "Humid", "Press"),               kernel = "k1PowExp",               dep = c(range = "cst", shape = "cst"),               value = c(shape = 1.8, range = 1.1))n <- 100; X <- matrix(runif(n*3), nrow = n, ncol = 3)try(V1 <- varVec(myCov, X = X)) ## bad colnamescolnames(X) <- inputNames(myCov)V2 <- varVec(myCov, X = X)Xnew <- matrix(runif(n * 3), nrow = n, ncol = 3)colnames(Xnew) <- inputNames(myCov)V2 <- varVec(myCov, X = X)

Warpings for Ordinal Inputs

Description

Given warpings for ordinal inputs.

Usage

warpNormwarpUnormwarpPowerwarpSpline1warpSpline2

Format

The format is a list of 6:

$ fun : the warping function. The second argument is the vector ofparameters. The function returns a numeric vector with an attribute"gradient" giving the derivative w.r.t. the parameters.

$ parNames : names of warping parameters (character vector).

$ parDefault: default values of warping parameters (numeric vector).

$ parLower : lower bounds of warping parameters (numeric vector).

$ parUpper : upper bounds of warping parameters (numeric vector).

$ hasGrad : a boolean equal toTRUE ifgradient issupplied as an attribute offun.

Details

SeecovOrd for the definition of a warping in thiscontext. At this stage, two warpings corresponding to cumulativedensity functions (cdf) are implemented:

• Normal distribution, truncated to[0,1]:

  F(x) = [N(x) - N(0)] / [N(1) - N(0)]

  whereN(x) = \Phi([x - \mu] / \sigma) is the cdf of the normaldistribution with mean\mu and standard deviation\sigma.

• Power distribution on[0, 1]:F(x) = x^{pow}.

Furthermore, a warping corresponding to unnormalized Normal cdf is implemented, as well as spline warpings of degree 1 and 2. Splines are defined by a sequence ofk knots between 0 and 1. The first knot is 0, and the last is 1.A spline warping of degree 1 is a continuous piecewise linear function. It is parameterized by a positive vector of lengthk-1, representing the increments at knots.A spline warping of degree 2 is a non-decreasing quadratic spline. It is obtained by integrating a spline of degree 1.Its parameters form a positive vector of lengthk, representing the derivatives at knots. The implementation relies on the functionscalingFun1d ofDiceKriging package.