diff --git a/NAMESPACE b/NAMESPACE
index 390bc12d5..df1154ad4 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -58,6 +58,7 @@ export(renderExampleRunPlot)
 export(setMBOControlInfill)
 export(setMBOControlMultiObj)
 export(setMBOControlMultiPoint)
+export(setMBOControlNoisy)
 export(setMBOControlTermination)
 export(trafoLog)
 export(trafoSqrt)
diff --git a/R/OptState.R b/R/OptState.R
index d0f676b49..11a4cc3ee 100644
--- a/R/OptState.R
+++ b/R/OptState.R
@@ -50,7 +50,8 @@ NULL
 
 makeOptState = function(opt.problem, loop = 0L, tasks = NULL, models = NULL,
   time.model = NULL, opt.result = NULL, state = "init", opt.path = NULL,
-  time.last.saved = Sys.time(), loop.starttime = Sys.time(), time.used = 0L, progress = 0, time.created = Sys.time()) {
+  time.last.saved = Sys.time(), loop.starttime = Sys.time(), time.used = 0L, progress = 0, time.created = Sys.time(),
+  identification.time.used = 0L) {
 
   opt.state = new.env()
 
@@ -68,6 +69,7 @@ makeOptState = function(opt.problem, loop = 0L, tasks = NULL, models = NULL,
   opt.state$loop.starttime = loop.starttime
   opt.state$time.used = time.used
   opt.state$progress = progress
+  opt.state$identification.time.used = identification.time.used
 
   opt.state$random.seed = getRandomSeed()
   opt.state$time.created = time.created
diff --git a/R/OptState_getter.R b/R/OptState_getter.R
index 6fa09db18..2395eefe3 100644
--- a/R/OptState_getter.R
+++ b/R/OptState_getter.R
@@ -144,3 +144,38 @@ getOptStateValidTerminationStates = function() {
   c("term.iter", "term.time", "term.exectime", "term.yval", "term.feval", "term.custom")
 }
 
+getOptStateIntensification = function(opt.state) {
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+  return(!is.null(control$noisy.method) && control$noisy.method %in% c("incumbent", "ocba"))
+}
+
+getOptStateIdentification = function(opt.state) {
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+  return(control$identification.time.budget > 0 | control$identification.max.evals > 0)
+}
+
+getOptStatePCS = function(opt.state) {
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+  return(control$noisy.identification.pcs)
+}
+
+getOptStateTerminationIdentification = function(opt.state) {
+  terminate = shouldTerminateIdentification.OptState(opt.state)
+  setOptStateProgress(opt.state, terminate$progress)
+  # update only if termination condition is met
+  if (terminate$term) {
+    setOptStateState(opt.state, terminate$code)
+  }
+  terminate
+}
+
+getOptStateStartTimeIdentification = function(opt.state) {
+  opt.state$identification.starttime
+}
+
+getOptStateTimeUsedIdentification = function(opt.state) {
+  opt.state$identification.time.used
+}
\ No newline at end of file
diff --git a/R/OptState_setter.R b/R/OptState_setter.R
index d417ff96f..b14941fca 100644
--- a/R/OptState_setter.R
+++ b/R/OptState_setter.R
@@ -41,6 +41,11 @@ setOptStateLoopStarttime = function(opt.state) {
   invisible()
 }
 
+setOptStateIdentificationStarttime = function(opt.state) {
+  opt.state$identification.starttime = Sys.time()
+  invisible()
+}
+
 setOptStateTimeUsed = function(opt.state, time.used = NULL, time.add = NULL) {
   if (!is.null(time.used)) {
     opt.state$time.used = time.used
@@ -64,3 +69,15 @@ setOptStateProgress = function(opt.state, progress) {
   opt.state$progress = progress
   invisible()
 }
+
+setOptStateTimeUsedIdentification = function(opt.state, time.used = NULL, time.add = NULL) {
+  if (!is.null(time.used)) {
+    opt.state$identification.time.used = time.used
+  } else if (!is.null(time.add)) {
+    opt.state$identification.time.used = getOptStateTimeUsedIdentification(opt.state) + time.add
+  } else {
+    opt.state$identification.time.used = getOptStateTimeUsedIdentification(opt.state) + difftime(Sys.time(), getOptStateStartTimeIdentification(opt.state), units = "secs")
+    setOptStateIdentificationStarttime(opt.state)
+  }
+  invisible()
+}
diff --git a/R/SMBO.R b/R/SMBO.R
index e8989adc2..ab86eb2a7 100644
--- a/R/SMBO.R
+++ b/R/SMBO.R
@@ -49,6 +49,8 @@ initSMBO = function(par.set, design, learner = NULL, control, minimize = rep(TRU
 #'   Outcome of the optimization.
 #'   For multiple results use a list.
 #'   For a result of a multi-objective function use a numeric vector.
+#'   For multiple results of for noisy instances use a list.
+#'   Each list element should correspond to one x value.
 #'
 #' @return [\code{\link{OptState}}]
 #' @export
diff --git a/R/evalTargetFun.R b/R/evalTargetFun.R
index 3f79fb118..125fe33ee 100644
--- a/R/evalTargetFun.R
+++ b/R/evalTargetFun.R
@@ -25,13 +25,30 @@ evalTargetFun.OptState = function(opt.state, xs, extras) {
   # short names and so on
   nevals = length(xs)
   ny = control$n.objectives
+
+  # trafo X points
+  xs.trafo = lapply(xs, trafoValue, par = par.set)
+
+  # handle noisy instances
+  if (isTRUE(control$noisy.instances > 1L)) {
+    nevals = nevals * control$noisy.instances
+    xs = rep(xs, each = control$noisy.instances)
+    extras = rep(extras, each = control$noisy.instances)
+    if (!control$noisy.self.replicating) {
+      xs.trafo = rep(xs.trafo, each = control$noisy.instances)
+      if (!is.na(control$noisy.instance.param)) {
+        inst.param = lapply(seq_len(control$noisy.instances), function(x) setNames(list(x), control$noisy.instance.param))
+        xs.trafo = Map(c, xs.trafo, inst.param)
+      }
+    }
+  }
+
+
   num.format = control$output.num.format
   num.format.string = paste("%s = ", num.format, sep = "")
   dobs = ensureVector(asInteger(getOptStateLoop(opt.state)), n = nevals, cl = "integer")
   imputeY = control$impute.y.fun
 
-  # trafo X points
-  xs.trafo = lapply(xs, trafoValue, par = par.set)
 
   # function to measure of fun call
     wrapFun = function(x) {
@@ -43,6 +60,9 @@ evalTargetFun.OptState = function(opt.state, xs, extras) {
         user.extras = attr(y, "extras")
         y = setAttribute(y, "extras", NULL)
       }
+      if (!is.null(control$noisy.instance.param) && !is.na(control$noisy.instance.param) && !control$noisy.self.replicating) {
+        user.extras = c(user.extras, x[control$noisy.instance.param])
+      }
       st = proc.time() - st
       list(y = y, time = st[3], user.extras = user.extras)
     }
@@ -56,6 +76,21 @@ evalTargetFun.OptState = function(opt.state, xs, extras) {
   res = parallelMap(wrapFun, xs.trafo, level = "mlrMBO.feval",
     impute.error = if (is.null(imputeY)) NULL else identity)
 
+  # handle noisy instances of self.replicating functions
+  if (isTRUE(control$noisy.instances > 1L) && control$noisy.self.replicating) {
+    xs.trafo = rep(xs.trafo, each = control$noisy.instances)
+    res = lapply(res, function(r) {
+      if (is.error(r)) {
+        rep(list(r), control$noisy.instances)
+      } else {
+        lapply(seq_along(r$y), function(i) {
+          list(y = r$y[i], time = r$time / length(r$y), user.extras = c(r$user.extras, setNames(list(i), control$noisy.instance.param)))
+        })
+      }
+    })
+    res = unlist(res, recursive = FALSE)
+  }
+
   # loop evals and to some post-processing
   for (i in seq_len(nevals)) {
     r = res[[i]]; x = xs[[i]]; x.trafo = xs.trafo[[i]]; dob = dobs[i]
diff --git a/R/identifyFinalPoints.R b/R/identifyFinalPoints.R
new file mode 100644
index 000000000..572a7f1c8
--- /dev/null
+++ b/R/identifyFinalPoints.R
@@ -0,0 +1,68 @@
+identifyFinalPoints = function(opt.state, min.pcs = NULL, time.budget = NULL) {
+
+  # some initialization
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+  par.set = getOptProblemParSet(opt.problem)
+  op = as.data.table(getOptStateOptPath(opt.state))
+  par.names = colnames(op)[1:(which(colnames(op) == "y") - 1)] #FIXME: This sucks
+  min.pcs = min.pcs %??% getOptStatePCS(opt.state)
+
+  # calculate summary of the design
+  ds = getOptPathSummary(opt.state, par.names)
+  nds = nrow(ds)
+
+  # set start time for identification here
+
+  # make sure that initially, each point is evaluated at least minrep times
+  xinit = rep(seq_len(nds), pmax(2 - ds$runs, 0))
+  opt.state = replicatePoint(opt.state, x = ds[xinit, ..par.names], type = paste("initeval"))
+
+  setOptStateTimeUsedIdentification(opt.state)
+  terminate = getOptStateTerminationIdentification(opt.state)
+
+  pcs = calculatePCS(opt.state)
+  while(pcs < min.pcs) {
+      ds = getOptPathSummary(opt.state, par.names)
+      add = distributeOCBA(ds, budget = 3)
+      reps = rep(seq_len(nrow(ds)), add)
+      replicatePoint(opt.state, x = ds[reps, ..par.names], type = paste("identification"))
+      setOptStateTimeUsedIdentification(opt.state)
+      terminate = getOptStateTerminationIdentification(opt.state)
+      showInfo(getOptProblemShowInfo(opt.problem), "[mbo] identification: P(CS) %.3f / %.3f", pcs, min.pcs)
+      pcs = calculatePCS(opt.state)
+
+      if(terminate$term)
+        break
+  }
+
+  return(opt.state)
+}
+
+
+calculatePCS = function(opt.state) {
+  # calculate the (approximate) probability of correct selection
+  # best observed design
+  op = as.data.table(getOptStateOptPath(opt.state))
+  par.names = colnames(op)[1:(which(colnames(op) == "y") - 1)] #FIXME: This sucks
+  ds = getOptPathSummary(opt.state, par.names)
+
+  ds = ds[order(ds$y), ]
+  b = 1
+
+  # mean anc covariance for vector (y1 - yb, y2 - yb, ..., yn - yb)
+  vf = rep(1, nrow(ds) - 1)
+  trafo.mat = cbind(vf, diag(- vf)) # transformation matrix A
+
+  m = ds[b, ]$y - ds[- b, ]$y
+
+  # covariance for vector (y1 - yb, y2 - yb, ..., yn - yb)
+  # is obtained by multiplication with matrix A = [(1 -1, 0), (1, 0, -1)]
+  # rule is cov(A * Y) = A * cov(Y) * t(A)
+  sigma = trafo.mat %*% diag(ds$ysd^2 / ds$runs) %*% t(trafo.mat)
+ #  sigma = diag(sqrt(ds$ysd[-b]^2 / ds$runs[-b]^2 + ds$ysd[b]^2 / ds$runs[b]^2))
+
+  pcs = pmvnorm(upper = 0, mean = m, sigma = sigma)
+
+  return(pcs[1])
+}
diff --git a/R/intensifyOptState.R b/R/intensifyOptState.R
new file mode 100644
index 000000000..3afd6916f
--- /dev/null
+++ b/R/intensifyOptState.R
@@ -0,0 +1,178 @@
+intensifyOptState = function(opt.state) {
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+  
+  switch(control$noisy.method, 
+      "incumbent" = intensifyIncumbent(opt.state),
+      "ocba" = intensifyOCBA(opt.state)
+    )
+}
+
+intensifyIncumbent = function(opt.state) {
+
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+  op = as.data.table(getOptStateOptPath(opt.state))
+  par.names = colnames(op)[1:(which(colnames(op) == "y") - 1)] #FIXME: This sucks
+  
+  # get a summary of the design  
+  ds = getOptPathSummary(opt.state, par.names)
+  nds = nrow(ds)
+
+  # incumbent: current best point w. r. t. mean over all function evaluations
+  # the newest point cannot be the incumbent, it is always a challenger
+  # DOES NOT WORK FOR MULTIPOINT PROPOSAL YET
+  inc = which.min(ds[- nds, ]$y) 
+  # incumbent is replicated once in each iteration
+  replicatePoint(opt.state, x = ds[inc, ..par.names], type = "incumbent", reps = 1L)
+
+  # determine a set of challengers
+  if (control$noisy.incumbent.nchallengers == 0L) {
+    cls = c(nds)
+  } else {
+    # determine set of p points to be challenged against incumbent
+    # incumbent is excluded (cannot be challenged against itself)
+    # and new point is always set as a challenger
+    # points are drawn randomly without replacement with probability prop. to their function value
+    cls = setdiff(seq_len(nds), c(inc, nds))
+    p = min(control$noisy.incumbent.nchallengers, nds - 2)
+    probs = exp(- ds[cls, ]$y) / sum(exp(- ds[cls, ]$y))
+    cls = sample(cls, size = p, prob = probs, replace = FALSE)
+    cls = c(cls, nds)
+  }
+
+  # start the race
+  for (cl in cls) {
+
+    r = 1L
+    replicatePoint(opt.state, x = ds[cl, ..par.names], type = paste("challenger"), reps = r)
+    ds = getOptPathSummary(opt.state, par.names)
+
+    # proceed as long as challenger has less runs than incumbent and is better than incumbent
+    while((ds[cl, "runs"] < ds[inc, "runs"]) && (ds[cl, "y"] < ds[inc, "y"])) {
+      r = 2L * r
+      replicatePoint(opt.state, x = ds[cl, ..par.names], type = paste("challenger"), reps = r)
+      ds = getOptPathSummary(opt.state, par.names)
+    }
+
+  }
+  return(opt.state)
+}
+
+replicatePoint = function(opt.state, x, type, reps = 1L) {
+
+  # replicate rows according to the number of desired replicates
+  xs = seq_len(nrow(x))
+  xrep = x[rep(xs, reps), ]
+   
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+
+  prop = makeProposal(control, xrep, prop.type = rep(type, nrow(xrep)))
+  evalProposedPoints.OptState(opt.state, prop)
+
+  return(opt.state)
+}
+
+getOptPathSummary = function(opt.state, par.names) {
+    op = as.data.table(getOptStateOptPath(opt.state))
+    ds = op[, .(y = mean(y), ysd = sd(y), runs = .N), by = par.names]
+    return(ds)
+}
+
+
+intensifyOCBA = function(opt.state) {
+  
+  # some intialization
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+  par.set = getOptProblemParSet(opt.problem)
+
+  op = as.data.table(getOptStateOptPath(opt.state))
+  par.names = colnames(op)[1:(which(colnames(op) == "y") - 1)] #FIXME: This sucks
+  
+  # minimum number of replicates at each point 
+  minrep = max(control$noisy.ocba.initial, 2L)
+
+  # calculate summary of the dsign 
+  ds = getOptPathSummary(opt.state, par.names)
+  nds = nrow(ds)
+
+  # make sure that initially, each point is evaluated at least minrep times 
+  xinit = rep(seq_len(nds), pmax(minrep - ds$runs, 0))
+  opt.state = replicatePoint(opt.state, x = ds[xinit, ..par.names], type = paste("initeval"))
+
+  ds = getOptPathSummary(opt.state, par.names)
+  add = distributeOCBA(ds, budget = control$noisy.ocba.budget)
+  reps = rep(seq_len(nds), add)
+  
+  replicatePoint(opt.state, x = ds[reps, ..par.names], type = paste("OCBA"))
+
+  return(opt.state)
+}
+
+
+distributeOCBA = function(ds, budget = 3L) {
+
+  nds = nrow(ds)
+
+  # TODO: until now only minimization possible 
+  tbudget = budget + sum(ds$runs)
+  
+  # search for the best and second-best design
+  b = order(ds$y)[1]
+  or = order(ds$y)[2:nds]
+  s = or[ds[or, ]$ysd > 0][1]
+
+  
+  # vector of ratios
+  ratio = rep(0, nds)
+  ratio[s] = 1
+  
+  # calculate ratios
+  tmp = (ds[b, ]$y - ds[s, ]$y) / (ds[b, ]$y - ds[- c(s, b), ]$y)
+  tmp[!is.finite(tmp)] = min(tmp[is.finite(tmp)], 10e-10, na.rm = TRUE)
+  ratio[- c(s, b)] = tmp^2 * ds[- c(s, b), ]$ysd^2 / ds[s, ]$ysd^2
+
+  # numerically ysd is sometimes 0 (in fact, this has prob. 0 in the noisy case)
+  # these indices need to be excluded
+  ysdval = which(ds$ysd > 0)
+  ratio[b] = ds[b, ]$ysd * sqrt(sum(ratio[ysdval]^2 / ds$ysd[ysdval]^2))
+  
+  # additional replications
+  add = rep(0, nds)
+  
+  # disable designs where ratio is zero
+  disabled = (ratio == 0)
+  
+  more_alloc = TRUE
+  
+  while (more_alloc) {
+    
+    add[!disabled] = roundPreserveSum(tbudget / sum(ratio[!disabled]) * ratio[!disabled])
+    
+    # disable designs that have been run too much
+    disabled = disabled | (ds$runs > add)
+    more_alloc = any(ds$runs > add)
+    
+    # set additional replications s.t. already run replications are set
+    add[disabled] = ds[disabled, ]$runs
+    
+    # decrease total budget correspondingly
+    tbudget = budget + sum(ds$runs) - sum(add[disabled])
+  }
+  
+  add = add - ds$runs
+  
+  return(add)
+}
+
+roundPreserveSum = function(x, digits = 0) {
+  x[!is.finite(x)] = 0 # happens because of 0 / 0 and should be 0
+  up = 10 ^ digits
+  x = x * up
+  y = floor(x)
+  indices = tail(order(x - y), round(sum(x)) - sum(y))
+  y[indices] = y[indices] + 1
+  y / up
+}
\ No newline at end of file
diff --git a/R/makeTaskSingleObj.R b/R/makeTaskSingleObj.R
index 6f7aca2cc..4416fb89a 100644
--- a/R/makeTaskSingleObj.R
+++ b/R/makeTaskSingleObj.R
@@ -21,5 +21,13 @@ makeTaskSingleObj = function(opt.path, control) {
     data[[y.name]] = trafo.y.fun(data[[y.name]])
   }
 
+  agg = control$noisy.instance.aggregation
+
+  if (!is.null(agg)) {
+    par.names = colnames(data)[1:(which(colnames(data) == "y") - 1)] #FIXME: This sucks
+    data = setDT(data)[, .(y = agg(y)), by = par.names]
+    data = data.frame(data)
+  }
+
   makeRegrTask(target = control$y.name, data = data)
 }
diff --git a/R/mboTemplate.R b/R/mboTemplate.R
index dfb0a27c4..7c734e386 100644
--- a/R/mboTemplate.R
+++ b/R/mboTemplate.R
@@ -34,15 +34,28 @@ mboTemplate.OptState = function(obj) {
       was satisfied right after the creation of the initial design!", terminate$message)
     return(opt.state)
   }
+
+  # optimization phase
   repeat {
     prop = proposePoints(opt.state)
     evalProposedPoints.OptState(opt.state, prop)
+    intensify = getOptStateIntensification(opt.state)
+    if (intensify) {
+      intensifyOptState(opt.state)
+    }
     finalizeMboLoop(opt.state)
     terminate = getOptStateTermination(opt.state)
     if (terminate$term) {
       break
     }
   }
+
+  # identification phase
+  identify = getOptStateIdentification(opt.state)
+  if (identify) {
+    setOptStateIdentificationStarttime(opt.state)
+    identifyFinalPoints(opt.state) 
+  }
   opt.state
 }
 
diff --git a/R/setMBOControlNoisy.R b/R/setMBOControlNoisy.R
new file mode 100644
index 000000000..e5ae5a28b
--- /dev/null
+++ b/R/setMBOControlNoisy.R
@@ -0,0 +1,63 @@
+#' @title Set options for handling noisy functions.
+#' @description
+#' Extends an MBO control object with options for handling noisy functions. 
+#' @template arg_control
+#' @param method [\code{character(1)}]\cr
+#'   Which of the replication strategies should be used? Possible values are:
+#'   \dQuote{fixed}: Every point is evaluated \code{instances} times. \cr
+#'   \dQuote{incumbent}: Use an incumbent strategy as intensification strategy. 
+#'    The size of the set of additional challengers (apart from the incumbent) can be specified in \code{incumbent.nchallengers}. \cr
+#'   \dQuote{ocba}: Distribution replication budget according to OCBA. 
+#'    The replication budget per iteration is specified in \code{ocba.budget},
+#'    the initial number of iterations per parameter is specified in \code{ocba.initial}. \cr 
+#' @param instances [\code{integer(1)}]\cr
+#'   How many instances of one parameter will be calculated?
+#' @param instance.param [\code{character(1)}]\cr
+#'   What is the name of the function param that defines the instance?
+#' @param self.replicating [\code{logical(1)}]\cr
+#'   TRUE if the function returns a vector of noisy results for one input. Then \code{instances} specifies the length of the result we expect.
+#' @param incumbent.nchallengers [\code{integer(1)}]\cr
+#'   The size of the set of additional challengers (apart from the incumbent), defaults to \code{0}.
+#' @param ocba.budget [\code{integer(1)}]\cr
+#'   The budget that is allocated in each iteration per Optimal Computing Budget Allocation (OCBA) rule, defaults to 10.
+#' @param ocba.initial [\code{integer(1)}]\cr
+#'   The number of initial replications at each new point, defaults to \code{3}. 
+#'   This needs to be larger than 1, since OCBA requires an initial variance estimate at each point.
+#' @param instance.aggregation [\code{function}]\cr
+#'   Should data be aggregated per instance? If yes, a function (e. g. mean) needs to be specified.  
+#' @param identification.pcs [\code{numeric(1)}]\cr
+#'   Minimum probability of correct selection. 
+#'   If an identification budget is set, the algorithm spends this to achieve a minimum probability of correct selection of the final point.
+#'   Defaults to \code{0.7}.
+#' @return [\code{\link{MBOControl}}].
+#' @family MBOControl
+#' @export
+setMBOControlNoisy = function(control,
+  method = NULL,
+  instances = NULL,
+  instance.param = NULL,
+  instance.aggregation = NULL,
+  self.replicating = NULL,
+  incumbent.nchallengers = NULL, 
+  ocba.budget = NULL,
+  ocba.initial = NULL,
+  identification.pcs = NULL) {
+
+  assertClass(control, "MBOControl")
+  control$noisy.method = coalesce(method, control$noisy.method, "fixed")
+  assertChoice(control$noisy.method, choices = c("fixed", "incumbent", "ocba"))
+  control$noisy.instances = assertInt(instances, lower = 1L, null.ok = TRUE, na.ok = FALSE) %??% control$noisy.instances %??% 1L
+  control$noisy.self.replicating = assertFlag(self.replicating, null.ok = TRUE, na.ok = FALSE) %??% control$noisy.self.replicating %??% FALSE
+  control$noisy.instance.param = assertString(instance.param, null.ok = TRUE, na.ok = TRUE) %??% control$noisy.instance.param %??% ifelse(control$noisy.self.replicating, "noisy.repl", NA_character_)
+  control$noisy.instance.aggregation = assertClass(instance.aggregation, "function", null.ok = TRUE) %??% control$noisy.instance.aggregation 
+  control$noisy.ocba.budget = assertInt(ocba.budget, lower = 0L, null.ok = TRUE, na.ok = FALSE) %??% control$noisy.ocba.budget %??% 3L
+  control$noisy.ocba.initial = assertInt(ocba.initial, lower = 2L, null.ok = TRUE, na.ok = FALSE) %??% control$noisy.ocba.initial %??% 3L
+  control$noisy.incumbent.nchallengers = assertInt(incumbent.nchallengers, lower = 0L, null.ok = TRUE, na.ok = FALSE) %??% control$noisy.incumbent.nchallengers %??% 0L
+  control$noisy.identification.pcs = assertNumeric(identification.pcs, lower = 0, upper = 1, null.ok = TRUE) %??% control$noisy.identification.pcs %??% 0.7
+
+  if (control$noisy.self.replicating && control$noisy.instance.param != "noisy.repl") {
+    stop("You can not change the instance.param for self replicating functions.")
+  }
+
+  return(control)
+}
diff --git a/R/setMBOControlTermination.R b/R/setMBOControlTermination.R
index 7c587c991..c1abd3c71 100644
--- a/R/setMBOControlTermination.R
+++ b/R/setMBOControlTermination.R
@@ -36,6 +36,14 @@
 #'   The default is \code{NULL} which means, that the first supplied argument is taken, following the order of the function signature.
 #'   Other values can be \code{"iters"}, \code{"time.budget"}, etc.\cr
 #'   If you want to to use it together with a criterion you supplied in \code{more.termination.conds}, \code{more.termination.conds} has to be a named list and the function further has to return a list element \code{progress} with values between 0 and 1.
+#' @param identification.time.budget [\code{integer(1)} | NULL] \cr
+#'   Time budget for identification in seconds. This budget is put on top to the runtime budget and 
+#'   used for final identification of the best point in the noisy case. Note that this doesn't make sense if the function is not noisy.
+#'   The default \code{NULL} means: there is no budget spent for final identification.
+#' @param identification.max.evals [\code{integer(1)} | NULL] \cr
+#'   Maximum number of evaluations performed for identification. This budget is put on top to the budget spent for optimization and 
+#'   used for final identification of the best point in the noisy case. Note that this doesn't make sense if the function is not noisy.
+#'   The default \code{NULL} means: there are no evaluations spent for final identification.
 #' @return [\code{\link{MBOControl}}].
 #' @family MBOControl
 #' @export
@@ -61,13 +69,15 @@
 #' res = mbo(fn, control = ctrl)
 #' print(res)
 setMBOControlTermination = function(control,
-  iters = NULL, time.budget = NULL, exec.time.budget = NULL, target.fun.value = NULL, max.evals = NULL, more.termination.conds = list(), use.for.adaptive.infill = NULL) {
+  iters = NULL, time.budget = NULL, exec.time.budget = NULL, target.fun.value = NULL, max.evals = NULL, more.termination.conds = list(), use.for.adaptive.infill = NULL,
+  identification.time.budget = NULL, identification.max.evals = NULL) {
 
   assertClass(control, "MBOControl")
   assertList(more.termination.conds)
   assertCharacter(use.for.adaptive.infill, null.ok = TRUE)
 
   stop.conds = more.termination.conds
+  stop.conds.identification = list()
 
   if (is.null(iters) && is.null(time.budget) && is.null(exec.time.budget) && is.null(max.evals) && length(stop.conds) == 0L) {
     stopf("You need to specify a maximal number of iteration, a time budget or at least
@@ -96,6 +106,13 @@ setMBOControlTermination = function(control,
     stop.conds = c(stop.conds, max.evals = makeMBOTerminationMaxEvals(max.evals))
   }
 
+  if (!is.null(identification.time.budget)) {
+    stop.conds.identification = c(stop.conds.identification, time.budget = makeMBOTerminationIdentificationMaxBudget(identification.time.budget))
+  }
+
+  if (!is.null(identification.max.evals)) {
+    stop.conds.identification = c(stop.conds.identification, max.evals = makeMBOTerminationIdentificationMaxEvals(identification.max.evals))
+  }
 
   # sanity check termination conditions
   lapply(stop.conds, function(stop.on) {
@@ -110,6 +127,7 @@ setMBOControlTermination = function(control,
   }
 
   control$stop.conds = stop.conds
+  control$stop.conds.identification = stop.conds.identification
 
   # store stuff in control object since it is needed internally
   control$iters = coalesce(iters, control$iters, Inf)
@@ -117,6 +135,8 @@ setMBOControlTermination = function(control,
   control$exec.time.budget = coalesce(exec.time.budget, control$exec.time.budget, Inf)
   control$max.evals = coalesce(max.evals, Inf)
   control$use.for.adaptive.infill = use.for.adaptive.infill
+  control$identification.time.budget = coalesce(identification.time.budget, control$identification.time.budget, 0L)
+  control$identification.max.evals = coalesce(identification.max.evals, control$identification.max.evals, 0L)
 
   return(control)
 }
diff --git a/R/shouldTerminate.R b/R/shouldTerminate.R
index f9281d556..5b8c4c217 100644
--- a/R/shouldTerminate.R
+++ b/R/shouldTerminate.R
@@ -25,3 +25,25 @@ shouldTerminate.OptState = function(opt.state) {
   # "fallback"
   return(list(term = FALSE, message = NA_character_, code = NA_character_, progress = progress))
 }
+
+
+# Helper which checks identification termination criterion. 
+#
+# @param opt.state [\code{OptState}]\cr
+shouldTerminateIdentification.OptState = function(opt.state) {
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+  stop.conds = control$stop.conds.identification
+  progress = NULL
+
+  for (i in seq_along(stop.conds)) {
+    stop.cond = stop.conds[[i]]
+    stop.obj = stop.cond(opt.state)
+    if (stop.obj$term) {
+      return(stop.obj)
+    }
+  }
+
+  # "fallback"
+  return(list(term = FALSE, message = NA_character_, code = NA_character_, progress = progress))
+}
diff --git a/R/term_conds.R b/R/term_conds.R
index d2d164ee8..6cd39c9a6 100644
--- a/R/term_conds.R
+++ b/R/term_conds.R
@@ -101,3 +101,39 @@ makeMBOTerminationMaxEvals = function(max.evals) {
     return(list(term = term, message = message, code = "term.feval", progress = (evals-init.evals) / (max.evals-init.evals)))
   }
 }
+
+
+# @title
+# Time budget for identification condition.
+#
+# @param time.budget [numeric(1)]
+#   Time budget in seconds.
+makeMBOTerminationIdentificationMaxBudget = function(time.budget) {
+  assertNumber(time.budget, na.ok = FALSE)
+  force(time.budget)
+  function(opt.state) {
+    time.used = as.numeric(getOptStateTimeUsedIdentification(opt.state), units = "secs")
+    term = (time.used > time.budget)
+    message = if (!term) NA_character_ else sprintf("Time budged for identification %f reached.", time.budget)
+    return(list(term = term, message = message, code = "term.time", progress = time.used / time.budget))
+  }
+}
+
+
+# @title
+# Maximum evaluations for identification 
+#
+# @param time.budget [numeric(1)]
+#   Time budget in seconds.
+makeMBOTerminationIdentificationMaxEvals = function(max.evals) {
+  max.evals = asInt(max.evals, na.ok = FALSE, lower = 1L)
+  force(max.evals)
+  function(opt.state) {
+    opt.path = getOptStateOptPath(opt.state)
+    evals = getOptPathLength(opt.path)
+    id.evals = sum(as.data.frame(opt.path)$prop.type == "identification")
+    term = (id.evals > max.evals)
+    message = if (!term) NA_character_ else sprintf("Maximal number of function evaluations %i for identification reached.", max.evals)
+    return(list(term = term, message = message, code = "term.feval", progress = id.evals / max.evals))
+  }
+}
diff --git a/R/utils.R b/R/utils.R
index 1fcb463b3..14b0c6ea2 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -11,6 +11,8 @@ loadPackages = function(control) {
     requirePackages("cmaesr", why = "proposePoints")
   if (control$n.objectives == 1L && control$propose.points > 1L && control$multipoint.method == "moimbo")
     requirePackages("emoa", why = "proposePoints")
+  if (control$identification.time.budget > 0)
+    requirePackages("mvtnorm", why = "intensification")
 }
 
 
diff --git a/man/setMBOControlNoisy.Rd b/man/setMBOControlNoisy.Rd
new file mode 100644
index 000000000..d317da586
--- /dev/null
+++ b/man/setMBOControlNoisy.Rd
@@ -0,0 +1,59 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/setMBOControlNoisy.R
+\name{setMBOControlNoisy}
+\alias{setMBOControlNoisy}
+\title{Set options for handling noisy functions.}
+\usage{
+setMBOControlNoisy(control, method = NULL, instances = NULL,
+  instance.param = NULL, instance.aggregation = NULL,
+  self.replicating = NULL, incumbent.nchallengers = NULL,
+  ocba.budget = NULL, ocba.initial = NULL)
+}
+\arguments{
+\item{control}{[\code{\link{MBOControl}}]\cr
+Control object for mbo.}
+
+\item{method}{[\code{character(1)}]\cr
+Which of the replication strategies should be used? Possible values are:
+\dQuote{fixed}: Every point is evaluated \code{instances} times. \cr
+\dQuote{incumbent}: Use an incumbent strategy as intensification strategy. 
+ The size of the set of additional challengers (apart from the incumbent) can be specified in \code{incumbent.nchallengers}. \cr
+\dQuote{ocba}: Distribution replication budget according to OCBA. 
+ The replication budget per iteration is specified in \code{ocba.budget},
+ the initial number of iterations per parameter is specified in \code{ocba.initial}. \cr}
+
+\item{instances}{[\code{integer(1)}]\cr
+How many instances of one parameter will be calculated?}
+
+\item{instance.param}{[\code{character(1)}]\cr
+What is the name of the function param that defines the instance?}
+
+\item{instance.aggregation}{[\code{function}]\cr
+Should data be aggregated per instance? If yes, a function (e. g. mean) needs to be specified.}
+
+\item{self.replicating}{[\code{logical(1)}]\cr
+TRUE if the function returns a vector of noisy results for one input. Then \code{instances} specifies the length of the result we expect.}
+
+\item{incumbent.nchallengers}{[\code{integer(1)}]\cr
+The size of the set of additional challengers (apart from the incumbent), defaults to \code{0}.}
+
+\item{ocba.budget}{[\code{integer(1)}]\cr
+The budget that is allocated in each iteration per Optimal Computing Budget Allocation (OCBA) rule, defaults to 10.}
+
+\item{ocba.initial}{[\code{integer(1)}]\cr
+The number of initial replications at each new point, defaults to \code{3}. 
+This needs to be larger than 1, since OCBA requires an initial variance estimate at each point.}
+}
+\value{
+[\code{\link{MBOControl}}].
+}
+\description{
+Extends an MBO control object with options for handling noisy functions.
+}
+\seealso{
+Other MBOControl: \code{\link{makeMBOControl}},
+  \code{\link{setMBOControlInfill}},
+  \code{\link{setMBOControlMultiObj}},
+  \code{\link{setMBOControlMultiPoint}},
+  \code{\link{setMBOControlTermination}}
+}
diff --git a/tests/testthat/test_identification.R b/tests/testthat/test_identification.R
new file mode 100644
index 000000000..a71985c90
--- /dev/null
+++ b/tests/testthat/test_identification.R
@@ -0,0 +1,69 @@
+context("mbo identification of final points")
+
+test_that("final OCBA identification works with time constraint", {
+  ps = makeNumericParamSet("x", 1L, -3, 3)
+  fun = makeSingleObjectiveFunction(
+    fn = function(x) (x + 0.5)^2 + 5 * sin(3 * (x + 0.5)) + rnorm(1, sd = 2),
+    par.set = ps, 
+    noisy = TRUE
+  )
+
+  design = generateDesign(n = 3L, par.set = ps)
+
+  ctrl = makeMBOControl()
+  ctrl = setMBOControlTermination(ctrl, time.budget = 2L, identification.time.budget = 2L)
+  ctrl = setMBOControlNoisy(ctrl, method = "ocba", ocba.initial = 2L, ocba.budget = 0L, identification.pcs = 1)
+
+  # check if time works correctly 
+  time = Sys.time()
+  or = mbo(fun, design = design, control = ctrl, show.info = TRUE)
+  timeend = as.numeric(difftime(Sys.time(), time), units = "secs")
+
+  # check that time budget is correct 
+  expect_true(timeend <= ctrl$time.budget + ctrl$identification.time.budget + 5)
+
+  # check if each point is at least evaluated twice
+  ds = as.data.table(or$opt.path)
+  ds = ds[, runs := .N, by = "x"]
+  expect_true(all(ds$runs >= 2L))
+})
+
+test_that("final OCBA identification works with pcs constraint", {
+  ps = makeNumericParamSet("x", 1L, -3, 3)
+  fun = makeSingleObjectiveFunction(
+    fn = function(x) (x + 0.5)^2 + 5 * sin(3 * (x + 0.5)) + rnorm(1, sd = 2),
+    par.set = ps, 
+    noisy = TRUE
+  )
+
+  ctrl = makeMBOControl()
+  ctrl = setMBOControlTermination(ctrl, time.budget = 2L, identification.time.budget = 10000L)
+  ctrl = setMBOControlNoisy(ctrl, method = "ocba", ocba.initial = 2L, ocba.budget = 0L, identification.pcs = 0.01)
+
+  or = mbo(fun, control = ctrl, show.info = TRUE)
+
+  # check if time works correctly 
+  pcs = mlrMBO:::calculatePCS(or$final.opt.state)
+  expect_true(pcs > ctrl$noisy.identification.pcs)
+})
+
+
+test_that("identification works with budget constraint", {
+  ps = makeNumericParamSet("x", 1L, -3, 3)
+  fun = makeSingleObjectiveFunction(
+    fn = function(x) (x + 0.5)^2 + 5 * sin(3 * (x + 0.5)) + rnorm(1, sd = 2),
+    par.set = ps, 
+    noisy = TRUE
+  )
+
+  ctrl = makeMBOControl()
+  ctrl = setMBOControlTermination(ctrl, max.evals = 10L, identification.max.evals = 9L)
+  ctrl = setMBOControlNoisy(ctrl, method = "ocba", ocba.initial = 2L, ocba.budget = 0L, identification.pcs = 1)
+
+  or = mbo(fun, control = ctrl, show.info = TRUE)
+
+  # check if time works correctly 
+  ds = as.data.table(or$opt.path)
+  expect_true(sum(ds$prop.type == "identification") == 12L)
+})
+
diff --git a/tests/testthat/test_intensify.R b/tests/testthat/test_intensify.R
new file mode 100644
index 000000000..d4416d0ac
--- /dev/null
+++ b/tests/testthat/test_intensify.R
@@ -0,0 +1,96 @@
+context("mbo intensification")
+
+test_that("mbo works with incumbent strategy", {
+  ps = makeNumericParamSet("x", 1L, -3, 3)
+  fun = smoof::makeSingleObjectiveFunction(
+	  fn = function(x) (x + 0.5)^2 + 5 * sin(3 * (x + 0.5)) + rnorm(1, sd = 2),
+	  par.set = ps, 
+	  noisy = TRUE
+	)
+
+  ctrl = makeMBOControl()
+  ctrl = setMBOControlTermination(ctrl, iters = 5L)
+  ctrl = setMBOControlNoisy(ctrl, method = "incumbent", incumbent.nchallengers = 2L)
+  or = mbo(fun, control = ctrl)
+  opdf = as.data.frame(or$opt.path)
+
+  # incumbent and proposed point are evaluated once in each iteration
+  expect_true(all(table(opdf$prop.type)[c("incumbent", "infill_cb")] == 5L))
+  # exactly 3 challengers are evaluated in each iteration (1 new point + 2 old points)  
+  expect_true(all(setDT(opdf)[opdf$prop.type == "challenger", .(length(unique(x))), by = "dob"]$V1 == 3L))
+
+  # incumbent of iteration i + 1 should be the best point of iteration i
+  for (i in 1:5L) {
+    curropdf = opdf[opdf$dob < i, .(ymean = mean(y)), by = x]
+    xinc = opdf[opdf$prop.type == "incumbent", "x"][i]
+    expect_true(curropdf[which.min(curropdf$ymean), ]$x == xinc)
+  }
+})
+
+
+test_that("mbo works with ocba replication strategy", {
+  ps = makeNumericParamSet("x", 1L, -3, 3)
+  fun = makeSingleObjectiveFunction(
+	  fn = function(x) (x + 0.5)^2 + 5 * sin(3 * (x + 0.5)) + rnorm(1, sd = 2),
+	  par.set = ps, 
+	  noisy = TRUE
+	)
+
+	ctrl = makeMBOControl()
+	ctrl = setMBOControlTermination(ctrl, iters = 5L)
+	ctrl = setMBOControlNoisy(ctrl, method = "ocba", ocba.budget = 3L, ocba.initial = 2L)
+
+	or = mbo(fun, control = ctrl)
+	opdf = as.data.frame(or$opt.path)
+
+	# in each iteration we do 12L replications
+	opdfs = setDT(opdf)[, .(N = sum(prop.type == "OCBA")), by = dob]
+	expect_true(all(opdfs[opdfs$dob > 0, "N"] == 3L))
+
+	# each point is evaluated 3L times at least
+	opdfs = setDT(opdf)[, .N, by = "x"]
+	expect_true(all(opdfs$N >= 2L))
+
+})
+
+
+test_that("mbo with replication works in multidimensional case", {
+  f = makeAckleyFunction(5L)
+
+  fun = makeSingleObjectiveFunction(
+  	fn = function(x) f(x) + rnorm(1, sd = 2), 
+  	par.set = getParamSet(f)
+  	)
+
+  ctrl = makeMBOControl()
+  ctrl = setMBOControlTermination(ctrl, iters = 5L)
+  ctrl = setMBOControlNoisy(ctrl, method = "ocba", ocba.budget = 12L, ocba.initial = 3L)
+
+  or = mbo(fun, control = ctrl)
+  opdf = as.data.frame(or$opt.path)
+
+  # in each iteration we do 12L replications
+  opdfs = setDT(opdf)[, .(N = sum(prop.type == "OCBA")), by = dob]
+  expect_true(all(opdfs[opdfs$dob > 0, "N"] == 12L))
+
+  # each point is evaluated 3L times at least
+  opdfs = setDT(opdf)[, .N, by = eval(paste("x", 1:5, sep = ""))]
+  expect_true(all(opdfs$N >= 3L))
+})
+
+test_that("per instance aggregation works for different functions", {
+  f = makeAckleyFunction(1L)
+  fun = makeSingleObjectiveFunction(
+  	fn = function(x) f(x) + rnorm(1, sd = 2), 
+  	par.set = getParamSet(f)
+  )
+
+  ctrl = makeMBOControl()
+  ctrl = setMBOControlTermination(ctrl, iters = 3L)
+  ctrl = setMBOControlNoisy(ctrl, method = "incumbent", instance.aggregation = median)
+
+  or = mbo(fun, control = ctrl)
+  opdf = as.data.frame(or$opt.path)
+
+  expect_true(all(length(or$models[[1]]$subset) == length(unique(opdf$x))))
+})
diff --git a/tests/testthat/test_mbo_km.R b/tests/testthat/test_mbo_km.R
index ededfc171..49885174a 100644
--- a/tests/testthat/test_mbo_km.R
+++ b/tests/testthat/test_mbo_km.R
@@ -1,6 +1,6 @@
-context("mbo km")
+context("mbo noisy")
 
-test_that("mbo works with km", {
+test_that("mbo works with multiple instances of noisy problems", {
   des = testd.fsphere.2d
   des$y = apply(des, 1, testf.fsphere.2d)
   learner = makeLearner("regr.km", nugget.estim = TRUE)
diff --git a/tests/testthat/test_mbo_noisy.R b/tests/testthat/test_mbo_noisy.R
new file mode 100644
index 000000000..9cb433858
--- /dev/null
+++ b/tests/testthat/test_mbo_noisy.R
@@ -0,0 +1,52 @@
+context("mbo noisy")
+
+test_that("mbo works with multiple instances of noisy problems", {
+  ps = makeNumericParamSet("x", 1, -7, 7)
+  fun = smoof::makeSingleObjectiveFunction(
+    fn = function(x) x^2 + rnorm(1, 0, 0.5),
+    par.set = ps,
+    noisy = TRUE
+  )
+  ctrl = makeMBOControl()
+  ctrl = setMBOControlTermination(ctrl, iters = 5L)
+  ctrl = setMBOControlInfill(ctrl, crit = crit.aei, opt.focussearch.points = 100L)
+  ctrl = setMBOControlNoisy(ctrl, instances = 5L)
+  or = mbo(fun, control = ctrl)
+  opdf = as.data.frame(or$opt.path)
+  expect_true(all(table(opdf$x) == 5))
+})
+
+test_that("mbo works with multiple fixed instances of noisy problems", {
+  ps = makeNumericParamSet("x", 1, -7, 7)
+  fun = smoof::makeSingleObjectiveFunction(
+    fn = function(x) x$x^2 + rnorm(1, (x$i - 2)/100, 0.05),
+    par.set = ps,
+    noisy = TRUE,
+    has.simple.signature = FALSE
+  )
+  ctrl = makeMBOControl()
+  ctrl = setMBOControlTermination(ctrl, iters = 5L)
+  ctrl = setMBOControlInfill(ctrl, crit = crit.aei, opt.focussearch.points = 100L)
+  ctrl = setMBOControlNoisy(ctrl, instances = 5L, instance.param = "i")
+  or = mbo(fun, control = ctrl)
+  opdf = as.data.frame(or$opt.path)
+  expect_true(all(opdf$i %in% 1:5))
+  expect_true(all(table(opdf$x) == 5))
+})
+
+test_that("mbo works with self replicating instances of noisy problems", {
+  ps = makeNumericParamSet("x", 1, -7, 7)
+  fun = smoof::makeSingleObjectiveFunction(
+    fn = function(x) x^2 + rnorm(5, 0.01),
+    par.set = ps,
+    noisy = TRUE
+  )
+  ctrl = makeMBOControl()
+  ctrl = setMBOControlTermination(ctrl, iters = 5L)
+  ctrl = setMBOControlInfill(ctrl, crit = crit.aei, opt.focussearch.points = 100L)
+  ctrl = setMBOControlNoisy(ctrl, instances = 5L, self.replicating = TRUE)
+  or = mbo(fun, control = ctrl)
+  opdf = as.data.frame(or$opt.path)
+  expect_true(all(opdf$noisy.repl %in% 1:5))
+  expect_true(all(table(opdf$x) == 5))
+})