Kim, Kim, and Song

Selection of the Most Probable Best

Selection of the Most Probable Best

Taeho Kim Department of Industrial and Manufacturing Engineering, Pennsylvania
State University, PA 16802,
<EMAIL_ADDRESS>Kyoung-Kuk Kim School of Management Engineering, Korea
Advanced Institute of Science and Technology, Seoul, Republic of Korea,
<EMAIL_ADDRESS>Eunhye Song School of Industrial and Systems Engineering,
Georgia Institute of Technology, GA 30332,
<EMAIL_ADDRESS>
We consider an expected-value ranking and selection (R&S) problem where all
$k$ solutions’ simulation outputs depend on a common uncertain input model.
Given that the uncertainty of the input model is captured by a probability
simplex on a finite support, we define the most probable best (MPB) to be the
solution whose probability of being optimal is the largest. To devise an
efficient sampling algorithm to find the MPB, we first derive a lower bound to
the large deviation rate of the probability of falsely selecting the MPB, then
formulate an optimal computing budget allocation (OCBA) problem to find the
optimal static sampling ratios for all solution-input model pairs that
maximize the lower bound. We devise a series of sequential algorithms that
apply interpretable and computationally efficient sampling rules and prove
their sampling ratios achieve the optimality conditions for the OCBA problem
as the simulation budget increases. The algorithms are benchmarked against a
state-of-the-art sequential sampling algorithm designed for contextual R&S
problems and demonstrated to have superior empirical performances at finding
the MPB.

Ranking and selection; Input uncertainty; Large deviation theory; Optimal
computing budget allocation; Sequential sampling algorithm

##  1 Introduction

When randomness in a simulation model is driven by input models estimated from
finite data, the simulation output is subject to uncertainty caused by
estimation errors in the input models. This additional uncertainty,
distinguished from inherent stochastic simulation errors, is often referred to
as _input uncertainty_. If the simulation model is applied to find an optimal
design or policy for the target real-world system, then input uncertainty must
be accounted for to make a correct statistical inference on the performance of
a selected solution. In this work, we consider a ranking and selection (R&S)
problem under input uncertainty focusing on the case when all $k$ solutions in
comparison share a common input model estimated from data. Because the
simulation output mean is a functional of the input model, the configuration
of the R&S problem (and its optimum, correspondingly) depends on the input
model, which itself is uncertain. We assume that the uncertainty about the
input model is captured by a probability simplex defined on a finite number
of, say $B$, candidate input models estimated from data. Given the probability
simplex, we define _the most probable best_ (MPB) as the solution that has the
highest likelihood of being the best.

Finding the MPB leads to solving a nested R&S problem; for each candidate
input model, the inner-level R&S problem compares $k$ solutions to find the
conditional best given the input model. Once all $B$ inner problems are
solved, the probability of each solution being optimal given the probability
simplex on the input model can be computed. We refer to this probability as
the _preference probability_ of each solution. The outer-level R&S problem
selects the solution with the largest preference probability as the MPB.

The concept of the MPB has been applied in the medical decision-making
literature (Fenwick et al. 2001, Barton et al. 2008), where the goal is to
assess the probability each treatment is the most cost-effective given a
patient’s willingness-to-pay and uncertainties about a disease. These problems
can be considered as an R&S problem under input model uncertainty.

If a Bayesian approach is taken so that the probability simplex on the input
model is indeed a posterior distribution given data, then the MPB can be
considered the optimal Bayes or maximum a posteriori decision rule (Berger
1985, Robert 2007). This decision rule returns a solution that maximizes the
expected utility under the posterior. In particular, our utility is an
indicator function, which returns $1$, if each solution is the best decision
given a realized input model and $0$, otherwise. This decision rule has been
applied to problems in finance and econometrics (Harvey et al. 2010, Hirano
2010), marketing (Rossi et al. 2012), and classification in machine learning
(Bishop 2006).

Several frameworks have been proposed to account for input uncertainty in R&S
and optimization via simulation in general. We categorize them into three
groups according to the treatment of input uncertainty in their problem
formulation. The first approach is to apply a risk measure with respect to
input model estimation error to the performance measure of each solution, and
then to compare the risk measure values; Corlu and Biller (2015), Wu and Zhou
(2019), Pearce and Branke (2017), and Ungredda et al. (2022) use a risk-
neutral measure (mean) and Xie and Zhou (2015), Wu et al. (2018) and Zhu et
al. (2020) explore (conditional) value at risk. The second category takes the
distributionally robust optimization approach; it first assumes an ambiguity
set on the input model, finds the worst-case input model for each solution
within the ambiguity set, and selects the solution with the best worst-case
performance. Gao et al. (2017) and Fan et al. (2020) fall under this category.
The last focuses on providing a probability guarantee that the best solution
chosen under the best estimate of the real-world input distribution is in fact
optimal; Corlu and Biller (2013), Song et al. (2015) and Song and Nelson
(2019) take this view, however, point out that when estimation error of the
input model is large, the probability guarantee may be low. The MPB
formulation does not fit in any of the three categories above, which warrants
its exploration.

The MPB formulation also lets us assess the model risk in the R&S problem
caused by input uncertainty. If the MPB’s preference probability clearly
dominates all other solutions’, then the risk due to input uncertainty is low.
On the other hand, if there are several solutions whose preference
probabilities are close to the MPB’s, then collecting additional data to
reduce the uncertainty about the input model may reduce the risk of making a
suboptimal decision.

Our objective is to devise a sequential R&S procedure to find the MPB
efficiently. To estimate the mean performance of each solution-input model
pair at any point in the algorithm, we adopt the Bayesian machinery as in
Ryzhov (2016); we impose the normal-normal conjugate models and independent
priors for all pairs. The probability of false selection (PFS) is defined as
the probability that the sample MPB determined from the posterior means of the
Bayesian model is not equal to the true MPB. We take the fixed-budget approach
and formulate an optimal computing budget allocation (OCBA) problem that aims
to minimize the PFS given the budget. Solving the OCBA problem exactly is a
challenge even for a classical R&S problem (Chen et al. 2000). Instead, Glynn
and Juneja (2004) derive the large-deviation rate (LDR) of the PFS and find
the optimality conditions for the static sampling ratios that maximize the LDR
for the classical R&S problem. We take the same approach to derive the
optimality conditions for the asymptotically optimal static sampling ratios
for our nested R&S problem. These optimality conditions depend on the unknown
means of simulation outputs. To reconcile, we adopt a plug-in approach that
replaces the unknown means with the posterior means of the Bayesian model and
utilize the sample-version of the optimality conditions to propose several
sequential sampling algorithms.

Such a plug-in approach to design a sequential algorithm from asymptotic
optimality conditions has been widely adopted in the R&S literature (Fu et al.
2007, He et al. 2007, Chen et al. 2008). However, the static allocation is
only asymptotically optimal and the plug-in policy’s performance may be far
from optimal (Peng et al. 2016, 2018). Indeed, for our problem, a plug-in
algorithm (Algorithm 1) is not guaranteed to find the MPB almost surely
(a.s.). Instead, the version (Algorithm 2) that combines the plug-in estimates
and posterior sampling is shown to satisfy the optimality conditions for the
static allocation ratios.

Alternatively, a dynamic program (DP) can be formulated to maximize the
expected reward (e.g., $0$-$1$ reward or negative regret) after a finite
simulation budget is spent. However, solving the DP exactly is computationally
challenging even for a moderate problem size. Instead, myopic or one-step
optimal policies have been proposed and studied (see Chick and Inoue (2001),
Frazier et al. (2008), Ryzhov (2016), Peng and Fu (2017) for instance).

Notably, there are some recent work that bridge the gap between the two
abovementioned approaches by proving that a sequential sampling algorithm
based on a myopic policy indeed achieves the asymptotically optimal static
allocations as the simulation budget grows. Chen and Ryzhov (2019b) propose a
sequential sampling algorithm based on the complete expected improvement
(CEI), which is originally created as a myopic sampling policy (Salemi et al.
2019), and prove that its sampling ratios for the solutions asymptotically
satisfy the optimality conditions in Glynn and Juneja (2004). Avci et al.
(2021) derive the same result for the policy that applies the gradient of the
CEI with respect to the sample size as a sampling rule.

While these results are desirable, in this paper, we focus on studying the
optimal static allocation ratios and creating sequential policies from them as
this is the first paper that formalizes selection of the MPB as a nested R&S
problem. We believe that these results will provide a foundation to create and
analyze a DP-based sampling rule for the MPB formulation in our future study.

If one views a candidate input model as a context, our work can be connected
to the contextual R&S studied by Gao et al. (2019), Li et al. (2020), and Shen
et al. (2021). The contextual R&S solves an R&S problem conditional on a
covariate vector. However, at the outer level, it aims to maximize the mean
(or worst-case) probability of correct selection over all possible covariates.
Therefore, there is no outer-level R&S problem for the contextual R&S as in
the MPB formulation. Such a difference between two formulations manifests in
the optimality conditions for the respective OCBA problems as we discuss in
Section 4.

Because of the unique nested structure, deriving the LDR of the PFS for the
MPB formulation has a technical challenge that is not present in classical or
contextual R&S. The outer-level R&S problem’s correct selection depends on how
well the inner-level problems are solved under each candidate input model.
However, the importance of each inner-level problem for achieving the outer-
level correct selection differs; even if false selection occurs at some inner-
level problems, the MPB may still be found correctly. Hence, the MPB
formulation requires quantifying the importance of each inner-level problem
for outer-level correct selection; there is no straightforward extension of an
existing R&S algorithm that provides this feature. We introduce the _balance
weight_ to measure how important it is to sample each solution-input model
pair in achieving the outer-level correct selection. The balance weights
appear in the opitmality conditions for our OCBA problem.

As a byproduct, we also prove asymptotic convergence of the sequential
sampling algorithm for the contextual R&S problem stated in Gao et al. (2019)
without a full proof; only a sketch of proof is provided. This result is an
intermediate step to show our Theorem 5.3 and we could not find a rigorous
proof in the literature at the time of submission of this paper.

A preliminary version of this work published in a conference proceeding
discusses the optimality conditions for the static sampling allocations
without a proof and the sequential sampling algorithms are analyzed only
empirically (Kim et al. 2021). In this work, we provide a complete set of
theoretical analyses for both static and sequential sampling algorithms we
propose.

The remainder of the paper is organized as follows. In Section 2, we formally
define our problem of interest. Section 3 provides background on the Bayesian
learning model. Section 4 presents the OCBA formulation for the MPB problem
and derive the optimality conditions for the static optimal sampling ratios
for all solution-input model pairs. Based on this result, sequential sampling
algorithms are introduced and their convergence properties are discussed in
Section 5. All algorithms are benchmarked against a state-of-the-art
contextual R&S algorithm in empirical studies in Section 6. Conclusions are
given in Section 7. Due to the page limit, we omit the proofs of all
theoretical results from the paper and provide them in the Electronic
Companion of this paper.

##  2 Problem Definition

Suppose the input model of the simulator is parameterized by $\theta$ whose
value is uncertain. We assume that the distribution of $\theta$ can be modeled
by a probability simplex $\pi$ defined on a finite support. Note that this
framework can also accommodate a nonparametric input model with a finite
support. Let $y_{i}(\theta):=\mathrm{E}\left[Y_{i}(\theta)|\theta\right]$ be
the conditional mean of the $i$th solution’s simulation output given $\theta$.
The most probable best (MPB) is defined as

$i^{*}:=\argmax\nolimits_{1\leq i\leq
k}\mathrm{P}_{\pi}\left\\{y_{i}(\theta)\leq y_{j}(\theta),\forall j\neq
i\right\\}.$ (1)

We refer to the probability in (1) as the _preference probability_ of the
$i$th solution; it measures the probability that Solution $i$ performs no
worse than any other solutions under $\pi$. Therefore, the MPB is the solution
with the largest preference probability under $\pi$.

Although the MPB can be defined for general distribution $\pi$, in this paper,
we focus on finding the MPB when $\theta$ has finite support
$\Theta=\\{\theta_{1},\theta_{2},\ldots,\theta_{B}\\}$. Thus, $\pi$ can be
represented as probability simplex
$\pi=\\{p_{1},p_{2},\ldots,p_{B}|\sum_{b}p_{b}=1\\}$ on $\Theta$. Under this
assumption, (1) can be rewritten as

$i^{*}=\argmax\nolimits_{1\leq i\leq
k}\sum\nolimits_{b=1}^{B}p_{b}\mathbf{1}\left\\{y_{i}(\theta_{b})=y_{i^{b}}(\theta_{b})\right\\},$
(2)

where $i^{b}:=\argmin_{i}y_{i}(\theta_{b})$ is the conditional optimum at
$\theta_{b}$. Finding $i^{*}$ can be formulated as a nested R&S problem;
selecting $i^{b}$ given $\theta_{b}$ corresponds to an inner-level R&S
problem. The outer-level R&S problem then compares all $k$ solutions’
preference probabilities to find $i^{*}$.

To simplify the analysis, we assume (i) $i^{b}$ is unique for each
$\theta_{b}$; and (ii) $i^{*}$ is unique. With these assumptions, (2) further
simplifies to $i^{*}=\argmax_{1\leq i\leq
k}\sum_{b=1}^{B}p_{b}\mathbf{1}\left\\{i=i^{b}\right\\}$. Additionally, we
define the favorable set of $i$ to be $\Theta_{i}:=\\{\theta_{b}:i=i^{b},1\leq
b\leq B\\}$. In words, $\Theta_{i}$ is the set of $\theta$ at which $i$
outperforms other solutions. The preference probability of each $i$ can be
written as $\mathrm{P}_{\pi}(\Theta_{i})$.

Our objective is to create a simulation optimization algorithm to find the MPB
in (2) efficiently given simulation budget $n$. Suppose an algorithm returns
$i^{*}_{n}$ as an estimate for $i^{*}$ after spending $n$ simulation
replications. In case two or more solutions tie in the estimated preference
probabilities, $i^{*}_{n}$ may be a set of solutions. One measure of
efficiency of the algorithm is the _probability of correct selection_ (PCS),
where the correct selection event is defined as
CS$:=\\{i_{n}^{*}=\\{i^{*}\\}\\}$. Given $n$, an efficient algorithm would
maximize the PCS. Equivalently, one can define the _false selection_ event as
FS$:=\\{i_{n}^{*}\neq\\{i^{*}\\}\\}$ and attempt to minimize the _probability
of false selection_ (PFS), which is the approach we take in this paper. Note
that any event of a tie is included in FS as $i^{*}$ is assumed to be unique.
Furthermore, we define the conditional FS event at $\theta_{b}$ as
FS${}_{b}:=\\{i^{b}_{n}\neq\\{i^{b}\\}\\}$, where $i^{b}_{n}$ is the estimated
$i^{b}$ after spending $n$ replications. If at each $\theta_{b}$, conditional
correct selection occurs, then we have CS on the outer-level problem. However,
even if FSb occurs at some $\theta_{b}$, we may still achieve CS at the outer-
level as long as the estimated preference probability of $i^{*}$ is the
largest.

###  2.1 Relations to other formulations

Daskin et al. (1997) consider a value-at-risk (VaR) regret minimization
problem. Under our notation, this approach can be formulated as

$\min\nolimits_{1\leq i\leq
k}\text{VaR}_{\pi}^{1-\alpha}\\{y_{i}(\theta)-\min\nolimits_{1\leq j\leq
k}y_{j}(\theta)\\},$ (3)

where $\text{VaR}_{\pi}^{1-\alpha}$ is the $1-\alpha$ level VaR with respect
to $\pi$. Observe that for some $i$, if $\mathrm{P}_{\pi}(\Theta_{i})\geq
1-\alpha$, then $\text{VaR}_{\pi}^{1-\alpha}\\{y_{i}(\theta)-\min_{1\leq j\leq
k}y_{j}(\theta)\\}=0$. Hence, the MPB belongs to the solution set of (3) when
its preference probability is no less than the risk level, $1-\alpha$.

The contextual R&S regards $\theta$ as a context that parameterizes each R&S
problem and aims to maximize the average or worst-case PCS. Under our
notation, the two measures can be written as
$\mathrm{E}_{\pi}\left[\mathrm{P}(i^{b}_{n}=i^{b})\right]=\sum_{b=1}^{B}p_{b}\mathrm{P}(i^{b}_{n}=i^{b})$
and $\min_{1\leq b\leq B}\mathrm{P}(i^{b}_{n}=i^{b})$, respectively, which are
shown to be asymptotically equivalent regardless of $\\{p_{b}\\}$, if
$p_{b}>0$ for all $1\leq b\leq B$ (Gao et al. 2019). Contextual R&S puts
emphasis on solving all $\theta$-problems equitably, which clearly differs
from solving (2) to find $i^{*}$. Nevertheless, one can still apply a
contextual R&S algorithm to solve (2) hoping that if each $\theta_{b}$ problem
is solved correctly, then $i^{*}$ is found correctly. However, this would not
be as efficient as algorithms we tailor to solve (2) specifically. See Section
6 for empirical comparisons.

Lastly, the MPB formulation is closely connected to the multinomial R&S
(Bechhofer et al. 1959, Chen 1992, Tollefson et al. 2014). Their objective is
to select a solution among $k$ with the largest probability such that its
simulation output is the largest (or smallest) among all $k$ simulation
outputs. Our problem has a nested structure that the classical multinomial R&S
does not; we solve the mean-value inner R&S problems, then solve the
multinoimal R&S at the outer level.

##  3 Bayesian Learning Model

To facilitate sequential learning of $y_{i}(\theta)$, we introduce stochastic
process $\eta_{i}(\theta)$ to model uncertainty about $y_{i}(\theta)$. Instead
of taking a purely Bayesian stance, we view $y_{i}(\theta)$ to be fixed and
adopt $\eta_{i}(\theta)$ as a model to simply assist the learning process of
$y_{i}(\theta)$; a similar view is taken in Ryzhov (2016). Adopting the
normal-normal conjugate model for $\eta_{i}(\theta)$ and $Y_{i}(\theta)$, we
assume for each $(i,\theta)$

$\displaystyle Y_{i}(\theta)\sim
N(\eta_{i}(\theta),\lambda^{2}_{i}(\theta)),\hskip 20.0pt\eta_{i}(\theta)\sim
N({\mu}_{i,0}(\theta),\sigma^{2}_{i,0}(\theta)),$ (4)

where $\lambda_{i}(\theta)$ is the simulation error variance of
$Y_{i}(\theta)$ given $\theta$, and $\mu_{i,0}(\theta)$ and
$\sigma_{i,0}^{2}(\theta)$ are the mean and variance of the prior distribution
of $\eta_{i}(\theta)$, respectively. Model (4) assumes that the simulation
error at each $(i,\theta)$ is normally distributed. For a non-normal
simulation error, batching can be adopted (Kim and Nelson 2006). We
empirically demonstrate robustness of batching and the normality assumption
for non-normal simulation outputs in Section 6. Note that
$\lambda_{i}(\theta)$ is assumed to be known throughout the paper until
Section 6.2. For unknown variances, (4) can be replaced with a normal-gamma
model; See 8 for details.

After $n$ simulation replications are made, let $N^{n}_{i}(\theta)\geq 1$
denote the number of replications allocated to $(i,\theta)$. Conditional on
the simulation outputs,
$Y_{i1}(\theta),Y_{i2}(\theta),\ldots,Y_{iN^{n}_{i}(\theta)}(\theta)$, the
posterior mean and variance, $\mu_{i,n}(\theta)$ and
$\sigma_{i,n}^{2}(\theta)$, of $\eta_{i}(\theta)$ are updated as

$\displaystyle\mu_{i,n}(\theta)=\sigma_{i,n}^{2}(\theta)\left(\frac{\mu_{i,0}(\theta)}{\sigma_{i,0}^{2}(\theta)}+\frac{N^{n}_{i}(\theta)}{\lambda_{i}^{2}(\theta)}\frac{1}{N^{n}_{i}(\theta)}\sum\nolimits_{r=1}^{N^{n}_{i}(\theta)}Y_{ir}(\theta)\right),\;\;\sigma_{i,n}^{2}(\theta)=\left(\frac{1}{\sigma_{i,0}^{2}(\theta)}+\frac{N^{n}_{i}(\theta)}{\lambda_{i}^{2}(\theta)}\right)^{-1}.$
(5)

We adopt the noninformative prior, $\sigma_{i,0}^{2}(\theta)=\infty$, for all
$(i,\theta)$. For this choice, (5) simplifies to
$\mu_{i,n}(\theta)=\sum_{r=1}^{N^{n}_{i}(\theta)}Y_{ir}(\theta)/N^{n}_{i}(\theta)$
and $\sigma_{i,n}^{2}(\theta)=\lambda_{i}^{2}(\theta)/N^{n}_{i}(\theta)$
regardless of the choice for $\mu_{i,0}(\theta).$ Notice that
$\mu_{i,n}(\theta)$ and $\sigma_{i,n}^{2}(\theta)$ are same as the sample mean
and its variance in the frequentist’s setting. We adopt
$\mu_{i,n}(\theta_{b})$ as the estimator for $y_{i}(\theta_{b})$ for all
$1\leq i\leq k,1\leq b\leq B$.

