Title: | Autoregressive Cokriging Models for Multifidelity Codes |
---|---|
Description: | For emulating multifidelity computer models. The major methods include univariate autoregressive cokriging and multivariate autoregressive cokriging. The autoregressive cokriging methods are implemented for both hierarchically nested design and non-nested design. For hierarchically nested design, the model parameters are estimated via standard optimization algorithms; For non-nested design, the model parameters are estimated via Monte Carlo expectation-maximization (MCEM) algorithms. In both cases, the priors are chosen such that the posterior distributions are proper. Notice that the uniform priors on range parameters in the correlation function lead to improper posteriors. This should be avoided when Bayesian analysis is adopted. The development of objective priors for autoregressive cokriging models can be found in Pulong Ma (2020) <DOI:10.1137/19M1289893>. The development of the multivariate autoregressive cokriging models with possibly non-nested design can be found in Pulong Ma, Georgios Karagiannis, Bledar A Konomi, Taylor G Asher, Gabriel R Toro, and Andrew T Cox (2019) <arXiv:1909.01836>. |
Authors: | Pulong Ma [aut, cre] |
Maintainer: | Pulong Ma <[email protected]> |
License: | GPL (>= 2) |
Version: | 0.1.2 |
Built: | 2024-10-31 06:38:17 UTC |
Source: | CRAN |
This is a simple and high-level funciton to fit autoregressive cokriging models to multifidelity computer model outputs.
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() )
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() )
formula |
a list of |
output |
a list of |
input |
a list of |
cov.model |
a string indicating the type of covariance function in AR-cokriging models. Current covariance functions include
|
nugget.est |
a logical value indicating whether nugget parameter is included or not. Default value is |
input.new |
a matrix including new inputs for making prediction |
prior |
a list of arguments to setup the prior distributions
|
opt |
a list of arguments to setup the |
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
|
info |
a list that contains
|
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.
Pulong Ma <[email protected]>
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.
cokm
, cokm.param
, cokm.fit
, cokm.predict
############################################################## ############################################################## ############ 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 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)
############################################################## ############################################################## ############ 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 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)
This function constructs the cokm object in autogressive cokriging models
cokm( formula = list(~1, ~1), output, input, cov.model = "matern_5_2", nugget.est = FALSE, prior = list(), opt = list(), NestDesign = TRUE, tuning = list(), info = list() )
cokm( formula = list(~1, ~1), output, input, cov.model = "matern_5_2", nugget.est = FALSE, prior = list(), opt = list(), NestDesign = TRUE, tuning = list(), info = list() )
formula |
a list of |
output |
a list of |
input |
a list of |
cov.model |
a string indicating the type of covariance function in AR-cokriging models. Current covariance functions include
|
nugget.est |
a logical value indicating whether the nugget is included or not. Default value is |
prior |
a list of arguments to setup the prior distributions with the reference prior as default.
|
opt |
a list of arguments to setup the |
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
|
info |
a list that contains
|
Pulong Ma <[email protected]>
ARCokrig
, cokm.fit
, cokm.predict
This is an S4 class definition for cokm
in the ARCokrig
package
output
a list of elements, each of which contains a matrix of computer model outputs.
input
a list of elements, each of which contains a matrix of inputs.
param
a list of elements, each of which contains a vector of initial values for
correlation parameters (and nugget variance parameters if
nugget terms are included in AR-cokriging models).
cov.model
a string indicating the type of covariance function in AR-cokriging models. Current covariance functions include
product form of exponential covariance functions.
product form of Matern covariance functions with smoothness parameter 3/2.
product form of Matern covariance functions with smoothness parameter 5/2.
product form of Gaussian covariance functions.
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
.
prior
a list of arguments to setup the prior distributions with the reference prior as default
the name of the prior. Current implementation includes
JR
, Reference
, Jeffreys
, Ind_Jeffreys
hyperparameters in the priors.
For jointly robust (JR) prior, three parameters are included:
refers to the polynomial penalty to avoid singular correlation
matrix with a default value 0.2;
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
the maximum number of MCEM iterations.
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.
the number of Monte Carlo samples in the MCEM algorithm.
a logical value to show the MCEM iterations if it is true.
info
a list that contains
number of iterations used in the MCEM algorithm
parameter difference after the MCEM algorithm stops.
Pulong Ma <[email protected]>
This function simulate from predictive distributions in autogressive cokriging models
cokm.condsim(obj, input.new, nsample = 30)
cokm.condsim(obj, input.new, nsample = 30)
obj |
a |
input.new |
a matrix including new inputs for making prediction |
nsample |
a numerical value indicating the number of samples |
Pulong Ma <[email protected]>
cokm
, cokm.fit
, cokm.predict
, ARCokrig
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)) ## create the cokm object prior = list(name="Reference") obj = cokm(formula=list(~1,~1+x1), output=list(c(zc), c(zf)), input=list(as.matrix(Dc), as.matrix(Df)), prior=prior, cov.model="matern_5_2") ## update model parameters in the cokm object obj = cokm.fit(obj) cokrige = cokm.condsim(obj, input.new, nsample=30)
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)) ## create the cokm object prior = list(name="Reference") obj = cokm(formula=list(~1,~1+x1), output=list(c(zc), c(zf)), input=list(as.matrix(Dc), as.matrix(Df)), prior=prior, cov.model="matern_5_2") ## update model parameters in the cokm object obj = cokm.fit(obj) cokrige = cokm.condsim(obj, input.new, nsample=30)
This function estimates parameters in autogressive cokriging models
cokm.fit(obj)
cokm.fit(obj)
obj |
a |
Pulong Ma <[email protected]>
cokm
, cokm.param
, cokm.predict
, ARCokrig
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)) ## create the cokm object prior = list(name="JR") obj = cokm(formula=list(~1,~1+x1), output=list(c(zc), c(zf)), input=list(as.matrix(Dc), as.matrix(Df)), prior=prior, cov.model="matern_5_2") ## update model parameters in the cokm object obj = cokm.fit(obj)
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)) ## create the cokm object prior = list(name="JR") obj = cokm(formula=list(~1,~1+x1), output=list(c(zc), c(zf)), input=list(as.matrix(Dc), as.matrix(Df)), prior=prior, cov.model="matern_5_2") ## update model parameters in the cokm object obj = cokm.fit(obj)
This function compute estimates for regression and variance parameters given the correlation parameters are known. It is used to show all model parameters in one place.
cokm.param(obj)
cokm.param(obj)
obj |
a |
a list of model parameters including regression coefficients ,
scale discrepancy
, variance parameters
, and correlation parameters
in covariance functions.
If nugget parameters are included in the model, then nugget parameters are shown in
.
Pulong Ma <[email protected]>
cokm
, cokm.fit
, cokm.condsim
, ARCokrig
This function makes prediction in autogressive cokriging models. If a nested design is used, the predictive mean and predictive variance are computed exactly; otherwise, Monte Carlo simulation from the predictive distribution is used to approximate the predictive mean and predictive variance.
cokm.predict(obj, input.new)
cokm.predict(obj, input.new)
obj |
a |
input.new |
a matrix including new inputs for making prediction |
Pulong Ma <[email protected]>
cokm
, cokm.fit
, cokm.condsim
, ARCokrig
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)) ## create the cokm object prior = list(name="Reference") obj = cokm(formula=list(~1,~1+x1), output=list(c(zc), c(zf)), input=list(as.matrix(Dc), as.matrix(Df)), prior=prior, cov.model="matern_5_2") ## update model parameters in the cokm object obj = cokm.fit(obj) cokrige = cokm.predict(obj, input.new)
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)) ## create the cokm object prior = list(name="Reference") obj = cokm(formula=list(~1,~1+x1), output=list(c(zc), c(zf)), input=list(as.matrix(Dc), as.matrix(Df)), prior=prior, cov.model="matern_5_2") ## update model parameters in the cokm object obj = cokm.fit(obj) cokrige = cokm.predict(obj, input.new)
This function compute the continous rank probability score for normal distributions. It is mainly used to evaluate the validility of predictive distributions.
CRPS(x, mu, sig)
CRPS(x, mu, sig)
x |
a vector of true values (held-out data) |
mu |
a vector of predictive means |
sig |
a vector of predictive standard deviations |
Pulong Ma <[email protected]>
This function constructs the mvcokm object in autogressive cokriging models for multivariate outputs. The model is known as the parallel partial (PP) cokriging emulator.
mvcokm( formula = list(~1, ~1), output, input, cov.model = "matern_5_2", nugget.est = FALSE, prior = list(), opt = list(), NestDesign = TRUE, tuning = list(), info = list() )
mvcokm( formula = list(~1, ~1), output, input, cov.model = "matern_5_2", nugget.est = FALSE, prior = list(), opt = list(), NestDesign = TRUE, tuning = list(), info = list() )
formula |
a list of |
output |
a list of |
input |
a list of |
cov.model |
a string indicating the type of covariance function in the PP cokriging models. Current covariance functions include
|
nugget.est |
a logical value indicating whether the nugget is included or not. Default value is |
prior |
a list of arguments to setup the prior distributions with the jointly robust prior as default
|
opt |
a list of arguments to setup the |
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
|
info |
a list that contains
|
Pulong Ma <[email protected]>
ARCokrig
, mvcokm.fit
, mvcokm.predict
, mvcokm.condsim
This is an S4 class definition for mvcokm
in the ARCokrig
package
output
a list of elements, each of which contains a matrix of computer model outputs.
input
a list of elements, each of which contains a matrix of inputs.
param
a list of elements, each of which contains a vector of initial values for
correlation parameters (and nugget variance parameters if
nugget terms are included in AR-cokriging models).
cov.model
a string indicating the type of covariance function in AR-cokriging models. Current covariance functions include
product form of exponential covariance functions.
product form of Matern covariance functions with smoothness parameter 3/2.
product form of Matern covariance functions with smoothness parameter 5/2.
product form of Gaussian covariance functions.
product form of power-exponential covariance functions with roughness parameter fixed at 1.9.
anisotropic form of exponential covariance function.
anisotropic form of Matern covariance functions with smoothness parameter 3/2.
anisotropic form of Matern covariance functions with smoothness parameter 5/2.
nugget.est
a logical value indicating whether the nugget is included or not. Default value is FALSE
.
prior
a list of arguments to setup the prior distributions with the jointly robust prior as default
the name of the prior. Current implementation includes
JR
, Reference
, Jeffreys
, Ind_Jeffreys
hyperparameters in the priors.
For jointly robust (JR) prior, three parameters are included:
refers to the polynomial penalty to avoid singular correlation
matrix with a default value 0.2;
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
the maximum number of MCEM iterations.
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.
the number of Monte Carlo samples in the MCEM algorithm.
info
a list that contains
number of iterations used in the MCEM algorithm
parameter difference after the MCEM algorithm stops
Pulong Ma <[email protected]>
This function makes prediction based on conditional simulation in autogressive cokriging models for multivariate output
mvcokm.condsim(obj, input.new, nsample = 30)
mvcokm.condsim(obj, input.new, nsample = 30)
obj |
a |
input.new |
a matrix including new inputs for making prediction |
nsample |
a numerical value indicating the number of samples |
Pulong Ma <[email protected]>
mvcokm
, mvcokm.fit
, mvcokm.predict
, ARCokrig
This function estimates parameters in the parallel partial cokriging model
mvcokm.fit(obj)
mvcokm.fit(obj)
obj |
a |
Pulong Ma <[email protected]>
mvcokm
, mvcokm.predict
, mvcokm.condsim
, ARCokrig
This function computes estimates for regression and variance parameters given the correlation parameters are known. It is used to show all model parameters in one place.
mvcokm.param(obj)
mvcokm.param(obj)
obj |
a |
a list of model parameters including regression coefficients ,
scale discrepancy
, variance parameters
, and correlation parameters
in covariance functions.
If nugget parameters are included in the model, then nugget parameters are shown in
.
Pulong Ma <[email protected]>
mvcokm
, mvcokm.fit
, mvcokm.predict
, ARCokrig
This function makes prediction in the parallel partial cokriging model. If a nested design is used, the predictive mean and predictive variance are computed exactly; otherwise, Monte Carlo simulation from the predictive distribution is used to approximate the predictive mean and predictive variance.
mvcokm.predict(obj, input.new)
mvcokm.predict(obj, input.new)
obj |
a |
input.new |
a matrix including new inputs for making prediction |
Pulong Ma <[email protected]>
mvcokm
, mvcokm.fit
, mvcokm.condsim
, ARCokrig