corr_matrix function

Power Exponential or Matern Correlation Matrix

Computes the power exponential or Matern correlation matrix for a set of n design points in d-dimensional input region and a vector of d correlation hyper-parameters (beta).

corr_matrix(X, beta, corr = list(type = "exponential", power = 1.95))

Arguments

• X: the (n x d) design matrix
• beta: a (d x 1) vector of correlation hyper-parameters in $(-\infty, \infty)$
• corr: a list that specifies the type of correlation function along with the smoothness parameter. The default corresponds to power exponential correlation with smoothness parameter "power=1.95". One can specify a different power (between 1.0 and 2.0) for the power exponential, or use the Matern correlation function, specified as corr=list(type = "matern", nu=(2*k+1)/2), where c("$k \\in\n$", "$\\{0,1,2,...\\}$")

Returns

The (n x n) correlation matrix, R, for the design matrix (X) and the hyper-parameters (beta).

Details

The power exponential correlation function is given by

c("$\n$", "R_{ij} = \\prod exp(-10^\\beta*|x_{i}-x_{j}|^{power})").

The Matern correlation function given by Santner, Williams, and Notz (2003) is

c("$R_{ij} = \\prod\n$", "\\frac{1}{\\Gamma(\\nu)2^{\\nu-1}}(2\\sqrt{\\nu}|x_{ik} - x_{jk}|10^{\\beta_k})^\\nu\n", "$\\kappa_{\\nu}(2\\sqrt{\\nu}|x_{ik} -\n$", "$x_{jk}|10^{\\beta_k})$")c("$R_{ij} = \\prod\n$", "\\frac{1}{\\Gamma(\\nu)2^{\\nu-1}}(2\\sqrt{\\nu}|x_{ik} - x_{jk}|10^{\\beta_k})^\\nu\n", "$\\kappa_{\\nu}(2\\sqrt{\\nu}|x_{ik} -\n$", "$x_{jk}|10^{\\beta_k})$")c("$R_{ij} = \\prod\n$", "\\frac{1}{\\Gamma(\\nu)2^{\\nu-1}}(2\\sqrt{\\nu}|x_{ik} - x_{jk}|10^{\\beta_k})^\\nu\n", "$\\kappa_{\\nu}(2\\sqrt{\\nu}|x_{ik} - x_{jk}|10^{\\beta_k})$"),

where $\kappa_{\nu}$ is the modified Bessel function of order $\nu$.

Note

Both Matern and power exponential correlation functions use the new $\beta$ parametrization of hyper-parameters given by c("$\\theta_k =\n$", "$10^{\\beta_k}$") for easier likelihood optimization. That is, beta is a log scale parameter (see MacDonald et al. (2015)).

Examples

## 1D Example - 1
n = 5
d = 1
set.seed(3)
library(lhs)
x = maximinLHS(n,d)
beta =  rnorm(1)
corr_matrix(x,beta)

## 1D Example - 2
beta = rnorm(1)
corr_matrix(x,beta,corr = list(type = "matern"))

## 2D example - 1
n = 10
d = 2
set.seed(2)
library(lhs)
x = maximinLHS(n,d)
beta = rnorm(2)
corr_matrix(x, beta,
    corr = list(type = "exponential", power = 2))

## 2D example - 2
beta = rnorm(2)
R = corr_matrix(x,beta,corr = list(type = "matern", nu = 5/2))
print(R)

References

MacDonald, K.B., Ranjan, P. and Chipman, H. (2015). GPfit: An R Package for Fitting a Gaussian Process Model to Deterministic Simulator Outputs. Journal of Statistical Software, 64(12), 1-23. https://www.jstatsoft.org/v64/i12/

Ranjan, P., Haynes, R., and Karsten, R. (2011). A Computationally Stable Approach to Gaussian Process Interpolation of Deterministic Computer Simulation Data, Technometrics, 53(4), 366 - 378.

Santner, T.J., Williams, B., and Notz, W. (2003), The design and analysis of computer experiments, Springer Verlag, New York.

Author(s)

Blake MacDonald, Hugh Chipman, Pritam Ranjan