These R files contains an analytical implementation of the criterion “IECI” proposed for Optimization under unknown constraints (2010) by Robert Gramacy and Herbert Lee .

This implementation is based on DiceKriging R package (CRAN version):

if ( ! is.element ( "DiceKriging" , installed.packages ())) install.packages ( "DiceKriging" )
# More robust (and expensive) EI optimization
if ( ! is.element ( "rgenoud" , installed.packages ())) install.packages ( "rgenoud" )
# Also load DiceView for easy plot
if ( ! is.element ( "DiceView" , installed.packages ())) install.packages ( "DiceView" )
# IECI will need randtoolbox, some KrigInv routines
if ( ! is.element ( "randtoolbox" , installed.packages ())) install.packages ( "randtoolbox" )
if ( ! is.element ( "KrigInv" , installed.packages ())) install.packages ( "KrigInv" )

Objective function
fun = function ( x ){ - ( 1-1 / 2 * ( sin ( 12 * x ) / ( 1 + x ) +2 * cos ( 7 * x ) * x ^ 5+0.7 ))}
plot ( fun , type = 'l' , ylim = c ( -1.2 , 0.2 ), ylab = "y" )

X = data.frame ( X = matrix ( c ( .0 , .33 , .737 , 1 ), ncol = 1 ))
y = fun ( X )

Expected Improvement (EI)
par ( mfrow = c ( 2 , 1 ))
library ( DiceKriging )
set.seed ( 123 )
kmi <- km ( design = X , response = y , control = list ( trace = FALSE ), optim.method = 'gen' )
library ( DiceView )

`## Loading required package: DiceEval`

`## Loading required package: rgl`

```
## Warning in rgl.init(initValue, onlyNULL): RGL: unable to open X11
## display
```

`## Warning: 'rgl_init' failed, running with rgl.useNULL = TRUE`

sectionview.km ( kmi , col_surf = 'blue' , col_points = 'blue' , xlim = c ( 0 , 1 ), ylim = c ( -1.1 , 0 ), title = "" , Xname = "x" , yname = "Y(x)" )
plot ( fun , type = 'l' , add = T )
abline ( h = min ( y ), col = 'blue' , lty = 2 )
text ( labels = "m" , x = 0.05 , y = min ( y ) +0.05 , col = 'blue' )
legend ( bg = rgb ( 1 , 1 , 1 , .4 ), 0.3 , 0 , legend = c ( "(unknown) objective function" , "conditionnal gaussian process" , "observations" ), col = c ( 'black' , 'blue' , 'blue' ), lty = c ( 1 , 2 , 0 ), pch = c ( -1 , -1 , 20 ))
source ( "https://raw.githubusercontent.com/IRSN/IECI/master/EI.R" )
.xx = seq ( f = 0 , t = 1 , l = 501 )
plot ( .xx , EI ( .xx , kmi ), col = 'blue' , type = 'l' , xlab = "x" , ylab = "EI(x)" )

Expected Conditional Improvement (ECI) by Monte Carlo estimation (! slow !)
par ( mfrow = c ( 2 , 1 ))
xn = 0.4
sectionview.km ( kmi , col_surf = 'blue' , col_points = 'blue' , xlim = c ( xn -0.2 , xn +0.2 ), ylim = c ( -1.1 , -.7 ), title = "" , Xname = "x" )
plot ( fun , type = 'l' , add = T )
abline ( h = min ( y ), col = 'blue' , lty = 2 )
text ( labels = "m" , x = 0.35 , y = min ( y ) -0.05 , col = 'blue' )
abline ( v = xn , col = 'red' )
text ( labels = "xn" , x = xn +0.01 , y = -0.75 , col = 'red' )
Yn = predict.km ( kmi , newdata = xn , type = "UK" , checkNames = F )
.yy = seq ( f = Yn $ lower , t = Yn $ upper , l = 11 )
#lines(x=dnorm(.xxx,yn$mean,yn$sd)/100+xn,y=.xxx,col='red')
points ( x = rep ( xn , length ( .yy )), y = .yy , col = rgb ( 1 , 0 , 0 , 0.2+0.8 * dnorm ( .yy , Yn $ mean , Yn $ sd ) * sqrt ( 2 * pi ) * Yn $ sd ), pch = 20 )

