Off-Policy Evaluation of Bandit Algorithm from Dependent Samples under Batch Update Policy

The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.

50,117 characters

Off-Policy Evaluation of Bandit Algorithm from Dependent Samples under Batch Update Policy



\title{Off-Policy Evaluation of Bandit Algorithm from Dependent Samples under Batch Update Policy}


\author{
Masahiro Kato,\ \ \ \ \ Yusuke Kaneko\\
CyberAgent Inc.\\
\texttt{masahiro\[email removed]}\\
\texttt{kaneko\[email removed]}
}

\maketitle

\begin{abstract}
The goal of \emph{off-policy evaluation} (OPE) is to evaluate a new policy using historical data obtained via a \emph{behavior policy}. However, because the contextual bandit algorithm updates the policy based on past observations, the samples are not \emph{independent and identically distributed} (i.i.d.). This paper tackles this problem by constructing an estimator from a \emph{martingale difference sequence} (MDS) for the dependent samples. In the data-generating process, we do not assume the convergence of the policy, but the policy uses the same conditional probability of choosing an action during a certain period. Then, we derive an asymptotically normal estimator of the value of an \emph{evaluation policy}. As another advantage of our method, the batch-based approach simultaneously solves deficient support problem. Using benchmark and real-world datasets, we experimentally confirm the effectiveness of the proposed method.
\end{abstract}

\section{Introduction}
As an instance of sequential decision-making problems, the \emph{multi-armed bandit} (MAB) algorithms have attracted significant attention in various applications, such as ad optimization, personalized medicine, search engines, and recommendation systems.  Recently, various methods for evaluating a new policy using historical data obtained via the MAB algorithms \citep{kdd2009_ads,li2010contextual} have emerged. The goal of \emph{off-policy evaluation} (OPE) is to evaluate a new policy by estimating the expected reward obtained from the new policy \citep{dudik2011doubly,wang2017optimal,narita2019counterfactual,pmlr-v97-bibaut19a,Kallus2019IntrinsicallyES,Oberst2019}. Although an OPE algorithm estimates the expected reward from a new policy, most existing studies presume that the samples are \emph{independent and identically distributed} (i.i.d.). However, the MAB algorithm policy updates the probability of choosing an action based on past observations, and samples are not i.i.d. owing to this update. In this case, such existing studies do not guarantee that their estimators have asymptotic normality and $\sqrt{T}$-consistency for a sample size $T$. Therefore, there is a strong motivation to establish a novel method for OPE from dependent samples.

Several pioneering studies address OPE from dependent samples \citep{Laan2008TheCA,Hahn2011,Laan2016onlinetml,Luedtke2016,hadad2019,Kato2020}. We can group the methods for deriving asymptotic normality into the following three groups: (a) \citet{Laan2008TheCA}, \citet{Laan2016onlinetml}, \citet{hadad2019}, and \citet{Kato2020} derive the asymptotic normality with the central limit theorem (CLT) of a martingale difference sequence (MDS) by assuming that the probability of choosing an action converges to a time-invariant probability; (b) \citet{Luedtke2016} derives the asymptotic normality by standardizing a MDS, which is also used in the first group; (c) \citet{Hahn2011} apply asymptotic theory for the batched probability update process.

This paper focuses on an approach of the third group; that is, there are sufficiently large sample sizes in each batch. Compared with \citet{Hahn2011}, our proposed method is more general and applicable in practical applications. Our method has the following three advantages compared with existing studies: (i) it does not assume convergence of the probability of choosing an action; (ii) it allows the probability of choosing an action to be $0$ for some actions in some batches; (iii) we can also use non-Donsker nuisance estimators as well as \cite{Laan2016onlinetml}.

This paper has three main contributions. First, we provide a solution for OPE from dependent samples obtained via the MAB algorithms. Second, under the batch update policy, the proposed estimator achieves the asymptotic normality with fewer assumptions. Third, the estimator also experimentally shows a lower mean squared error (MSE) in some cases.

\section{Problem Setting}
Here, we formulate OPE under a batch update.

\subsection{Date-Generating Process}
Let $A_t$ be an action taking variable in $\mathcal{A}=\{1,2,\dots,K\}$, $X_t$ the \emph{covariate} observed by the decision maker when choosing an action , and $\mathcal{X}$ the space of covariate. Let us denote a random variable of a reward at period $t$ as $Y_t=\sum^K_{a=1}\mathbbm{1}[A_t = a]Y_t(a)$, where $Y_t:\mathcal{A}\to\mathbb{R}$ is a potential outcome\footnote{We can express the DGP without using the potential outcome variable \citep{kato_uehara_2020}.}. In this paper, we have access to a dataset $\{(X_t, A_t, Y_t)\}^{T}_{t=1}$ with the following data-generating process (DGP):
\begin{align}
\label{eq:DGP}
\big\{(X_t, A_t, Y_t)\big\}^{T}_{t=1}\sim p(x)\pi_t(a\mid x, \Omega_{t-1})p(y\mid a, x),
\end{align}
where $\Omega_{t-1}\in\mathcal{M}_{t-1}$ denotes the history until $t-1$ period defined as $\Omega_{t-1}=\{X_{t-1}, A_{t-1}, Y_{t-1}, \dots, X_{1}, A_1, Y_{1}\}$ with the space $\mathcal{M}_{t-1}$, $p(x)$ denotes the density of the covariate $X_t$, $\pi_t(a\mid x, \Omega_{t-1})$ denotes the probability of choosing an action $A_t$ conditioned on $X_t$ and $\Omega_{t-1}$, and $p(y\mid a, x)$ denotes the density of an outcome $Y_t$ conditioned on $A_t$ and $X_t$. We assume that $p(x)$ and $p(y\mid a, x)$ are invariant across periods, but $\pi_t(a\mid x, \Omega_{t-1})$ can take different values across periods. Let us call a policy inducing $\pi_t(a\mid x, \Omega_{t-1})$ a \emph{behavior policy}.

\subsection{Off-Policy Evaluation}
\label{sec:opeopl}
This paper considers estimating the value of an \emph{evaluation policy} using samples obtained under the behavior policy. Let an evaluation policy $\pi^\mathrm{e}:\mathcal{A}\times\mathcal{X}\to[0,1]$ be a probability of choosing an action $A_t$ conditioned on a covariate $X_t$. We are interested in estimating the expected reward from any pre-specified evaluation policy $\pi^{\mathrm{e}}(a \mid x)$. Then, we define the expected reward under an evaluation policy as $R(\pi^\mathrm{e}) := \mathbb{E}\left[\sum^K_{a=1}\pi^\mathrm{e}(a \mid x)Y_t(a)\right]$. For brevity, we also denote $R(\pi^\mathrm{e})$ as $\theta_0$. The goal of this paper is to estimate $\theta_0=R(\pi^{\mathrm{e}})$ using dependent samples under a batch update policy. To identify $\theta_0$, we assume overlaps in policy and the boundedness of the outcome.
\begin{assumption}\label{asm:overlap_pol}
There exists a constant $C_1$ such that $0\leq \frac{\pi^\mathrm{e}(a\mid x)}{\pi_t(a\mid x, \Omega_{t-1})}\leq C_1$.
\end{assumption}
\begin{assumption}\label{asm:overlap_outcome}
There exists a constant $C_2$ such that $|Y| \leq C_2$.
\end{assumption}

