An Upper Confidence Bound Approach to Estimating the Maximum Mean

Before developing estimators for $\mu^{*}$, we highlight some key properties of the GUCB policy. To convey the main idea, for a while we focus on a special case where the exploration rate $\nu_{n}=\log n$. In this case, it has been well known that (see Auer et al. 2002), for system $k$ ($k\neq k^{*}$),

$${\mathbb{E}}[T_{k}(n)]\leq C\log n,$$

for some constant $C$ that depends on $\Delta_{k}$. Then, it can be easily seen that ${\mathbb{E}}[T_{k^{*}}(n)]$, the expected number of times system $k^{*}$ is chosen, is of order $n$ during the first $n$ rounds, because the summation of $T_{k}(n)$’s equals $n$.

In other words, it is expected that among the first $n$ rounds, system $k^{*}$ is chosen in the majority of rounds. Recall that our objective is to estimate $\mu^{*}$, the mean of system $k^{*}$. It is, therefore, reasonable to use the grand average of the samples of all the $n$ rounds as an estimator of $\mu^{*}$, i.e., the GA estimator. Although it uses the samples that are not drawn from system $k^{*}$, the validity of the GA estimator can be justified by the fact that the estimation error due to such samples may phase out when $n$ is sufficiently large, as the number of such samples is negligible compared to those drawn from system $k^{*}$.

Specifically, the GA estimator is given by

$$\widetilde{M}_{n}=\frac{1}{n}\sum_{k=1}^{K}\sum_{j=1}^{T_{k}(n)}X_{k,j},$$

(2)

where $\{X_{k,j},j=1,\ldots,T_{k}(n)\}$ denotes the samples drawn from system $k$ during the first $n$ rounds using the GUCB policy in Algorithm 2.

The GA estimator is not new in the literature, which has served as a key ingredient for algorithms for MDPs and reinforcement learning under the special case that $\nu_{n}=\log n$. However, research on the statistical properties of the estimator has been underdeveloped. To the best of our knowledge, only the asymptotic unbiasedness of the estimator has been established; see, e.g., Kocsis and Szepesvári (2006), while other statistical properties seem to be missing. One of the contributions of this paper is to fill this gap. In particular, we establish strong consistency, asymptotic MSE, and CLT for the GA estimator, which shall be discussed in detail in Section 4.

By taking grand average of the samples of all the systems including those with means smaller than $\mu^{*}$, the GA estimator is naturally negatively biased. While this bias gradually dies out as the sample size $n$ increases, its impact cannot be neglected with finite samples, which has been considered as a main drawback of the GA estimator.

As an attempt to remedy the bias drawback of the GA estimator, we propose a new estimator based on a key observation that the number of samples drawn from system $k^{*}$ usually dominates those from other systems. We therefore define, at the $n$th round,

$$I_{n}^{*}\triangleq\mathop{\arg\max}_{k=1,...,K}T_{k}\left(n\right).$$

(3)

In other words, $I_{n}^{*}$ is the index of the system from which the GUCB policy has sampled most frequently so far. It is reasonable to expect that it is very likely to observe $I_{n}^{*}=k^{*}$. Recall that $\mu^{*}$ is the mean of the performance of system $k^{*}$. It is, therefore, reasonable to propose a new estimator of $\mu^{*}$ by taking sample average of the system indexed by $I_{n}^{*}$, i.e., the LSA estimator. Specifically, the LSA estimator is given by

$$M_{I_{n}^{*}}=\frac{1}{T_{I_{n}^{*}}(n)}\sum_{j=1}^{T_{I_{n}^{*}}(n)}X_{I_{n}^{*},j},$$

(4)

where the samples are drawn using the GUCB policy in Algorithm 2.

The LSA estimator and the maximum estimator share a similar structure, both of them being sample averages of certain identified systems. However, the way of identifying the target system differs significantly. The maximum estimator chooses the system with the largest sample average and uses directly this sample average to estimate $\mu^{*}$. Instead, the LSA estimator chooses the system from which the GUCB policy has sampled most frequently so far, and uses the sample average of this system to estimate $\mu^{*}$. The latter seems to be more intuitively appealing in that it is virtually an unbiased sample-mean estimator conditioning on the event that system $k^{*}$ is correctly identified. Arguably, correctly identifying $k^{*}$ may be relatively easy as the number of rounds GUCB samples from system $k^{*}$ is usually overwhelmingly larger than those from other systems.

