Goal-Oriented Lower-Tail Calibration of Gaussian Processes for Bayesian

Optimization

Aurélien Pion 1 2 Emmanuel Vazquez 2

Abstract optimization (BO) (see, e.g., Mockus et al., 1978; Jones

et al., 1998; Villemonteix et al., 2009; Srinivas et al., 2010).
arXiv:2605.20145v1 [[Link]] 19 May 2026

Bayesian optimization (BO) selects evaluation

In the single-objective setting, BO aims at identifying a
points for expensive black-box objectives using
global minimizer x⋆ ∈ arg minx∈X f (x) using a limited
Gaussian process (GP) predictive distributions.
evaluation budget, modeling f a priori as a Gaussian process
Kernel choice and hyperparameter selection can
(GP), denoted by ξ ∼ GP(m, k), with mean function m
lead to miscalibrated predictive distributions and
and covariance kernel k. Given an initial set of evaluations,
an inappropriate exploration–exploitation trade-
the GP posterior yields predictive distributions for f : for
off. For minimization, sampling criteria such as
each x ∈ X, it provides a predictive cumulative distribution
expected improvement (EI) depend on the predic-
function (CDF) F̂n (· | x) for the unknown value f (x) at
tive distribution below the current best value, so
iteration n. At each iteration, the next evaluation point is
lower-tail miscalibration directly affects the sam-
selected by maximizing a sampling criterion, also referred
pling decision. This article studies goal-oriented
to as an acquisition function, that depends on this predictive
calibration of GP predictive distributions below
distribution, such as expected improvement (EI) (Mockus
a low threshold t in the noiseless setting, for
et al., 1978; Schonlau & Welch, 1996).
standard GP models with hyperparameters se-
lected by maximum likelihood. A framework To avoid unnecessary evaluations of f , the predictive dis-
for predictive reliability below t is introduced, tributions should reflect uncertainty at each iteration n in
based on two notions of spatial calibration: occur- order to guide the selection of new evaluation points toward
rence calibration over the design space and thresh- the minimum. In practice, GP predictive distributions are
olded µ-calibration on sublevel sets of the form sensitive to kernel choice and hyperparameter misspecifica-
{x ∈ X, f (x) ≤ t}. Building on this framework, tion, which can lead to over- or under-estimated uncertainty
we propose tcGP, a post-hoc method that cali- (Pion & Vazquez, 2025). Such misspecification can mis-
brates GP predictive distributions below t, and we lead the sampling criterion and substantially degrade the
show that the resulting EI-based global optimiza- performance of BO (Bogunovic & Krause, 2021).
tion algorithm remains dense in the design space.
While improving GP predictive performance over the whole
Experiments on standard benchmarks show im-
domain is desirable (Bogunovic & Krause, 2021; Guo et al.,
proved lower-tail calibration and BO performance
2021; Tuo & Wang, 2022; Deshpande et al., 2024; Tom
relative to standard GP models and globally cali-
et al., 2025), BO can benefit from goal-oriented predic-
brated GP models.
tion, in particular from accurate predictions below a low
threshold t ∈ R. Here, t is a data-driven threshold chosen
in the lower tail of the observations, e.g., as an empirical
1. Introduction δ-quantile of the observed responses for some δ ∈ (0, 1].
This article considers the minimization of an expensive- Along these lines, Petit et al. (2025) propose the relaxed
to-evaluate function f defined over a design space X ⊂ Gaussian process (RE GP), a GP-based model designed to
Rd , with noiseless evaluations (a query at x returns f (x) improve predictive accuracy below t. RE GP builds on GP in-
exactly). A widely used approach for this setting is Bayesian terpolation and relaxes the interpolation constraints above t
to improve predictions below this level. They report im-
1
Transvalor S.A., Biot, France 2 Univ. Paris-Saclay, CNRS, proved BO performance compared with standard GP-based
CentraleSupélec, L2S, Gif-sur-Yvette, France. Correspondence to: methods.
Emmanuel Vazquez <[Link]@[Link]>.
Following the idea of goal-oriented prediction below a
Proceedings of the 43 rd International Conference on Machine threshold t, this work investigates whether BO can benefit
Learning, Seoul, South Korea. PMLR 306, 2026. Copyright 2026
by the author(s). from probabilistic models that are calibrated below a low

1
 Goal-Oriented Calibration of GP for BO

Forrester 1d. Probablity of improvement. Scaled expected improvement.

GP interval 1.0 Proba. (GP Model) 1.0 EI (GP Model)
15 calibrated interval Proba. (Calibrated Model) EI (Calibrated Model)
posterior mean 0.8 f (x) < mn 0.8 f (x) < mn

Prob. improvement
10 f
minimum observed mn
0.6 0.6
evaluations

EI
z

5 f (x) < mn
0.4 0.4

0
0.2 0.2

−5
0.0 0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

Figure 1. Comparison between a standard GP and a model calibrated below qδ,n (TC GP), with δ = 0.3, on a standard one-dimensional test
function and with n = 10 evaluations. Left: observations with GP predictions; the current design does not include a point near the global
(0) (1)
minimizer. Middle: probabilities of improvement F̂n (mn | x) (GP, black) and F̂n (mn | x) (TC GP, red), where mn = mini≤n f (Xi ).
Right: EI for the GP and TC GP, each rescaled by its maximum over X. Here, calibration raises the predicted probability of improvement
in the sparsely sampled region near the global minimizer, where the current design has no evaluation, and EI increases there accordingly.
This exploration behavior seems appropriate in this example.

threshold t, and compares them with standard GP models A post-hoc calibration method for GP predictive distribu-
and with post-hoc globally calibrated GP predictive distribu- tions is also introduced, inspired by Pion & Vazquez (2025),
tions over X. By calibration below t, we mean that the lower termed calibrated prediction below t for GPs (TC GP). Stan-
tail of the predictive CDFs F̂n (· | x) is reliable for inputs dardized residuals are modeled with a parametric general-
drawn from µ. The average probability assigned to the event ized normal family, with parameters selected by optimizing
f (X) ≤ t, with X ∼ µ, matches its frequency, and, on the a criterion targeting spatial calibration below t. The re-
sublevel set {f (X) ≤ t}, the transformed lower-tail values sulting predictive distributions are used in BO with the EI
have the correct distribution. criterion. Analytical expressions for EI with TC GP are de-
rived, together with a convergence result. We also report
Consider Figure 1, which compares a GP model with hyper-
additional experiments using the UCB sampling criterion in
parameters selected by maximum likelihood and a model
the appendix.
calibrated below the empirical δ-quantile qδ,n (δ = 0.3), on
a standard one-dimensional test function (from Forrester We evaluate TC GP and RE GP empirically and find improved
et al., 2008). The middle panel displays the predicted proba- lower-tail calibration and BO performance relative to stan-
bilities of improvement dard GP models, globally calibrated GP models (Pion &
Vazquez, 2025), sequential calibration (Deshpande et al.,
F̂n(j) (mn | x), j ∈ {0, 1}, (1) 2024), and variants that target only one of the two calibra-
(j) tion notions.
where F̂n (· | x) is the predictive CDF at x, with j = 0
for the GP and j = 1 for the calibrated predictive distri- Section 2 reviews related work. Sections 3–5 introduce the
bution, and mn = mini≤n f (Xi ). In the shaded region, setting, TC GP, and the TC GP-based EI algorithm. Section 6
the GP assigns low lower-tail probabilities despite sparse reports experiments.
sampling, while the calibrated model assigns larger values,
which raises EI in that region. This example is constructed 2. Related Work
to illustrate the issue in a minimal setting. In higher dimen-
sions, similar behavior can arise when the selected hyperpa- Calibration and goal-oriented prediction Classical no-
rameters lead to local overconfidence in sparsely sampled tions of calibration, such as probabilistic calibration and
regions, with underestimated lower-tail probabilities and proper scoring rules, typically address global calibration of
reduced exploration. This motivates lower-tail calibration of predictive distributions over the domain X (Gneiting et al.,
predictive distributions, since improvement-based sampling 2007; Gneiting & Resin, 2023; Pion & Vazquez, 2025). Cal-
criteria are sensitive to tail miscalibration. ibration in BO with non-i.i.d. data is studied by Deshpande
et al. (2024), who propose online recalibration and confor-
This work makes two contributions. Inspired by the no-
mal procedures to maintain quantile calibration within the
tion of probabilistic tail calibration introduced by Allen
BO loop. Sahoo et al. (2021) introduced threshold cali-
et al. (2025), we define two notions of spatial calibration
bration for decision-making applications and proposed a
below t: occurrence calibration over X and thresholded
method for obtaining reliable estimates of threshold-based
µ-calibration on the sublevel set {x ∈ X, f (x) ≤ t}.

2
 Goal-Oriented Calibration of GP for BO

losses. Their notion is closely related to thresholded µ- We focus on the EI sampling criterion in this work. (We also
calibration considered here. This article is formulated in consider the UCB policy and report complementary results
a spatial setting over X and considers two notions of cali- in Appendix F.6.2.) For the family of predictive densities
bration below t: occurrence calibration over X and thresh- p̂n (· | x), with x ∈ X, the EI is defined as
olded µ-calibration on the sublevel set {x ∈ X, f (x) ≤ t}. Z
Friedli et al. (2025) introduced a GP-based sequential de- ρn (x) = (mn − z)+ p̂n (z | x)dz,
sign method built from the threshold-weighted continuous R
ranked probability score (twCRPS), and applied it to excur- where mn = min(f (X1 ), . . . , f (Xn )). EI is primarily con-
sion set estimation. trolled by the predicted probability mass below mn and
the shape of the predictive CDF on (−∞, mn ], motivating
Conformal prediction. Conformal prediction (CP) pro- calibration specifically in the lower tail rather than globally.
vides finite-sample coverage guarantees for a new response
Consider a fixed reference measure µ on X used to define
when the data are exchangeable (Vovk et al., 2005). Confor-
and assess spatial calibration (e.g., uniform on X). The mea-
mal predictive systems (CPS) extend CP to full predictive
sure µ may coincide with µ0 , but it is fixed independently
distributions (Vovk et al., 2019). CPS have been applied to
of the BO policy. For fixed n and Dn , let X ∼ µ be an aux-
GP interpolation by Pion & Vazquez (2025). Goal-oriented
iliary draw, independent of Dn . We write Pn (·) := PX (·)
CP methods target coverage on selected subsets of inputs
and En (·) := EX (·) for probability and expectation with
(Zhang & Candès, 2024; Jin & Ren, 2025), but typically
respect to this auxiliary draw X ∼ µ, with Dn fixed. The
yield piecewise-constant predictive distributions and still
subscript indicates that the integrands depend on Dn through
rely on exchangeability. CP has also been used for sequen-
F̂n (· | X).
tial optimization with noisy evaluations (Stanton et al., 2023;
Kim et al., 2025). These approaches are not tailored to lower-
3.2. A Generalized Normal Distribution Model for
tail BO criteria and can be computationally demanding. In
Prediction Errors
this work, we adopt a GP-specific post-hoc calibration ap-
proach inspired by Pion & Vazquez (2025) and adapt it to In this section, we depart from the Bayesian interpreta-
prediction below a threshold. tion of the GP model. The function f is treated as fixed
and deterministic, and the GP formalism is used only as
Alternative BO formulations. Misspecification in BO a device to construct, from Dn , a predictive mean fn (·)
and sequential decision-making has been studied by and a predictive variance σn2 (·). We retain from this con-
Neiswanger & Ramdas (2021), who analyzed the impact struction the associated Gaussian predictive

distributions
of global variance misspecification. Several works also F̂nGP (z | x) = Φ (z − fn (x))/σn (x) , x ∈ X, z ∈ R,
modify the probabilistic model or the sampling criterion where Φ is the standard normal CDF, with the convention
to improve BO, without targeting calibration of predictive F̂nGP (z | x) = 1{z ≥ fn (x)} when σn (x) = 0.
lower-tail probabilities. Picheny et al. (2019) considers BO To assess deviations from this Gaussian predictive distribu-
with ordinal observations (e.g., pairwise preferences or rank tion for X ∼ µ, define the standardized prediction error
information) and builds on variational ordinal GP regression
(Chu & Ghahramani, 2005). Picheny et al. (2022) focuses f (x) − fn (x)
Rn (x, f (x)) = , (2)
on quantiles and expectiles of the predictive distribution. σn (x)
with Rn (x, f (x)) = 0 if σn (x) = 0. Under the GP model,
3. Background and Problem Setup for any fixed x with σn (x) > 0, the conditional distribu-
tion of (ξ(x) − fn (x))/σn (x) given Dn is standard normal.
3.1. Problem Setting and Notation
In contrast, since f is deterministic, Rn (X, f (X)) is ran-
Bayesian optimization proceeds sequentially. We start from dom only through X ∼ µ and its distribution need not be
an initial dataset Dn0 = {(Xi , Zi )}ni=1
0
, where Zi = f (Xi ) Gaussian.
and the inputs Xi are drawn according to an initial design
Following Pion & Vazquez (2025), we model Rn (X, f (X))
measure µ0 on X. We then run n − n0 additional BO it-
with a generalized normal (GN) family, denoted by
erations, n > n0 . For i > n0 , the design points Xi are
GN (β, l, λ) with β > 0, λ > 0, and l ∈ R (density recalled
selected by maximizing a sampling criterion, so their distri-
in Appendix A). The parameter l is a location parameter,
bution is determined by the BO strategy. Using the dataset
λ controls dispersion, and β controls tail decay; β = 2
Dn we construct a family of predictive CDFs F̂n (· | x) for
corresponds to the Gaussian case and β = 1 to the Laplace
x ∈ X. This family is used to define the sampling criterion
case.
and thereby governs the choice of the next point Xn+1 to
evaluate f . For any x ∈ X and parameters β, l, λ, we define the predic-