\begin{remark}[Existing Methods for OPE]
We review three types of standard estimators of $R(\pi^{\mathrm{e}})$ under the case where $\pi_1(a\mid x, \Omega_{0})=\pi_2(a\mid x, \Omega_{1})=\cdots=\pi_T(a\mid x, \Omega_{T-1})=p(a\mid x)$ in the DGP defined in (\ref{eq:DGP}). The first estimator is an inverse probability weighting (IPW) estimator given by $\frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\frac{\pi^\mathrm{e}(a\mid X_t)\mathbbm{1}[A_t=a]Y_t}{p(a\mid X_t)}$ \citep{rubin87,hirano2003efficient,swaminathan15a}. Although this estimator is unbiased when the behavior policy is known, it suffers from high variance. The second estimator is a direct method (DM) estimator $\frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\hat{f}_{T}(a, X_t)$, where $\hat{f}_{T}(a, X_t)$ is an estimator of $f^*(a, X_t)$ \citep{HahnJinyong1998OtRo}. This estimator is known to be weak against model misspecification for $f^*(a, X_t)$. The third estimator is an augmented IPW (AIPW) defined as $\frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\Bigg(\frac{\pi^\mathrm{e}(a\mid X_t)\mathbbm{1}[A_t=a]\big(Y_t - \hat{f}_{T}(a, X_t)\big)}{p(a\mid X_t)} + \pi^\mathrm{e}(a\mid X_t)\hat{f}_{T}(a, X_t)\Bigg)$ \citep{robins94,ChernozhukovVictor2018Dmlf}.
Under certain conditions, it is known that this estimator achieves the efficiency bound (a.k.a semiparametric lower bound), which is the lower bound of the asymptotic MSE of OPE among regular $\sqrt{T}$-consistent estimators \citep{VaartA.W.vander1998As}.
\end{remark}

\begin{remark}[Semiparametric Lower Bound]
\label{rem:semi_low}
The lower bound of the variance is defined for an estimator of $\theta_0$ under some posited models of the DGP. If this posited model is a parametric model, it is equal to the Cram\'{e}r-Rao lower bound. When this posited model is a semiparametric model, we can define a corresponding Cram\'{e}r-Rao lower bound \citet{bickel98}. \citet{narita2019counterfactual} gives the semiparametric lower bound of the DGP (\ref{eq:DGP}) under $\pi_1(a\mid x, \Omega_{0})=\cdots=\pi_T(a\mid x, \Omega_{T-1})=p(a\mid x)$ as $\mathbb{E}\Big[\sum^{K}_{a=1}\frac{\big(\pi^\mathrm{e}(a\mid X)\big)^2v^*(a, X_t)}{p(a\mid X)}+\left(\sum^{K}_{a=1}\pi^\mathrm{e}(a\mid X)f^*(a, X_t) - \theta_0\right)^2\Big]$.
\end{remark}

\paragraph{Notations:} Let us denote $\mathbb{E}[Y_t(a)\mid x]$ and $\mathrm{Var}(Y_t(a)\mid x)$ as $f^*(a, x)$ and $v^*(a, x)$, respectively. Let $\mathcal{F}$ be the class of $f^*(a, x)$. Let $\hat{f}_{t}(a, x\mid \Omega_{t-1})$ be an estimator of $f^*(a, x)$ constructed from $\Omega_{t-1}$, respectively. Let $\mathcal{N}(\mu, \mathrm{var})$ be the normal distribution with the mean $\mu$ and the variance $\mathrm{var}$. For a random variable $Z$ and function $\mu$, let $\|\mu(Z)\|_2=\int |\mu(z)|^2 p(z) dz $ be the $L^{2}$-norm.

\subsection{Patterns of Probability Update}
In this paper, based on the update of $\pi_t(a\mid x, \Omega_{t-1})$, we classify the policies into two patterns, \emph{sequential update policy} and \emph{batch update policy}. For the sequential update policy, the policy updates $\pi_t(a\mid x, \Omega_{t-1})$ at each period \citep{Laan2008TheCA,Laan2016onlinetml,Kato2020}. Under the batch update policy, after the policy continues using a fixed probability $\pi_t(a\mid x, \Omega_{t-1})$ for some periods without updates, the policy updates $\pi_t(a\mid x, \Omega_{t-1})$ \citep{Hahn2011,narita2019counterfactual}. Although the sequential update is standard in the MAB problem, we often apply batch updates in industrial applications such as ad-optimization \citep{narita2019counterfactual}. For OPE under the sequential update, \citet{Laan2008TheCA}, \citet{Laan2016onlinetml}, \citet{Luedtke2016}, \citet{hadad2019}, and \citet{Kato2020} proposed estimators with the asymptotic normality. For instance, an adaptive AIPW (A2IPW) estimator \citep{Laan2016onlinetml,hadad2019,Kato2020} has asymptotic normality if the behavior policy converges. On the other hand, we consider OPE under a batch update. Let $M$ denote the number of updates and $\tau\in I = \{1,2,\dots,M\}$ denotes the batch index. For $\tau\in I$, the probability is updated at a period $t_\tau$, where $t_\tau - t_{\tau-1} = Tr_\tau$, using samples $\{(X_t, Y_t, A_t)\}^{t_{\tau}}_{t=t_{\tau-1}}$, where $r_1 + r_2 + \cdots + r_M = 1$ and $t_0 = 0$. Thus, in addition to the DGP (\ref{eq:DGP}), we assume
\begin{align*}
\big\{(X_t, A_t, Y_t)\big\}^{t_\tau}_{t=t_{\tau-1}} \stackrel{i.i.d.}{\sim}  p(x)\pi_\tau(a\mid x, \Omega_{t_{\tau-1}})p(y\mid a, x),
\end{align*}
where $\pi_\tau(a\mid x, \Omega_{t_{\tau-1}})$ denotes the assignment probability updated based on samples until the period $t_{\tau-1}$.

\begin{comment}
\begin{remark}[Policy and Probability]
To transform a function $\pi_\tau(a\mid x, \Omega_{t_{\tau-1}})$ to a probability, we consider the following process. First, we draw a random variable $\xi_t$ following the uniform distribution on $[0,1]$ at the $t$-th round. Then, at each round $t$, we select an action $A_t$ such that $\sum^{A_t-1}_{a=0}\pi_\tau(a\mid X_t, \Omega_{t_{\tau-1}}) \leq \xi_t \leq \sum^{A_t}_{a=0}\pi_\tau(a\mid X_t, \Omega_{t_{\tau-1}})$, where $\pi_t(0 \mid X_t, \Omega_{t-1})=0$. Under this process, we regard the probability $p_t(a\mid x)$ as $\pi_\tau(a\mid x, \Omega_{t_{\tau-1}})$.
\end{remark}
\end{comment}

\subsection{Related Work}
For the sequential update, \citet{Laan2008TheCA}, \citet{Laan2016onlinetml}, \citet{hadad2019}, and \citet{Kato2020} assume that the probability of choosing an action converges to a time-invariant function almost certainly; that is, $\pi_t(a\mid x, \Omega_{t-1})\xrightarrow{\mathrm{p}} \alpha(a\mid x)$, where $\alpha:\mathcal{X}\to (0,1)$. This assumption enables us to apply the CLT for MDS. \citet{Laan2016onlinetml} proposed constructing step-wise nuisance estimators, which enables us to derive asymptotic normality without Donsker's conditions of nuisance estimators. This technique is a generalization of sample-splitting, which is also called cross-fitting in a context of double/debiased machine learning \citep{klaassen1987,ZhengWenjing2011CTME,ChernozhukovVictor2018Dmlf}. For the A2IPW estimator, \citet{hadad2019} proposed using an adaptive weight for stabilizing the behavior, and \citet{Kato2020} derived concentration inequality based on the law of iterated logarithms. On the other hand, we also construct an MDS and apply the CLT, but do not assume $\pi_t(a\mid x, \Omega_{t-1})\xrightarrow{\mathrm{p}} \alpha(a\mid x)$ by using batch update policy. Instead, we assume a sufficient sample size for each batch.

