Package 'kergp'

Description

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

Details

Warning

As a lab, kergp may strongly evolve in its future life. Users interested in stable software for the Analysis of Computer Experiments are encouraged to use other packages such as DiceKriging instead.

Note

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

Author(s)

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

Maintainer: Olivier Roustant, <[email protected]>

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 and Analysis. 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 of Computer 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 <- 2
xGrid <- 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 <- 25
tab <- 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) -- 
covSymGauss
symgp <- 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)

Description

Coerce a covTP object representing a Tensor-Product covariance kernel on the $d$-dimensional Euclidean space into a list containing $d$ one-dimensional kernels.

Usage

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

Arguments

Value

A list with length d or d + 1 where d is the "dimension" slot x@d of the object x. The first d elements of the list are one-dimensional correlation kernel objects with class "covTP". When x is a covariance kernel (as opposed to a correlation kernel), the list contains one more element which gives the variance.

Caution

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

See Also

covTP and covTP-class.

Examples

set.seed(123)
d <- 6
myCov1 <- 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 <- 10
X <- 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))))
}

Description

Check the gradient provided in a covMan 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 matrices x1 and x2 with n1 and n2 rows can be given or instead be drawn at random. The matrix of kernel values with dimension c(n1, n2) is computed, together with its gradient with dimension c(n1, n2, npar) where npar is the number of parameters of the kernel. A numerical differentiation w.r.t. the kernel parameters is performed for the kernel value at x1 and x2, and the result is compared to that provided by the kernel function (the function described in the slot named "kernel" of object). Note that the value of the parameter vector is the value provided by coef(object) and it can be changed by using the replacement method `coef<-` if needed.

Value

A list of results related to the Jacobians

Caution

For now the function only works when object has class "covMan".

Note

As a rule of thumb, a gradient coded without error gives a value of test less than 1e-4, and usually the value is much smaller than that.

Author(s)

Yves Deville

Description

Check length/names for a vector of values for parameters or bounds.

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)

Description

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

Usage

checkX(object, X, ...)

Arguments

Value

A matrix with columns taken from X and with column names identical to inputNames(object).

See Also

The inputNames method.

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 matrix X must have the number of columns expected from the covariance kernel object description, and it must have named columns conforming to the kernel input names as returned by the inputNames method. If the two sets of names are identical but the names are in a different order, the columns are permuted in order to 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 with the kernel object, i.e. inputNames(object). The columns are copies of those found under the same names in X, but are put in the order of inputNames(object). When an input name does not exist in colnames(X) an error occurs.

See Also

The inputNames method.

Description

Extract some of or all the coefficients of a covariance kernel object as vector, 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 by as containing the coefficients selected by type.

See Also

The coef<- replacement method which takes a vector of replacement values.

Examples

d <- 3
myCov1 <- covTS(d = d, kernel = "k1Exp", dep = c(range = "input"),
                value = c(range = 1.1))
myCov1
## versatile 'coef' method
coef(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")

Description

Generic function for the replacement of coefficient values.

Usage

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

Arguments

Value

The modified object.

Description

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

Usage

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

Arguments

Value

The lower or upper bounds on the covariance kernel parameters.

Description

Modified Helmert contrast (or coding) matrix.

Usage

contr.helmod(n)

Arguments

Details

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

Value

An orthogonal matrix with n rows and n - 1 columns. The columns form a basis of the subspace orthogonal to a vector of n ones.

Examples

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

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 the optional gradient attribute.

Note

When lowerSQRT is FALSE, the implementation used is always in R because no gain would then result from an implementation in C.

Author(s)

Yves Deville

Examples

checkGrad <- TRUE
lowerSQRT <- FALSE
nlevels <- 12
set.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))
## time
rbind(R = tR, C = tC, C2 = tC2)

## results
max(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")
}

Description

Compute the correlation or covariance matrix for a diagonal structure.

Usage

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

Arguments

Value

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

Examples

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

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 size $m$ is the general symmetric correlation matrix with rank $\leq r$ where $r$ is given, as described by Rapisarda et al. It depends on $(r - 1) \times (m - r / 2) / 2$ parameters $\theta_{ij}$ where the indices $i$ and $j$ are such that $1 \leq j < i$ for $i \leq r$ or such that $1 \leq j < r$ for $r < i \leq n$. The parameters $\theta_{ij}$ are angles and are to be taken to be in $[0, 2\pi)$ if $j = 1$ and in $[0, \pi)$ otherwise.

