# Two-Sample Testing in Reinforcement Learning

Martin Waltz  Chair of Econometrics and Statistics, esp. in the Transport
Sector,
Technische Universität Dresden, 01062 Dresden, Germany  Ostap Okhrin  Chair of
Econometrics and Statistics, esp. in the Transport Sector,
Technische Universität Dresden, 01062 Dresden, Germany

###### Abstract

Value-based reinforcement-learning algorithms have shown strong performances
in games, robotics, and other real-world applications. The most popular
sample-based method is $Q$-Learning. It subsequently performs updates by
adjusting the current $Q$-estimate towards the observed reward and the maximum
of the $Q$-estimates of the next state. The procedure introduces maximization
bias with approaches like Double $Q$-Learning. We frame the bias problem
statistically and consider it an instance of estimating the maximum expected
value (MEV) of a set of random variables. We propose the $T$-Estimator (TE)
based on two-sample testing for the mean, that flexibly interpolates between
over- and underestimation by adjusting the significance level of the
underlying hypothesis tests. A generalization, termed $K$-Estimator (KE),
obeys the same bias and variance bounds as the TE while relying on a nearly
arbitrary kernel function. We introduce modifications of $Q$-Learning and the
Bootstrapped Deep $Q$-Network (BDQN) using the TE and the KE. Furthermore, we
propose an adaptive variant of the TE-based BDQN that dynamically adjusts the
significance level to minimize the absolute estimation bias. All proposed
estimators and algorithms are thoroughly tested and validated on diverse tasks
and environments, illustrating the bias control and performance potential of
the TE and KE.

Keywords:  maximum expected value, two-sample testing, reinforcement learning,
Keywords:  $Q$-learning, estimation bias

##  1 Introduction

Estimating the maximum expected value (MEV) of a set of random variables is a
long-standing statistical problem, including early contributions of Blumenthal
and Cohen (1968), Dudewicz (1971), and Dhariyal et al. (1985). These works
show that for various true underlying distributions an unbiased estimator does
not exist. The problem has recently attained increased attention since it also
arises in modern machine learning domains, most notably in reinforcement
learning (RL). RL aims at finding a policy - a mapping from states to actions
- that maximizes a numerical reward signal (Sutton and Barto, 2018).
Frequently used approaches define a policy-dependent action-value, also called
$Q$-value, for each state-action pair. This value represents the expected sum
of discounted rewards when executing the given action in the given state and
following a specific policy afterward. In particular, the update rule of the
$Q$-Learning (Watkins and Dayan, 1992) algorithm is based on adjusting the
$Q$-estimate for a given state-action pair towards the observed reward and the
maximum of the estimated $Q$-values of the next state. However, this use of
the Maximum Estimator (ME) of the MEV leads to overestimations of action-
values, which are transmitted throughout the update routine (Van Hasselt,
2010). These can damage the learning performance or even lead to failure of
the algorithm (Thrun and Schwartz, 1993), especially when function
approximation is used (Van Hasselt et al., 2016).

Van Hasselt (2010) proposed the Double Estimator (DE), which splits the data
into independent sets, thereby separating the selection and evaluation of a
maximizing value. The corresponding Double $Q$-Learning is a popular choice
among practitioners. Although the DE introduces underestimation bias, Double
$Q$-Learning offers improved robustness and strong performances, especially in
highly stochastic environments. Another crucial contribution is D’Eramo et al.
(2016), in which the Weighted Estimator (WE) alongside Weighted $Q$-Learning
is proposed. From a bias perspective, the estimator builds a compromise
between the overestimating ME and the underestimating DE. However, the WE does
not offer additional flexibility in selecting the level of bias and is
computationally demanding since it requires numerical integration or Monte
Carlo approximation. Notably, those estimators led to modifications of the
Deep $Q$-Networks (DQN, Mnih et al. 2015), which expanded $Q$-Learning to the
deep neural network (DNN) setting and paved the path for the striking success
of RL in recent years (Silver et al., 2017; Vinyals et al., 2019).

D’Eramo et al. (2016), Lan et al. (2020), among others, have shown that both
over- and underestimation of $Q$-values might not always be harmful, depending
on, e.g., the stochasticity of the environment, the difference of the action-
values, the size of the action space, or the time horizon. We argue that a
competitive estimator of the MEV should thus be able to interpolate between
over- and underestimation via an interpretable hyperparameter, enabling it to
deal with a diverse set of environments. Furthermore, the estimator should
obey a variance bound and, for practical application, should be fast and
stable to compute. Fulfilling these criteria, we propose an estimator based on
two-sample testing for the mean, named $T$-Estimator (TE). The idea is to get
a statistically significant statement of whether one mean is truly larger than
others. Consequently, the hyperparameter is the level of significance
$\alpha$, a familiar quantity for researchers and practitioners from diverse
application domains. The ME is shown to be a special case of the TE with
$\alpha=0.5$. Building on the two-sample test statistic, we further consider a
generalization termed $K$-Estimator (KE), which is characterized by a suitable
kernel function and can smooth the discontinuities around testing decisions of
the TE. We theoretically and empirically analyze the TE and KE regarding their
biases and variances, for which general sharp bounds are derived. Using those
newly defined estimators, we propose RL algorithms for the table-based case
and with DNNs as function approximators. Since the two-sample testing
procedure incorporates variance estimates of the involved variables, we employ
an online variance update routine (D’Eramo et al., 2019) in the table-based
scenario, and the framework of the Boostrapped DQN (BDQN, Osband et al. 2016)
in the DNN setting.

The empirical evidence that over- and underestimation of action-values is not
necessarily detrimental to learning performance might be explained by the
connection of $Q$-estimates to the exploration procedure of algorithms (Fox et
al., 2016). However, Fox (2019) and Liang et al. (2021) argue that these
topics should be addressed separately by focusing firstly on unbiased value-
estimation and secondly on improved exploration schemes. We acknowledge this
perspective by additionally proposing an adaptive tuning mechanism for the
significance level $\alpha$ of the TE in the DNN setting with the objective of
minimizing the absolute estimation bias.111We are thankful to two anonymous
reviewers that pushed us to think in this direction with their valuable
comments. The approach complements recent proposals of Dorka et al. (2021) and
Wang et al. (2021). The dynamic adjustment of $\alpha$ is realized by running
partial greedy episodes and comparing $n$-step returns (Sutton and Barto,
2018) with the action-value estimates for the visited state-action pairs.
Furthermore, through learning $\alpha$, we avoid the computational demanding
tuning process of this environment-specific hyperparameter. Finally, we
demonstrate the performance potential of all newly proposed estimators and
algorithms by extensively testing them in various tasks and environments, with
and without function approximation.

The paper is organized as follows: Section 2 formalizes the problem of
estimating the MEV. Section 3 details the proposed estimators and analyzes
them with and without fulfilling the assumption of independently and
identically distributed (iid) data. Section 4 introduces the RL setup and
presents the new temporal-difference algorithms, while Section 5 details the
measurement of estimation bias and introduces the adaptive update mechanism of
$\alpha$. The experiments are shown and thoroughly discussed in Section 6,
with the code being available at: https://github.com/MarWaltz/TUD_RL. Section
7 provides further literature on the state-of-the-art and Section 8 concludes.

##  2 Estimating the Maximum Expected Value

###  2.1 Problem Definition

Let us consider $M\geq 2$ independent random variables $X_{1},\ldots,X_{M}$
with finite expectations
$\mu_{1}=\operatorname{E}(X_{1}),\ldots,\mu_{M}=\operatorname{E}(X_{M})$ and
variances
$\sigma_{1}^{2}=\operatorname{Var}(X_{1}),\ldots,\sigma_{M}^{2}=\operatorname{Var}(X_{M})$.
The corresponding probability density functions (pdfs) and cumulative
distributions functions (cdfs) are denoted $f_{X_{1}},\ldots,f_{X_{M}}$ and
$F_{X_{1}},\ldots,F_{X_{M}}$, respectively. The quantity of interest is the
_maximum expected value_ : $\mu_{*}=\max_{i}\mu_{i}$. Estimation is performed
based on samples $S=\\{S_{1},\ldots,S_{M}\\}$ without knowing moments or
imposing distributional assumptions. The realizations in a sample $S_{i}$ are
assumed to be iid. The unbiased sample mean of $S_{i}$ is denoted
$\hat{\mu}_{i}(S_{i})$, while an estimator of the MEV is referred to as
$\hat{\mu}_{*}(S)$. Throughout the paper, we frequently abbreviate these
notations via $\hat{\mu}_{i}=\hat{\mu}_{i}(S_{i})$,
$\hat{\mu}_{*}=\hat{\mu}_{*}(S)$, and similar, for conciseness. Primary
evaluation criteria of an estimator are its bias
$\operatorname{Bias}(\hat{\mu}_{*})=\operatorname{E}(\hat{\mu}_{*})-\mu_{*}$,
and variance
$\operatorname{Var}(\hat{\mu}_{*})=\operatorname{E}\left\\{[\hat{\mu}_{*}-\operatorname{E}(\hat{\mu}_{*})]^{2}\right\\}$.
These can be aggregated to the mean squared error
$\operatorname{MSE}(\hat{\mu}_{*})=\operatorname{Bias}(\hat{\mu}_{*})^{2}+\operatorname{Var}(\hat{\mu}_{*})$.

###  2.2 Maximum Estimator

The ME $\hat{\mu}^{ME}_{*}$ is the classic approach and takes the maximum of
unbiased mean estimates:

$\hat{\mu}^{ME}_{*}=\max_{i}\hat{\mu}_{i}.$

Denoting the pdf of $\hat{\mu}_{i}$ as $\hat{f}_{i}$ and the corresponding cdf
as $\hat{F}_{i}$, it holds:

$\operatorname{E}\left(\hat{\mu}^{ME}_{*}\right)=\sum_{i=1}^{M}\int_{-\infty}^{\infty}x\hat{f}_{i}(x)\prod\limits_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{M}\hat{F}_{j}(x)dx.$ (1)

The ME is positively biased:
$\operatorname{E}\left(\hat{\mu}^{ME}_{*}\right)\geq\mu_{*}$, see Van Hasselt
(2013). More precisely, following Aven (1985), a general upper bound for the
bias can be given:

$0\leq\operatorname{Bias}\left(\hat{\mu}^{ME}_{*}\right)\leq\sqrt{\frac{M-1}{M}\sum_{i=1}^{M}\operatorname{Var}\left(\hat{\mu}_{i}\right)}.$
(2)

The bias is particularly large when $\mu_{1}\approx\ldots\approx\mu_{M}$.
Furthermore, it can be shown that the variance of $\hat{\mu}^{ME}_{*}$ is
bounded from above:
$\operatorname{Var}\left(\hat{\mu}^{ME}_{*}\right)\leq\sum_{i=1}^{M}\frac{\sigma_{i}^{2}}{|S_{i}|}$,
where $|S_{i}|$ is the sample size of $S_{i}$, see Van Hasselt (2013).

###  2.3 Double Estimator

Van Hasselt (2010) introduced the DE, which is thoroughly analyzed in Van
Hasselt (2013). The key idea is to separate the selection of the maximizing
random variable and the evaluation of its sample mean, which is performed
simultaneously in the ME. The DE splits $S$ randomly into disjoint subsets
$S^{A}=\\{S^{A}_{1},\ldots,S^{A}_{M}\\}$ and
$S^{B}=\\{S^{B}_{1},\ldots,S^{B}_{M}\\}$, guaranteeing that means based on the
subsets are still unbiased. Afterwards, one selects an index which maximizes
the sample mean in $S^{A}$:
$a^{*}\in\\{i\mid\hat{\mu}_{i}(S^{A})=\max_{j}\hat{\mu}_{j}(S^{A})\\}$. The DE
is defined by evaluating $a^{*}$ on $S^{B}$:
$\hat{\mu}^{DE}_{*}(S)=\hat{\mu}_{a^{*}}(S^{B})$. Consequently, one can
perform the same procedure with $S^{A}$ and $S^{B}$ switched to get a second
DE estimate. Averaging both DE estimates yields the 2-fold Cross-Validation
estimator (CVE) $\hat{\mu}^{CVE}_{*}$, which has a reduced variance in
comparison to a single DE estimate. The expectations of the DE and the CVE are
equal since both DE estimates (for $S^{A}$ and $S^{B}$) have identical
expectations:

$\displaystyle\operatorname{E}\left(\hat{\mu}^{CVE}_{*}\right)=\operatorname{E}\left(\hat{\mu}^{DE}_{*}\right)$
$\displaystyle=\sum_{i=1}^{M}\operatorname{E}\left[\hat{\mu}_{i}(S^{B})\right]P(i=a^{*})$
$\displaystyle=\sum_{i=1}^{M}\operatorname{E}\left[\hat{\mu}_{i}(S^{B})\right]\int_{-\infty}^{\infty}\hat{f}_{i}^{A}(x)\prod\limits_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{M}\hat{F}_{j}^{A}(x)dx,$ (3)

where $\hat{f}_{i}^{A}$ and $\hat{F}_{i}^{A}$ are the cdf and pdf of
$\hat{\mu}_{i}(S^{A})$, respectively. Comparing (3) with (1), it appears that
(1) performs integration over $x$, while the corresponding term
$\operatorname{E}[\hat{\mu}_{i}(S^{B})]$ in (3) is outside the integral, so
that probability and value are independent. Van Hasselt (2010) showed that the
DE is prone to underestimation:
$\operatorname{E}\left(\hat{\mu}^{DE}_{*}\right)\leq\mu_{*}$, because it might
attribute non-zero selection probability to non-maximum variables.
Furthermore, Van Hasselt (2013) conjectures the following lower bound for the
bias:

$-\frac{1}{2}\left(\sqrt{\sum_{i=1}^{M}\frac{\sigma_{i}^{2}}{|S_{i}^{A}|}}+\sqrt{\sum_{i=1}^{M}\frac{\sigma_{i}^{2}}{|S_{i}^{B}|}}\right)<\operatorname{Bias}(\hat{\mu}^{DE}_{*})\leq
0,$

while the variance of the CVE is shown to be bounded as the ME:
$\operatorname{Var}\left(\hat{\mu}^{CVE}_{*}\right)\leq\sum_{i=1}^{M}\frac{\sigma_{i}^{2}}{|S_{i}|}$.
Worth mentioning is that the variance of the CVE is not necessarily half the
variance of the DE, as there is non-zero covariance between the two DE
estimates, see the example in Appendix A. Throughout the experiments, we
follow D’Eramo et al. (2021) and use the CVE instead of the DE whenever
possible.

###  2.4 Weighted Estimator