For such a non-stationary setting, \citet{Luedtke2016} also proposed an estimator with asymptotically normality for sequential update policy without using batch update policy. For deriving the asymptotic normality, \citet{Luedtke2016} used standardization for an MDS. Although the method enables us to construct an asymptotically normal estimator for various estimators, the proposed estimator only has $\sqrt{T - \ell}$-consistency for another sample size $\ell > 0$, not $\sqrt{T}$, to estimate the variances of the MDS.

As other related work, in the MAB problem, \citet{perchet2016} considered the setting of batch policy updates. In OPE, \citet{narita2019counterfactual} also discuss a similar problem setting, but they assume that samples are i.i.d. Independently, \citet{kelly2020} provided a method for deriving a confidence interval of an ordinary least squares estimator, which is a different parameter of what we want to estimate.

\section{OPE under Batch Update Policy}
This section introduces a concept for conducting OPE under the batch update policy and a method based on the concept with its theoretical properties.

\subsection{Strategy for OPE}
For OPE under batch update policy, we consider asymptotic properties based on the assumption of $t_\tau - t_{\tau-1} \to \infty$ as $T\to \infty$ for fixed $\tau$. Because $\big\{(X_t, A_t, Y_t)\big\}^{t_\tau}_{t=t_{\tau-1}}$ is i.i.d., we can use the standard limit theorems for the partial sum of the samples to obtain an asymptotically normal estimator of $\theta_0 = R(\pi^\mathrm{e})$. However, we also have the motivation to use all samples together to increase the efficiency of the estimator. Therefore, based on the idea of \emph{generalized method of moments} (GMM), we propose an estimator of $\theta_0$ considering the sample averages of each block as an empirical moment conditions. The main difference from the standard GMM is the assumption that the samples are not i.i.d. However, for the case under the batch update, we can apply the central limit theorems (CLT) for the martingale difference sequences (MDS) by appropriately constructing an estimator. We describe the proposed method as follows.

\subsection{Estimator for OPE}
We propose an estimator of $\theta_0$ based on a idea of GMM. For an index of batch $\tau\in I$, a function $f\in\mathcal{F}$ such that $f:\mathcal{A}\times \mathcal{X} \to \mathbb{R}$ and an evaluation policy $\pi^\mathrm{e} \in \Pi$, we define $h^{\mathrm{OPE}}_t:\mathcal{X}\times\mathcal{A}\times\mathbb{R}\times I \times \Theta\times\mathcal{F}\times\Pi \to \mathbb{R}$ as $h^{\mathrm{OPE}}_t(x, k, y; \tau, \theta, f, \pi^\mathrm{e}) = \frac{1}{r_\tau}\eta_t(x, k, y; \tau, \theta, f, \pi^\mathrm{e})\mathbbm{1}\big[t_{\tau - 1} < t \leq t_{\tau} \big]$, where $\eta_t(x, k, y; \tau, \theta, f, \pi^\mathrm{e}) := \phi_t(x, k, y; \tau, f, \pi^\mathrm{e}) - \theta$ and $\phi_t(x, k, y; \tau, f, \pi^\mathrm{e}):=$
\begin{align*}
&\sum^{K}_{a=1}\pi^\mathrm{e}(a\mid x)\Bigg\{\frac{\mathbbm{1}[k=a]\big\{y - f(a, x)\big\}}{\pi_{\tau}(a\mid x, \Omega_{t_{\tau - 1}})}+f(a, x)\Bigg\}.
\end{align*}
Let us note that, for $\tau\in I$, $\theta_0\in \Theta$, $f_{t-1}\in \mathcal{F}$, and $\pi^\mathrm{e}\in\Pi$, the sequence $\big\{h^{\mathrm{OPE}}_t(X_t, A_t, Y_t; \tau, \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e})\big\}^{T}_{t=1}$ is an MDS: for $h^{\mathrm{OPE}}_t(X_t, A_t, Y_t; \tau, \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e})$, by  $\mathbb{E}[\mathbbm{1}[A_t=a]\mid X_t, \Omega_{t-1}]=\pi_{\tau}(a\mid X_t, \Omega_{t_{\tau - 1}})$, we have
\begin{align*}
&\mathbb{E}\left[h^{\mathrm{OPE}}_t(X_t, A_t, Y_t; \tau, \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e})\mid \Omega_{t-1}\right]\\
& = \mathbb{E}\left[\frac{\mathbbm{1}\big[t_{\tau - 1} < t \leq t_{\tau} \big]}{r_\tau}\eta_t(x, k, y; \tau, \theta, \hat{f}_{t-1}, \pi^\mathrm{e})\mid \Omega_{t-1}\right]\\
& = \frac{\mathbbm{1}\big[t_{\tau - 1} < t \leq t_{\tau} \big]}{r_\tau}\mathbb{E}\left[\eta_t(x, k, y; \tau, \theta, \hat{f}_{t-1}, \pi^\mathrm{e})\mid \Omega_{t-1}\right]\\
& = \frac{\mathbbm{1}\big[t_{\tau - 1} < t \leq t_{\tau} \big]}{r_\tau}\times 0 = 0.
\end{align*}
Let us also define $\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta, \hat{f}_{t-1}, \pi^\mathrm{e}\right):=$
\begin{align*}
\begin{pmatrix}
h^{\mathrm{OPE}}_t(X_t, A_t, Y_t; 1, \theta, \hat{f}_{t-1}, \pi^\mathrm{e}) \\
h^{\mathrm{OPE}}_t(X_t, A_t, Y_t; 2, \theta, \hat{f}_{t-1}, \pi^\mathrm{e}) \\
 \vdots \\
h^{\mathrm{OPE}}_t(X_t, A_t, Y_t; M, \theta, \hat{f}_{t-1}, \pi^\mathrm{e})
\end{pmatrix}.
\end{align*}
Then, the sequence $\left\{\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e}\right)\right\}^{T}_{t=1}$ is an MDS with respect to $\big\{\Omega_t\big\}^{T-1}_{t=0}$; that is, $\mathbb{E}\left[\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e}\right)\mid \Omega_{t-1}\right] = \bm{0}$. Using the sequence $\left\{\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta, \hat{f}_{t-1}, \pi^\mathrm{e}\right)\right\}^{T}_{t=1}$, we define an estimator of OPE as $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e}) := $
\begin{align}
\label{gmmdm_ope}
\operatorname*{arg\,min}_{\theta\in\Theta} \left(\hat{\bm{q}}^{\mathrm{OPE}}_T(\theta)\right)^\top \hat{W}_T \left(\hat{\bm{q}}^{\mathrm{OPE}}_T(\theta)\right),
\end{align}
where $\hat{\bm{q}}^{\mathrm{OPE}}_T(\theta) = \frac{1}{T}\sum^{T}_{t=1}\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta, \hat{f}_{t-1}, \pi^\mathrm{e}\right)$ and $\hat{W}_T$ is a data-dependent $(M\times M)$-dimensional positive semi-definite matrix. Let us note that the estimator defined in Eq.~(\ref{gmmdm_ope}) is an application of GMM with the moment condition $\bm{q}^{\mathrm{OPE}}(\theta_0) = \mathbb{E}\left[\frac{1}{T}\sum^T_{t=1}\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e}\right)\right]=0$. For the minimization problem defined in Eq.~(\ref{gmmdm_ope}), we can analytically calculate the minimizer as $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e}) = w^\top_T D_T(\pi^\mathrm{e})$, where $w_T = (w_{T, 1}\ \cdots\ w_{T, M})^\top$ is an $M$-dimensional vector such that $\sum^M_{\tau=1}w_{T, \tau}=1$, and $D_T(\pi^\mathrm{e})$ is
\begin{align*}
\begin{pmatrix}
\frac{1}{t_{1}}\sum^{t_{1}}_{t=1} \phi_t(X_t, A_t, Y_t; 1, \hat{f}_{t-1}, \pi^\mathrm{e}) \\
\frac{1}{t_{2} - {t_{1}}}\sum^{t_2}_{t=t_1+1} \phi_t(X_t, A_t, Y_t; 2, \hat{f}_{t-1}, \pi^\mathrm{e})\\
\vdots\\
\frac{1}{T-t_{M-1}}\sum^{T}_{t=t_{M-1}+1} \phi_t(X_t, A_t, Y_t; M, \hat{f}_{t-1}, \pi^\mathrm{e})\\
\end{pmatrix}.
\end{align*}
We call the estimator a \emph{Batch-based Adaptive AIPW} (BA2IPW) estimator. In Appendix~\ref{appdx:sec:gmmdm}, we discuss the GMM perspective in more detail.