3
 Goal-Oriented Calibration of GP for BO

tive CDF Recall the probability integral transform (PIT), which mo-
tivates the definition below: if a random variable Z has
z − fn (x)
 
F̂nβ,l,λ (z | x) = Θβ,l,λ , z ∈ R, (3) a continuous and strictly increasing CDF F , then the PIT
σn (x) F (Z) is uniformly distributed on [0, 1]. Hence, if a pre-
dictive CDF F̂ matches the distribution of Z, then the PIT
with the convention F̂nβ,l,λ (z | x) = 1{z ≥ fn (x)} when F̂ (Z) is uniform: departures from uniformity indicate mis-
σn (x) = 0, where Θβ,l,λ is the CDF of GN (β, l, λ). calibration. Pion & Vazquez (2025) extend this idea to a
Pion & Vazquez (2025) fix l = 0 and select (β, λ) via spatial (design-marginal) setting by introducing the µ-PIT,
a Bayesian procedure inspired by tolerance-interval con- based on an auxiliary draw X ∼ µ. Here we apply the same
structions (Meeker et al., 2017), targeting µ-probabilistic idea conditionally on {f (X) ≤ t} through the truncated
calibration of the predictive distribution over X. The result- CDF F̂n,t (· | X).
ing method is the Bayesian calibration of residuals for GPs Assume that µ({x ∈ X : f (x) ≤ t}) > 0 and µ({x ∈ X :
(BCR GP). f (x) = t}) = 0. The second condition avoids an atom
at one in the thresholded PIT. The µ–probability integral
4. Spatial Calibration below a Threshold transform restricted to the thresholded region (µ-tPIT) is
defined by
4.1. Statistical Formulation (
F̂n,t (f (X) | X), f (X) ≤ t,
In this section, we introduce notions of calibration below a Ut := (6)
threshold t ∈ R with respect to a fixed measure µ on X. For 1, f (X) > t.
a given n, consider predictive CDFs F̂n (· | x), x ∈ X, such
that for µ-almost all x the map z 7→ F̂n (z | x) is continuous We say that F̂n is µ-probabilistically calibrated below t
and strictly increasing on R. We assess calibration of the if Ut | {f (X) ≤ t} is uniform on [0, 1]. This is the
predictive distributions in the lower tail up to t. thresholded analogue of µ–probabilistic calibration (Pion
& Vazquez, 2025): among points such that f (X) ≤ t, the
For z ≤ t, the predictive CDF can be decomposed as values F̂n,t (f (X) | X) should be uniformly distributed on
[0, 1].
F̂n (z | x) = F̂n (t | x) F̂n,t (z | x), (4)
For u ∈ [0, 1], define
where F̂n (t | x) is the predicted probability mass below t,
Gµ,t (u) := Pn (Ut ≤ u | f (X) ≤ t) . (7)
and F̂n,t (· | x) is the predictive CDF truncated below t,
such that, for all z ≤ t, Then µ-probabilistic calibration below t is equivalent to
Gµ,t (u) = u for all u ∈ [0, 1].
F̂n (z | x)
F̂n,t (z | x) = PF̂n (·|x) (Z ≤ z | Z ≤ t) = , (5)
F̂n (t | x) Occurrence calibration. Thresholded µ-calibration con-
strains the distribution of F̂n,t (f (X) | X) conditional on
with Z ∼ F̂n (· | x). We therefore introduce two comple- {f (X) ≤ t}. It does not constrain the total probability mass
mentary notions: thresholded µ-calibration, targeting the assigned to {f (X) ≤ t}, which is governed by F̂n (t | x).
shape of the predictive distribution within the thresholded re-
Let pt denote the excursion probability for X ∼ µ,
gion through F̂n,t , and occurrence calibration, targeting the
probability mass assigned to the region through F̂n (t | x) pt = P (f (X) ≤ t) . (8)
averaged over X.
The occurrence discrepancy is
This decomposition parallels the severity–occurrence split
in probabilistic tail calibration (Allen et al., 2025). In that  
rt,n = pt − En F̂n (t | X) , X ∼ µ. (9)
setting, severity refers to calibration of the conditional dis-
tribution given Z ≤ t, and occurrence refers to calibration
Equivalently, pt = µ({x ∈ X : f (x) ≤ t}) and En [F̂n (t |
of the predictive probability of the event {Z ≤ t}. The def- R
X)] = X F̂n (t | x) µ(dx). We say that F̂n is occurrence-
initions introduced here adapt this split to a spatial setting
calibrated at threshold t if rt,n = 0.
over X with a fixed threshold t.
4.2. Metrics for Prediction below a Threshold
Thresholded µ-calibration. Thresholded µ-calibration is
a tail version of µ-probabilistic calibration (Pion & Vazquez, We assess thresholded µ–probabilistic calibration by com-
2025), obtained by restricting attention to the thresholded paring the distribution of the µ–tPIT to the uniform distribu-
region {f (X) ≤ t} for an auxiliary draw X ∼ µ. tion using the Kolmogorov–Smirnov distance. This yields

4
 Goal-Oriented Calibration of GP for BO

the tKS–PIT metric: Remark 4.1. Calibration alone does not guarantee useful
uncertainty quantification. A well-calibrated predictive dis-
JtKS-PIT (F̂n | t) := sup |Gµ,t (u) − u| . (10)
u∈[0,1] tribution can be overly diffuse. Proper scoring rules assess
calibration and sharpness jointly (see Appendix C.1).
Weighted LOO µ-tPIT and tKS–PIT. The tKS–PIT and
occurrence discrepancy rt are defined with respect to the 4.3. Calibration below a Threshold in BO
reference measure µ on X. In principle, they should be
At iteration n, EI depends on the predictive lower tail below
evaluated on an independent dataset (Xj′ , f (Xj′ ))L
j=1 with the current best value mn , through F̂n (mn | x) and the
Xj′ ∼ µ. In BO, such an evaluation set is typically unavail-
restriction of F̂n (· | x) to (−∞, mn ]. However, calibration
able. We therefore rely on the observed BO data Dn and
at level mn cannot be assessed from Dn : by definition, no
use leave-one-out (LOO) predictive CDFs. Since the BO
observation satisfies f (Xi ) < mn . We therefore impose
design points are not distributed according to µ in general,
calibration below a higher threshold tn , chosen from the
we also allow for reweighting to approximate µ–averages.
empirical δ-quantile qδ,n of (f (X1 ), . . . , f (Xn )) for some
Let Dn = {(Xi , Zi )}ni=1 with Zi = f (Xi ). For each δ ∈ (0, 1]. In the BO algorithm, this threshold can be kept
i ∈ {1, . . . , n}, let F̂n,−i (· | x) denote the predictive CDF fixed when too few observations fall below the new quantile.
constructed from Dn \ {(Xi , Zi )}. Fix t ∈ R. A formal consistency statement showing that occurrence
Let P(wi )ni=1 be nonnegative weights and set w̃i = calibration at mn forces the average predicted improvement
n probability to vanish as mn approaches the optimum is
wi / j=1 wj . The weights are chosen so that weighted av-
erages over (Xi ) approximate µ–averages. The unweighted given in Appendix D.
case corresponds to wi = 1. In this work, (wi ) is obtained
by density-ratio estimation based on a density estimate of 5. A Goal-Oriented Calibration Method
the BO design distribution; details are given in Appendix E.3
(paragraph Density-ratio estimation and weights). 5.1. Description of the Method
Assume that F̂n,−i (t | Xi ) > 0 whenever Zi ≤ t. Define In this section, we introduce a procedure to select the param-
the LOO empirical µ-tPIT by eters β and λ of the GN predictive model. We target thresh-
(  olded µ-calibration below a fixed threshold t, assessed via
LOO
F̂n,−i (Zi | Xi ) F̂n,−i (t | Xi ), Zi ≤ t, the truncated CDF F̂n,t on {f (X) ≤ t}, and occurrence cal-
Ut,i =  
1, Zi > t. ibration at t, assessed by matching pt with En F̂n (t | X) .
(11) Following Pion & Vazquez (2025), we fix l = 0, so that
calibration changes the shape of the predictive distribution
We estimate pt using the weighted empirical frequency
without shifting its mean away from the GP mean fn . Ac-
n
w
X cordingly, we write F̂nβ,λ (· | x) instead of F̂nβ,0,λ (· | x).
p̂t,n = w̃i 1{Zi ≤ t}. (12)
i=1
Parameter selection criterion. Let Utβ,λ be the µ-tPIT
w
Assume that > 0 and define
p̂t,n β,λ
associated with the truncated predictive CDF F̂n,t (· | X).
Pn LOO
w̃i 1{Zi ≤ t} 1{Ut,i ≤ u} Thresholded µ–probabilistic calibration below t is equiva-
GLOO,w
t,n (u) = i=1 w . (13) lent to
p̂t,n
 
The weighted empirical tKS–PIT metric is Pn Utβ,λ ≤ u f (X) ≤ t = u, u ∈ [0, 1]. (16)
LOO,w
JtKS-PIT,n (F̂n | t) = sup GLOO,w
t,n (u) − u . (14)
u∈[0,1] To incorporate occurrence calibration at t, define the occur-
rence ratio h i
Weighted estimation of pt and occurrence discrepancy. En F̂nβ,λ (t | X)
The corresponding weighted LOO occurrence discrepancy κβ,λ
t = . (17)
pt
estimator is
n
X The numerator is the µ-average predicted probability of
LOO,w w {f (X) ≤ t}, while the denominator is the µ-probability of
rt,n = p̂t,n − w̃i F̂n,−i (t | Xi ) . (15)
i=1 {f (X) ≤ t}, so κβ,λt = 1 is equivalent to occurrence cali-
bration at t. Note that κβ,λ
t > 1 indicates an overestimation
When an independent test set is available (as in our synthetic β,λ
of pt and κt < 1 an underestimation.
benchmarks), we also report direct test-set versions of these
metrics; see Appendix E.4. We compare u 7→ Pn (Utβ,λ ≤ u | f (X) ≤ t) to u 7→

5
 Goal-Oriented Calibration of GP for BO

u κβ,λ
t and define Predictive distribution and EI. Fix (β, λ). Conditional
  on Dn and for x ∈ X, TC GP uses the predictive CDF
J(β, λ) = sup Pn Utβ,λ ≤ u f (X) ≤ t − u κβ,λ
t . F̂nβ,λ (· | x) defined in (3). Let Zn,x be a predictive random
u∈[0,1] variable with this CDF. With the GN model,
(18)
By construction, J(β, λ) = 0 when F̂nβ,λ is µ–


Zn,x ∼ GN β, fn (x), λ σn (x) , (22)
probabilistically calibrated below t and occurrence-
calibrated at t. When κβ,λt = 1, J(β, λ) reduces to the with the convention that Zn,x = fn (x) when σn (x) = 0.
tKS–PIT metric. Appendix F.2 reports alternative parameter
selection criteria. We use J in the experiments. The EI criterion induced by F̂nβ,λ (· | x) is
 
ρn (x) = E (mn − Zn,x )+ | Dn , (23)
LOO approximation. The criterion J(β, λ) involves µ–
expectations and the excursion probability pt , which are where the expectation is taken with respect to the predictive
unknown given only Dn . We use a weighted leave-one-out distribution of Zn,x conditional on Dn . Proposition 5.1
(LOO) approximation. Let GLOO,w t,n (u; β, λ) be defined by yields the closed form
β,λ
(13), with F̂n,−i (· | x) replaced by F̂n,−i (· | x). Define

ρn (x) = γ mn − fn (x), λ σn (x), β , (24)
Pn β,λ
i=1 w̃i F̂n,−i (t | Xi )
bβ,λ
κ t,n = w , (19) where γ is defined in (26).
p̂t,n
Proposition 5.1. If Z ∼ GN (β, l, λ), then for any a ∈ R,
w
assuming p̂t,n > 0. The resulting approximation of (18) is  
E (a − Z)+ = γ(a − l, λ, β), (25)
LOO,w
Jt,n (β, λ) = sup GLOO,w
t,n bβ,λ
(u; β, λ) − u κt,n .
u∈[0,1] where, for λ > 0,
(20)  
z λ 2 z β
The parameters β and λ are selected by solving γ(z, λ, β) = z Θβ + Γ , , (26)
λ 2 Γ(1/β) β λ
LOO,w
(β ∗ , λ∗ ) = arg min Jt,n (β, λ), (21) and γ(z, 0, β) = max(z, 0). Here Θβ denotes the CDF
(β,λ)∈[β0 ,β1 ]×[λ0 ,λ1 ]
of GN (β, 0, 1) and Γ(·, ·) is the upper incomplete gamma
where 0 < λ0 < λ1 and 0 < β0 < β1 . We refer to the function. For λ > 0, γ(·, λ, β) is continuous and satisfies
method as calibrated prediction below t for GPs (TC GP). γ(z, λ, β) > 0 for all z ∈ R.
The complete TC GP procedure, including its integration
into BO, is described in Appendix E.3. Proof. See Appendix B.1.
Remark 5.2. For β = 2, the GN distribution coincides with
5.2. EI with TC GP a Gaussian distribution, and (26) reduces to the standard EI
expression; see Appendix B.2.
Section 3 defined EI and the BO loop. We now describe the
EI algorithm obtained by using TC GP predictive distribu-
BO update. At iteration n, we set t = tn and obtain
tions.
(βn , λn ) by solving (21). The next evaluation point is se-
lected as
Choice of the threshold. The objective is to obtain predic-
Xn+1 ∈ arg max ρn (x). (27)
tive distributions that are well calibrated below the current x∈X
best value mn . Calibrating at level mn is not feasible from Additional implementation details are given in Ap-
Dn , since no observation satisfies Zi < mn by definition. pendix E.3.
We therefore calibrate below a higher threshold tn chosen
in the lower tail of the observed responses. Let qδ,n be
5.3. Convergence of EI with TC GP and Fixed GP
the empirical δ-quantile of (Z1 , . . . , Zn ). At the initial it-
Hyperparameters
eration, we set tn0 = qδ,n0 . We use qδ,n as the candidate
threshold when p̂qwδ,n ,n ≥ pmin , in which case tn = qδ,n . We prove a convergence result for EI when predictive dis-
If not, we keep the previous threshold. Here, pmin is the tributions are obtained with TC GP from a zero-mean GP
minimum estimated excursion probability below a candidate model with fixed covariance kernel k. This is a sanity check:
threshold that we allow. This avoids calibrating on too few with fixed hyperparameters and bounded TC GP parameters,
points when the design concentrates near a minimizer and EI still produces a dense sequence of evaluation points in X.
the lower-tail excursion probability becomes small. Proposition 5.3 extends the exploration result of Vazquez

6
 Goal-Oriented Calibration of GP for BO

& Bect (2010) by replacing the Gaussian predictive distri- All methods use a GP model with constant mean and an
bution with the GN family, with parameters (βn , λn ) con- anisotropic Matérn kernel (Appendix E.1), with hyperpa-
strained to a compact set. rameters selected by maximum likelihood at each iteration.
Experiments are implemented with gpmp (Vazquez, 2026).
We assume that k satisfies the no-empty ball (NEB) property
The benchmark uses standard deterministic test functions
(Vazquez & Bect, 2010): for any sequence (Xn )n≥1 in X
(Appendix E.7). Each run starts from an initial design of size
and any x ∈ X, x is an adherent point of {Xn : n ≥ 1} iff
n0 = 10d drawn uniformly on X, followed by EI-based BO
σn2 (x) → 0 as n → ∞.
iterations. EI is maximized using a sequential Monte Carlo
Proposition 5.3. Assume that X is compact and that k is procedure available in gpmp (Appendix E.2). For each test
continuous, stationary, strictly positive definite, and satisfies function, we generate 100 independent initial datasets.
the NEB property. Let H be the associated RKHS and as-
sume f ∈ H. Let fn and σn be the zero-mean noiseless ker- BO performance is summarized by the excursion probability
nel interpolant and kriging standard deviation constructed below the current best value mn : for each run and each
from Dn . Assume that there exists M ≥ 1 such that for all iteration n, we consider pmn = P (f (X) ≤ mn ), X ∼
n ≥ M, U(X), and estimate pmn by subset simulation (Bect et al.,
(βn , λn ) ∈ [β0 , β1 ] × [λ0 , λ1 ], 2017) as the main performance metric. We then report, as a
function of n, the median and the 10% and 90% quantiles
with 0 < β0 < β1 and 0 < λ0 < λ1 . Let (Xn )n≥1 be of pmn across the 100 runs.
generated by EI, that is
6.2. Thresholded and Occurrence Calibration:
Xn+1 ∈ arg max ρn (x), Empirical Comparison
x∈X

where ρn is EI computed from the TC GP predictive CDFs We compare the roles of occurrence calibration, thresh-
with parameters (βn , λn ). Then (Xn ) is dense in X. olded tail-rank calibration, and their combination in BO.
We compare a standard GP model with three TC GP variants.
The maximizer exists because fn , σn , and γ are continuous, The first variant, TC GP-thres, selects (β, λ) by minimiz-
so ρn is continuous on the compact set X. ing the weighted LOO tKS–PIT, targeting thresholded µ–
calibration. The second, TC GP-occ, minimizes the weighted
Proof. See Appendix B.3. LOO occurrence discrepancy rt . The joint variant, TC GP,
uses the selection criterion of Section 5.
Proposition 5.3 implies that (Xn ) intersects every nonempty At each BO iteration, we report BO performance and test-
open subset of X. Since f is continuous on the compact set set calibration metrics evaluated at the current best value
X, for any global minimizer x∗ ∈ arg minx∈X f (x) there mn : twCRPS, tKS–PIT below mn , and rt (details about the
exists a subsequence (Xnj )j such that Xnj → x∗ , hence metrics provided in Appendices C.2 and E.4). The twCRPS
mn → minx∈X f (x). and rt use 1000 test points sampled in X, while tKS–PIT
Remark 5.4. In practice, GP hyperparameters are re- uses 900 test points conditioned on f (x) ≤ mn .
estimated at each iteration. Extending Proposition 5.3 to Figure 2 reports results for Ackley (d = 4) and Goldstein–
sequential hyperparameter selection is nontrivial (see, e.g., Price. On Ackley, TC GP-thres reduces tKS–PIT but does
Bull, 2011) and left for future work. not improve rt or twCRPS, and BO progresses more slowly.
TC GP-occ substantially reduces rt and improves optimiza-
6. Experiments tion performance despite a degraded tKS–PIT, suggesting
that correcting tail mass can matter more for EI than im-
6.1. Experimental Setup proving tail ranks alone. Neither single-component variant
We study how calibration below a threshold affects BO. We improves upon the standard GP model. The joint method
report results for variants targeting occurrence calibration, TC GP performs best. It rapidly reduces rt and improves
thresholded µ-calibration, or both, and we compare TC GP tKS–PIT early, leading to faster decreases of mn . As the de-
with RE GP (relaxed interpolation above tn ), BCR GP with sign concentrates near a minimizer and mn approaches the
density-ratio reweighting (Appendix E.5) for calibration optimum, calibration below mn relies on fewer test points
over X, and the online conformal recalibration method of and optimization slows down.
Deshpande et al. (2024) based on the quantile pinball loss On Goldstein–Price, all TC GP-based variants improve over
(onGP). For RE GP, we use δ = 0.25, and for TC GP, δ = the GP and quickly reduce twCRPS and rt . TC GP and
0.05 (see also Appendices F.3 and F.4). We fix pmin = TC GP-occ behave similarly and outperform TC GP-thres,
0.015. Additional diagnostics at tn and at the current best again indicating that occurrence calibration has a stronger
value mn are reported in Appendix F.5.

7
 Goal-Oriented Calibration of GP for BO

Excursion Prob. twCRPS Occurence discrepancy tKS–PIT

100
10−1
10−1
10−3 10−2 0.8
10−2
−3
10−3
10
10−5 0.6
10−4
10−4
10−5 10−5
10−7
10−6 0.4
10−6
10−9 10−7 10−7
0.2
0 10−8

1 50 100 150 200 1 50 100 150 200 1 50 100 150 200 1 50 100 150 200
steps steps steps steps
GP tcGP tcGP-occ tcGP-thres

(a) Ackley with d = 4

Excursion Prob. twCRPS Occurence discrepancy tKS–PIT
10−1 103 1.0
102 10−2
10−3 101 0.9
100
10−5 10−4
10−1
0.8
10−2
10−7 10−6
10−3
10−4 0.7
10−9 10−5 10−8
10−6
0.6
10−11 10−7 10−10
0
1 50 100 150 200 1 50 100 150 200 1 50 100 150 200 1 50 100 150 200
steps steps steps steps
GP tcGP tcGP-occ tcGP-thres

(b) Goldstein–Price
Figure 2. BO performance and calibration metrics. From left to right: median and 10%/90% quantiles across runs of the estimated
excursion probability pmn = P(f (X) ≤ mn ) with X ∼ U(X); median twCRPS; median occurrence discrepancy rt ; and median
tKS–PIT. Calibration metrics are evaluated on a test set at the current best value mn . Results are shown for a standard GP and three
TC GP variants using J, tKS–PIT, or rt as the parameter selection criterion (δ = 0.05).

Goldstein–Price Dixon–Price4 Rosenbrock6 Ackley4

10−1
10−2 10−2
10−3 10−3
10−4 10−4
−5
10 10−6
10−6 10−5
10−7 10−8
−8
10
10−9 10−10 10−7

−12
10−10
−11 10
10
10−9
10−14 10−12
10−13

0 50 100 150 200 0 50 100 150 200 0 50 100 150 200 0 50 100 150 200
steps steps steps steps
GP onGP bcrGP reGP δ = 0.25 tcGP δ = 0.05

Figure 3. BO performance summarized by the excursion probability below the current best value. For each run and iteration n, we estimate
pmn = P(f (X) ≤ mn ) with X ∼ U(X), and report the median and 10%/90% quantiles of pmn across runs for Goldstein–Price,
Dixon–Price (d = 4), Rosenbrock (d = 6), and Ackley (d = 4). Methods: GP, BCR GP, onGP, RE GP (δ = 0.25), TC GP (δ = 0.05).

effect on EI-driven BO than thresholded µ-calibration alone. 6.3. Comparison across Methods
Overall, these results suggest that accurate occurrence cali- We compare EI-based BO with a standard GP model, RE GP,
bration at mn is often important for BO, with thresholded TC GP, onGP, and BCR GP on four test functions: Ackley
µ-calibration providing an additional gain. Focusing only (d = 4), Dixon–Price (d = 4), Rosenbrock (d = 6), and
on the thresholded notion can yield poor optimization per- Goldstein–Price. Figure 3 reports the evolution of the me-
formance because it does not control systematic over- or dian and the 10%/90% quantiles of pmn = P(f (X) ≤ mn ),
under-estimation of improvement probabilities over X. where mn denotes the best observed value so far. To com-
plement these results, Appendix F.6 reports additional EI ex-

8
 Goal-Oriented Calibration of GP for BO

periments on Ackley, Crossintray, Dixon–Price, Goldstein– Impact Statement

Price, Hartmann, Michalewicz, Rosenbrock, Shekel, and
Perm, in dimensions ranging from d = 2 to d = 20, to- Improved lower-tail calibration can reduce the number of ex-
gether with additional results for the upper confidence bound pensive evaluations required in Bayesian optimization. The
(UCB) criterion. impact is application-dependent, since the same method-
ology can be used for beneficial or harmful optimization
On Goldstein–Price, Dixon–Price, and Rosenbrock, RE GP goals.
and TC GP improve over the standard GP model, with RE GP
performing best on these smoother objectives. On Ackley,
TC GP provides the largest gains.
Conflict of Interest Disclosure
BCR GP does not consistently improve over the standard GP Aurélien Pion is affiliated with TRANSVALOR S . A . The
model, suggesting that global calibration over X alone is article does not evaluate a proprietary product or model
not sufficient to improve EI. The online conformal method developed by TRANSVALOR S . A . The authors declare no
onGP is close to the standard GP model on these noiseless other financial conflicts of interest.
benchmarks, except on Ackley.
Appendix F.6 reports similar qualitative patterns for both
References
EI and UCB sampling criteria. On functions where the Allen, S., Bhend, J., Martius, O., and Ziegel, J. Weighted
standard GP model is already competitive, such as Ackley, verification tools to evaluate univariate and multivariate
Michalewicz, or Shekel, TC GP and RE GP provide only probabilistic forecasts for high-impact weather events.
small improvements. In these cases, calibration has little Weather and Forecasting, 38(3):499 – 516, 2023.
effect on BO performance.
Allen, S., Koh, J., Segers, J., and Ziegel, J. Tail calibration
Calibrating below a low threshold improves BO in these of probabilistic forecasts. J. Amer. Statist. Assoc., 120
settings. RE GP and TC GP show complementary strengths, (552):2796–2808, 2025.
with RE GP incurring a larger computational overhead than
TC GP (Appendix F.1). Auer, P., Cesa-Bianchi, N., and Fischer, P. Finite-time
analysis of the multiarmed bandit problem. Mach. Learn.,
47:235–256, 2002.
7. Discussion and Limitations
Bect, J., Li, L., and Vazquez, E. Bayesian subset simulation.
We introduced a goal-oriented calibration framework for SIAM/ASA Journal on Uncertainty Quantification, 5(1):
Bayesian optimization, built on two design-marginal no- 762–786, 2017.
tions: occurrence calibration at a threshold and thresholded
µ-calibration below that threshold. The ablation study indi- Bogunovic, I. and Krause, A. Misspecified Gaussian process
cates that the joint criterion is more consistent than variants bandit optimization. In Ranzato, M., Beygelzimer, A.,
targeting only one of the two notions. In the deterministic Dauphin, Y., Liang, P., and Vaughan, J. W. (eds.), Adv.
GP-based experiments, TC GP changes the behavior of EI Neural Inf. Process. Syst., volume 34, pp. 3004–3015.
and UCB and improves performance on several benchmarks, Curran Associates, Inc., 2021.
while remaining cheaper to train than RE GP.
Bull, A. D. Convergence rates of efficient global optimiza-
The study is restricted to the exact-evaluation setting. In tion algorithms. J. Mach. Learn. Res., 12(88):2879–2904,
noisy problems, the latent objective is no longer directly 2011.
observed, so the calibration targets considered here are not
directly accessible from the data and would need to be re- Chu, W. and Ghahramani, Z. Gaussian processes for ordinal
defined. The experiments are also mostly low-to-moderate regression. J. Mach. Learn. Res., 6(35):1019–1041, 2005.
dimensional, where the weighted leave-one-out estimates Deshpande, S., Marx, C., and Kuleshov, V. Online cali-
are more stable. Future work can study extensions to noisy brated and conformal prediction improves Bayesian op-
BO and to higher-dimensional problems. timization. In Dasgupta, S., Mandt, S., and Li, Y. (eds.),
Proceedings of The 27th International Conference on Ar-
Acknowledgements tificial Intelligence and Statistics, volume 238 of Proc.
Mach. Learn. Res., pp. 1450–1458. PMLR, 02–04 May
This work was supported by TRANSVALOR S . A . through a 2024.
CIFRE PhD agreement. We thank TRANSVALOR S . A . for
its support. Feliot, P., Bect, J., and Vazquez, E. A Bayesian approach to
constrained single- and multi-objective optimization. J.
Global Optim., 67(1-2):97–133, April 2016.

9
 Goal-Oriented Calibration of GP for BO

Forrester, A. I. J., Sóbester, A., and Keane, A. J. Engineering Neiswanger, W. and Ramdas, A. Uncertainty quantification
Design via Surrogate Modelling: A Practical Guide. John using martingales for misspecified Gaussian processes.
Wiley & Sons, 2008. In Feldman, V., Ligett, K., and Sabato, S. (eds.), Proceed-
ings of the 32nd International Conference on Algorithmic
Friedli, L., Gautier, A., Broccard, A., and Ginsbourger, Learning Theory, volume 132 of Proc. Mach. Learn. Res.,
D. CRPS-based targeted sequential design with applica- pp. 963–982. PMLR, 16–19 Mar 2021.
tion in chemical space, 2025. URL [Link]
org/abs/2503.11250. Nelder, J. A. and Mead, R. A simplex method for function
minimization. The Computer Journal, 7(4):308–313, 01
Gneiting, T. and Resin, J. Regression diagnostics meets 1965.
forecast evaluation: conditional calibration, reliability
diagrams, and coefficient of determination. Electron. J. Petit, S. J., Bect, J., and Vazquez, E. Relaxed Gaussian pro-
Stat., 17(2), January 2023. cess interpolation: a goal-oriented approach to Bayesian
optimization. Journal of Machine Learning Research, 26
Gneiting, T., Balabdaoui, F., and Raftery, A. E. Probabilistic (195):1–70, 2025.
forecasts, calibration and sharpness. J. R. Stat. Soc. Ser.
B Stat. Methodol., 69(2):243–268, 2007. Picheny, V., Vakili, S., and Artemev, A. Ordinal Bayesian
optimisation, 2019. URL [Link]
Gneiting, T. G. and Raftery, A. E. Strictly proper scoring
abs/1912.02493.
rules, prediction, and estimation. J. Am. Stat. Assoc., 102
(477):359–378, 2007. Picheny, V., Moss, H., Torossian, L., and Durrande, N.
Bayesian quantile and expectile optimisation. In Cussens,
Guo, Z., Ong, Y. S., and Liu, H. Calibrated and recalibrated
J. and Zhang, K. (eds.), Proceedings of the Thirty-Eighth
expected improvements for Bayesian optimization. Struct.
Conference on Uncertainty in Artificial Intelligence, vol-
Multidiscip. Optim., 64:3549–3567, 2021.
ume 180 of Proc. Mach. Learn. Res., pp. 1623–1633.
Jin, Y. and Ren, Z. Confidence on the focal: conformal PMLR, 01–05 Aug 2022.
prediction with selection-conditional coverage. J. R. Stat.
Pion, A. and Vazquez, E. Design-marginal calibration of
Soc. Ser. B. Stat. Methodol., 87(4):1239–1259, 04 2025.
Gaussian process predictive distributions: Bayesian and
Jones, D., Schonlau, M., and Welch, W. Efficient global conformal approaches, 2025. URL [Link]
optimization of expensive black-box functions. J. Global org/abs/2512.05611.
Optim., 13:455–492, 12 1998.
Sahoo, R., Zhao, S., Chen, A., and Ermon, S. Reliable
Kim, D., Zecchin, M., Park, S., Kang, J., and Simeone, decisions with threshold calibration. In Ranzato, M.,
O. Robust Bayesian optimization via localized online Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan,
conformal prediction. IEEE Trans. Signal Process., 73: J. W. (eds.), Adv. Neural Inf. Process. Syst., volume 34,
2039–2052, 2025. pp. 1831–1844. Curran Associates, Inc., 2021.

Lai, T. L. and Robbins, H. Asymptotically efficient adaptive Schonlau, M. and Welch, W. J. Global optimization with
allocation rules. Adv. Appl. Math., 6(1):4–22, 1985. nonparametric function fitting. In Proceedings of the
ASA, Section on Physical and Engineering Sciences, pp.
Matheson, J. E. and Winkler, R. L. Scoring rules for con- 183–186. Amer. Statist. Assoc., 1996.
tinuous probability distributions. Manag. Sci., 22(10):
1087–1096, 1976. SciPy Developers. SciPy Project, 2026. URL https:
//[Link]/doc/scipy/reference/
Meeker, W. Q., Hahn, G. J., and Escobar, L. A. Statistical generated/[Link].
Intervals: A Guide for Practitioners and Researchers. SciPy API Reference, accessed 2026-01.
John Wiley & Sons, Hoboken, New Jersey, second edition,
2017. ISBN 978-0-471-68717-7. Scott, D. W. Multivariate Density Estimation: Theory, Prac-
tice, and Visualization. John Wiley & Sons, New York,
Mockus, J., Tiesis, V., and Zilinskas, A. The application of 1992.
Bayesian methods for seeking the extremum. In Dixon,
L. C. W. and Szegö, G. P. (eds.), Towards Global Optimi- Shimodaira, H. Improving predictive inference under co-
sation, volume 2, pp. 117–129. North-Holland, Amster- variate shift by weighting the log-likelihood function. J.
dam, 1978. Statist. Plann. Inference, 90(2):227–244, 2000.

Nadarajah, S. A generalized normal distribution. J. Appl. Srinivas, N., Krause, A., Kakade, S. M., and Seeger, M.
Stat., 32(7):685–694, 2005. Gaussian process optimization in the bandit setting: No

10
 Goal-Oriented Calibration of GP for BO

regret and experimental design. In Proc. 27th Interna-

tional Conference on Machine Learning (ICML 2010),
pp. 1015–1022, 2010.
Stanton, S., Maddox, W., and Wilson, A. G. Bayesian
optimization with conformal prediction sets. In Ruiz, F.,
Dy, J., and van de Meent, J. W. (eds.), Proceedings of The
26th International Conference on Artificial Intelligence
and Statistics, volume 206 of Proc. Mach. Learn. Res.,
pp. 959–986. PMLR, 25–27 Apr 2023.
Stein, M. L. Interpolation of Spatial Data: Some Theory for
Kriging. Springer Ser. Stat. Springer New York, 1999.

Surjanovic, S. and Bingham, D. Virtual library of simula-

tion experiments: Test functions and datasets, 2013. URL
[Link] Accessed
November 2025.
Tom, G., Lo, S., Corapi, S., Aspuru-Guzik, A., and Sanchez-
Lengeling, B. Ranking over regression for Bayesian opti-
mization and molecule selection. APL Machine Learning,
3(3):036113, 08 2025.
Tuo, R. and Wang, W. Uncertainty quantification for
Bayesian optimization. In Camps-Valls, G., Ruiz, F. J. R.,
and Valera, I. (eds.), Proceedings of the 25th Interna-
tional Conference on Artificial Intelligence and Statistics
(AISTATS), volume 151 of Proc. Mach. Learn. Res., pp.
2862–2884. PMLR, 2022.
Vazquez, E. Gpmp: the Gaussian process micro package,
2026. URL [Link]
gpmp.
Vazquez, E. and Bect, J. Convergence properties of the
expected improvement algorithm with fixed mean and
covariance functions. J. Statist. Plann. Inference, 140:
3088–3095, 11 2010.

Villemonteix, J., Vazquez, E., and Walter, E. An informa-

tional approach to the global optimization of expensive-
to-evaluate functions. J. Global Optim., 44(4):509–534,
2009.
Vovk, V., Gammerman, A., and Shafer, G. Algorithmic
Learning in a Random World. Springer, 2005.
Vovk, V., Shen, J., Manokhin, V., and Xie, M. Nonparamet-
ric predictive distributions based on conformal prediction.
Mach. Learn., 108(3):445–474, 2019.
Zhang, Y. and Candès, E. J. Posterior conformal predic-
tion, 2024. URL [Link]
19712.

11
 Goal-Oriented Calibration of GP for BO

Appendices
The appendix is organized as follows. Appendix A recalls the generalized normal distribution. Appendix B gathers the
proofs of the main results. Appendix C briefly reviews scoring rules. Appendix D presents the consistency statement
discussed in Section 4.3. Appendix E provides implementation details, and Appendix F reports additional experimental
results.

A. Generalized Normal Distribution

The generalized normal distribution GN (β, l, λ) is a three-parameter family on R, with shape parameter β > 0, location
parameter l ∈ R, and scale parameter λ > 0. Its density is
β !
|z − l|

β
g(z) = exp − , z ∈ R. (28)
2 Γ(1/β) λ λ

It is symmetric about l. The mean is E[Z] = l and the variance is

Γ(3/β)
var(Z) = λ2 . (29)
Γ(1/β)
The parameter β controls tail decay: β = 2 yields a Gaussian distribution, β < 2 corresponds to heavier tails, and β > 2 to
lighter tails. In particular,  2
λ
GN (2, l, λ) = N l, , (30)
2
since Γ(3/2)/Γ(1/2) = 1/2.
The cumulative distribution function can be written using the lower incomplete gamma function
Z x
γ(a, x) = ta−1 e−t dt. (31)
0

For z ∈ R, the CDF Θβ,l,λ of GN (β, l, λ) may be written as

 
l−z β

 21 − 2 Γ(1/β)
1
γ β1 ,


λ , z ≤ l,
Θβ,l,λ (z) = P(Z ≤ z) = 
 (32)
1 + 1 1 z−l β
2 2 Γ(1/β) γ β, λ , z ≥ l.

We also write Θβ := Θβ,0,1 .

Further properties are given in Nadarajah (2005). The distribution GN (β, l, λ) is available in SciPy as
[Link]; see SciPy Developers (2026).

B. Proofs
B.1. Proof of Proposition 5.1
Proof. Let β > 0 and λ > 0. Write Z = l + λT with T ∼ GN (β, 0, 1), whose density is
β
gβ (t) = exp(−|t|β ). (33)
2 Γ(1/β)
Set z = a − l and s = z/λ. Then
     
Z s
E (a − Z)+ = E (z − λT )+ = λ E (s − T )+ = λ (s − t) gβ (t) dt. (34)
−∞

Therefore,
 
 Z s 
E (a − Z)+ = λ s Θβ (s) − t gβ (t) dt . (35)
−∞

12
 Goal-Oriented Calibration of GP for BO
Rs
It remains to compute −∞
t gβ (t) dt. Since t 7→ t exp(−|t|β ) is odd, for any s ∈ R,
Z s Z ∞
t exp(−|t|β ) dt = − u exp(−uβ ) du. (36)
−∞ |s|

By the change of variables y = uβ ,

Z ∞
1 ∞ 2/β−1
 
1 2
Z
u exp(−uβ ) du = y exp(−y) dy = Γ , |s|β . (37)
|s| β |s|β β β

Hence,
s s  
1 2
Z Z
β
t gβ (t) dt = t exp(−|t|β ) dt = − Γ , |s|β . (38)
−∞ 2 Γ(1/β) −∞ 2 Γ(1/β) β
Substituting back yields
    
1 2 z λ 2 z β
, |s|β
 
E (a − Z)+ = λ s Θβ (s) + Γ = z Θβ + Γ , . (39)
2 Γ(1/β) β λ 2 Γ(1/β) β λ
This is (26).
For λ = 0, Z = l almost surely, so E[(a − Z)+ ] = max(a − l, 0), which defines γ(z, 0, β). Continuity at λ = 0 follows
by considering three cases. If z > 0, then Θβ (z/λ) → 1 and the incomplete-gamma term vanishes. If z < 0, then
Θβ (z/λ) → 0 and the same incomplete-gamma term vanishes. If z = 0, the first term is identically zero and the second
term is λΓ(2/β)/(2Γ(1/β)), hence converges to zero. Strict positivity for λ > 0 holds because Z has a continuous density
that is strictly positive on R, hence P(Z < a) > 0 for any finite a, and (a − Z)+ > 0 on {Z < a}.

B.2. Connection between EI in the Gaussian Case and the Generalized Normal Model
For β = 2, the generalized normal distribution coincides with a Gaussian:
 2
λ
GN (2, l, λ) = N l, . (40)
2
√ √
Let σ = λ/ 2 and τ = (a − l)/σ. For β = 2 we have Θ2 (s) = Φ( 2 s), and Γ(1, x) = exp(−x). Substituting into (26)
yields
γ(a − l, λ, 2) = (a − l)Φ(τ ) + σ φ(τ ), (41)
where Φ and φ denote the CDF and PDF of N (0, 1). This is the usual EI formula for a Gaussian predictive distribution with
mean l and standard deviation σ.

B.3. Proof of Proposition 5.3

Recall that, at iteration n, EI with TC GP is


ρn (x) = γ mn − fn (x), λn σn (x), βn , mn = min f (Xi ), (42)
1≤i≤n

and we set
νn (f ) = max ρn (x) = ρn (Xn+1 ), Xn+1 ∈ arg max ρn (x), (43)
x∈X x∈X

with (βn , λn ) ∈ [β0 , β1 ] × [λ0 , λ1 ] and λ0 > 0 for all sufficiently large n.
Lemma B.1. Let 0 < β0 < β1 . The map (z, s, β) 7→ γ(z, s, β), with γ(z, 0, β) = z+ , is continuous on compact subsets of
R × [0, ∞) × [β0 , β1 ].

Proof. For s > 0, continuity follows from the closed-form expression (26). It remains to check continuity at s = 0,
uniformly for β ∈ [β0 , β1 ]. Let Tβ ∼ GN (β, 0, 1). Since the map x 7→ x+ is Lipschitz,
 
|γ(z, s, β) − z+ | = E (z − sTβ )+ − z+ ≤ s E|Tβ |. (44)

13
 Goal-Oriented Calibration of GP for BO

Moreover,
Γ(2/β)
E|Tβ | = (45)
Γ(1/β)
is bounded on [β0 , β1 ] by continuity. Hence γ(z, s, β) → z+ as s ↓ 0, uniformly for β ∈ [β0 , β1 ]. If also z ′ → z, then
′
|z+ − z+ | → 0, so continuity holds at s = 0 as well. This proves the claim.
Lemma B.2. For all f ∈ H, lim inf n→∞ νn (f ) = 0.

Proof. Fix f ∈ H. Since X is compact, (Xn ) has a cluster point x̃ ∈ X and there exists a strictly increasing map ϕ : N → N
such that Xϕ(n) → x̃. Since k is continuous and X is compact, every f ∈ H is continuous on X, hence f (Xϕ(n) ) → f (x̃).
By Vazquez & Bect (2010, Prop. 10),

fϕ(n)−1 (Xϕ(n) ) → f (x̃), σϕ(n)−1 (Xϕ(n) ) → 0. (46)

Moreover, since ϕ is increasing,

mϕ(n)−1 ≤ mϕ(n−1) ≤ f (Xϕ(n−1) ) → f (x̃), (47)
hence lim supn→∞ mϕ(n)−1 ≤ f (x̃). Set

zϕ(n)−1 := mϕ(n)−1 − fϕ(n)−1 (Xϕ(n) ). (48)

Then lim supn→∞ zϕ(n)−1 ≤ 0. Define also

tϕ(n)−1 := λϕ(n)−1 σϕ(n)−1 (Xϕ(n) ), (49)

so that tϕ(n)−1 → 0 since 0 < λϕ(n)−1 ≤ λ1 .

Fix ε > 0. For n large enough, zϕ(n)−1 ≤ ε. For λ > 0, z 7→ γ(z, λ, β) = E[(z − Z)+ ] is nondecreasing, hence

νϕ(n)−1 (f ) = ρϕ(n)−1 (Xϕ(n) ) (50)

= γ zϕ(n)−1 , tϕ(n)−1 , βϕ(n)−1 (51)


≤ γ ε, tϕ(n)−1 , βϕ(n)−1 . (52)

Since tϕ(n)−1 → 0 and, for all large n, βϕ(n)−1 ∈ [β0 , β1 ], every subsequence of γ(ε, tϕ(n)−1 , βϕ(n)−1 ) has a further
subsequence along which βϕ(n)−1 → β̄ ∈ [β0 , β1 ]. Along this further subsequence, Lemma B.1 gives

γ(ε, tϕ(n)−1 , βϕ(n)−1 ) −→ γ(ε, 0, β̄) = ε. (53)

Thus the full sequence converges to ε. Therefore,

lim sup νϕ(n)−1 (f ) ≤ ε. (54)

n→∞

Since ε > 0 is arbitrary, lim inf n→∞ νn (f ) = 0.

Proof of Proposition 5.3. Assume by contradiction that (Xn ) is not dense in X. Then there exists x̃ ∈ X such that x̃ is not
adherent to {Xn : n ≥ 1}. By the NEB property, σn2 (x̃) does not converge to 0. Since (σn2 (x̃))n≥1 is nonincreasing, it
converges to some c2 > 0, hence inf n≥1 σn (x̃) ≥ c.
Since fn is the H-projection of f onto the span of kernel sections, ∥fn ∥H ≤ ∥f ∥H , hence
p p
|fn (x̃)| ≤ ∥fn ∥H k(x̃, x̃) ≤ ∥f ∥H k(x̃, x̃). (55)

Also, (mn ) is bounded because f is continuous

p on the compact set X. Therefore, zn∗ := mn − fn (x̃) ranges in a compact
interval I ⊂ R and σn (x̃) ranges in [c, k(x̃, x̃)].
Define
ψ(z, σ, β, λ) = γ(z, λσ, β). (56)

14
 Goal-Oriented Calibration of GP for BO

The map ψ is continuous by Lemma B.1. Consider the compact set

p
K = I × [c, k(x̃, x̃)] × [β0 , β1 ] × [λ0 , λ1 ]. (57)

Since λ0 > 0 and γ(·, λ, β) is strictly positive for λ > 0, we have ψ > 0 on K, hence

η := inf ψ(z, σ, β, λ) > 0. (58)

(z,σ,β,λ)∈K

For all sufficiently large n,

ρn (x̃) = γ(zn∗ , λn σn (x̃), βn ) ≥ η. (59)
Since Xn+1 maximizes ρn , we obtain

νn (f ) = ρn (Xn+1 ) ≥ ρn (x̃) ≥ η for all sufficiently large n, (60)

which contradicts Lemma B.2. Therefore, (Xn ) is dense in X.

C. Scoring Rules
C.1. General Definition
To assess both sharpness and calibration, we use proper scoring rules (Gneiting & Raftery, 2007). Scoring rules assign a
numerical score S(F̂ , z) to a predictive forecast F̂ and an outcome z ∈ R. A scoring rule is strictly proper if, for the true
distribution F ,
EZ∼F [S(F, Z)] ≤ EZ∼F [S(F̂ , Z)], (61)
with equality if and only if F̂ = F .
For goal-oriented prediction, Matheson & Winkler (1976) introduces the weighted continuous ranked probability score
(twCRPS), defined for a weight function w, for instance w(u) = 1{u ≤ t}, by
Z
StwCRPS (F̂ , z) = (F̂ (u) − 1{u ≥ z})2 w(u) du. (62)

For nonconstant weights, twCRPS is a proper weighted score and does not identify the full predictive distribution outside
the weighted region.
Given a true CDF F , a forecasted CDF F̂ , and a scoring rule S the predictive performance is evaluated using the expected
score
JS (F̂ ) = EZ∼F (S(F̂ , Z)). (63)

An empirical version is given below. In the spatial setting the expected score of the family of predicted forecast CDFs
F̂n (· | x), indexed by x ∈ X is given by,
h i
JS,µ (F̂n ) = En S(F̂n (· | X), f (X)) . (64)

C.2. Empirical twCRPS

Let weight function w : R → R+ and F̂ a forecasted CDF. In the following consider the chaining function v defined as any
R x′
function that verifies v(x′ ) − v(x) = x w(u)du (see, e.g., Allen et al., 2023). The twCRPS at z ∈ R can be rewritten as

 1 
StwCRPS (F̂ , z) = EF̂ |v(Y ) − v(z)| − EF̂ |v(Y ) − v(Y ′ )| ,
 
(65)
2

where Y and Y ′ are independent random variables with CDF F̂ .

For a threshold t ∈ R, let w(z) = 1{z ≤ t}, then v(z) = min(z, t).

15
 Goal-Oriented Calibration of GP for BO

For the empirical forecast-observation pair (Z, F̂ ) and a sample (Yℓ )K

ℓ=1 from F̂ , the twCRPS at Z can be approximated by
the Monte Carlo estimator
K K K
1 X 1 XX
SbtwCRPS,K (F̂ , Z) = |v(Yℓ ) − v(Z)| − |v(Yℓ ) − v(Yr )|. (66)
K 2K 2 r=1
ℓ=1 ℓ=1

For multiple forecast-observation pairs (Zi , F̂i )m

i=1 , the empirical score is

m
1 Xb
JbtwCRPS,m = StwCRPS,K (F̂i , Zi ). (67)
m i=1

D. Impact of Calibration at the Current Best Value

This appendix provides the consistency statement mentioned in Section 4.3. It shows that, if one could enforce occurrence
calibration at the current best value mn , then the predicted probability mass below mn must collapse as mn approaches the
global minimum.
a.s.
Proposition D.1 (Vanishing tail mass with occurrence calibration at mn ). Assume that mn −→ f (x∗ ) for some global
minimizer x∗ ∈ arg minx∈X f (x). If F̂n is occurrence-calibrated at threshold mn , that is,
 

En F̂n (mn | X) = Pn f (X) ≤ mn , X ∼ µ, (68)

then for every ε ∈ (0, 1], we have

lim sup Pn F̂n (mn | X) > ε ≤ ε−1 µ ({x ∈ X : f (x) = f (x∗ )})


a.s. (69)
n→∞

In particular, if µ ({x ∈ X : f (x) = f (x∗ )}) = 0, then

Pn F̂n (mn | X) > ε → 0 a.s. (70)

Proof. Fix ε ∈ (0, 1] and set An := {F̂n (mn | X) > ε}. By Markov’s inequality,

Pn (An ) ≤ ε−1 En F̂n (mn | X) .

 

Using (68), this yields

Pn F̂n (mn | X) > ε ≤ ε−1 Pn f (X) ≤ mn .



(71)

On any sample path such that mn → f (x∗ ), we have f (x) ≥ f (x∗ ) for all x ∈ X, and therefore

1{f (x) ≤ mn } → 1{f (x) = f (x∗ )} for all x ∈ X. (72)

Moreover, 0 ≤ 1{f (x) ≤ mn } ≤ 1, so dominated convergence with respect to µ gives

Pn f (X) ≤ mn = En 1{f (X) ≤ mn } → En 1{f (X) = f (x∗ )} = µ ({x ∈ X : f (x) = f (x∗ )})

   
a.s. (73)

Combining (71) with (73) and taking lim sup yields (69). The last statement follows when

µ ({x ∈ X : f (x) = f (x∗ )}) = 0. (74)

Interpretation. Proposition D.1 is a population statement with an auxiliary draw X ∼ µ (independent of Dn ). It implies
that, when the set of global minimizers has µ-measure 0, enforcing occurrence calibration at mn forces F̂n (mn | X) to be
small for µ-almost all X at large n.
The result does not describe how the predictive lower tail behaves in the region explored by the BO policy. This motivates
calibrating below a higher threshold tn and using notions based on the truncated CDF F̂n,tn (· | x).

16
 Goal-Oriented Calibration of GP for BO

E. Implementation Details
E.1. GP Matérn Model and Parameters
In the experiment reported Section 6, the GP model GP(m, k) is such that the mean function is constant over X, m(x) = µ
for x ∈ X. The kernel k is selected from the anisotropic Matérn kernel family. The anisotropic Matérn kernel family is
defined, for x, y ∈ Rd , as
d
2 2
X (x[i] − y[i] )2
kσ,ν,ρ (x, y) = σ κν (h), h = , (75)
i=1
ρ2i

where σ 2 is the variance parameter, ρ = (ρ1 , . . . , ρd ) are component-wise lengthscales, and κν denotes the Matérn
correlation function (Stein, 1999, Chapter 2.7). The smoothness parameter is fixed to ν = p + 1/2 with p ∈ N⋆ .
The parameter (µ, ρ, σ) is selected by maximum likelihood. In the Experiments in Section 6, p = 2.

E.2. Maximization of EI
EI is maximized with the sequential Monte Carlo (SMC) procedure of Feliot et al. (2016), available in gpmp. A population
of particles is iteratively reweighted with EI-proportional weights and propagated through mutation steps, yielding a discrete
approximation of the EI maximizers. The best particle is then used to initialize a local maximization with the sequential
least squares quadratic programming (SLSQP) algorithm. This hybrid global–local procedure is applied for all methods. In
the Experiments in Section 6 we used 1000 particles.

E.3. Bayesian Optimization with TC GP

Algorithm 1 summarizes EI-based BO with TC GP. It follows the procedure described in the main text. At iteration n, the
empirical quantile qδ,n is used as a candidate threshold when there is sufficient data support in the lower tail, and the previous
threshold is kept otherwise. Concretely, we estimate the excursion probability below qδ,n and enforce p̂qwδ,n ,n ≥ pmin . If the
condition fails, the threshold is frozen to avoid calibrating on too few points. Since the BO design is adaptive, this support
rule is only a safeguard.

Density-ratio estimation and weights. The weighted LOO quantities in Section 4.2 aim at approximating µ–averages
using the adaptively sampled BO locations (Xi )ni=1 , which are not distributed according to µ. We use an empirical density-
ratio correction as a heuristic for this adaptive design. If νn denotes the smoothed empirical design density of the BO
locations, the idealized weights have the form
dµ
wi ∝ (Xi ), (76)
dνn
and are used to approximate, for the integrands φ appearing in the weighted LOO criteria,
Z n
X wi
φ(x) µ(dx) ≈ w̃i φ(Xi ), w̃i = Pn . (77)
X i=1 j=1 wj

In our experiments, µ is uniform on X, hence

1
wi ∝ . (78)
νn (Xi )

In practice, we use a KDE estimate ν̂n of this design density. Concretely, we rescale the design points Xi to [0, 1]d using
the bounds of X, and fit a Gaussian kernel density estimator (KDE) to the rescaled points (Scott, 1992). The bandwidth is
chosen according to Scott’s rule (Scott, 1992). We then set

dµ
wi ∝ (Xi ), (79)
dν̂n
which reduces to wi ∝ 1/ν̂n (Xi ) for uniform µ, and normalize to obtain (w̃i )ni=1 .
Large density ratios can lead to highly variable weights. A standard safeguard is weight clipping, for instance wi ←
min(wi , wmax ) for a fixed wmax > 0, followed by renormalization. We did not use clipping in the reported experiments.

17
 Goal-Oriented Calibration of GP for BO

Algorithm 1 BO with TC GP
Require: Initial dataset Dn0 = {(Xi , Zi )}ni=1
0
with Zi = f (Xi ), design space X, budget Nmax , quantile level δ ∈ (0, 1],
minimum excursion probability pmin > 0.
Ensure: Best value mNmax and associated point Xmin .
1: Set n ← n0 .
2: Set mn ← min1≤i≤n Zi and choose any Xmin ∈ arg min1≤i≤n Zi .
3: Set the current threshold t ← qδ,n .
4: while n < Nmax do
5: Fit a GP model to Dn by maximum likelihood and obtain functions x 7→ fn (x) and x 7→ σn (x).
6: Estimate the importance weights (w̃i )ni=1 (Section 4.2 and paragraph below).
7: Set t̃ ← qδ,n .
8: Compute the weighted excursion frequency
n
X
w
p̂t̃,n = w̃i 1{Zi ≤ t̃}.
i=1

w
9: if p̂t̃,n ≥ pmin then
10: Update t ← t̃.
11: end if
LOO,w
12: Select (βn , λn ) by minimizing Jt,n on [β0 , β1 ] × [λ0 , λ1 ] (Section 5.1 and paragraph below).
13: Define EI ρn using the TC GP predictive CDFs with parameters (βn , λn ) (Section 5.2).
14: Choose Xn+1 ∈ arg maxx∈X ρn (x) (paragraph below).
15: Evaluate Zn+1 = f (Xn+1 ) and set Dn+1 ← Dn ∪ {(Xn+1 , Zn+1 )}.
16: if Zn+1 < mn then
17: Set mn+1 ← Zn+1 and Xmin ← Xn+1 .
18: else
19: Set mn+1 ← mn .
20: end if
21: Set n ← n + 1.
22: end while
23: Return (Xmin , mNmax ).

Selection of (βn , λn ). We minimize JtLOO,w

n ,n
with a two-stage procedure. First, we sample (β, λ) uniformly on [β0 , β1 ] ×
[λ0 , λ1 ] (900 candidates) and select the best pair (β̃, λ̃) according to JtLOO,w
n ,n
. Second, we refine with the Nelder–Mead
algorithm (Nelder & Mead, 1965) initialized at (β̃, λ̃). The supremum over u ∈ [0, 1] in JtLOO,w
n ,n
is evaluated on a fixed grid
of u values.

E.4. Test-Set Evaluation of µ–Calibration Metrics

When an independent evaluation set is available, µ–calibration metrics can be computed for X ∼ µ without leave-one-out
or reweighting. Let (Xj′ , Zj′ )L ′ ′ ′
j=1 with Xj ∼ µ i.i.d. and Zj = f (Xj ), and assume that the predictive CDFs F̂n (· | x) are
constructed from Dn only.

Direct test-set estimators. Fix t ∈ R. Define

(
F̂n (Zj′ | Xj′ ) F̂n (t | Xj′ ), Zj′ ≤ t,

′
Ut,j = (80)
1, Zj′ > t,
PL
with the convention that F̂n (t | Xj′ ) > 0 when Zj′ ≤ t. Assume that j=1 1{Zj′ ≤ t} > 0 and set
PL
j=1 1{Zj′ ≤ t} 1{Ut,j
′
≤ u}
G′t,L (u) = PL , ′
JtKS-PIT,L (F̂n | t) = sup G′t,L (u) − u . (81)
′
j=1 1{Zj ≤ t} u∈[0,1]

18
 Goal-Oriented Calibration of GP for BO

The excursion probability pt = P(f (X) ≤ t) is estimated by

L
1X
p̂t,L = 1{Zj′ ≤ t}, (82)
L j=1

and the corresponding occurrence discrepancy estimator is

L
′ 1X
rt,L = p̂t,L − F̂n (t | Xj′ ) . (83)
L j=1

SMC-based estimators for small thresholds. For small t (in particular t = mn at later BO iterations), direct sampling
from µ yields too few points with f (X) ≤ t to estimate Gµ,t and pt accurately. We therefore use subset simulation
implemented via sequential Monte Carlo (SMC) (Bect et al., 2017) to obtain particles approximately distributed according
to µ(· | f (X) ≤ t) and an estimate of pt = P(f (X) ≤ t).
At the initial iteration with n0 observations, we run subset simulation down to t = mn0 to compute an estimate p̂mn0 and
(0) (0)
a particle cloud (X̃i )m m
i=1 ≈ µ(· | f (X) ≤ mn0 ). We evaluate tKS–PIT at mn0 using (X̃i )i=1 , and we evaluate the
occurrence discrepancy using p̂mn0 .
When BO produces a new best value mn < mn−1 , we do not rerun subset simulation from scratch. We reuse the
particles available at level mn−1 as an initial population and perform additional SMC steps targeting the new excursion set
(n)
{x : f (x) ≤ mn }. This yields updated particles (X̃i )m i=1 ≈ µ(· | f (X) ≤ mn ) together with an updated estimate p̂mn .
The particle cloud is used to form the empirical distribution of the µ-tPIT below mn (and hence tKS–PIT), while p̂mn is
used in the occurrence discrepancy.
In all experiments, we use m = 1000 particles and Metropolis–Hastings moves for the Markov transition steps.

E.5. Reweighting for BCR GP

The BCR GP approach of Pion & Vazquez (2025) specifies a generalized normal model for standardized residuals and places
a uniform prior on θ = (β, λ). We recall the resulting posterior construction and describe how we incorporate importance
weights to target µ-probabilistic calibration over X. In our experiments, global calibration diagnostics for BCR GP are
reported using KS–PIT, i.e., the limit of tKS–PIT as t → ∞.
Let Rn,−i (Xi ) denote the leave-one-out standardized residual at Xi ,
Zi − fn,−i (Xi )
Rn,−i (Xi ) = , (84)
σn,−i (Xi )
where fn,−i and σn,−i are the GP predictive mean and standard deviation obtained from Dn \ {(Xi , Zi )}. Let p(r | θ) be
the density of GN (β, 0, λ) evaluated at r.
With a uniform prior on (0, a) × (0, b), BCR GP defines
n
Y
p (θ | {Rn,−i (Xi )}ni=1 ) ∝ p(Rn,−i (Xi ) | θ) 1{0 < β < a} 1{0 < λ < b}. (85)
i=1

Let ν̂n be the density estimate used above. To incorporate importance weights, we replace the standard likelihood by
a weighted likelihood, following Shimodaira (2000). Let (ωi )ni=1 be nonnegative weights proportional to dµ/dν̂n (Xi )
(Section 4.2). For the weighted likelihood, the global scale of the exponents affects posterior concentration, so we fix the
normalization X n Xn
w̄i = n ωi ωj , w̄i = n. (86)
j=1 i=1

The resulting pseudo-posterior is

n
Y w̄i
pw (θ | {Rn,−i (Xi )}ni=1 ) ∝ p(Rn,−i (Xi ) | θ) 1{0 < β < a} 1{0 < λ < b}. (87)
i=1

19
 Goal-Oriented Calibration of GP for BO

∗ ∗
Given Monte Carlo samples (βi , λi )m i=1 , Pion & Vazquez (2025) derive a Bayesian selection rule for (β , λ ) targeting
probabilistic calibration, quantified by the KS–PIT up to a prescribed quantile. To limit computational cost, we set m = 200
in this work; increasing m further had negligible impact on BO performance.

E.6. Experimental Parameter Settings

We summarize below the main parameters used in all numerical experiments.

Test-set evaluation parameters

• Test-set evaluation of calibration metrics: The subset simulation procedure uses k0 = 1000 particles to estimate
P(f (X) ≤ mn ) and to generate samples from µ(· | f (X) ≤ mn ).
• Occurrence discrepancy: The occurrence discrepancy rt is computed using an additional k1 = 1000 points sampled
uniformly over X.
• tKS–PIT computation: To compute the tKS–PIT, we subsample k2 = 900 particles from the subset simulation output,
ensuring that all corresponding function values lie below the current threshold mn at each iteration.

Model parameters

• GP prior regularity: We use a Matérn covariance kernel with smoothness parameter ν = p + 12 , with p = 2,
corresponding to twice mean-square differentiable sample paths. All GP hyperparameters are estimated by maximum
likelihood.
• Parameter bounds for calibration: The calibration parameters are constrained to

λ ∈ [λ0 , λ1 ] = [5 × 10−3 , 10], β ∈ [β0 , β1 ] = [0.1, 10].

• pmin = 0.015.

E.7. Test Functions Used in the Experiments

Name Domain Expression of f (x)

2
Goldstein–Price [−2, 2] (1 + (x1 + x2 + 1)2 (19 − 14x1 + 3x21 − 14x2 + 6x1 x2 + 3x22 ))(30 + (2x1 − 3x2 )2 (18 −
32x1 + 12x21 + 48x2 − 36x1 x2 + 27x 2
 2 ))
d Pd−1 2 2 2
Rosenbrock [−5, 10] i=1 100(x
 i+1q− xP i ) +(xi − 1)
d
[−32.768, 32.768] −20 exp −0.2 d i x2i − exp d1 i cos(2πxi ) + 20 + e
1
P

Ackley
Dixon–Price–d [−10, 10]d (x1 − 1)2 + di=2 i(2x2i − xi−1 )2
P
(d) (d) (d)

− i=1 ci exp − dj=1 aij (xj − pij )2
d P4 P
Hartmann–d [0, 1]
Michalewicz–d [0, π]d − di=1 sin(xi ) sin2m (ix2i /π)
P
  √ 2 2  0.1
x1 +x2
Cross-In-Tray [−10, 10]2 −0.0001 sin(x1 ) sin(x2 ) exp 100 − π
+1
P −1
[0, 10]4 − m 4 2
P
Shekel–m i=1 j=1 (xj − Cji ) + βi
Pd Pd  2
Perm–d [−d, d]d i=1
i
j=1 (j + βP ) xj − j i
1

Table 1. Test functions used in the numerical experiments. Each function is continuous, deterministic, and defined on the domain X.
Expressions and domains follow the standard definitions provided in Surjanovic & Bingham (2013), with parameter values specified
below.

Hartmann parameters.
3.0 10.0 30.0 1.0 3689 1170 2673
     
(3) 0.1 10.0 35.0 (3) 1.2 (3) −4 4699 4387 7470
a = , c =  , p = 10 .
3.0 10.0 30.0 3.0 1091 8732 5547
0.1 10.0 35.0 3.2 381 5743 8828

20
 Goal-Oriented Calibration of GP for BO

10 3 17 3.5 1.7 8 1.0 1312 1696 5569 124 8283 5886

     
(6) 0.05 10 17 0.1 8 14 (6) 1.2 (6) −4 2329 4135 8307 3736 1004 9991
a = , c =  , p = 10 .
3 3.5 1.7 10 17 8 3.0 2348 1451 3522 2883 3047 6650
17 8 0.05 10 0.1 14 3.2 4047 8828 8732 5743 1091 381

Michalewicz parameter.
m = 10

Perm parameter.
βP = 1

Shekel parameters.
1 T
β= (1, 2, 2, 4, 4, 6, 3, 7, 5, 5)
10

4 1 8 6 3 2 5 8 6 7
 
4 1 8 6 7 9 3 1 2 3.6
C= .
4 1 8 6 3 2 5 8 6 7 
4 1 8 6 7 9 3 1 2 3.6

F. Additional Experiments
F.1. Running Time of the Methods
In this section, we report the per-iteration time required to train the GP or calibrated model (excluding sampling criteria
maximization). All calibration methods have a similar computation time relatively to GP, but RE GP is more costly (by a
factor 10).

102 102
102 102

101
101
101
101
0
10

100 100
10−1 100

0 50 100 150 200 0 50 100 150 200 0 50 100 150 200 0 50 100 150 200
Goldstein–Price Rosenbrock6 Ackley4 Dixon–Price4
GP reGP δ = 0.25 tcGP δ = 0.05
onGP bcrGP

Figure 4. Training time for RE GP, GP, TC GP, onGP, BCR GP at each iteration of BO.

F.2. Alternative Selection Criteria for TC GP

We compare several objectives for selecting (β, λ) in TC GP, with the aim of controlling thresholded µ-calibration and
occurrence calibration below t = qδ,n . Rule 0 is the criterion J defined in (18). Figure 5 reports BO results over 100
iterations with δ = 0.1.
A natural variant of (18) replaces the target u 7→ u κβ,λ
t by the truncated map
 
u 7→ min u κβ,λ t , 1 . (88)

When κβ,λ
t > 1, truncation saturates the target at 1. Large values of κβ,λ
t then become weakly penalized, and the objective
can become insensitive to overestimation of lower-tail mass. We do not use this variant.

Rule 1. Multiply the weighted LOO tKS–PIT metric by the occurrence discrepancy:
 
LOO,w
J (1) (β, λ) = JtKS-PIT,n F̂nβ,λ | t rt . (89)

21
 Goal-Oriented Calibration of GP for BO

Excursion Prob. Proportion Excursion Prob. Proportion

100
10−1
10−2
10−3
10−4
10−5
10−6
10−7
10−8 10−9

10−10 10−11

10−12 10−13

1 20 40 60 80 100 1 0.2 0.4 0.6 0.8 1 1 20 40 60 80 100 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion

rule 0, J rule 1, J (1) rule 2, J (2) rule 3 J (3) rule 0, J rule 1, J (1) rule 2, J (2) rule 3 J (3)

(a) Ackley with d = 4 (b) Dixon–Price with d = 4

Excursion Prob. Proportion Excursion Prob. Proportion
100
10−1
10−2
10−3
10−4
10−5
10−6
10−7
10−8
10−9
10−10
10−11
10−12 10−13

1 20 40 60 80 100 1 0.2 0.4 0.6 0.8 1 1 20 40 60 80 100 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion

rule 0, J rule 1, J (1) rule 2, J (2) rule 3 J (3) rule 0, J rule 1, J (1) rule 2, J (2) rule 3 J (3)

(c) Goldstein–Price (d) Rosenbrock with d = 6

Figure 5. Comparison of TC GP variants. Left: median and 10%/90% quantiles of pmn = P(f (X) ≤ mn ), where mn is the best
observed value after n evaluations and X ∼ µ. Right: fraction of runs reaching the prescribed target value.

Rule 2. Following Allen et al. (2025), define

 
Pn Utβ,λ ≤ u | f (X) ≤ t
ct (u) = . (90)
κβ,λ
t

Then use the KS-type deviation

J (2) (β, λ) = sup |ct (u) − u| . (91)
u∈[0,1]

Rule 3. Compare the conditional tPIT curve to the rescaled map u 7→ u/κβ,λ
t :
  u
J (3) (β, λ) = sup Pn Utβ,λ ≤ u | f (X) ≤ t − β,λ . (92)
u∈[0,1] κt

Choice of objective. Figure 5 reports BO performance for the different selection rules. Performance is summarized by
pmn = P(f (X) ≤ mn ), where mn is the best value observed after n evaluations and X ∼ µ (with µ uniform), and by the
fraction of runs reaching the prescribed target.
Empirically, the rules yield similar results overall. Rule 0 is used in the main experiments because it is the most consistent
across the considered test functions. In particular, it behaves well on both highly multimodal landscapes (e.g., Ackley with
d = 4) and on difficult smooth landscapes (e.g., Goldstein–Price). The alternative rules can either overweight tail-rank
calibration, which can degrade optimization on smoother problems, or penalize overly optimistic lower-tail mass too weakly
on multimodal problems. Rule 0 keeps a direct penalty on deviations of κβ,λ t from 1 without additional tuning.

F.3. Effect of the Choice of δ for RE GP on BO Performance

Figure 6 studies the sensitivity of RE GP to the threshold parameter δ. BO performance is summarized by pmn = P(f (X) ≤
mn ), where mn is the best value observed after n evaluations and X ∼ µ (with µ uniform), and by the fraction of runs

22
 Goal-Oriented Calibration of GP for BO

Excursion Prob. Proportion Excursion Prob. Proportion

100 100

10−2 10−2

10−4 10−4

−6
10 10−6

10−8
10−8
−10
10
10−10

1 20 40 60 80 100 1 0.2 0.4 0.6 0.8 1 1 20 40 60 80 100 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion

reGP δ = 0.05 reGP δ = 0.1 reGP δ = 0.25 reGP δ = 0.5 reGP δ = 0.05 reGP δ = 0.1 reGP δ = 0.25 reGP δ = 0.5

(a) Ackley with d = 4 (b) Dixon–Price with d = 4

Excursion Prob. Proportion Excursion Prob. Proportion
100
−1
10
10−2
10−3
10−4
10−5

10−7 10−6

10 −9 10−8

10 −11 10−10

10−13 10−12

1 20 40 60 80 100 1 0.2 0.4 0.6 0.8 1 1 20 40 60 80 100 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion

reGP δ = 0.05 reGP δ = 0.1 reGP δ = 0.25 reGP δ = 0.5 reGP δ = 0.05 reGP δ = 0.1 reGP δ = 0.25 reGP δ = 0.5

(c) Goldstein–Price (d) Rosenbrock with d = 6

Figure 6. Sensitivity of RE GP to δ (results shown after n = 100 BO iterations). Left: median and 10%/90% quantiles of pmn =
P(f (X) ≤ mn ), where mn is the best value observed after n evaluations and X ∼ µ. Right: fraction of runs reaching the prescribed
target level.

reaching the prescribed target.

Empirically, a large value of δ (e.g., δ = 0.5) consistently leads to degraded performance across the considered problems. In
contrast, δ = 0.25 achieves the best overall performance and is consistently comparable to, or better than, smaller values
such as δ = 0.05 and δ = 0.1. We therefore use δ = 0.25 for all RE GP experiments.

F.4. Effect of the Choice of δ for TC GP on BO Performance

Figure 7 studies the influence of the threshold parameter δ on TC GP. Performance is summarized by pmn = P(f (X) ≤ mn ),
where mn is the best value observed after n evaluations and X ∼ µ (with µ uniform), and by the fraction of runs reaching
the prescribed target.
Empirically, δ = 0.05 provides the most robust behavior across the considered problem classes. Larger values (e.g.,
δ = 0.25) can yield the best performance on difficult but convex functions such as Goldstein–Price, although δ = 0.05
remains competitive there. On highly multimodal, nonconvex problems such as Ackley, smaller values of δ perform
substantially better. We therefore use δ = 0.05 for all TC GP experiments.

F.5. Calibration on a Fixed Dataset

F.5.1. C ALIBRATION BELOW THE P RESCRIBED T HRESHOLD
We evaluate whether TC GP, either with the combined objective or with objectives targeting only thresholded or occurrence
calibration, improves calibration at the operational threshold tn = qn,δ used to select (β, λ), with δ ∈ {0.25, 0.1, 0.05}.
The GP setup matches Section 6, except that we work with a fixed training dataset Dn of size n = 30d obtained from
uniform sampling on X. For each test function, we generate 100 independent datasets and report averages over the replicates.
For a given δ, calibration and predictive quality below qn,δ are assessed using three metrics: twCRPS (jointly measuring

23
 Goal-Oriented Calibration of GP for BO

Excursion Prob. Proportion Excursion Prob. Proportion

100
10−1
10−2
10−3
10−4
10−5
10−6
10−7
10−8
10−9
10−10
10−11
10−12
10−13

1 20 40 60 80 100 1 0.2 0.4 0.6 0.8 1 1 20 40 60 80 100 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion

tcGP δ = 0.05 tcGP δ = 0.1 tcGP δ = 0.25 tcGP δ = 0.5 tcGP δ = 0.05 tcGP δ = 0.1 tcGP δ = 0.25 tcGP δ = 0.5

(a) Ackley with d = 4 (b) Dixon–Price with d = 4

Excursion Prob. Proportion Excursion Prob. Proportion
100
10−1
10−2
10−3
10−4
10−5
10−6
10−7
10−8

10−10 10−9

10 −12 10−11

1 20 40 60 80 100 1 0.2 0.4 0.6 0.8 1 1 20 40 60 80 100 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion

tcGP δ = 0.05 tcGP δ = 0.1 tcGP δ = 0.25 tcGP δ = 0.5 tcGP δ = 0.05 tcGP δ = 0.1 tcGP δ = 0.25 tcGP δ = 0.5

(c) Goldstein–Price (d) Rosenbrock with d = 6

Figure 7. Sensitivity of TC GP to δ (results shown after n = 100 BO iterations). Left: median and 10%/90% quantiles of pmn =
P(f (X) ≤ mn ), where mn is the best value observed after n evaluations and X ∼ µ. Right: fraction of runs reaching the prescribed
target level.

calibration and sharpness below qn,δ ), rt (occurrence calibration at tn ), and tKS–PIT (thresholded µ–calibration below tn ).
The tKS–PIT is computed using 4000 points sampled conditional on f (x) ≤ qn,δ by subset simulation (see Section 2 of
Bect et al., 2017). In addition, 4000 points are sampled uniformly on X to estimate rt and twCRPS.
Results are reported in Table 2. Overall, TC GP improves both rt and tKS–PIT for all values of δ, with smaller gains in
twCRPS. The local variant TC GP-thres consistently improves tKS–PIT, but can worsen rt on some problems (e.g., Hartman
with d = 6 and Ackley with d = 4). The marginal variant TC GP-occ improves rt but degrades tKS–PIT on several test
functions. For Ackley, TC GP shows no improvement over the GP baseline for δ ∈ {0.25, 0.1}, and only a limited gain for
δ = 0.05.

F.5.2. C OMPARISON BETWEEN M ULTIPLE M ETHODS

We evaluate whether RE GP and TC GP improve calibration below the current best value mn relative to a standard GP model.
As above, we consider a fixed training dataset Dn of size n = 30d obtained from uniform sampling on X. Results are
reported in Table 3. Overall, both methods achieve improved or comparable calibration across the test functions, with RE GP
showing the largest gains on smooth (nearly convex) objectives and TC GP on the more challenging functions.
On Ackley, we observe a marked improvement in both tKS–PIT and rt evaluated at mn relative to the other methods, while
improvements at the calibration threshold t = qn,δ (with δ = 0.05) were more limited. This suggests that calibration below
an operational threshold tn can still translate into substantial gains at more extreme levels t ≤ tn . A similar pattern is
observed on Hartman6.

F.6. Additional Results for EI and UCB

F.6.1. A DDITIONAL R ESULTS FOR EI
In this section, we report results on additional test functions in Figures 8 and 9. BO performance is summarized using two
complementary views: the attained-level probability pmn = P(f (X) ≤ mn ), where mn is the best observed value after n

24
 Goal-Oriented Calibration of GP for BO

Method Goldstein Price Rosenbrock6 Hartman6 Dixon Price 4 Ackley4

twCRPS rt tKS–PIT twCRPS rt tKS–PIT twCRPS rt tKS–PIT twCRPS rt tKS–PIT twCRPS rt tKS–PIT

n = 30, t = qn,δ , δ = 0.25

GP 6.9e2 0.04 0.64 6.0e3 0.01 0.16 0.049 0.095 0.15 2.7e2 0.0059 0.16 0.12 0.093 0.14
TC GP 5.7e2 0.036 0.27 6.0e3 0.01 0.14 0.052 0.06 0.28 2.6e2 0.0059 0.13 0.12 0.027 0.24
TC GP-occ 6.1e2 0.036 0.52 8.9e3 0.0095 0.36 0.059 0.046 0.79 4.3e2 0.0071 0.42 0.13 0.018 0.51
TC GP-thres 5.8e2 0.036 0.27 6e3 0.01 0.14 0.05 0.11 0.15 2.6e2 0.0059 0.13 0.12 0.078 0.13
n = 30, t = qn,δ , δ = 0.1
GP 6.0e2 0.11 0.83 2.8e3 0.011 0.3 0.027 0.018 0.27 1.1e2 0.0047 0.22 0.063 0.042 0.19
TC GP 5.1e2 0.065 0.48 2.8e3 0.011 0.27 0.027 0.021 0.24 1.0e2 0.0051 0.19 0.066 0.021 0.25
TC GP-occ 5.2e2 0.065 0.68 3.4e3 0.011 0.58 0.031 0.01 0.67 1.6e2 0.0057 0.44 0.071 0.012 0.52
TC GP-thres 5.1e2 0.065 0.41 2.9e3 0.011 0.27 0.031 0.096 0.12 1.0e2 0.0053 0.19 0.085 0.12 0.19
n = 30, t = qn,δ , δ = 0.05
GP 5.8e2 0.16 0.91 1.8e3 0.023 0.39 0.016 0.0088 0.39 56 0.005 0.32 0.042 0.014 0.29
TC GP δ = 0.05 5e2 0.089 0.65 1.8e3 0.02 0.34 0.016 0.016 0.18 55 0.0054 0.26 0.044 0.025 0.28
TC GP-occ 4.9e2 0.09 0.75 2.1e3 0.02 0.61 0.017 0.011 0.39 64 0.0056 0.48 0.045 0.011 0.42
TC GP-thres 4.9e2 0.089 0.62 1.8e3 0.02 0.31 0.023 0.09 0.15 56 0.0059 0.27 0.14 0.16 0.22
Table 2. Comparison of calibration using the twCRPS, rt , and tKS–PIT below t = qn,δ for δ ∈ {0.25, 0.1, 0.05} for TC GP with the
combined objective, or with tKS–PIT (TC GP-thres), or rt (TC GP-occ) to select the parameters, and standard GP, all calibrated below
qn,δ . The twCRPS and rt are evaluated on a test grid (Xi , f (Xi ))mi=1 , with m = 4000, and the Xi s are uniformly distributed. tKS–PIT
is approximated using another test grid (X̃i , f (X̃i ))m
i=1 , with m = 4000, where f (X̃i ) ≤ qn,δ and the X̃i are sampled using a subset
simulation algorithm (Bect et al., 2017).

Method Goldstein Price Rosenbrock6 Hartman6 Dixon Price 4 Ackley4

twCRPS rt tKS–PIT twCRPS rt tKS–PIT twCRPS rt tKS–PIT twCRPS rt tKS–PIT twCRPS rt tKS–PIT

n = 30
GP 5.7e2 0.2 0.99 6.8e2 0.036 0.77 0.0016 0.0043 0.61 16 0.015 0.83 0.012 0.0061 0.72
TC GP δ = 0.05 4.9e2 0.12 0.78 6.2e2 0.027 0.56 0.0015 0.0033 0.33 14 0.011 0.68 0.012 0.0047 0.49
TC GP δ = 0.1 4.8e2 0.12 0.72 6.3e2 0.03 0.62 0.0015 0.004 0.59 15 0.012 0.74 0.012 0.0059 0.66
TC GP δ = 0.25 4.5e2 0.13 0.87 6.2e2 0.031 0.67 0.0016 0.0044 0.83 14 0.012 0.73 0.012 0.0063 0.64
RE GP δ = 0.05 1.9 0.021 0.77 39 0.0039 0.66 0.0042 0.0048 0.89 5.3 0.0047 0.74 0.013 0.0046 0.95
RE GP δ = 0.1 2.9 0.019 0.61 54 0.0035 0.43 0.0038 0.0038 0.84 11 0.0055 0.66 0.013 0.0043 0.89
RE GP δ = 0.25 6.7 0.029 0.69 4.1e2 0.018 0.79 0.0018 0.0032 0.72 25 0.015 0.82 0.012 0.006 0.76
Table 3. Comparison of calibration using the twCRPS, rt , and tKS–PIT below t = mn for TC GP with multiple δ, RE GP, and the standard
GP. The calibration methods use thresholds based on qn,δ , while the metrics are evaluated at mn . The twCRPS and rt are evaluated on
a test grid (Xi , f (Xi ))m i=1 , with m = 4000, and the Xi s are uniformly distributed. tKS–PIT is approximated using another test grid
(X̃i , f (X̃i ))m
i=1 , with m = 4000, where f (X̃i ) ≤ mn and the X̃i are sampled using a subset simulation algorithm (Bect et al., 2017).

25
 Goal-Oriented Calibration of GP for BO

iterations and X ∼ µ (here µ is uniform on X), and the fraction of runs that reach a prescribed target level.
Overall, RE GP and TC GP consistently improve upon the standard GP baseline and typically outperform BCR GP and onGP.
An exception is Ackley (d = 10), where onGP achieves the best performance. On the smoother (nearly convex) objectives,
RE GP and TC GP are the top-performing methods; on the Perm function, TC GP yields only modest gains over the GP
baseline compared to RE GP. On highly non-convex functions, TC GP most often provides the largest improvements.
Figure 10 reports diagnostic experiments on moderate-dimensional test functions (d = 10 to d = 20) over 150 BO iterations.
These results extend the controlled benchmark evidence beyond the main text, but do not constitute a high-dimensional BO
evaluation.

Excursion Prob. Proportion Excursion Prob. Proportion

100
10−1
10−2
10−4
10−4
10−7
10−6
10−10
10−8 10−13
10−10 10−16

10−12 10−19

1 50 100 150 200 1 0.2 0.4 0.6 0.8 1 1 50 100 150 200 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion
GP bcrGP reGP δ = 0.25 tcGP δ = 0.05 GP bcrGP reGP δ = 0.25 tcGP δ = 0.05
onGP onGP

(a) Ackley with d = 4 (b) Ackley with d = 10

Excursion Prob. Proportion Excursion Prob. Proportion
100

10−2
10−3
−5
10
10−6

10−8 10−9

10−11 10−12

10−14 10−15

1 50 100 150 200 1 0.2 0.4 0.6 0.8 1 1 50 100 150 200 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion
GP bcrGP reGP δ = 0.25 tcGP δ = 0.05 GP bcrGP reGP δ = 0.25 tcGP δ = 0.05
onGP onGP

(c) Crossintray (d) Dixon–Price with d = 4

Excursion Prob. Proportion Excursion Prob. Proportion
100 101

10−2 10−4

10−4 10−9

10−6 10−14
10−8
10−19
10−10
10−24
10−12
10−29
10−14

1 50 100 150 200 1 0.2 0.4 0.6 0.8 1 1 50 100 150 200 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion
GP bcrGP reGP δ = 0.25 tcGP δ = 0.05 GP bcrGP reGP δ = 0.25 tcGP δ = 0.05
onGP onGP

(e) Goldstein–Price (f) Hartman with d = 6

Figure 8. Comparison of the BO performance with EI as sampling criteria for GP, RE GP with δ = 0.25, onGP and TC GP with δ = 0.05.
Left: median and 10%/90% quantile of pmn = P(f (X) ≤ mn ), where mn is the best observed value so far. Right: fraction of successful
runs reaching a prescribed target level.

26
 Goal-Oriented Calibration of GP for BO

Excursion Prob. Proportion Excursion Prob. Proportion

100 100

10−4 10−3

10 −8 10−6

10−12 10−9

10−16 10−12

10−20 10−15

1 50 100 150 200 1 0.2 0.4 0.6 0.8 1 1 50 100 150 200 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion
GP bcrGP reGP δ = 0.25 tcGP δ = 0.05 GP bcrGP reGP δ = 0.25 tcGP δ = 0.05
onGP onGP

(a) Michalewicz with d = 4 (b) Rosenbrock with d = 6

Excursion Prob. Proportion Excursion Prob. Proportion
1
10
100
10−3
10−2
10−7
10−4
10−11

10−15 10−6

10−19 10−8

10−23
10−10

1 50 100 150 200 1 0.2 0.4 0.6 0.8 1 1 50 100 150 200 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion
GP bcrGP reGP δ = 0.25 tcGP δ = 0.05 GP bcrGP reGP δ = 0.25 tcGP δ = 0.05
onGP onGP

(c) Shekel with m = 10 (d) Perm with d = 6

Figure 9. Comparison of the BO performance with EI as sampling criteria for GP, RE GP with δ = 0.25, onGP and TC GP with δ = 0.05.
Left: median and 10%/90% quantile of pmn = P(f (X) ≤ mn ), where mn is the best observed value so far. Right: fraction of successful
runs reaching a prescribed target level.

27
 Goal-Oriented Calibration of GP for BO

Excursion Prob. Proportion Excursion Prob. Proportion

10−1 10−2
10−3
10−5
10−5
10−8
10−7
10−11
10−9
10−14
10−11
10−17
10−13

1 50 100 150 1 0.2 0.4 0.6 0.8 1 1 50 100 150 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion
GP bcrGP reGP δ = 0.25 tcGP δ = 0.05 GP bcrGP reGP δ = 0.25 tcGP δ = 0.05
onGP onGP

(a) Dixon–Price with d = 10 (b) Dixon–Price with d = 20

Excursion Prob. Proportion Excursion Prob. Proportion
100

10−2 10−2

−4
10
10−5
10−6
10−8
10−8
−10
10 10−11
10−12
10−14
10−14

1 50 100 150 1 0.2 0.4 0.6 0.8 1 1 50 100 150 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion
GP bcrGP reGP δ = 0.25 tcGP δ = 0.05 GP bcrGP reGP δ = 0.25 tcGP δ = 0.05
onGP onGP

(c) Rosenbrock with d = 10 (d) Rosenbrock with d = 15

Figure 10. Comparison of the BO performance in moderate dimension d = 10 to d = 20 for 150 iterations with EI as sampling criteria
for GP, RE GP with δ = 0.25, onGP and TC GP with δ = 0.05. Left: median and 10%/90% quantile of pmn = P(f (X) ≤ mn ), where
mn is the best observed value so far. Right: fraction of successful runs reaching a prescribed target level.

28
 Goal-Oriented Calibration of GP for BO

F.6.2. R ESULTS FOR UCB

We compare a standard GP model and TC GP using the upper confidence bound (UCB) sampling criterion (Lai & Robbins,
1985; Auer et al., 2002; Srinivas et al., 2010), to assess whether tail calibration also affects a criterion other than EI.
As above, BO performance is summarized by the probability of excursion pmn = P(f (X) ≤ mn ), where mn is the
best observed value after n iterations, and by the fraction of runs reaching a prescribed target level, which provides a
complementary view of method reliability.

UCB with a GP model. Fix ε ∈ (0, 1]. With a Gaussian predictive distribution, the (1 − ε) lower confidence bound used
for minimization is
fn (x) − Φ−1 (1 − ε) σn (x), (93)
where Φ denotes the standard normal CDF.

UCB with TC GP. With TC GP parameters (β, λ), the predictive distribution is generalized normal, and the corresponding
(1 − ε) lower confidence bound becomes

fn (x) − Θ−1
β,0,1 (1 − ε) λ σn (x), (94)

where Θβ,0,1 is the CDF defined in (3).

Results are reported in Figures 11 and 12 with UCB and ε = 0.1. Overall, TC GP outperforms the GP baseline on most test
functions. On Ackley (d = 4) the gains are modest and results are similar for Ackley (d = 10); in these cases UCB with a
standard GP is already strong, making improvements harder to obtain.

29
 Goal-Oriented Calibration of GP for BO

Excursion Prob. Proportion Excursion Prob. Proportion

−1
10 10−2

10−4 10−6

10−7 10−10

10−10 10−14
−13
10 10−18

10−16 10−22

1 50 100 150 200 1 0.2 0.4 0.6 0.8 1 1 50 100 150 200 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion
GP bcrGP reGP δ = 0.25 tcGP δ = 0.05 GP bcrGP reGP δ = 0.25 tcGP δ = 0.05
onGP onGP

(a) Ackley with d = 4 (b) Ackley with d = 10

Excursion Prob. Proportion Excursion Prob. Proportion
100
10−1

10−3 10−5

10−6 10−9

10−9 10−13

10 −12 10−17

10−15 10−21

1 50 100 150 200 1 0.2 0.4 0.6 0.8 1 1 50 100 150 200 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion
GP bcrGP reGP δ = 0.25 tcGP δ = 0.05 GP bcrGP reGP δ = 0.25 tcGP δ = 0.05
onGP onGP

(c) Crossintray (d) Dixon–Price with d = 4

Excursion Prob. Proportion Excursion Prob. Proportion
100
−2
10
10−2
10−7
10−4
10−12
10−6
10−17
10−8
10−22
10−10
10−27
10−12
−14
10−32
10

1 50 100 150 200 1 0.2 0.4 0.6 0.8 1 1 50 100 150 200 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion
GP bcrGP reGP δ = 0.25 tcGP δ = 0.05 GP bcrGP reGP δ = 0.25 tcGP δ = 0.05
onGP onGP

(e) Goldstein–Price (f) Hartman with d = 6

Figure 11. Comparison of the BO performance with UCB with ε = 0.1 as sampling criteria for GP, and TC GP with δ = 0.05. Left:
median and 10%/90% quantile of pmn = P(f (X) ≤ mn ), where mn is the best observed value so far. Right: fraction of successful runs
reaching a prescribed target level.

30
 Goal-Oriented Calibration of GP for BO

Excursion Prob. Proportion Excursion Prob. Proportion

101
−1
10
10−3
10−4
10−7

10−11 10−7

10−15 10−10

10−19 10−13

10−23 10−16

1 50 100 150 200 1 0.2 0.4 0.6 0.8 1 1 50 100 150 200 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion
GP bcrGP reGP δ = 0.25 tcGP δ = 0.05 GP bcrGP reGP δ = 0.25 tcGP δ = 0.05
onGP onGP

(a) Michalewicz with d = 4 (b) Rosenbrock with d = 6

Excursion Prob. Proportion Excursion Prob. Proportion
100 100

10−4
10−2
10−8
10−4
10−12
10−6
10−16

10−20 10−8

10−24
10−10

1 50 100 150 200 1 0.2 0.4 0.6 0.8 1 1 50 100 150 200 1 0.2 0.4 0.6 0.8 1
steps proportion steps proportion
GP bcrGP reGP δ = 0.25 tcGP δ = 0.05 GP bcrGP reGP δ = 0.25 tcGP δ = 0.05
onGP onGP

(c) Shekel with m = 10 (d) Perm with d = 6

Figure 12. Comparison of the BO performance with UCB with ε = 0.1 as sampling criteria for GP, and TC GP with δ = 0.05. Left:
median and 10%/90% quantile of pmn = P(f (X) ≤ mn ), where mn is the best observed value so far. Right: fraction of successful runs
reaching a prescribed target level.

31