D’Eramo et al. (2016) and D’Eramo et al. (2021) introduced the Weighted
Estimator (WE) for the MEV, which is a weighted mean of all sample averages.
Each weight corresponds to the probability of $\hat{\mu}_{i}$ being larger
than all other means:

$\hat{\mu}^{WE}_{*}=\sum_{i=1}^{M}w_{i}\hat{\mu}_{i}=\sum_{i=1}^{M}P\left(\hat{\mu}_{i}=\max_{j}\hat{\mu}_{j}\right)\hat{\mu}_{i}.$

Since the probabilities depend on the unknown mean distributions
$\hat{f_{i}}$, the authors propose a Gaussian approximation based on the
central limit theorem:

$\hat{\mu}^{WE}_{*}=\sum_{i=1}^{M}\hat{\mu}_{i}\int_{-\infty}^{\infty}\tilde{f}_{i}(x)\prod\limits_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{M}\tilde{F}_{j}(x)dx,$ (4)

where $\tilde{f}_{i}$ is the Gaussian pdf with mean $\hat{\mu}_{i}$ and
variance $\frac{\hat{\sigma}_{i}^{2}}{|S_{i}|}$. The unbiased estimate of
${\sigma}_{i}^{2}$ is denoted $\hat{\sigma}_{i}^{2}$ and $|S_{i}|$ refers to
the sample size. Crucially, the bias of the WE is bounded by the ME and DE:

$\operatorname{Bias}(\hat{\mu}^{DE}_{*})\leq\operatorname{Bias}(\hat{\mu}^{WE}_{*})\leq\operatorname{Bias}(\hat{\mu}^{ME}_{*}),$

while it exhibits the same variance bound:
$\operatorname{Var}\left(\hat{\mu}^{WE}_{*}\right)\leq\sum_{i=1}^{M}\frac{\sigma_{i}^{2}}{|S_{i}|}$.
Thus, the bias of the WE might be positive or negative, depending on the
distribution of the random variables. A drawback of this estimator lies in
increased computation time since calculating the integrals in (4) is a
demanding process. Tackling this issue, D’Eramo et al. (2021) propose to use
Monte Carlo approximations instead, and we follow this approach when computing
the WE in the experiments. The Monte Carlo sample sizes in those cases are set
to 100.

##  3 Two-Sample Testing-based Estimators

###  3.1 T-Estimator

To create a flexible estimator which:

  1. (a)

is able to interpolate between over- and underestimation,

  2. (b)

obeys a variance bound similar to the ME,

  3. (c)

has an interpretable hyperparameter,

  4. (d)

is fast and easy to compute,

we propose a procedure based on one-sided two-sample testing for the mean.
Generally, for two random variables $X_{1},X_{2}$, we consider the hypothesis
$H_{0}:\mu_{1}\geq\mu_{2}$. The test statistic is constructed as follows
(Wackerly et al., 2008):

$T=\frac{\hat{\mu}_{1}-\hat{\mu}_{2}}{\sqrt{\frac{\hat{\sigma}^{2}_{1}}{|S_{1}|}+\frac{\hat{\sigma}^{2}_{2}}{|S_{2}|}}},$

where $\hat{\sigma}^{2}_{1}$ and $\hat{\sigma}^{2}_{2}$ are unbiased estimates
of the variances $\sigma^{2}_{1}$ and $\sigma^{2}_{2}$ of $X_{1}$ and $X_{2}$,
respectively. $|S_{1}|$ and $|S_{2}|$ denote the corresponding sample sizes.
If the realization of $T$ is smaller than the $\alpha$-quantile of the
standard normal distribution, $z_{\alpha}=\Phi^{-1}(\alpha)$, hypothesis
$H_{0}$ is rejected. The use of the normal distribution as the asymptotic
distribution of the test statistics for $H_{0}$ can be justified via the
central limit theorem, since it holds
$\sqrt{|S_{i}|}\frac{\hat{\mu}_{i}-\mu_{i}}{\sigma_{i}}\xrightarrow[|S_{i}|\rightarrow\infty]{d}\mathcal{N}(0,1)$
for $i=1,2$, and using Slutsky’s theorem (Casella and Berger, 2002), as
$\hat{\sigma}_{i}$ converges almost surely to $\sigma_{i}$.

Based on this test, the following procedure for estimating the MEV is
proposed: First, we consider the complete set of indices
$\mathcal{L}=\\{1,\ldots,M\\}$ and select an index that corresponds to a
variable with the value of the ME:
$i^{*}\in\\{i\mid\hat{\mu}_{i}=\max_{j}\hat{\mu}_{j}\\}$. Second, we test for
all $i\in\mathcal{L}$ the $H_{0}$: $\mu_{i}\geq\mu_{i^{*}}$. If $H_{0}$ is
rejected for some $i^{\prime}$, we assert $\mu_{i^{\prime}}<\mu_{i^{*}}$ and
remove variable index $i^{\prime}$ from the index set:
$\mathcal{L}\leftarrow\mathcal{L}\char 92\relax\\{i^{\prime}\\}$. Third, we
average the remaining $\\{\hat{\mu}_{i}\mid i\in\mathcal{L}\\}$. Compactly
written:

$\hat{\mu}^{TE}_{*}(\alpha)=\left[\sum_{i=1}^{M}\mathcal{I}\left(\frac{\hat{\mu}_{i}-\hat{\mu}^{ME}_{*}}{\sqrt{\frac{\hat{\sigma}^{2}_{i}}{|S_{i}|}+\frac{\hat{\sigma}^{2}_{i^{*}}}{|S_{i^{*}}|}}}\geq
z_{\alpha}\right)\right]^{-1}\sum_{i=1}^{M}\mathcal{I}\left(\frac{\hat{\mu}_{i}-\hat{\mu}^{ME}_{*}}{\sqrt{\frac{\hat{\sigma}^{2}_{i}}{|S_{i}|}+\frac{\hat{\sigma}^{2}_{i^{*}}}{|S_{i^{*}}|}}}\geq
z_{\alpha}\right)\hat{\mu}_{i},$ (5)

where $\mathcal{I}(\cdot)$ is the indicator function. We refer to (5) as _T-
Estimator_ (TE). In simple words: TE averages the means of all variables,
which are statistically not smaller than the one of the ME. Consequently, the
selection is a binary decision of rejection or non-rejection of the underlying
hypothesis. A key aspect of the TE is the consideration of the values of the
sample means together with their uncertainties, expressed by variances.
Asymptotically, meaning $|S_{i}|\rightarrow\infty$ for $i=1,\ldots,M$, the TE
follows a normal distribution since it is an average of asymptotically normal
distributed variables.

The hyperparameter is the significance level $\alpha$, which is an
interpretable quantity for practitioners and researchers, and is naturally
restricted to $\alpha\in(0,0.5]$. One can directly determine the extreme case
on the upper domain limit:
$\hat{\mu}^{TE}_{*}(\alpha=0.5)=\hat{\mu}^{ME}_{*}$. Thus, the ME is a special
case of the TE, being prone to overestimation bias with the bounds given in
Section 2.2. By reducing $\alpha$, we reduce the bias since we tend to non-
reject $H_{0}$ for smaller sample means. If one would consider a significance
level of zero, the TE would collapse into the Average Estimator (AE):
$\hat{\mu}^{AVG}_{*}=M^{-1}\sum_{i=1}^{M}\hat{\mu}_{i}$. Imagaw and Kaneko
(2017) provide a similar definition in a multi-armed bandit context. The AE
has low variance:
$\operatorname{Var}(\hat{\mu}^{AVG}_{*})=M^{-2}\sum_{i=1}^{M}\frac{\sigma_{i}^{2}}{|S_{i}|}$,
but severe negative bias:
$\operatorname{Bias}(\hat{\mu}^{AVG}_{*})=-M^{-1}\sum_{i=1}^{M}(\max_{j}\mu_{j}-\mu_{i})$.
However, we do not include $\alpha=0$ in our definition domain for the TE
since 1) it is statistically not reasonable to consider such hypothesis tests
and 2) the uncertainties of the sample means, quantified through variances,
are not present anymore.

In general, we flexibly interpolate between under- and overestimation by
selecting a lower or larger level of significance. The following lemma
contains the sharp bias bounds of the TE.

######  Lemma 1.

For $\alpha\in(0,0.5]$, it holds:

$\frac{1}{2}\left[\min_{i}\mu_{i}-\max_{i}\mu_{i}-\sqrt{\frac{M-1}{M}\sum_{i=1}^{M}\mathrm{Var}\left(\hat{\mu}_{i}\right)}\right]\leq\mathrm{Bias}\left[\hat{\mu}^{TE}_{*}(\alpha)\right]\leq\mathrm{Bias}(\hat{\mu}^{ME}_{*}).$

###### Proof.

The upper bound is straightforward since the TE is a weighted average of
sample means, while the ME is the extreme case of weighting the maximum sample
mean with one. Regarding the lower bound, we use that per construction:

$\hat{\mu}^{TE}_{*}(\alpha)\geq\frac{\max_{i}\hat{\mu}_{i}+\min_{i}\hat{\mu}_{i}}{2}.$
(6)

To see this, we first note that the numerator of the test statistics in (5) is
always zero for the ME, leading to value one of the corresponding indicator
function for all $\alpha\in(0,0.5]$. However, since the test statistics for
index $i$ positively correlates with $\hat{\mu}_{i}$ and
$\hat{\sigma}_{i}^{2}$, extreme variance scenarios are possible in which the
index with the minimum mean yields the only other non-rejected hypothesis.

Building expectations of (6), we have:

$\displaystyle\operatorname{E}\left[\hat{\mu}^{TE}_{*}(\alpha)\right]$
$\displaystyle\geq\frac{1}{2}\left[\operatorname{E}\left(\max_{i}\hat{\mu}_{i}\right)+\operatorname{E}\left(\min_{i}\hat{\mu}_{i}\right)\right]$
$\displaystyle\geq\frac{1}{2}\left[\max_{i}\mu_{i}+\min_{i}\mu_{i}-\sqrt{\frac{M-1}{M}\sum_{i=1}^{M}\mathrm{Var}\left(\hat{\mu}_{i}\right)}\right],$

where the last line uses (2) and the bound for the minimum sample average of
Aven (1985). The bias follows immediately. ∎

Moreover, a general expression for the expectation of the TE for arbitrary $M$
can be derived for known true variances.

######  Lemma 2.

The expectation of the TE is:

$\displaystyle\operatorname{E}$
$\displaystyle[\hat{\mu}^{TE}_{*}(\alpha)]=\sum_{i=1}^{M}\int_{-\infty}^{\infty}\int_{-\infty}^{x_{i}}\cdots\int_{-\infty}^{x_{i}}$
$\displaystyle\left[\sum_{j=1}^{M}\mathcal{I}\left(\frac{x_{j}-x_{i}}{\theta_{ij}}\geq
z_{\alpha}\right)\right]^{-1}\left[\sum_{j=1}^{M}\mathcal{I}\left(\frac{x_{j}-x_{i}}{\theta_{ij}}\geq
z_{\alpha}\right)x_{j}\right]\left[\prod_{j=1}^{M}\hat{f}_{j}(x_{j})\right]\left[\prod\limits_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{M}dx_{j}\right]dx_{i},$ (7)

where
$\theta_{ij}=\sqrt{\frac{\sigma_{i}^{2}}{|S_{i}|}+\frac{\sigma_{j}^{2}}{|S_{j}|}}$,
$\hat{f_{j}}$ is the pdf of $\hat{\mu}_{j}$, being asymptotically normal, and
$\sigma_{i}^{2}$ are known.

That this expression holds true appears by generalizing the example with $M=2$
detailed in Appendix A to higher dimensions. Furthermore, the TE is consistent
for the MEV since, with increasing sample size, each sample mean approaches
its population mean, the ME approaches the true MEV, and the variances of the
means tend to zero. Consequently, all tests except for the true MEV variable
will reject the $H_{0}$, and only the MEV will be left. Note that the TE
depends on performing multiple hypothesis tests, which can increase type I
error. However, for simplicity, we stick to the procedure as presented above
and do not pursue multiple testing corrections, e.g., a Bonferroni correction
(Armstrong, 2014). Regarding variance, the TE shares the common overly
pessimistic variance bound, while the proof relies on the TE being a weighted
average of means and is similar to D’Eramo et al. (2016):

######  Lemma 3.

For $\alpha\in(0,0.5]$, it holds:
$\mathrm{Var}\left[\hat{\mu}^{TE}(\alpha)\right]\leq\sum_{i=1}^{M}\frac{\sigma_{i}^{2}}{|S_{i}|}$.

###  3.2 K-Estimator

The TE is a flexible estimator fulfilling our initial requirements, but can be
put in an even broader context. Crucially, the TE involves indicator
functions, whose derivatives with relation to $\alpha$ are either zero or not
existent, leading to the binary decision of exclusion and non-exclusion of a
particular mean. Avoiding this behaviour, we propose to apply the standard
Gaussian cdf $\Phi$ directly to the test statistics and use the resulting
values as a smoothed weighting:

$\hat{\mu}^{\Phi}_{*}=\left[\sum_{i=1}^{M}\Phi\left(T_{i}\right)\right]^{-1}\sum_{i=1}^{M}\Phi\left(T_{i}\right)\hat{\mu}_{i},\quad\text{where}\quad
T_{i}=\frac{\hat{\mu}_{i}-\hat{\mu}^{ME}_{*}}{\sqrt{\frac{\hat{\sigma}^{2}_{i}}{|S_{i}|}+\frac{\hat{\sigma}^{2}_{i^{*}}}{|S_{i^{*}}|}}}.$
(8)

In fact, it can be generalized even further by considering a weighting kernel
$\kappa(\cdot)$:

$\hat{\mu}^{KE}_{*}=\left[\sum_{i=1}^{M}\kappa\left(T_{i}\right)\right]^{-1}\sum_{i=1}^{M}\kappa\left(T_{i}\right)\hat{\mu}_{i}.$
(9)