\subsection{Asymptotic Properties}
\label{sec:asymp_ba2ipw}
Here, we show the consistency and asymptotic normality of the proposed BA2IPW estimator $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e})$.

\begin{theorem}[Consistency of the BA2IPW Estimator]
\label{thm:consistency}
Suppose that there exists a constant $C_f > 0$ such that $\big|f_{t_{\tau-1}}(a, x)\big| < C_f$ for $\tau \in I$. Then, under Assumptions~\ref{asm:overlap_pol} and \ref{asm:overlap_outcome}, $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e})\xrightarrow{\mathrm{p}}\theta_0$.
\end{theorem}
\begin{proof}
We use the law of large numbers for an MDS from the boundedness of $\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e}\right)$, and we have $\frac{1}{T}\sum^T_{t=1}\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e}\right)\xrightarrow{\mathrm{p}}0$ (Proposition~\ref{prp:mrtgl_WLLN} in Appendix~\ref{sec:prelim}). This result means that $D_T(\pi^\mathrm{e})\xrightarrow{\mathrm{p}}I\theta_0$, where $I = (1\ 1\ \cdots\ 1)^\top$ is an $M$-dimensional vector. Therefore, $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e}) = w_T D_T(\pi^\mathrm{e}) \xrightarrow{\mathrm{p}} w_T I \theta_0 = \theta_0$.
\end{proof}
\begin{theorem}[Asymptotic Distribution of the BA2IPW Estimator]
Suppose that (i) $w_T = (w_{T,1}\ \cdots\ w_{T,M})^\top \xrightarrow{\mathrm{p}} w = (w_{1}\ \cdots\ w_{M})^\top$; (ii) $w_{T,\tau} > 0$ and $\sum^M_{\tau=1}w_{T,\tau} = 1$; (iii) $\hat{f}_{t-1}(a, x) \xrightarrow{\mathrm{p}} f^*(a, x)$ for all $a \in \mathcal{A}$ and $x \in \mathcal{X}$; (iv) There exists a constant $C_f > 0$ such that $\big|\hat{f}_{t-1}(a, x)\big| < C_f$. Then, under Assumptions~\ref{asm:overlap_pol} and \ref{asm:overlap_outcome}, $\sqrt{T}\big(\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e}) - \theta_0\big)  \xrightarrow{\mathrm{d}}\mathcal{N}\big(0, \sigma^2\big)$, where $\sigma^2 = \sum^M_{\tau=1}w_{\tau}\sigma^2_{\tau}$ and $\sigma^2_{\tau}=\frac{1}{r_\tau}\mathbb{E}\Big[\sum^{K}_{a=1}\frac{\big(\pi^\mathrm{e}(a\mid X)\big)^2\nu^*(a, X)}{\pi_{\tau}(a\mid X, \Omega_{t_{\tau - 1}})} + \left(\sum^{K}_{a=1}\pi^\mathrm{e}(a\mid X)f^*(a, X) - \theta_0\right)^2\Big]$.
\end{theorem}
The proof is shown in Appendix~\ref{appdx:main}. Readers might consider that the use of MDS for deriving the asymptotic normality is unnecessary. We discuss the necessity of MDS in Appendix~\ref{appdx:nec_mds}. We can also define a corresponding \emph{Batch-based Adaptive IPW} (BAdaIPW). For the BAdaIPW estimator, the variance of a batch $\tau$ is  $\sigma^2_{\mathrm{IPW}, \tau}=\frac{1}{r_\tau}\mathbb{E}\left[\sum^{K}_{a=1}\frac{\big(\pi^\mathrm{e}(a\mid X)\big)^2\mathbb{E}[Y^2_t\mid X_t]}{\pi_{\tau}(a\mid X, \Omega_{t_{\tau - 1}})} - \theta^2_0\right]$.

\begin{remark}[Construction of $f_{t_m}$ and Donsker Condition]
As well as the cross-fitting of double/debiased machine learning proposed by \citet{klaassen1987},\citet{ZhengWenjing2011CTME},\citet{Laan2016onlinetml}, and \citet{ChernozhukovVictor2018Dmlf}, the proposed estimator does not require Donsker's condition for asymptotic normality. This property comes from the MDS as pointed by \citet{Laan2016onlinetml}. On the other hand, because the samples are not independent, we cannot use the standard regression to obtain a consistent estimator of $f^*$. For example, \citet{yang2002} propose a nonparametric method for the bandit process under some mild conditions.
\end{remark}

\subsection{Weight of the Proposed Estimator}
\label{sec:efficiency}
Next, we discuss the choice of weight $w_T$.

\paragraph{Equal Weight:} A naive choice is weighting the moment conditions equality; that is, $w_{T, \tau} = \frac{1}{M}$. In this case, the proposed estimator boils down to $I^\top D_T(\pi^\mathrm{e})$, which is almost the same as the A2IPW estimator. Although the estimator itself is similar to the A2IPW estimator, the theoretical guarantee for the asymptotic normality is different. While the A2IPW estimator uses the assumption that the policy converges to a time-invariant policy, the proposed BA2IPW estimator uses the assumption of the batch update. We call the BA2IPW estimator with the equal weight a Plain BA2IPW (PBA2IPW) estimator.

\paragraph{Efficient Weight:} First, we consider an efficient weight $w_T$ that minimizes the asymptotic variance of $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e})$. As well as the standard GMM, the $\tau$-th element of the efficient weight is given as $w^*_{\tau}=\frac{1}{\sigma^2_\tau}/\sum^M_{\tau'=1}\frac{1}{\sigma^2_{\tau'}}$ \citep{GVK126800421}. Here, we use the orthogonality among moment conditions; that is, zero covariance. In this case, the asymptotic variance becomes $1/\sum^M_{\tau'=1}\frac{1}{\sigma^2_{\tau'}}$. Therefore, for gaining efficiency, we use a weight $\hat{w}_{T,\tau} = \frac{1}{ \hat{\sigma}^2_{T, \tau}}/\sum^M_{\tau'=1}\frac{1}{\hat{\sigma}^2_{T, \tau'}}$, where $\hat{\sigma}^2_{T, \tau}$ is an estimator of $\sigma^2_\tau$. We call the BA2IPW estimator with the efficient weight an Efficient BA2IPW (EBA2IPW) estimator.

\begin{remark}[Estimation of $w^*$]
If $\hat{\sigma}^2_{T, \tau}\xrightarrow{\mathrm{p}}\sigma^2_\tau$, we also have $\hat{w}_{T,\tau} \xrightarrow{\mathrm{p}}w^*_\tau$ from the continuous mapping theorem. In this paper, we propose two estimators defined as $\hat{\sigma}^2_{T, \tau}=\frac{1}{r_\tau T}\sum^{T}_{t=1}\left\{h^{\mathrm{OPE}}_t(x, a, y; \tau, \widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e}), \hat{f}_{t-1}, \pi^\mathrm{e})\right\}^2\times$ $\mathbbm{1}\big[t_{\tau - 1} < t \leq t_{\tau} \big]$ and $\tilde{\sigma}^2_{T, \tau}=\frac{1}{r_\tau T}\sum^{T}_{t=1}\left\{h^{\mathrm{OPE}}_t(x, a, y; \tau, \widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e}), \hat{f}_{T}, \pi^\mathrm{e})\right\}^2\times$ $\mathbbm{1}\big[t_{\tau - 1} < t \leq t_{\tau} \big]$.
\end{remark}

