kergp (version 0.5.7)

covMan: Creator Function for covMan Objects

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

kernel

A (semi-)positive definite function. This must be an object of class "function" with formal arguments named "x1", "x2" and "par". The first two formal arguments are locations vectors or matrices. The third formal is for the vector \(\boldsymbol{\theta}\) of all covariance parameters. An analytical gradient can be computed and returned as an attribute of the result with name "gradient". See Details.

hasGrad

Logical indicating whether the kernel function returns the gradient w.r.t. the vector of parameters as a "gradient" attribute of the result. See Details

acceptMatrix

Logical indicating whether kernel admits matrices as arguments. Default is FALSE. See Examples below.

inputs

Character vector giving the names of the inputs used as arguments of kernel. Optional if d is given.

d

Integer specifying the spatial dimension (equal to the number of inputs). Optional if inputs is given.

parNames

Vector of character strings containing the parameter names.

par, parLower, parUpper

Optional numeric vectors containing the parameter values, lower bounds and upper bounds.

label

Optional character string describing the kernel.

...

Not used at this stage.

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\).

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! 

if (FALSE) {
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
}

## 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