with $\kappa:\mathbb{R}^{-}\rightarrow\mathbb{R}^{+}$, where
$\mathbb{R}^{-}=(-\infty;0]$ and
$\mathbb{R}^{+}=\mathbb{R}\setminus\mathbb{R}^{-}\cup\left\\{0\right\\}$. We
require that $\kappa(\cdot)$ is monotonically increasing to build a reasonable
kernel since $T_{i}\leq 0$, $\forall i=1,\ldots,M$, and that
$\lim\limits_{T_{i}\rightarrow-\infty}\kappa(T_{i})=0$ for consistency.
Similar kernel functions are considered in Mammen (1991) for isotonic
regressions. We refer to (9) as the _K-Estimator_ (KE). Crucially, the
hyperparameter of the KE is not a fixed scalar anymore, but the specification
of $\kappa(\cdot)$. For example, setting
$\kappa(T_{i};\alpha)=\mathcal{I}(T_{i}\geq z_{\alpha})$ results in obtaining
the TE. Further options for $\kappa(\cdot)$ are listed in Table 1.
Additionally, more flexible parametrized specifications are available.
Consider, for example, the cdf of the beta distribution
$\mathcal{B}_{\mathfrak{a},\mathfrak{b}}$ with two shape parameters
$\mathfrak{a}$, $\mathfrak{b}$. Although the latter is naturally defined on
$[0,1]$, one could simply scale and shift it to, e.g., $[-5,0]$ to generate a
more valid specification. This particular case is used in the example of
Section 3.3. Apart from that, we primarily apply the standard Gaussian cdf
throughout the paper.

Kernel |  $\kappa(T)$ for $T\leq 0$
---|---
cdf: Gaussian $\Phi_{\lambda}$ | $\int_{-\infty}^{T}\frac{1}{\sqrt{2\pi\lambda^{2}}}\exp{\left\\{-\frac{1}{2}\left(\frac{t}{\lambda}\right)^{2}\right\\}}dt$
cdf: $t$-distribution $t_{\nu}$ | $\int_{-\infty}^{T}\frac{\Gamma(\frac{\nu+1}{2})}{\sqrt{\nu\pi}\Gamma(\frac{\nu}{2})}\left(1+\frac{t^{2}}{\nu}\right)^{-\frac{\nu+1}{2}}dt$
Epanechnikov | $\frac{3}{4}(1-T^{2})\mathcal{I}(|T|\leq 1)$
Laplace | $\frac{1}{2}\exp{(-|T|)}$
Triangle | $(1-|T|)\mathcal{I}(|T|\leq 1)$
Table 1: Exemplary kernel functions. $\lambda$ is the standard deviation of
the Gaussian cdf, $\nu$ are the degrees of freedom of the $t$-distribution,
and $\Gamma(x)=\int_{0}^{\infty}t^{x-1}\exp(-t)dt$ denotes the gamma function.
We abbreviate the standard Gaussian kernel $\Phi_{\lambda=1}$ with $\Phi$.

Bias and variance of the KE depend on the chosen kernel specification, but the
bounds of the TE are still valid as long as the kernel function fulfills the
requirements stated above.

######  Corollary 1.

For the KE, it holds:

$\frac{1}{2}\left[\min_{i}\mu_{i}-\max_{i}\mu_{i}-\sqrt{\frac{M-1}{M}\sum_{i=1}^{M}\mathrm{Var}\left(\hat{\mu}_{i}\right)}\right]\leq\mathrm{Bias}\left(\hat{\mu}^{KE}_{*}\right)\leq\mathrm{Bias}(\hat{\mu}^{ME}_{*}).$

###### Proof.

The KE cannot exceed the ME, thus the upper bound holds. For the lower bound,
we note the same relationship as for the TE:

$\hat{\mu}^{KE}_{*}\geq\frac{\max_{i}\hat{\mu}_{i}+\min_{i}\hat{\mu}_{i}}{2}.$

This follows since the sample mean corresponding to the ME is per construction
weighted with $\kappa(0)$ (before normalization). Simultaneously, for extreme
variance scenarios, it is possible that the weight of the minimum sample mean
tends to $\kappa(0)$ as well, while the weights of the remaining means tend to
$\lim\limits_{T_{i}\rightarrow-\infty}\kappa(T_{i})$, which is required to be
zero. The remaining steps are identical to the proof of Lemma 1. ∎

######  Corollary 2.

For the KE, it holds:
$\mathrm{Var}\left(\hat{\mu}^{KE}\right)\leq\sum_{i=1}^{M}\frac{\sigma_{i}^{2}}{|S_{i}|}$.

The proof is again similar to D’Eramo et al. (2016). Finally, we can
generalize (2) of the TE for the KE, describing the expectation for known
variances.

######  Corollary 3.

The expectation of the KE is:

$\displaystyle\operatorname{E}[\hat{\mu}^{KE}_{*}]$
$\displaystyle=\sum_{i=1}^{M}\int_{-\infty}^{\infty}\int_{-\infty}^{x_{i}}\cdots\int_{-\infty}^{x_{i}}$
$\displaystyle\left[\sum_{j=1}^{M}\kappa\left(\frac{x_{j}-x_{i}}{\theta_{ij}}\right)\right]^{-1}\left[\sum_{j=1}^{M}\kappa\left(\frac{x_{j}-x_{i}}{\theta_{ij}}\right)x_{j}\right]\left[\prod_{j=1}^{M}\hat{f}_{j}(x_{j})\right]\left[\prod\limits_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{M}dx_{j}\right]dx_{i},$

where
$\theta_{ij}=\sqrt{\frac{\sigma_{i}^{2}}{|S_{i}|}+\frac{\sigma_{j}^{2}}{|S_{j}|}}$,
$\hat{f_{j}}$ is the pdf of $\hat{\mu}_{j}$, being asymptotically normal, and
$\sigma_{i}^{2}$ are known.

Finally, we highlight three main differences between the KE (including the TE
as a special case) and the WE of Section 2.4 since both are constructed as a
weighted sum of sample means. First, the weights of the WE are probabilities,
while the weights of the KE do not have a probabilistic interpretation.
Second, while the KE allows for a multitude of specifications, the WE is not
tunable and thus cannot be adjusted to a given scenario in a practical
problem. Third, the computation of the WE requires integration or Monte Carlo
approximation schemes, while the KE’s computation is extremely fast.

###  3.3 Example (iid): Bias, Variance, and MSE

To further analyze the behaviour of the proposed estimators, we assume a
similar setup to D’Eramo et al. (2016) with $M=2$. Precisely, we consider
Gaussian random variables $X_{1}\sim\mathcal{N}(\mu_{1},\sigma^{2})$,
$X_{2}\sim\mathcal{N}(\mu_{2},\sigma^{2})$, where $\sigma^{2}=100$ is the
common known variance, and we assume to have $|S_{1}|=|S_{2}|=100$
observations of each variable. As assumed throughout prior sections, the
realizations in a sample $S_{i}$ are iid. We fix $\mu_{2}=0$ and compute bias,
variance, and MSE for different $\mu_{1}\in[0,5]$. For completeness, we report
the analytic forms for the expectation and variance of the estimators in this
case in Appendix A. For the TE, we select significance levels
$\alpha\in\\{0.05,0.10,0.15\\}$ and for the KE, we analyze the standard
Gaussian kernel $\Phi$, the Epanechnikov kernel, and the shifted and scaled
$\mathcal{B}_{\mathfrak{a},\mathfrak{b}}$ cdf kernel with $\mathfrak{a}=2$,
$\mathfrak{b}=0.5$ as described above. Results are displayed in Figures 1 and
2.

Figure 1: Comparison of the ME, DE, and TE with level of significance in
parentheses. Figure 2: Comparison of the ME, DE, and KE with kernel in
parentheses.

Regarding the TE, we see how the bias decreases with a smaller significance
level. In general, the estimator operates between the biases of ME and DE,
although the bias of the DE is not necessarily a lower bound (see
$\alpha=0.05$). In the mean equality case, $\mu_{1}=\mu_{2}$, the TE can avoid
the significant overestimation of the ME while having only slightly increased
variance. Consequently, the TE outperforms the conventional competitors for
all considered significance levels under the MSE criterion for
$\mu_{1}=\mu_{2}$. However, if the difference of the true expectations of the
random variables is large, all estimators become unbiased. In this scenario,
the ME is the best choice due to its low variance. Regarding the KE in Figure
2, we see that the standard Gaussian and Epanechnikov kernels can achieve a
desirable balance between under- and overestimation while maintaining a
smaller variance than the considered TE.

However, the chosen specification of the beta distribution appears sub-optimal
for this particular problem because of its rather strong underestimation for
large $\mu_{1}-\mu_{2}$. To find a better fitting parametrization, we
numerically solved the optimization problem of minimizing the squared bias for
the depicted range of $\mu_{1}-\mu_{2}$ over the parameters $\mathfrak{a}$,
$\mathfrak{b}$ of the $\mathcal{B}_{\mathfrak{a},\mathfrak{b}}$ kernel. To
enable comparability between the estimators, we run the identical optimization
for the parameter $\lambda$ of the Gaussian kernel $\Phi_{\lambda}$ (deviating
from the unit variance specification) and the level of significance of the TE.
The optimized kernel functions alongside specifications from Figures 1, 2 are
depicted in Figure 3. The functions are normalized to $[0,1]$ by division
through $\kappa(0)$ of the respective kernel. Optimizing the standard
deviation of the Gaussian cdf yields $\lambda\approx 0.83$, which is close to
the unit variance specification. On the other hand, the bias-optimal value for
the significance level of the TE is $\approx 0.14$. Regarding the
$\mathcal{B}_{\mathfrak{a},\mathfrak{b}}$ specification, one needs to recall
that the beta kernel is capable of approximating both the optimized TE and the
optimized KE with the (non-standard) Gaussian cdf kernel. Following Figure 3,
the optimized TE is favorable in this scenario since the optimized beta cdf
yields a non-unique, zero-variance solution, which is in line with the
optimized TE.

Overall, we have seen through this investigation that both the TE and the KE
can achieve flexible trade-offs between bias and variance in the estimation of
the MEV. Considering typical levels of significance like $0.05$, $0.10$, or
$0.15$ in the TE builds a robust estimator, further enhanced by accurately
specifying a suited KE.

Figure 3: Original kernel functions and optimized specification for minimizing
the squared bias in Figures 1 and 2.

###  3.4 Example (non-iid): Bias and Variance

The estimators of the MEV will be transferred to solving sequential decision
problems in Section 4, and thus auto-correlations, time dependencies, and
decision-making play an essential role in the algorithms. We illustrate the
limitations of the estimators by applying them to auto-correlated data. First,
consider the following two processes:

$\displaystyle X_{i,t}$ $\displaystyle=(1-\rho)\mu_{i}+\rho
X_{i,t-1}+\varepsilon_{i,t},$ (10) $\displaystyle\varepsilon_{i,t}$
$\displaystyle\sim\mathcal{N}\left\\{0,(1-\rho^{2})\sigma^{2}\right\\},$

with $X_{i,0}=\mu_{i}$ for $i=1,2$ and time steps $t=0,1,\ldots,T$. The given
specification yields the unconditional moments:
$\operatorname{E}[X_{i,t}]=\mu_{i}$ and
$\operatorname{Var}[X_{i,t}]=\sigma^{2}$. Thus, by setting $\mu_{1}=1$,
$\mu_{2}=0$, and $\sigma^{2}=100$, we obtain a special case of the scenarios
studied in Section 3.3 while introducing auto-correlation in the samples. The
latter can be tuned via the parameter $\rho\in[0,1]$, and further details on
such time series can be found in Tsay (2010).

Second, in analogy to the temporal-difference algorithms considered later in
Section 4, we compute exponentially weighted mean estimates based on (10):

$\hat{\mu}_{i,t}=\hat{\mu}_{i,t-1}+\tau\left(X_{i,t}-\hat{\mu}_{i,t-1}\right),$
(11)

where we initialize $\hat{\mu}_{i,0}=\mu_{i}$ for $i=1,2$ and consider time
steps $t=0,1,\ldots,T$ with learning rate $\tau$. Transferred to RL, the
$X_{i,t}$ are - potentially correlated - returns, while $\hat{\mu}_{i,t}$
corresponds to the estimate of the $Q$-value. Figure 4 illustrates exemplary
realizations of (10) and (11).

Figure 4: Realizations of the processes in (10) and (11) for varying $\rho$
and $\tau$. While $\rho=0.0$ represents independent random noise, the
dependence structure is visible for larger $\rho$. The parameter $\tau$
controls the smoothness of the mean estimate. Figure 5: Estimator comparison
in the non-iid setting based on kernel-density estimates for $10^{6}$ runs.
The black dashed line is the true MEV.

We consider $\rho=\\{0.0,0.3,0.7,0.95\\}$ and $\tau=\\{0.1,0.3\\}$ with time
horizon $T=100$, and compare the ME, DE, TE, and KE. After simulating (10) and
(11), we use the final mean estimate to compute the MEV, where we use the
sample variance of the time series $\hat{\mu}_{i,t}$ as the required variance
estimate for the TE and KE. Regarding the DE, we simulate for each variable
($i=1,2$) two mean estimate processes with $T=50$, where we use one for index
selection and one for evaluation as described in Section 2.3. The simulation
and estimation procedure is repeated $10^{6}$ times and we compute a kernel
density estimate (Gaussian kernel, see Sheather 2004) on the MEV estimates.
Figure 5 displays the results.

We see the behavior of the iid-scenario reflected for $\rho=0.0$. The ME
overestimates the MEV in all cases (with the peak of the density being
slightly to the right), while the DE tends to underestimate it (with the peak
being at the left). Similar to Section 3.3, the TE and KE can balance the two
conventional approaches, and the relative ordering between the estimators
remain with increasing auto-correlation. However, all estimators exhibit large
variances for the extreme case $\rho=0.95$, and a general tendency towards
positive bias is visible. This observation is an immediate consequence of the
construction in (10) since it is likely to observe extreme mean estimates
paired with large differences between the two mean processes; see column
$\rho=0.95$ in Figure 4. In summary, when moving to the RL scenario in the
following section, the principal behavior of the estimators remains. However,
one should keep in mind that the MEV estimators are constructed for the ideal
iid-case, and a violation of this assumption can impact the behavior of the
estimators when facing sequential data.

##  4 Application to Reinforcement Learning

###  4.1 Tabular Version