\paragraph{Weight for Numerical Stability:}
As explained above, we can obtain an efficient weight that minimizes the asymptotic variance. However, as \citet{hadad2019} pointed out, when using time-variance estimators $f_t$, early variable $\phi_t(X_t, A_t, Y_t; \tau, f_t, \pi^\mathrm{e})$ might be unstable due to the existence of an inaccurate estimators $f_t$ at early stages. \citet{hadad2019} proposed an adaptive weight that puts more weight on later variables $\phi_t(X_t, A_t, Y_t; \tau, f_t, \pi^\mathrm{e})$ in the estimator. In our problem setting, we can also introduce such weights. However, unlike the adaptive weight proposed by \citet{hadad2019}, which must be martingale, we do not require martingales on the weights. This property is a benefit of the batch update. For instance, for stabilization, we define a weight $\frac{1}{\ddot{\sigma}^2_{T, \tau}}/\sum^M_{\tau'=1}\frac{1}{\ddot{\sigma}^2_{T, \tau'}}$, where $\ddot{\sigma}^2_{T, \tau} = \hat{\sigma}^2_{T, \tau}+ \alpha\frac{1}{t_\tau}\sum^{t_{\tau}}_{t=t_{\tau-1}+1} \big(f_t(A_t, X_t) - f_T(A_t, X_t)\big)^2$ and $\alpha > 0$ is a constant. The first term is an efficient weight described above. The second term reflects the deviation between $f_t$ and $f_T$, which would be more accurate because it uses more samples.

\subsection{Main Algorithm}
\label{sec:main_algorithm}
As discussed in Section~\ref{sec:efficiency}, we can minimize the asymptotic variance of the proposed estimator $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e})$ by choosing $w_T$ appropriately. However, to obtain $w_{T,\tau}\xrightarrow{\mathrm{p}} \frac{1}{\sigma^2_\tau}/\sum^M_{\tau'=1}\frac{1}{\sigma^2_{\tau'}}$, which is the optimal weight that minimizes the asymptotic variance, we need a consistent estimator of $\theta_0$, which is what we want to estimate. On the other hand, we have a consistent estimator of $\theta_0$ without using an optimal weight matrix $w_{T,\tau}\xrightarrow{\mathrm{p}} \frac{1}{\sigma^2_\tau}/\sum^M_{\tau'=1}\frac{1}{\sigma^2_{\tau'}}$. Based on these properties, we propose \emph{two-step estimation}. First, using an arbitrary positive definite weight $w^{(0)}_T$, such as the identity matrix, we obtain an initial estimate $\widehat{R}^{\mathrm{BA2IPW}, (1)}_T(\pi^\mathrm{e})$. Then, using $\widehat{R}^{\mathrm{BA2IPW}, (1)}_T(\pi^\mathrm{e})$, we construct $w^{(1)}_{T,\tau}\xrightarrow{\mathrm{p}} \frac{1}{\sigma^2_\tau}/\sum^M_{\tau'=1}\frac{1}{\sigma^2_{\tau'}}$. We can obtain an efficient estimator $\widehat{R}^{\mathrm{BA2IPW}, (2)}_T(\pi^\mathrm{e})$ of $\theta_0$, as discussed. More generally, we consider an algorithm with iteration such that after obtaining $\widehat{R}^{\mathrm{BA2IPW}, (i-1)}_T(\pi^\mathrm{e})$, we estimate $w^{(i-1)}_{T,\tau}\xrightarrow{\mathrm{p}} \frac{1}{\sigma^2_\tau}/\sum^M_{\tau'=1}\frac{1}{\sigma^2_{\tau'}}$ and obtain a next estimator $\widehat{R}^{\mathrm{BA2IPW}, (i)}_T(\pi^\mathrm{e})$ by using $w^{(i-1)}_{T,\tau}$. We refer to this algorithm with $N$-iterations as \emph{$N$-step BA2IPW estimation}. We can use sufficiently large $N$ because, at each iteration, we only calculate the weighted average of the moment conditions using $\widehat{R}^{\mathrm{BA2IPW}, (i)}_T(\pi^\mathrm{e})$, which is not time-consuming. Although the asymptotic properties of the iterated estimator are the same as those of the two-step estimator, we report that the iteration improves the empirical performance in some cases. For $N\geq 2$, we summarize the $N$-Step BA2IPW Estimation in Algorithm~\ref{alg}. We can use any method to construct $f_{t_\tau}$ and $\pi_\tau$ as long as they are consistent for bandit data and satisfy some regularity conditions needed for Theorem~\ref{thm:main}.

\begin{algorithm}[tb]
   \caption{$N$-step efficient BA2IPW estimation}
   \label{alg}
\begin{algorithmic}
   \STATE {\bfseries Input:} $\big\{(X_t, A_t, Y_t)\big\}^{T}_{t=1}$ and  $\big\{f_{t}\big\}^{T-1}_{t=0}$.
   \STATE {\bfseries Initialization:} Let $w^{(0)}_T$ be a positive definite matrix such as the identity matrix.
   \FOR{$i=1$ to $N$}
   \STATE Using $w^{(i-1)}_T$, obtain $\widehat{R}^{\mathrm{BA2IPW}, (i)}_T(\pi^\mathrm{e})  = w^{(i-1)}_T D_T(\pi^\mathrm{e})$.
   \STATE Using $\widehat{R}^{\mathrm{BA2IPW}, (i)}_T(\pi^\mathrm{e})$, construct weights $\{\hat{w}^{(i)}_T\}$.
   \ENDFOR
   \STATE{\bfseries Output:} Estimator $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e}) = \widehat{R}^{\mathrm{BA2IPW}, (N)}_T(\pi^\mathrm{e})$
\end{algorithmic}
\end{algorithm}

\section{Deficient Support Problem}
As an application of BA2IPW, we consider an OPE without Assumption~\ref{asm:overlap_pol}, which assumes that there exists $C_1$ such that $0\leq \frac{\pi^\mathrm{e}(a\mid x)}{p_t(a\mid x)}\leq C_1$. Instead of Assumption~\ref{asm:mds}, we consider a situation in which we are allowed to change the support of actions in each batch. For example, in the first batch, we choose an action from a set $\{1,2,3\}$ with a probability larger than $0$, but we choose an action from a set $\{1,2,4\}$ with a probability larger than $0$ in the second batch. In this case, the probability of choosing the action $4$ is $0$ in the first batch, while the probability of choosing the  action $3$ is $0$ in the second batch. This situation is a common in practice and called deficient support problem \citep{Sachdeva2020}. For this problem, instead of Assumption~\ref{asm:overlap_pol}, we use the following assumption.
\begin{assumption}\label{asm:overlap_pol2}
For $a\in\{1,2,\dots,K\}$, there exist $\tau\in\{1,2,\dots,M\}$ and $C_1$ such that $0\leq \frac{\pi^\mathrm{e}(a\mid x)}{\pi_{\tau}(a\mid x, \Omega_{t_{\tau - 1}})}\leq C_1$.
\end{assumption}
Under this assumption, if $\pi_{\tau}(a\mid x, \Omega_{t_{\tau - 1}})>0$ for at least one batch, we are allowed to use $\pi_{\tau'}(a\mid x, \Omega_{t_{\tau' - 1}})=0$ for $\tau'\neq \tau$. With this assumption, we derive the asymptotic normality in Appendix~\ref{appdx:incomp_support_of_arms}. Thus, our approach provides a new solution to this problem.

