| GPdiffMesh {diffMeshGP} | R Documentation |
This R function implements the nonstationary Kriging model proposed by Tuo, Wu and Yu (2014) for analyzing multi-fidelity computer outputs. This function computes the maximum likelihood estimates for the model parameters as well as the predictive means and variances of the exact solution (i.e., the conceptually highest fidelity).
GPdiffMesh(x,X,meshT,Y,regFunX = function(x){ return(0*matrix(x[,1]))},
regFunT = function(x){ return(1*matrix(x[,1]))},
phi1 = 1,
sigma12 = 1,
sigma22 = 1,
phi2 = 1,
mybeta = FALSE,
l = 4)
x |
An (n_0 x d) matrix of control points on which the values of the exact solution will be predicted. Each row represents one input point. |
X |
An (n x d) matrix of control variables. |
meshT |
An (n x 1) mesh density matrix. A same row of X and meshT represents the inputs from one trial of the computer experiment. |
Y |
An (n x 1) computer output matrix corresponding to (X, meshT). |
regFunX |
A scalar or a vector-valued regression function f_1(x). The default value is f_1(x)=1. See Details below. |
regFunT |
A scalar or a vector-valued regression function f_2(x,t). The default value is f_2(x,t) = t. See Details below. |
phi1,sigma12,sigma22,phi2,mybeta |
Initial values of the parameters for the maximum likelihood estimation. The default values are phi1 = 1, sigma12 = 1, sigma22 = 1, phi2 = 1, mybeta = 1. See Details below. |
l |
The value of parameter l in the covariance function of the nonstationary Kriging model. (We used fixed l in this R function.) The default value is l=4. See Details below. |
This R function implements the nonstationary Kriging model proposed by Tuo, Wu and Yu (2014) for the modeling and analysis of multi-fidelity computer experiments. Denote an input-output pair from the computer simulation by (x, t, y), where x is the vector of input variables, t is the mesh density, and y is the corresponding computer output. Tuo, Wu and Yu (2014) use the following Gaussian process model:
y = η(x, t) = η(x, 0) + δ(x, t),
where η(x, 0) and δ(x, t) are realizations of two mutually independent Gaussian processes V(x) and Z(x, t). Assume
E(V(x))=f_{1}^{T}(x)β_1, E(Z(x, t))=f^{T}_{2}(x,t)β_2,
where f_{1}(x) can be set in regFunX, and f_{2}(x,t) can be set in regFunT. Both f_{1}(x) and f_{2}(x,t) can be vector-valued function. f_{2}(x,t) should satisfy for any x, \lim_{t\rightarrow 0}f_{2}(x,t) = 0. The default functions are f_{1}(x) = 1 and f_{2}(x,t) = t. Assume the covariance of V(x) is
Cov(V(x_{1}),V(x_{2}))=σ_{1}^{2}e^{-φ^2_1\|x_{1}-x_{2}\|^{2}_{2}},
and the covariance of Z(x, t) is
Cov(Z(x_{1}, t_{1}),Z(x_{2}, t_{2}))=σ_{2}^{2}\min(t_{1},t_{2})^{l}e^{-φ_{2}^{2}\|x_{1}-x_{2}\|^{2}_{2}},
where l is a fixed parameter and is not estimated in this function.
The goal of the Tuo-Wu-Yu model is to predict for the exact solution \varphi(.):=η(.,0). For a set of control points x=(x_{1},…,x_{n_0}), this R function predicts η(x_{i},0)=\varphi(x_{i}) for i=1,…,n_{0}.
This R function uses maximum likelihood method to estimate the model parameters. Nelder-Mead method is used to maximize the likelihood function. The solution may depend on the choice of initial values. Users may specify the initial values or use the default values.
outy |
An (n_0 x 1) matrix of predictive means corresponding to control points. A same row of x and outy represents a pair of one control point and corresponding predictive mean. |
sigy |
An (n_0 x 1) matrix of predictive variances corresponding to control points. A same row of x and sigy represents a pair of one control point and corresponding predictive variances. |
estipar |
A list of parameter estimates. |
hig02 <- function(s)
{
#The test function is from [2].
y <- s*sin(s) / 10
return(y)
}
myX <- matrix(c(seq(from = 0,to = 10, by = 1),
seq(from = 0,to = 10, by = 1)),ncol = 2)
myy <- hig02(matrix(myX[,1]))
myT <- matrix(c(0.01,0.5,0.01,0.02,0.02,0.01,0.01,0.02,0.002,
0.003,0.03))
myregf <- function(x){
return(x)
}
myregfn <- function(s){
return(cbind((matrix(s[,1])^2*matrix(s[,2])),
(matrix(s[,1])*matrix(s[,2]))))
}
#Here s=cbind(t,x), where x is a matrix of input variables and
# t is the corresponding mesh density matrix.
x <- matrix(c(seq(from = 0,to = 10, by = 0.1),
seq(from = 0,to = 10, by = 0.1)),ncol = 2)
myploty <- hig02(matrix(x[,1]))
y <- GPdiffMesh(x, myX, myT, myy, regFunX = myregf, regFunT = myregfn)
# The regression function is beta_0 + beta_1 x + (t^2x, tx)' beta_2
y$outy
y$sigy
y$estipar
plot(x[,1], myploty,"l")
lines(x[,1],y$outy, type="o", pch=22, lty=2, col="red")