##  4 Asymptotic Analysis of Sampling Allocation

In this section, we study the asymptotically optimal sampling ratios for all
$kB$ solution-parameter pairs that maximizes the convergence rate of PFS for
Problem (2) as $n$ increases.

From the Bayesian learning model in Section 3, $i^{b}$ for each $\theta_{b}$
and $i^{*}$ can be estimated by $i^{b}_{n}=\argmin\nolimits_{1\leq i\leq
k}\mu_{i,n}(\theta_{b})$ and $i^{*}_{n}=\argmax\nolimits_{1\leq i\leq
k}\sum\nolimits_{b=1}^{B}p_{b}\mathbf{1}\left\\{i=i^{b}_{n}\right\\},$
respectively. From the normality assumption on the simulation output,
$i^{b}_{n}$ is unique with probability $1$; in fact, even if the simulation
output is non-normal, as long as it is continuously distributed, no tie occurs
at each inner-level problem. On the other hand, there may be ties at the
outer-level problem for finite $n$ in particular when $p_{b}=1/B$ for all
$1\leq b\leq B$.

Characterizing the PFS for finite $n$ is challenging in general for an R&S
problem, but especially difficult for our setting due to the nested structure
of (2); the complex relationship between FS event at the outer-level problem
and FSb events at the inner levels makes the analysis more challenging. The
approach we take is to formulate and solve the OCBA problem when $n\to\infty$
to find the asymptotically optimal static sampling ratios and use them to
guide our sequential sampling algorithms in Section 5.

If simulation budget $n$ is spent judiciously, the FS event becomes a rare
event as $n$ increases. In general, the large-deviation rate (LDR) of sequence
of rare events $\\{E_{n}\\}$ is defined as
$\lim_{n\rightarrow\infty}-\frac{1}{n}\log\mathrm{P}(E_{n})=I,$ provided that
$I>0$ exists. Namely, $I$ represents the exponential decay rate of the
probability of the sequence of events, $\left\\{E_{n}\right\\}$. When such a
limit does not exist, one can consider
$\lim\inf_{n\rightarrow\infty}-\frac{1}{n}\log\mathrm{P}(E_{n})$, instead.

We first focus on characterizing the static allocation rule that maximizes the
LDR. Let $\bm{\alpha}$ be the vector of $\\{\alpha_{i}(\theta_{b})\\}$, where
each $\alpha_{i}(\theta_{b})$ represents the asymptotic sampling ratio of $n$
allocated to $(i,\theta_{b})$. We assume that
$N_{i}^{n}(\theta_{b})\rightarrow\infty$ even if $\alpha_{i}(\theta_{b})=0$;
i.e., $N_{i}^{n}(\theta_{b})$ diverges sublinearly in $n$.

Let
$\mathrm{LDR}_{\text{FS}}:=\liminf_{n\rightarrow\infty}-\frac{1}{n}\log\mathrm{P}\left(i_{n}^{*}\neq\\{i^{*}\\}\right)$.
Then, the OCBA problem of our interest is

$\displaystyle\max\nolimits_{\bm{\alpha}}\mathrm{LDR}_{\text{FS}}\;\;\;\;\textrm{subject
to}\;\;\sum\nolimits_{i=1}^{k}\sum\nolimits_{b=1}^{B}\alpha_{i}(\theta_{b})=1,\alpha_{i}(\theta_{b})\geq
0,1\leq i\leq k,1\leq b\leq B.$ (6)

Problem (6) finds an optimal allocation $\bm{\alpha}$ that maximizes the limit
infimum on the exponential decay rate of the PFS in the limit. Clearly,
$\mathrm{LDR}_{\text{FS}}$ depends on $\bm{\alpha}$, but finding the explicit
dependence is challenging. In Section 4.1, we analyze a lower bound of
$\mathrm{LDR}_{\text{FS}}$ as a function of $\bm{\alpha}.$ In Section 4.2, we
present a relaxation of Problem (6) by replacing its objective with
(maximization of) the lower bound on $\mathrm{LDR}_{\text{FS}}$ and derive its
optimality conditions.

###  4.1 Optimal Computing Budget Allocation Formulation

The FS occurs when a subset of $\\{$FS${}_{b}\\}$ occurs. As pointed out in
Section 2, not all of {FSb} are equally critical for CS. Hence, the key to
characterizing $\text{LDR}_{\text{FS}}$ is to understand which subset of {FSb}
leads to FS at the outer-level problem.

Let $d_{j}$ be the difference in preference probabilities between $i^{*}$ and
$j$:

$\displaystyle d_{j}$
$\displaystyle:=\sum\nolimits_{b=1}^{B}p_{b}\mathbf{1}\left\\{i^{*}=i^{b}\right\\}-\sum\nolimits_{b=1}^{B}p_{b}\mathbf{1}\left\\{j=i^{b}\right\\}.$

Namely, $d_{j}$ measures how much inferior $j$ is to $i^{*}$ at the outer-
level problem; the smaller $d_{j}$ is the more competitive $j$ is against
$i^{*}$. Since $i^{*}$ is unique by assumption, $d_{j}>0$ for all $j\neq
i^{*}$ and $d_{i^{*}}=0$. When some FSb leads to overestimation of the
preference probability of Solution $j$ whose $d_{j}$ is small, the chance of
FS in the outer level increases. We formalize this idea in the following.

For each inner problem given $\theta_{b}$, there are $k$ possible solutions
that may be returned as $i^{b}_{n}$, whereas only one of them is correct. To
consider all inner problems jointly, we define

$\mathcal{M}:=\left\\{\BFM|\BFM\in\\{0,1\\}^{k\times
B},\BFM^{\top}\BFone_{k}=\BFone_{B}\right\\},$ (7)

where $\BFone_{k}$ denotes the $k$-dimensional column vector of ones. For each
$k\times B$ matrix $\BFM\in\mathcal{M}$, let $m_{i,b}$ denote its $(i,b)$th
element. The equality constraint in (7) implies that only one $m_{i,b}$ in
each column equals $1$. If we regard $m_{i,b}$ to be an indicator variable to
specify whether $i$ is the conditional optimum at $\theta_{b}$, then
$\mathcal{M}$ is the collection of all possible conditional optima mappings
for a generic $k\times B$ problem while only one element of $\mathcal{M}$
corresponds to the true configuration of Problem (2). Let $\BFM^{*}$ denote
the true mapping. Then, $\mathcal{M}\backslash\left\\{\BFM^{*}\right\\}$
enumerates all possible ways one or more FSbs jointly occur. The question is,
which of these mappings lead to FS in the outer level? To answer this, we
define $d_{j}(\BFM)$ and $I(\BFM)$ for each $\BFM\in\mathcal{M}$ as

$d_{j}(\BFM):=\sum\nolimits_{b=1}^{B}p_{b}m_{i^{*},b}-\sum\nolimits_{b=1}^{B}p_{b}m_{j,b}\;\;\;\mbox{and}\;\;\;I(\BFM):=\left\\{(i,\theta_{b})|m_{i,b}=1,i\neq
i^{b}\right\\}.$ (8)

In words, $d_{j}(\BFM)$ measures the difference in preference probabilities of
$i^{*}$ and $j$ under $\BFM$, and $I(\BFM)$ is the set of misspceified
conditional optimum-parameter pairs under $\BFM$. By definition,
$d_{j}(\BFM^{*})=d_{j}$ for all $1\leq j\leq k,$ and $I(\BFM^{*})=\emptyset.$
After $n$ simulation replications, consider realized conditional optima
$\left\\{i_{n}^{b}\right\\}_{b=1}^{B}$. Define $k\times B$ matrix $\BFM_{n}$
whose $(i,b)$th element is $m_{i,b}^{n}=\mathbf{1}\\{i=i_{n}^{b}\\}.$ Then,
there exists a unique element of $\mathcal{M}$ that coincides with $\BFM_{n}$.
Recall that we define FS as the event such that there exists $j\neq i^{*}$
that ties with or dominates $i^{*}$ at the outer-level problem. This implies
that $\BFM_{n}$ leads to FS, if and only if, $d_{j}(\BFM_{n})\leq 0$ for some
$j\neq i^{*}$. To this end, we define

$\mathcal{A}_{j}:=\left\\{\BFM\in\mathcal{M}|d_{j}(\BFM)\leq 0\right\\}$ (9)

for each $j\neq i^{*}.$ In words, $\mathcal{A}_{j}$ is a collection of
mappings that make $j$ as good as or better than $i^{*}$ in the outer-level
problem. Thus, if $\BFM_{n}\in\mathcal{A}_{j}$ for some $j\neq i^{*}$, then FS
occurs. Equivalently, FS$=\bigcup_{j\neq
i^{*}}\left\\{\BFM_{n}\in\mathcal{A}_{j}\right\\}$. Observe that $\max_{j\neq
i^{*}}\mathrm{P}(\BFM_{n}\in\mathcal{A}_{j})\leq$ PFS $\leq\sum_{j\neq
i^{*}}\mathrm{P}(\BFM_{n}\in\mathcal{A}_{j})$, which implies

$\displaystyle\liminf\nolimits_{n\rightarrow\infty}-\tfrac{1}{n}\log\text{PFS}=\min\nolimits_{j\neq
i^{*}}\liminf\nolimits_{n\rightarrow\infty}-\tfrac{1}{n}\log\mathrm{P}(\BFM_{n}\in\mathcal{A}_{j}).$
(10)

Recall that the left-hand side of (10) is the definition of
$\mathrm{LDR}_{\mathrm{FS}}$. Thus, (10) stipulates that characterizing
$\mathrm{LDR}_{\mathrm{FS}}$ boils down to finding the LDR of
$\mathrm{P}\left(\BFM_{n}\in\mathcal{A}_{j}\right)$ for all $j\neq i^{*}$.

Before introducing our main result, we first present an auxiliary result in
Proposition 4.1 that shows the LDR of falsely selecting $i\neq i^{b}$ as
$i_{n}^{b}$ at $\theta_{b}$. Proposition 4.1 and subsequent analyses require
the following rate function to be defined for each $i\neq i^{b}$:

$G_{i}(\theta_{b}):=\frac{(y_{i}(\theta_{b})-y_{i^{b}}(\theta_{b}))^{2}}{2\left(\lambda_{i}^{2}(\theta_{b})/\alpha_{i}(\theta_{b})+\lambda_{i^{b}}^{2}(\theta)/\alpha_{i^{b}}(\theta_{b})\right)}.$

We let $G_{i}(\theta_{b})=0$, if either $\alpha_{i}(\theta_{b})$ or
$\alpha_{i^{b}}(\theta_{b})$ is $0$. Note that $G_{i}(\theta_{b})$ is the LDR
of the probability of making an incorrect pairwise comparison between $i$ and
$i^{b}$ under the normality assumption on $Y_{i}(\theta_{b})$ for the inner-
level R&S problem at $\theta_{b}$ (Glynn and Juneja 2004).

###### Proposition 4.1

Given fixed $\bm{\alpha}$,
$\liminf\nolimits_{n\rightarrow\infty}-\frac{1}{n}\log$ for each $i\neq
i^{b}$.

Proposition 4.1 is a step towards finding a lower bound for (10). Roughly
speaking, the LDR of is determined by $\BFM$ in $\mathcal{A}_{j}$ such that
decays at the slowest rate. Theorem 4.2 below formalizes this idea to obtain a
lower bound to $\mathrm{LDR}_{\text{FS}}$ for any fixed $\bm{\alpha}$.

###### Theorem 4.2

Given fixed $\bm{\alpha}$, $\mathrm{LDR}_{\textup{FS}}\geq\min_{j\neq
i^{*}}\mathrm{LDR}_{j,i^{*}}$, where

$\mathrm{LDR}_{j,i^{*}}:=\min\nolimits_{\BFM\in\mathcal{A}_{j}}\sum\nolimits_{(i,\theta_{b})\in
I(\BFM)}G_{i}(\theta_{b}).$ (11)

From Theorem 4.2, the following approximate OCBA framework can be formulated
to find the MPB replacing the objective function in (6) with the lower bound
in Theorem 4.2:

$\displaystyle\mathrm{LDR}^{*}:=\max_{\bm{\alpha}}\min_{j\neq
i^{*}}\mathrm{LDR}_{j,i^{*}}\;\;\textrm{subject
to}\;\;\sum\nolimits_{i=1}^{k}\sum\nolimits_{b=1}^{B}\alpha_{i}(\theta_{b})=1,\alpha_{i}(\theta_{b})\geq
0,1\leq i\leq k,1\leq b\leq B.$ (12)

###  4.2 Relaxed OCBA and its Optimality Conditions

The expression for $\mathrm{LDR}_{j,i^{*}}$ in (11) involves $\mathcal{A}_{j}$
defined in (9). Unfortunately, computational complexity of characterizing
$\mathcal{A}_{j}$ increases combinatorially in $k$ and $B$. To avoid this, in
this section, we relax (12) to obtain a concave program in $\bm{\alpha}$ whose
optimality conditions can be easily derived by the Karush-Kuhn-Tucker (KKT)
conditions.

For ease of exposition, we first define two solution-parameter sets:
$\Xi^{\text{adv}}:=\left\\{(i^{*},\theta_{b}):\theta_{b}\in\Theta^{c}_{i^{*}}\right\\}$
includes all parameters in the adversarial set of $i^{*}$ paired with $i^{*}$;
and $\Xi:=\left\\{(i,\theta_{b})|i\neq i^{*},i\neq i^{b}\right\\}$. Therefore,
$\Xi\cup\Xi^{\text{adv}}$ includes all possible parameter-suboptimal pairs.

Observe that (11) has an embedded minimization problem. In words,
$\mathrm{LDR}_{j,i^{*}}$ is determined by $\BFM\in\mathcal{A}_{j}$ whose sum
of rates $G_{i}(\theta_{b})$ of the misspecified conditional optimum pairs
$(i,\theta_{b})$ in $I(\BFM)$ is the smallest. Conveniently, (11) can be
formulated as a (minimization version of a) knapsack problem (Martello and
Toth 1990) with an additional constraint as follows:

$\min\;\BFg^{\top}\BFm\hskip 15.0pt\text{ subject to
}\;\BFv_{j}^{\top}\BFm\geq d_{j},\hskip 10.0pt\BFM\in\mathcal{M},$ (13)

where $\BFm$ is the vectorized version of $\BFM$ by taking its elements
corresponding to $\Xi\cup\Xi^{\text{adv}}$ only; $\mathbf{g}$ is the vector of
$G_{i}(\theta_{b})$ for all $(i,\theta_{b})\in\Xi\cup\Xi^{\text{adv}}$;
$\mathbf{v}_{j}$ is the vector of
$v_{j}[(i,\theta_{b})]^{+}=\max\\{v_{j}[(i,\theta_{b})],0\\}$, and

$v_{j}[(i,\theta_{b})]=p_{b}\left(\mathbf{1}\left\\{j=i\right\\}-\mathbf{1}\left\\{j=i^{b}\right\\}-\mathbf{1}\left\\{i^{*}=i\right\\}+\mathbf{1}\left\\{i^{*}=i^{b}\right\\}\right)$

for all $(i,\theta_{b})\in\Xi\cup\Xi^{\text{adv}}$. Elements of both
$\mathbf{g}$ and $\mathbf{v}$ are sorted in the same order as they appear in
$\mathbf{m}$. Recall from (7) that $\BFM\in\mathcal{M}$ implies that the
elements of $\BFM$ are binary and each column sum is one; other than the
latter condition, (13) has the exact structure of a knapsack problem. While
detailed derivation of (13) can be found in Appendix 9, we provide an
intuition here. The objective function is obtained from that (11) minimizes
the sum of $G_{i}(\theta_{b})$’s. The inequality constraint of (13) is
equivalent to requiring $\BFM\in\mathcal{A}_{j}$. From (9),
$\BFM\in\mathcal{A}_{j}$, if and only if, $d_{j}(\BFM)\leq 0$ and
$v_{j}[(i,\theta_{b})]$ measures how much misspecifying $i$ to be $i^{b}$ at
$\theta_{b}$ affects $d_{j}-d_{j}(\BFM)$. So, the knapsack analogy is: (13)
finds the set of misspecified $(i,\theta_{b})$ pairs with the smallest sum of
LDRs that is just enough to cause $j$ to (incorrectly) outperform or tie with
$i^{*}$ in the outer-level problem.

There is no closed-form expression for the optimal solution of the knapsack
problem in general. Since our goal is to obtain easy-to-compute optimality
conditions for $\bm{\alpha}$ to guide a sequential sampling algorithm, we drop
the integrality constraints in (13) and compute the optimal objective function
value of the relaxed problem analytically. This provides an analytical
expression for the lower bound for $\mathrm{LDR}_{j,i^{*}}$ for each $j\neq
i^{*}$. Theorem 4.3 provides the analytical expressions.

###### Theorem 4.3

Let $W_{i}(\theta_{b})=\infty$ for all $(i,\theta_{b})\in\Xi^{c}$. For each
$(i,\theta_{b})\in\Xi,$ let

$W_{i}(\theta_{b})=\begin{cases}\max\left\\{\min\left(\min_{j\neq
i^{*}}d_{j},\frac{d_{i}}{2}\right)/p_{b},1\right\\},&\mbox{if
}\theta_{b}\in\Theta_{i^{*}}\\\ \max\left\\{d_{i}/p_{b},1\right\\},&\mbox{if
}\theta_{b}\notin\Theta_{i^{*}}\end{cases}.$ (14)

Given fixed $\bm{\alpha},$ $\min_{j\neq
i^{*}}\mathrm{LDR}_{j,i^{*}}\geq\underline{\mathrm{LDR}}:=\min_{(i,\theta_{b})\in\Xi}W_{i}(\theta_{b})G_{i}(\theta_{b})$.

We refer to $W_{i}(\theta_{b})$ as the _balance weight_ at each
$(i,\theta_{b})$. Intuitively, $W_{i}(\theta_{b})$ reflects the _importance of
misspecifying $(i,\theta_{b}$) to the FS event_. If $G_{i}(\theta_{b})$ is the
same for all $(i,\theta_{b})\in\Xi$, then the smallest $W_{i}(\theta_{b})$
determines $\underline{\mathrm{LDR}}$. Given $\theta_{b}$, the smaller $d_{i}$
is, the smaller $W_{i}(\theta_{b})$ is; this makes sense as smaller $d_{i}$
implies that Solution $i$ is more competitive against $i^{*}$ and therefore,
misspecifying $i$ as $i^{b}$ is more critical to FS. Given $i$,
$W_{i}(\theta_{b})$ is smaller when $\theta_{b}\in\Theta_{i^{*}}$; again, this
is intuitive as we would like to make fewer mistakes at specifying $i^{*}$ to
be the optima at the parameters in $\Theta_{i^{*}}$ to avoid underestimating
its preference probability.

Using $\underline{\mathrm{LDR}}$ defined in Theorem 4.3, an approximate
version of (12) is formulated as

$\displaystyle\underline{\mathrm{LDR}}^{*}:=\max_{\bm{\alpha}}\underline{\mathrm{LDR}}\;\;\;\;\textrm{subject
to}\;\;\sum\nolimits_{i=1}^{k}\sum\nolimits_{b=1}^{B}\alpha_{i}(\theta_{b})=1,\alpha_{i}(\theta_{b})\geq
0,1\leq i\leq k,1\leq b\leq B.$ (15)

Theorem 4.3 implies $\underline{\mathrm{LDR}}^{*}\leq\mathrm{LDR}^{*}$. Unlike
(12), the optimality conditions for (15) can be derived without fully
characterizing $\left\\{\mathcal{A}_{j}\right\\}_{j\neq i^{*}}$. Since
$\underline{\mathrm{LDR}}$ is a concave function of $\bm{\alpha}$ and the
constraints of (15) are linear, the KKT conditions provide its optimality
conditions as stated in Theorem 4.4.

######  Theorem 4.4 (Optimality conditions for (15))

Any ${\bm{\alpha}}$ is optimal for (15), if and only if, ${\bm{\alpha}}$
satisfies the following system of equations:

  * •

(Global balance condition) For all $1\leq b\leq B$,
$\;{{\alpha}^{2}_{i^{b}}(\theta_{b})}/{\lambda^{2}_{i^{b}}(\theta_{b})}=\sum_{i\neq
i^{b},i_{n}^{*}}{{\alpha}^{2}_{i}(\theta_{b})}/{\lambda^{2}_{i}(\theta_{b})}.$

  * •