Compared to the GA estimator, the LSA estimator avoids taking averages of those systems with smaller means with a high probability. It is thus expected that the LSA estimator may have smaller bias, which shall be shown theoretically in Section 4.

In what follows, we demonstrate how the maximum mean estimators may lead to a new approach to a multiple comparison problem that finds important application in clinical trials.

In a clinical trial, suppose that one is interested in comparing $K-1$ treatments with unknown mean effects $\mu_{1},\ldots,\mu_{K-1}$ to a control treatment with known mean effect $\mu_{0}$, where $K\geq 2$ and a higher mean effect implies a better treatment. In the current practice, a commonly used approach to this comparison problem is to test simultaneously $K-1$ null hypotheses $H_{0,k}$’s against alternatives $H_{1,k}$’s, where

$\displaystyle H_{0,k}:\mu_{k}\leq\mu_{0}\quad v.s.\quad H_{1,k}:\mu_{k}>\mu_{0},\quad k=1,\ldots,K-1.$

(5)

With the above multiple hypotheses framework, a meta null hypothesis is set as $H_{0,k}$’s being true simultaneously for all $k=1,\ldots,K-1$, which is rejected if any of $H_{0,k}$’s is rejected, implying that at least one of the $K-1$ treatments outperforms the control with certain confidence. Suppose that the overall type I error of this test with multiple hypotheses is set to be $\alpha$. In this case, type I errors associated with individual hypotheses cannot be set as $\alpha$, in order to control the overall type I error. Instead, Bonferroni correction is often applied, which assigns a type I error of $\alpha/(K-1)$ to the $K-1$ individual tests. By doing so, it can be verified that the family-wise error rate (FWER), defined as the probability of rejecting at least one true $H_{0,k}$ is not larger than $\alpha$, provided that these individual tests are independent. To see this, suppose that there are $m_{0}$ true null hypotheses that are assumed to be $\{H_{0,1},\ldots,H_{0,m_{0}}\}$ without loss of generality, where $m_{0}$ takes values in $\{1,\ldots,K-1\}$ and is unknown to the researcher. Then,

$\displaystyle\mbox{FWER}=\Pr\left(\cup_{i=1}^{m_{0}}\{\mbox{Rejecting $H_{0,i}$}\}\right)\leq\sum_{i=1}^{m_{0}}\Pr\left(\mbox{Rejecting $H_{0,i}$}\right)=m_{0}{\alpha\over K-1}\leq\alpha.$

While it is easy to implement, Bonferroni correction has been criticized, due to its conservativeness, i.e., a large gap between FWER and $\alpha$ especially when $K$ is large, and the drawback that it usually increases the probability of type II errors. It is thus desirable to develop other procedures for the comparison of treatments that avoid Bonferroni correction.

We note that with the estimators of the maximum mean, a more appealing approach can be developed for this multiple comparison problem. More specifically, in contrast to the multiple hypotheses formulation, we propose a single hypothesis test, with

$\displaystyle H_{0}:\max_{k=0,\ldots,K-1}\mu_{k}=\mu_{0}\quad v.s.\quad H_{1}:\max_{k=0,\ldots,K-1}\mu_{k}>\mu_{0}.$

(6)

Relabelling $\mu_{0}$ as $\mu_{K}$, we note that the GA or LSA estimator, $\widetilde{M}_{n}$ or $M_{I_{n}^{*}}$, can be used to estimate $\max_{k=0,\ldots,K-1}\mu_{k}\equiv\max_{k=1,\ldots,K}\mu_{k}$. Essentially, by comparing $\widetilde{M}_{n}$ (or $M_{I_{n}^{*}}$) to $\mu_{0}$, one may be able to decide whether to reject $H_{0}$ or not. Nevertheless, one needs to know the asymptotic distribution of $\widetilde{M}_{n}$ (or $M_{I_{n}^{*}}$), in order to formalize the test. We shall show in Section 4, more specifically, Theorems 5 and 8, that both $\widetilde{M}_{n}$ and $M_{I_{n}^{*}}$ follow asymptotic normal distributions, thus providing theoretical support for the proposed single hypothesis test. Exact forms of the testing statistics resulting from the GA and LSA estimators are delayed to (11) and (14) in the subsequent section.

