saEI: Saddlepoint Approximate Expected Improvement Criterion for...

In DynamicGP: Modelling and Analysis of Dynamic Computer Experiments

Saddlepoint Approximate Expected Improvement Criterion for the Sequential Design for Inverse Problems

Description

This function performs the sequential design procedure for the inverse problem. It starts from an initial design set xi and selects the follow-up design points from the candidate set candei as per the expected improvement (EI) criterion which is numerically approximated by the saddlepoint approximation technique in Huang and Oosterlee (2011). The surrogate is refitted using the augmented data via svdGP. After the selection of nadd follow-up points, the solution of the inverse problem is estimated either by the ESL2D approach or by the SL2D approach. Details are provided in Chapter 4 of Zhang (2018).

Usage

saEI(xi,yi,yobs,nadd,candei,candest,func,...,
     mtype=c("zmean","cmean","lmean"),
     estsol=c("ESL2D","SL2D"),
     frac=.95, nstarts=5, gstart=0.0001,
     nthread=1, clutype="PSOCK")

Arguments

xi	An N0 by d matrix of N0 initial design points.
yi	An L by N0 response matrix of xi, where L is the length of the time series outputs, N0 is the number of design points.
yobs	A vector of length L of the time-series valued field observations or the target response.
nadd	The number of the follow-up design points selected by this function.
candei	An M1 by d matrix of M1 candidate points on which the follow-up design points are selected.
candest	An M2 by d matrix of M2 candidate points on which the (final) estimated solution to the inverse problem is extracted.
func	An R function of the dynamic computer simulator. The first argument of func should be a vector of d-dimensional inputs. The simulator func should return a vector of length L as the output.
...	The remaining arguments of the simulator func.
mtype	The type of mean functions for the GP models. The choice "zmean" denotes zero-mean, "cmean" indicates constant-mean, "lmean" indicates linear-mean. The default choice is "zmean".
estsol	The method for estimating the final solution to the inverse problem after all follow-up design points are included, "ESL2D" denotes the ESL2D approach, "SL2D" denotes the SL2D approach. The default choice is "ESL2D".
frac	The threshold in the cumulative percentage criterion to select the number of SVD bases. The default value is 0.95.
nstarts	The number of starting points used in the numerical maximization of the posterior density function. The larger nstarts will typically lead to more accurate prediction but longer computational time. The default value is 5.
gstart	The starting number and upper bound for estimating the nugget parameter. If gstart = sqrt(.Machine$double.eps), the nugget parameter will be fixed at sqrt(.Machine$double.eps), since sqrt(.Machine$double.eps) is the lower bound of the nugget term. The default value is 0.0001.
nthread	The number of threads (processes) used in parallel execution of this function. nthread=1 implies no parallelization. The default value is 1.
clutype	The type of cluster in the R package "parallel" to perform parallelization. The default value is "PSOCK". Required only if nthread>1.

Value

xx	The design set selected by the sequential design approach, which includes both the initial and the follow-up design points.
yy	The response matrix collected on the design set xx.
xhat	The estimated solution to the inverse problem obtained on the candidate set candest from the final fitted surrogate.
maxei	A vector of length nadd, it collects the maximum value of the EI criterion in each iteration of the sequential design approach.

Author(s)

Ru Zhang [e-mail],

C. Devon Lin [e-mail],

Pritam Ranjan [e-mail]

References

Huang, X. and Oosterlee, C. W. (2011) Saddlepoint approximations for expectations and an application to CDO pricing, SIAM Journal on Financial Mathematics, 2(1) 692-714.

Zhang, R. (2018) Modeling and Analysis of Dynamic Computer Experiments, PhD thesis, Queen's University, ON, Canada.

See Also

ESL2D, SL2D, svdGP.

Examples

  library("lhs")
  forretal <- function(x,t,shift=1)
  {
    par1 <- x[1]*6+4
    par2 <- x[2]*16+4
    par3 <- x[3]*6+1
    t <- t+shift
    y <- (par1*t-2)^2*sin(par2*t-par3)
  }
  timepoints <- seq(0,1,len=200)
  xi <- lhs::randomLHS(30,3)
  candei <- lhs::randomLHS(500,3)
  candest <- lhs::randomLHS(500,3)
  candest <- rbind(candest, xi)

  ## evaluate the response matrix on the design matrix
  yi <- apply(xi,1,forretal,timepoints)
  x0 <- runif(3)
  y0 <- forretal(x0,timepoints)
  yobs <- y0+rnorm(200,0,sd(y0)/sqrt(50))
  ret <- saEI(xi,yi,yobs,1,candei,candest,forretal,timepoints,
              nstarts=1, nthread=1)
  yhat <- forretal(ret$xhat,timepoints)

  ## draw a figure to illustrate
  plot(y0,ylim=c(min(y0,yhat),max(y0,yhat)))
  lines(yhat,col="red")

DynamicGP documentation built on Nov. 10, 2022, 5:15 p.m.