R/CGGP_corr_fs.R In CollinErickson/CGGP: Composite Grid Gaussian Processes

Documented in CGGP_internal_CorrMatCauchyCGGP_internal_CorrMatCauchySQCGGP_internal_CorrMatCauchySQTCGGP_internal_CorrMatGaussianCGGP_internal_CorrMatMatern32CGGP_internal_CorrMatMatern52CGGP_internal_CorrMatPowerExpCGGP_internal_CorrMatWendland0CGGP_internal_CorrMatWendland1CGGP_internal_CorrMatWendland2

#' Cauchy correlation function
#'
#' Calculate correlation matrix for two sets of points in one dimension.
#' Note that this is not the correlation between two vectors.
#'
#' @param x1 Vector of coordinates from same dimension
#' @param x2 Vector of coordinates from same dimension
# ' @param LS Log of parameter that controls lengthscale
# ' @param FD Logit of 0.5*parameter  that controls the fractal dimension
# ' @param HE Log of parameter that controls the hurst effect
#' @param theta Correlation parameters:
#' \itemize{
#'     \item LS Log of parameter that controls lengthscale
#'     \item FD Logit of 0.5*parameter  that controls the fractal dimension
#'     \item HE Log of parameter that controls the hurst effect
#' }
#' @param return_dCdtheta Should dCdtheta be returned?
#' @param return_numpara Should it just return the number of parameters?
#' @param returnlogs Should log of correlation be returned?
#'
#' @return Matrix of correlation values between x1 and x2
#' @export
#' @family correlation functions
#'
#' @examples
#' CGGP_internal_CorrMatCauchy(c(0,.2,.4),c(.1,.3,.5), theta=c(-1,.9,.1))
CGGP_internal_CorrMatCauchy <- function(x1, x2, theta, return_dCdtheta=FALSE,
return_numpara=FALSE,
returnlogs=FALSE) {
if(return_numpara){
return(3)
}else{
if (length(theta) != 3) {stop("CorrMatCauchy theta should be length 3")}
diffmat =abs(outer(x1,x2,'-'));

expLS = exp(3*(theta[1]))
expHE = exp(3*(theta[2]))
h = diffmat/expLS
alpha = 2*exp(3*theta[3]+2)/(1+exp(3*theta[3]+2))
halpha = h^alpha
pow = -expHE/alpha

if (!returnlogs) {
C = (1+halpha)^pow
} else {
C = pow * log(1+halpha)
}
if(return_dCdtheta){
if (!returnlogs) {
dCdtheta = cbind(3*expHE*((1+halpha)^(pow-1))*(halpha),
3*C*pow*log(1+halpha),
(C*(expHE*log(1+halpha)/alpha^2 -
expHE*halpha*log(h)/alpha/(1+halpha))) *
6*exp(3*theta[3]+2)/(1+exp(3*theta[3]+2))^2)
} else {
dCdtheta = cbind(3*expHE*halpha/(1+halpha),
3*pow*log(1+halpha),
((expHE*log(1+halpha)/alpha^2 -
expHE*halpha*log(h)/alpha/(1+halpha))) *
6*exp(3*theta[3]+2)/(1+exp(3*theta[3]+2))^2)
}
dCdtheta[is.na(dCdtheta)] = 0
out <- list(C=C,dCdtheta=dCdtheta)
return(out)
}else{
return(C)
}
}
}