(Pairwise balance condition) For all
$(i,\theta_{b}),(j,\theta_{b^{\prime}})\in\Xi$, $\\{W_{i}(\theta_{b})\\}$
given in (14) satisfies

$W_{i}(\theta_{b})G_{i}(\theta_{b})=W_{j}(\theta_{b^{\prime}})G_{j}(\theta_{{b^{\prime}}}).$
(16)

  * •

${\alpha}_{i^{*}}(\theta_{b})=0$ for all $\theta_{b}\in\Theta^{c}_{i^{*}}$.

The last condition states that the asymptotically optimal sampling ratio for
$i^{*}$ at any adversarial $\theta_{b}$ is $0$. This may be surprising as the
optimal sampling ratios for the classical R&S problem are strictly positive
(Glynn and Juneja 2004). This stark difference has implications in our
setting. Suppose $i^{*}$ is correctly identified as the conditional optimum at
all $\theta_{b}\in\Theta_{i^{*}}$. Then, even if $i^{*}$ is incorrectly
identified as the conditional optimum at some
$\theta_{b}\in\Theta^{c}_{i^{*}}$, CS still occurs as long as the best among
$i\neq i^{*}$ is correctly identified not to overestimate the preference
probability of any $i\neq i^{*}$.

The global balance condition ensures that $i^{b}$ at each $\theta_{b}$ is
correctly identified by allocating enough replications to each $i^{b}$
relative to other solutions. One can observe that this condition is identical
to the global balance condition for the classical R&S problem (Glynn and
Juneja 2004) as the optimal ${\alpha}_{i^{*}}(\theta_{b})=0$ for all
$\theta_{b}\in\Theta^{c}_{i^{*}}$ in our setting according to the last
condition in Theorem 4.4.

The pairwise balance condition matches the product of $W_{i}(\theta_{b})$ and
$G_{i}(\theta_{b})$ for all suboptimal solution-parameter pairs. The balance
weight for each $G_{i}(\theta_{b})$ depends on the performance of $i$. Large
$W_{i}(\theta_{b})$ means sampling at $(i,\theta_{b})$ is less important
although its LDR $G_{i}(\theta_{b})$ may be small. Also note that (14) assigns
a smaller weight when $\theta_{b}$ is in the favorable set of $i^{*}$, which
implies that solving those problems correctly is more important for finding
$i^{*}$.

Let $\BFW$ and $\BFG$ be $k\times B$ matrices whose $(i,b)$th entries are
$W_{i}(\theta_{b})$ and $G_{i}(\theta_{b})$, respectively. For two matrices
$A$ and $B$ of the same size, the Hadamard product, $A\circ B$, is defined as
$(A\circ B)_{ij}=A_{ij}B_{ij}$. For ease of exposition, we define
$0\cdot\infty=\infty$. Thus, (16) is equivalent to that the elements of
$\mathbf{W}\circ\mathbf{G}$ corresponding to all $(i,\theta_{b})\in\Xi$ are
identical, and the rest of elements are all $\infty$.

Recall that the objective of contextual R&S is to minimize the PFS averaged
over the contexts. The optimality conditions for C-OCBA in Gao et al. (2019)
have the same global balance conditions as ours. However, our pairwise balance
conditions are the same as theirs if $W_{i}(\theta_{b})$’s were one for all
$(i,\theta_{b})\in\Xi\cup\Xi^{\text{adv}}$. This makes sense since all
suboptimal pairs are equally important for their objective function. The third
condition in Theorem 4.4 is unique to our problem.

##  5 Sequential Learning Procedures

The optimality conditions in Theorem 4.4 depend on unknown
$\left\\{y_{i}(\theta_{b})\right\\}$. In this section, we present four
sequential sampling algorithms that simultaneously learn
$\left\\{y_{i}(\theta_{b})\right\\}$ and asymptotically optimal allocations
from estimated balance weights along with their convergence properties. To
facilitate the discussion, we define some notation for sample statistics. Let
$\Xi_{n}$, $d_{i,n}$, $\mathbf{W}_{n}$, and $\mathbf{G}_{n}$ be the plug-in
versions of $\Xi$, $d_{i}=d_{i}(\BFM^{*})$, $\mathbf{W}$, and $\BFG$ defined
by replacing $\BFM^{*}$ with $\BFM_{n}$, $i^{*}$ with $i_{n}^{*}$, and $i^{b}$
with $i_{n}^{b}$ for all $b$, respectively. We further denote the estimated
favorable set of $i_{n}^{*}$ by
$\Theta_{n}^{*}:=\\{\theta_{b}:i_{n}^{*}=i_{n}^{b}\\}$ and the fraction of $n$
allocated to $(i,\theta_{b})$ by $\alpha_{i,n}(\theta_{b})$.

Algorithm  ‣ 5 presents the basic structure that all four algorithms follow.
In Step 5, two or more solutions may be tied when determining $i^{*}_{n}$.
Should a tie occur, we may select $i^{*}_{n}$ by applying a tie-breaking rule
such as the first index rule or random sampling rule.

Algorithm 0 Algorithmic framework for selection of the most probable best

1: Given total simulation budget $N$, allocate $n_{0}$ initial replications at
all $(i,\theta_{b})$. Let $n=n_{0}kB$.

2: while $n<N$ do

3: Calculate $\mu_{i,n}(\theta_{b})$ in (5) for all $(i,\theta_{b})$.

4: Inner-level problem : For each $\theta_{b}$, find
$i_{n}^{b}:=\argmin_{1\leq i\leq k}\mu_{i,n}(\theta_{b})$.

5: Outer-level problem : Find $i^{*}_{n}=\argmax_{1\leq i\leq
k}\sum_{b=1}^{B}p_{b}\mathbf{1}\left\\{i=i^{b}_{n}\right\\}.$

6: Make the next sampling decision, $(i,\theta_{b})$, and run a replication at
$(i,\theta_{b})$; $n\leftarrow n+1$.

7: end while

8: return $i_{N}^{*}$ as the estimated MPB.

Our first sequential sampling algorithm, Algorithm 1, applies the plug-in
versions of the balance conditions in Theorem 4.4 to make sampling decisions
in Steps 4–9. In the event of a tie for $i^{*}_{n}$, we apply the same tie-
breaking rule in Algorithm 1 and all other variants: we choose $\argmax_{j\in
i_{n}^{*}}\min_{b:i_{n}^{b}\neq j}G_{j,n}(\theta_{b})$ as $i^{*}_{n}$. Note
that the inner minimization term represents the empirical LDR of Solution $j$
identified as one of the tied $i^{*}_{n}$ incorrectly beating a solution in
its (estimated) adversarial set. Then, the outer maximization problem picks
$j\in i_{n}^{*}$ for which such an event is least likely to occur (i.e, its
adversarial set is most likely to be correct). If the algorithm samples all
$(i,\theta_{b})$ infinitely, the event of a tie happens only finitely many
times with probability $1$ under the uniqueness assumption of $i^{*}$ by the
strong law of large numbers.

Algorithm 1 Sequential Sampling Algorithm for Selection of the MPB

1: Initialization with Step 1 of Algorithm  ‣ 5.

2: while $n<N$ do

3: Run Steps 3–5 of Algorithm  ‣ 5 and update $G_{i,n}(\theta_{b})$ for all
$i\neq i_{n}^{b}$.

4: Construct $\BFW_{n}$ and $\BFG_{n}$. Find
$(i,\theta_{b})=\argmin\BFW_{n}\circ\BFG_{n}$.

5: if
$\Large({N^{n}_{i^{b}_{n}}(\theta_{b})}/{\lambda_{i_{n}^{b}}(\theta_{b})}\Large)^{2}<\sum_{j\neq
i_{n}^{b},i_{n}^{*}}\left({N^{n}_{j}(\theta_{b})}/{\lambda_{j}(\theta_{b})}\right)^{2}$
then

6: Let $j=i_{n}^{b}$.

7: else

8: Let $j=i$.

9: end if

10: Run a replication at $(j,\theta_{b})$; $n\leftarrow n+1$.

11: end while

12: return $i_{N}^{*}$ as the estimated MPB.

Algorithm 1 is easy to implement as we only need to compute $\mathbf{G}_{n}$,
$\mathbf{W}_{n}$, and
$\left(N^{n}_{i}(\theta_{b})/\lambda_{i}(\theta_{b})\right)^{2}$ at each
iteration. Unfortunately, Algorithm 1 does not guarantee that
$N^{n}_{i}(\theta_{b})\rightarrow\infty$ for all $(i,\theta_{b})$ a.s. because
it tends to stop simulating $\left\\{(i_{n}^{*},\theta_{b}):i_{n}^{b}\neq
i_{n}^{*}\right\\}$ once it correctly specifies the minimum number of
$\Theta_{i^{*}}$ needed to distinguish $i_{n}^{*}$ from the second best. For
instance, suppose that $B=50$, $|\Theta_{i^{*}}|=15$, and the second best
solution is the conditional optima at $10$ $\theta_{b}$s. Then, it is enough
to correctly specify $i^{*}$ as the conditional optima at $13$ out of $15$
$\theta_{b}$s in $\Theta_{i^{*}}$ to achieve CS provided that the conditional
optima are correctly specified at all other parameters in
$\Theta^{c}_{i^{*}}$. Even if at two remaining parameters in $\Theta_{i^{*}}$
the second best is incorrectly specified as the conditional optima, CS occurs.
Hence, Algorithm 1 may not fully characterize $\Theta_{i^{*}}$ in the limit
even if $i_{n}^{*}\to i^{*}$. This behavior may also impede empirical
convergence of the algorithm; underestimating $|\Theta_{i^{*}}|$ means that
the algorithm perceives the problem to be harder than it actually is.

The next variant, Algorithm 2, guarantees that $N^{n}_{i^{*}}(\theta_{b})$
increases sublinearly in $n$ for all $\theta_{b}\in\Theta_{i^{*}}^{c}$ so that
$N^{n}_{i^{*}}(\theta_{b})\rightarrow\infty$ and
$\alpha_{i^{*},n}(\theta_{b})\rightarrow 0$ a.s. The only difference from
Algorithm 1 is Step 4: we estimate $\eta_{i^{*}_{n}}(\theta_{b})$ for any
$\theta_{b}$ in the adversarial set of $i^{*}_{n}$ with a posterior sample
from Model (5), not its mean. Define event
$A=\\{\bar{\mu}_{i_{n}^{*},n}(\theta_{b})<\mu_{i_{n}^{b},n}(\theta_{b})\\}$
for such $\theta_{b}$. When $A$ occurs, then the updated $i_{n}^{b}$ in Step 4
must be equal to $i_{n}^{*}$ since all other plug-in estimates at $\theta_{b}$
are unchanged. This increases the estimated preference probability of
$i^{*}_{n}$. Consequently, $i_{n}^{*}$ remains unchanged from Step 3. Note
that the posterior samples are discarded after the current iteration.

Algorithm 2 simulates $(i_{n}^{*},\theta_{b})$ such that
$\theta_{b}\in(\Theta^{*}_{n})^{c}$, if event $A$ occurs and Step 5 selects
$(i_{n}^{*},\theta_{b})$ as the next sampling pair. Hence, the probability of
sampling $(i_{n}^{*},\theta_{b})$ is strictly positive and bounded from above
by , where $\Phi(\cdot)$ is the standard normal CDF. This contrasts with
Algorithm 1, which never simulates such $(i_{n}^{*},\theta_{b})$. Note that
vanishes as $N_{i_{n}^{*}}^{n}(\theta_{b})\rightarrow\infty$. Thus, Algorithm
2 allows sampling $(i_{n}^{*},\theta_{b})$ for
$\theta_{b}\in(\Theta^{*}_{n})^{c}$ recognizing that $\Theta^{*}_{n}$ may not
equal $\Theta_{i^{*}}$ for finite $n$, whereas as $n$ increases, its behavior
becomes increasingly similar to that of Algorithm 1 reckoning that
$\Theta^{*}_{n}$ is estimated more precisely. Formalizing this idea, we show
that Algorithm 2 leads to the optimality conditions in Theorem 4.4 as stated
below.

Algorithm 2 Hybrid of Posterior Sample for Selection of the MPB

1: Initialization with Step 1 of Algorithm  ‣ 5.

2: while $n<N$ do

3: Run Step 3 of Algorithm 1.

4: For all $\theta_{b}\in\left(\Theta_{n}^{*}\right)^{c}$, replace
$\mu_{i^{*}_{n},n}(\theta_{b})$ with $\bar{\mu}_{i_{n}^{*},n}(\theta_{b})\sim
N(\mu_{i^{*}_{n},n}(\theta_{b}),\sigma_{i^{*}_{n},n}^{2}(\theta_{b}))$ in (5);
update $i_{n}^{b}$.

5: Run Steps 4–10 of Algorithm 1.

6: end while

7: return $i_{N}^{*}$ as the estimated MPB.

###### Theorem 5.1

The sampling allocations made by Algorithm 2 satisfy the followings a.s.:

$\lim_{n\rightarrow\infty}\Big{\\{}\left(\tfrac{\alpha_{i^{b},n}(\theta_{b})}{\lambda_{i^{b}}(\theta_{b})}\right)^{2}-\sum\nolimits_{j\neq
i^{b},i^{*}}\left(\tfrac{\alpha_{j,n}(\theta_{b})}{\lambda_{j}(\theta_{b})}\right)^{2}\Big{\\}}=0,\forall\theta_{b}\mbox{
and
}\lim_{n\rightarrow\infty}\tfrac{W_{i}(\theta_{b})G_{i,n}(\theta_{b})}{W_{j}(\theta_{c})G_{j,n}(\theta_{c})}=1,\forall(i,\theta_{b})\neq(j,\theta_{c})\in\Xi.$

Also,
$\limsup_{n\rightarrow\infty}N^{n}_{i}(\theta_{b})/N^{n}_{j}(\theta_{c})<\infty$
for all $(i,\theta_{b}),(j,\theta_{c})\in\left(\Xi^{\textup{adv}}\right)^{c}$
and $\lim_{n\rightarrow\infty}N^{n}_{i^{*}}(\theta_{b})/n=0$ for all
$(i^{*},\theta_{b})\in\Xi^{\textup{adv}}$.

###  5.1 Algorithms for Learning Both MPB and its Favorable Set

As mentioned in Section 1, the preference probability of the MPB can be used
as a measure of how much model risk input uncertainty poses for the underlying
R&S problem. This motivates us to devise algorithms that learn both $i^{*}$
and $\Theta_{i^{*}}$ correctly. Assuming $i^{*}$ is learned correctly,
identifying $\Theta_{i^{*}}$ can be viewed as a binary classification task for
each $\theta_{b}$ to determine whether $\theta_{b}\in\Theta_{i^{*}}$ or not.
From this viewpoint, we can construct a confusion matrix for identifying
$\Theta_{i^{*}}$ as follows:

$\blockarray{ccc}\text{Actual $\backslash$
Predicted}&\theta_{b}\in\Theta_{n}^{*}\theta_{b}\in(\Theta_{n}^{*})^{c}\\\
\block{c(cc)}\theta_{b}\in\Theta_{i^{*}}\text{True Positive}\text{False
Negative}\\\ \theta_{b}\in\Theta_{i^{*}}^{c}\text{False Positive}\text{True
Negative}\\\ $ (17)

From (17), we define the _false negative rate_ (FNR) and _accuracy_ (ACC) as

$\text{FNR}:=\frac{\prob{}_{\pi}\left(\Theta_{i^{*}}\cap\left(\Theta^{*}_{n}\right)^{c}\right)}{\prob{}_{\pi}\left(\Theta_{i^{*}}\right)};\hskip
15.0pt\text{ACC}:=$ (18)

ACC is maximized ($=1$) when $\Theta_{n}^{*}=\Theta_{i^{*}}$, i.e., the
estimated preference probability of $i^{*}$ is exact. This motivates us to
design Algorithm 3 that maximizes the probability of event
$\text{ACC}_{*}:=\\{i_{n}^{*}=i^{*}\\}\cap\\{\Theta_{n}^{*}=\Theta_{i^{*}}\\}$.
Recall that we apply the same tie-breaking rule as in Algorithm 1 if
$i_{n}^{*}$ is not a singleton. ACC penalizes both false negative and false
positive. Meanwhile, recall that underestimation of $\Theta_{i^{*}}$ (and
preference probability of $i^{*}$) causes performance degradation of Algorithm
1. Thus, we also consider Algorithm 4 that penalizes only false negatives by
minimizing FNR. Since FNR$=0$ when
$\\{\Theta^{*}_{n}\supseteq\Theta_{i^{*}}\\}$, Algorithm 4 maximizes the
probability of event
$\text{FN}_{*}:=\\{i_{n}^{*}=i^{*}\\}\cap\\{\Theta^{*}_{n}\supseteq\Theta_{i^{*}}\\}.$
In the following, we discuss the optimality conditions for the OCBA
formulation when each of and is adopted in lieu of PFS.

First, the LDR of $\text{ACC}_{*}^{c}$ can be bounded from below as

$\displaystyle\liminf\nolimits_{n\rightarrow\infty}-\tfrac{1}{n}\log$
$\displaystyle=$
$\displaystyle\liminf\nolimits_{n\rightarrow\infty}-\tfrac{1}{n}\log$ (19)
$\displaystyle\geq$
$\displaystyle\min\left\\{\underline{\mathrm{LDR}},\liminf\nolimits_{n\rightarrow\infty}-\tfrac{1}{n}\log\right\\}.$

Thus, it only remains to find the LDR of . To guarantee
$\left\\{\Theta^{*}_{n}=\Theta_{i^{*}}\right\\}$, we need to ensure (a) for
each $\theta_{b}$ in $\Theta_{i^{*}}$, the conditional optimum is correctly
identified to be $i^{*}$; and (b) for each $\theta_{b}$ in
$\Theta_{i^{*}}^{c}$, $i^{*}$ is not chosen as the conditional optimum.
Consequently, we have

$\liminf_{n\rightarrow\infty}-\tfrac{1}{n}\log$ (20)

Combining (19) with (20), $\liminf_{n\rightarrow\infty}-\frac{1}{n}\log$ is
bounded from below by

$\underline{\mathrm{LDR}}_{\textup{ACC}}:=\min\left\\{\underline{\mathrm{LDR}},\min\nolimits_{\left\\{i\neq
i^{*},\theta_{b}\in\Theta_{i^{*}}\right\\}}G_{i}(\theta_{b}),\min\nolimits_{\theta_{b}\in\Theta_{i^{*}}^{c}}G_{i^{*}}(\theta_{b})\right\\}.$

Then, the OCBA for Algorithm 3 can be formulated as

$\displaystyle\underline{\mathrm{LDR}}_{\textup{ACC}}^{*}:=\max_{\bm{\alpha}}\underline{\mathrm{LDR}}_{\textup{ACC}}\;\;\textrm{subject
to}\;\;\sum\nolimits_{i=1}^{k}\sum\nolimits_{b=1}^{B}\alpha_{i}(\theta_{b})=1,\alpha_{i}(\theta_{b})\geq
0,1\leq i\leq k,1\leq b\leq B.$ (21)

To maximize , note that
$\left\\{\Theta_{i^{*}}\subseteq\Theta^{*}_{n}\right\\}=\left\\{\mu_{i^{b}_{n},n}(\theta_{b})\leq\mu_{i^{*},n}(\theta_{b}),\forall\theta_{b}\in(\Theta^{*}_{n})^{c}\right\\}$.
Thus, the LDR of $\textup{FN}_{*}^{c}$ is bounded from below as

$\displaystyle\liminf_{n\rightarrow\infty}-\tfrac{1}{n}\log\prob\left(\textup{FN}_{*}^{c}\right)$
$\displaystyle=\liminf_{n\rightarrow\infty}-\tfrac{1}{n}\log$
$\displaystyle\geq\min\left\\{\underline{\mathrm{LDR}},\min\nolimits_{\theta_{b}\in\Theta^{c}_{i^{*}}}G_{i^{*}}(\theta_{b})\right\\}.$

Defining
$\underline{\mathrm{LDR}}_{\textup{FN}}:=\min\left\\{\underline{\mathrm{LDR}},\min_{\theta_{b}\in\Theta^{c}_{i^{*}}}G_{i^{*}}(\theta_{b})\right\\}$,
the OCBA for Algorithm 4 can be written as

$\displaystyle\underline{\mathrm{LDR}}_{\textup{FN}}^{*}:=\max_{\bm{\alpha}}\underline{\mathrm{LDR}}_{\textup{FN}}\;\;\textrm{subject
to}\;\;\sum\nolimits_{i=1}^{k}\sum\nolimits_{b=1}^{B}\alpha_{i}(\theta_{b})=1,\alpha_{i}(\theta_{b})\geq
0,1\leq i\leq k,1\leq b\leq B.$ (22)

