##########################################################################
### EXAMPLE 1: 3 OPTIMIZATION STEPS USING EQI WITH KNOWN NOISE ###
### AND KNOWN COVARIANCE PARAMETERS FOR THE BRANIN FUNCTION ###
##########################################################################
set.seed(10)
library(DiceDesign)
# Set test problem parameters
doe.size <- 9
dim <- 2
test.function <- get("branin2")
lower <- rep(0,1,dim)
upper <- rep(1,1,dim)
noise.var <- 0.1
# Build noisy simulator
funnoise <- function(x)
{ f.new <- test.function(x) + sqrt(noise.var)*rnorm(n=1)
return(f.new)}
# Generate DOE and response
doe <- as.data.frame(lhsDesign(doe.size, dim)$design)
y.tilde <- funnoise(doe)
# Create kriging model
model <- km(y~1, design=doe, response=data.frame(y=y.tilde),
covtype="gauss", noise.var=rep(noise.var,1,doe.size),
lower=rep(.1,dim), upper=rep(1,dim), control=list(trace=FALSE))
# Optimisation with noisy.optimizer (n.ite can be increased)
n.ite <- 2
optim.param <- list()
optim.param$quantile <- .9
optim.result <- noisy.optimizer(optim.crit="EQI", optim.param=optim.param, model=model,
n.ite=n.ite, noise.var=noise.var, funnoise=funnoise, lower=lower, upper=upper,
NoiseReEstimate=FALSE, CovReEstimate=FALSE)
new.model <- optim.result$model
best.x <- optim.result$best.x
new.doe <- optim.result$history.x
## Not run:
# ##### DRAW RESULTS #####
# # Compute actual function on a grid
# n.grid <- 12
# x.grid <- y.grid <- seq(0,1,length=n.grid)
# design.grid <- expand.grid(x.grid, y.grid)
# names(design.grid) <- c("V1","V2")
# nt <- nrow(design.grid)
# func.grid <- rep(0,1,nt)
#
# for (i in 1:nt)
# { func.grid[i] <- test.function(x=design.grid[i,])}
#
# # Compute initial and final kriging on a grid
# pred <- predict(object=model, newdata=design.grid, type="UK", checkNames = FALSE)
# mk.grid1 <- pred$m
# sk.grid1 <- pred$sd
#
# pred <- predict(object=new.model, newdata=design.grid, type="UK", checkNames = FALSE)
# mk.grid2 <- pred$m
# sk.grid2 <- pred$sd
#
# # Plot initial kriging mean
# z.grid <- matrix(mk.grid1, n.grid, n.grid)
# filled.contour(x.grid,y.grid, z.grid, nlevels=50, color = topo.colors,
# plot.axes = {title("Initial kriging mean");
# points(model@X[,1],model@X[,2],pch=17,col="black");
# axis(1); axis(2)})
#
# # Plot initial kriging variance
# z.grid <- matrix(sk.grid1^2, n.grid, n.grid)
# filled.contour(x.grid,y.grid, z.grid, nlevels=50, color = topo.colors,
# plot.axes = {title("Initial kriging variance");
# points(model@X[,1],model@X[,2],pch=17,col="black");
# axis(1); axis(2)})
#
# # Plot final kriging mean
# z.grid <- matrix(mk.grid2, n.grid, n.grid)
# filled.contour(x.grid,y.grid, z.grid, nlevels=50, color = topo.colors,
# plot.axes = {title("Final kriging mean");
# points(new.model@X[,1],new.model@X[,2],pch=17,col="black");
# axis(1); axis(2)})
#
# # Plot final kriging variance
# z.grid <- matrix(sk.grid2^2, n.grid, n.grid)
# filled.contour(x.grid,y.grid, z.grid, nlevels=50, color = topo.colors,
# plot.axes = {title("Final kriging variance");
# points(new.model@X[,1],new.model@X[,2],pch=17,col="black");
# axis(1); axis(2)})
#
# # Plot actual function and observations
# z.grid <- matrix(func.grid, n.grid, n.grid)
# tit <- "Actual function - Black: initial points; red: added points"
# filled.contour(x.grid,y.grid, z.grid, nlevels=50, color = topo.colors,
# plot.axes = {title(tit);points(model@X[,1],model@X[,2],pch=17,col="black");
# points(new.doe[1,],new.doe[2,],pch=15,col="red");
# axis(1); axis(2)})
# ## End(Not run)
##########################################################################
### EXAMPLE 2: 3 OPTIMIZATION STEPS USING EQI WITH UNKNOWN NOISE ###
### AND UNKNOWN COVARIANCE PARAMETERS FOR THE BRANIN FUNCTION ###
##########################################################################
# Same initial model and parameters as for example 1
n.ite <- 2 # May be changed to a larger value
res <- noisy.optimizer(optim.crit="min.quantile",
optim.param=list(type="quantile",quantile=0.01),
model=model, n.ite=n.ite, noise.var=noise.var, funnoise=funnoise,
lower=lower, upper=upper,
control=list(print.level=0),CovReEstimate=TRUE, NoiseReEstimate=TRUE)
# Plot actual function and observations
plot(model@X[,1], model@X[,2], pch=17,xlim=c(0,1),ylim=c(0,1))
points(res$history.x[1,], res$history.x[2,], col="blue")
# Restart: requires the output estim.model of the previous run
# to deal with potential repetitions
res2 <- noisy.optimizer(optim.crit="min.quantile",
optim.param=list(type="quantile",quantile=0.01),
model=res$model, n.ite=n.ite, noise.var=noise.var, funnoise=funnoise,
lower=lower, upper=upper, estim.model=res$estim.model,
control=list(print.level=0),CovReEstimate=TRUE, NoiseReEstimate=TRUE)
# Plot new observations
points(res2$history.x[1,], res2$history.x[2,], col="red")
Run the code above in your browser using DataLab