#' CauchySQT correlation function
#'
#' Calculate correlation matrix for two sets of points in one dimension.
#' Note that this is not the correlation between two vectors.
#'
#' @inheritParams CGGP_internal_CorrMatCauchy
#'
#' @return Matrix of correlation values between x1 and x2
#' @export
#' @family correlation functions
#'
#' @examples
#' CGGP_internal_CorrMatCauchySQT(c(0,.2,.4),c(.1,.3,.5), theta=c(-.1,.3,-.7))
CGGP_internal_CorrMatCauchySQT <- function(x1, x2,theta, return_dCdtheta = FALSE,
return_numpara=FALSE,
returnlogs=FALSE) {
if(return_numpara){
return(3);
} else{
if (length(theta) != 3) {stop("CorrMatCauchySQT theta should be length 3")}
expTILT = exp((theta[3]))
expLS = exp(3*(theta[1]))
x1t = (x1+10^(-2))^expTILT
x2t = (x2+10^(-2))^expTILT
x1ts = x1t/expLS
x2ts = x2t/expLS

diffmat =abs(outer(x1ts,x2ts,'-'));
expHE = exp(3*(theta[2]))
h = diffmat
alpha = 2*exp(5)/(1+exp(5))
halpha = h^alpha
pow = -expHE/alpha

if (!returnlogs) {
C = (1+halpha)^pow
} else {
C = pow * log(1+halpha)
}

if(return_dCdtheta){
Q = ((1+halpha)^(pow-1))
gt1 = x1t*log(x1+10^(-2))
gt2 = x2t*log(x2+10^(-2))
lh =outer(gt1,gt2,'-')
hnabs = outer(x1ts,x2ts,'-')
LO = alpha*expTILT*(pow/expLS)*(abs(h)^(alpha-1)*lh*sign(hnabs))
if (!returnlogs) {
dCdtheta = cbind(3*expHE*((1+halpha)^(pow-1))*(halpha),3*C*pow*log(1+halpha),LO*Q)
} else {
dCdtheta = cbind(3*expHE*halpha/(1+halpha), 3*pow*log(1+halpha), LO/(1+halpha))
}
out <- list(C=C,dCdtheta=dCdtheta)
return(out)
}else{
return(C)
}
}
}

#' CauchySQ correlation function
#'
#' Calculate correlation matrix for two sets of points in one dimension
#' Note that this is not the correlation between two vectors.
#'
#' @inheritParams CGGP_internal_CorrMatCauchy
#'
#' @return Matrix of correlation values between x1 and x2
#' @export
#' @family correlation functions
#'
#' @examples
#' CGGP_internal_CorrMatCauchySQ(c(0,.2,.4),c(.1,.3,.5), theta=c(-.7,-.5))
CGGP_internal_CorrMatCauchySQ <- function(x1, x2,theta, return_dCdtheta = FALSE,
return_numpara =FALSE,returnlogs = FALSE) {
if(return_numpara){
return(2);
}else{
if (length(theta) != 2) {stop("CorrMatCauchySQ theta should be length 2")}
diffmat =abs(outer(x1,x2,'-'));

expLS = exp(3*theta[1])
expHE = exp(3*theta[2])
h = diffmat/expLS
alpha = 2*exp(0+6)/(1+exp(0+6))
halpha = h^alpha
pow = -expHE/alpha

if(!returnlogs){
C = (1+halpha)^pow
}else{
C = pow*log(1+halpha)
}
if(return_dCdtheta){
if(!returnlogs){
dCdtheta = cbind(3*expHE*((1+halpha)^(pow-1))*(halpha),3*C*pow*log(1+halpha))
}else{
dCdtheta = cbind(3*expHE*halpha/(1+halpha),3*C)
}
dCdtheta[is.na(dCdtheta)] = 0
out <- list(C=C,dCdtheta=dCdtheta)
return(out)
}else{
return(C)
}
}
}

#' Gaussian correlation function
#'
#' Calculate correlation matrix for two sets of points in one dimension
#' Note that this is not the correlation between two vectors.
#'
#' WE HIGHLY ADVISE NOT USING THIS CORRELATION FUNCTION.
#' Try Power Exponential, CauchySQT, Cauchy, or Matern 3/2 instead.
#'
#' @inheritParams CGGP_internal_CorrMatCauchy
#'
#' @return Matrix of correlation values between x1 and x2
#' @export
#' @family correlation functions
#'
#' @examples
#' CGGP_internal_CorrMatGaussian(c(0,.2,.4),c(.1,.3,.5), theta=c(-.7))
CGGP_internal_CorrMatGaussian <- function(x1, x2,theta, return_dCdtheta = FALSE,
return_numpara=FALSE,
returnlogs=FALSE) {
if(return_numpara){
return(1);
}else{
if (length(theta) != 1) {stop("CorrMatGaussian theta should be length 1")}
diffmat =abs(outer(x1,x2,'-'))
diffmat2 <- diffmat^2

expLS = exp(3*theta[1])
h = diffmat2/expLS

# Gaussian corr is awful, always needs a nugget
nug <- 1e-10
if (!returnlogs) {
C = (1-nug)*exp(-h) + nug*(diffmat<10^(-4))
# C = exp(-h)
} else {
# C = -h
C = (1-nug)*exp(-h) + nug*(diffmat<10^(-4))
C <- log(C)
}
if(return_dCdtheta){
if (!returnlogs) {
dCdtheta <- 3*C*diffmat2 / expLS
} else {
dCdtheta <- 3*diffmat2 / expLS
}
dCdtheta[is.na(dCdtheta)] = 0
out <- list(C=C,dCdtheta=dCdtheta)
return(out)
}else{
return(C)
}
}
}