`## Error in rgb(1, 0, 0, 0.2 + 0.8 * dnorm(.yy, Yn$mean, Yn$sd) * sqrt(2 * : alpha level 1, not in [0,1]`

points ( x = xn , y = .yy [ 8 ], col = rgb ( 0 , 1 , 0 ), pch = 1 )
points ( x = xn , y = .yy [ 4 ], col = rgb ( 0 , 1 , 0 ), pch = 1 )
plot ( .xx , EI ( .xx , kmi ), col = 'blue' , type = 'l' , xlab = "x" , ylab = "ECI(xn,x)" , xlim = c ( xn -0.2 , xn +0.2 ), ylim = c ( 0 , 0.05 ))
text ( labels = "EI(x)" , x = .4 , y = 0.03 , col = 'blue' )
for ( yn in .yy ) {
kmii <- km ( design = rbind ( X , xn ), response = rbind ( y , yn ), control = list ( trace = FALSE ))
lines ( .xx , 5 * EI ( .xx , kmii ), col = rgb ( 1 , 0 , 0 , 0.2+0.8 * dnorm ( yn , Yn $ mean , Yn $ sd ) * sqrt ( 2 * pi ) * Yn $ sd ), lty = 2 )
}

`## Error in rgb(1, 0, 0, 0.2 + 0.8 * dnorm(yn, Yn$mean, Yn$sd) * sqrt(2 * : alpha level 1, not in [0,1]`

yn = .yy [ 8 ]
kmii <- km ( design = rbind ( X , xn ), response = rbind ( y , yn ), control = list ( trace = FALSE ))
lines ( .xx , 5 * EI ( .xx , kmii ), col = rgb ( 0 , 1 , 0 ), lty = 2 )
yn = .yy [ 4 ]
kmii <- km ( design = rbind ( X , xn ), response = rbind ( y , yn ), control = list ( trace = FALSE ))
lines ( .xx , 5 * EI ( .xx , kmii ), col = rgb ( 0 , 1 , 0 ), lty = 6 )
source ( "https://raw.githubusercontent.com/IRSN/IECI/master/IECI_mc.R" , echo = FALSE )

`## Loading required package: rngWELL`

`## This is randtoolbox. For overview, type 'help("randtoolbox")'.`

lines ( .xx , Vectorize ( function ( .x ) 5 * ECI_mc_vec.o2 ( .x , xn , kmi ))( .xx ), col = 'red' , type = 'l' )

`## Error in predict_update_km(newXvalue = y0, newXmean = my0, newXvar = sdy0^2, : could not find function "predict_update_km"`

text ( labels = "ECI(xn,x) = E[ EI(x) | yn~Y(xn)] (x50)" , x = .3 , y = 0.045 , col = 'red' )

Expected Conditional Improvement (ECI) by exact calculation (! fast !)
par ( mfrow = c ( 2 , 1 ))
sectionview.km ( kmi , col_surf = 'blue' , col_points = 'blue' , xlim = c ( xn -0.2 , xn +0.2 ), ylim = c ( -1.1 , -.7 ), title = "" , Xname = "x" )
plot ( fun , type = 'l' , add = T )
abline ( h = min ( y ), col = 'blue' , lty = 2 )
text ( labels = "m" , x = 0.35 , y = min ( y ) -0.05 , col = 'blue' )
abline ( v = xn , col = 'red' )
text ( labels = "xn" , x = xn +0.01 , y = -0.75 , col = 'red' )
Yn = predict.km ( kmi , newdata = xn , type = "UK" , checkNames = F )
.yy = seq ( f = Yn $ lower , t = Yn $ upper , l = 11 )
#lines(x=dnorm(.xxx,yn$mean,yn$sd)/100+xn,y=.xxx,col='red')
plot ( .xx , EI ( .xx , kmi ), col = 'blue' , type = 'l' , xlab = "x" , ylab = "ECI(xn,x)" , xlim = c ( xn -0.2 , xn +0.2 ), ylim = c ( 0 , 0.05 ))
text ( labels = "EI(x)" , x = .4 , y = 0.03 , col = 'blue' )
source ( "https://raw.githubusercontent.com/IRSN/IECI/master/IECI.R" , echo = FALSE )

