saEI {DynamicGP}

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