CaliCo (version 0.1.1)

sequentialDesign: Calibration with a sequential design

Description

The aim is to reduce the error produced by the initial estimation of the Gaussian process by fortifying the initial DOE. The method consists in proposing new points based on the expectancy improvement criterion. The method and the algorithm are detailed in [Damblin et al. 2018]

Usage

sequentialDesign(md, pr, opt.estim, k)

Arguments

md

the model to improve (model2 or model4)

pr

list of priors to use for calibration

opt.estim

estimation options

• Nmh Number of iteration of the Metropolis Hastings algorithm

• thetaInit Initial point

• r regulation percentage in the modification of the k in the Metropolis Hastings

• sig Covariance matrix for the proposition distribution (\(k*sig\))

• Nchains Number of MCMC chains to run (if Nchain>1 an output is created called mcmc which is a coda object codamenu)

• burnIn Number of iteration to withdraw

k

number of iteration in the algorithm

Value

a seqDesign.class

References

DAMBLIN, Guillaume, BARBILLON, Pierre, KELLER, Merlin, et al. Adaptive numerical designs for the calibration of computer codes. SIAM/ASA Journal on Uncertainty Quantification, 2018, vol. 6, no 1, p. 151-179.

See Also

model, prior, calibrate, sequentialDesign

Examples

# NOT RUN {
###### The code to calibrate
X <- cbind(seq(0,1,length.out=5),seq(0,1,length.out=5))
code <- function(X,theta)
{
  return((6*X[,1]*theta[2]-2)^2*theta[1]*sin(theta[3]*X[,2]-4))
}
Yexp <- code(X,c(1,1,11))+rnorm(5,0,0.1)

###### For the second model
### code function is available, no DOE generated upstream
binf <- c(0.9,0.9,10.5)
bsup <- c(1.1,1.1,11.5)
opt.gp <- list(type="matern5_2", DOE=NULL)
opt.emul <- list(p=3,n.emul=150,binf=binf,bsup=bsup,type="maximinLHS")
model2 <- model(code,X,Yexp,"model2",
                opt.gp=opt.gp,
                opt.emul=opt.emul)
model2 %<% list(theta=c(1,1,11),var=0.1)

pr <- prior(type.prior=c("gaussian","gaussian","gaussian","gamma"),opt.prior=
list(c(1,0.01),c(1,0.01),c(11,3),c(2,0.1)))
###### Definition of the calibration options
opt.estim=list(Ngibbs=200,Nmh=400,thetaInit=c(1,1,11,0.1),r=c(0.3,0.3),
sig=diag(4),Nchains=1,burnIn=100)
###### Run the sequential calibration
mdfit <- sequentialDesign(model2,pr,opt.estim,2)
#plot(mdfit,X[,1])
# }