#EGIparallel
set.seed(9)
N <- 20 #number of observations
T <- c(20,60) #thresholds
testfun <- branin
lower <- c(0,0)
upper <- c(1,1)
#a 20 points initial design
design <- data.frame( matrix(runif(2*N),ncol=2) )
response <- testfun(design)
#km object with matern3_2 covariance
#params estimated by ML from the observations
model <- km(formula=~., design = design,
response = response,covtype="matern3_2")
optimcontrol <- list(method="genoud",pop.size=50)
integcontrol <- list(distrib="sur",n.points=50)
iter <- 1
batchsize <- 6
if (FALSE) {
obj <- EGIparallel(T=T,model=model,method="sur",batchsize=batchsize,
fun=testfun,iter=iter,lower=lower,upper=upper,
optimcontrol=optimcontrol,integcontrol=integcontrol)
par(mfrow=c(1,2))
print_uncertainty_2d(model=model,T=T,main="probability of excursion",
type="pn",new.points=0,cex.points=2)
print_uncertainty_2d(model=obj$lastmodel,T=T,
main="probability of excursion, parallel sur sampling",
type="pn",new.points=iter*batchsize,col.points.end="red",cex.points=2)
}
##############
#same example with noisy initial observations and noisy new observations
branin.noise <- function(x) return(branin(x)+rnorm(n=1,sd=30))
set.seed(9)
N <- 20;T <- c(20,60)
testfun <- branin.noise
lower <- c(0,0);upper <- c(1,1)
design <- data.frame( matrix(runif(2*N),ncol=2) )
response.noise <- apply(design,1,testfun)
response.noise - response
model.noise <- km(formula=~., design = design, response = response.noise,
covtype="matern3_2",noise.var=rep(30*30,times=N))
optimcontrol <- list(method="genoud",pop.size=50)
integcontrol <- list(distrib="sur",n.points=50)
iter <- 1
batchsize <- 6
if (FALSE) {
obj <- EGIparallel(T=T,model=model.noise,method="sur",batchsize=batchsize,
fun=testfun,iter=iter,lower=lower,upper=upper,
optimcontrol=optimcontrol,integcontrol=integcontrol,
new.noise.var=10*10)
par(mfrow=c(1,2))
print_uncertainty_2d(model=model.noise,T=T,
main="probability of excursion, noisy obs.",
type="pn",new.points=0,cex.points=2)
print_uncertainty_2d(model=obj$lastmodel,T=T,
main="probability of excursion, parallel sur sampling, noisy obs.",
type="pn",new.points=iter*batchsize,col.points.end="red",cex.points=2)
}
##############
# Conservative estimates with non-noisy initial observations
if (FALSE) {
testfun <- branin
# The conservative sampling strategies
# only work with 1 threshold
T <- 20
## Minimize Type II error sampling
# The list method.param contains all parameters for the
# conservative estimate and the conservative sequential
# strategy. Below are parameters for a type II strategy
# with conservative estimates at 0.95
method.param = list(penalization=0, # Type II strategy
typeEx=">", consLevel = 0.95,
n_discrete_design=500*model@d)
# If the CE for the initial model is already computed
# it is possible to pass the conservative Vorob'ev quantile
# level with method.param$consVorbLevel
obj_T2 <- EGIparallel(T=T,model=model,method="vorobCons",batchsize=batchsize,
fun=testfun,iter=iter,lower=lower,upper=upper,
optimcontrol=optimcontrol,
integcontrol=integcontrol,method.param=method.param)
par(mfrow=c(1,2))
print_uncertainty_2d(model=model,T=T,main="probability of excursion",
type="pn",new.points=0,cex.points=2,consQuantile = obj_T2$allCE_lvs[1])
print_uncertainty_2d(model=obj_T2$lastmodel,T=T,
main="probability of excursion, parallel Type II sampling",
type="pn",new.points=iter*batchsize,col.points.end="red",
cex.points=2,consQuantile = obj_T2$allCE_lvs[2])
## Maximize conservative estimate's volume
# Set up method.param
# Here we pass the conservative level computed
# in the previous step for the initial model
method.param = list(typeEx=">", consLevel = 0.95,
n_discrete_design=500*model@d,
consVorbLevel=obj_T2$allCE_lvs[1]
)
obj_consVol <- EGIparallel(T=T,model=model,method="vorobVol",batchsize=batchsize,
fun=testfun,iter=iter,lower=lower,upper=upper,
optimcontrol=optimcontrol,
integcontrol=integcontrol,method.param=method.param)
par(mfrow=c(1,2))
print_uncertainty_2d(model=model,T=T,main="probability of excursion",
type="pn",new.points=0,cex.points=2,consQuantile = obj_consVol$allCE_lvs[1])
print_uncertainty_2d(model=obj_consVol$lastmodel,T=T,
main="probability of excursion, parallel consVol sampling",
type="pn",new.points=iter*batchsize,col.points.end="red",
cex.points=2,consQuantile = obj_consVol$allCE_lvs[2])
}
Run the code above in your browser using DataLab