# R/corGaussian.R In OakleyJ/MUCM: Gaussian process emulator methods based on the MUCM toolkit

#' @title Correlation functions
#'
#' @description Functions to calculate the correlation between two points.
#' The correlation functions give the correlation between the simulator output at any given vector of input values and
#' the output at any other given input vector. Within this package, when these two vectors are the same then the correlation
#' function \code{C} at input vector \code{x} is \code{C(x,x) = 1}
#' In principal, the correlation function can take a very wide variety of forms. Only a few are included in this package, however,
#' a user is able to define his own function. A user defined function must follow the same signature as those included in this package.

#' @param inputs A data frame, matrix or vector containing the input values of the data where each row gives the observations.
#' @param inputs2 An additional data frame, matrix or vector containing the input values of the data where each row gives the observations.
#' @param phi Estimate of the roughness parameter (a vector).

#' @return Returns the correlation matrix \eqn{A} calculated by the specific function.

#' @details \code{corGaussian} uses the Gaussian correlation function defined as
#' \deqn{\prod _{i=1} ^ {n} \exp [- (\frac{(x_i - x_i')}{\delta _i}) ^2]}{\prod _{i=1}^{n} exp(-((x_i - x_i') / \delta_i)^2)}
#'

#' @note These functions use \code{\link{rdist}} from package \pkg{fields} for speed.

#' @importFrom fields rdist

#' @author Sajni Malde and Jeremy Oakley
#' @export
corGaussian <- function(inputs, inputs2, phi) {

if (missing(inputs2) || is.null(inputs2))
return(corGaussianSquare(inputs, phi))

delta <- exp(phi)
exp(-(rdist(inputs / rep(delta, each = nrow(inputs)), inputs2 / rep(delta, each = nrow(inputs2))) ^ 2))
}

corGaussianSquare <- function(inputs, phi) {
delta <- exp(phi)
exp(-(rdist(inputs / rep(delta, each = nrow(inputs))) ^ 2))
}

############ ----------------   MATERN   ------------------- ################

#' @rdname corGaussian
#' @details \code{corMatern2.5} uses the Matern correlation function with \eqn{\nu = 2.5} defined as
#' \deqn{(1 + 5^{0.5}r + \frac{5}{3}r^2)\exp(-5^{0.5}r)}{(1 + (5^0.5)*r + (5*r^2)/3) * exp(-(5^0.5)*r)} where \eqn{r} is distance between inputs \eqn{x_i} and \eqn{x_j}, scaled by
#' delta: \deqn{[\frac{(x_{i1}-x_{j1})^2}{\delta _1^2} +...+ \frac{(x_{id}-x_{jd})^2}{\delta_d^2}]^{0.5}}{[(x_{i1}-x_{j1})^2 / \delta_1^2 +...+ (x_{id}-x_{jd})^2 / \delta_d^2 ]^{0.5}}
#' @export
corMatern2.5 <- function(inputs, inputs2, phi) {

if (missing(inputs2) || is.null(inputs2))
return(corMaternSquare2.5(inputs, phi))

delta <- exp(phi)

# calculate distance matrix
distances <- rdist(inputs / rep(delta, each = nrow(inputs)), inputs2/rep(delta, each = nrow(inputs2)))
squared.distances <- distances ^ 2

# return correlation matrix
(1 + sqrt(5) * distances + (5 / 3) * squared.distances) * exp(-sqrt(5) * distances)
}

corMaternSquare2.5 <- function(inputs, phi) {
n <- nrow(inputs)
delta <- exp(phi)

# calculate distance matrix
distances <- rdist(inputs / rep(delta, each = n))
squared.distances <- distances ^ 2

# return correlation matrix
(1 + sqrt(5) * distances + (5 / 3) * squared.distances) * exp(-sqrt(5) * distances)
}

OakleyJ/MUCM documentation built on May 7, 2019, 9:01 p.m.