corLevCompSymm: Correlation Matrix for the Compound Symmetry Structure
In kergp: Gaussian Process Laboratory

Description Usage Arguments Value Note Author(s) Examples

Compute the correlation matrix for a the compound symmetry structure.

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

`par`	Numeric vector of length `1` if `cov` is `TRUE` or with length `2` else. The first element is the correlation coefficient and the second one (when it exists) is the variance.
`nlevels`	Number of levels.
`levels`	Character representing the levels.
`lowerSQRT`	Logical. When `TRUE` the (lower) Cholesky root L of the correlation matrix C is returned instead of the correlation matrix.
`compGrad`	Logical. Should the gradient be computed?
`cov`	Logical. If `TRUE` the matrix is a covariance matrix (or its Cholesky root) rather than a correlation matrix and the last element in `par` is the variance.
`impl`	A character telling which of the C and R implementations should be chosen.

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

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

Yves Deville

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