In what follows, we establish strong consistency, asymptotic MSE and asymptotic normality for the GA estimator $\widetilde{M}_{n}$. Let $X_{I_{j},j}$ denote the sample drawn at the $j$th round. Note that

$\displaystyle\widetilde{M}_{n}$

$\displaystyle=$

$\displaystyle\frac{1}{n}\sum_{k=1}^{K}T_{k}(n)\bar{X}_{k}[T_{k}(n)]=\frac{1}{n}\sum_{j=1}^{n}X_{I_{j},j}=\frac{1}{n}\sum_{j=1}^{n}\left(X_{I_{j},j}-\mu_{I_{j}}\right)+\frac{1}{n}\sum_{j=1}^{n}\mu_{I_{j}}.$

Define

$\displaystyle Z_{n}\triangleq\sum_{j=1}^{n}\left(X_{I_{j},j}-\mu_{I_{j}}\right).$

Then,

$\displaystyle\widetilde{M}_{n}=\frac{Z_{n}}{n}+\frac{1}{n}\sum_{j=1}^{n}\mu_{I_{j}}.$

(7)

Let ${\cal F}_{n}$ be the $\sigma$-field generated by the first $n$ samples for $n\geq 1$, and ${\cal F}_{0}=\{\Omega,\emptyset\}$. Note that $I_{n}$ is ${\cal F}_{n-1}$-measurable. It follows that for $n\geq 1$,

$\displaystyle{\mathbb{E}}\left[\left.X_{I_{n},n}\right|{\cal F}_{n-1}\right]=\mu_{I_{n}},\ {\rm and\ then}\ {\mathbb{E}}X_{I_{n},n}={\mathbb{E}}\mu_{I_{n}}.$

Therefore, ${\mathbb{E}}Z_{n}=0$ and

$\displaystyle{\mathbb{E}}\left[\left.Z_{n}\right|{\cal F}_{n-1}\right]=Z_{n-1}+{\mathbb{E}}\left[\left.X_{I_{n},n}-\mu_{I_{n}}\right|{\cal F}_{n-1}\right]=Z_{n-1},$

i.e., $Z_{n}$ is a martingale with mean 0.

Note that $X_{k}-\mu_{k}$ is sub-Gaussian with parameter $\gamma_{k}$ for $k=1,...,K$. Because $I_{n+1}$ is ${\cal F}_{n}$-measurable,

$\displaystyle{\mathbb{E}}\left[\left.\exp\{s(X_{I_{n+1},n+1}-\mu_{I_{n+1}})\}\right|{\cal F}_{n}\right]\leq\exp\left\{\frac{\bar{\gamma}^{2}s^{2}}{2}\right\},\mbox{ for }s\in{\mathbb{R}},$

implies that $X_{I_{n+1},n+1}-\mu_{I_{n+1}}$ is sub-Gaussian with parameter $\bar{\gamma}$.

Recall that by the definition of $Z_{n}$,

$\displaystyle Z_{n+1}-Z_{n}=X_{I_{n+1},n+1}-\mu_{I_{n+1}}.$

Then, we have

$\displaystyle{\mathbb{P}}\left(\left.|Z_{n+1}-Z_{n}|>a\right|{\cal F}_{n}\right)\leq 2\exp\left\{-\frac{a^{2}}{2\bar{\gamma}^{2}}\right\},$

by the Hoeffding bound in Lemma A1.

Further note that $\{Z_{n+1}-Z_{n},n=1,2,...\}$ is a martingale difference sequence, the condition of Theorem 2 in Shamir (2011) is thus satisfied. Applying this theorem, for any $\delta>0$, it holds that with probability at least $1-\delta$,