\begin{table*}[t]
\begin{center}
\caption{Results of OPE under the RW policy. We highlight in bold the best two estimators in each dataset.}
\medskip
\label{tbl:exp_table1}
\vspace{-0.3cm}
\scalebox{0.70}[0.70]{
\begin{tabular}{l|rr|rr|rr|rr|rr|rr}
\toprule
Datasets &  \multicolumn{2}{c|}{satimage}& \multicolumn{2}{c|}{pendigits}& \multicolumn{2}{c|}{mnist}& \multicolumn{2}{c|}{letter}& \multicolumn{2}{c|}{sensorless}& \multicolumn{2}{c}{connect-4} \\
Metrics &      MSE &      SD &      MSE &      SD &      MSE &      SD &      MSE &      SD &      MSE &      SD &     MSE &     SD \\
\hline
PBA2IPW &  \textbf{0.038} &  0.002 &  \textbf{0.129} &  0.047 &  \textbf{0.173} &  0.104 &  \textbf{0.331} &  0.547 &  \textbf{0.146} &  0.089 &  \textbf{0.021} &  0.021 \\
EBA2IPW &  0.050 &  0.006 &  0.190 &  0.030 &  0.191 &  0.025 &  0.398 &  0.062 &  0.182 &  0.028 &  \textbf{0.025} &  0.022 \\
EBA2IPW' &  0.044 &  0.003 &  0.182 &  0.029 &  \textbf{0.102} &  0.012 &  0.389 &  0.064 &  0.177 &  0.027 &  \textbf{0.025} &  0.024 \\
BAdaIPW &  0.077 &  0.010 &  0.178 &  0.082 &  0.200 &  0.111 &  0.333 &  0.537 &  0.160 &  0.093 &  0.027 &  0.027 \\
AdaDM &  0.141 &  0.010 &  0.493 &  0.034 &  0.434 &  0.036 &  0.476 &  0.023 &  0.413 &  0.031 &  0.142 &  0.024 \\
AIPW &  \textbf{0.032} &  0.001 &  \textbf{0.110} &  0.030 &  0.244 &  0.028 &  \textbf{0.254} &  0.216 &  \textbf{0.128} &  0.064 &  0.055 &  0.022 \\
DM &  0.099 &  0.004 &  0.452 &  0.025 &  0.282 &  0.028 &  0.459 &  0.023 &  0.395 &  0.025 &  0.086 &  0.018 \\
\bottomrule
\end{tabular}
}
\end{center}
\vspace{-0.3cm}
\begin{center}
\caption{Results of OPE under the UCB policy. We highlight in bold the best two estimators in each dataset.}
\medskip
\label{tbl:exp_table2}
\vspace{-0.cm}
\scalebox{0.70}[0.70]{
\begin{tabular}{l|rr|rr|rr|rr|rr|rr}
\toprule
Datasets &  \multicolumn{2}{c|}{satimage}& \multicolumn{2}{c|}{pendigits}& \multicolumn{2}{c|}{mnist}& \multicolumn{2}{c|}{letter}& \multicolumn{2}{c|}{sensorless}& \multicolumn{2}{c}{connect-4} \\
Metrics &      MSE &      SD &      MSE &      SD &      MSE &      SD &      MSE &      SD &      MSE &      SD &     MSE &     SD \\
\hline
PBA2IPW &  0.050 &  0.005 &  0.088 &  0.015 &  0.240 &  0.415 &  \textbf{0.205} &  0.088 &  \textbf{0.162} &  0.057 &  0.032 &  0.032 \\
EBA2IPW &  \textbf{0.014} &  0.000 &  \textbf{0.029} &  0.002 &  0.240 &  0.070 &  0.434 &  0.036 &  0.239 &  0.054 &  0.030 &  0.022 \\
EBA2IPW' &  0.036 &  0.005 &  0.081 &  0.017 &  \textbf{0.142} &  0.030 &  0.422 &  0.032 &  0.202 &  0.038 &  0.037 &  0.033 \\
BAdaIPW &  0.087 &  0.037 &  0.122 &  0.030 &  0.275 &  0.406 &  0.219 &  0.090 &  0.183 &  0.067 &  0.057 &  0.057 \\
AdaDM &  0.076 &  0.002 &  0.230 &  0.008 &  0.372 &  0.020 &  0.451 &  0.018 &  0.327 &  0.021 &  0.071 &  0.023 \\
AIPW &  0.036 &  0.003 &  0.058 &  0.012 &  \textbf{0.136} &  0.007 &  \textbf{0.170} &  0.043 &  \textbf{0.134} &  0.038 &  \textbf{0.022} &  0.012 \\
DM &  \textbf{0.009} &  0.000 &  \textbf{0.049} &  0.001 &  0.161 &  0.008 &  0.371 &  0.019 &  0.214 &  0.013 &  \textbf{0.019} &  0.010 \\
\bottomrule
\end{tabular}
}
\end{center}
\vspace{-0.3cm}
\end{table*}

\section{Estimation of the Behavior Policy}
In the proposed BA2IPW method, we assume that the true behavior policy is known. However, in many real-world applications, the assumption does not hold. To solve this problem, by using an estimator $\hat{g}_{t-1}$ of $\pi_\tau$, which is constructed from $\Omega_{t-1}$ as well as $\hat{f}_{t-1}$, we also propose a \emph{Batch-based Adaptive Doubly Robust} (BADR) as $\hat{R}^{\mathrm{BADR}}_T(\pi^\mathrm{e})  = w^\top_T \widetilde{D}_T(\pi^\mathrm{e})$, where $w_T = (w_{T, 1}\ \cdots\ w_{T, M})^\top$ is an $M$-dimensional vector such that $\sum^M_{\tau=1}w_{T, \tau}=1$, $\tilde{D}_T(\pi^\mathrm{e})=$
\begin{align*}
\begin{pmatrix}
\frac{1}{t_{1}}\sum^{t_{1}}_{t=1} \tilde{\phi}_t(X_t, A_t, Y_t; 1, \hat{f}_{t-1}, \hat{g}_{t-1}, \pi^\mathrm{e}) \\
\frac{1}{t_{2} - {t_{1}}}\sum^{t_2}_{t=t_1+1} \tilde{\phi}_t(X_t, A_t, Y_t; 2, \hat{f}_{t-1}, \hat{g}_{t-1}, \pi^\mathrm{e})\\
\vdots\\
\frac{1}{T-t_{M-1}}\sum^{T}_{t=t_{M-1}+1} \tilde{\phi}_t(X_t, A_t, Y_t; M, \hat{f}_{t-1}, \hat{g}_{t-1}, \pi^\mathrm{e})\\
\end{pmatrix},
\end{align*}
and $\tilde{\phi}_t(x, k, y; \tau, f, g, \pi^\mathrm{e}):=$
\begin{align*}
&\sum^{K}_{a=1}\pi^\mathrm{e}(a\mid x)\Bigg\{\frac{\mathbbm{1}[k=a]\big\{y - f(a, x)\big\}}{g(a\mid x)}+f(a, x)\Bigg\}.
\end{align*}
For the BADR estimator, we show the asymptotic normality as follows. The proof is shown in Appendix~\ref{appdx:badr_dist}.