Since $\underline{\mathrm{LDR}}_{\textup{ACC}}$ and
$\underline{\mathrm{LDR}}_{\textup{FN}}$ are concave in $\bm{\alpha}$, (21)
and (22) are concave programs. Therefore, we can derive their optimality
conditions from the KKT conditions as summarized in Theorem 5.2.

######  Theorem 5.2 (Optimality conditions for (21) and (22))

Any ${\bm{\alpha}}$ is optimal for (21) (respectively (22)), if and only if,
${\bm{\alpha}}$ satisfies the following system of equations:

  * •

(Global balance condition) For all $1\leq b\leq B$,
$\;{\alpha}^{2}_{i^{b}}(\theta_{b})/{\lambda^{2}_{i^{b}}(\theta_{b})}=\sum_{i\neq
i^{b}}{\alpha}^{2}_{i}(\theta_{b})/{\lambda^{2}_{i}(\theta_{b})}.$

  * •

(Pairwise balance condition) For all
$(i,\theta_{b}),(j,\theta_{b^{\prime}})\in\Xi\cup\Xi^{\textup{adv}}$,
$W^{\textup{ACC}}_{i}(\theta_{b})G_{i}(\theta_{b})=W^{\textup{ACC}}_{j}(\theta_{b^{\prime}})G_{j}(\theta_{{b^{\prime}}})$
(resp.
$W^{\textup{FN}}_{i}(\theta_{b})G_{i}(\theta_{b})=W^{\textup{FN}}_{j}(\theta_{b^{\prime}})G_{j}(\theta_{{b^{\prime}}})$),
where

$W^{\textup{ACC}}_{i}(\theta_{b}):=\begin{cases}1,&\mbox{if
}\theta_{b}\in\Theta_{i^{*}}\mbox{ or }(i,\theta_{b})\in\Xi^{\textup{adv}}\\\
W_{i}(\theta_{b}),&\mbox{if }\theta_{b}\notin\Theta_{i^{*}}\end{cases},\hskip
15.0ptW^{\textup{FN}}_{i}(\theta_{b})=\begin{cases}W_{i}(\theta_{b}),&\mbox{if
}(i,\theta_{b})\in\Xi\\\ 1,&\mbox{if
}(i,\theta_{b})\in\Xi^{\textup{adv}}\end{cases},$ (23)

and $W^{\textup{ACC}}_{i}(\theta_{b})=W^{\textup{FN}}_{i}(\theta_{b})=\infty$
for all $(i,\theta_{b})\in\left(\Xi\cup\Xi^{\textup{adv}}\right)^{c}$.

Let $\BFW_{n}^{\textup{ACC}}$ and $\BFW_{n}^{\textup{FN}}$ be the sample
version of $\BFW^{\textup{ACC}}$ and $\BFW^{\textup{FN}}$ after $n$
replications are made. Algorithm 3 (Algorithm 4) is identical to Algorithm 1
except that $\mathbf{W}_{n}$ in Step 3 and $\sum_{j\neq
i_{n}^{b},i_{n}^{*}}\left({N^{n}_{j}(\theta_{b})}/{\lambda_{j}(\theta_{b})}\right)^{2}$
in Step 5 are replaced with $\mathbf{W}^{\textup{ACC}}_{n}$
($\mathbf{W}^{\textup{FN}}_{n}$) and $\sum_{j\neq
i_{n}^{b}}\left({N^{n}_{j}(\theta_{b})}/{\lambda_{j}(\theta_{b})}\right)^{2}$,
respectively. Theorem 5.3 states that Algorithm 3 and 4 indeed satisfy the
optimality conditions in Theorem 5.2, asymptotically.

###### Theorem 5.3

The sampling allocations made by Algorithm 3 satisfy the followings a.s.:

$\lim_{n\rightarrow\infty}\Big{\\{}\left(\tfrac{\alpha_{i^{b},n}(\theta_{b})}{\lambda_{i^{b}}(\theta_{b})}\right)^{2}-\sum_{j\neq
i^{b}}\left(\tfrac{\alpha_{j,n}(\theta_{b})}{\lambda_{j}(\theta_{b})}\right)^{2}\Big{\\}}=0,\forall\theta_{b}\mbox{
and
}\lim_{n\rightarrow\infty}\tfrac{W^{\textup{ACC}}_{i}(\theta_{b})G_{i,n}(\theta_{b})}{W^{\textup{ACC}}_{j}(\theta_{c})G_{j,n}(\theta_{c})}=1,\forall(i,\theta_{b}),(j,\theta_{c})\in\Xi\cup\Xi^{\textup{adv}}.$

Also,
$\limsup_{n\rightarrow\infty}{N^{n}_{i}(\theta_{b})}/{N^{n}_{j}(\theta_{b^{\prime}})}<\infty$
for all $(i,\theta_{b})\neq(j,\theta_{b^{\prime}})$. For Algorithm 4, the same
statements hold when $W^{\textup{ACC}}_{i}(\theta_{b})$ is replaced with
$W^{\textup{FN}}_{i}(\theta_{b})$.

Note that
$\limsup_{n\rightarrow\infty}{N^{n}_{i}(\theta_{b})}/{N^{n}_{j}(\theta_{b^{\prime}})}<\infty$
for all $(i,\theta_{b})\neq(j,\theta_{b^{\prime}})$ implies that Algorithms 3
and 4 distribute $n$ to each $(i,\theta_{b})$ at a linear rate. This contrasts
with Algorithm 2, which allocates sampling effort sublinear in $n$ to
$(i,\theta_{b})$ in $\Xi^{\text{adv}}$.

##  6 Empirical Analysis

In this section, we present numerical performance analyses using synthetic and
market simulation examples; the latter is a simplified version of General
Motors’ vehicle content optimization problem (Song et al. 2020). We compare
Algorithms 1–4 with Equal Allocation (EA) that samples
$\left\\{(i,\theta_{b})\right\\}$ uniformly given a budget and C-OCBA proposed
by Gao et al. (2019) to solve a contextual R&S problem. For all algorithms,
initial sample size $n_{0}=5$ is adopted for the synthetic example and
$n_{0}=10$ for the market simulation example. We evaluate four metrics for the
performance analysis: (i) PFS: ; (ii) FNR; and (iii) $1-$ACC, where FNR and
ACC are defined in (18). All metrics are reported in the base-10 logarithmic
scale.

###  6.1 Synthetic example

We consider a set of synthetic examples to test robustness of the six
algorithms in comparison under different problem configurations. We adopt
$k=10$ and $B=50$, and assume $p_{b}=1/B$ for all $1\leq b\leq 50$ for all
experiments; the unequal probability case can be found in Appendix 12. Both
cases exhibit the same conclusion. For each macro run, we randomize the
simulation output standard deviations, $\\{\lambda_{i}(\theta_{b})\\}$, by
sampling them from uniform distribution $U[4,6]$.

The conditional optimum at each $\theta_{b}$ is set to be

$i^{b}=\ell\;\;(\mbox{if }5\ell-4\leq b\leq 5\ell,\mbox{ for some
}1\leq\ell\leq 7),\;\;i^{b}=8\;\;(\mbox{if }36\leq b\leq
41),\;\;i^{b}=10\;\;(\mbox{if }42\leq b\leq 50).$ (24)

This implies $i^{*}=10$ with preference probability $0.18$. The second best is
$i=8$ whose preference probability is $0.12$. All solutions except for $i=9$
are conditional optima at some $\theta_{b}$. For each macro run, we set
$y_{i^{b}}(\theta_{b})=1$ and fill in
$\left\\{y_{i}(\theta_{b})\right\\}_{1\leq i\leq k,i\neq i^{b}}$ randomly
without replacement from $\left\\{2,3,\cdots,k\right\\}$ for each $1\leq b\leq
50$. In the baseline experiment, the simulation output distribution is normal
for all $(i,\theta_{b})$ and their variances are assumed to be known.

Three modifications to the baseline instance are included in the experiment,
each of which is designed to examine the algorithms’ sensitivity to a problem
feature:

  * •

Scenario 1 (more dominant $i^{*}$): $i^{*}=10$ is made more dominant by
modifying (24) so that $i^{b}=10$ for all $36\leq b\leq 50$. This makes any
$1\leq i\leq 7$ be the second best and the difference in preference
probability compared to $i^{*}$ is $0.2$, which is increased from $0.06$ in
the baseline.

  * •

Scenario 2 (larger simulation error variance): We draw
$\lambda_{i}(\theta_{b})$ from $U[8,12]$ for each macro run. In other words,
we amplify the simulation variance twice as much as in the baseline on
average.

  * •

Scenario 3 (non-normal simulation output): To test robustness when the
normality assumption fails to hold, the simulation output is set to follow a
shifted exponential distribution, i.e.,
$Y_{i}(\theta_{b})\overset{D}{=}\mu_{i}(\theta_{b})-\lambda_{i}(\theta_{b})+\text{Exp}(\lambda_{i}(\theta_{b}))$,
where $\overset{D}{=}$ means equivalence in distribution and
$\lambda_{i}(\theta)$ is the mean of the exponential distribution.

Figure 1: Performance measure estimates from 5,000 macro runs for each
scenario. The vertical dashed line in each plot represents the end of the
initial sampling phase.

Figure 1 summarizes numerical results for all scenarios. The baseline is
presented in the first column. Notice that EA and C-OCBA show slower
convergence in PFS than our algorithms, where the gaps widen as the simulation
budget increases. Algorithm 1 falls behind Algorithms 3 and 4 in PFS. As
mentioned in Section 5, the shortfall is caused by that Algorithm 1
underestimates the MPB’s preference probability. This can be confirmed from
the FNR plot. Algorithm 2 corrects this behavior and achieves the fastest PFS
convergence rate while showing good performances in FNR and $1-$ACC. On the
other hand, Algorithms 3 and 4 exhibit the best FNR and $1-$ACC, respectively,
as they are designed to maximize the respective measures.

For Scenario 1, finding the MPB is expected to be easier than the baseline.
While all algorithms show improved PFS, notice that C-COBA outperforms
Algorithm 1. This is because Algorithm 1’s underestimation of the MPB’s
preference probability is more severe when $i^{*}$ is dominant. Algorithms 2,
4, and 3 continue to show the best performances in PFS, FNR, and $1-$ACC,
respectively.

Scenario 2 has higher simulation variances, which explains why all algorithms
show worse performances than the baseline given the same budget. The overall
trends across the algorithms are similar to the baseline. Finally, the last
column (Scenario 3) shows that even if the simulation error is non-Gaussian,
similar trends can be observed among all other algorithms while Algorithm 2
clearly outperforms in all three metrics.

###  6.2 Market simulation for product portfolio optimization

We consider a product portfolio optimization example formulated as a R&S
problem. Each solution represents a product mix (portfolio) to be sold at a
store and the objective is to maximize the expected sales. A product has
several features and each feature is offered at different levels. For
instance, consider a vehicle model with features such as seat type, engine,
power train, and etc. For seats, there may be two levels for the material:
fabric or leather. In our example, each product has five features and each
feature is offered in three levels. Consequently, there are total $3^{5}=243$
unique products. The marginal price of the $\ell$th level of the $m$th feature
is $p_{m\ell}=50\ell m+10\ell^{2}$. The price of a product is the sum of the
marginal prices of all features of the product.

We assume that the store’s capacity allows $40$ products to be displayed and
the product portfolio stipulates which products fill those slots. Since not
all products are equally popular, predicting what combinations of feature-
levels are most desired by the customers is critical for maximizing the sales.
Given a product portfolio, the market simulator simulates customers’ arrivals
to the store and purchase decisions; aggregating them, the expected sales is
estimated. Although the store capacity is $40$, each customer may not see all
$40$ products due to temporary stock-outs. To simulate this, we generate
$M_{n}\sim$ Poisson($4$) and sample $M_{n}$ products from a given portfolio to
construct a set of products that the $n$th customer compares before making a
purchase decision. The customer may choose to purchase one of the products in
the set or none. Given this choice set, a customer’s decision is simulated
based on a discrete choice model (Train 2009) as we describe briefly below.

Suppose the $n$th customer’s simulated choice set has size $J$; first $J-1$
elements are products and the last is no-purchase. Let $K$ be the number of
feature-levels plus one, i.e., $K=5\times 3+1=16$. The $j$th product in the
choice set can be represented by
$x_{nj}=(0,1,\cdots,0,-p_{j})\in\mathbb{R}^{K}$, where each of the first $K-1$
elements is an indicator denoting whether the corresponding feature-level is
included in the product. The last element corresponds to the negative price of
the $j$th product. For the no-purchase option, $x_{nJ}$ is a vector of zeroes.
The following linear model to represent the $n$th customer’s utility for the
$j$th option in the choice set:

$U_{nj}=\beta_{n}^{\top}x_{nj}+\epsilon_{nj},\mbox{ where
}\\{\epsilon_{nj}\\}\;\;\text{is i.i.d and}\;\;$ (25)

where $\beta_{n}$ is a $K$-dimensional vector that represents the $n$th
customer’s utilities for $(K-1)$ feature-levels and price, and $\epsilon_{nj}$
is random noise that cannot be captured by the linear model. Let $y_{n}$
represent the $n$th customer’s choice. Under Model (25), the probability that
the $n$th customer selects the $j$th option can be derived as for $1\leq j\leq
J$. Once the choice set for the $n$th simulated customer is given, the
expected sales for the customer is equal to the probability of purchasing any
products in the choice set,
$1-\\{\sum_{j=1}^{J}\exp(\beta_{n}^{\top}x_{nj})\\}^{-1}$. Let $N$ be the
total number of consumers we simulate in each replication. Then, the expected
sales can be estimated by sampling $N$ choice sets and computing:
$\sum_{n=1}^{N}(1-\\{\sum_{j=1}^{J}\exp(\beta_{n}^{\top}x_{nj})\\}^{-1})$.
Hence, the key input parameter for this simulation model is the utility
vector, $\\{\beta_{n}\\}$.

We adopt a Bayesian conjoint analysis to model uncertainty about $\beta_{n}$
from synthetic survey data. For details of the model, data generation and
posterior sampling algorithm of $\beta_{n}$, see 13.

Let $\mathcal{B}=\left\\{\beta_{1},\beta_{2},\cdots,\beta_{N}\right\\}$, where
each $\beta_{n}$ is sampled from its posterior distribution. Thus,
$\mathcal{B}$ represents the set of all input vectors required to set up a
simulation model. To construct the finite support for $\mathcal{B}$, we
generate $\mathcal{B}_{1},\mathcal{B}_{2},\cdots,\mathcal{B}_{B}$ from the
posterior sampling algorithm for $\beta_{n}$.

Table 1: Portfolio description based on the utilities at prior mean.

Portfolio | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---|---|---|---|---
Description | Top 40 | Top 20$\times 2$ | Top 10$\times 4$ | Top 5$\times 8$ | Top 4$\times 10$ |  | Top 5$\times 3$
---
+Top 6$\sim$30
| Top 10$\times 2$
---
+Top 11$\sim$30
| Top 10$\times 3$
---
+Top 11$\sim$20
| Top 4$\times 6$
---
+Top 5$\sim$20
Preference Prob. | 0 | 0 | 0.30 | 0.44 | 0.02 | 0.08 | 0 | 0.04 | 0.12

Table 1 presents nine portfolios constructed by first sorting all $243$
products by utilities computed from the prior mean for $\beta_{n}$. The
rationale for duplicating the popular items is to ensure they are less likely
to be out of stock. Since the exact expected sales cannot be computed
analytically, we report Monte Carlo estimates in Table 4. The preference
probabilities in Table 1 are computed from Table 4. Observe that the MPB is
Portfolio 4 whose preference probability is $0.44$ followed by Portfolio 3.
Since their difference in preference probabilities is $0.14$, a risk-sensitive
decision maker may prefer to reduce the posterior uncertainty by increasing
the survey data size. Interestingly, if one takes the risk-neutral view and
evaluates each portfolio by its mean sales averaged over all $\mathcal{B}$s,
then Portfolio 3 is the best. Furthermore, if one takes the set of sampled
$\mathcal{B}$ as the ambiguity set and adopts the distributionally robust
approach, then the best is still Portfolio 3. In neither approach, Portfolio 4
is chosen as the best, yet it has the largest preference probability. On the
other hand, our formulation can also detect Portfolio 3 as a closely
contending solution.

Figure 2: Numerical results of market portfolio optimization averaged over
5,000 macro runs; each replication simulates $N=20$ consumers.

Finally, Figure 2 shows the performances of the six algorithms applied to this
example. Unlike the synthetic examples in Section 6.1, the simulation error
variances are unknown in this case; thus, Model 26 in the Electronic Companion
is adopted instead of Model (4). In addition, a batch of $10$ replications is
run every time a sampling decision is made. From Figure 2, observe that
Algorithm 2 clearly outperforms the other algorithms in PFS. In terms of FNR
and $1-$ACC, Algorithms 4 and 3 exhibit best performances, respectively. This
reaffirms the conclusion in Section 6.1.

##  7 Conclusions

We introduced the MPB as a solution approach for a R&S problem subject to
input uncertainty. From the theoretical aspect, we have established the LDR of
PFS for the MPB selection problem. A lower bound of the true LDR has been
derived to retain interpretable optimality conditions for the asymptotic
sampling ratios, which involve balance weights that measure the importance of
sampling each solution-parameter pair to reduce the FS rate. We have also
proposed several sequential sampling algorithms some of which are shown to
achieve a set of optimality conditions in the limit. Our numerical experiments
support the efficiency of the proposed algorithms.

Further study is required to create a DP-based algorithm to find the MPB. We
will also consider more general utility functions other than the identity
function adopted in this work to define the preference probability. Moreover,
when a user can provide the largest negligible optimality gap, the preference
probability may be modified to reflect this information to control the level
of conservatism. Incorporating a continuous support for the input parameter is
another remaining research question.