#' Matern 3/2 correlation function
#'
#' Calculate correlation matrix for two sets of points in one dimension.
#' Note that this is not the correlation between two vectors.
#'
#' @inheritParams CGGP_internal_CorrMatCauchy
#'
#' @return Matrix of correlation values between x1 and x2
#' @export
#' @family correlation functions
#'
#' @examples
#' CGGP_internal_CorrMatMatern32(c(0,.2,.4),c(.1,.3,.5), theta=c(-.7))
CGGP_internal_CorrMatMatern32 <- function(x1, x2,theta, return_dCdtheta=FALSE,
return_numpara=FALSE,
returnlogs=FALSE) {
if(return_numpara){
return(1);
}else{
if (length(theta) != 1) {stop("CorrMatMatern32 theta should be length 1")}
diffmat =abs(outer(x1,x2,'-'))

expLS = exp(3*theta[1])
h = diffmat/expLS

if (!returnlogs) {
# C = (1-10^(-10))*(1+sqrt(3)*h)*exp(-sqrt(3)*h) + 10^(-10)*(diffmat<10^(-4))
C = (1+sqrt(3)*h)*exp(-sqrt(3)*h)
} else {
C <- log(1+sqrt(3)*h) - sqrt(3)*h
}
if(return_dCdtheta){
if (!returnlogs) {
dCdtheta <- (sqrt(3)*diffmat*exp(-sqrt(3)*h) - sqrt(3)*C*diffmat) * (-3/expLS)
} else {
dCdtheta <- (sqrt(3)*diffmat/(1+sqrt(3)*h) - sqrt(3)*diffmat) * (-3/expLS)
}
dCdtheta[is.na(dCdtheta)] = 0
out <- list(C=C,dCdtheta=dCdtheta)
return(out)
}else{
return(C)
}
}
}

#' Matern 5/2 correlation function
#'
#' Calculate correlation matrix for two sets of points in one dimension.
#' Note that this is not the correlation between two vectors.
#'
#' @inheritParams CGGP_internal_CorrMatCauchy
#'
#' @return Matrix of correlation values between x1 and x2
#' @export
#' @family correlation functions
#'
#' @examples
#' CGGP_internal_CorrMatMatern52(c(0,.2,.4),c(.1,.3,.5), theta=c(-.7))
CGGP_internal_CorrMatMatern52 <- function(x1, x2,theta, return_dCdtheta=FALSE,
return_numpara=FALSE,
returnlogs=FALSE) {
if(return_numpara){
return(1);
}else{
if (length(theta) != 1) {stop("CorrMatMatern52 theta should be length 1")}
diffmat =abs(outer(x1,x2,'-'))
expLS = exp(3*theta[1])
h = diffmat/expLS
if (!returnlogs) {
# C = (1-10^(-10))*(1+sqrt(5)*h+5/3*h^2)*exp(-sqrt(5)*h) + 10^(-10)*(diffmat<10^(-4))
C = (1+sqrt(5)*h+5/3*h^2)*exp(-sqrt(5)*h)
} else {
C = log(1+sqrt(5)*h+5/3*h^2) - sqrt(5)*h
}
if(return_dCdtheta){
if (!returnlogs) {
dCdtheta <- ((sqrt(5)*diffmat+10/3*diffmat*h)*exp(-sqrt(5)*h) -
sqrt(5)*C*diffmat) * (-3/expLS)
} else {
dCdtheta <- ((sqrt(5)*diffmat+10/3*diffmat*h)/(1+sqrt(5)*h+5/3*h^2) -
sqrt(5)*diffmat) * (-3/expLS)
}
dCdtheta[is.na(dCdtheta)] = 0
out <- list(C=C,dCdtheta=dCdtheta)
return(out)
}else{
return(C)
}
}
}