Reinforcement learning describes a collection of learning techniques for
sequential decision processes, in which an agent aims to maximize its reward
while interacting with an environment, see Sutton and Barto (2018) and
Bertsekas (2019). The problem is modeled as a Markov Decision Process (MDP,
Puterman 1994), consisting of a state space $\mathcal{S}$, a finite action
space $\mathcal{A}$, a state transition probability distribution
$\mathcal{P}:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow[0,1]$, a
bounded reward function
$\mathcal{R}:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$, and a
discount factor $\gamma\in[0,1]$. If $\gamma=1$, we assume there is a zero-
reward absorbing state and that the probability of reaching this state
converges to one as time tends to infinity, see Lan et al. (2020). At each
time step $t$, the agent takes an action $a_{t}\in\mathcal{A}$ based on state
information $s_{t}\in\mathcal{S}$, receives a reward
$r_{t}=\mathcal{R}(s_{t},a_{t})$, and transitions with probability
$\mathcal{P}(s_{t+1}\mid s_{t},a_{t})$ to a new state $s_{t+1}\in\mathcal{S}$.
Objective is to optimize for a policy
$\pi:\mathcal{S}\times\mathcal{A}\rightarrow[0,1]$, a mapping from states to
distributions over actions, that maximizes the expected return
$\operatorname{E}_{\pi}\left[\sum_{t=0}^{\infty}\gamma^{t}r_{t}\right]$.
Common practice are value-based methods that define action-values
$Q^{\pi}(s,a)=\operatorname{E}_{\pi}\left[\sum_{t=0}^{\infty}\gamma^{t}r_{t}\mid
s_{0}=s,a_{0}=a\right]$ for a certain policy. Thus, $Q^{\pi}(s,a)$ is the
expected return when starting in state $s$, executing $a$, and following
policy $\pi$ afterwards. There exists an optimal deterministic stationary
policy $\pi^{*}(s)=\operatorname*{arg\,max}_{a\in\mathcal{A}}Q^{*}(s,a)$, that
is connected with optimal action-values $Q^{*}(s,a)=\max_{\pi}Q^{\pi}(s,a)$
for all $s\in\mathcal{S}$ and $a\in\mathcal{A}$ if the state space is finite
or countably infinite (Puterman, 1994, Theorem 6.2.10). To optimize for
$Q^{*}(s,a)$, one uses a recursive relationship known as Bellman (1954)
optimality equation:

$Q^{*}(s,a)=\mathcal{R}(s,a)+\gamma\sum_{s^{\prime}\in\mathcal{S}}\mathcal{P}(s^{\prime}\mid
s,a)\max_{a^{\prime}\in\mathcal{A}}Q^{*}(s^{\prime},a^{\prime}),$ (12)

where $s^{\prime}$ is a successor state after performing action $a$ in state
$s$. Since $Q^{*}(s^{\prime},a^{\prime})$ is the expected return from
executing $a^{\prime}$ in $s^{\prime}$ and following the optimal policy
afterward, the problem immediately appears as an instance of the estimation of
the MEV of a set of random variables, namely the stochastic returns.
Consequently, the methodology of Section 2 applies. $Q$-Learning, from Watkins
and Dayan (1992), translates (12) into a sample-based algorithm by using the
ME:

$\hat{Q}^{*}(s,a)\leftarrow\hat{Q}^{*}(s,a)+\tau\left[y^{Q}-\hat{Q}^{*}(s,a)\right],$

with target
$y^{Q}=r+\gamma\max_{a^{\prime}\in\mathcal{A}}\hat{Q}^{*}(s^{\prime},a^{\prime})$
and learning rate $\tau$ (recall Section 3.4). An estimate of $Q^{*}(s,a)$ is
denoted $\hat{Q}^{*}(s,a)$. The algorithm is known to converge to the optimal
action-values if the conditions of Robbins and Monro (1951) on the learning
rate are fulfilled, and each state-action pair is visited infinitely often
(Tsitsiklis, 1994). However, especially in early stages of the training when
the $Q$-estimates are very imprecise, the algorithm tends to transmit
overestimated values. Overcoming this issue, Van Hasselt (2010) uses the DE
and stores two separate $Q$-tables with estimates $\hat{Q}^{*}_{A}$ and
$\hat{Q}^{*}_{B}$, leading to the target:
$y^{DQ}=r+\gamma\hat{Q}^{*}_{B}\left[s^{\prime},\operatorname*{arg\,max}_{a^{\prime}\in\mathcal{A}}\hat{Q}^{*}_{A}(s^{\prime},a^{\prime})\right]$.
Finally, to apply the TE and KE in a $Q$-Learning setup, we propose to replace
the target by:

$y^{KQ}=r+\gamma\left\\{\sum_{a^{\prime}\in\mathcal{A}}\kappa\left[T_{\hat{Q}^{*}}(s^{\prime},a^{\prime})\right]\right\\}^{-1}\sum_{a^{\prime}\in\mathcal{A}}\kappa\left[T_{\hat{Q}^{*}}(s^{\prime},a^{\prime})\right]\hat{Q}^{*}(s^{\prime},a^{\prime}),$
(13)

where

$T_{\hat{Q}^{*}}(s^{\prime},a^{\prime})=\frac{\hat{Q}^{*}(s^{\prime},a^{\prime})-\max_{a^{\prime\prime}\in\mathcal{A}}\hat{Q}^{*}(s^{\prime},a^{\prime\prime})}{\sqrt{\widehat{\operatorname{Var}}\left[\hat{Q}^{*}(s^{\prime},a^{\prime})\right]+\widehat{\operatorname{Var}}\left[\hat{Q}^{*}(s^{\prime},a^{*})\right]}},$

for a currently maximizing action
$a^{*}\in\\{a\in\mathcal{A}\mid\hat{Q}^{*}(s^{\prime},a)=\max_{a^{\prime\prime}\in\mathcal{A}}\hat{Q}^{*}(s^{\prime},a^{\prime\prime})\\}$.
The variance estimate of $\hat{Q}^{*}(s^{\prime},a^{*})$ is denoted
$\widehat{\operatorname{Var}}\left[\hat{Q}^{*}(s^{\prime},a^{*})\right]$ and
will be generated following the proposal of D’Eramo et al. (2019). First, the
variance of the underlying return is estimated via an exponentially-weighted
online update:

$\widehat{\sigma^{2}}_{\rm
process}(s,a)\leftarrow(1-\tau)\left\\{\widehat{\sigma^{2}}_{\rm
process}(s,a)+\tau\left[y^{KQ}-\hat{Q}^{*}(s,a)\right]^{2}\right\\}.$

Second, a normalization by Kish (1965) effective sample size $n_{\rm
eff}(s,a)$ is performed:

$\widehat{\operatorname{Var}}\left[\hat{Q}^{*}(s,a)\right]=\frac{\widehat{\sigma^{2}}_{\rm
process}(s,a)}{n_{\rm eff}(s,a)}.$

The effective sample size weights each sample depending on the learning rate
and is computed via $n_{\rm
eff}(s,a)=\frac{\omega(s,a)^{2}}{\omega^{2}(s,a)}$, in which numerator and
denominator are incrementally updated:

$\displaystyle\omega(s,a)$ $\displaystyle\leftarrow(1-\tau)\omega(s,a)+\tau,$
$\displaystyle\omega^{2}(s,a)$
$\displaystyle\leftarrow(1-\tau)^{2}\omega^{2}(s,a)+\tau^{2}.$

With this approach, we can introduce TE-$Q$-Learning (TE-$Q$) and
KE-$Q$-Learning (KE-$Q$), respectively, being summarized in Algorithm 1.

initialize $\forall s\in\mathcal{S},a\in\mathcal{A}:$
$\hat{Q}^{*}(s,a)=0,\widehat{\sigma^{2}}_{\rm
process}(s,a)>0,w(s,a)=0,w^{2}(s,a)=0$

repeat

Initialize $s$

repeat

Choose action $a$ from state $s$ with policy derived from $Q$ (e.g.,
$\epsilon$-greedy)

Take action $a$, observe reward $r$ and next state $s^{\prime}$

[0.75ex] _Update effective sample size:_

$\omega(s,a)\leftarrow(1-\tau)\omega(s,a)+\tau$

$\omega^{2}(s,a)\leftarrow(1-\tau)^{2}\omega^{2}(s,a)+\tau^{2}$

$n_{\rm eff}(s,a)\leftarrow\frac{\omega(s,a)^{2}}{\omega^{2}(s,a)}$

[0.75ex] Calculate target $y_{KQ}$ via (13)

[0.75ex] _Update process variance:_

$\widehat{\sigma^{2}}_{\rm
process}(s,a)\leftarrow(1-\tau)\left\\{\widehat{\sigma^{2}}_{\rm
process}(s,a)+\tau\left[y^{KQ}-\hat{Q}^{*}(s,a)\right]^{2}\right\\}$

[0.75ex] _Update $Q$-estimate:_

$\hat{Q}^{*}(s,a)\leftarrow\hat{Q}^{*}(s,a)+\tau\left[y^{KQ}-\hat{Q}^{*}(s,a)\right]$

$s\leftarrow s^{\prime}$

until  _$s$ is terminal_

until  __

Algorithm 1 TE-$Q$-Learning/KE-$Q$-Learning

###  4.2 Deep Version

Table-based algorithms cannot deal with continuous state spaces, and function
approximators like DNNs are used to parametrize the function
$\hat{Q}^{*}(s,a;\theta)$, with $\theta$ being the parameter set of the neural
network. The resulting DQN (Mnih et al., 2015) and their extensions (Van
Hasselt et al., 2016; Schaul et al., 2016; Hessel et al., 2018) have shown
great performances on a variety of challenging tasks. The optimization
procedure of these algorithms is still based on the Bellman optimality
equation (12), but uses gradient descend to update $\theta$:

$\theta\leftarrow\theta+\tau\left[y^{DQN}-\hat{Q}^{*}(s,a;\theta)\right]\nabla_{\theta}\hat{Q}^{*}(s,a;\theta),$

where
$y^{DQN}=r+\gamma\max_{a^{\prime}\in\mathcal{A}}\hat{Q}^{*}(s^{\prime},a^{\prime};\theta^{-})$.
The set $\theta^{-}$ refers to the parameters of the target network, which is
a time-delayed copy of the main network with parameter $\theta$. Instead of
updating fully online, DQN samples minibatches of past experiences from a
replay buffer $D$ to stabilize training. Van Hasselt et al. (2016) proposed
the Double DQN (DDQN) and uses the DE to compute the target:
$y^{DDQN}=r+\gamma\hat{Q}^{*}[s^{\prime},\operatorname*{arg\,max}_{a^{\prime}\in\mathcal{A}}\hat{Q}^{*}(s^{\prime},a^{\prime};\theta);\theta^{-}]$.
The action selection is performed via the main network, while the evaluation
uses the target network like the regular DQN.

To translate the TE and KE to the function approximation case, we require a
variance estimate of the $Q$-estimates. We follow D’Eramo et al. (2019) and
use the framework of the BDQN (Osband et al., 2016) to accomplish this task.
Generally, the bootstrap is a method for computing measures of accuracy for
statistical estimates (Efron and Tibshirani, 1994). The method trains
different regressors of the target function based on bootstrap samples
generated by sampling with replacement from the original dataset. The BDQN
transfers this idea to the DQN algorithm by maintaining $K\in\mathbb{N}$
differently initialized $Q$-networks $\hat{Q}^{*}_{k}(s,a;\theta_{k})$ with
parameter $\theta_{k}$, $k=1,\ldots,K$, each equipped with its own target
network $\hat{Q}^{*}_{k}(s,a;\theta_{k}^{-})$. This DQN modification was
proposed to improve over the usual $\epsilon$-greedy exploration strategy. At
the beginning of each episode, one $\hat{Q}^{*}_{k}$ is selected randomly, and
the agent acts greedy with relation to this $\hat{Q}^{*}_{k}$. At test time,
the majority vote of the function approximators is used. Generally, the BDQN
can be implemented using $K$ different networks or maintaining a common
network body and specifying $K$ different heads. We pursue the latter
approach. Diversification across the heads is achieved by different random
initialization of the parameters and random generation of binary masks
$m^{1},\ldots,m^{K}\in[0,1]$, to indicate which head should be trained on
which experience sample. A Bernoulli distribution with parameter $p$ is a
possible choice for the generating distribution $M$. Crucially, through the
$K$ heads, we can directly use the sample variance of the $Q$-estimates and
transfer (13) to compute target $y^{KDQN,k}$ for the $k$-th network
$\hat{Q}^{*}_{k}$:

$y^{KDQN,k}=r+\gamma\left\\{\sum_{a^{\prime}\in\mathcal{A}}\kappa\left[T_{\hat{Q}^{*}_{k}}(s^{\prime},a^{\prime})\right]\right\\}^{-1}\sum_{a^{\prime}\in\mathcal{A}}\kappa\left[T_{\hat{Q}^{*}_{k}}(s^{\prime},a^{\prime})\right]\hat{Q}^{*}_{k}(s^{\prime},a^{\prime};\theta_{k}^{-}),$
(14)

where

$T_{\hat{Q}^{*}_{k}}(s^{\prime},a^{\prime})=\frac{\hat{Q}^{*}_{k}(s^{\prime},a^{\prime};\theta^{-})-\max_{a^{\prime\prime}\in\mathcal{A}}\hat{Q}^{*}_{k}(s^{\prime},a^{\prime\prime};\theta_{k}^{-})}{\sqrt{\widehat{\operatorname{Var}}\left[\hat{Q}^{*}_{k}(s^{\prime},a^{\prime};\theta_{k}^{-})\right]+\widehat{\operatorname{Var}}\left[\hat{Q}^{*}_{k}(s^{\prime},a^{*};\theta_{k}^{-})\right]}},$

for a currently maximizing action
$a^{*}\in\\{a\in\mathcal{A}\mid\hat{Q}^{*}_{k}(s^{\prime},a;\theta_{k}^{-})=\max_{a^{\prime\prime}\in\mathcal{A}}\hat{Q}^{*}_{k}(s^{\prime},a^{\prime\prime};\theta_{k}^{-})\\}$.
The resulting gradient $g_{i}^{k}$ for the $i$-th tuple from the replay buffer
is:

$g_{i}^{k}=m_{i}^{k}\left[y_{i}^{KDQN,k}-\hat{Q}^{*}_{k}(s_{i},a_{i};\theta_{k})\right]\nabla_{\theta}\hat{Q}^{*}_{k}(s_{i},a_{i};\theta_{k}).$
(15)

The full procedure is termed TE-BDQN and KE-BDQN, respectively, being detailed
in Algorithm 2.

initialize Action-value estimate networks with $K$ outputs
$\left\\{\hat{Q}^{*}_{k}\right\\}^{K}_{k=1}$, masking distribution $M$, empty
replay buffer $D$

repeat

Initialize $s$

Pick a value function to act: $k\sim\text{Uniform}\\{1,\ldots,K\\}$

repeat

Choose action $a$ from state $s$ with greedy policy derived from
$\hat{Q}^{*}_{k}$

Take action $a$, observe reward $r$ and next state $s^{\prime}$

Sample bootstrap mask $m\sim M$

Add $(s,a,r,s^{\prime},m)$ to replay buffer $D$

Sample random minibatch of transitions
$\left\\{(s_{i},a_{i},s^{\prime}_{i},r_{i},m_{i})\right\\}_{i=1}^{B}$ from $D$

Perform gradient descent step based on (15)

Every $C$ steps reset $\theta_{k}=\theta_{k}^{-}$ for $k=1,\ldots,K$