## References

  * Avci et al. (2021) Avci H, Nelson BL, Waechter A (2021) Getting to “rate-optimal” in ranking & selection. Kim S, Feng B, Smith K, Masoud S, Zheng Z, eds., _Proceedings of the 2021 Winter Simulation Conference_ (Piscataway, New Jersey: IEEE).
  * Barton et al. (2008) Barton GR, Briggs AH, Fenwick EAL (2008) Optimal cost-effectiveness decisions:the role of the cost-effectiveness acceptability curve (ceac), the cost-effectiveness acceptability frontier (ceaf), and the expected value of perfection information (evpi). _Value in Health_ 11(5):886–897.
  * Bechhofer et al. (1959) Bechhofer RE, Elmaghraby S, Morse N (1959) A single-sample multiple-decision procedure for selecting the multinomial event which has the highest probability. _The Annals of Mathematical Statistics_ 102–119.
  * Berger (1985) Berger JO (1985) _Statistical Decision Theory and Bayesian Analysis_ (Springer Science & Business Media).
  * Bishop (2006) Bishop CM (2006) Pattern recognition. _Machine learning_ 128(9).
  * Chen et al. (2008) Chen CH, He D, Fu M, Lee LH (2008) Efficient simulation budget allocation for selecting an optimal subset. _INFORMS Journal on Computing_ 20(4):579–595.
  * Chen et al. (2000) Chen CH, Lin J, Yücesan E, Chick SE (2000) Simulation budget allocation for further enhancing the efficiency of ordinal optimization. _Discrete Event Dynamic Systems_ 10(3):251–270.
  * Chen (1992) Chen P (1992) Truncated selection procedures for the most probable event and the least probable event. _Annals of the Institute of Statistical Mathematics_ 44(4):613–622.
  * Chen and Ryzhov (2019a) Chen Y, Ryzhov IO (2019a) Balancing optimal large deviations in ranking and selection, working paper.
  * Chen and Ryzhov (2019b) Chen Y, Ryzhov IO (2019b) Complete expected improvement converges to an optimal budget allocation. _Advances in Applied Probability_ 51(1):209–235.
  * Chick and Inoue (2001) Chick SE, Inoue K (2001) New two-stage and sequential procedures for selecting the best simulated system. _Operations Research_ 49(5):732–743.
  * Corlu and Biller (2013) Corlu C, Biller B (2013) A subset selection procedure under input parameter uncertainty. Pasupathy R, Kim SH, Tolk A, Hill R, Kuhl ME, eds., _Proceedings of the 2013 Winter Simulation Conference_ , 463–473 (Piscataway, New Jersey: IEEE).
  * Corlu and Biller (2015) Corlu CG, Biller B (2015) Subset selection for simulations accounting for input uncertainty. Yilmaz L, Chan WKV, I Moon TMKR, Macal C, Rossetti MD, eds., _Proceedings of the 2015 Winter Simulation Conference_ , 437–446 (Piscataway, New Jersey: IEEE).
  * Daskin et al. (1997) Daskin MS, Hesse SM, Revelle CS (1997) $\alpha$-reliable p-minimax regret: A new model for strategic facility location modeling. _Location science_ 5(4):227–246.
  * Durrett (2019) Durrett R (2019) _Probability: theory and examples_ , volume 49 (Cambridge university press).
  * Fan et al. (2020) Fan W, Hong LJ, Zhang X (2020) Distributionally robust selection of the best. _Management Science_ 66(1):190–208.
  * Fenwick et al. (2001) Fenwick E, Claxton K, Sculpher M (2001) Representing uncertainty: the role of cost-effectiveness acceptability curves. _Health Economics_ 10(8):779–787.
  * Frazier et al. (2008) Frazier PI, Powell WB, Dayanik S (2008) A knowledge-gradient policy for sequential information collection. _SIAM Journal on Control and Optimization_ 47(5):2410–2439.
  * Fu et al. (2007) Fu MC, Hu JQ, Chen CH, Xiong X (2007) Simulation allocation for determining the best design in the presence of correlated sampling. _INFORMS Journal on Computing_ 19(1):101–111.
  * Gao et al. (2019) Gao S, Du J, Chen CH (2019) Selecting the optimal system design under covariates. _2019 IEEE 15th International Conference on Automation Science and Engineering (CASE)_ , 547–552 (Piscataway, New Jersey: IEEE).
  * Gao et al. (2017) Gao S, Xiao H, Zhou E, Chen W (2017) Robust ranking and selection with optimal computing budget allocation. _Automatica_ 81:30–36.
  * Glynn and Juneja (2004) Glynn P, Juneja S (2004) A large deviations perspective on ordinal optimization. Ingalls RG, Rossetti MD, Smith JS, Peters BA, eds., _Proceedings of the 2004 Winter Simulation Conference_ , 577–585 (Piscataway, New Jersey: IEEE).
  * Harvey et al. (2010) Harvey CR, Liechty JC, Liechty MW, Müller P (2010) Portfolio selection with higher moments. _Quantitative Finance_ 10(5):469–485.
  * He et al. (2007) He D, Chick SE, Chen CH (2007) Opportunity cost and ocba selection procedures in ordinal optimization for a fixed number of alternative systems. _IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews)_ 37(5):951–961.
  * Hirano (2010) Hirano K (2010) Decision theory in econometrics. _Microeconometrics_ , 29–35 (Springer).
  * Kim et al. (2021) Kim KK, Kim T, Song E (2021) Selection of the most probable best under input uncertainty. Kim S, Feng B, Smith K, Masoud S, Zheng Z, eds., _Proceedings of the 2021 Winter Simulation Conference_ (Piscataway, New Jersey: IEEE).
  * Kim and Nelson (2006) Kim SH, Nelson BL (2006) Selecting the best system. Henderson SG, Nelson BL, eds., _Simulation_ , volume 13 of _Handbooks in Operations Research and Management Science_ , 501–534 (Elsevier).
  * Li et al. (2020) Li H, Lam H, Liang Z, Peng Y (2020) Context-dependent ranking and selection under a Bayesian framework. Bae KHG, Feng B, Kim S, Lazarova-Molnar S, Zheng Z, Roeder T, Thiesing R, eds., _Proceedings of the 2020 Winter Simulation Conference_ , 2060–2070 (Piscataway, New Jersey: IEEE).
  * Martello and Toth (1990) Martello S, Toth P (1990) _Knapsack problems: algorithms and computer implementations_ (John Wiley & Sons, Inc.).
  * Pearce and Branke (2017) Pearce M, Branke J (2017) Efficient expected improvement estimation for continuous multiple ranking and selection. Chan WKV, D’Ambrogio A, Zacharewicz G, Mustafee N, Wainer G, Page E, eds., _Proceedings of the 2017 Winter Simulation Conference_ , 2161–2172 (Piscataway, New Jersey: IEEE).
  * Peng et al. (2016) Peng Y, Chen CH, Fu MC, Hu JQ (2016) Dynamic sampling allocation and design selection. _INFORMS Journal on Computing_ 28(2):195–208.
  * Peng et al. (2018) Peng Y, Chong EK, Chen CH, Fu MC (2018) Ranking and selection as stochastic control. _IEEE Transactions on Automatic Control_ 63(8):2359–2373.
  * Peng and Fu (2017) Peng Y, Fu M (2017) Myopic allocation policy with asymptotically optimal sampling rate. _IEEE Transactions on Automatic Control_ 62(4):2041–2047.
  * Powell and Ryzhov (2012) Powell WB, Ryzhov IO (2012) _Optimal learning_ , volume 841 (John Wiley & Sons).
  * Robert (2007) Robert C (2007) _The Bayesian choice: from decision-theoretic foundations to computational implementation_ (Springer Science & Business Media).
  * Rossi et al. (2012) Rossi PE, Allenby GM, McCulloch R (2012) _Bayesian statistics and marketing_ (John Wiley & Sons).
  * Ryzhov (2016) Ryzhov IO (2016) On the convergence rates of expected improvement methods. _Operations Research_ 64(6):1515–1528.
  * Salemi et al. (2019) Salemi P, Song E, Nelson BL, Staum J (2019) Gaussian markov random fields for discrete optimization via simulation: Framework and algorithms. _Operations Research_ 67(1):250–266.
  * Shen et al. (2021) Shen H, Hong LJ, Zhang X (2021) Ranking and selection with covariates for personalized decision making. _INFORMS Journal on Computing_ 33(4):1500–1519.
  * Song and Nelson (2019) Song E, Nelson BL (2019) Input–output uncertainty comparisons for discrete optimization via simulation. _Operations Research_ 67(2):562–576.
  * Song et al. (2015) Song E, Nelson BL, Hong LJ (2015) Input uncertainty and indifference-zone ranking & selection. Yilmaz L, Chan WKV, I Moon TMKR, Macal C, Rossetti MD, eds., _Proceedings of the 2015 Winter Simulation Conference_ , 414–424 (Piscataway, New Jersey: IEEE).
  * Song et al. (2020) Song E, Wu-Smith P, Nelson BL (2020) Uncertainty quantification in vehicle content optimization for general motors. _INFORMS Journal on Applied Analytics_ 50(4):225–238.
  * Tollefson et al. (2014) Tollefson E, Goldsman D, Kleywegt A, Tovey C (2014) Optimal selection of the most probable multinomial alternative. _Sequential Analysis_ 33(4):491–508.
  * Train (2009) Train KE (2009) _Discrete choice methods with simulation_ (Cambridge university press).
  * Ungredda et al. (2022) Ungredda J, Pearce M, Branke J (2022) Bayesian optimisation vs. input uncertainty reduction. _ACM Transactions on Modeling and Computer Simulation_ In press.
  * Wu and Zhou (2019) Wu D, Zhou E (2019) Ranking and selection under input uncertainty: Fixed confidence and fixed budget, https://arxiv.org/abs/1708.08526.
  * Wu et al. (2018) Wu D, Zhu H, Zhou E (2018) A Bayesian risk approach to data-driven stochastic optimization: formulations and asymptotics. _SIAM Journal on Optimization_ 28(2):1588–1612.
  * Xie and Zhou (2015) Xie W, Zhou E (2015) Simulation optimization when facing input uncertainty. Yilmaz L, Chan WKV, Moon I, Roeder TMK, Macal C, Rossetti MD, eds., _Proceedings of the 2015 Winter Simulation Conference_ , 3714–3724 (Piscataway, New Jersey: IEEE).
  * Zhu et al. (2020) Zhu H, Liu T, Zhou E (2020) Risk quantification in stochastic simulation under input uncertainty. _ACM Transactions on Modeling and Computer Simulation_ 30(1):1–24.

Electronic Companion of “Selection of the Most Probable Best” This electronic
companion provides technical proofs and some experiment details on market
simulation not presented in the main body of the paper.

##  8 Bayesian Learning Model for Unknown Variance

In this section, we consider the extension of Bayesian model in (4) to the
unknown variance case. Suppose both mean and variance of simulation model are
unknown in advance. The normal-gamma conjugate pair is frequently adopted in
the Bayesian learning literature (e.g., see Section 5 in Ryzhov (2016) and
Section 2.3.5 in Powell and Ryzhov (2012)) to model the case of unknown
simulation variance. The prior of this model can be written as

$\displaystyle Y_{i}(\theta)$ $\displaystyle\sim
N(\eta_{i}(\theta),\lambda^{2}_{i}(\theta)),$ (26)
$\displaystyle(\eta_{i}(\theta),\lambda^{-2}_{i}(\theta))$ $\displaystyle\sim
NG({\mu}_{i,0}(\theta),\tau_{i,0}(\theta),\gamma_{i,0}(\theta),\omega_{i,0}(\theta)),$

where $NG$ represents the normal-gamma distribution. From (26), one can derive
$\lambda^{-2}_{i}(\theta)\sim\text{Gamma}(\gamma_{i,0}(\theta),\omega_{i,0}(\theta))$
where $\omega_{i,0}(\theta)$ is a rate parameter and
$\eta_{i}(\theta)|\lambda_{i}^{-2}(\theta)\sim
N(\mu_{i,0}(\theta),1/(\tau_{i,0}(\theta)\lambda_{i}^{-2}(\theta)))$. From the
conjugacy of Model (26), the joint posterior distribution of
$(\eta_{i}(\theta),\lambda_{i}^{-2}(\theta))$ given the simulation outputs,
$Y_{i1}(\theta),Y_{i2}(\theta),\ldots,Y_{iN^{n}_{i}(\theta)}(\theta)$, at
$(i,\theta)$ is
$NG(\mu_{i,n}(\theta),\tau_{i,n}(\theta),\gamma_{i,n}(\theta),\omega_{i,n}(\theta)),$
where

$\displaystyle\mu_{i,n}(\theta)$
$\displaystyle=\frac{\mu_{i,0}(\theta)\tau_{i,0}(\theta)+\sum_{r=1}^{N_{i}^{n}(\theta)}Y_{ir}(\theta)}{\tau_{i,0}(\theta)+N_{i}^{n}(\theta)},\tau_{i,n}(\theta)=\tau_{i,0}(\theta)+N_{i}^{n}(\theta),\gamma_{i,n}(\theta)=\gamma_{i,0}(\theta)+\frac{1}{2}N_{i}^{n}(\theta),$
$\displaystyle\omega_{i,n}(\theta)$
$\displaystyle=\omega_{i,0}(\theta)+\frac{1}{2}\sum_{r=1}^{N_{i}^{n}(\theta)}\left(Y_{ir}(\theta)-m_{i}(\theta)\right)^{2}+\frac{\tau_{i,0}(\theta)N^{n}_{i}(\theta)}{2\left(\tau_{i,0}(\theta)+N^{n}_{i}(\theta)\right)}(m_{i}(\theta)-\mu_{i,0}(\theta))^{2},$

and $m_{i}(\theta)$ is the sample mean of simulation outputs made at
$(i,\theta)$. We introduce
$S^{2}_{i}(\theta):=\mathrm{E}\left[\lambda_{i}^{-2}(\theta)|\mathcal{E}_{n}\right]^{-1}$
as an estimator for the simulation error variance at $(i,\theta)$. Then,
$S^{2}_{i}(\theta)={\omega_{i,n}(\theta)}/{\gamma_{i,n}(\theta)}$ from that
the marginal distribution of $\lambda^{-2}_{i}(\theta)$ is
Gamma$(\gamma_{i,n}(\theta),\omega_{i,n}(\theta))$. Since $S^{-2}_{i}(\theta)$
is a nonnegative martingale, it converges to
$\mathrm{E}[\lambda^{-2}_{i}(\theta)|\mathcal{E}_{\infty}]=\lambda^{-2}_{i}(\theta)$
a.s. from the martingale convergence theorem. Adopting $S_{i}^{2}(\theta_{b})$
as a plug-in estimator of $\lambda^{2}_{i}(\theta_{b})$, the plug-in sample
LDR, $G_{i,n}(\theta_{b})$, is given as

$G_{i,n}(\theta_{b})=\dfrac{(\mu_{i,n}(\theta_{b})-\mu_{i^{b}_{n},n}(\theta_{b}))^{2}}{2\left({S^{2}_{i}(\theta_{b})}/{\alpha_{i,n}(\theta_{b})}+{S^{2}_{i_{n}^{b}}(\theta_{b})}/{\alpha_{i_{n}^{b},n}(\theta_{b})}\right)}.$

If noninformative prior
$(\tau_{i,0}(\theta),\gamma_{i,0}(\theta),\omega_{i,0}(\theta))=(0,0,0)$ is
adopted, then $\mu_{i,n}(\theta)$ and $S_{i}^{2}(\theta_{b})$ coincide with
the sample mean and sample variance of the simulation outputs, respectively.

##  9 Proofs of Proposition 4.1, Theorems 4.2–4.4, and 5.2

###### Proof 9.1

Proof of Proposition 4.1. Let us fix $\theta_{b}.$ Suppose both
$\alpha_{i}(\theta_{b})$ and $\alpha_{i^{b}}(\theta_{b})$ are positive. Then,
we have $\lim_{n\rightarrow\infty}-\frac{1}{n}\log$ (Glynn and Juneja 2004).
If either $\alpha_{i}(\theta_{b})$ or $\alpha_{i^{b}}(\theta_{b})$ is $0$,
then there is no guarantee for to converge to $0$, thus
$\liminf_{n\rightarrow\infty}-\frac{1}{n}\log$ according to the definition of
$G_{i}(\theta_{b})$. Observe that
$\left\\{m_{i,b}^{n}=1\right\\}\subseteq\left\\{\mu_{i,n}(\theta_{b})<\mu_{i^{b},n}(\theta_{b})\right\\}$,
therefore, we have $\liminf_{n\rightarrow\infty}-\frac{1}{n}\log$, as desired.

###### Proof 9.2

Proof of Theorem 4.2. From the definition of $\mathcal{A}_{j}$, for each
$j\neq i^{*}$,

$\displaystyle\max_{\BFM\in\mathcal{A}_{j}}$

Therefore, we have

$\displaystyle\liminf_{n\rightarrow\infty}-\frac{1}{n}\log$ $\displaystyle=$
$\displaystyle\min_{\BFM\in\mathcal{A}_{j}}\liminf_{n\rightarrow\infty}-\frac{1}{n}\log$
$\displaystyle=$
$\displaystyle\min_{\BFM\in\mathcal{A}_{j}}\liminf_{n\rightarrow\infty}-\frac{1}{n}\log$
$\displaystyle=$
$\displaystyle\min_{\BFM\in\mathcal{A}_{j}}\liminf_{n\rightarrow\infty}-\frac{1}{n}\log\prod_{(i,\theta_{b})\in
I(\BFM)}$ $\displaystyle=$
$\displaystyle\min_{\BFM\in\mathcal{A}_{j}}\sum_{(i,\theta_{b})\in
I(\BFM)}\liminf_{n\rightarrow\infty}-\frac{1}{n}\log$ $\displaystyle\geq$
$\displaystyle\min_{\BFM\in\mathcal{A}_{j}}\sum_{(i,\theta_{b})\in
I(\BFM)}G_{i}(\theta_{b}).$

The second equality holds from that only $(i,\theta_{b})\in I(\BFM)$ affects
the LDR of . The third equality holds because simulation outputs are
independent across $\theta_{b}$ given fixed $\bm{\alpha}$. The last inequality
follows from Proposition 4.1.

Discussions on $v_{j}[(i,\theta_{b})]$: Before we proceed to the proofs of
Theorems 4.3, 4.4, and 5.2, we present some discussions on
$v_{j}[(i,\theta_{b})]$. We define
$v_{j}:\Xi\cup\Xi^{\text{adv}}\rightarrow\mathbb{R}$ as a function where
$v_{j}[(i,\theta_{b})]:=d_{j}-d_{j}(\BFM)\text{ such that
}I(\BFM)=\left\\{(i,\theta_{b})\right\\}$. The quantity,
$v_{j}[(i,\theta_{b})]$, measures importance of misspecifying $(i,\theta_{b})$
for the comparison of $i^{*}$ and $j$ at the outer-level problem. For any
$\BFM\in\mathcal{M}$, we have

$\displaystyle d_{j}-d_{j}(\BFM)$ $\displaystyle=$
$\displaystyle\sum_{b=1}^{B}p_{b}\left(m_{j,b}-\mathbf{1}\left\\{j=i^{b}\right\\}-m_{i^{*},b}+\mathbf{1}\\{i^{*}=i^{b}\\}\right)$
(27) $\displaystyle=$
$\displaystyle\sum_{b=1}^{B}\sum_{i=1}^{k}p_{b}\left[\left(m_{j,b}-\mathbf{1}\left\\{j=i^{b}\right\\}-m_{i^{*},b}+\mathbf{1}\\{i^{*}=i^{b}\\}\right)\right]m_{i,b}$
$\displaystyle=$
$\displaystyle\sum_{(i,\theta_{b})\in\Xi\cup\Xi^{\text{adv}}}v_{j}[(i,\theta_{b})]m_{i,b},$

where the first equality holds from the definition of $d_{j}(\BFM)$ in (8) and
the second follows as only one of $m_{1,b},m_{2,b},\ldots,m_{k,b}$ equals one
while all else are zeroes; (27) holds by definition. From (27), Definition (9)
can be rewritten as
$\mathcal{A}_{j}=\left\\{\BFM:d_{j}-\sum\nolimits_{(i,\theta_{b})\in\Xi\cup\Xi^{\text{adv}}}v_{j}[(i,\theta_{b})]m_{i,b}\leq
0\right\\}$ for $j\neq i^{*}$. Recall
$v_{j}[(i,\theta_{b})]^{+}:=\max\\{0,v_{j}[(i,\theta_{b})]\\}$. As
$v_{j}[(i,\theta_{b})]\leq v_{j}[(i,\theta_{b})]^{+}$, we have for any $j\neq
i^{*}$

$\displaystyle\mathcal{A}_{j}\subseteq\widetilde{\mathcal{A}}_{j}:=\bigg{\\{}\BFM:d_{j}-\sum_{(i,\theta_{b})\in\Xi\cup\Xi^{\text{adv}}}v_{j}[(i,\theta_{b})]^{+}m_{i,b}\leq
0\bigg{\\}},$
$\displaystyle\mathrm{LDR}_{j,i^{*}}\underset{\text{\eqref{eq:LDR}}}{=}\;\min_{\BFM\in\mathcal{A}_{j}}\sum_{(i,\theta_{b})\in\Xi\cup\Xi^{\text{adv}}}G_{i}(\theta_{b})m_{i,b}\geq\min_{\BFM\in\widetilde{\mathcal{A}}_{j}}\sum_{(i,\theta_{b})\in\Xi\cup\Xi^{\text{adv}}}G_{i}(\theta_{b})m_{i,b}.$
(28)

Interestingly, it turns out that the inequality in (28) holds in fact as
equality; hence we can replace $\mathcal{A}_{j}$ with
$\widetilde{\mathcal{A}}_{j}$ to obtain the same $\mathrm{LDR}_{j,i^{*}}$.

To see this, we introduce a partial order relation $\preceq$ between two
matrices in $\mathcal{M}$ defined as

$\BFM_{1}\preceq\BFM_{2}\text{ if $I(\BFM_{1})\subseteq I(\BFM_{2})$ holds}.$
(29)

In words, $\preceq$ orders by counting which matrix has more misspecfied
conditional optima. All required conditions (reflexivity, antisymmetry, and
trasitivity) are inherited from the subset relation in (29). Furthermore,
suppose $\BFM_{1}\preceq\BFM_{2}$ with
$\BFM_{1},\BFM_{2}\in\widetilde{\mathcal{A}}_{j}$. Then, we have
$\sum\nolimits_{(i,\theta_{b})\in
I(\BFM_{1})}G_{i}(\theta_{b})\leq\sum\nolimits_{(i,\theta_{b})\in
I(\BFM_{2})}G_{i}(\theta_{b})$. This implies that $\mathrm{LDR}_{j,i^{*}}$
remains unchanged when $\widetilde{\mathcal{A}}_{j}$ is replaced with
$\widetilde{\mathcal{A}}_{j}-\left\\{\BFM_{2}\right\\}$; that is, we can
remove any non-minimal elements with respect to $\preceq$ from
$\widetilde{\mathcal{A}}_{j}$.

The following lemma stipulates that the inequality in (28) indeed holds as
equality.

###### Lemma 9.3

The relation $\min_{\BFM\in\mathcal{A}_{j}}\sum_{(i,\theta_{b})\in
I(\BFM)}G_{i}(\theta_{b})=\min_{\BFM\in\widetilde{\mathcal{A}}_{j}}\sum_{(i,\theta_{b})\in
I(\BFM)}G_{i}(\theta_{b})$ holds.

###### Proof 9.4