$\displaystyle\frac{Z_{n}}{n}\leq\sqrt{\frac{56\log(1/\delta)}{n/(2\bar{\gamma}^{2})}}=\sqrt{\frac{112\bar{\gamma}^{2}\log(1/\delta)}{n}}.$

Let $\varepsilon=\sqrt{\frac{112\bar{\gamma}^{2}\log(1/\delta)}{n}}$. Then, $\delta=\exp\left\{-\frac{n\varepsilon^{2}}{112\bar{\gamma}^{2}}\right\}$ and

$\displaystyle{\mathbb{P}}\left(\frac{Z_{n}}{n}>\varepsilon\right)\leq\delta=\exp\left\{-\frac{n\varepsilon^{2}}{112\bar{\gamma}^{2}}\right\}.$

Applying the same argument to the sequence $\{-Z_{n},n=1,2,...\}$, symmetrically we have

$\displaystyle{\mathbb{P}}\left(\frac{Z_{n}}{n}<-\varepsilon\right)\leq\delta=\exp\left\{-\frac{n\varepsilon^{2}}{112\bar{\gamma}^{2}}\right\}.$

Therefore,

$\displaystyle{\mathbb{P}}\left(\left|\frac{Z_{n}}{n}\right|>\varepsilon\right)\leq 2\exp\left\{-\frac{n\varepsilon^{2}}{112\bar{\gamma}^{2}}\right\},$

and moreover,

$\displaystyle\sum_{n=1}^{\infty}{\mathbb{P}}\left(\left|Z_{n}/n\right|>\varepsilon\right)<\infty.$

By Borel-Cantelli lemma, $Z_{n}/n\overset{a.s.}{\longrightarrow}0$ as $n\to\infty$. i.e., the first term on the RHS of (7) converges to 0 almost surely.

Moreover, the second term on the RHS of (7) satisfies

$\displaystyle\frac{1}{n}\sum_{j=1}^{n}\mu_{I_{j}}=\frac{1}{n}\sum_{k=1}^{K}T_{k}(n)\mu_{k}\overset{a.s.}{\longrightarrow}\mu_{k^{*}},$

where the convergence follows from Theorem 1 and the continuous mapping theorem. Therefore, by (7), we establish strong consistency of $\widetilde{M}_{n}$ that is summarized as follows.

Theorem 3 ensures that the GA estimator, $\widetilde{M}_{n}$, converges to $\mu^{*}$ with probability 1, as $n\rightarrow\infty$, which serves as an important theoretical guarantee for the estimator. To further understand its probabilistic behavior, it is desirable to conduct error analysis, especially on its bias and variance. In the following theorem, we establish upper bounds for its bias, variance and MSE. The proof of the theorem is provided in Section A.5 of the appendix.

Theorem 4 shows that the MSE of the GA estimator converges to $0$ at a rate of $n^{-1}$. In particular, the variance is the dominant term in the MSE, compared to the square of the bias, which is of order $\nu_{n}^{2}/n^{2}$. This result implies that when the exploration rate takes value in the range $\nu_{n}\in\left[\alpha\log n,n^{1/2-\delta}\right]$ with $0<\delta<1/2$, the bias of the GA estimator is negligible compared to its variance, making it possible to construct CIs by ignoring its bias. Furthermore, as suggested by Theorem 4, the leading term of the variance does not depend on $\nu_{n}$. It is thus asymptotically optimal to choose logarithmic $\nu_{n}$, the slowest possible rate, so as to reduce the asymptotic bias as far as possible. The superiority of this choice of $\nu_{n}$ is also confirmed by the numerical experiments in Section 5.

To provide a formal theoretical support for asymptotically valid CIs and the single hypothesis test in (6), we establish a CLT for $\widetilde{M}_{n}$ in the following theorem, whose proof is provided in Section A.6 of the appendix.

Theorem 5 shows that the GA estimator is asymptotically normally distributed with mean $\mu^{*}$ and variance $\sigma_{k^{*}}^{2}/n$. Based on Theorem 5, an asymptotically valid CI can be constructed for $\mu^{*}$. To do so, a remaining issue is on how to estimate the unknown $\sigma_{k^{*}}^{2}$. In light of the proof of Theorem 5, we estimate $\sigma_{k^{*}}^{2}$ using