Value

A correlation matrix (or its root) with the optional gradient attribute.

Note

This function is essentially for internal use and the corresponding correlation or covariance kernels are created as covQual objects by using the q1LowRank creator.

Here the parameters $\theta_{ij}$ are used in row order rather than in the column order. This order simplifies the computation 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". Linear Algebra and its Applications, 422: 629-653.

See Also

The q1LowRank creator of a corresponding kernel object with class "covQual", and the similar corLevSymm function for the full-rank case.

Description

Compute the correlation matrix for a general symmetric correlation structure.

Usage

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

Arguments

Details

The correlation matrix with dimension $n$ is the general symmetric correlation matrix as described by Pinheiro and Bates and implemented in the nlme package. It depends on $n \times (n - 1) / 2$ parameters $\theta_{ij}$ where the indices $i$ and $j$ are such that $1 \leq j < i \leq n$. The parameters $\theta_{ij}$ are angles and are to be taken to be in $[0, \pi)$ for a one-to-one parameterisation.

Value

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

Note

This function is essentially for internal use and the corresponding correlation or covariance kernels are created as covQual objects by using the q1Symm creator.

The parameters $\theta_{ij}$ are used in row order rather than in the column order as in the reference or in the nlme package. This order simplifies the computation of the gradients.

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-Effects Models in S and S-PLUS, Springer.

See Also

The corSymm correlation structure in the nlme package.

Examples

checkGrad <- TRUE
nlevels <- 12
npar <- nlevels * (nlevels - 1) / 2
par <- 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))
## time
rbind(R = tR, C = tC, C2 = tC2)

## results
max(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")
}

Description

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

Methods

checkX

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

checkX

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

inputNames

signature(object = "covAll"): returns the character vector 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")

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 the $d$-dimensional Euclidean space takes 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, and $r_\ell$ is given by $r_\ell := [x_\ell - x'_\ell] / \theta_\ell$ for $\ell = 1$ to $d$.

In this default form, the ANOVA kernel depends on $2d + 1$ parameters: the ranges $\theta_\ell >0$, the variance ratios $\tau_\ell^2$, and the variance parameter $\delta^2$.

An isotropic form uses the same range $\theta$ for all inputs, i.e. sets $\theta_\ell = \theta$ for all $\ell$. This is obtained by using iso = TRUE.

A correlation version uses $\delta^2 = 1$. This is obtained by using cov = "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 on a "shape" parameter, e.g. have the form $\kappa(r;\,\alpha)$. The extra shape parameter $\alpha$ will be considered then as a parameter of the resulting ANOVA kernel, making it possible to estimate it by 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)

Description

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

Objects from the Class

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

Slots

k1Fun1:

Object of class "function" A function of a scalar numeric variable.

k1Fun1Char:

Object of class "character" describing the function in the slot k1Fun1.

hasGrad:

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

cov:

Object of class "integer". The value 1L corresponds to a general covariance kernel. The value of 0L sets the variance parameter to 1, which does not correspond to a correlation kernel. See Section 'details' of covANOVA.

iso:

Object of class "integer". The value 1L corresponds to an isotropic covariance, with all the inputs sharing the same range value.

iso1:

Object of class "integer" used only when the function in the slot k1Fun1 depends on parameters i.e. has more than one formal argument. NOT IMPLEMENTED YET.

label:

Object of class "character". Short description of the object.

d:

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

inputNames:

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

parLower:

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

parUpper:

Object of class "numeric". Numeric values for the upper bounds on the parameters. Can be Inf.

par:

Object of class "numeric". Numeric values for the parameters. Can be NA.

kern1ParN1:

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

parN1:

Object of class "integer". Number of parameters of the function kern1Fun (such as a shape).

parN:

Object of class "integer". Number of parameters for the object. The include: direct parameters in the function kern1Fun, ranges, and variance.

kern1ParNames:

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

kernParNames:

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

Extends

Class "covAll", directly.

Methods

coef

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

coef<-

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

coefLower

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

coefLower<-

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

coefUpper

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

coefUpper<-

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

covMat

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

npar

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

scores

signature(object = "covANOVA"): Compute the scores i.e. the derivatives w.r.t. the parameters of the contribution of the covariance in the log-likelihood of a gp.

show

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

varVec

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

See Also

covTP.

Examples

showClass("covANOVA")

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 using formula which involves kernel objects inheriting from the class "covAll" and possibly of other scalar numeric parameters called top parameters. The formula can be thought of as involving the covariance matrices rather than the kernel objects, each kernel object say obj being replaced by covMat(obj, X) for some design matrix or data frame X. Indeed, the sum or the product of two kernel objects lead to a covariance which is simply the sum or product of the kernel covariances. The top parameters are considered as parameters of the covariance structure, as well as the parameters of the covariance objects used in the formula. Their value at the creation time will be used and thus will serve as initial value in estimation.

Value

An object with S4 class "covComp".

Caution

The class definition and its creator are to regarded as a DRAFT, many changes being necessary until a stable implementation will be reached. The functions relating to this class are not for final users of GP models, but rather to those interested in the conception and specification in view of a future release of the kergp 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 <- 100
sigma2 <- 4

myCovComp <- covComp(formula = ~ tau2 * myCovGauss2() * myCovExp3() + sigma2 * k())

myCovComp1 <- covComp(formula = ~ myCovGauss2() * myCovExp3() + k())

inputNames(myCovComp)
coef(myCovComp)

n <- 5
set.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)