Proof. Suppose $\BFM\in\widetilde{\mathcal{A}}_{j}\setminus\mathcal{A}_{j}$.
Then, it suffices to show that $\BFM$ is non-minimal of
$\widetilde{\mathcal{A}}_{j}$. For such $\BFM$, there exists
$\BFM_{1}\in\mathcal{M}$ such that $I(\BFM_{1})=\left\\{(i,\theta_{b})\in
I(\BFM)|v_{j}[(i,\theta_{b})]>0\right\\}$. From the definition of $\preceq$,
we have $\BFM_{1}\preceq\BFM$. Also, $\BFM_{1}\in\mathcal{A}_{j}$ as
$d_{j}(\BFM_{1})=d_{j}-\sum_{(i,\theta_{b})\in\Xi\cup\Xi^{\text{adv}}}v_{j}[(i,\theta_{b})]^{+}m_{i,b}\leq
0$ from (27) and $\BFM\in\mathcal{A}_{j}$. This implies $\BFM$ is non-minimal,
as desired.

Table 2 summarizes classification of each $(i,\theta_{b})$ in
$\Xi\cup\Xi^{\text{adv}}$ according to $v_{j}[(i,\theta_{b})]^{+}$. Since
$v_{j}[(i,\theta_{b})]/p_{b}=\mathbf{1}\left\\{j=i\right\\}-\mathbf{1}\left\\{j=i^{b}\right\\}-\mathbf{1}\left\\{i^{*}=i\right\\}+\mathbf{1}\left\\{i^{*}=i^{b}\right\\}$,
each $v_{j}/p_{b}$ ranges from -2 to 2. Thus, the value in Table 2,
$v_{j}[(i,\theta_{b})]^{+}/p_{b}$, must be 0, 1, or 2. From Table 2, one can
confirm that $v_{j}[(i^{*},\theta_{b})]^{+}=0$ holds for all $j$ and
$(i^{*},\theta_{b})\in\Xi^{\text{adv}}$. Similar to Lemma 9.3, for all $j\neq
i^{*}$, we can rule out $\BFM\in\widetilde{\mathcal{A}}_{j}$ such that
$m_{i^{*},b}=1$ holds for some $(i^{*},\theta_{b})\in\Xi^{\text{adv}}$. In
this regards, we can restrict our attention to $\Xi$ only.

Table 2: Given $j\neq i^{*}$, $\Xi$ is partitioned into five subsets according to the row/column conditions. All $(i,\theta_{b})$ in each subset have an identical value of $v_{j}[(i,\theta_{b})]^{+}/p_{b}$ given in each cell. The first cell is marked with N/A as no pair in $\Xi$ satisfies the two conditions (i.e., null set).  $i$ $\theta_{b}$ |  $i^{b}=i$ |  $i^{b}\neq i$ & $i^{b}\neq i^{*}$ |  $i^{b}=i^{*}$
---|---|---|---
$i=j$ | N/A | $1$ | $2$
$i\neq j$ | $0$ | $0$ | $1$

###### Proof 9.5

Proof of Theorem 4.3.

Recall that we vectorize $(v_{j}[(i,\theta_{b})]^{+})$, $\BFM$, and $\BFG$ by
taking elements in $\Xi$ only, and denote them as $\BFv_{j}$, $\BFm$, and
$\BFg$. From Lemma 9.3, we have the following integer programming (IP) as a
relaxation of (11):

$\min\;\BFg^{\top}\BFm\hskip 15.0pt\text{ subject to }\;\BFv_{j}^{\top}\BFm-
d_{j}\geq 0,\;\BFM\in\mathcal{M}.$ (30)

Recall the definition, $\mathcal{M}=\left\\{\BFM:\BFM\in\\{0,1\\}^{k\times
B},\BFM^{\top}\BFone_{k}=\BFone_{B}\right\\}$. Since the IP in (30) has no
closed-form solution in general, we further drop the integrality and equality
constraints in $\mathcal{M}$ and take a nonnegative constraint only to bound
its optimal objective function value.

$\min\;\BFg^{\top}\BFm\hskip 15.0pt\text{ subject to }\;\BFv_{j}^{\top}\BFm-
d_{j}\geq 0,\;\BFm\geq 0.$ (31)

The resulting linear program in (31) can be easily solved by deriving its KKT
conditions; one can show that the optimal objective function value of (30) is
lower-bounded by
$d_{j}\left\\{{\min_{(i,\theta_{b})\in\Xi}G_{i}(\theta_{b})}/{v_{j}[(i,\theta_{b})]^{+}}\right\\}$,
which in turn implies

$\mathrm{LDR}_{j,i^{*}}\geq
d_{j}\left\\{{\min_{(i,\theta_{b})\in\Xi}G_{i}(\theta_{b})}/{v_{j}[(i,\theta_{b})]^{+}}\right\\}.$
(32)

Let us define

$\widetilde{W}_{i}(\theta_{b}):=\min_{j\neq
i^{*}}\left({d_{j}}/{v_{j}[(i,\theta_{b})]^{+}}\right).$ (33)

From (32), we have

$\min_{j\neq i^{*}}\mathrm{LDR}_{j,i^{*}}\geq\min_{j\neq
i^{*}}\left[d_{j}\left\\{{\min_{(i,\theta_{b})\in\Xi}G_{i}(\theta_{b})}/{v_{j}[(i,\theta_{b})]^{+}}\right\\}\right]\stackrel{{\scriptstyle\eqref{def:Wib}}}{{=}}\min_{(i,\theta_{b})\in\Xi}\widetilde{W}_{i}(\theta_{b})G_{i}(\theta_{b}).$
(34)

Next, we show that $\widetilde{W}_{i}(\theta_{b})$ in (34) can be replaced
with $W_{i}(\theta_{b})$. Without loss of generality, let us consider a
particular $j\neq i^{*}$ and define
$\widetilde{\Xi}_{j}:=\\{(i,\theta_{b})\in\Xi|d_{j}(\BFM)\leq
0,I(\BFM)=\left\\{(i,\theta_{b})\right\\}\\}.$ Namely, misspecifying a single
element of $\widetilde{\Xi}_{j}$ makes Solution $j$ outperform or tie with
$i^{*}$. We further define two disjoint subsets of $\mathcal{A}_{j}$ as
$\mathcal{B}_{1}:=\\{\BFM\in\mathcal{A}_{j}|I(\BFM)\cap\widetilde{\Xi}_{j}=\emptyset\\}$
and
$\mathcal{B}_{2}:=\\{\BFM\in\mathcal{A}_{j}|I(\BFM)\subseteq\widetilde{\Xi}_{j}\\}.$
Suppose there exists $\BFM\in\mathcal{A}_{j}$ such that
$\BFM\notin\mathcal{B}_{1}\cup\mathcal{B}_{2}$ and define
$I(\BFM_{1})=I(\BFM)\cap\widetilde{\Xi}_{j}$. From the definition of $\preceq$
in (29), we have $\BFM_{1}\preceq\BFM$. Hence, any element in
$\mathcal{A}_{j}\setminus\mathcal{B}_{1}\cup\mathcal{B}_{2}$ can be excluded
from the LDR derivation. Since $\mathcal{B}_{1}$ and $\mathcal{B}_{2}$ are
disjoint, we have

$\displaystyle\mathrm{LDR}_{j,i^{*}}$ $\displaystyle=$
$\displaystyle\min_{\BFM\in\mathcal{A}_{j}}\sum_{(i,\theta_{b})\in
I(\BFM)}G_{i}(\theta_{b})=\min_{\BFM\in\mathcal{B}_{1}\cup\mathcal{B}_{2}}\sum_{(i,\theta_{b})\in
I(\BFM)}G_{i}(\theta_{b})$ (35) $\displaystyle=$
$\displaystyle\min\left\\{\min_{\BFM\in\mathcal{B}_{1}}\sum_{(i,\theta_{b})\in
I(\BFM)}G_{i}(\theta_{b}),\min_{\BFM\in\mathcal{B}_{2}}\sum_{(i,\theta_{b})\in
I(\BFM)}G_{i}(\theta_{b})\right\\}$ $\displaystyle=$
$\displaystyle\min\left\\{\min_{\BFM\in\mathcal{B}_{1}}\sum_{(i,\theta_{b})\in
I(\BFM)}G_{i}(\theta_{b}),\min_{(i,\theta_{b})\in\widetilde{\Xi}_{j}}G_{i}(\theta_{b})\right\\}.$

We first consider the first inner minimization problem in (35). Applying the
linear relaxation similar to that in (31), we obtain the following lower bound

$\min_{\BFM\in\mathcal{B}_{1}}\sum_{(i,\theta_{b})\in
I(\BFM)}G_{i}(\theta_{b})\geq
d_{j}\left\\{{\min_{(i,\theta_{b})\in\Xi\setminus\widetilde{\Xi}_{j}}G_{i}(\theta_{b})}/{v_{j}[(i,\theta_{b})]^{+}}\right\\}.$
(36)

Note that $(i,\theta_{b})\in\widetilde{\Xi}_{j}$, if and only if,
$d_{j}/v_{j}[(i,\theta_{b})]^{+}=d_{j}/(d_{j}-d_{j}(\BFM))^{+}\leq 1$, where
$I(\BFM)=\left\\{(i,\theta_{b})\right\\}$. Combining this observation with
(36), (35) can be bounded from below by

$\mathrm{LDR}_{j,i^{*}}\geq\min_{(i,\theta_{b})\in\Xi}\left\\{\max\left[1,\frac{d_{j}}{v_{j}[(i,\theta_{b})]^{+}}\right]G_{i}(\theta_{b})\right\\}.$

Putting all $\mathrm{LDR}_{j,i^{*}}$ terms together, we obtain

$\displaystyle\min_{j\neq i^{*}}\mathrm{LDR}_{j,i^{*}}\geq\min_{j\neq
i^{*}}\min_{(i,\theta_{b})\in\Xi}\left\\{\max\left[1,\frac{d_{j}}{v_{j}[(i,\theta_{b})]^{+}}\right]G_{i}(\theta_{b})\right\\}=\min_{(i,\theta_{b})\in\Xi}\left\\{\max\left[1,\min_{j\neq
i^{*}}\frac{d_{j}}{v_{j}[(i,\theta_{b})]^{+}}\right]G_{i}(\theta_{b})\right\\}.$

The equality holds since $\min_{1\leq n\leq N}\max(1,a_{n})=\max(1,\min_{1\leq
n\leq N}a_{n})$ for all finite $N$.

Finally, we show that $\max\left[1,\min_{j\neq
i^{*}}{d_{j}}/{v_{j}[(i,\theta_{b})]^{+}}\right]=W_{i}(\theta_{b})$ holds for
all $(i,\theta_{b})\in\Xi$. For $\theta_{b}\in\Theta_{i^{*}}$ (or
$i^{b}=i^{*})$, we have $v_{i}[(i,\theta_{b})]^{+}/p_{b}=~{}2$ and
$v_{j}[(i,\theta_{b})]^{+}/p_{b}=1$ for all $j\neq i,i^{*}$ from Table 2.
Then, from Definition (14), we have
$W_{i}(\theta_{b})=\max\left[1,\min\left\\{{d_{i}}/{2},\min_{j\neq
i^{*}}d_{j}\right\\}/p_{b}\right]=\max\left[1,\min_{j\neq
i^{*}}{d_{j}}/{v_{j}[(i,\theta_{b})]^{+}}\right]$. For
$\theta_{b}\notin\Theta_{i^{*}}$ (or $i^{b}\neq i^{*}$), we have
$v_{i}[(i,\theta_{b})]^{+}/p_{b}=1$ and $v_{j}[(i,\theta_{b})]^{+}/p_{b}=0$
for all $j\neq i,i^{*}$ from Table 2. Therefore, again from Definition (14),
$W_{i}(\theta_{b})=\max[1,d_{i}/p_{b}]=\max\left[1,\min_{j\neq
i^{*}}{d_{j}}/{v_{j}[(i,\theta_{b})]^{+}}\right]$, as desired.

In the following, we present proofs of Theorems 4.4 and 5.2. For Theorem 5.2,
it suffices to see how balance weights in (23) are derived since the remaining
parts are straightforward from the proof of Theorem 4.4.

###### Proof 9.6

Proof of Theorem 4.4.

From the definition of $\underline{\mathrm{LDR}}$, (15) is equivalent to

$\displaystyle\max$ $\displaystyle C$ (37) s.t. $\displaystyle C\leq
W_{i}(\theta_{b})G_{i}(\theta_{b}),\forall(i,\theta_{b})\in\Xi,$
$\displaystyle\sum_{i=1}^{k}\sum_{b=1}^{B}\alpha_{i}(\theta_{b})=1,$
$\displaystyle\alpha_{i}(\theta_{b})\geq 0.$

Note that harmonic mean function $f(x,y)=1/(1/x+1/y)$ is concave in both $x$
and $y$, so each $G_{i}(\theta_{b})$ is concave; thus, (37) is a concave
program. On the one hand, $G_{i}(\theta_{b})$ for
$(i,\theta_{b})\in\Xi^{\text{adv}}$ does not appear in this program. For any
such $(i,\theta_{b})$, suppose $\alpha_{i}(\theta_{b})>0$. Then, we can always
strictly improve the objective function value by setting
$\alpha_{i}(\theta_{b})=0$ and increasing
$\alpha_{i^{\prime}}(\theta_{b^{\prime}})$ at
$(i^{\prime},\theta_{b^{\prime}})$ pairs for which the first constraint of
(37) is binding. Hence ${\alpha}_{i^{*}}(\theta_{b})=0$ for all
$\theta_{b}\in\Theta^{c}_{i^{*}}$ at optimality. This proves the last
optimality condition.

To solve (37), we set up a Lagrangian function
$\mathcal{L}(\bm{\rho},\gamma,\bm{\kappa};C,\bm{\alpha})$ defined as

$\begin{split}\mathcal{L}(\bm{\rho},\gamma,\bm{\kappa};C,\bm{\alpha})=-C+&\sum\nolimits_{(i,\theta_{b})\in\Xi}\rho_{i}(\theta_{b})(C-W_{i}(\theta_{b})G_{i}(\theta_{b}))\\\
&+\gamma\left(\sum\nolimits_{(i,\theta_{b})\in(\Xi^{\text{adv}})^{c}}\alpha_{i}(\theta_{b})-1\right)-\sum\nolimits_{(i,\theta_{b})\in(\Xi^{\text{adv}})^{c}}\kappa_{i}(\theta_{b})\alpha_{i}(\theta_{b}).\end{split}$
(38)

where $\bm{\rho}=\left\\{\rho_{i}(\theta_{b})\right\\}$, $\gamma$, and
$\bm{\kappa}=\left\\{\kappa_{i}(\theta_{b})\right\\}$ are the Lagrangian
multipliers of the first, second, and third constraints in (37), respectively.
Note that in the third and fourth terms summations are over
$(i,\theta_{b})\in(\Xi^{\text{adv}})^{c}$ thanks to the last optimality
condition. Next, we argue that $\kappa_{i}(\theta_{b})=0$ for all
$(i,\theta_{b})\in(\Xi^{\text{adv}})^{c}$. Suppose $\alpha_{i}(\theta_{b})=0$
for some $(i,\theta_{b})\in\left(\Xi^{\text{adv}}\right)^{c}$. Then,
$G_{i}(\theta_{b})=0$ or $G_{j}(\theta_{b})$ = 0 for some $j$ (when
$i=i^{b}$). This implies $C\leq 0$. However, if we set
$\alpha_{i}(\theta_{b})=1/kB$ for all $1\leq i\leq k,1\leq b\leq B$, the
optimal $C$ must be positive; for all $(i,\theta_{b})\in\Xi$,
$W_{i}(\theta_{b})>0$ and positive $\alpha_{i}(\theta_{b})$’s make
$G_{i}(\theta_{b})$ positive. Therefore, optimal $\alpha_{i}(\theta_{b})$ must
be positive for all $(i,\theta_{b})\in\left(\Xi^{\text{adv}}\right)^{c}$. As a
result, Lagrangian multipliers associated with the nonnegativity constraint
should be zero as they are not binding, i.e., $\kappa_{i}(\theta_{b})=0$.
Further, the uniform allocation, $\alpha_{i}(\theta_{b})=1/kB$, with $C=0$
becomes an interior point, which is strictly feasible. Hence, Slater’s
condition implies that the strong duality holds; thus, Karush-Kuhn-Tucker
(KKT) conditions become necessary and sufficient conditions for optimality.
Consequently, we derive the KKT conditions for the following Lagrangian
function instead of (38)

$\mathcal{L}(\bm{\rho},\gamma;C,\bm{\alpha}):=-C+\sum\nolimits_{(i,\theta_{b})\in\Xi}\rho_{i}(\theta_{b})\left(C-W_{i}(\theta_{b})G_{i}(\theta_{b})\right)+\gamma\left(\sum\nolimits_{(i,\theta_{b})\in\left(\Xi^{\text{adv}}\right)^{c}}\alpha_{i}(\theta_{b})-1\right)$
(39)

as follows:

$\displaystyle\rho_{i}(\theta_{b})$ $\displaystyle\geq$ $\displaystyle
0,\forall(i,\theta_{b})\in\Xi,$ (40)
$\displaystyle\frac{\partial\mathcal{L}}{\partial C}$ $\displaystyle=$
$\displaystyle-1+\sum_{(i,\theta_{b})\in\Xi}\rho_{i}(\theta_{b})=0,$ (41)
$\displaystyle\frac{\partial\mathcal{L}}{\partial\alpha_{i}(\theta_{b})}$
$\displaystyle=$
$\displaystyle-\rho_{i}(\theta_{b})W_{i}(\theta_{b})\frac{\partial
G_{i}(\theta_{b})}{\partial\alpha_{i}(\theta_{b})}+\gamma=0,\forall(i,\theta_{b})\in\Xi,$
(42)
$\displaystyle\frac{\partial\mathcal{L}}{\partial\alpha_{i^{b}}(\theta_{b})}$
$\displaystyle=$ $\displaystyle-\sum_{i\neq
i^{b},i^{*}}\rho_{i}(\theta_{b})W_{i}(\theta_{b})\frac{\partial
G_{i}(\theta_{b})}{\partial\alpha_{i^{b}}(\theta_{b})}+\gamma=0,\forall b,$
(43)
$\displaystyle\rho_{i}(\theta_{b})\left(C-W_{i}(\theta_{b})G_{i}(\theta_{b})\right)$
$\displaystyle=$ $\displaystyle 0,\forall(i,\theta_{b})\in\Xi.$ (44)

From (40) and (41), there exists $(i,\theta_{b})\in\Xi$ such that
$\rho_{i}(\theta_{b})>0$. Therefore, we further obtain $\gamma~{}>~{}0$ from
(42), which in turn implies that $\rho_{i}(\theta_{b})>0$ for all
$(i,\theta_{b})\in\Xi$. Finally, (44) implies that
$C=W_{i}(\theta_{b})G_{i}(\theta_{b})$ for all $(i,\theta_{b})\in\Xi$, which
corresponds to the pairwise balance condition. Combining (42) with (43), we
have

$\displaystyle 1$ $\displaystyle=\frac{1}{\gamma}\sum_{i\neq
i^{b},i^{*}}\rho_{i}(\theta_{b})W_{i}(\theta_{b})\frac{\partial
G_{i}(\theta_{b})}{\partial\alpha_{i^{b}}(\theta_{b})}=\sum_{i\neq
i^{b},i^{*}}\frac{\rho_{i}(\theta_{b})W_{i}(\theta_{b})\frac{\partial
G_{i}(\theta_{b})}{\partial\alpha_{i^{b}}(\theta_{b})}}{\rho_{i}(\theta_{b})W_{i}(\theta_{b})\frac{\partial
G_{i}(\theta_{b})}{\partial\alpha_{i}(\theta_{b})}}$ (45)
$\displaystyle=\sum_{i\neq i^{b},i^{*}}\left(\frac{\partial
G_{i}(\theta_{b})}{\partial\alpha_{i^{b}}(\theta_{b})}\middle/\frac{\partial
G_{i}(\theta_{b})}{\partial\alpha_{i}(\theta_{b})}\right)=\sum_{i\neq
i^{b},i^{*}}\frac{\lambda^{2}_{i^{b}}(\theta_{b})/\alpha^{2}_{i^{b}}(\theta_{b})}{\lambda^{2}_{i}(\theta_{b})/\alpha^{2}_{i}(\theta_{b})}.$

For the third equality, we exploit that $\dfrac{\partial
G_{i}(\theta_{b})}{\partial\alpha_{i}(\theta_{b})}=G_{i}(\theta_{b})\dfrac{\lambda^{2}_{i}(\theta_{b})/\alpha_{i}^{2}(\theta_{b})}{\lambda^{2}_{i}(\theta_{b})/\alpha_{i}^{2}(\theta_{b})+\lambda^{2}_{i^{b}}(\theta_{b})/\alpha_{i^{b}}^{2}(\theta_{b})}.$

