ARCokrig (version 0.1.2)

ARCokrig: Fit the AR-Cokriging model and make predictions

Description

This is a simple and high-level funciton to fit autoregressive cokriging models to multifidelity computer model outputs.

Usage

ARCokrig(
  formula = list(~1, ~1),
  output,
  input,
  cov.model = "matern_5_2",
  nugget.est = FALSE,
  input.new,
  prior = list(),
  opt = list(),
  NestDesign = TRUE,
  tuning = list(),
  info = list()
)

Arguments

formula

a list of \(s\) elements, each of which contains the formula to specify fixed basis functions or regressors.

output

a list of \(s\) elements, each of which contains a matrix of computer model outputs.

input

a list of \(s\) elements, each of which contains a matrix of inputs.

cov.model

a string indicating the type of covariance function in AR-cokriging models. Current covariance functions include

exp

product form of exponential covariance functions.

matern_3_2

product form of Matern covariance functions with smoothness parameter 3/2.

matern_5_2

product form of Matern covariance functions with smoothness parameter 5/2.

Gaussian

product form of Gaussian covariance functions.

powexp

product form of power-exponential covariance functions with roughness parameter fixed at 1.9.

nugget.est

a logical value indicating whether nugget parameter is included or not. Default value is FALSE.

input.new

a matrix including new inputs for making prediction

prior

a list of arguments to setup the prior distributions

name

the name of the prior. Current implementation includes JR, Reference, Jeffreys, Ind_Jeffreys

hyperparam

hyperparameters in the priors. For jointly robust (JR) prior, three parameters are included: \(a\) refers to the polynomial penalty to avoid singular correlation matrix with a default value 0.2; \(b\) refers to the exponenetial penalty to avoid diagonal correlation matrix with a default value 1; nugget.UB is the upper bound of the nugget variance with default value 1, which indicates that the nugget variance has support (0, 1).

opt

a list of arguments to setup the optim routine.

NestDesign

a logical value indicating whether the experimental design is hierarchically nested within each level of the code.

tuning

a list of arguments to control the MCEM algorithm for non-nested design. It includes the arguments

maxit

the maximum number of MCEM iterations.

tol

a tolerance to stop the MCEM algorithm. If the parameter difference between any two consecutive MCEM algorithm is less than this tolerance, the MCEM algorithm is stopped.

n.sample

the number of Monte Carlo samples in the MCEM algorithm.

verbose

a logical value to show the MCEM iterations if it is true.

info

a list that contains

iter

number of iterations used in the MCEM algorithm

eps

parameter difference after the MCEM algorithm stops

Value

The main call inside ARCokrig consists of cokm, cokm.fit, and cokm.predict. Thus, the function returns the cokm object and predictions over new inputs.

References

• Ma, P. (2019). ``Objective Bayesian Analysis of a Cokriging Model for Hierarchical Multifidelity Codes." arXiv:1910.10225. https://arxiv.org/abs/1910.10225.

• Ma, P., Karagiannis, G., Konomi, B., Asher, T., Toro, G., and Cox, A. (2019) ``Multifidelity Computer Model Emulation with High-Dimensional Output: An Application to Storm Surge." arXiv:1909.01836. https://arxiv.org/abs/1909.01836.

See Also

cokm, cokm.param, cokm.fit, cokm.predict

Examples

# NOT RUN {
##############################################################
##############################################################
############ Example
Funcc = function(x){
  return(0.5*(6*x-2)^2*sin(12*x-4)+10*(x-0.5)-5)
}

Funcf = function(x){
  z1 = Funcc(x)
  z2 = 2*z1-20*x+20 + sin(10*cos(5*x))
  return(z2)
}

#####################################################################
###### Nested design 
#####################################################################
Dc <- seq(-1,1,0.1)
indDf <- c(1, 3, 6, 8, 10, 13, 17, 21)
zc <- Funcc(Dc)
Df <- Dc[indDf]
zf <- Funcf(Df)

input.new = as.matrix(seq(-1,1,length.out=200))


## fit and predict with the AR-Cokriging model

out = ARCokrig(formula=list(~1,~1+x1), output=list(c(zc), c(zf)),
              input=list(as.matrix(Dc), as.matrix(Df)),
              cov.model="matern_5_2",
              input.new=input.new)

## plot results
# }
# NOT RUN {
library(ggplot2)
cokrig = out$cokrig
df.l1 = data.frame(x=c(Dc), y=c(zc))
df.l2 = data.frame(x=c(Df), y=c(zf))
CI.lower = cokrig$lower95[[2]]
CI.upper = cokrig$upper95[[2]]
df.CI = data.frame(x=c(input.new),lower=CI.lower, upper=CI.upper)
df.pred = data.frame(x=c(input.new), y=cokrig$mu[[2]])

g = ggplot(data.frame(x=c(-1,1)), aes(x)) + 
 stat_function(fun=Funcc, geom="line", aes(colour="level 1"), n=500) +
 stat_function(fun=Funcf, geom="line", aes(colour="level 2"), n=500) 

g = g + geom_point(data=df.l1, mapping=aes(x=x, y=y), shape=16, size=2, color="black") + 
 geom_point(data=df.l2, mapping=aes(x=x, y=y), shape=17, size=2, color="black")

g = g + geom_line(data=df.pred, aes(x=x, y=y, colour="cokriging"), inherit.aes=FALSE) +
 geom_ribbon(data=df.CI, mapping=aes(x=x,ymin=lower, ymax=upper), fill="gray40",
             alpha=0.3, inherit.aes=FALSE)
g = g + scale_colour_manual(name=NULL, values=c("red","blue", "turquoise3"), 
                           breaks=c("cokriging","level 1", "level 2"))  

g = g + ggtitle("A Two-Level Example") + 
 theme(plot.title=element_text(size=14),
       axis.title.x=element_text(size=14),
       axis.text.x=element_text(size=14),
       axis.title.y=element_text(size=14),
       axis.text.y=element_text(size=14),
       legend.text = element_text(size=12),
       legend.direction = "horizontal",
       legend.position = c(0.6, 0.1)) + xlab("") + ylab("")
print(g)

# }
# NOT RUN {


# }