sequentialDesign: Calibration with a sequential design

In mathieucarmassi/calibrationcode: Code Calibration in a Bayesian Framework

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 NgibbsNumber of iteration of the algorithm Metropolis within Gibbs 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])

## End(Not run)

mathieucarmassi/calibrationcode documentation built on Aug. 14, 2019, 12:35 a.m.