$s\leftarrow s^{\prime}$

until  _$s$ is terminal_

until  __

Algorithm 2 TE-BDQN/KE-BDQN

##  5 Adaptive Absolute Bias Minimization

While the TE can interpolate between under- and overestimation by selecting a
smaller or larger $\alpha$, it is a priori not known which $\alpha$ is
adequate for an unknown environment. The method of choice for practitioners in
such situations is a grid search to select a suitable value empirically. Next
to the selection issue, a fixed parameter might not be sufficient for
controlling the estimation bias, and, e.g., an ascending strategy would be
favourable. Leveraging these considerations, we propose an adaptive
modification of TE-BDQN to adjust $\alpha$ under the objective of minimizing
the absolute estimation bias during training.

###  5.1 Bias Estimation

Before introducing the adaptive mechanism, we first outline how to estimate
the bias of given action-value estimates at a certain point in training.
Following Chen et al. (2021a), we consider a current policy $\pi$, which is
connected with true action-values $Q^{\pi}(s,a)$. The aggregated bias of
estimates $\hat{Q}^{\pi}(s,a)$ of $Q^{\pi}(s,a)$ for all
$s\in\mathcal{S},a\in\mathcal{A}$ is defined as:

$\operatorname{Bias}(\hat{Q}^{\pi},\pi)=\operatorname{E}_{s\sim\rho^{\pi},a\sim\pi}[\hat{Q}^{\pi}(s,a)-Q^{\pi}(s,a)],$

where $\rho^{\pi}$ is the state-visitation distribution of $\pi$. Chen et al.
(2021a) proposed to repeatedly run analysis episodes from random initial
states while following $\pi$. The observed Monte Carlo return of an
encountered state-action pair serves as an unbiased estimate of its true
$Q$-value. Averaging over all encountered state-action pairs yields the
estimate of the estimation bias:

$\widehat{\operatorname{Bias}}(\hat{Q}^{\pi},\pi)=\frac{1}{\lvert\mathcal{T}\rvert}\sum_{(s,a,R)\in\mathcal{T}}[\hat{Q}^{\pi}(s,a)-R],$
(16)

where $\mathcal{T}$ is a set of encountered $(s,a,R)$-tuples, where $s$ is the
state, $a$ the executed action, and $R$ the (later) observed Monte Carlo
return. $\lvert\mathcal{T}\rvert$ is the cardinality of $\mathcal{T}$. Chen et
al. (2021a) applied this procedure to the Soft-Actor Critic algorithm
(Haarnoja et al., 2018), which uses an actor dictating the policy $\pi$ and a
critic providing the estimates $\hat{Q}^{\pi}$. In standard $Q$-Learning, no
explicit actor is providing the policy, and the algorithm directly
approximates the optimal action-values $Q^{*}(s,a)$, leading to estimates
$\hat{Q}^{*}(s,a)$ (Sutton and Barto, 2018). To still generate insights into
the action-value estimation accuracy of the algorithm, Van Hasselt et al.
(2016) compare the $\hat{Q}^{*}(s,a)$ generated _during_ training with Monte
Carlo returns of the final greedy policy _after_ training. Although this
approach is certainly instructive, it does not enable an assessment without
having a converged baseline.

We instead propose to use (16) with
$\pi(s)=\operatorname*{arg\,max}_{a}\hat{Q}^{*}(s,a)$ for all
$s\in\mathcal{S}$ already _during_ training, being briefly summarized in
Algorithm 3. The method is justified since $Q$-Learning uses a greedy target
policy in its update. Thus, the algorithm _evaluates_ the greedy policy with
respect to its own action-value estimates, and we can assess whether the
$Q$-estimates are too optimistic or pessimistic. Transferred to the BDQN and
its modifications, the procedure can be similarly applied by assessing each
head separately since they constitute different estimates $\hat{Q}^{*}_{k}$
for $k=1,\ldots,K$. Thus, we run Algorithm 3 for each head and average the
output to receive an aggregated bias estimate for the bootstrap-ensemble.

Input Estimates $\hat{Q}^{*}(s,a)$ for all $s\in\mathcal{S},a\in\mathcal{A}$

Set $\mathcal{T}=\emptyset$ and
$\pi(s)=\operatorname*{arg\,max}_{a}\hat{Q}^{*}(s,a)$ for all
$s\in\mathcal{S}$

for  _number of episodes_ do

Randomly initialize $s$

Play episode following $\pi$ and append encountered state-action-return tuples
to $\mathcal{T}$

end for

return
_$\frac{1}{\lvert\mathcal{T}\rvert}\sum_{(s,a,R)\in\mathcal{T}}[\hat{Q}^{*}(s,a)-R]$_

Algorithm 3 Bias estimation for $Q$-Learning like algorithms

###  5.2 Adaptive TE-BDQN

Algorithm 3 will be used to assess the algorithms in the experiments of
Section 6. Furthermore, the approach serves as a basis for dynamically
adjusting the $\alpha$ of the TE-BDQN. Intuitively, since larger $\alpha$
leads to larger $Q$-estimates, $\alpha$ is reduced if the $Q$-estimates are
too high. Vice versa, we increase $\alpha$ if the $Q$-estimates are too small.
Precisely, we perform the following update:

$\alpha\leftarrow\alpha+\frac{\tau_{\rm
Ada}}{K}\sum_{k=1}^{K}\sum_{t=1}^{T_{\rm
Ada}}\left[R_{k}(s_{t,k},a_{t,k})-\hat{Q}^{*}_{k}(s_{t,k},a_{t,k};\theta_{k})\right],$
(17)

with step size $\tau_{\rm Ada}$ and roll-out length $T_{\rm Ada}$.
Importantly, we use $n$-steps returns (Sutton and Barto, 2018):

$R_{k}(s_{t,k},a_{t,k})=r_{t,k}+\gamma
r_{t+1,k}+\gamma^{2}r_{t+2,k}+\ldots+\gamma^{T_{\rm Ada}-t}r_{T_{\rm
Ada},k}+\gamma^{T_{\rm Ada}-t+1}\max_{a^{\prime}}\hat{Q}^{*}_{k}(s_{T_{\rm
Ada}+1,k},a^{\prime};\theta_{k}),$

where $t=1,\ldots,T_{\rm Ada}$. We denote state and action at time $t$ under
head $k$ as $s_{t,k}$ and $a_{t,k}$, respectively. The resulting immediate
reward is $r_{t,k}$ and the initial states $s_{1,k}$ for $k=1,\ldots,K$ are
randomly sampled from the replay buffer. As motivated in Section 5.1, the
actions $a_{t,k}$ under head $k$ are selected by acting greedy with relation
to $\hat{Q}^{*}_{k}$. Although $n$-step returns are generally not an unbiased
estimate of the expected return of a policy like a complete Monte Carlo roll-
out, we found them empirically much more practicable since they do not require
running full episodes while still allowing us to judge the accuracy of the
current value estimates. Consequently, this approach can also be applied to
non-episodic problems.

We update $\alpha$ immediately after the target networks, avoiding
instabilities in the learning process. A similar proposal to (17) in an
episodic context with continuous action spaces based on full Monte Carlo roll-
outs has recently been made by Dorka et al. (2021). Note that it is possible
to maintain a separate $\alpha$ for each bootstrap head, enabling a tailored
parametrization for the bias of each approximator. However, we only consider
one parameter for the whole ensemble for simplicity in the following. The
resulting algorithm is called Ada-TE-BDQN and is shown in Appendix B.

##  6 Experiments

We analyze the proposed estimators of the MEV on a statistically motivated
real-world example before considering two tabular environments that serve as a
proof-of-concept for TE/KE-$Q$-Learning. The experiments with function
approximation are carried out in the MinAtar environments of Young and Tian
(2019), which allow for a thorough algorithmic comparison.

###  6.1 Internet Ads

We consider the internet ad problem previously studied by Van Hasselt (2013),
D’Eramo et al. (2021), and Jiang et al. (2021). There are $M$ different ads,
and each has the same return per click. Consequently, the click rate is the
only quantity of interest and is modeled as a Bernoulli variable (click or no
click). The true expectations $\mu_{1},\ldots,\mu_{M}$ of these $M$ variables
equal the respective click probability and are equally spaced in a specific
interval $\mu_{int}$. We consider $N$ customers and assume that the ads are
presented equally often to have $N/M$ samples per ad. The objective is to
estimate the maximum true mean accurately, thus finding the best ad based on
the given samples. We compare the TE ($\alpha=0.1$) and the KE (standard
Gaussian kernel) with the ME, DE, and WE based on bias, variance, and MSE. Six
configurations of the problem are considered by varying the number of
customers $N$, the number of ads $M$, or the upper limit of the sampling
interval $\mu_{\rm int}$, while the lower limit is always fixed at 0.02.
Figure 6 displays the results.

Figure 6: Comparison of ME, DE, WE, TE, and KE on the internet ad problem.
Results are averaged over 10 000 runs.

Both TE and KE yield lower MSE than their competitors in most scenarios. We
emphasize that $\alpha=0.1$ was not cherry-picked for this problem, and the
TE’s performance could thus be further increased by tailoring $\alpha$ for
each experiment. In general, all estimators’ MSEs decrease with an increasing
number of ads, while they are more accurate with a higher number of customers
$N$, constituting reasonable observations. The DE often yields large
variances, while the ME provides biased estimates. Despite producing a higher
MSE than TE and KE in most cases, the WE outperforms the conventional
competitors ME and DE.

###  6.2 Maximization Bias Example

We consider the example in Figure 6.5 of Sutton and Barto (2018). A simple MDP
with two non-terminal states A and B is given. The agent starts in A. If it
goes ’right’ from A, the episode ends, and zero reward is received. If action
’left’ is selected in A, the agent deterministically gets to state B and
receives zero reward. There are eight actions to choose from B, but all lead
to a terminal state and yield a reward sampled from $\mathcal{N}(-0.1,1)$. The
parameters are $\gamma=1$, $\epsilon=0.1$, and $\tau=0.1$. Since this is an
undiscounted task, the expected return starting with action ’left’ is $-0.1$,
and the agent should always prefer going ’right.’ However, due to the random
selection of the $\epsilon$-greedy strategy, the ’left’ action will always be
picked at least $5\%$ in expectation.

Figure 7 depicts the results. The upper-left part displays the percentage of
selecting action ’left’ in A; the upper-right plot contains the same
percentage after 500 training episodes. The lower-right graph shows the
estimate of $Q^{*}$(A, left) over training. Finally, the lower-left plot
displays the number of averaged means of TE-$Q$-Learning when updating
$\hat{Q}^{*}$(A, left). $Q$-Learning is the same algorithm as TE-$Q$-Learning
with $\alpha=0.5$ and only included for comparison.

$Q$-Learning initially overestimates the value of the ’left’ action in state A
and selects it nearly twice as optimal after 500 episodes. Double $Q$-Learning
performs better and achieves a final rate of roughly $6\%$. On the other hand,
TE-$Q$-Learning can modulate the overestimation bias through its significance
level $\alpha$ and reaches a near-optimal selection percentage for $\alpha\leq
0.10$. Interestingly, TE-$Q$ ($\alpha=0.4$) performs worse after 500 episodes
than $Q$-Learning. Although the initial overestimation is not as large as for
$\alpha=0.5$, the effect is more persistent, and more interactions are needed
to reduce the estimate. For additional insights, we display the number of
means averaged by the TE. TE-$Q$-Learning with $\alpha=0.5$ naturally
considers only the maximum sample mean, which is non-unique in the first
several episodes since all action-value estimates are initialized with zero.
The lower the significance level of TE-$Q$ gets, the more means are averaged
until nearly all sample means are considered for $\alpha\leq 0.1$.
Furthermore, KE-$Q$-Learning with the standard Gaussian kernel performs
reasonably well and achieves a final selection rate below Double $Q$-Learning.

Figure 7: Maximization Bias example from Sutton and Barto (2018) with
parameters $\gamma=1$, $\epsilon=0.1$, and $\tau=0.1$. Q and DQ refer to
$Q$-Learning and Double $Q$-Learning, respectively. Results are averaged over
$100\,000$ runs and 95% confidence intervals are included. The action-value
estimate in the lower-right graph for DQ is generated by averaging over both
$Q$-tables.

###  6.3 Cliff Walking

We examine the Cliff Walking task from Example 6.6 in Sutton and Barto (2018),
which is an undiscounted, episodic task with start and goal states. Our
environment is a grid of width 10 and height 5. Start state S is the lower-
left grid point; goal state G is the lower-right grid point. All transitions
are rewarded with $-1$, except those which lead to the grid points directly
between S and G. Those are referred to as ’Cliff,’ yield reward $-100$, and
send the agent back to S. Actions are the four movement directions up, down,
right, and left. Performance is measured via the return during an episode.

Figure 8 follows the setup of Zhu and Rigotti (2021) and contains results for
constant $\epsilon=0.1$ and annealing $\epsilon=1/\sqrt{n(s)}$, where $n(s)$
is the number of times state $s$ has been visited. The learning rate is
$\tau=0.1(100+1)/\left[100+n(s,a)\right]$, with $n(s,a)$ being the number of
updates for the state-action pair. Next to $Q$\- and Double $Q$-Learning, we
consider Weighted $Q$-Learning (WQ, D’Eramo et al. 2016) and Self-correcting
$Q$-Learning (SCQ, Zhu and Rigotti 2021) with $\beta\in\\{2,3\\}$, following
the recommendation of the authors. We run each algorithm for 3000 episodes and
average the results over 500 independent runs. Additionally, we display the
maximum action-value estimate of the start state S over training. For
comparison, since at least eleven steps are necessary to walk across our map,
it holds for the optimal policy: $\max_{a^{\prime}}Q^{*}(S,a^{\prime})=-11$.

We see the strong performance of the newly proposed algorithms for both
exploration strategies. Similar to the example in Section 6.2, especially
TE-$Q$ with $\alpha=0.05$ and KE-$Q$ are appropriate for this task and achieve
the highest returns together with WQ. Furthermore, the higher action-value
estimates for $Q$-Learning are apparent, while Double $Q$-Learning leads to
severe underestimation. Finally, the returns are higher for all algorithms
with $\epsilon=1/\sqrt{n(s)}$ than with a constant exploration rate, which is
reasonable since action selection yields a higher probability of selecting
greedy in the long-term.

Figure 8: Cliff Walking example from Sutton and Barto (2018) with parameters
$\gamma=1$, $\tau=0.1(100+1)/\left[100+n(s,a)\right]$, and two different
$\epsilon$-greedy strategies. Results are averaged over 500 runs,
exponentially smoothed for visualization purposes, and 95% confidence
intervals are included. The maximum action-value estimate of the start state
for DQ is computed by averaging this quantity over both $Q$-tables.

