Package 'diffMeshGP'

Multi-Fidelity Computer Experiments Using The Tuo-Wu-Yu Model

Description

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).

Usage

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)

Arguments

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.

Details

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 = \eta(x, t) = \eta(x, 0) + \delta(x, t),$

where $\eta(x, 0)$ and $\delta(x, t)$ are realizations of two mutually independent Gaussian processes $V(x)$ and $Z(x, t)$. Assume

$E(V(x))=f_{1}^{T}(x)\beta_1, E(Z(x, t))=f^{T}_{2}(x,t)\beta_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}))=\sigma_{1}^{2}e^{-\phi^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}))=\sigma_{2}^{2}\min(t_{1},t_{2})^{l}e^{-\phi_{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(.):=\eta(.,0)$. For a set of control points $x=(x_{1},\ldots,x_{n_0})$, this R function predicts $\eta(x_{i},0)=\varphi(x_{i})$ for $i=1,\ldots,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.

Value

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.

Examples

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")

---