\begin{theorem}[Asymptotic Distribution of the BADR Estimator]
\label{thm:badr_asymp}
Suppose that (i) $w_T = (w_{T,1}\ \cdots\ w_{T,M})^\top \xrightarrow{\mathrm{p}} w = (w_{1}\ \cdots\ w_{M})^\top$; (ii) $w_{T,\tau} > 0$ and $\sum^M_{\tau=1}w_{T,\tau} = 1$; (iii) for $\alpha\beta=\mathrm{o}_{p}((t-t_{\tau-1})^{-1/2}),\alpha=\mathrm{o}_{p}(1), \beta=\mathrm{o}_{p}(1)$, and each $\tau\in I$, the nuisance estimators satisfy $\|\hat{g}_{t-1}(a\mid X_t) - \pi_{\tau}(a\mid X_t, \Omega_{t_{\tau-1}})\|_{2}=\alpha$, and $\|\hat{f}_{t-1}(a,X_t)-f^*(a,X_t)\|_2=\beta$, where the expectation of the norm is over $X_t$; (iv) there exit constants $C_f$ and $C_g$ such that $|\hat{f}_{t-1}(a, x)| \leq C_f$ and $0 < \left|\frac{\pi^\mathrm{e}(a\mid x)}{\hat{g}_{t-1}(a\mid x)}\right| \leq C_g$ for all $a\in\mathcal{A}$ and $x\in\mathcal{X}$. Then, under Assumptions~\ref{asm:overlap_pol} and \ref{asm:overlap_outcome}, $\sqrt{T}\big(\hat{R}^{\mathrm{BADR}}_T(\pi^\mathrm{e})  - \theta_0\big)  \xrightarrow{\mathrm{d}}\mathcal{N}\big(0, \sigma^2\big)$.
\end{theorem}

\section{Off-Policy Learning}
An important application of OPE is Off-Policy Learning (OPL), which attempts to determine the optimal policy maximizing the expected reward. Let us define the optimal policy $\pi^*$ as $\pi^* = \operatorname*{arg\,max}_{\pi\in \Pi} R(\pi)$, where $\Pi$ is a policy class. By applying each OPE estimator, we estimate the optimal policy as $\hat{\pi} = \operatorname*{arg\,max}_{\pi\in \Pi}\widehat{R}^{\mathrm{BA2IPW}}_T(\pi)$.

\begin{table*}[t]
\begin{center}
\caption{The coverage ratios (CRs) are shown. The left graph shows the results with the RW policy. The right graph shows the results with the UCB policy.}
\medskip
\label{tbl:exp_table3}
\vspace{-0.3cm}
\begin{tabular}{cc}

\begin{minipage}{0.5\hsize}
\begin{center}
\scalebox{0.70}[0.70]{
\begin{tabular}{l|r|r|r|r|r|r}
\toprule
Datasets &  satimage& pendigits& mnist& letter& sensorless& connect-4 \\
\hline
PBA2IPW &  1.00 &  0.96 &  1.00 &  0.85 &  0.94 &  1.00 \\
EBA2IPW &  0.66 &  0.22 &  0.15 &  0.02 &  0.19 &  0.88 \\
EBA2IPW' &  0.67 &  0.27 &  0.47 &  0.02 &  0.19 &  0.89 \\
BAdaIPW &  0.93 &  0.90 &  0.90 &  0.83 &  0.91 &  0.99 \\
\bottomrule
\end{tabular}
}
\end{center}
\end{minipage}

\begin{minipage}{0.5\hsize}
\begin{center}
\scalebox{0.70}[0.70]{
\begin{tabular}{l|r|r|r|r|r|r}
\toprule
Datasets &  satimage& pendigits& mnist& letter& sensorless& connect-4 \\
\hline
PBA2IPW &  1.00 &  1.00 &  1.00 &  0.94 &  0.98 &  1.00 \\
EBA2IPW &  0.92 &  0.74 &  0.01 &  0.00 &  0.02 &  0.78 \\
EBA2IPW' &  0.77 &  0.58 &  0.17 &  0.00 &  0.03 &  0.67 \\
BAdaIPW &  0.91 &  0.88 &  0.91 &  0.89 &  0.82 &  0.93 \\
\bottomrule
\end{tabular}
}
\end{center}
\end{minipage}
\end{tabular}
\end{center}
\vspace{-0.4cm}
\end{table*}

\begin{table*}[t]
\begin{center}
\caption{Results of OPL under the RW policy. We highlight in bold the best two estimators in each dataset.}
\medskip
\label{tbl:exp_table4}
\vspace{-0.3cm}
\scalebox{0.73}[0.73]{
\begin{tabular}{l|rr|rr|rr|rr|rr|rr}
\toprule
Datasets &  \multicolumn{2}{c|}{satimage}& \multicolumn{2}{c|}{pendigits}& \multicolumn{2}{c|}{mnist}& \multicolumn{2}{c|}{letter}& \multicolumn{2}{c|}{sensorless}& \multicolumn{2}{c}{connect-4} \\
Metrics &      RWD &      SD &      RWD &      SD &      RWD &      SD &      RWD &      SD &      RWD &      SD &     RWD &     SD \\
\hline
PBA2IPW &  0.812 &  0.020 &  0.690 &  0.072 &  0.493 &  0.298 &  \textbf{0.172} &  0.090 &  0.264 &  0.131 &  0.665 &  0.030 \\
EBA2IPW &  0.813 &  0.022 &  \textbf{0.717} &  0.063 &  \textbf{0.519} &  0.313 &  0.135 &  0.073 &  \textbf{0.289} &  0.126 &  \textbf{0.679} &  0.024 \\
BAdaIPW &  \textbf{0.815} &  0.023 &  \textbf{0.697} &  0.089 &  \textbf{0.515} &  0.312 &  0.150 &  0.078 &  \textbf{0.313} &  0.145 &  0.677 &  0.024 \\
AdaDM &  0.777 &  0.033 &  0.478 &  0.076 &  0.191 &  0.151 &  0.046 &  0.026 &  0.190 &  0.091 &  0.654 &  0.023 \\
AIPW &  \textbf{0.819} &  0.020 &  0.698 &  0.062 &  \textbf{0.515} &  0.312 &  \textbf{0.154} &  0.087 &  0.287 &  0.141 &  \textbf{0.678} &  0.024 \\
DM &  0.791 &  0.034 &  0.544 &  0.071 &  0.247 &  0.174 &  0.057 &  0.036 &  0.210 &  0.090 &  0.654 &  0.023 \\
\bottomrule
\end{tabular}
}
\end{center}
\vspace{-0.5cm}
\end{table*}