$\displaystyle\widetilde{\sigma}_{n}^{2}=\frac{1}{n}\sum_{k=1}^{K}T_{k}(n)\hat{\sigma}_{k}^{2},$

(8)

where for $k=1,...,K$,

$\displaystyle\hat{\sigma}_{k}^{2}=\frac{1}{T_{k}(n)}\sum_{j=1}^{T_{k}(n)}X_{k,j}^{2}-\left[\frac{1}{T_{k}(n)}\sum_{j=1}^{T_{k}(n)}X_{k,j}\right]^{2}.$

(9)

It can be shown that $\widetilde{\sigma}_{n}^{2}$ converges to $\sigma_{k^{*}}^{2}$ in probability, as $n\rightarrow\infty$. This result is summarized in the following proposition, whose proof is provided in Section A.7 of the appendix.

From Theorem 5 and Proposition 2, it follows that

$\displaystyle\sqrt{n}\widetilde{\sigma}_{n}^{-1}\left(\widetilde{M}_{n}-\mu^{*}\right)\Rightarrow N(0,1).$

Then, an asymptotically valid 100$(1-\beta)\%$ CI of $\mu^{*}$ is given by

$\displaystyle\left(\widetilde{M}_{n}-z_{1-\beta/2}\widetilde{\sigma}_{n}/\sqrt{n},\quad\widetilde{M}_{n}+z_{1-\beta/2}\widetilde{\sigma}_{n}/\sqrt{n}\right),$

(10)

where $z_{1-\beta/2}$ is the $1-\beta/2$ quantile of the standard normal distribution.

The CLT in Theorem 5 and the variance estimator $\widetilde{\sigma}_{n}^{2}$ in Proposition 2 also serve as the theoretical foundation for the proposed single hypothesis test in Section 3.3. More specifically, the following testing statistics can be constructed:

$$\widetilde{U}_{n}\triangleq{\sqrt{n}\left(\widetilde{M}_{n}-\mu_{0}\right)\over\widetilde{\sigma}_{n}}.$$

(11)

Note that under either the null or alternative hypothesis, $\max_{k=0,\ldots,K-1}\mu_{k}\geq\mu_{0}$. Therefore, a one-side test can be conducted using the GA-based statistics $\widetilde{U}_{n}$. In particular, the null hypothesis is rejected when $\widetilde{U}_{n}>z_{1-\alpha}$, suggesting that at least one of the considered treatments outperform the control treatment, where $\alpha$ denotes the overall type I error.

Strong consistency of the LSA estimator is stated in the following theorem, whose proof is provided in Section A.8 of the appendix.

Theorem 6 shows that the LSA estimator converges to $\mu^{*}$ as $n\rightarrow\infty$. To better understand its asymptotic error, we establish upper bounds for its asymptotic bias, variance and MSE, which are summarized in the following theorem with proof provided in Section A.9 of the appendix.

Theorem 7 shows that the bias of LSA estimator decays at a rate at least as fast as $\nu_{n}/n^{3/2}$, which is much faster than that of GA estimator (see Theorem 4). Interestingly, the rate of convergence of its bias can be faster than $n^{-1}$ if $\nu_{n}=o(n^{1/2})$. Furthermore, as suggested by Theorem 7, the leading term of the variance does not depend on $\nu_{n}$. It is thus asymptotically optimal to choose logarithmic $\nu_{n}$, the slowest possible rate, so as to reduce the asymptotic bias as far as possible. The superiority of this choice of $\nu_{n}$ is also confirmed by the numerical experiments in Section 5.

It should be remarked that in Theorems 7 and 4, the dominant coefficients of the variance bounds for the LSA and GA estimators, i.e., $\sigma_{k^{*}}^{2}K^{2}$ and $2K\sigma_{k^{*}}^{2}$, are not tight and can be further reduced with more elaborate analysis. These coefficients do not imply that the GA estimator has a smaller variance, and in fact, numerical results in Section 5 suggest the opposite.