#' Power exponential correlation function
#'
#' Calculate correlation matrix for two sets of points in one dimension.
#' Note that this is not the correlation between two vectors.
#'
#' @inheritParams CGGP_internal_CorrMatCauchy
#'
#' @return Matrix of correlation values between x1 and x2
# @rdname CGGP_internal_CorrMatCauchy
#' @export
#' @family correlation functions
#'
#' @examples
#' CGGP_internal_CorrMatPowerExp(c(0,.2,.4),c(.1,.3,.5), theta=c(-.7,.2))
CGGP_internal_CorrMatPowerExp <- function(x1, x2,theta,
return_dCdtheta = FALSE,
return_numpara=FALSE,
returnlogs=FALSE) {
if(return_numpara){
return(2);
}else{
if (length(theta) != 2) {stop("CorrMatPowerExp theta should be length 2")}
diffmat =abs(outer(x1,x2,'-'))
tmax <- 3
expLS = exp(tmax*theta[1])
minpower <- 1
maxpower <- 1.95
alpha <- minpower + (theta[2]+1)/2 * (maxpower - minpower)
h = diffmat/expLS
if (!returnlogs) {
# C = (1-nug)*exp(-(h)^alpha) + nug*(diffmat<10^(-4))
C = exp(-(h)^alpha)
} else {
C = -(h^alpha)
}
if(return_dCdtheta){
if (!returnlogs) {
dCdtheta <- cbind(tmax*alpha*C*diffmat^alpha/expLS^alpha,
-C*h^alpha*log(h)/2 * (maxpower - minpower))
} else {
dCdtheta <- cbind(tmax*alpha*diffmat^alpha/expLS^alpha,
-h^alpha*log(h)/2 * (maxpower - minpower))
}
dCdtheta[is.na(dCdtheta)] = 0
out <- list(C=C,dCdtheta=dCdtheta)
return(out)
}else{
return(C)
}
}
}

#' Wendland0 (Triangle) correlation function
#'
#' Calculate correlation matrix for two sets of points in one dimension.
#' Note that this is not the correlation between two vectors.
#'
#' @inheritParams CGGP_internal_CorrMatCauchy
#'
#' @return Matrix of correlation values between x1 and x2
# @rdname CGGP_internal_CorrMatCauchy
#' @export
#' @family correlation functions
#'
#' @examples
#' CGGP_internal_CorrMatWendland0(c(0,.2,.4),c(.1,.3,.5), theta=-.7)
CGGP_internal_CorrMatWendland0 <- function(x1, x2,theta,
return_dCdtheta = FALSE,
return_numpara=FALSE,
returnlogs=FALSE) {
if(return_numpara){
return(1)
}else{
if (length(theta) != 1) {stop("CorrMatWendland0 theta should be length 1")}
diffmat =abs(outer(x1,x2,'-'))
tmax <- 3

expLS = exp(tmax*theta[1])

wherecov = which(diffmat<expLS)
h = matrix(0,dim(diffmat)[1],dim(diffmat)[2])
h[wherecov] = 1-diffmat[wherecov] / expLS
if (!returnlogs) {
C = matrix(0,dim(h)[1],dim(h)[2])
C[wherecov] <- h[wherecov]
} else {
C = -Inf * matrix(1,dim(h)[1],dim(h)[2])
C[wherecov] = log(h[wherecov])
}
if(return_dCdtheta){
if (!returnlogs) {
dCdtheta = matrix(0,dim(diffmat)[1],dim(diffmat)[2])
dCdtheta[wherecov] <- tmax * (1-h[wherecov])
} else {
dCdtheta = matrix(0,dim(diffmat)[1],dim(diffmat)[2])
dCdtheta[wherecov] <- tmax * ((1-h[wherecov])/h[wherecov])
}
dCdtheta[is.na(dCdtheta)] = 0
out <- list(C=C,dCdtheta=dCdtheta)
return(out)
}else{
return(C)
}
}
}