```
## Warning in file(filename, "r", encoding = encoding): cannot open file
## 'int_Phi.R': No such file or directory
```

`## Error in file(filename, "r", encoding = encoding): cannot open the connection`

lines ( .xx , 5 * ECI ( matrix ( .xx , ncol = 1 ), xn , kmi ), col = 'red' , type = 'l' )

`## Error in ECI(matrix(.xx, ncol = 1), xn, kmi): could not find function "ECI"`

text ( labels = "ECI(xn,x) = E[ EI(x) | yn~Y(xn)] (x50)" , x = .3 , y = 0.045 , col = 'red' )

And finally, Integrated Expected Conditional Improvement (IECI) calculated at the point xn:

```
> IECI_mc(xn,kmi)
[1] 0.002014966
> IECI(xn,kmi)
[1] 0.002115793
```

Overall criterions values (EI and IECI):

par ( mfrow = c ( 1 , 1 ))
sectionview.km ( kmi , ylim = c ( -1 , 0 ))
abline ( v = X $ X , lty = 2 )
n = function ( x ){( x - min ( x )) / ( max ( x ) - min ( x )) -1 }
xx = seq ( f = 0 , t = 1 , l = 1000 )
lines ( xx ,( fun ( xx )), type = 'l' )
lines ( xx , n ( IECI ( x 0 = xx , model = kmi , lower = 0 , upper = 1 )), type = 'l' , col = 'blue' )

`## Error in IECI(x0 = xx, model = kmi, lower = 0, upper = 1): could not find function "IECI"`

text ( labels = "IECI(x)" , x = 0.6 , y = 0.9-1 , col = 'blue' )
lines ( xx , n ( EI ( xx , model = kmi )), col = 'red' )
text ( labels = "EI(x)" , x = 0.6 , y = 0.1-1 , col = 'red' )

Analytical development of IECI:
ECI is viewed in two main parts : updated realizations where the new conditional point $y_n$ is below current minimum $m$, and updated realizations where the new conditional point $y_n$ is over current minimum $m$ (which then does not change this current minimum for EI calculation):

Change variable: $z := \frac {y_n-\mu(x_n)}{\sigma(x_n)}$

$y_n=\sigma(x_n) z + \mu(x_n)$
$dz = dy_n/\sigma(x_n)$
$y_n \in [m,+\infty[ \sim z \in [m’ = (m-\mu(x_n))/\sigma(x_n), + \infty[$
So,

Considering that ([@KrigingUpdate]),

Finally,

Given,

Change variable: $z := \frac {y_n-\mu(x_n)}{\sigma(x_n)}$

$y_n=\sigma(x_n) z + \mu(x_n)$
$dz = dy_n/\sigma(x_n)$
$y_n \in [m,+\infty[ \sim z \in [m’ = (m-\mu(x_n))/\sigma(x_n), + \infty[$
So,

Finally,

Given,

These integrals need the following expressions to become tractable:

$\int_{x}^{y} (a+b z) \Phi(a+b z) e^{-\frac{z^2}{2}} dz$
$\int_{x}^{y} \phi(a+b z) e^{-\frac{z^2}{2}} dz$
We have (see below):

And

It should be noticed that this exact calculation of ECI is vectorizable, which means that synchronized computations may be performed in an efficient manner.

Detailed calculation of integrals in $ECI^y_n$ and $ECI^m_n$
$\int_{x}^{y} \phi(a+b z) e^{-\frac{z^2}{2}} dz$

$\int_{x}^{y} (a+b z) \Phi(a+b z) e^{-\frac{z^2}{2}} dz$
Changing variable: $t:=a+b z$,

We use the bi-normal cumulative density approximation [@Genz1992]:

Changing variable $t:=v-(a+b u)$,

and thus,