Similar to the case for GA, we establish a CLT for the LSA estimator, which is stated in the following theorem. Proof of the theorem is provided in Section A.10 of the appendix.

Theorem 8 shows the asymptotic normality of the LSA estimator, whose rate of convergence is $n^{-1/2}$. From its proof, it can also be seen that if replacing $\sqrt{n}$ by $\sqrt{T_{I^{*}_{n}}(n)}$, the CLT still holds.

While both CLTs for the GA and LSA estimators (in Theorems 5 and 8, respectively) have the same rate of convergence, it should be emphasized that they require different conditions on the exploration rate. Specifically, CLT holds for the LSA estimator so long as the exploration rate $\nu_{n}$ is within the range of $\left[\alpha\log n,n^{1-\delta}\right]$, while CLT of the GA estimator can only be guaranteed under a more restricted condition that $\nu_{n}\in\left[\alpha\log n,n^{{1/2}-\delta}\right]$.

In what follows, we construct an asymptotically valid CI for $\mu^{*}$ based on Theorem 8. To do this, we propose to estimate the unknown variance $\sigma_{k^{*}}^{2}$ by

$\displaystyle\bar{\sigma}_{n}^{2}=\frac{1}{T_{I^{*}_{n}}(n)}\sum_{j=1}^{T_{I^{*}_{n}}(n)}X_{I^{*}_{n},j}^{2}-\left[\frac{1}{T_{I^{*}_{n}}(n)}\sum_{j=1}^{T_{I^{*}_{n}}(n)}X_{I^{*}_{n},j}\right]^{2}.$

Strong consistency of the estimator $\bar{\sigma}_{n}^{2}$ is established in the following proposition, whose proof is provided in Section A.11 of the appendix.

Similar to (10), an asymptotically valid 100$(1-\beta)\%$ CI of $\mu^{*}$ is given by

$\displaystyle\left(M_{I^{*}_{n}}-z_{1-\beta/2}\bar{\sigma}_{n}/\sqrt{n},\quad M_{I^{*}_{n}}+z_{1-\beta/2}\bar{\sigma}_{n}/\sqrt{n}\right).$

(13)

Similar to the discussion in Section 4.1, for the proposed single hypothesis test in Section 3.3, we construct the following LSA-based testing statistics:

$$U_{n}^{*}\triangleq{\sqrt{n}\left(M_{I_{n}^{*}}-\mu_{0}\right)\over\bar{\sigma}_{n}},$$

(14)

and reject the null hypothesis when $U_{n}^{*}>z_{1-\alpha}$, where $\alpha$ denotes the overall type I error.

In the numerical experiments, we observed that the bias of the GA estimator with finite samples could be quite large, especially when the variances of the systems are large, or the exploration rate $\nu_{n}$ is large. This phenomenon is intuitively explained by that there may be too many samples drawn from systems other than $k^{*}$, leading to a significant underestimation of the maximum mean. However, this adverse effect dies out gradually in the long run. As a means to alleviate this adverse effect in finite-sample setting, we propose to choose a certain amount of samples as a warm-up period and discard the warm-up samples in the estimation. Numerical experiments suggest that a small amount of warm-up samples, e.g. $10\%$ of the total sampling budget, may greatly improve the estimation quality of the GA estimator. Throughout our experiments, we set 10% of the budget as the warm-up period for the GA estimator.

By contrast, the LSA estimator is insensitive to the warm-up samples. Whether to keep or to discard these samples does not have much impact on the quality of the LSA estimator. During the implementation, we discard the warm-up period samples for the LSA estimator as well, to be parallel to the setting for the GA estimator.

From an implementation perspective, it has also been documented that slightly tuning the UCB policy by taking into consideration of the variances of the systems may substantially improve the performances of algorithms for multi-armed bandit problems; see, e.g., Auer et al. (2002). We observe the same phenomenon when estimating the maximum mean. Following the same idea of Auer et al. (2002), we propose to replace the UCB, i.e., $\sqrt{2\nu_{n}/T_{k}(n-1)}$, by

