blockCor: Correlation Matrix without or with Derivatives

Description

Calculation of a correlation matrix with or without derivatives according to the specification of a correlation structure.

Usage

blockCor(x, ...)
# S3 method for default
blockCor(x, theta, covtype = c("matern5_2", "matern3_2", "gaussian"), 
	derivatives = FALSE, ...)
# S3 method for gekm
blockCor(x, ...)

Value

blockCor returns a list with the following components:

K: The correlation matrix without derivatives.
R: If derivatives = TRUE, the correlation matrix with first partial derivatives, otherwise NULL.
S: If derivatives = TRUE, the correlation matrix with second partial derivatives, otherwise NULL.

The components of the list can be combined in the following form to constructed the complete correlation matrix with derivatives: cbind(rbind(K, R), rbind(t(R), S)).

Arguments

x: a numeric matrix or an object of class gekm.
theta: numeric vector of length d for the hyperparameters.
covtype: character specifying the correlation function to be used. Must be one of "matern5_2", "matern3_2" or "gaussian". See ‘Details’.
derivatives: logical, if TRUE the first and second partial derivatives of the correlation matrix are calculated, otherwise not.
...: further arguments passed to or from other methods.

Author

Carmen van Meegen

Details

The correlation matrix with derivatives is defined as a block matrix of the form $$\begin{pmatrix} K & R^{\rm T} \\ R & S \end{pmatrix}.$$

Three correlation functions are currently implemented in blockCor:

Matérn 5/2 kernel with covtype = "matern5_2": $$\phi(x, x'; \theta) = \prod_{k = 1}^d \left(1 + \frac{\sqrt{5}|x_k - x_k'|}{\theta_k} + \frac{5 |x_k - x_k'|^2}{3 \theta_k^2}\right) \exp\left( -\frac{\sqrt{5}|x_k - x_k'|}{\theta_k}\right)$$
Matérn 3/2 kernel with covtype = "matern3_2": $$\phi(x, x'; \theta) = \prod_{k = 1}^d \left(1 + \frac{\sqrt{3} |x_k - x_k'|}{\theta_k} \right) \exp\left( -\frac{\sqrt{3} |x_k - x_k'|}{\theta_k}\right)$$
Gaussian kernel with covtype = "gaussian": $$\phi(x, x'; \theta) = \prod_{k = 1}^d \exp\left( -\frac{(x_k - x_k')^2}{2\theta_k^2}\right)$$

References

Koehler, J. and Owen, A. (1996). Computer Experiments. In Ghosh, S. and Rao, C. (eds.), Design and Analysis of Experiments, volume 13 of Handbook of Statistics, pp. 261–308. Elsevier Science. tools:::Rd_expr_doi("10.1016/S0169-7161(96)13011-X").

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. The MIT Press. https://gaussianprocess.org/gpml/.

Stein, M. L. (1999). Interpolation of Spatial Data: Some Theory for Kriging. Springer Series in Statistics. Springer-Verlag. tools:::Rd_expr_doi("10.1007/978-1-4612-1494-6").

Examples

Run this code

# Some examples of correlation matrices:

x <- matrix(1:4, ncol = 1)
blockCor(x, theta = 1)
blockCor(x, theta = 1, covtype = "gaussian")
blockCor(x, theta = 1, covtype = "gaussian", derivatives = TRUE)
blockCor(x, theta = 1, covtype = "matern3_2", derivatives = TRUE)

res <- blockCor(x, theta = 2, covtype = "gaussian", derivatives = TRUE)
cbind(rbind(res$K, res$R), rbind(t(res$R), res$S))


# Illustration of correlation functions and their derivatives:

x <- seq(0, 1, length = 100)
X <- matrix(x, ncol = 1)
gaussian <- blockCor(X, theta = 0.25, covtype = "gaussian", derivatives = TRUE)
matern5_2 <- blockCor(X, theta = 0.25, covtype = "matern5_2", derivatives = TRUE)
matern3_2 <- blockCor(X, theta = 0.25, covtype = "matern3_2", derivatives = TRUE)

# Correlation functions and first partial derivatives:
index <- c(10, 20, 40, 80)
par(mar = c(5.1, 5.1, 4.1, 2.1))

# Matern 3/2
plot(x, matern3_2$K[1, ], type = "l", xlab = expression(group("|", x - x*minute, "|")),
	ylab = expression(phi(x, x*minute, theta == 0.25)), lwd = 2)
tangents(x[index], matern3_2$K[index, 1], matern3_2$R[index, 1], 
	length = 0.15, lwd = 2, col = 2)
points(x[index], matern3_2$K[index, 1], pch = 16)

# Matern 5/2
lines(x, matern5_2$K[1, ], lwd = 2, col = 3)
tangents(x[index], matern5_2$K[index, 1], matern5_2$R[index, 1], 
	length = 0.15, lwd = 2, col = 2)
points(x[index], matern5_2$K[index, 1], pch = 16)

# Gaussian
lines(x, gaussian$K[1, ], lwd = 2, col = 4)
tangents(x[index], gaussian$K[index, 1], gaussian$R[index, 1], 
	length = 0.15, lwd = 2, col = 2)
points(x[index], gaussian$K[index, 1], pch = 16)

legend("topright", lty = 1, lwd = 2, col = c(1, 3, 4), bty = "n",
	legend = c("Matern 3/2", "Matern 5/2", "Gaussian"))


# First and second partial derivatives of correlation functions:
index <- c(5, 10, 20, 40, 80)
par(mar = c(5.1, 6.1, 4.1, 2.1))

# Gaussian
plot(x, matern3_2$R[1, ], type = "l", xlab = expression(group("|", x - x*minute, "|")),
	ylab = expression(frac(partialdiff * phi(x, x*minute, theta == 0.25), 
		partialdiff * x * minute)), lwd = 2)
tangents(x[index], matern3_2$R[1, index], matern3_2$S[index, 1], 
	length = 0.15, lwd = 2, col = 2)
points(x[index], matern3_2$R[1, index], pch = 16)

# Matern 5/2
lines(x, matern5_2$R[1, ], lwd = 2, col = 3)
tangents(x[index], matern5_2$R[1, index], matern5_2$S[index, 1], 
	length = 0.15, lwd = 2, col = 2)
points(x[index], matern5_2$R[1, index], pch = 16)

# Matern 3/2
lines(x, gaussian$R[1, ], lwd = 2, col = 4)
tangents(x[index], gaussian$R[1, index], gaussian$S[index, 1], 
	length = 0.15, lwd = 2, col = 2)
points(x[index], gaussian$R[1, index], pch = 16)

legend("topright", lty = 1, lwd = 2, col = c(1, 3, 4), bty = "n",
	legend = c("Matern 3/2", "Matern 5/2", "Gaussian"))

Run the code above in your browser using DataLab