###### Proof 9.7

Proof of Theorem 5.2. Consider (21) first. The formulation (21) is equivalent
to

$\displaystyle\max$ $\displaystyle C$ s.t. $\displaystyle C\leq
W_{i}(\theta_{b})G_{i}(\theta_{b}),\forall(i,\theta_{b})\in\Xi,$
$\displaystyle C\leq G_{i}(\theta_{b}),\forall i\neq
i^{*},\theta_{b}\in\Theta_{i^{*}},$ $\displaystyle C\leq
G_{i^{*}}(\theta_{b}),\forall\theta_{b}\in\Theta_{i^{*}}^{c},$
$\displaystyle\sum_{i=1}^{k}\sum_{b=1}^{B}\alpha_{i}(\theta_{b})=1,$
$\displaystyle\alpha_{i}(\theta_{b})\geq 0.$

As in (37), we rewrite the first two constraints to obtain

$\displaystyle\max$ $\displaystyle C$ (46) s.t. $\displaystyle
C\leq\min\\{W_{i}(\theta_{b}),1\\}G_{i}(\theta_{b}),\forall(i,\theta_{b})\in\Xi\;\;\textup{where}\;\;\theta_{b}\in\Theta_{i^{*}}$
$\displaystyle C\leq
W_{i}(\theta_{b})G_{i}(\theta_{b}),\forall(i,\theta_{b})\in\Xi\;\;\textup{where}\;\;\theta_{b}\in\Theta^{c}_{i^{*}}$
$\displaystyle C\leq
G_{i^{*}}(\theta_{b}),\forall\theta_{b}\in\Theta_{i^{*}}^{c},$
$\displaystyle\sum_{i=1}^{k}\sum_{b=1}^{B}\alpha_{i}(\theta_{b})=1,$
$\displaystyle\alpha_{i}(\theta_{b})\geq 0.$

Unlike Theorem 4.4, all $G_{i^{*}}(\theta_{b})$ with
$\theta_{b}\in\Theta^{c}_{i^{*}}$ appear in (46) from the third constraint in
(46). Therefore, the sampling ratios
$\left\\{\alpha_{i^{*}}(\theta_{b}),\theta_{b}\in\Theta^{c}_{i^{*}}\right\\}$
are positive in this case. From the definition of
$\left\\{W^{\textup{ACC}}_{i}(\theta_{b})\right\\}$ in (23), (46) is
equivalent to

$\displaystyle\max$ $\displaystyle C$ (47) s.t. $\displaystyle C\leq
W^{\textup{ACC}}_{i}(\theta_{b})G_{i}(\theta_{b}),\forall(i,\theta_{b})\in\Xi\cup\Xi^{\text{adv}},$
$\displaystyle\sum_{i=1}^{k}\sum_{b=1}^{B}\alpha_{i}(\theta_{b})=1,$
$\displaystyle\alpha_{i}(\theta_{b})\geq 0.$

Again, from the KKT conditions, we can derive the pairwise and global balance
conditions; note that the only difference here from the proof of Theorem 4.4
is the Lagrangian function for (47), which has the same from as (39) except
that $\left\\{W_{i}(\theta_{b})\right\\}$ are substituted with
$\left\\{W^{\textup{ACC}}_{i}(\theta_{b})\right\\}$. For the FN case in (22),
the proof proceeds similarly.

##  10 Proofs of Theorems 5.1 and 5.3

We reverse the order to show Theorem 5.3 first. Showing Theorem 5.1 is more
involved as Algorithm 2 combines both posterior mean plug in and posterior
sampling to decide which solution-parameter pair to simulate next.

###### Proof 10.1

Proof of Theorem 5.3.

For two positive random sequences $\left\\{a_{n}\right\\}_{n\in\mathbb{N}}$
and $\left\\{b_{n}\right\\}_{n\in\mathbb{N}}$, we adopt the notation,
$a_{n}={O_{\text{a.s.}}}(b_{n})$, if
$\limsup_{n\rightarrow\infty}{a_{n}}/{b_{n}}<\infty$ a.s. We further define
$a_{n}={\Theta_{\text{a.s.}}}(b_{n})$, if $a_{n}={O_{\text{a.s.}}}(b_{n})$ and
$b_{n}={O_{\text{a.s.}}}(a_{n})$. Equivalently,
$a_{n}={\Theta_{\text{a.s.}}}(b_{n})$ holds, if and only if,
$0<\liminf_{n\rightarrow\infty}{a_{n}}/{b_{n}}\leq\limsup_{n\rightarrow\infty}{a_{n}}/{b_{n}}<\infty$
a.s. Do not confuse ${\Theta_{\text{a.s.}}}$ with the set notation defined for
$\theta$.

For each sample path, we define two partitions of the set of positive integers
$\mathbb{N}$, $\left\\{\mathcal{N}_{i}(\theta_{b})\right\\}_{1\leq i\leq
k,1\leq b\leq B}$ and $\left\\{\mathcal{N}(\theta_{b})\right\\}_{1\leq b\leq
B}$, where

$\mathcal{N}_{i}(\theta_{b}):=\\{n\in\mathbb{N}|\text{$(i,\theta_{b})$ is
sampled at the $(n+1)$th simulation run}\\}\text{ and
}\mathcal{N}(\theta_{b}):=\bigcup_{i=1}^{k}\mathcal{N}_{i}(\theta_{b}).$

Equivalently, $\mathcal{N}(\theta_{b})=\\{n\in\mathbb{N}|\min\nolimits_{1\leq
i\leq k,i\neq
i_{n}^{b}}W_{i,n}(\theta_{b})G_{i,n}(\theta_{b})<\min\nolimits_{1\leq i\leq
k,i\neq
i_{n}^{b^{\prime}}}W_{i,n}(\theta_{b^{\prime}})G_{i,n}(\theta_{{b^{\prime}}}),\forall{b^{\prime}}\neq
b\\}$ from the sampling rule. We also have
$\bigcup_{i=1}^{k}\bigcup_{b=1}^{B}\mathcal{N}_{i}(\theta_{b})=\bigcup_{b=1}^{B}\mathcal{N}(\theta_{b})=\mathbb{N}$.
Further, the total allocation to $\theta_{b}$ by the $n$th replication is
denoted by
$N_{\text{tot}}^{n}(\theta_{b}):=\sum_{i=1}^{k}N_{i}^{n}(\theta_{b})$.

For each $\theta_{b}$, define
$S_{1}(\theta_{b})=\\{i:\lim_{n\rightarrow\infty}N^{n}_{i}(\theta_{b})=\infty\\}$
and
$S_{2}(\theta_{b})=\\{i:\lim_{n\rightarrow\infty}N^{n}_{i}(\theta_{b})<\infty\\}$;
together, they form a partition of the solution set. Observe that both
$\argmin_{i\in S_{1}(\theta_{b})}\mu_{i,n}(\theta_{b})$ and $\argmin_{i\in
S_{2}(\theta_{b})}\mu_{i,n}(\theta_{b})$ converge given the sample path; the
former is by the SLLN and latter is due to the finiteness of
$N^{n}_{i}(\theta_{b})$ for all $i\in S_{2}(\theta_{b})$. Hence,
$i_{n}^{b}=i^{b}_{0}$ holds for some $i^{b}_{0}$ for all sufficiently large
$n$. Corresponding to $\\{i_{0}^{b}\\}$, we define
$\Xi_{0},\Xi_{0}^{\text{adv}}$ and $\\{W_{i,0}(\theta_{b})\\}$ so that for
sufficiently large $n$, $\Xi_{n}=\Xi_{0}$,
$\Xi_{n}^{\text{adv}}=\Xi_{0}^{\text{adv}}$, and
$W_{i,n}(\theta_{b})=W_{i,0}(\theta_{b})$. Note that we cannot guarantee
$i^{b}_{0}=i^{b}$ since $i^{b}_{0}$ depends on sample paths assigned to
$(i,\theta_{b})$ with $i\in S_{2}(\theta_{b})$ as long as
$S_{2}(\theta_{b})\neq\emptyset$. As a first step, we show that for all
$\theta_{b}$, $i_{0}^{b}=i^{b}$ holds almost surely.

Step 1. $N^{n}_{i}(\theta_{b})\rightarrow\infty$ a.s. as $n\rightarrow\infty$
for all $1\leq i\leq k,1\leq b\leq B$.

Since $\sum_{b=1}^{B}N^{n}_{\text{tot}}(\theta_{b})=n$, there exists at least
one $\theta_{b}$ such that $N^{n}_{\text{tot}}(\theta_{b})\rightarrow\infty$
as $n\rightarrow\infty$. For such $\theta_{b}$, suppose there exists $1\leq
j\leq k$ such that $\lim_{n\rightarrow\infty}N^{n}_{j}(\theta_{b})<\infty$. We
first argue that $N_{i_{0}^{b}}^{n}(\theta_{b})\rightarrow\infty$. By means of
contradiction, suppose otherwise. Because
$\lim_{n\rightarrow\infty}N^{n}_{i_{0}^{b}}(\theta_{b})<\infty$,

$\left(\dfrac{N^{n}_{i_{0}^{b}}(\theta_{b})}{\lambda_{i_{0}^{b}}(\theta_{b})}\right)^{2}<\sum\nolimits_{i\neq
i_{0}^{b}}\left(\dfrac{N^{n}_{i}(\theta_{b})}{\lambda_{i}(\theta_{b})}\right)^{2}$
(48)

must hold at finitely many $n\in\mathcal{N}(\theta_{b})$ because otherwise the
algorithm would allocate infinite simulation effort to $i_{0}^{b}$. This
implies that $\lim_{n\rightarrow\infty}N^{n}_{i}(\theta_{b})<\infty$ for all
$1\leq i\leq k$, which contradicts that
$N^{n}_{\text{tot}}(\theta_{b})\rightarrow\infty$. Consequently, we have
$j\neq i_{0}^{b}$, and $N^{n}_{i_{0}^{b}}(\theta_{b})\rightarrow\infty$
follows.

To achieve $N^{n}_{i_{0}^{b}}(\theta_{b})\rightarrow\infty$, (48) should hold
infinitely often. Hence, there exists $i_{1}\neq i_{0}$ such that
$N^{n}_{i_{1}}(\theta_{b})\rightarrow\infty$ that makes the right-hand side of
(48) diverge. For such $i_{1}$, we can see

$nG_{i_{1},n}(\theta_{b})=\dfrac{(\mu_{i_{1},n}(\theta_{b})-\mu_{i_{0}^{b},n}(\theta_{b}))^{2}}{2\left(\dfrac{\lambda^{2}_{i_{1}}(\theta_{b})}{N_{i_{1}}^{n}(\theta_{b})}+\dfrac{\lambda^{2}_{i_{0}^{b}}(\theta_{b})}{N_{i_{0}^{b}}^{n}(\theta_{b})}\right)}\rightarrow\infty.$
(49)

Observe that we have $\lim_{n\rightarrow\infty}nG_{j,n}(\theta_{b})<\infty$ as
$\lim_{n\rightarrow\infty}N_{j}^{n}(\theta_{b})<\infty$, which in turn implies
that
$W_{j,0}(\theta_{b})G_{j,n}(\theta_{b})<W_{i_{1},0}(\theta_{b})G_{i_{1},n}(\theta_{b})$
holds for all sufficiently large $n$. Then, Steps 4 and 5 of Algorithm 1 would
never select $(i_{1},\theta_{b})$ as the next sampling pair, which contradicts
$N_{i_{1}}^{n}(\theta_{b})\rightarrow\infty$. Hence, we conclude that
$N^{n}_{i}(\theta_{b})\rightarrow\infty$ for all $1\leq i\leq k$.

Now, it remains to show that $\\{N_{i}(\theta_{b})\\}_{1\leq i\leq k}$ go to
infinity for all $\theta_{b}$. By Step 4 of Algorithm 1 again, we claim that
$nG_{i^{\prime},n}(\theta_{b^{\prime}})\rightarrow\infty$ for all
$(i^{\prime},\theta_{{b^{\prime}}})\in\Xi_{0}\cup\Xi_{0}^{\text{adv}}$ from
the existence of $\theta_{b}$ such that
$N_{\text{tot}}^{n}(\theta_{b})\rightarrow\infty$; otherwise, the existence of
$(i^{\prime},\theta_{b^{\prime}})$ satisfying
$\limsup_{n\rightarrow\infty}nW_{i^{\prime},0}(\theta_{b^{\prime}})G_{i^{\prime},n}(\theta_{b^{\prime}})<~{}\infty$
implies
$W_{i^{\prime},0}(\theta_{b^{\prime}})G_{i^{\prime},n}(\theta_{b^{\prime}})<\min_{i\neq
i_{0}^{b}}W_{i,0}(\theta_{b})G_{i,n}(\theta_{b})$ for all sufficiently large
$n$, and this makes $n\in\mathcal{N}(\theta_{b})$ hold for only finitely many
$n$, which contradicts $N_{\text{tot}}^{n}(\theta_{b})\rightarrow\infty$.
Observe that in order to have
$nG_{i^{\prime},n}(\theta_{b^{\prime}})\rightarrow\infty$ a.s.,
$N_{i^{\prime}}^{n}(\theta_{b^{\prime}})$ must diverge as seen in (49). ∎
(Step 1)

Since $N^{n}_{i}(\theta_{b^{\prime}})\rightarrow\infty$ for any
$(i,\theta_{b})$, the SLLN implies that we can characterize all conditional
optima for sufficiently large $n$. That is, we always have $i^{b}_{n}=i^{b}$,
$i^{*}_{n}=i^{*}$, and $W_{i,n}(\theta_{b})=W_{i}(\theta_{b})$ for
sufficiently large $n$. Thus, we drop $n$ from these quantities in the
remainder of the proof.

Note that for any sequence $\left\\{a_{n}\right\\}_{n\in\mathbb{N}}$
satisfying $\limsup_{n\rightarrow\infty}a_{n}=a$ with
$a\in\mathbb{R}\cup\left\\{\infty\right\\}$, there exists subsequence
$\left\\{m_{\ell}\right\\}_{\ell\in\mathbb{N}}$ such that
$\lim_{\ell\rightarrow\infty}a_{m_{\ell}}=a$. In the following, we repeatedly
make this _convergent subsequence argument_ to verify the limit results.

Step 2. For each $b$,
$\limsup_{n\rightarrow\infty}{N^{n}_{i}(\theta_{b})}/{N^{n}_{j}(\theta_{b})}<\infty$
holds for all $i\neq j$.

First, let us consider when $j=i^{b}$ and $i\neq i^{b}$. Suppose
$\limsup_{n\rightarrow\infty}{N^{n}_{i}(\theta_{b})}/{N^{n}_{i^{b}}(\theta_{b})}=\infty$
by means of contradiction. Then, there exists subsequence
$\left\\{n_{\ell}\right\\}_{\ell\in\mathbb{N}}$ such that
$\lim_{\ell\rightarrow\infty}{N^{n_{\ell}}_{i}(\theta_{b})}/{N^{n_{\ell}}_{i^{b}}(\theta_{b})}=\infty$.
We further define another subsequence
$\left\\{m_{\ell}\right\\}_{\ell\in\mathbb{N}}$ such that
$m_{\ell}:=\sup\left\\{m\leq
n_{\ell}:m\in\mathcal{N}_{i}(\theta_{b})\right\\}$. Note that
$\left\\{m_{\ell}\right\\}_{\ell\in\mathbb{N}}$ is nondecreasing and
$m_{\ell}\rightarrow\infty$ since $N^{n}_{i}(\theta_{b})\rightarrow\infty$. We
also have
$\lim_{\ell\rightarrow\infty}{N^{m_{\ell}}_{i}(\theta_{b})}/{N^{m_{\ell}}_{i^{b}}(\theta_{b})}\geq\lim_{\ell\rightarrow\infty}{N^{n_{\ell}}_{i}(\theta_{b})}/{N^{n_{\ell}}_{i^{b}}(\theta_{b})}=\infty$
since $N_{i^{b}}^{n_{\ell}}(\theta_{b})\geq N_{i^{b}}^{m_{\ell}}(\theta_{b})$
and $N_{i}^{n_{\ell}}(\theta_{b})=N_{i}^{m_{\ell}}(\theta_{b})$ from the
definition of $m_{\ell}$. Thus, for large $\ell$, we obtain that

$\left(\dfrac{N_{i^{b}}^{m_{\ell}}(\theta_{b})}{\lambda_{i^{b}}(\theta_{b})}\right)^{2}<\left(\dfrac{N_{i}^{m_{\ell}}(\theta_{b})}{\lambda_{i}(\theta_{b})}\right)^{2}<\sum_{j\neq
i^{b}}\left(\dfrac{N^{m_{\ell}}_{j}(\theta_{b})}{\lambda_{j}(\theta_{b})}\right)^{2},$

where the last inequality holds since $i\neq i^{b}$. However, Step 4 in
Algorithm 1 forces $(i^{b},\theta_{b})$ to be sampled next should $\theta_{b}$
is sampled at all. This contradicts that
$m_{\ell}\in\mathcal{N}_{i}(\theta_{b})$. Therefore, we have
$\limsup_{n\rightarrow\infty}{N^{n}_{i}(\theta_{b})}/{N^{n}_{i^{b}}(\theta_{b})}<\infty$.

Let us consider when $i,j\neq i^{b}$. Again, suppose
$\limsup_{n\rightarrow\infty}{N^{n}_{i}(\theta_{b})}/{N^{n}_{j}(\theta_{b})}=\infty$
by means of contradiction. Similarly, we can take subsequence
$\left\\{m_{\ell}\right\\}_{\ell\in\mathbb{N}}\subseteq\mathcal{N}_{i}(\theta_{b})$
such that
$\lim_{\ell\rightarrow\infty}{N^{m_{\ell}}_{i}(\theta_{b})}/{N^{m_{\ell}}_{j}(\theta_{b})}=\infty$.
Furthermore, we have

$\displaystyle nG_{i,m_{\ell}}(\theta_{b})$
$\displaystyle=\dfrac{(\mu_{i,m_{\ell}}(\theta_{b})-\mu_{i^{b},m_{\ell}}(\theta_{b}))^{2}}{2\left(\lambda^{2}_{i}(\theta_{b})+\lambda^{2}_{i^{b}}(\theta_{b})N^{m_{\ell}}_{i}(\theta_{b})/N^{m_{\ell}}_{i^{b}}(\theta_{b})\right)}N^{m_{\ell}}_{i}(\theta_{b}),$
(50) $\displaystyle nG_{j,m_{\ell}}(\theta_{b})$
$\displaystyle\leq\dfrac{(\mu_{j,m_{\ell}}(\theta_{b})-\mu_{i^{b},m_{\ell}}(\theta_{b}))^{2}}{2\lambda^{2}_{j}(\theta_{b})}N^{m_{\ell}}_{j}(\theta_{b}),$

where the inequality follows because
$\lambda^{2}_{j}(\theta_{b})/N^{m_{\ell}}_{j}(\theta_{b})+\lambda^{2}_{i^{b}}(\theta_{b})/N^{m_{\ell}}_{i^{b}}(\theta_{b})\geq\lambda^{2}_{j}(\theta_{b})/N^{m_{\ell}}_{j}(\theta_{b})$
holds. Combining (50) with the product rule of limit superior, we can derive
that

$\displaystyle\limsup_{\ell\rightarrow\infty}\dfrac{G_{j,m_{\ell}}(\theta_{b})}{G_{i,m_{\ell}}(\theta_{b})}$
$\displaystyle\leq$
$\displaystyle\frac{\limsup_{\ell\rightarrow\infty}\dfrac{(\mu_{j,m_{\ell}}(\theta_{b})-\mu_{i^{b},m_{\ell}}(\theta_{b}))^{2}}{2\lambda^{2}_{j}(\theta_{b})}}{\liminf_{\ell\rightarrow\infty}\dfrac{(\mu_{i,m_{\ell}}(\theta_{b})-\mu_{i^{b},m_{\ell}}(\theta_{b}))^{2}}{2\left(\lambda^{2}_{i}(\theta_{b})+\lambda^{2}_{i^{b}}(\theta_{b})N^{m_{\ell}}_{i}(\theta_{b})/N^{m_{\ell}}_{i^{b}}(\theta_{b})\right)}}\limsup_{\ell\rightarrow\infty}\dfrac{N^{m_{\ell}}_{j}(\theta_{b})}{N^{m_{\ell}}_{i}(\theta_{b})}$
$\displaystyle=$ $\displaystyle
C_{1}\limsup_{\ell\rightarrow\infty}\dfrac{N^{m_{\ell}}_{j}(\theta_{b})}{N^{m_{\ell}}_{i}(\theta_{b})}=0,$