To investigate the interim and asymptotic behavior of the algorithms for
different learning rates, we run an analysis similar to Van Seijen et al.
(2009). Precisely, we consider learning rates in
$\left\\{0.1,0.2,\ldots,0.9,1.0\right\\}$ and employ again the two exploration
strategies $\epsilon=0.1$ and $\epsilon=1/\sqrt{n(s)}$. We analyze the average
return over the first 100 episodes and average the results over 5000 runs for
the interim performance. For the asymptotic scenario, we run each algorithm
for $50\,000$ episodes and average the results over 5 runs.

Figure 9: Cliff Walking example adapted from Van Seijen et al. (2009). The
algorithms’ interim (dotted lines) and asymptotic (solid lines) return
averages are analyzed for different learning rates. The number of episodes is
$n$.

The TE-$Q$ and KE-$Q$ algorithms offer the most robust interim progress across
learning rates for both exploration strategies, while the DQ and SCQ expose a
severe performance drop when $\tau=1$. This might be because both algorithms
rely on two different $Q$-tables, and complete replacement of the entries
yields instabilities in this case. Regarding the asymptotic analysis, the SCQ
and DQ algorithms improve on $Q$-Learning and are marginally above WQ, TE-$Q$,
and KE-$Q$ for $\epsilon=0.1$, while the return differences are close to zero
for the annealing exploration strategy due to long-term greedy action
selection.

###  6.4 MinAtar

We select the MinAtar (Young and Tian, 2019) environments to test the proposed
Deep RL algorithms. MinAtar is a testbed incorporating several Atari games
from the Arcade Learning Environment (Bellemare et al., 2013) with a reduced
state-representation. The platform incorporates sticky actions (Machado et
al., 2018) and is designed to enable thorough algorithmic comparisons due to
reduced computation times. Following Young and Tian (2019), the network
structure consists of a convolutional and a fully-connected layer. The
remaining hyperparameters match Young and Tian (2019), except that we use the
Adam (Kingma and Ba, 2014) optimizer, which led to much more stable results
during our experiments. Appendix C contains the full list of specifications.

The compared algorithms are the DQN, DDQN, Self-Correcting DQN (SCDQN, Zhu and
Rigotti 2021), BDQN, TE-BDQN, KE-BDQN (with standard Gaussian cdf), and Ada-
TE-BDQN. For the parametrization of the BDQN and its modifications, we follow
Osband et al. (2016) by using $K=10$ bootstrap heads, each corresponding to
one fully-connected layer, and setting $p=1$ for the masking Bernoulli
distribution. The BDQN uses the target computation of the DDQN, which we apply
consequently. Furthermore, we scale the gradients of the convolutional core
part for the bootstrap-based algorithms by $1/K$, which was also recommended
by Osband et al. (2016). We consider $\beta\in\\{2,3,4\\}$ for the SCDQN and
$\alpha\in\\{0.1,0.2,0.3,0.4\\}$ for the TE-BDQN. The bias parameter of the
Ada-TE-BDQN is initialized with $\alpha=0.25$ and we consider two step sizes
$\tau_{\rm Ada}\in\\{10^{-4},10^{-5}\\}$ with horizon $T_{\rm Ada}=32$.

Figure 10: Algorithm comparison on Asterix. The top three rows show the
training for $\tau=10^{-4.5}$, while the plots in the last row compare the
final return over different learning rates. DQN, DDQN, and SCDQN are in the
left column; the BDQN-based algorithms are in the right column. Figure 11:
Algorithm comparison on Breakout. The top three rows show the training for
$\tau=10^{-4.5}$, while the plots in the last row compare the final return
over different learning rates. DQN, DDQN, and SCDQN are in the left column;
the BDQN-based algorithms are in the right column. Figure 12: Algorithm
comparison on Freeway. The top three rows show the training for
$\tau=10^{-4.5}$, while the plots in the last row compare the final return
over different learning rates. DQN, DDQN, and SCDQN are in the left column;
the BDQN-based algorithms are in the right column. Figure 13: Algorithm
comparison on Seaquest. The top three rows show the training for
$\tau=10^{-4.5}$, while the plots in the last row compare the final return
over different learning rates. DQN, DDQN, and SCDQN are in the left column;
the BDQN-based algorithms are in the right column. Figure 14: Algorithm
comparison on SpaceInvaders. The top three rows show the training for
$\tau=10^{-4.5}$, while the plots in the last row compare the final return
over different learning rates. DQN, DDQN, and SCDQN are in the left column;
the BDQN-based algorithms are in the right column.

To check the robustness of the algorithms, we analyze three different learning
rates for each environment and algorithm:
$\tau\in\\{10^{-5},10^{-4.5},10^{-4}\\}$. Every $10\,000$ steps during an
experiment, we average the return of 10 test episodes. For the BDQN and its
variants, the majority vote of the ensemble is applied. Additionally, we run
for all algorithms bias estimation episodes from random initial states sampled
from the replay buffer, following Algorithm 3. The number of those episodes
are 10 for the DQN, DDQN, and SCDQN, while we run only 3 episodes for each
head of the BDQN-based algorithms due to computation time. We repeat all
experiments for ten independent runs, exponentially smooth the results for
clarity, and include $95\%$ point-wise confidence intervals over the runs.

Figures 10 \- 14 depict the results. We show the bias and return training
plots for $\tau=10^{-4.5}$ in the upper two rows of each figure, while the
last row contains the final return across learning rates. The BDQN-based
algorithms generally outperform their competitors, and especially the KE-BDQN
and Ada-TE-BDQN show a robust performance across environments, although the
algorithms’ variances are relatively high in Seaquest. As expected, the DQN is
affected by massive overestimations, while the DDQN can reduce the
$Q$-estimates in comparison. Although the DE theoretically underestimates the
MEV, the DDQN still offers a positive bias in the given set of experiments.
This observation is in line with Van Hasselt et al. (2016) and illustrates the
effect of time-dependence and function approximation on the analysis, compare
Section 3.4.

As theoretically discussed in prior sections, a larger $\alpha$ in the TE-BDQN
yields larger $Q$-estimates and, consequently, a larger estimation bias. The
adaptive mechanism of Ada-TE-BDQN, especially for $\tau_{\rm Ada}=10^{-4}$,
results in approximately unbiased action-value estimates. Throughout
environments, the $\alpha$ of the Ada-TE-BDQN mostly increases during training
but stays on moderate values in approximately $[0.2,0.4]$. Large values with
$\alpha\approx 0.5$ are not achieved even in later stages of training,
indicating the criticality of the ME combined with function approximation.

###  6.5 Discussion

The experimental results confirm that we can embed statistical thinking in
form of the TE/KE into value-based RL methods to control their estimation
bias. In addition, the experiments support the finding of Lan et al. (2020)
that unbiased $Q$-estimation does not necessarily translate into the best
return performance. For example, in SpaceInvaders, the KE-BDQN is return-wise
the strongest algorithm despite its severe negative bias. However, the TE-BDQN
($\alpha=0.1$) offers an even lower bias but cannot match the return of the
KE-BDQN. There seems to be a critical level - or path over time - of
estimation bias for a given MDP, which yields maximum performance. Careful
selection of a bias control parameter like $\alpha$ for the TE or the kernel
function for the KE thus constitutes a crucial component in designing
temporal-difference algorithms.

Analyzing the behavior of the algorithms in more detail, we see that the
estimated bias _changes over time_. With random initialization of the networks
and a couple of zero-return episodes, all algorithms’ bias is shortly
approximately zero. As soon as some non-zero rewards are observed, the
different target specifications affect the update routine and result in
severely different bias plots over time. Besides the Ada-TE-BDQN, each
algorithm reveals its tendency towards over- or underestimation, although
exceptions are possible. For example, the TE-BDQN ($\alpha=0.3$) offers slight
overestimations during the first three million steps in Breakout before
shifting towards underestimation. Importantly, none of the non-adaptive
algorithms shows reliable convergence to zero-bias as training proceeds, which
agrees with the observations of Van Hasselt et al. (2016). Finally, we
summarize the core findings of our investigation:

  1. 1.

_Absolute bias minimization does not equal return maximization._ In order to
maximize performance in a real application, different bias control
configurations should be considered, for which the TE/KE-BDQN build a flexible
framework.

  2. 2.

_Approximately unbiased estimation offers a robust baseline across tasks._
Although it is not always the return-maximizing choice, using a scheme for
approximately unbiased value-estimation appears more robust across tasks than
fixing a particular bias control parameter. The Ada-TE-BDQN is a powerful
candidate for such a scheme.

  3. 3.

_The compatibility between bias control algorithms and exploration schemes
requires systematic analysis._ The impact on exploration most likely
constitutes an essential factor in the occasional return-improvement of a
biased procedure over an unbiased one (Liang et al., 2021). Further study
needs to generate insights on how these components interact and, crucially,
whether the algorithmic approaches are compatible. Can we use the Ada-TE-BDQN
with a modified exploration scheme to boost return performance? Do we maintain
unbiased action-value estimation in this process? Can we achieve or even
improve on the return peaks of a fine-tuned bias control configuration in this
manner? Exploration and action-value estimation are even in off-policy RL not
necessarily orthogonal and thus constitute a crucial path for future research.

##  7 Related Works

Next to Van Hasselt (2010), D’Eramo et al. (2016), and Zhu and Rigotti (2021),
several proposals have been made to further tackle the issue of estimation
bias in temporal-difference algorithms. Zhang et al. (2017) proposed a hybrid
between the ME and the DE called Weighted Double Estimator. It relies on a
hyperparameter on the positive real axis, for which the authors propose a
heuristic based on empirical experiences. Lee et al. (2013) proposed Bias-
corrected $Q$-Learning, which incorporates a correction term depending on the
reward variance. Ensemble-based methods like Maxmin $Q$-Learning (Lan et al.,
2020) and Randomized Ensembled Double $Q$-Learning (REDQ, Chen et al. 2021a)
apply a minimization operator over the ensemble or a subset of the latter. The
Action-Candidate based Clipped Double Estimator (Jiang et al., 2021) extends
the DE by creating a so-called candidate set of indices of which the
maximizing one will be picked. Furthermore, the Clipped Double Estimator of
Fujimoto et al. (2018) and the Truncated Quantile Critic (TQC, Kuznetsov et
al. 2020) algorithm are relevant contributions to addressing the
overestimation issue in actor-critic frameworks. Finally, Lee et al. (2021)
pursuit a re-weighting strategy of sampled transitions based on uncertainty
estimates from an ensemble. Apart from methodological extensions, Chen et al.
(2021b) recently reported that a lower learning rate or an adequate schedule
could also avoid the massive overestimations of $Q$-Learning. However,
lowering the learning rate can come at the expense of impractically slow
learning, as seen in our Breakout experiments, and constitutes thus not a
practical option to address the issue of action-value estimation.

Recently, some proposals have been made to minimize the estimation bias of
temporal-difference algorithms through online parameter adjustments in the
spirit of the Ada-TE-BDQN. Liang et al. (2021) expand the work of Fox et al.
(2016), Fox (2019) by using an ensemble to adjust the temperature parameter in
a maximum entropy framework. Kuznetsov et al. (2021) and Dorka et al. (2021)
introduce adaptive variants of the TQC by adjusting the number of quantiles to
drop based on recent near on-policy trajectories. Finally, Wang et al. (2021)
generalize MaxMin $Q$-Learning and REDQ by changing the size of the subset of
the ensemble on which the minimization operator is performed. The metric
driving the adjustment is the ensemble’s function approximation error since it
is argued that high approximation error is connected with the overestimation
of action-values.

##  8 Conclusion

Reinforcement learning is a domain of artificial intelligence with significant
breakthroughs in a diverse set of real-world applications, particularly in the
last decade. A key issue of frequently applied temporal-difference algorithms
is the propagation of biased action-value estimates. We address this topic by
proposing the $T$-Estimator and the $K$-Estimator for the underlying problem
of estimating the maximum expected value of random variables. Both estimators
are easy to compute and allow to flexibly interpolate between over- and
underestimation bias, leading to promising modifications of $Q$-Learning and
the Bootstrapped DQN algorithm. Coupled with the dynamic selection procedure
of the significance level of TE, our work constitutes an important step
towards unbiased estimation of action-values with function approximation.

In future research, we will analyze the discussed interplay of action-value
estimation and exploration. As methodological extensions, we will investigate
possibilities to extend the two-sample testing procedures into continuous
actions spaces to modify policy gradient methods because the latter constitute
an elementary class of methods in several application domains. Furthermore,
next to the considered procedure, there are alternative approaches for
uncertainty quantification in the neural network scenario. For example, the
regularization technique dropout (Srivastava et al., 2014) will be applied
similar to D’Eramo et al. (2021) to obtain the required variance estimate for
the newly proposed algorithms, and the Bootstrapped DQN will be enhanced by
adding random prior functions (Osband et al., 2018). Finally, our work focuses
on the estimation bias over the whole state-action distribution of a policy,
and we aggregate this quantity into one scalar estimate. While we analyze how
this scalar changes throughout training, we do not differentiate how the bias
is distributed over the state-action space at a particular time in training.
We acknowledge that assessing complex MDPs in this fashion might result in an
over-simplification. In the future, we will work on more tailored solutions
and consider, e.g., an individual significance level of the TE for different
regions of the state-action space.

### Acknowledgements

We would like to thank Niklas Paulig for fruitful discussions in the early
stages of this work. Furthermore, we want to thank two anonymous reviewers for
their constructive feedback, which improved this manuscript thoroughly.