$$\sqrt{{2\nu_{n}\over T_{k}(n-1)}V_{k,n}}$$

for $k=1,\ldots,K$ during the implementation, where

$$V_{k,n}=\frac{1}{T_{k}(n-1)}\sum_{i=1}^{T_{k}(n-1)}X_{k,j}^{2}-\bar{X}_{k}^{2}[T_{k}(n-1)]+\sqrt{2\nu_{n}\over T_{k}(n-1)}$$

is an estimate of the UCB of the variance of system $k$.

Consider $K$ systems, with samples from normal distributions with unknown means $(\mu_{1},\mu_{2},\dots,\mu_{K})$. Parameter configuration is similar to that of Fan et al. (2016), who focused on selecting the system with the largest mean rather than estimating the maximum mean. More specifically, the means $(\mu_{1},\mu_{2},\dots,\mu_{K})$ of systems are set to be $\mu_{k}=1.5+0.5k$, $k=1,\ldots,K$, and the standard deviations are set to be $\sigma_{k}=\mu_{k}$ for all systems.

We examine the performances of different estimators when estimating $\mu^{*}$, by considering their biases, standard deviations, MSEs, and RRMSEs, where RRMSE is defined as the percentage of the root MSE to the true maximum mean. Results reported are based on $10^{4}$ independent replications unless otherwise stated.

In the first set of experiments, let $K=20$ with the total sampling budget $n$ ranging between $10^{4}$ and $10^{7}$. We compare biases, MSEs, and coverage probabilities for GA and LSA estimators with respect to (w.r.t.) $\nu_{n}=\log n$, $n^{1/2}$, and $n^{2/3}$. The results are summarized in Figures 1, 2, and 3, illustrating the rates of convergence of absolute biases and MSEs, and coverage probabilities of the 90% CIs, for the GA (left panel) and LSA (right panel) estimators.

• •

  Figure 1 shows that for different choices of $\nu_{n}$, the bias of the LSA estimator is substantially smaller than that of the GA estimator, consistent with Theorems 4 and 7 that the former converges to $0$ at a rate of at least $\nu_{n}/n^{3/2}$, much faster than the latter (with a rate of $\nu_{n}/n$). Furthermore, it is observed that logarithmic $\nu_{n}$ leads to the smallest biases among the three choices, for both the GA and LSA estimators.

• •

  Figure 2 shows that the MSE of the LSA estimator converges approximately to $0$ at a rate of $n^{-1}$ for all choice of $\nu_{n}$. For the GA estimator, the same rate is observed for logarithmic $\nu_{n}$, while the rate tends to be slower when $\nu_{n}=n^{1/2}$ and $\nu_{n}=n^{2/3}$. This coincides with the requirement in Theorems 4 that $\nu_{n}$ must be slower than $n^{1/2}$ for the GA estimator.

• •

  Figure 3 shows the LSA CIs outperform the GA ones substantially in terms of coverage probabilities, especially when sample sizes are relatively small, due to the significantly smaller bias of the LSA estimator. When $\nu_{n}=\log n$, both the GA and LSA CIs achieve the nominal coverage probability with sufficiently large sample sizes. However, the GA CIs fail to do so for $\nu_{n}=n^{1/2}$ and $\nu_{n}=n^{2/3}$ that violate the requirement in Theorem 5.

Given the superiority of logarithmic exploration rate, we set $\nu_{n}=\log n$ in the rest of the experiments in Section 5.

In the second set of experiments, we fix $K=20$, and compare the GA and LSA estimators to four competing estimators given as follows. Consider two maximum average estimators that take the maximum of the sample means of the $K$ systems. In particular, the first is a maximum average estimator based on the GUCB adaptive sampling policy, referred to as adaptive maximum average (AMA) estimator, while the other is the maximum average estimator based on a simple static sampling policy that randomly allocates each sample to one of the systems with equal probabilities, referred to as simple maximum average (SMA) estimator. Moreover, we consider two heuristic estimators in the form of the sample average of certain selected system resulting from existing best arm identification methods. In particular, we consider two Bayesian best arm identification methods, including the top-two probability sampling (TTPS) method of Russo (2020), and the top-two expected improvement (TTEI) of Qin et al. (2017). For these two estimators, we assume known variances of the systems as required by the TTPS and TTEI methods, and set the tuning parameter $\beta$ to be the default value of $0.5$, following Russo (2020) and Qin et al. (2017).