#' Wendland1 1 correlation function
#'
#' Calculate correlation matrix for two sets of points in one dimension.
#' Note that this is not the correlation between two vectors.
#'
#' @inheritParams CGGP_internal_CorrMatCauchy
#'
#' @return Matrix of correlation values between x1 and x2
# @rdname CGGP_internal_CorrMatCauchy
#' @export
#' @family correlation functions
#'
#' @examples
#' CGGP_internal_CorrMatWendland1(c(0,.2,.4),c(.1,.3,.5), theta=-.7)
CGGP_internal_CorrMatWendland1 <- function(x1, x2,theta,
return_dCdtheta = FALSE,
return_numpara=FALSE,
returnlogs=FALSE) {
if(return_numpara){
return(1)
}else{
if (length(theta) != 1) {stop("CorrMatWendland1 theta should be length 1")}
diffmat =abs(outer(x1,x2,'-'))
tmax <- 3
expLS = exp(tmax*theta[1])

wherecov = which(diffmat<expLS)
h = matrix(0,dim(diffmat)[1],dim(diffmat)[2])
h[wherecov] = 1-diffmat[wherecov] / expLS
if (!returnlogs) {
C = matrix(0,dim(h)[1],dim(h)[2])
C[wherecov] <- h[wherecov]^3 * (4 - 3*h[wherecov])
} else {
C = -Inf * matrix(1,dim(h)[1],dim(h)[2])
C[wherecov] = 3* log(h[wherecov]) + log(4 - 3*h[wherecov])
}

if(return_dCdtheta){
h2 = 1-h
if (!returnlogs) {
# dCdtheta <- ifelse(1-h > 0,
#                    # tmax*expLS * (3*(1-h2)^2*(h2/expLS)*(3*h+1) + (1-h2)^3*(-3*h2/expLS)),
#                    tmax*expLS * (3*(1-h2)^2*(h2/expLS)*(4-3*h) + (1-h2)^3*(-3*1/expLS)),
#                    0)
dCdtheta = matrix(0,dim(diffmat)[1],dim(diffmat)[2])
# dCdtheta[wherecov] <- 12*tmax*h[wherecov]^2*(1-h[wherecov])
dCdtheta[wherecov] <- 12*tmax*h[wherecov]^2*(1-h[wherecov]) * diffmat[wherecov] / expLS
} else {
dCdtheta = matrix(0,dim(diffmat)[1],dim(diffmat)[2])
# dCdtheta[wherecov] <- 12*tmax*(1-h[wherecov])/(h[wherecov] * (4 - 3*h[wherecov]))
# dCdtheta[wherecov] <- 12*tmax*h[wherecov]^2*(1-h[wherecov]) * diffmat[wherecov] / expLS / C[wherecov]
dCdtheta[wherecov] <- tmax * 3 * (1/h[wherecov] - 1/(4-3*h[wherecov])) * diffmat[wherecov] / expLS
}
dCdtheta[is.na(dCdtheta)] = 0
out <- list(C=C,dCdtheta=dCdtheta)
return(out)
}else{
return(C)
}
}
}

#' Wendland2 2 correlation function
#'
#' Calculate correlation matrix for two sets of points in one dimension.
#' Note that this is not the correlation between two vectors.
#'
#' @inheritParams CGGP_internal_CorrMatCauchy
#'
#' @return Matrix of correlation values between x1 and x2
# @rdname CGGP_internal_CorrMatCauchy
#' @export
#' @family correlation functions
#'
#' @examples
#' CGGP_internal_CorrMatWendland2(c(0,.2,.4),c(.1,.3,.5), theta=-.7)
CGGP_internal_CorrMatWendland2 <- function(x1, x2,theta,
return_dCdtheta = FALSE,
return_numpara=FALSE,
returnlogs=FALSE) {
if(return_numpara){
return(1)
}else{
if (length(theta) != 1) {stop("CorrMatWendland2 theta should be length 1")}
diffmat =abs(outer(x1,x2,'-'))
tmax <- 3
expLS = exp(tmax*theta[1])

wherecov = which(diffmat<expLS)
h = matrix(0,dim(diffmat)[1],dim(diffmat)[2])
h[wherecov] = 1-diffmat[wherecov] / expLS
if (!returnlogs) {
C = matrix(0,dim(h)[1],dim(h)[2])
C[wherecov] <- h[wherecov]^5 * (8*h[wherecov]^2- 21*h[wherecov] + 14)
} else {
C = -Inf * matrix(1,dim(h)[1],dim(h)[2])
C[wherecov] = 5* log(h[wherecov]) + log(8*h[wherecov]^2- 21*h[wherecov] + 14)
}
if(return_dCdtheta){
if (!returnlogs) {
dCdtheta = matrix(0,dim(h)[1],dim(h)[2])
dCdtheta[wherecov] <- tmax*14*h[wherecov]^4*(1-h[wherecov])^2*(5-4*h[wherecov])
} else {
dCdtheta = matrix(0,dim(h)[1],dim(h)[2])
dCdtheta[wherecov] <- tmax*14*(1-h[wherecov])^2*(5-4*h[wherecov])/(h[wherecov]*(8*h[wherecov]^2- 21*h[wherecov] + 14))
}
dCdtheta[is.na(dCdtheta)] = 0
out <- list(C=C,dCdtheta=dCdtheta)
return(out)
}else{
return(C)
}
}
}

CollinErickson/CGGP documentation built on Feb. 6, 2024, 2:24 a.m.