corLevSymm: Correlation Matrix for a General Symmetric Correlation...
In kergp: Gaussian Process Laboratory

Description Usage Arguments Details Value Note References See Also Examples

Compute the correlation matrix for a general symmetric correlation structure.

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

`par`	A numeric vector with length `npCor + npVar` where `npCor = nlevels *` `(nlevels - 1) / 2` is the number of correlation parameters, and `npVar` is the number of variance parameters, which depends on the value of `cov`. The value of `npVar` is `0`, `1` or `nlevels` corresponding to the values of `cov`: `0`, `1` and `2`. The correlation parameters are assumed to be located at the head of `par` i.e. at indices `1` to `npCor`. The variance parameter(s) are assumed to be at the tail, i.e. at indices `npCor + 1` to `npCor + npVar`.
`nlevels`	Number of levels.
`levels`	Character representing the levels.
`lowerSQRT`	Logical. When `TRUE` the (lower) Cholesky root L of the correlation or covariance matrix C is returned instead of the correlation matrix.
`compGrad`	Logical. Should the gradient be computed? This is only possible for the C implementation.
`cov`	Integer `0`, `1` or `2`. If `cov` is `0`, the matrix is a correlation matrix (or its Cholesky root). If `cov` is `1` or `2`, the matrix is a covariance (or its Cholesky root) with constant variance vector for `code = 1` and with arbitrary variance for `code = 2`. The variance parameters `par` are located at the tail of the `par` vector, so at locations `npCor + 1` to `npCor + nlevels` when `code = 2` where `npCor` is the number of correlation parameters, i.e. `nlevels * (nlevels - 1) / 2`.
`impl`	A character telling which of the C and R implementations should be chosen.

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 * (n - 1) / 2 parameters theta[i, j] where the indices i and j are such that 1 <= j < i <= n. The parameters theta[i, j] are angles and are to be taken to be in [0, pi) for a one-to-one parameterisation.

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

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[i, j] 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.

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.

The corSymm correlation structure in the nlme package.

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")
}