add probability of improvement

NAMESPACE

@@ -22,7 +22,9 @@ export(crit.dib1)
 export(crit.ei)
 export(crit.eqi)
 export(crit.mr)
+export(crit.poi)
 export(crit.se)
+export(crit.sei)
 export(exampleRun)
 export(exampleRunMultiObj)
 export(finalizeSMBO)
@@ -46,6 +48,8 @@ export(makeMBOInfillCritDIB)
 export(makeMBOInfillCritEI)
 export(makeMBOInfillCritEQI)
 export(makeMBOInfillCritMeanResponse)
+export(makeMBOInfillCritPOI)
+export(makeMBOInfillCritSEI)
 export(makeMBOInfillCritStandardError)
 export(makeMBOLearner)
 export(makeMBOTrafoFunction)

R/infill_crits.R

@@ -81,7 +81,6 @@ makeMBOInfillCritStandardError = function() {
 #' @rdname infillcrits
 makeMBOInfillCritEI = function(se.threshold = 1e-6) {
   assertNumber(se.threshold, lower = 1e-20)
-  force(se.threshold)
   makeMBOInfillCrit(
     fun = function(points, models, control, par.set, designs, iter, progress, attributes = FALSE) {
       model = models[[1L]]
@@ -93,12 +92,9 @@ makeMBOInfillCritEI = function(se.threshold = 1e-6) {
       p = predict(model, newdata = points)$data
       p.mu = maximize.mult * p$response
       p.se = p$se
-      y.min = min(y)
-      d = y.min - p.mu
+      d = min(y) - p.mu
       xcr = d / p.se
-      xcr.prob = pnorm(xcr)
-      xcr.dens = dnorm(xcr)
-      ei = d * xcr.prob + p.se * xcr.dens
+      ei = d * pnorm(xcr) + p.se * dnorm(xcr)
       res = ifelse(p.se < se.threshold, 0, -ei)
       if (attributes) {
         res = setAttribute(res, "crit.components", data.frame(se = p$se, mean = p$response))
@@ -114,6 +110,78 @@ makeMBOInfillCritEI = function(se.threshold = 1e-6) {
   )
 }
 
+#' @export
+#' @rdname infillcrits
+makeMBOInfillCritSEI = function(se.threshold = 1e-6) {
+  # On a New Improvement-Based Acquisition Function for Bayesian Optimization
+  # Noè et.al. 2018
+  # http://arxiv.org/abs/1808.06918
+
+  assertNumber(se.threshold, lower = 1e-20)
+  makeMBOInfillCrit(
+    fun = function(points, models, control, par.set, designs, iter, progress, attributes = FALSE) {
+      model = models[[1L]]
+      design = designs[[1]]
+      maximize.mult = if (control$minimize) 1 else -1
+      assertString(control$y.name)
+      y = maximize.mult * design[, control$y.name]
+      assertNumeric(y, any.missing = FALSE)
+      p = predict(model, newdata = points)$data
+      p.mu = maximize.mult * p$response
+      p.se = p$se #s(x)
+      d = min(y) - p.mu
+      xcr = d / p.se #u(x) in paper
+      p.xcr = pnorm(xcr)
+      d.xcr = dnorm(xcr)
+      ei = d * p.xcr + p.se * d.xcr
+      vi = p.se^2 * ((xcr^2 + 1) * p.xcr + xcr * d.xcr) - ei^2
+      vi = pmax(0, vi)
+      sei = ei / sqrt(vi)
+      res = ifelse(p.se < se.threshold | is.infinite(sei) | is.nan(sei), 0, -sei)
+      if (attributes) {
+        res = setAttribute(res, "crit.components", data.frame(se = p$se, mean = p$response))
+      }
+      return(res)
+    },
+    name = "Scaled expected improvement",
+    id = "sei",
+    components = c("se", "mean"),
+    params = list(se.threshold = se.threshold),
+    opt.direction = "maximize",
+    requires.se = TRUE
+  )
+}
+
+#' @export
+#' @rdname infillcrits
+makeMBOInfillCritPOI = function(se.threshold = 1e-6) {
+  # https://www.cse.wustl.edu/~garnett/cse515t/spring_2015/files/lecture_notes/12.pdf
+  assertNumber(se.threshold, lower = 1e-20)
+  makeMBOInfillCrit(
+    fun = function(points, models, control, par.set, designs, iter, progress, attributes = FALSE) {
+      model = models[[1L]]
+      design = designs[[1]]
+      maximize.mult = if (control$minimize) 1 else -1
+      assertString(control$y.name)
+      y = maximize.mult * design[, control$y.name]
+      assertNumeric(y, any.missing = FALSE)
+      p = predict(model, newdata = points)$data
+      poi = pnorm(min(y) , mean = maximize.mult * p$response, sd = p$se)
+      res = ifelse(p$se < se.threshold, 0, -poi)
+      if (attributes) {
+        res = setAttribute(res, "crit.components", data.frame(se = p$se, mean = p$response))
+      }
+      return(res)
+    },
+    name = "Probability of improvement",
+    id = "poi",
+    components = c("se", "mean"),
+    params = list(se.threshold = se.threshold),
+    opt.direction = "maximize",
+    requires.se = TRUE
+  )
+}
+
 #' @export
 #' @rdname infillcrits
 makeMBOInfillCritCB = function(cb.lambda = NULL) {

R/zzz.R

@@ -96,3 +96,17 @@ crit.eqi = makeMBOInfillCritEQI()
 #' @format NULL
 #' @keywords NULL
 crit.dib1 = makeMBOInfillCritDIB(cb.lambda = 1)
+#' @rdname MBOInfillCrit
+#' @export
+#' @usage NULL
+#' @docType NULL
+#' @format NULL
+#' @keywords NULL
+crit.poi = makeMBOInfillCritPOI()
+#' @rdname MBOInfillCrit
+#' @export
+#' @usage NULL
+#' @docType NULL
+#' @format NULL
+#' @keywords NULL
+crit.sei = makeMBOInfillCritSEI()

tests/testthat/test_infillcrits.R

@@ -47,10 +47,18 @@ test_that("infill crits", {
       expect_true(!anyMissing(opdf$infill_cb[, c("se","mean","lambda")]))
       expect_true(all(opdf$infill_cb$lambda == 2))
     }
-    if (!is.null(opdf$infill_aei))
+    if (!is.null(opdf$infill_aei)) {
       expect_true(!anyMissing(opdf$infill_aei[, c("se","mean","tau")]))
-    if (!is.null(opdf$infill_eqi))
+    }
+    if (!is.null(opdf$infill_eqi)) {
       expect_true(!anyMissing(opdf$infill_eqi[, c("se","mean","tau")]))
+    }
+    if (!is.null(opdf$infill_poi)) {
+      expect_true(!anyMissing(opdf$infill_poi[, c("poi","se","mean")]))
+    }
+    if (!is.null(opdf$infill_sei)) {
+      expect_true(!anyMissing(opdf$infill_poi[, c("sei","se","mean")]))
+    }
     expect_true(nrow(opdf$final_eval) == 10L)
   }
 
@@ -65,7 +73,7 @@ test_that("infill crits", {
   # we have converged and just waste time. we need to detect this somehow, or cope with it
   for (noisy in c(TRUE, FALSE)) {
     for (minimize in c(TRUE, FALSE)) {
-      crits = if (noisy) list(crit.aei, crit.eqi) else list(crit.mr, crit.se, crit.ei, crit.cb2)
+      crits = if (noisy) list(crit.aei, crit.eqi) else list(crit.mr, crit.se, crit.ei, crit.cb2, crit.poi, crit.sei)
       for (lrn in learners) {
         if (inherits(lrn, "regr.km"))
           lrn = setHyperPars(lrn, nugget.estim = noisy)