Description

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

Objects from the Class

Objects can be created by calls of the form new("covComp", ...) or by using covComp.

Slots

def:

Object of class "expression" defining the This is a parsed and cleaned version of the value of the formula formal in covComp.

covAlls:

Object of class "list" containing the kernel objects used by the formula. The coefficients of these kernels can be changed.

hasGrad:

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

label:

Object of class "character" A label attached to the kernel to 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. Its length is in slot parN.

inputNames:

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

topParN:

Object of class "integer": number of top 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 top parameters.

topParUpper:

Object of class "numeric". Upper bounds for the top parameters.

parsedFormula:

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

Extends

Class "covAll", directly.

Methods

as.list

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

checkX

signature(object = "covComp", X = "data.frame"): check that the inputs exist with suitable column names and suitable factor content. The levels should match the prescribed levels. Returns a matrix with the input columns in the order prescribed by object.

coef, coef<-

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

coefLower, coefUpper

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

scores

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

See Also

The covComp creator.

Examples

showClass("covComp")

Description

Creator function for covMan objects representing a covariance kernel 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 the kernel function must be in accordance with the value of acceptMatrix.

• When acceptMatrix is FALSE, the formal arguments x1 and x2 of kernel are numeric vectors with length $d$. The returned result is a numeric vector of length $1$. The attribute named "gradient" of the returned value (if provided in accordance with the value of hasGrad) must then be a numeric vector with length equal to the number of covariance parameters. It must contain the derivative of the kernel value $K(\mathbf{x}_1,\,\mathbf{x}_2;\,\boldsymbol{\theta})$ with respect to the parameter vector $\boldsymbol{\theta}$.

• When acceptMatrix is TRUE, the formals x1 and x2 are matrices with $d$ columns and with $n_1$ and $n_2$ rows. The result is then a covariance matrix with $n_1$ rows and $n_2$ columns. The gradient attribute (if provided in accordance with the value of hasGrad) must be a list with length equal to the number of covariance parameters. The list element $\ell$ must contain a numeric matrix with dimension $(n_1, n_2)$ which is the derivative of the covariance matrix w.r.t. the covariance parameter $\theta_\ell$.

Note

The kernel function must be symmetric with respect to its first two arguments, and it must be positive definite, which is not checked. If the function returns an object with a "gradient" attribute but hasGrad was set to FALSE, the gradient will not be used in optimization.

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

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 C
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;
     "
     
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 arguments
myCov <- 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 points
X <- matrix(c(0, 0), ncol = 2)
Xnew <- matrix(c(0.5, 1), ncol = 2)
colnames(X) <- colnames(Xnew) <- inputNames(myCov)
covMat(myCov, X)            ## same points
covMat(myCov, X, Xnew)      ## two different points

# computing covariances between sets of given locations
X <- 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 matrix
covMat(myCov, X, Xnew)   ## covariances between design and new data

Description

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

Objects from the Class

Objects can be created by calling new("covMan", ...) or by using the covMan function.