where $C_{1}$ is a finite constant obtained by computing the limit. Note that
$C_{1}$ is positive; the liminf term in (10.1) is positive since
$\limsup_{\ell\rightarrow\infty}N^{m_{\ell}}_{i}(\theta_{b})/N^{m_{\ell}}_{i^{b}}(\theta_{b})<\infty$
by the first case and the other terms are also positive and finite by the
SLLN. Hence, we get
$\lim_{\ell\rightarrow\infty}{G_{j,m_{\ell}}(\theta_{b})}/{G_{i,m_{\ell}}(\theta_{b})}=0$.
Then, for sufficiently large $\ell$,
$W_{j}(\theta_{b})G_{j,m_{\ell}}(\theta_{b})<W_{i}(\theta_{b})G_{i,m_{\ell}}(\theta_{b})$
holds and it implies that $m_{\ell}\notin\mathcal{N}_{i}(\theta_{b})$, which
contradicts the definition of $\left\\{m_{\ell}\right\\}_{\ell\in\mathbb{N}}$.
Thus, we obtain
$\limsup_{n\rightarrow\infty}{N^{n}_{i}(\theta_{b})}/{N^{n}_{j}(\theta_{b})}<\infty$
for all $i,j\neq i^{b}$.

Finally, it only remains to show the case when $i=i^{b}$ and $j\neq i^{b}$. We
argue again by contradiction; suppose
$\limsup_{n\rightarrow\infty}N^{n}_{i^{b}}(\theta_{b})/{N^{n}_{j}(\theta_{b})}=\infty$.
Proceeding similarly, we can take subsequence
$\left\\{m_{\ell}\right\\}_{\ell\in\mathbb{N}}\subseteq\mathcal{N}_{i^{b}}(\theta_{b})$
such that
$\lim_{\ell\rightarrow\infty}N^{m_{\ell}}_{i^{b}}(\theta_{b})/N^{m_{\ell}}_{j}(\theta_{b})=\infty$.
For any $i^{\prime},j^{\prime}\neq i^{b}$, we have shown that
$\limsup_{n\rightarrow\infty}N^{n}_{i^{\prime}}(\theta_{b})/N^{n}_{j^{\prime}}(\theta_{b})<\infty$.
From this, we obtain
$\limsup_{\ell\rightarrow\infty}{N^{m_{\ell}}_{i^{b}}(\theta_{b})}/{N^{m_{\ell}}_{j^{\prime}}(\theta_{b})}=\infty$
for all $j^{\prime}\neq i^{b}$. Repeating the same logic, we can take
subsequence
$\left\\{m^{\prime}_{\ell}\right\\}_{\ell\in\mathbb{N}}\subseteq\mathcal{N}_{i^{b}}(\theta_{b})$
satisfying
$\lim_{\ell\rightarrow\infty}N^{m^{\prime}_{\ell}}_{i^{b}}(\theta_{b})/N^{m^{\prime}_{\ell}}_{j^{\prime}}(\theta_{b})=\infty$
for any $j^{\prime}\neq i^{b}$. Then, for sufficiently large $\ell$, the
global balance condition,
$\left({N_{i^{b}}^{m^{\prime}_{\ell}}(\theta_{b})}/{\lambda_{i^{b}}(\theta_{b})}\right)^{2}>\sum_{j^{\prime}\neq
i^{b}}\left({N^{m^{\prime}_{\ell}}_{j^{\prime}}(\theta_{b})}/{\lambda_{j^{\prime}}(\theta_{b})}\right)^{2}$,
always holds. Then, from Step 4 of Algorithm 1,
$m^{\prime}_{\ell}\notin\mathcal{N}_{i^{b}}(\theta_{b})$, which contradicts
the definition of $m^{\prime}_{\ell}$. ∎ (Step 2)

Step 3.
$\limsup_{n\rightarrow\infty}{N^{n}_{i}(\theta_{b})}/{N^{n}_{j}(\theta_{c})}<\infty$
holds for all $(i,\theta_{b})\neq(j,\theta_{c})$.

From Step 2, we have
$N^{n}_{i}(\theta_{b})={\Theta_{\text{a.s.}}}(N^{n}_{\text{tot}}(\theta_{b}))$.
Hence, it suffices to show
$\limsup_{n\rightarrow\infty}{N^{n}_{\text{tot}}(\theta_{b})}/{N^{n}_{\text{tot}}(\theta_{c})}<\infty$.
Suppose otherwise. Then, we have subsequence
$\left\\{m_{\ell}\right\\}_{\ell\in\mathbb{N}}\subseteq\mathcal{N}(\theta_{b})$
such that
$\lim_{\ell\rightarrow\infty}N^{m_{\ell}}_{\text{tot}}(\theta_{b})/N^{m_{\ell}}_{\text{tot}}(\theta_{c})=\infty$.
Step 2 implies that
$m_{\ell}G_{i,m_{\ell}}(\theta_{b})={\Theta_{\text{a.s.}}}(N^{m_{\ell}}_{\text{tot}}(\theta_{b}))$
for all $1\leq i\leq k,1\leq b\leq B$. As a byproduct, we derive

$\frac{\min_{1\leq i\leq k,i\neq i^{c}}G_{i,m_{\ell}}(\theta_{c})}{\min_{1\leq
i\leq k,i\neq
i^{b}}G_{i,m_{\ell}}(\theta_{b})}=\frac{{\Theta_{\text{a.s.}}}(N^{m_{\ell}}_{\text{tot}}(\theta_{c}))}{{\Theta_{\text{a.s.}}}(N^{m_{\ell}}_{\text{tot}}(\theta_{b}))}={\Theta_{\text{a.s.}}}\left(\frac{N^{m_{\ell}}_{\text{tot}}(\theta_{c})}{N^{m_{\ell}}_{\text{tot}}(\theta_{b})}\right)\rightarrow
0.$ (52)

Hence, for sufficiently large $\ell$, (52) implies that $\min_{1\leq i\leq
k,i\neq i^{c}}W_{i}(\theta_{c})G_{i,m_{\ell}}(\theta_{c})<\min_{1\leq i\leq
k,i\neq i^{b}}W_{i}(\theta_{b})G_{i,m_{\ell}}(\theta_{b})$, that is,
$m_{\ell}\notin\mathcal{N}(\theta_{b})$, which contradicts the definition of
$m_{\ell}$. ∎ (Step 3)

From Step 3, we have $N_{i}^{n}(\theta_{b})={\Theta_{\text{a.s.}}}(n)$ for all
$(i,\theta_{b})$ and
$N^{n}_{\text{tot}}(\theta_{b})={\Theta_{\text{a.s.}}}(n)$ for all
$\theta_{b}$.

Step 4. The global balance condition holds in the limit.

In this step, we aim to show that for each $1\leq b\leq B$,

$\lim_{n\rightarrow\infty}\Delta_{n}(\theta_{b})=0,$ (53)

where auxiliary variable $\Delta_{n}(\theta_{b})$ is defined as

$\displaystyle\Delta_{n}(\theta_{b}):=\left(\frac{n}{N^{n}_{\text{tot}}(\theta_{b})}\right)^{2}\left\\{\left(\dfrac{\alpha_{i^{b},n}(\theta_{b})}{\lambda_{i^{b}}(\theta_{b})}\right)^{2}-\sum_{\ell\neq
i^{b}}\left(\dfrac{\alpha_{\ell,n}(\theta_{b})}{\lambda_{\ell}(\theta_{b})}\right)^{2}\right\\}.$

Combining with $N^{n}_{\text{tot}}(\theta_{b})={\Theta_{\text{a.s.}}}(n)$,
(53) implies that

$\lim\nolimits_{n\rightarrow\infty}\left\\{\left({\alpha_{i^{b},n}(\theta_{b})}/{\lambda_{i^{b}}(\theta_{b})}\right)^{2}-\sum_{\ell\neq
i^{b}}\left({\alpha_{\ell,n}(\theta_{b})}/{\lambda_{\ell}(\theta_{b})}\right)^{2}\right\\}=~{}0.$

To show (53), we extend Lemma 4.1 in Chen and Ryzhov (2019b) in Lemma 10.2
below. The proof of Lemma 10.2 can be found in Section 11.

######  Lemma 10.2 (Extension of Lemma 4.1 in Chen and Ryzhov (2019b))

We have

$\Delta_{n+1}(\theta_{b})-\Delta_{n}(\theta_{b})\begin{cases}>0,&\mbox{ if
$n\in\mathcal{N}_{i^{b}}(\theta_{b})$,}\\\ <0,&\mbox{ if
$n\in\mathcal{N}(\theta_{b})\setminus\mathcal{N}_{i^{b}}(\theta_{b})$,}\\\
=0,&\mbox{ otherwise.}\end{cases}$ (54)

Let $\lambda_{\min}=\min_{1\leq i\leq k}\lambda_{i}(\theta_{b})$. For
arbitrary $\epsilon>0$, there exists $n_{1}$ such that
$N^{n}_{\text{tot}}(\theta_{b})>2/\lambda_{\min}^{2}\epsilon-1$ for all $n\geq
n_{1}$. If $n\in\mathcal{N}_{i^{b}}(\theta_{b})$ and $n\geq n_{1}$, then
$\Delta_{n}(\theta_{b})<0$ and we get

$\displaystyle\Delta_{n+1}(\theta_{b})$ $\displaystyle=$
$\displaystyle\left(\dfrac{N^{n}_{i^{b}}(\theta_{b})+1}{(N^{n}_{\text{tot}}(\theta_{b})+1)\lambda_{i^{b}}(\theta_{b})}\right)^{2}-\sum_{j\neq
i^{b}}\left(\dfrac{N^{n}_{j}(\theta_{b})}{(N^{n}_{\text{tot}}(\theta_{b})+1)\lambda_{j}(\theta_{b})}\right)^{2}$
(55) $\displaystyle=$
$\displaystyle\left(\dfrac{N^{n}_{i^{b}}(\theta_{b})+1}{(N^{n}_{\text{tot}}(\theta_{b})+1)\lambda_{i^{b}}(\theta_{b})}\right)^{2}-\left(\dfrac{N^{n}_{i^{b}}(\theta_{b})}{(N^{n}_{\text{tot}}(\theta_{b})+1)\lambda_{i^{b}}(\theta_{b})}\right)^{2}+\underbrace{\left(\frac{N^{n}_{\text{tot}}(\theta_{b})}{N^{n}_{\text{tot}}(\theta_{b})+1}\right)^{2}\Delta_{n}(\theta_{b})}_{\text{$\leq
0$}}$ $\displaystyle\leq$
$\displaystyle\dfrac{2N^{n}_{i^{b}}(\theta_{b})+1}{(N^{n}_{\text{tot}}(\theta_{b})+1)^{2}\lambda^{2}_{i^{b}}(\theta_{b})}\leq\dfrac{2}{(N^{n}_{\text{tot}}(\theta_{b})+1)\lambda^{2}_{\min}}<\epsilon.$

Similarly, when
$n\in\mathcal{N}(\theta_{b})\setminus\mathcal{N}_{i^{b}}(\theta_{b})$ and
$n\geq n_{1}$, we can verify that $\Delta_{n+1}(\theta_{b})>-\epsilon$.

We argue there exists $n_{2}$ $(\geq n_{1})$ such that
$-\epsilon<\Delta_{n_{2}}(\theta_{b})<\epsilon$. To see this, suppose
$\Delta_{n_{1}}(\theta_{b})<0$. Selecting $n_{2}=\inf\left\\{n\geq
n_{1}:\Delta_{n}(\theta_{b})\geq 0\right\\}$, we have
$n_{2}-1\in\mathcal{N}_{i^{b}}(\theta_{b})$ as
$\Delta_{n_{2}-1}(\theta_{b})<0$ and $\Delta_{n_{2}}(\theta_{b})>0$. Hence,
(55) shows that $0<\Delta_{n_{2}}(\theta_{b})<\epsilon$. A similar argument
can be made when $\Delta_{n_{1}}(\theta_{b})\geq 0$ to show
$-\epsilon<\Delta_{n_{2}}(\theta_{b})<0$ by choosing $n_{2}=\inf\left\\{n\geq
n_{1}:\Delta_{n}(\theta_{b})<0\right\\}$.

If $-\epsilon<\Delta_{n_{2}}(\theta_{b})<0$, then
$n_{2}\in\mathcal{N}_{i^{b}}(\theta_{b})$ or $n_{2}\in\mathcal{N}(\theta_{c})$
for some $\theta_{c}\neq\theta_{b}$ since $\Delta_{n_{2}}(\theta_{b})<0$
implies that $(i^{b},\theta_{b})$ should be sampled at $n_{2}$ if
$n_{2}\in\mathcal{N}(\theta_{b})$. When
$n_{2}\in\mathcal{N}_{i^{b}}(\theta_{b})$, Lemma 10.2 and (55) imply that we
obtain
$-\epsilon<\Delta_{n_{2}}(\theta_{b})<\Delta_{n_{2}+1}(\theta_{b})<\epsilon$.
When $n_{2}\in\mathcal{N}(\theta_{c})$ for some $\theta_{c}\neq\theta_{b}$,
Lemma 10.2 shows
$-\epsilon<\Delta_{n_{2}}(\theta_{b})=\Delta_{n_{2}+1}(\theta_{b})<0$. In
other words, $\Delta_{n_{2}+1}(\theta_{b})\in(-\epsilon,\epsilon)$ always
holds. We can show an analogous result for when
$0<\Delta_{n_{2}}(\theta_{b})<\epsilon$. Hence, we obtain
$\Delta_{n}(\theta_{b})\in(-\epsilon,\epsilon)$ for all $n\geq n_{2}$ by
induction, which implies $\lim_{n\rightarrow\infty}\Delta_{n}(\theta_{b})=0$.
∎ (Step 4)

Step 5. The pairwise balance condition holds in the limit.

Remark: This proof is applicable for more general algorithms with positive
balance weights not limited to Algorithms 3 and 4. For example, notice that
our sampling rule is identical to the C-OCBA discussed in Gao et al. (2019) if
$W_{i}(\theta_{b})=1$ for all $(i,\theta_{b})\in\Xi\cup\Xi^{\text{adv}}$.
Theorem 2 in Gao et al. (2019) states that the sampling ratios from their
plug-in algorithm achieves convergence of global and pairwise balance
conditions. However, this is stated without a proof (only a brief sketch of
proof is provided). Our proof of Step 5 presented below also proves Theorem 2
in Gao et al. (2019) by simply setting $W_{i}(\theta_{b})=1$ for all
$(i,\theta_{b})\in\Xi\cup\Xi^{\text{adv}}$, which involves expending the
technical results in Chen and Ryzhov (2019a, b).

First, we aim to show that

$\limsup_{n\rightarrow\infty}r_{n}(\theta_{b},\theta_{b^{\prime}})\leq
1,\;\;\text{ for all ${b^{\prime}}\neq b$,}$ (56)

where $r_{n}(\theta_{b},\theta_{{b^{\prime}}}):={\min_{i\neq
i^{b}}W_{i}(\theta_{b})G_{i,n}(\theta_{b})}/{\min_{j\neq
i^{{b^{\prime}}}}W_{j}(\theta_{b^{\prime}})G_{j,n}(\theta_{{b^{\prime}}})}$,
which is the ratio between the smallest LDR functions of $\theta_{b}$ and
$\theta_{{b^{\prime}}}$. For convenience, we define
$\Gamma_{i,n}(\theta_{b})=nW_{i}(\theta_{b})G_{i,n}(\theta_{b})$; note that
$\Gamma_{i,n}(\theta_{b})={\Theta_{\text{a.s.}}}(n)$ as
$0<\liminf_{n\rightarrow\infty}G_{i,n}(\theta_{b})\leq\limsup_{n\rightarrow\infty}G_{i,n}(\theta_{b})<\infty$.
Once (56) holds, we further have
$\liminf_{n\rightarrow\infty}r_{n}(\theta_{b},\theta_{{b^{\prime}}})=1/\limsup_{n\rightarrow\infty}r_{n}(\theta_{{b^{\prime}}},\theta_{b})\geq
1$ by symmetry, which in turn implies
$\lim_{n\rightarrow\infty}r_{n}(\theta_{b},\theta_{b^{\prime}})=1$ for all
$b\neq{b^{\prime}}$ and shows the pairwise balance condition holds.

To show (56), the following two lemmas are introduced. They provide uniform
upper bounds on the increment of scaled LDR, $\Gamma_{i,n}(\theta_{b})$, and
the number of replications allocated to each input parameter, respectively.
Their proofs are presented in Section 11. To state these lemmas, suppose
$n\in\mathcal{N}(\theta_{b})$ and define
$m_{n}(\theta_{b}):=\inf\left\\{\ell>0:n+\ell\in\mathcal{N}(\theta_{b})\right\\}$,
the number of replications until the next time algorithm allocates a
replication to any of $\left\\{(i,\theta_{b})\right\\}_{1\leq i\leq k}$.

###### Lemma 10.3

$\left|\min_{i\neq
i^{b}}\Gamma_{i,n+m_{n}(\theta_{b})}(\theta_{b})-\min_{i\neq
i^{b}}\Gamma_{i,n}(\theta_{b})\right|={O_{\text{a.s.}}}(\sqrt{n\log\log n})$
for $n\in\mathcal{N}(\theta_{b})$.

###### Lemma 10.4

Given $n\in\mathcal{N}(\theta_{b})$, for each ${b^{\prime}}\neq b$, the number
of samples allocated to $\theta_{b^{\prime}}$ is bounded as
$\left|\mathcal{N}(\theta_{b^{\prime}})\cap(n,n+m_{n}(\theta_{b}))\right|={O_{\text{a.s.}}}(\sqrt{n\log\log
n}).$

Because Lemma 10.4 holds for all ${b^{\prime}}\neq b$, it implies that
$m_{n}(\theta_{b})={O_{\text{a.s.}}}(\sqrt{n\log\log n})$. From this, for each
$0\leq t\leq m_{n}(\theta_{b})$, we have
$\Gamma_{i,n+t}(\theta_{b^{\prime}})={\Theta_{\text{a.s.}}}(n+t)={\Theta_{\text{a.s.}}}(n)$
and $\min_{i\neq
i^{b^{\prime}}}\Gamma_{i,n+t}(\theta_{b^{\prime}})={\Theta_{\text{a.s.}}}(n)$.
With these observations and Lemma 10.3, we further obtain

$\displaystyle
r_{n+t}(\theta_{b},\theta_{b^{\prime}})-r_{n}(\theta_{b},\theta_{b^{\prime}})$
$\displaystyle=\frac{\min_{i\neq i^{b}}\Gamma_{i,n+t}(\theta_{b})}{\min_{i\neq
i^{{b^{\prime}}}}\Gamma_{i,n+t}(\theta_{{b^{\prime}}})}-\frac{\min_{i\neq
i^{b}}\Gamma_{i,n}(\theta_{b})}{\min_{i\neq
i^{{b^{\prime}}}}\Gamma_{i,n}(\theta_{{b^{\prime}}})}$
$\displaystyle=\frac{\min_{i\neq i^{b}}\Gamma_{i,n+t}(\theta_{b})\min_{i\neq
i^{{b^{\prime}}}}\Gamma_{i,n}(\theta_{{b^{\prime}}})-\min_{i\neq
i^{{b^{\prime}}}}\Gamma_{i,n+t}(\theta_{{b^{\prime}}})\min_{i\neq
i^{b}}\Gamma_{i,n}(\theta_{b})}{\min_{i\neq
i^{{b^{\prime}}}}\Gamma_{i,n+t}(\theta_{{b^{\prime}}})\min_{i\neq
i^{{b^{\prime}}}}\Gamma_{i,n}(\theta_{{b^{\prime}}})}$
$\displaystyle\leq\frac{\left|\min_{i\neq
i^{b}}\Gamma_{i,n+t}(\theta_{b})-\min_{i\neq
i^{b}}\Gamma_{i,n}(\theta_{b})\right|\min_{i\neq
i^{{b^{\prime}}}}\Gamma_{i,n}(\theta_{{b^{\prime}}})}{\min_{i\neq
i^{{b^{\prime}}}}\Gamma_{i,n+t}(\theta_{{b^{\prime}}})\min_{i\neq
i^{{b^{\prime}}}}\Gamma_{i,n}(\theta_{{b^{\prime}}})}$
$\displaystyle\quad\quad+\frac{\left|\min_{i\neq
i^{{b^{\prime}}}}\Gamma_{i,n+t}(\theta_{{b^{\prime}}})-\min_{i\neq
i^{{b^{\prime}}}}\Gamma_{i,n}(\theta_{{b^{\prime}}})\right|\min_{i\neq
i^{b}}\Gamma_{i,n}(\theta_{b})}{\min_{i\neq