## References

  * Armstrong (2014) Armstrong, R. A. When to use the Bonferroni correction. Ophthalmic and Physiological Optics, 2014, 34, 502–508.
  * Aven (1985) Aven, T. Upper (lower) bounds on the mean of the maximum (minimum) of a number of random variables. Journal of Applied Probability, 1985, 22, 723–728.
  * Bellemare et al. (2013) Bellemare, M. G., Naddaf, Y., Veness, J., and Bowling, M. The arcade learning environment: An evaluation platform for general agents. Journal of Artificial Intelligence Research, 2013, 47, 253–279.
  * Bellman (1954) Bellman, R. The theory of dynamic programming. Bulletin of the American Mathematical Society, 1954, 60, 503–515.
  * Bertsekas (2019) Bertsekas, D., Reinforcement Learning and Optimal Control, 2019, Belmont: Athena Scientific.
  * Blumenthal and Cohen (1968) Blumenthal, S. and Cohen, A. Estimation of the larger of two normal means. Journal of the American Statistical Association, 1968, 63, 861–876.
  * Casella and Berger (2002) Casella, G. and Berger, R. L., Statistical Inference (2nd ed.), 2002, Belmont, CA: Brooks/Cole Cengage Learning.
  * Chen et al. (2021a) Chen, X., Wang, C., Zhou, Z., and Ross, K. W. 2021a, Randomized Ensembled Double Q-Learning: Learning Fast Without a Model. In International Conference on Learning Representations.
  * Chen et al. (2021b) Chen, Y., Schomaker, L., and Wiering, M. A. 2021b, An Investigation Into the Effect of the Learning Rate on Overestimation Bias of Connectionist Q-learning. In International Conference on Agents and Artificial Intelligence, pp. 107–118.
  * D’Eramo et al. (2021) D’Eramo, C., Cini, A., Nuara, A., Pirotta, M., Alippi, C., Peters, J., and Restelli, M. Gaussian Approximation for Bias Reduction in Q-Learning. Journal of Machine Learning Research, 2021, 22, 1–51.
  * Dhariyal et al. (1985) Dhariyal, I., Sharma, D., and Krishnamoorthy, K. Non-existence of unbiased estimators of ordered parameters. Statistics: A Journal of Theoretical and Applied Statistics, 1985, 16, 89–95.
  * Dorka et al. (2021) Dorka, N., Boedecker, J., and Burgard, W. 2021, Adaptively Calibrated Critic Estimates for Deep Reinforcement Learning. In Deep RL Workshop NeurIPS 2021.
  * Dudewicz (1971) Dudewicz, E. J. Maximum likelihood estimators for ranked means. Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, 1971, 19, 29–42.
  * D’Eramo et al. (2019) D’Eramo, C., Cini, A., and Restelli, M. 2019, Exploiting Action-Value uncertainty to drive exploration in reinforcement learning. In International Joint Conference on Neural Networks, IEEE, pp. 1–8.
  * D’Eramo et al. (2016) D’Eramo, C., Restelli, M., and Nuara, A. 2016, Estimating maximum expected value through gaussian approximation. In International Conference on Machine Learning, PMLR, pp. 1032–1040.
  * Efron and Tibshirani (1994) Efron, B. and Tibshirani, R. J., An Introduction to the Bootstrap, 1994, CRC press.
  * Fox (2019) Fox, R. 2019, Toward provably unbiased temporal-difference value estimation. In Optimization Foundations for Reinforcement Learning Workshop at NeurIPS.
  * Fox et al. (2016) Fox, R., Pakman, A., and Tishby, N. 2016, Taming the noise in reinforcement learning via soft updates. In Conference on Uncertainty in Artificial Intelligence, pp. 202–211.
  * Fujimoto et al. (2018) Fujimoto, S., Hoof, H., and Meger, D. 2018, Addressing function approximation error in actor-critic methods. In International Conference on Machine Learning, PMLR, pp. 1587–1596.
  * Haarnoja et al. (2018) Haarnoja, T., Zhou, A., Abbeel, P., and Levine, S. 2018, Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor. In International conference on machine learning, PMLR, pp. 1861–1870.
  * Hessel et al. (2018) Hessel, M., Modayil, J., Van Hasselt, H., Schaul, T., Ostrovski, G., Dabney, W., Horgan, D., Piot, B., Azar, M., and Silver, D. 2018, Rainbow: Combining improvements in deep reinforcement learning. In AAAI Conference on Artificial Intelligence, vol. 32.
  * Imagaw and Kaneko (2017) Imagaw, T. and Kaneko, T. 2017, Estimating the maximum expected value through upper confidence bound of likelihood. In Conference on Technologies and Applications of Artificial Intelligence, IEEE, pp. 202–207.
  * Jiang et al. (2021) Jiang, H., Xie, J., and Yang, J. 2021, Action Candidate Based Clipped Double Q-learning for Discrete and Continuous Action Tasks. In AAAI Conference on Artificial Intelligence, vol. 35, pp. 7979–7986.
  * Kingma and Ba (2014) Kingma, D. P. and Ba, J. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
  * Kish (1965) Kish, L., Survey Sampling, 1965, New York: John Wiley & Sons.
  * Kuznetsov et al. (2021) Kuznetsov, A., Grishin, A., Tsypin, A., Ashukha, A., and Vetrov, D. Automating Control of Overestimation Bias for Continuous Reinforcement Learning. arXiv preprint arXiv:2110.13523, 2021.
  * Kuznetsov et al. (2020) Kuznetsov, A., Shvechikov, P., Grishin, A., and Vetrov, D. 2020, Controlling overestimation bias with truncated mixture of continuous distributional quantile critics. In International Conference on Machine Learning, PMLR, pp. 5556–5566.
  * Lan et al. (2020) Lan, Q., Pan, Y., Fyshe, A., and White, M. 2020, Maxmin Q-learning: Controlling the Estimation Bias of Q-learning. In International Conference on Learning Representations.
  * Lee et al. (2013) Lee, D., Defourny, B., and Powell, W. B. 2013, Bias-corrected q-learning to control max-operator bias in q-learning. In Symposium on Adaptive Dynamic Programming and Reinforcement Learning, IEEE, pp. 93–99.
  * Lee et al. (2021) Lee, K., Laskin, M., Srinivas, A., and Abbeel, P. 2021, SUNRISE: A simple unified framework for ensemble learning in deep reinforcement learning. In International Conference on Machine Learning, PMLR, pp. 6131–6141.
  * Liang et al. (2021) Liang, L., Xu, Y., McAleer, S. M., Hu, D., Ihler, A., Abbeel, P., and Fox, R. 2021, Temporal-Difference Value Estimation via Uncertainty-Guided Soft Updates. In Deep RL Workshop NeurIPS 2021.
  * Machado et al. (2018) Machado, M. C., Bellemare, M. G., Talvitie, E., Veness, J., Hausknecht, M., and Bowling, M. Revisiting the arcade learning environment: Evaluation protocols and open problems for general agents. Journal of Artificial Intelligence Research, 2018, 61, 523–562.
  * Mammen (1991) Mammen, E. Estimating a smooth monotone regression function. The Annals of Statistics, 1991, 19, 724–740.
  * Mnih et al. (2015) Mnih, V., Kavukcuoglu, K., Silver, D., Rusu, A. A., Veness, J., Bellemare, M. G., Graves, A., Riedmiller, M., Fidjeland, A. K., Ostrovski, G., et al. Human-level control through deep reinforcement learning. Nature, 2015, 518, 529–533.
  * Nadarajah and Kotz (2008) Nadarajah, S. and Kotz, S. Exact distribution of the max/min of two Gaussian random variables. IEEE Transactions on very large scale integration systems, 2008, 16, 210–212.
  * Osband et al. (2018) Osband, I., Aslanides, J., and Cassirer, A. Randomized Prior Functions for Deep Reinforcement Learning. Advances in Neural Information Processing Systems, 2018, 31.
  * Osband et al. (2016) Osband, I., Blundell, C., Pritzel, A., and Van Roy, B. Deep exploration via bootstrapped DQN. Advances in Neural Information Processing Systems, 2016, 29, 4026–4034.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Köpf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S. Pytorch: An imperative style, high-performance deep learning library. Advances in Neural Information Processing Systems, 2019, 32, 8026–8037.
  * Puterman (1994) Puterman, M. L., Markov Decision Processes: Discrete Stochastic Dynamic Programming, 1994, John Wiley & Sons.
  * Robbins and Monro (1951) Robbins, H. and Monro, S. A stochastic approximation method. The Annals of Mathematical Statistics, 1951, 400–407.
  * Schaul et al. (2016) Schaul, T., Quan, J., Antonoglou, I., and Silver, D. Prioritized experience replay. International Conference on Learning Representations, 2016.
  * Sheather (2004) Sheather, S. J. Density estimation. Statistical Science, 2004, 19, 588–597.
  * Silver et al. (2017) Silver, D., Schrittwieser, J., Simonyan, K., Antonoglou, I., Huang, A., Guez, A., Hubert, T., Baker, L., Lai, M., Bolton, A., et al. Mastering the game of go without human knowledge. Nature, 2017, 550, 354–359.
  * Srivastava et al. (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R. Dropout: a simple way to prevent neural networks from overfitting. Journal of Machine Learning Research, 2014, 15, 1929–1958.
  * Sutton and Barto (2018) Sutton, R. S. and Barto, A. G., Reinforcement Learning: An Introduction, 2018, Cambridge: The MIT Press.
  * Thrun and Schwartz (1993) Thrun, S. and Schwartz, A. 1993, Issues in using function approximation for reinforcement learning. In Proceedings of the Fourth Connectionist Models Summer School, Hillsdale, NJ, pp. 255–263.
  * Tsay (2010) Tsay, R. S., Analysis of Financial Time Series, 2010, New Jersey: John Wiley & Sons.
  * Tsitsiklis (1994) Tsitsiklis, J. N. Asynchronous stochastic approximation and Q-learning. Machine learning, 1994, 16, 185–202.
  * Van Hasselt (2010) Van Hasselt, H. Double Q-learning. Advances in Neural Information Processing Systems, 2010, 23, 2613–2621.
  * Van Hasselt (2013) —Estimating the maximum expected value: an analysis of (nested) cross validation and the maximum sample average. arXiv preprint arXiv:1302.7175, 2013.
  * Van Hasselt et al. (2016) Van Hasselt, H., Guez, A., and Silver, D. 2016, Deep reinforcement learning with double Q-learning. In AAAI Conference on Artificial Intelligence, vol. 30.
  * Van Seijen et al. (2009) Van Seijen, H., Van Hasselt, H., Whiteson, S., and Wiering, M. 2009, A theoretical and empirical analysis of Expected Sarsa. In Symposium on Adaptive Dynamic Programming and Reinforcement Learning, IEEE, pp. 177–184.
  * Vinyals et al. (2019) Vinyals, O., Babuschkin, I., Czarnecki, W. M., Mathieu, M., Dudzik, A., Chung, J., Choi, D. H., Powell, R., Ewalds, T., Georgiev, P., et al. Grandmaster level in StarCraft II using multi-agent reinforcement learning. Nature, 2019, 575, 350–354.
  * Wackerly et al. (2008) Wackerly, D. D., Mendenhall, W., and Scheaffer, R. L., Mathematical Statistics with Applications (7th ed.), 2008, Belmont, CA: Thomson Brooks/Cole.
  * Wang et al. (2021) Wang, H., Lin, S., and Zhang, J. Adaptive Ensemble Q-learning: Minimizing Estimation Bias via Error Feedback. Advances in Neural Information Processing Systems, 2021, 34.
  * Watkins and Dayan (1992) Watkins, C. J. and Dayan, P. Q-learning. Machine learning, 1992, 8, 279–292.
  * Young and Tian (2019) Young, K. and Tian, T. Minatar: An atari-inspired testbed for thorough and reproducible reinforcement learning experiments. arXiv preprint arXiv:1903.03176, 2019.
  * Zhang et al. (2017) Zhang, Z., Pan, Z., and Kochenderfer, M. J. 2017, Weighted Double Q-learning. In International Joint Conference on Artificial Intelligence, pp. 3455–3461.
  * Zhu and Rigotti (2021) Zhu, R. and Rigotti, M. 2021, Self-correcting Q-learning. In AAAI Conference on Artificial Intelligence, vol. 35, pp. 11185–11192.

##  Appendix Appendix A Analytic forms for Section 3.3

###  A.1 Maximum Estimator

Consider the following setup: $M=2$,
$X_{1}\sim\mathcal{N}(\mu_{1},\sigma^{2})$,
$X_{2}\sim\mathcal{N}(\mu_{2},\sigma^{2})$, and given sample sizes $|S_{1}|$,
$|S_{2}|$, from which follows:
$\hat{\mu}_{1}\sim\mathcal{N}(\mu_{1},\frac{\sigma^{2}}{|S_{1}|})$ and
$\hat{\mu}_{2}\sim\mathcal{N}(\mu_{2},\frac{\sigma^{2}}{|S_{2}|})$. Regarding
the ME, we use the expectation in (1) to compute the bias. However, we can
alternatively use the following closed-form solutions (Nadarajah and Kotz,
2008):

$\displaystyle\operatorname{E}\left[\max(\hat{\mu}_{1},\hat{\mu}_{2})\right]$
$\displaystyle=\mu_{1}\Phi\left(\frac{\mu_{1}-\mu_{2}}{\theta}\right)+\mu_{2}\Phi\left(\frac{\mu_{2}-\mu_{1}}{\theta}\right)+\theta\phi\left(\frac{\mu_{1}-\mu_{2}}{\theta}\right),$
$\displaystyle\operatorname{E}\left\\{\left[\max(\hat{\mu}_{1},\hat{\mu}_{2})\right]^{2}\right\\}$
$\displaystyle=\left(\frac{\sigma^{2}}{|S_{1}|}+\mu_{1}^{2}\right)\Phi\left(\frac{\mu_{1}-\mu_{2}}{\theta}\right)+\left(\frac{\sigma^{2}}{|S_{2}|}+\mu_{2}^{2}\right)\Phi\left(\frac{\mu_{2}-\mu_{1}}{\theta}\right)$
$\displaystyle+(\mu_{1}+\mu_{2})\theta\phi\left(\frac{\mu_{1}-\mu_{2}}{\theta}\right),$

where $\phi$ is the standard Gaussian pdf and
$\theta=\sqrt{\frac{\sigma^{2}}{|S_{1}|}+\frac{\sigma^{2}}{|S_{2}|}}$. The
expectation of the squared ME can be used to compute the variance:

$\operatorname{Var}\left[\max(\hat{\mu}_{1},\hat{\mu}_{2})\right]=\operatorname{E}\left\\{\left[\max(\hat{\mu}_{1},\hat{\mu}_{2})\right]^{2}\right\\}-\operatorname{E}\left[\max(\hat{\mu}_{1},\hat{\mu}_{2})\right]^{2}.$

###  A.2 Double Estimator

The expectation of the DE is given in (3), which directly yields the bias. As
already mentioned, what we refer to as the DE throughout the paper is actually
the CVE whenever possible, thus we compute the variance of the latter for this
example. For notation, we use $\hat{\mu}_{i}^{A}=\hat{\mu}_{i}(S^{A})$ and
$\hat{f}_{i}^{A}$, $\hat{F}_{i}^{A}$ for the pdf and cdf of
$\hat{\mu}_{i}^{A}$, respectively, and similar for $S^{B}$. We assume that the
sample $S$ is split evenly between $S_{A}$ and $S_{B}$, so that the
theoretical mean distribution $\hat{f}_{i}^{A}$ equals $\hat{f}_{i}^{B}$. The
DE estimate when index selection is performed on subsample $S^{A}$ is denoted
with $\hat{\mu}^{DE,A}_{*}$, and when selecting based on $S^{B}$ with
$\hat{\mu}^{DE,B}_{*}$. It follows:

$\displaystyle\operatorname{Var}\left(\hat{\mu}^{CVE}_{*}\right)$
$\displaystyle=\operatorname{Var}\left(\frac{\hat{\mu}^{DE,A}_{*}+\hat{\mu}^{DE,B}_{*}}{2}\right)$
$\displaystyle=\frac{1}{4}\operatorname{Var}\left(\hat{\mu}^{DE,A}_{*}\right)+\frac{1}{4}\operatorname{Var}\left(\hat{\mu}^{DE,B}_{*}\right)+\frac{1}{2}\operatorname{Cov}\left(\hat{\mu}^{DE,A}_{*},\hat{\mu}^{DE,B}_{*}\right)$
$\displaystyle=\frac{1}{2}\operatorname{Var}\left(\hat{\mu}^{DE,A}_{*}\right)+\frac{1}{2}\operatorname{Cov}\left(\hat{\mu}^{DE,A}_{*},\hat{\mu}^{DE,B}_{*}\right),$
(18)

because
$\operatorname{Var}\left(\hat{\mu}^{DE,A}_{*}\right)=\operatorname{Var}\left(\hat{\mu}^{DE,B}_{*}\right)$.
Using definition:

$\operatorname{Var}\left(\hat{\mu}^{DE,A}_{*}\right)=\operatorname{E}\left[\left(\hat{\mu}^{DE,A}_{*}\right)^{2}\right]-\operatorname{E}\left[\hat{\mu}^{DE,A}_{*}\right]^{2},$
(19)

in which:

$\operatorname{E}\left[\left(\hat{\mu}^{DE,A}_{*}\right)^{2}\right]=\operatorname{E}\left[\left(\hat{\mu}^{B}_{1}\right)^{2}\right]\int_{-\infty}^{\infty}\hat{f}_{1}^{A}(x)\hat{F}_{2}^{A}(x)dx+\operatorname{E}\left[\left(\hat{\mu}^{B}_{2}\right)^{2}\right]\int_{-\infty}^{\infty}\hat{f}_{2}^{A}(x)\hat{F}_{1}^{A}(x)dx,$

where we compute:
$\operatorname{E}\left[\left(\hat{\mu}^{B}_{1}\right)^{2}\right]=\operatorname{Var}\left(\hat{\mu}^{B}_{1}\right)+\operatorname{E}\left[\hat{\mu}^{B}_{1}\right]^{2}$;
and $\operatorname{E}\left[\left(\hat{\mu}^{B}_{2}\right)^{2}\right]$
analogously, so that (19) is complete. To compute the covariance in (A.2), we
have:

$\operatorname{Cov}\left(\hat{\mu}^{DE,A}_{*},\hat{\mu}^{DE,B}_{*}\right)=\operatorname{E}\left[\hat{\mu}^{DE,A}_{*}\hat{\mu}^{DE,B}_{*}\right]-\operatorname{E}\left[\hat{\mu}^{DE,A}_{*}\right]\operatorname{E}\left[\hat{\mu}^{DE,B}_{*}\right],$

with the expectation of the product being

$\begin{aligned}
\operatorname{E}\left[\hat{\mu}^{DE,A}_{*}\hat{\mu}^{DE,B}_{*}\right]&=\operatorname{E}\left\\{\left[\mathcal{I}\left(\hat{\mu}_{1}^{A}>\hat{\mu}_{2}^{A}\right)\hat{\mu}_{1}^{B}+\mathcal{I}\left(\hat{\mu}_{1}^{A}\leq\hat{\mu}_{2}^{A}\right)\hat{\mu}_{2}^{B}\right]\left[\mathcal{I}\left(\hat{\mu}_{1}^{B}>\hat{\mu}_{2}^{B}\right)\hat{\mu}_{1}^{A}+\mathcal{I}\left(\hat{\mu}_{1}^{B}\leq\hat{\mu}_{2}^{B}\right)\hat{\mu}_{2}^{A}\right]\right\\}\\\
&=\operatorname{E}\left[\mathcal{I}\left(\hat{\mu}_{1}^{A}>\hat{\mu}_{2}^{A}\right)\hat{\mu}_{1}^{A}\right]\operatorname{E}\left[\mathcal{I}\left(\hat{\mu}_{1}^{B}>\hat{\mu}_{2}^{B}\right)\hat{\mu}_{1}^{B}\right]+\operatorname{E}\left[\mathcal{I}\left(\hat{\mu}_{1}^{A}>\hat{\mu}_{2}^{A}\right)\hat{\mu}_{2}^{A}\right]\operatorname{E}\left[\mathcal{I}\left(\hat{\mu}_{1}^{B}\leq\hat{\mu}_{2}^{B}\right)\hat{\mu}_{1}^{B}\right]\\\
&+\operatorname{E}\left[\mathcal{I}\left(\hat{\mu}_{1}^{A}\leq\hat{\mu}_{2}^{A}\right)\hat{\mu}_{1}^{A}\right]\operatorname{E}\left[\mathcal{I}\left(\hat{\mu}_{1}^{B}>\hat{\mu}_{2}^{B}\right)\hat{\mu}_{2}^{B}\right]+\operatorname{E}\left[\mathcal{I}\left(\hat{\mu}_{1}^{A}\leq\hat{\mu}_{2}^{A}\right)\hat{\mu}_{2}^{A}\right]\operatorname{E}\left[\mathcal{I}\left(\hat{\mu}_{1}^{B}\leq\hat{\mu}_{2}^{B}\right)\hat{\mu}_{2}^{B}\right].\end{aligned}$

This expression is simplified using
$\mathcal{I}\left(\hat{\mu}_{1}^{A}>\hat{\mu}_{2}^{A}\right)=1-\mathcal{I}\left(\hat{\mu}_{1}^{A}\leq\hat{\mu}_{2}^{A}\right)$
to get

$\operatorname{E}\left[\hat{\mu}^{DE,A}_{*}\hat{\mu}^{DE,B}_{*}\right]=\mu_{1}^{2}+2I_{1}(\mu_{2}-\mu_{1})+(I_{1}-I_{2})^{2},$

where

$\displaystyle I_{1}$
$\displaystyle=\operatorname{E}\left[\mathcal{I}(\hat{\mu}_{1}^{A}\leq\hat{\mu}_{2}^{A})\hat{\mu}_{1}^{A}\right]=\mu_{1}-\int_{-\infty}^{\infty}x\hat{f}_{1}^{A}(x)\hat{F}_{2}^{A}(x)dx,$
$\displaystyle I_{2}$
$\displaystyle=\operatorname{E}\left[\mathcal{I}(\hat{\mu}_{1}^{A}\leq\hat{\mu}_{2}^{A})\hat{\mu}_{2}^{A}\right]=\int_{-\infty}^{\infty}x\hat{f}_{2}^{A}(x)\hat{F}_{1}^{A}(x)dx.$

###  A.3 T-Estimator and K-Estimator

Regarding the expectation of the KE, we have:

$\displaystyle\operatorname{E}\left[\hat{\mu}^{KE}_{*}\right]$
$\displaystyle=\operatorname{E}\left[\hat{\mu}^{KE}_{*}\mathcal{I}(\hat{\mu}_{1}>\hat{\mu}_{2})\right]+\operatorname{E}\left[\hat{\mu}^{KE}_{*}\mathcal{I}(\hat{\mu}_{1}\leq\hat{\mu}_{2})\right]$
$\displaystyle=\operatorname{E}\left[\left\\{\left[\sum_{j=1}^{2}\kappa\left(\frac{\hat{\mu}_{j}-\hat{\mu}_{1}}{\theta_{1j}}\right)\right]^{-1}\sum_{j=1}^{2}\kappa\left(\frac{\hat{\mu}_{j}-\hat{\mu}_{1}}{\theta_{1j}}\right)\hat{\mu}_{j}\right\\}\mathcal{I}\left(\hat{\mu}_{1}>\hat{\mu}_{2}\right)\right]$
$\displaystyle+\operatorname{E}\left[\left\\{\left[\sum_{j=1}^{2}\kappa\left(\frac{\hat{\mu}_{j}-\hat{\mu}_{2}}{\theta_{2j}}\right)\right]^{-1}\sum_{j=1}^{2}\kappa\left(\frac{\hat{\mu}_{j}-\hat{\mu}_{2}}{\theta_{2j}}\right)\hat{\mu}_{j}\right\\}\mathcal{I}\left(\hat{\mu}_{1}\leq\hat{\mu}_{2}\right)\right]$
$\displaystyle=\operatorname{E}\left\\{\frac{1}{\kappa(0)+\kappa\left(\frac{\hat{\mu}_{2}-\hat{\mu}_{1}}{\theta_{12}}\right)}\left[\kappa(0)\hat{\mu}_{1}+\kappa\left(\frac{\hat{\mu}_{2}-\hat{\mu}_{1}}{\theta_{12}}\right)\hat{\mu}_{2}\right]\mathcal{I}(\hat{\mu}_{1}>\hat{\mu}_{2})\right\\}$
$\displaystyle+\operatorname{E}\left\\{\frac{1}{\kappa\left(\frac{\hat{\mu}_{1}-\hat{\mu}_{2}}{\theta_{21}}\right)+\kappa(0)}\left[\kappa\left(\frac{\hat{\mu}_{1}-\hat{\mu}_{2}}{\theta_{21}}\right)\hat{\mu}_{1}+\kappa(0)\hat{\mu}_{2}\right]\mathcal{I}(\hat{\mu}_{1}\leq\hat{\mu}_{2})\right\\}$
$\displaystyle=\int_{-\infty}^{\infty}\int_{-\infty}^{x_{1}}\frac{1}{\kappa(0)+\kappa\left(\frac{x_{2}-x_{1}}{\theta_{12}}\right)}\left[\kappa(0)x_{1}+\kappa\left(\frac{x_{2}-x_{1}}{\theta_{12}}\right)x_{2}\right]\hat{f}_{1}(x_{1})\hat{f}_{2}(x_{2})dx_{2}dx_{1}$
$\displaystyle+\int_{-\infty}^{\infty}\int_{-\infty}^{x_{2}}\frac{1}{\kappa\left(\frac{x_{1}-x_{2}}{\theta_{21}}\right)+\kappa(0)}\left[\kappa\left(\frac{x_{1}-x_{2}}{\theta_{21}}\right)x_{1}+\kappa(0)x_{2}\right]\hat{f}_{1}(x_{1})\hat{f}_{2}(x_{2})dx_{1}dx_{2},$

where
$\theta_{ij}=\sqrt{\frac{\sigma^{2}}{|S_{i}|}+\frac{\sigma^{2}}{|S_{j}|}}$ and
$\hat{f}_{i}$ is the pdf of $\hat{\mu}_{i}$. For the variance, we can compute
$\operatorname{E}\left\\{[\hat{\mu}^{KE}_{*}]^{2}\right\\}$ analogously. Since
the TE is a special case of the KE with $\kappa(T)=\mathcal{I}(T\geq
z_{\alpha})$, the above formula is also applicable for TE.

##  Appendix Appendix B Adaptive TE-BDQN Algorithm

initialize Action-value estimate networks with $K$ outputs
$\left\\{\hat{Q}^{*}_{k}\right\\}^{K}_{k=1}$, masking distribution $M$, empty
replay buffer $D$

repeat

Initialize $s$

Pick a value function to act: $k\sim\text{Uniform}\\{1,\ldots,K\\}$

repeat

Choose action $a$ from state $s$ with greedy policy derived from
$\hat{Q}^{*}_{k}$

Take action $a$, observe reward $r$ and next state $s^{\prime}$

Sample bootstrap mask $m\sim M$

Add $(s,a,r,s^{\prime},m)$ to replay buffer $D$

Sample random minibatch of transitions
$\left\\{(s_{i},a_{i},s^{\prime}_{i},r_{i},m_{i})\right\\}_{i=1}^{B}$ from $D$

Perform gradient descent step based on (15)

Every $C$ steps:

Reset $\theta_{k}=\theta_{k}^{-}$ for $k=1,\ldots,K$

Run partial episodes to update $\alpha$ via:
$\alpha\leftarrow\alpha+\frac{\tau_{\rm
Ada}}{K}\sum_{k=1}^{K}\sum_{t=1}^{T_{\rm
Ada}}\left[R_{k}(s_{t,k},a_{t,k})-\hat{Q}^{*}_{k}(s_{t,k},a_{t,k};\theta_{k})\right]$

$s\leftarrow s^{\prime}$

until  _$s$ is terminal_

until  __

Algorithm 4 Ada-TE-BDQN

##  Appendix Appendix C Hyperparameters in MinAtar

Hyperparameter | Value
---|---
Batch size ($B$) | 32
Discount factor ($\gamma$) | 0.99
Loss function | MSE
Min. replay buffer size | $5\,000$
Max. replay buffer size | $100\,000$
Optimizer | Adam
Target network update frequency ($C$) | $1\,000$
Initial exploration rate* ($\epsilon_{\rm initial}$) | 1.0
Final exploration rate* ($\epsilon_{\rm final}$) | 0.1
Test exploration rate* ($\epsilon_{\rm test}$) | 0.0
Exploration steps* | $100\,000$
Bernoulli mask probability† ($p$) | 1.0
Number of bootstrap heads† ($K$) | 10
Initial bias parameter‡ ($\alpha$) | 0.25
Time horizon‡ ($T_{\rm Ada})$ | 32
Table 2: List of hyperparameters used in the MinAtar experiments. Parameters
with a * are used by DQN, DDQN, and SCDQN, while the ones with a † are
relevant for BDQN, TE-BDQN, KE-BDQN, and Ada-TE-BDQN. An ‡ exclusively refers
to Ada-TE-BDQN.

Table 2 details the settings for the experiments in MinAtar (Young and Tian,
2019). All algorithms were implemented using PyTorch (Paszke et al., 2019) and
the computation was performed on Intel(R) Xeon(R) CPUs E5-2680 v3 (12 cores) @
2.50GHz. The source code is available at: https://github.com/MarWaltz/TUD_RL.
Note that we replaced some extreme outlier seeds for the bootstrap-based
algorithms in Breakout and Seaquest environments. For example, the algorithm
TE-BDQN led to an astonishing peak performance with a test return of over 200
in a Breakout run, while it got stuck in a rare occasion on Seaquest.
Including those exceptions would paint an unrealistic picture of the actual
capabilities of the algorithm. Similar was observed for all bootstrap-based
algorithms, and we argue that those rare instabilities are due to the
algorithm’s dependence on the initialization of the bootstrap heads.