Slots

kernel:

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

hasGrad:

logical indicating whether kernel returns the gradient (w.r.t. the vector of parameters) as "gradient" attribute of the result.

acceptMatrix:

logical indicating whether kernel admits matrix arguments. Default is FALSE.

label:

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

d:

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

inputNames:

object of class "character", vector of input names. Length d.

parLower:

,

parUpper:

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

par:

object of class "numeric", numeric vector of parameter values.

parN:

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

kernParNames:

object of class "character", name of the kernel parameters.

Methods

coef<-

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

coefLower<-

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

coefLower

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

coef

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

coefUpper<-

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

coefUpper

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

covMat

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

scores

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

show

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

Note

While the coef<- replacement method is typically intended for internal use during likelihood maximization, the coefLower<- and coefUpper<- replacement methods can be used when some information is available on the possible values of the parameters.

Author(s)

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

See Also

The covMan function providing a creator.

Examples

showClass("covMan")

Description

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

Usage

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

Arguments

Value

A rectangular matrix with nrow(X) rows and nrow(Xnew) columns containing the covariances $K(\mathbf{x}_1, \mathbf{x}_2)$ for all the couples of sites $\mathbf{x}_1$ and $\mathbf{x}_2$.

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 the C code using eval.

Value

A $n_1 \times n_2$ matrix with general element $C_{ij} := K(\mathbf{x}_{1,i},\,\mathbf{x}_{2,j};\,\boldsymbol{\theta})$ where $K(\mathbf{x}_1,\,\mathbf{x}_2;\,\boldsymbol{\theta})$ is the covariance kernel function.

Note

The value of the parameter $\boldsymbol{\theta}$ can be extracted from the object with the coef 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 colnames
colnames(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))

Description

Creator function for the class covOrd-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.

Let $u$ be an ordered factor with levels $u_1, \dots, u_L$. Let $k_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. Then $k$ is defined by:

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

where $z_1, \dots, z_L$ form a regular sequence from $0$ to $1$ (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 unnormalized F, we set $\theta$ to 1, in order to avoid overparameterization.

Value

An object of class covOrd-class, inheriting from covQual-class.

See Also

covOrd-class, warpFun

Examples

u <- ordered(1:6, labels = letters[1:6])

myCov <- covOrd(ordered = u, cov = "homo", intAsChar = FALSE)
myCov
coef(myCov) <- c(mean = 0.5, sd = 1, theta = 3, sigma2 = 2)
myCov

checkX(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) 
myCov2

plot(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))
myCov3

plot(myCov3, type = "cor", method = "ellipse")
plot(myCov3, type = "warp", col = "blue", lwd = 2)


str(warpPower)  # details on the list describing the Power cdf
str(warpNorm)   # details on the list describing the Normal cdf

Description

Covariance kernel for qualitative ordered inputs obtained by warping.

Let $u$ be an ordered factor with levels $u_1, \dots, u_L$. Let $k_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. Then $k$ is defined by:

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