It should be pointed out that, in spite of being reasonable estimators, asymptotic properties of the AMA, TTPS and TTEI estimators, such as their biases, MSEs and central limit theorems, are not known in the literature. The comparison of these estimators in our experiments serves merely to numerically examine their merits, while development of theoretical underpinnings is a topic that deserves future investigation.

Relative biases and standard deviations as percentages of the true maximum mean, and RRMSEs are summarized in Table 1⁴⁴4For cases of $n=10^{4}$ and $10^{5}$, reported errors are estimated based on $10^{5}$ independent replications.. From Table 1, it can be seen that all the considered estimators converge as sample size increases. Among them, LSA and AMA have almost the same biases and variances, and perform the best with the smallest MSEs. Their biases are also smaller or comparable to other estimators. GA and SMA suffer from a disadvantage of significantly larger biases compared to other estimators. TTPS and TTEI have biases comparable to those of LSA and AMA, but with larger variances. It is also interesting to observe that while SMA is positively biased, AMA is negatively biased, suggesting that the dependence structure among the samples drawn from the UCB policy may lead to a negative bias.

Consider a risk measurement application in which we are interested in estimating a coherent risk measure that is the maximum of the expected discounted portfolio loss taken over $K$ probability distributions. Lesnevski et al. (2007) proposed a simulation-based procedure to construct a fixed-width CI for the risk measure, and used an example with $K=256$ to demonstrate the performance of the constructed CIs. For the sake of completeness, we briefly describe this example as follows.

A coherent risk measure can be written as

$\displaystyle\sup_{{\mathbb{P}}\in{\cal P}}{\mathbb{E}}_{\mathbb{P}}[-Y/r],$

where $Y$ is the value of a portfolio at a future time horizon, $1/r$ is a stochastic discount factor which represents the time value of money, and ${\cal P}$ is a set of probability measures simplified by ${\cal P}=\{{\mathbb{P}}_{1},...,{\mathbb{P}}_{K}\}$. Denote $X=-Y/r$ and $\mu_{k}={\mathbb{E}}_{{\mathbb{P}}_{k}}[X]$. Let $X_{i}$ be a random variable whose distribution under a probability measure ${\mathbb{P}}$ is the same as that of $X$ under ${\mathbb{P}}_{k}$, i.e., ${\mathbb{P}}(X_{k}\leq x)={\mathbb{P}}_{k}(X\leq x)$. Then the coherent risk measure is exactly the maximum mean:

$\displaystyle\max_{k=1,...,K}{\mathbb{E}}X_{k}.$

(15)

Consider estimating the coherent risk measure of a portfolio consisting of European-style call and put options written on three assets with price $S_{j}(t)$ at time $t$ for $j=1,2,3$. Let $S_{0}(t)$ denote a market index. All options in the portfolio have a common terminal time $T$. For each of $j=0,1,2,3$, $S_{j}(t)$ follows geometric Brownian motion with drift $d_{j}$ and volatility $\zeta_{j}$, and therefore $\log S_{j}(t)=\log S_{j}(0)+(d_{j}-\zeta_{j}^{2}/2)t+\zeta_{j}W_{j}\sqrt{t}$, where $W_{j}$ is a standard normal random variable. The assets are dependent in the sense that $W_{0}=Z_{0}$, and $W_{j}=\lambda_{j}Z_{0}+\sqrt{1-\lambda_{j}^{2}}Z_{j}$ for $j=1,2,3$, where under probability measure ${\mathbb{P}}$, $Z_{0}$, $Z_{1}$, $Z_{2}$, and $Z_{3}$ are independent standard normal random variables. In this model, $Z_{0}$ corresponds to the market factor common to all assets, while $Z_{1}$, $Z_{2}$, and $Z_{3}$ are idiosyncratic factors corresponding to each individual asset.