\section{Experiments}
Using benchmark datasets and real-world logged data, we demonstrate the effectiveness of the BA2IPW estimator with an equal weight (PBA2IPW) and efficient weight using variance estimators $\hat{\sigma}^2_{T,\tau}$ (EBA2IPW) and $\tilde{\sigma}^2_{T,\tau}$ (EBA2IPW'), and BAdaIPW estimator with an equal weight (BAdaIPW). Note that although the forms of the several estimators are the same as the existing studies, the theoretical guarantees are different.

\subsection{Experiments with Benchmark Dataset}
Following \citet{dudik2011doubly} and \citet{Chow2018}, we evaluate the proposed estimators using classification datasets by transforming them into contextual bandit data. From the LIBSVM repository, we use the satimage, pendigits, mnist, letter, sensorless, and connect-4 datasets \footnote{\url{https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/}}. For a batched update behavior policy, we use the random walk (RW) and LinearUCB (UCB) \citep{sutton1998reinforcement,li2010contextual,Wei2011} policies. When using the RW policy, we first decide the probability of choosing an action from uniform distribution. Then, we add a noise $0.01\mathcal{N}(0,1)$ at each batch, i.e., the policy is a random walk. At each batch, we standardize the values of random walk to be probability, i.e., all values are positive and the sum is $1$. When using the UCB policy, we choose an estimated best arm $\hat{A}_t$ firstly. Then, we create an adaptive policy. Then, we construct a behavior policy as a policy that chooses $\hat{A}_t$ with probability $0.8$ and the other arms with equal probability. While the probability of choosing an action of the UCB policy converges, that of the RW policy does not converge.

For each dataset, we compare the performances of the following estimators of policy value: PBA2IPW, EBA2IPW, EBA2IPW', BAdaIPW, Adaptive DM estimator (AdaDM) defined as $\frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\pi^\mathrm{e}(a\mid X_t)\hat{f}_{t}(a, X_t)$, an AIPW defined as $\frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\frac{\pi^\mathrm{e}(a\mid X_t)\mathbbm{1}[A_t=a]\big(Y_t - \hat{f}_{T}(a, X_t)\big)}{\pi_t(a\mid X_t, \Omega_{t-1})} + \frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\pi^\mathrm{e}(a\mid X_t)\hat{f}_{T}(a, X_t)$, and DM estimator (DM) defined as $\frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\pi^\mathrm{e}(a\mid X_t)\hat{f}_{T}(a, X_t)$. When estimating $f^*$, we use the Nadaraya-Watson regression (NW) estimator \citep{yang2002}.

\paragraph{MSEs:} To construct an evaluation policy, we create a deterministic policy $\pi^d$ by training a logistic regression classifier on historical data and set the output as $\pi^d$. Through experiments, the behavior policy $\pi_t$ is assumed to be known. More details are in Appendix~\ref{appdx:det_exp}. Let us construct the evaluation policy $\pi^\mathrm{e}$ as a mixture of $\pi^d$ and the uniform random policy $\pi^u$, defined as $\pi^\mathrm{e} = 0.9\pi^d+0.1\pi^u$. We construct the evaluation policy $\pi^\mathrm{e}$ as a mixture of $\pi^d$ and the uniform random policy $\pi^u$ defined as $\pi^\mathrm{e} = 0.9\pi^d+0.1\pi^u$. We compare the MSEs of six estimators, the PBA2IPW, EBA2IPW, EBA2IPW', BAdaIPW, AdaDM, AIPW, and DM estimators. For estimating the weights of EBA2IPW and EBA2IPW' estimators, we iteratively estimate the wights and the value $\theta_0$ $10$ times. In each experiment, we have historical data with a sample size $T=1500$. When estimating $f^*$, we use the Nadaraya-Watson regression (NW) estimator \citep{yang2002}.

The resulting MSEs and their standard deviations (SDs) over $100$ replications of each experiment are shown in Tables~\ref{tbl:exp_table1} and \ref{tbl:exp_table2}. In many cases, the proposed methods show the preferable the existing methods. When using the RW policy, the policy does not converge to a time-invariant policy. Therefore, the proposed method is theoretically preferable for the situation. More importantly, we can construct confidence intervals from the proposed methods, but cannot construct it from the A2IPW estimator. On the other hand, when using the UCB policy, the policy converges to a time-invariant policy, but the proposed methods still show higher performance for various datasets in some cases. We consider that this result is based on the fact that, even though the policy approaches to a time-invariant policy, the update is only allowed in batch, and does not converge sufficiently. In Appendix~\ref{appdx:det_exp}, we show the additional results.

\paragraph{Coverage Ratio of Confidence Interval:} In Table~\ref{tbl:exp_table3}, we show the coverage ratio of the confidence intervals derived in the previous experiments together with the MSEs. The coverage ratio of the confidence interval is a percentage at which it covers the true value $\theta_0$ in the confidence interval. For the $100$ trials of the previous experiment, we calculate the coverage ratio (CR) of $95\%$ confidence interval, which is constructed as $\left[\hat{\theta}_T - 1.96\sqrt{\frac{\hat{\sigma}^2}{1500}}, \hat{\theta}_T + 1.96\sqrt{\frac{\hat{\sigma}^2}{1500}}\right]$, where $\hat{\theta}_T$ is an estimator of $\theta_0$ and $\hat{\sigma}^2$ is its estimated asymptotic variance. In the results, the PBA2IPW and BAdaIPW estimator shows CR close to $0.95$ in many cases. The BAdaIPW estimator does not require an estimator of $f^*$. Therefore, compared with the other estimators, it shows more preferable performances. EBA2IPW estimator requires several variance estimators, and their estimation error worsen the results compared with the PBA2IPW and BAdaIPW estimators.

\paragraph{OPL:}
In the experiments of OPL, we compare the performances of estimated policy maximizing expected reward obtained from the PBA2IPW, EBA2IPW, BAdaIPW, AdaDM, AIPW, and DM estimators. We conducted $5$ trials for each experiment with $T=1500$. The resulting expected rewards over the evaluation data (RWDs) and the SDs are shown in Table~\ref{tbl:exp_table4}, where we highlight in bold the best two estimators.

\subsection{Experiments with Real-World Data}
We apply our estimators to evaluate a policy using the real-world dataset in CyberAgent Inc., which is the second-largest Japanese advertisement company with about $7$ billion USD market capitalization (as of August 2020). This company simultaneously runs Thompson sampling and uniformly random sampling to determine the design of advertisements. The Thompson sampling updates the parameter every $30$ minute; therefore, there are batches with samples obtained during the $30$ minute. We use the logged data produced by the algorithms to confirm the empirical performance of the proposed estimators. To check the performance, we calculate the estimation error between the estimates of the value of the uniformly random sampling policy estimated from the dataset obtained by the Thompson sampling and the observed average reward of the uniformly random policy. More details are shown in Appendix~\ref{appdx:det_exp2}. We apply the PBA2IPW, EBA2IPW, EBA2IPW', BAdaIPW, and AdaDM estimators. The results are shown in Table~\ref{tbl:exp_table5} and Figure~\ref{fig1:box_plot}. We show the Bias, MSE, and averaged confidence intervals. While the PBA2IPW and BAdaIPW estimators suffer the high variance, the EBA2IPW, EBA2IPW', and AdaDM estimators show the preferable performances. Although the AdaDM estimator also shows the effectiveness for this dataset, the BA2IPW estimator is theoretically more robust because it is consistent even if $\hat{f}_t$ does not converge to $f^*$.

\begin{table}[h]
\begin{center}
\caption{Results of the CyberAgent dataset. }
\medskip
\label{tbl:exp_table5}
\scalebox{0.73}[0.73]{
\begin{tabular}{lrrrrrrrrrrrrrrrr}
\toprule
{} &      Bias &      MSE &   $\frac{1.96\times \mathrm{variance}}{\mathrm{sample}\ \mathrm{size}}$\\
\midrule
PBA2IPW &  -0.01766 &  8.25830e-02 &  8.25830e-02  \\
EBA2IPW &  0.00019 &  4.52574e-05 &  3.58261e-03  \\
EBA2IPW' &  0.00095 &  4.03248e-05 &  2.92494e-03  \\
BAdaIPW &  0.14650 &  5.70524e-02 &  7.68610e-02  \\
AdaDM &  0.00714 &  2.26107e-04 &  2.23604e-05  \\
\bottomrule
\end{tabular}
}
\end{center}
\vspace{-0.3cm}
\end{table}

\begin{figure}[h]
\begin{center}
 \includegraphics[width=55mm]{fig_boxplot.pdf}
\end{center}
\caption{Box plot of the CyberAgent dataset.}
\label{fig1:box_plot}
\vspace{-0.3cm}
\end{figure}

\section{Conclusion}
This study presented solutions for causal inference from dependent samples obtained via a batch-based bandit algorithm. By using the asymptotic property in batch, we applied the CLT for an MDS to obtain the asymptotic normality without requiring the convergence of the assignment probability. Additionally, we showed that the proposed method is applicable when the support of arms is incomplete, which is a notorious problem in OPE. In experiments, the proposed batch-based estimators showed theoretically expected performances for benchmark and real-world datasets.

\bibliographystyle{icml2019}
\bibliography{arxiv}

\clearpage

\onecolumn