where $z_1, \dots, z_L$ form a regular sequence from $0$ to $1$ (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 form new("covOrd", ...).

Slots

covLevels:

Same as for covQual-class.

covLevMat:

Same as for covQual-class.

hasGrad:

Same as for covQual-class.

acceptLowerSQRT:

Same as for covQual-class.

label:

Same as for covQual-class.

d:

Same as for covQual-class. Here equal to 1.

inputNames:

Same as for covQual-class.

nlevels:

Same as for covQual-class.

levels:

Same as for covQual-class.

parLower:

Same as for covQual-class.

parUpper:

Same as for covQual-class.

par:

Same as for covQual-class.

parN:

Same as for covQual-class.

kernParNames:

Same as for covQual-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 value 0L corresponds to a correlation kernel while 1L is for a covariance kernel.

parNk1:

Object of class "integer". Number of parameters of k1Fun1. Equal to 0 at this stage.

parNwarp:

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

k1ParNames:

Object of class "character". Parameter names of k1Fun1.

warpParNames:

Object of class "character". Parameter names of warpFun.

warpKnots:

Object of class "numeric". Parameters of warpFun.

ordered:

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

intAsChar:

Object of class "logical". If TRUE (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 that the inputs exist with suitable column names and suitable factor content. The levels should match the prescribed levels. Returns a matrix with the input columns in the order prescribed by object.

signature(object = "covOrd", X = "matrix"): check that the inputs exist with suitable column names and suitable numeric content for coercion into a factor with the prescribed levels. Returns a data frame with the input columns in the order prescribed by object.

coef<-

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

coefLower<-

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

coefLower

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

coef

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

coefUpper<-

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

coefUpper

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

covMat

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

npar

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

scores

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

simulate

signature(object = "covOrd"): simulate nsim paths from a Gaussian Process having the covariance structure. The paths are indexed by the finite set of levels of factor inputs, and they are returned as columns of a matrix.

varVec

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

Note

This class is to be regarded as experimental. The slot names or list may be changed in the future. The methods npar, inputNames or `inputNames<-` should provide a more robust access to some slot values.

See Also

See covMan for a comparable structure dedicated to kernels with continuous inputs.

Examples

showClass("covOrd")

Description

Covariance kernel for qualitative inputs.

Objects from the Class

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

Slots

covLevels:

Object of class "function". This function has arguments 'par' and optional arguments lowerSQRT and compGrad. 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 function covLevels (former slot) with lowerSQRT = FALSE and gradient = FALSE.

hasGrad:

Object of class "logical". When TRUE, the covariance matrix returned by the function in slot covLevels must compute the gradients. The returned covariance matrix must have a "gradient" attribute; this must be an array with dimension c(m, m, np) where m stands for the number of levels and $np$ is the number of parameters.

acceptLowerSQRT:

Object of class "logical". When TRUE, the function in slot covLevels must have a formal lowerSQRT which can receive a logical value. When the value is TRUE the Cholesky (lower) root of the covariance is returned instead of the covariance.

label:

Object of class "character". A description of the kernel which 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 frame when the kernel will be evaluated.

nlevels:

Object of class "integer". A vector with length d giving the number of levels for each of the d inputs.

levels:

Object of class "list". A list of length d containing the d character vectors of levels for the d (qualitative) inputs.

parLower:

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

parUpper:

Object of class "numeric". Vector of parN upper values for the parameters of the structure. The value Inf can be used when needed.

par:

Object of class "numeric". Vector of parN current values for the structure.

parN:

Object of class "integer". Number of parameters for the structure, as returned by the npar method.

kernParNames:

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

ordered:

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

intAsChar:

Object of class "logical" indicating how to cope with an integer input. When intAsChar is TRUE the input is coerced into a character; the values taken by this character vector should then match the levels in the covQual object as given by levels(object)[[1]]. If instead intAsChar is FALSE, the integer values are assumed to correspond to the levels of the covQual object in the same order.

Methods

checkX

signature(object = "covQual", X = "data.frame"): check that the inputs exist with suitable column names and suitable factor content. The levels should match the prescribed levels. Returns a matrix with the input columns in the order prescribed by object.

signature(object = "covQual", X = "matrix"): check that the inputs exist with suitable column names and suitable numeric content for coercion into a factor with the prescribed levels. Returns a data frame with the input columns in the order prescribed by object.

coef<-

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

coefLower<-

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

coefLower

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

coef

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

coefUpper<-

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

coefUpper

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

covMat

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

npar

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

plot

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

scores

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

simulate

signature(object = "covQual"): simulate nsim paths from a Gaussian Process having the covariance structure. The paths are indexed by the finite set of levels of factor inputs, and they are returned as columns of a matrix.

varVec

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

Note

This class is to be regarded as experimental. The slot names or list may be changed in the future. The methods npar, inputNames or `inputNames<-` should provide a more robust access to some slot values.

See Also

See covMan for a comparable structure dedicated to kernels with continuous inputs.

Examples

showClass("covQual")

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

When covBet and covWith are zero, the resulting matrix is not a correlation matrix, due to the mode of construction. The "between" covariance matrix is a correlation but diagonal blocks are added to the extended matrix obtained by re-sizing the "between" covariance into a $n \times n$ matrix.

Note

For now the replacement method such as 'coef<-' are inherited from the class covQuall. Consequently when these methods are used they do not update the covariance structure in the between slot nor those in the within (list) slot.

This covariance kernel involves two categorical (i.e. factor) inputs, but these are nested. It could be aliased in the future as q1Nested or q2Nested.

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 invocated
nest2
plot(nest2, type = "cor")

## check that the covariance matrices match for nest1 and nest2
max(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"))

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 form new("covQualNested", ...).

Slots

covLevels:

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

covLevMat:

Object of class "matrix". The matrix returned by the function in slot covLevels. Since this matrix is often needed, it can be stored rather than recomputed.

hasGrad:

Object of class "logical". If TRUE, the analytical gradient can be computed.

acceptLowerSQRT:

Object of class "logical". If TRUE, the lower square root of the matrix can be returned

label:

Object of class "character". A label to describe the kernel.

d:

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

inputNames:

Object of class "character" Names of the inputs.

nlevels:

Object of class "integer" with length d give the number of levels for the factors.

levels:

Object of class "list" with length d. Gives the levels 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 final level.

groupLevels:

Object of class "character". Vector of labels for the groups.

between:

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

within:

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

parNCum:

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

contrasts:

Object of class "function". A contrast function like contr.helmod. This function must return a contrast matrix with 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")

Description

Creator for the class "covRadial", which describes radial kernels.

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 the $d$-dimensional Euclidean space takes the form

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

where $k_1(r)$ is a suitable correlation kernel for a one-dimensional input, and $r$ 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 on $d + 1$ parameters: the ranges $\theta_\ell >0$ and the variance $\sigma^2$.

An isotropic form uses the same range $\theta$ for all inputs, i.e. sets $\theta_\ell = \theta$ for all $\ell$. This is obtained by using iso = TRUE.

A correlation version uses $\sigma^2 = 1$. This is obtained by using cov = "corr".

Finally, the correlation kernel $k_1(r)$ can depend on a "shape" parameter, e.g. have the form $k_1(r;\,\alpha)$. The extra shape parameter $\alpha$ will be considered then as a parameter of the resulting radial kernel, making it possible to estimate it by ML along with the range(s) and the variance.

Value

An object with class "covRadial".

Note

When k1Fun1 has more than one formal argument, its arguments with position > 1 are assumed to be "shape" parameters of the model. Examples are functions with formals function(x, shape = 1.0) or function(x, alpha = 2.0, beta = 3.0), corresponding to vector of parameter names c("shape") and c("alpha", "beta"). Using more than one shape parameter has not been tested yet.

Remind that using a one-dimensional correlation kernel $k_1(r)$ here does not warrant that a positive semi-definite kernel will result for any dimension $d$. This question relates to Schoenberg's theorem and the concept of completely monotone functions.

References

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

See Also

k1Fun1Exp, k1Fun1Matern3_2, k1Fun1Matern5_2 or k1Fun1Gauss for examples of functions that can be used as values for the k1Fun1 formal.

Examples

set.seed(123)
d <- 2; ng <- 20
xg <- 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 <- 2
myCovRadial <- 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))

Description

Class of radial covariance kernels.

Objects from the Class

Objects can be created by calls of the form covRadial(...) of new("covRadial", ...).

Slots

k1Fun1:

Object of class "function" A function of a scalar numeric variable. Note that using a one-dimensional kernel here does not warrant that a positive semi-definite kernel results for any dimension $d$.

k1Fun1Char:

Object of class "character" describing the function in the slot k1Fun1.

hasGrad:

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

cov:

Object of class "integer". The value 0L corresponds to a correlation kernel while 1L is for a covariance kernel.

iso:

Object of class "integer". The value 1L corresponds to an isotropic covariance, with all the inputs sharing the same range value.

label:

Object of class "character". Short description of the object.

d:

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

inputNames:

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

parLower:

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

parUpper:

Object of class "numeric". Numeric values for the upper bounds on the parameters. Can be Inf.

par:

Object of class "numeric". Numeric values for the parameters. Can be NA.

parN1:

Object of class "integer". Number of parameters of the function kern1Fun (such as a shape).

parN:

Object of class "integer". Number of parameters for the object. The include: direct parameters in the function kern1Fun, ranges, and variance.

kern1ParNames:

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

kernParNames:

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

Extends

Class "covAll", directly.

Methods

coef<-

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

coefLower<-

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

coefLower

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

coef

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

coefUpper<-

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

coefUpper

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

covMat

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

npar

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

scores

signature(object = "covRadial"): Compute the scores i.e. the derivatives w.r.t. the parameters of the contribution of the covariance in the log-likelihood of a gp.

show

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

varVec

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