Statistical Inference of the Value Function for Reinforcement Learning in Infinite Horizon Settings

Introduction

Reinforcement learning (RL) is a general technique that allows an agent to learn and interact with an environment. A policy defines the agent’s way of behaving. It maps the states of environments to a set of actions to be chosen from. RL algorithms have made tremendous achievements and found extensive applications in video games (Silver et al. 2016), robotics (Kormushev et al. 2013), bidding (Jin et al. 2018), ridesharing (Xu et al. 2018), etc. In particular, a number of RL methods have been proposed in precision medicine, to derive an optimal policy as a set of sequential treatment decision rules that optimize patients’ clinical outcomes over a fixed period of time (finite horizon). References include Murphy 2003; Zhang et al. 2013; Zhao et al. 2015; Shi et al. 2018a; Shi et al. 2018b; Zhang et al. 2018, to name a few.

Mobile health (or mHealth) technology has recently emerged due to the use of mobile devices such as mobile phones, tablet computers or wearable devices in health care. It allows health-care providers to communicate with patients and manage their illness in real time. It also collects rich longitudinal data (e.g., through mobile health apps) that can be used to estimate the optimal policy. Data from mHealth applications differ from those in finite horizon settings in that the number of treatment decision points for each patient is not necessarily fixed (infinite horizon) while the total number of patients could be limited. Take the OhioT1DM dataset (Marling and Bunescu 2018) as an example. It contains data for six patients with type 1 diabetes. For all patients, their continuous glucose monitoring (CGM) blood glucose levels, insulin doses including bolus and basal rates, self-reported times of meals and exercises are continually measured and recorded for eight weeks. Developing an optimal policy as functions of these time-varying covariates could potentially assist these patients in improving their health status.

In this paper, we focus on the infinite horizon setting where the data generating process is modeled by a Markov decision process (Puterman 1994, MDP,). Specifically, at each time point, the agent selects an action based on the observed state. The system responds by giving the decision maker a corresponding outcome and moving into a new state in the next time step. This model is generally applicable to sequential decision making, including applications from mHealth, games, robotics, ridesharing, etc. After a policy is being proposed, it is important to examine its benefit prior to recommending it for practical use. The goodness of a policy is quantified by its (state) value function, corresponding to the discounted cumulative reward that the agent receives on average, starting from some initial state. The inference of the value function helps a decision maker to evaluate the impact of implementing a policy when the environment is in a certain state. In some applications, it is also important to evaluate the integrated value of a policy aggregated over different initial states. For example, in medical studies, one might wish to know the mean outcome of patients in the population. The integrated value could thus be used as a criterion for comparing different policies.

In statistics literature, a few methods have been proposed to estimate the optimal policy in infinite horizons. Ertefaie and Strawderman 2018 proposed a variant of gradient Q-learning method. Luckett et al. 2019 proposed a V-learning to directly search the optimal policy among a restricted class of policies. Inference of the value function under a generic (data-dependent) policy has not been studied in these papers. In the computer science literature, Thomas et al. 2015 and Jiang and Li 2016 proposed (augmented) inverse propensity-score weighted ((A)IPW) estimators for the the value function in infinite horizons and derived their associated CIs. However, these methods are not suitable for settings where only a limited number of trajectories (e.g., plays of a game or patients in medical studies) are available, since (A)IPW estimators become increasingly unstable as the number of decision points diverges to infinity. Recently, Kallus and Uehara 2019 proposed a double reinforcement learning (DRL) method that achieves consistent estimation of the value under a fixed policy even with limited number of trajectories. Their method computes a Q-function and a marginalized density ratio. Learning the density ratio is challenging in general and it remains difficult to investigate the goodness-of-fit of the estimated density ratio in practice.

The focus of this paper is to construct confidence intervals (CIs) for a (possibly data-dependent) policy’s value function at a given state as well as its integrated value with respect to a given reference distribution. Our proposed CI is derived by estimating the state-action value function (Q-function) under the target policy. Similar to the value, the Q-function measures the discounted cumulative reward that the agent receives on average, starting from some initial state-action pair. We use series/sieve method to approximate the Q-function based on $L$ basis functions, where $L$ grows with the total number of observations. The advances of our proposed method are summarized as follows. First, the proposed inference method is generally applicable. Specifically, it can be applied to any fixed policy (either deterministic or random) and any data-dependent policy whose value converges at a certain rate. The latter includes policies estimated by general Q-learning type algorithms that learns an optimal Q-function from the observed data, such as gradient Q-learning (Maei et al. 2010; Ertefaie and Strawderman 2018), fitted Q-iteration (Ernst et al. 2005; Riedmiller 2005, see for example,), etc. See Section 3.2.4 for detailed illustrations.

Second, when applied to data-dependent policies, our method is valid in nonregular cases where the optimal policy is not uniquely defined. Inference without requiring the uniqueness of the optimal policy is extremely challenging even in the simpler finite-horizon settings (Luedtke and van der Laan 2016, see the related discussions in). The major challenge lies in that the estimated policy may not stabilize as sample size grows, making the variance of the value estimator difficult to estimate (see Section 3.2.1 for details). We achieve valid inference by proposing a SequentiAl Value Evaluation (SAVE) method that splits the data into several blocks and recursively update the estimated policy and its value estimator. It is worth mentioning that the data-splitting rule cannot be arbitrarily determined since the observations are time dependent in infinite horizon settings (see Section 3.2.2 for details).

Third, our CI is valid as long as either the number of trajectories $n$ in the data, or the number of decision points $T$ per trajectory diverges to infinity. It can thus be applied to a wide variety of real applications in infinite horizons ranging from the Framingham heart study (Tsao and Vasan 2015) with over two thousand patients to the OhioT1DM dataset that contains eight weeks’ worth of data for six people. We also allow both $n$ and $T$ to approach infinity, which is the case in applications from video games. In contrast, CIs proposed by Thomas et al. 2015 and Jiang and Li 2016 require $n$ to grow to infinity to achieve nominal coverage.

Lastly, we consider both off-policy and on-policy learning methods. In off-policy settings, CIs are derived based on historical data collected by a potentially different behavior policy. Off-policy evaluation is critical in situations where running the target policy could be expensive, risky or unethical. In on-policy settings, the estimated policy is recursively updated as batches of new observations arrive. To our knowledge, this is the first work on statistical inference of a data-dependent policy in on-policy settings in sequential decision making with infinite horizons.

To study the asymptotic properties of our proposed CI, we focus on tensor-product spline and wavelet series estimators. Our technical contributions are described as follows. First, we introduce a bidirectional-asymptotic framework that allows either $n$ or $T$ to approach infinity. Our major technical contribution is to derive a nonasymptotic error bound for the spectral norm of sums of mean zero random matrices formed by the data transactions from MDP as a function of $n$, $T$ and $L$ (see e.g., Lemma 3). This result is important in studying the limiting distribution of series estimators under such a theoretical framework.

Second, for policies that are estimated by Q-learning type algorithms such as the greedy gradient Q-learning, fitted Q-iteration and deep Q-network (Mnih et al. 2015), we relate the convergence rate of their values to the prediction error of the corresponding estimated Q-functions. We show in Theorems 3 and 4 that the values can converge at faster rates than the estimated Q-functions under certain margin type conditions on the optimal Q-function. To the best of our knowledge, these findings have not been discovered in the reinforcement-learning literature. Our theorems form a basis for researchers to study the value properties of Q-learning type algorithms. Moreover, our theoretical results are consistent with findings in point treatment studies where there is only one single decision point (Qian and Murphy 2011; Luedtke and van der Laan 2016, see e.g.,). However, the derivation of Theorems 3 and 4 is more involved since the value function in our settings is an infinite series involving both immediate and future rewards.

Third, when these basis functions are used, we mathematically characterize the approximation error of the Q-function as a function of $L$, the dimension of the state variables, and the smoothness of the Markov transition function and the conditional mean of the immediate reward as a function of the state-action pair. This offers some guidance to practitioners on the choice of the number of basis functions $L$, when some prior knowledge on the degree of smoothness of the aforementioned functions are available.

The rest of the paper is organized as follows. We introduce the model setup in Section 2. In Sections 3 and 4, we present the proposed off-policy and on-policy evaluation methods, respectively. Simulation studies are conducted to evaluate the empirical performance of the proposed inference methods in Section 5. We apply the proposed inference method to the OhioT1DM dataset in Section 6, Finally, we conclude our paper by a discussion section.

Optimal policy in infinite-horizon settings

for some transition function $\mathcal{P}$. Here, $\mathcal{P}$ defines the next state distribution conditional on the current state-action pair. Moreover, suppose the following conditional mean independence assumption (CMIA) holds

for some reward function $r$. By MA, CMIA automatically holds when $Y_{0,t}$ is a deterministic function of $X_{0,t},A_{0,t}$ and $X_{0,t+1}$ that measures the system’s status at time $t+1$. The latter is satisfied in our real data application (see Section 6 for details) and is commonly assumed in the reinforcement learning literature. CMIA is thus weaker than this condition. MA and CMIA are important to guarantee the existence of an optimal policy (see (2.1)) and derive the bidirectional-asymptotic theory of the proposed CI (see the discussions below Theorem 1). We assume both assumptions hold throughout this paper.

Similar to Theorem 6.2.12 of Puterman 1994, we can show under the given conditions that there exists at least one optimal policy $\pi^{\tiny{opt}}$ that satisfies

To better understand $\pi^{\tiny{opt}}$, we introduce the state-action function (Q-function) under a policy $\pi$ as

Let $Q^{\tiny{opt}}$ denote the optimal Q-function, i.e, $Q^{\tiny{opt}}(\cdot,\cdot)=\sup_{\pi}Q(\pi;\cdot,\cdot)$. It can be shown that $\pi^{\tiny{opt}}$ satisfies

where $\sargmax$ denotes the smallest maximizer when the argmax is not unique. Such a deterministic optimal policy may be appealing in medical studies. For example, in optimal dose studies, it is preferred to assign each patient the smallest optimal dose level to avoid toxicity.

Off-policy evaluation

Let $n$ denote the number of trajectories in the dataset. For the $i$-th trajectory, let $\{A_{i,t}\}_{t\geq 0}$, $\{X_{i,t}\}_{t\geq 0}$ and $\{Y_{i,t}\}_{t\geq 0}$ denote the sequence of actions, states and rewards, respectively. It is worth mentioning that the time points are not necessarily homogeneous across different trajectories. Suppose the data are generated according to a fixed policy $b(\cdot|\cdot)$, better known as the behavior policy such that

are i.i.d copies of $\{(X_{0,t},A_{0,t},Y_{0,t})\}_{t\geq 0}$. The observed data can thus be summarized as $\{(X_{i,t},A_{i,t},Y_{i,t},X_{i,t+1})\}_{0\leq t<T_{i},1\leq i\leq n}$, where $T_{i}$ is the termination time of the $i$-th trajectory. The goal of off-policy evaluation is to learn the value under a target policy $\pi(\cdot|\cdot)$, possibly different from $b(\cdot|\cdot)$.

Luckett et al. 2019 showed that the value function satisfies

Based on (3.5), they directly modelled the value function, constructed an estimator for the integrated value under their estimated optimal policy and proved that it is asymptotically normal (Luckett et al. 2019, see Theorem 4.3,).

Following their procedure, for a fixed policy $\pi$, one might estimate $V(\pi;\cdot)$ nonparametrically and construct the CI using the resulting estimate. However, such an approach might not be appropriate for polices that are discontinuous functions of the covariates. To better illustrate this, notice that $V(\pi,\cdot)$ satisfies the following Bellman equation

When $\mathcal{P}$ satisfies certain smoothness conditions (see Condition A1 below), we have

for any $\pi$. Suppose $r(\cdot,a)$ is continuous for any $a\in\mathcal{A}$. Then $C(\pi;x,a)$ is continuous in $x$ for any $\pi$ and $a$. When $\pi$ is a non-continuous function of $x$, it follows from (3.6) that $V(\pi;\cdot)$ is not continuous either. However, many nonparametric methods, such as kernel smoothers and series estimation, require the underlying function to possess certain degree of smoothness in order to achieve estimation consistency. Notice that any non-constant deterministic policy has jumps and is not continuous at certain points (such as the optimal policy $\pi_{0}$ given in (2)). This poses significant challenges in performing inference to these policies.

To allow valid inference for both deterministic and random policies, we consider modelling the Q-function. Under CMIA, we have

As a result, the Q-function satisfies the following Bellman equation

Similar to (3.7), we can show the second term on the right-hand-side (RHS) of (3.7) is a smooth function of $x$ for any $\pi$ and $a$. When $r(\cdot,a)$ is smooth, so is $Q(\pi,\cdot,a)$. To formally establish these results, we introduce the notion of $p$-smoothness (also known as Hölder smoothness with exponent $p$) below.

Here, $x_{j}$ denotes the $j$-th element of $x$. For any $p>0$, let $\lfloor p\rfloor$ denote the largest integer that is smaller than $p$. Define the class of $p$-smooth functions as follows:

When $0<p\leq 1$, we have $\lfloor p\rfloor=0$. It is equivalent to require $h$ to satisfy $\sup_{x,y}|h(x)-h(y)|/\|x-y\|_{2}^{p}\leq c$. The notion of $p$-smoothness is thus reduced to the Hölder continuity.

Under A1, there exists some constant $c^{\prime}>0$ such that $Q(\pi;\cdot,a)\in\Lambda(p,c^{\prime})$ for any policy $\pi$ and $a\in\mathcal{A}$.

Lemma 1 implies the Q-function has bounded derivatives up to order $\lfloor p\rfloor$. This motivates us to first estimate the Q-function and then derive the corresponding value estimators based on the relation $V(\pi;x)=\sum_{a\in\mathcal{A}}\pi(a|x)Q(\pi;x,a)$. By the Bellman equation (3.8), we can show the Q-function satisfies

The above equation forms a basis of our methods to learn $Q(\pi;\cdot,\cdot)$ (see details in the next section). In contrast to Equation (3.5), the sampling ratio $\pi(a|x)/b(a|x)$ does not appear in (3.9). This is because $A_{i,t}$ is the only sampling action and no further actions are involved in (3.9). As a result, our method does not require correct specification of the behavior policy. Nor do we need to estimate it from the observed dataset. This is another advantage of modelling the Q-function over the value.

1.2 Method

We describe our procedure in this section. We propose to approximate $Q(\pi;\cdot,\cdot)$ based on linear sieves, which takes the form

where $\Phi_{L}(\cdot)=\{\phi_{L,1}(\cdot),\cdots,\phi_{L,L}(\cdot)\}^{\top}$ is a vector consisting of $L$ sieve basis functions, such as splines or wavelet bases (see for example, Huang 1998, for choices of basis functions). We allow $L$ to grow with the sample size to reduce the bias of the resulting estimates. Under certain mild conditions, there exist some $\{\beta_{\pi,a}^{*}\}_{a\in\mathcal{A}}$ that satisfy

for any $a^{\prime}\in\mathcal{A}$. Recall that $\mathcal{A}=\{0,1,\dots,m-1\}$. Define $\bm{\beta}^{*}_{\pi}=(\beta_{\pi,1}^{*T},\cdots,\beta_{\pi,m}^{*T})^{\top}$,

Let $\widehat{\bm{\beta}}_{\pi}=(\widehat{\beta}_{\pi,1}^{\top},\cdots,\widehat{\beta}_{\pi,m}^{\top})^{\top}$, we propose to estimate $V(\pi;x)$ by

where $z_{\alpha}$ denotes the upper $\alpha$-th quantile of a standard normal distribution, and

1.3 Theory

Following the behavior policy $b(\cdot|\cdot)$, the set of variables $\{X_{0,t}\}_{t\geq 0}$ forms a time-homogeneous Markov chain. Its transition kernel $\mathcal{P}_{X}$ is given by

(A2.) The Markov chain $\{X_{0,t}\}_{t\geq 0}$ has an unique invariant distribution with some density function $\mu(\cdot)$. The density functions $\mu$ and $\nu_{0}$ are uniformly bounded away from $0$ and $\infty$.

(ii) The Markov chain $\{X_{0,t}\}_{t\geq 0}$ is geometrically ergodic.

We make a few remarks. First, we do not require the limiting density function $\mu$ to be equal to the initial state density $\nu_{0}$.

Finally, when $\nu_{0}=\mu$, $\{X_{0,t}\}_{t\geq 0}$ is stationary. Under Condition A3(ii), it follows from Theorem 3.7 of Bradley 2005 that $\{X_{0,t}\}_{t\geq 0}$ is exponentially $\beta$-mixing (see the proof of Lemma 3 for details). When $T\to\infty$, A3(ii) enables us to derive matrix concentration inequalities for $\widehat{\bm{\Sigma}}_{\pi}$. This together with A3(i) implies that $\widehat{\bm{\Sigma}}_{\pi}$ is invertible, with probability approaching $1$ (wpa1). We remark that A3(ii) is not needed when $T$ is bounded.

By MA, CMIA and (3.9), the leading term on the RHS of (3.14) forms a mean-zero martingale (details can be found in Section E.5). As either $n$ or $T$ grows to infinity, the asymptotic normality follows from the martingale central limit theorem.

2 Inference of the value under an (estimated) optimal policy

For simplicity, we assume $T_{1}=T_{2}=\cdots=T_{n}=T$ throughout this section. Consider an estimated policy $\widehat{\pi}$, computed based on the data $\{(X_{i,t},A_{i,t},Y_{i,t},X_{i,t+1})\}_{0\leq t<T,1\leq i\leq n}$. The integrated value under $\widehat{\pi}$ is given by

We will require the value of $\widehat{\pi}$ to converge to some fixed policy $\pi^{*}$ (possibly different from $\pi^{\tiny{opt}}$), i.e,

We begin by outlining the challenge of obtaining inference in the nonregular cases. Suppose $\pi^{*}\in\Pi^{\tiny{opt}}$. When the optimal policy is not unique, $\widehat{\pi}$ might not converge to a fixed policy, despite that its value converges (see (3.15)). To better illustrate this, suppose $\widehat{\pi}$ is computed by some Q-learning type algorithms, i.e,

To allow for valid inference, we use a sequential value evaluation procedure to construct the CI. That is, we propose sequentially estimating the optimal policy and evaluating its value using different data subsets. This allows us to treat the estimated optimal policy as known conditional on past observations (see Equation (E.37) in Appendix E.2). The martingale CLT can thus be applied to obtain the limiting distribution for our estimator (see (E.39) and the related discussions). We detail our procedure in the next section.

2.2 SAVE for the value under an (estimated) optimal policy

We begin by dividing $\mathcal{I}_{0}=\{(i,t):1\leq i\leq n,0\leq t<T\}$ into $K$ non-overlapping subsets, denoted by $\mathcal{I}_{1},\mathcal{I}_{2},\cdots,\mathcal{I}_{K}$. At the $k$-th step, we use the sub-dataset

As commented in the introduction, the data-splitting rule cannot be arbitrary. For any of the two tuples $(i_{1},t_{1})$ and $(i_{2},t_{2})$, define an order $(i_{2},t_{2})\succ(i_{1},t_{1})$ if either $t_{2}>t_{1}$ or $i_{2}>i_{1}$. For any $(i_{2},t_{2})\in\mathcal{I}_{k+1}$, we require the following:

It remains to specify $\mathcal{I}_{1},\mathcal{I}_{2},\cdots,\mathcal{I}_{K}$ that satisfy (3.20). Consider some positive integers $n_{\min}\leq n$, $T_{\min}\leq T$. Assume $n$ and $T$ are divisible by $n_{\min}$ and $T_{\min}$, respectively. Let $K_{n}=n/n_{\min}$ and $K_{T}=T/T_{\min}$. We set $K=K_{n}K_{T}$. For any $1\leq k_{n}\leq K_{n}$, $1\leq k_{T}\leq K_{T}$, define a set $\mathcal{I}(k_{n},k_{T})$ by

Thus, each block $\mathcal{I}(k_{n},k_{T})$ contains data from $n_{\min}$ trajectories with $T_{\min}$ decision time points. Below, we introduce two special examples.

When only a few trajectories are available, we may set $n_{\min}=n$. Then, the blocks are constructed according to the times that decisions were being made.

When each trajectory contains a very short time period, we may set $T_{\min}=T$. Then, the observations are divided according to the trajectories they belong to.

Based on this order, we set $\mathcal{I}_{k}=\mathcal{I}(n(k),T(k))$ where $n(k)$ and $T(k)$ are the unique positive integers that satisfy $k=n(k)+(T(k)-1)K_{n}$. For any $k_{2}>k_{1}$, we have either $n(k_{2})>n(k_{1})$ or $T(k_{2})>T(k_{1})$. Thus, the proposed data-splitting rule guarantees (3.20) holds for any $k$.

In Theorem 2 below, we establish the validity of our CI in (3.21). It relies on Condition A3* and A4. A3* is very similar to A3 and we present the detailed definition in Appendix A to save space.

We provide a sketch for the proof of Theorem 2 in Appendix E.2.

2.3 Convergence of the value under an estimated optimal policy

For any $\mathcal{I}\subseteq\mathcal{I}_{0}$, we use $\widehat{\pi}_{\mathcal{I}}$ to denote an estimated optimal policy based on observations in $\mathcal{I}$. Let $\widehat{Q}_{\mathcal{I}}(\cdot,\cdot)$ denote some consistent estimator for $Q^{\tiny{opt}}(\cdot,\cdot)$ and $\widehat{\pi}_{\mathcal{I}}$ denote the greedy policy with respect to $\widehat{Q}_{\mathcal{I}}(\cdot,\cdot)$ (see Equation (3.2.1)).

(A5) Assume there exist some constants $\alpha,\delta_{0}>0$ such that

where $\lambda$ denotes the Lebesgue measure, the big-$O$ terms are uniform in $0<\varepsilon\leq\delta_{0}$, and $\max_{a^{\prime}\in\mathcal{A}-\argmax_{a}Q^{\tiny{opt}}(x,a)}Q^{\tiny{opt}}(x,a^{\prime})=-\infty$ if the set $\mathcal{A}-\argmax_{a}Q^{\tiny{opt}}(x,a)=\emptyset$.

For each $x$, the quantity $\max_{a}Q^{\tiny{opt}}(x,a)-\max_{a^{\prime}\in\mathcal{A}-\argmax_{a}Q^{\tiny{opt}}(x,a)}Q^{\tiny{opt}}(x,a^{\prime})$ measures the difference in value between $\pi^{\tiny{opt}}$ and the policy that assigns the best suboptimal treatment(s) at the first decision point and follows $\pi^{\tiny{opt}}$ subsequently. In point treatment studies, Qian and Murphy 2011 imposed a similar condition (Qian and Murphy 2011, see Equation (3.3),) to derive sharp convergence rate for the value under an estimated optimal individualized treatment regime. Here, we generalize their condition in infinite-horizon settings. A5 is also closely related to the margin condition commonly used to bound the excess misclassification error (Tsybakov 2004; Audibert and Tsybakov 2007).

The margin-type condition is mild. In Appendix A.3, we present detailed examples and show the condition holds under these examples. The following theorems summarize our results.

Assume A1, (3.22) and (3.23) hold. Suppose the following event occurs with probability at least $1-O(|\mathcal{I}|^{-\kappa})$ for any finite $\kappa>0$,

In Theorem 3, we require the estimated Q-function to satisfy certain uniform convergence rate. In Theorem 4 below, we relax this condition by assuming that the integrated loss converges to zero at certain rate.

It can be seen from Theorems 3 and 4 that the integrated value converges faster compared to the Q-function. We provide a sketch for the proofs of both theorems in Appendix E.3.

2.4 Applications

In this section, we provide several examples to illustrate the convergence rate of $\widehat{Q}_{\mathcal{I}}$. The proposed methods can be applied to evaluating the values under these estimated policies. The algorithm in Example 1 requires to impose a linear model assumption for the optimal Q-function. The algorithm in Example 2 allows more general nonlinear and nonparametric models for the optimal Q-function.

Suppose we model $Q^{\tiny{opt}}(x,a)$ by linear sieves $\Phi_{L}^{\top}(x)\theta_{a}$. Then we can compute $\{\widehat{\theta}_{a,\mathcal{I}}\}_{a\in\mathcal{A}}$ by minimizing the following projected Bellman error:

where $\delta_{i,t}(\{\theta_{a}\}_{a\in\mathcal{A}})=Y_{i,t}+\gamma\max_{a^{\prime}\in\mathcal{A}}\Phi_{L}^{\top}(X_{i,t+1})\theta_{a^{\prime}}-\Phi_{L}^{\top}(X_{i,t})\theta_{A_{i,t}}$. The above loss is non-smooth and non-convex as a function of $\{\theta_{a}\}_{a\in\mathcal{A}}$. The estimator $\{\widehat{\theta}_{a,\mathcal{I}}\}_{a\in\mathcal{A}}$ can be computed based on the greedy gradient Q-learning algorithm.

In fitted $Q$-iteration (FQI), the optimal Q-function is approximated by some nonparametric models $Q(\cdot,\cdot,\theta)$ indexed by $\theta$. The parameter $\theta$ is iteratively updated by

Extensions to on-policy evaluation

We now extend our methodology in Section 3 to on-policy settings. The proposed CI is similar to that presented in Section 3.2.2 and applies to any reinforcement learning algorithms that iteratively update the estimated policy based on batches of observations. Let $\{T(k)\}_{k\geq 1}$ be a monotonically increasing sequence that diverges to infinity. At the $k$-th iteration, define $\bar{\mathcal{I}}_{k}=\{(i,t):1\leq i\leq n,0\leq t<\sum_{j=1}^{k}T(j)\}$. The data observed so far can be summarized as $\{(X_{i,t},A_{i,t},Y_{i,t},X_{i,t+1})\}_{1\leq i\leq n,0\leq t<\sum_{j=1}^{k}T(j)}$. We compute the estimated policy $\widehat{\pi}_{\bar{\mathcal{I}}_{k}}$ based on these data. Then we determine the behavior policy $\widehat{b}_{\bar{\mathcal{I}}_{k}}$ as a function of $\widehat{\pi}_{\bar{\mathcal{I}}_{k}}$ and generate new observations

according to $\widehat{b}_{\bar{\mathcal{I}}_{k}}$. To balance the exploration-exploitation trade-off, a common choice of $\widehat{b}_{\bar{\mathcal{I}}_{k}}$ is the $\epsilon$-greedy policy with respect to $\widehat{\pi}_{\bar{\mathcal{I}}_{k}}$.

Simulations

We consider three scenarios. In Scenarios (A) and (B), the system dynamics are given by

for $t\geq 0$, where $\{z_{t}\}_{t\geq 0}\stackrel{{\scriptstyle iid}}{{\sim}}N(0_{2},I_{2}/4)$ and $X_{0,0}\sim N(0_{2},I_{2})$. In Scenario (A), we consider a completely randomized study and set $\{A_{0,t}\}_{t\geq 0}$ to i.i.d Bernoulli random variables with expectation $0.5$. In Scenario (B), we allow the treatment assignment mechanism to depend on the observed state. Specifically, we set $m=2$ and ${\mbox{Pr}}(A_{0,t}=1|X_{0,t})=0.5\textrm{sigmoid}(X_{0,t,1})+0.5\textrm{sigmoid}(X_{0,t,2})$ where $X_{0,t,i}$ denotes the $i$th element in $X_{0,t}$. The target policy we consider is designed as follows,

In Scenario (C), we consider a standard RL setting included in OpenAI Gym (Brockman et al. 2016): Cliff Walking. This RL example is detailed in Example 6.6 in Sutton and Barto 2018. The objective is to identify the optimal path from the starting point S to the destination point G without falling off the cliff (see Figure 1). This scenario corresponds to an episodic task where the agent will be sent instantly to the starting point wherever it steps into the cliff or arrives at the destination. We manually add some noises to the immediate rewards simulated by the OpenAI Gym to ensure that the system dynamics are not deterministic. We remark that this task is considered in Kallus and Uehara 2020 as well. The target policy is the optimal policy and the behavior policy is a 50-50 mixture of the optimal and uniform random policies.

The DRL estimator has been shown to be much more efficient than AIPW or IPW estimators (Thomas et al. 2015; Jiang and Li 2016). So we focus on comparing our approach with DRL. DRL requires the calculation of the Q-function, the marginalized density ratio and the behavior policy. Here, we treat the behavior policy as known and estimate the Q-function and the density ratio based on nonparametric sieve regression.

In Figure 2 and Table 1, we report the empirical coverage probabilities (ECPs) and average lengths (ALs) of CIs constructed by the proposed method, with different choices of $n$ and $T$. It can be seen that our CI achieves nominal coverage in all cases. Its length decreases as $nT$ increases. This is consistent with our theoretical findings where we show the proposed value estimator converges at a rate of $n^{-1/2}T^{-1/2}$ under certain conditions (see the discussions below Theorem 1).

Comparing our method with DRL, it is clear that our CIs are in general narrower than those constructed by DRL. In addition, MSEs of the proposed value estimates are smaller than those based on DRL. This is consistent with our theoretical analysis in Appendix C.2 where we show the variance of our value estimator is strictly smaller than that based on DRL under certain conditions. In addition, it can be seen from Table 1 that ECPs of DRL are below 90% in Scenario (C).

2 Off-policy evaluation with an (estimated) optimal policy

In addition, we design a non-regular setting Scenario (D) where the actions do not have effects on the transition dynamics or the immediate rewards. Specifically, for any $t\geq 0$, we set

3 On-policy evaluation with an (estimated) optimal policy

Application to the OhioT1DM dataset

As commented in the introduction, this dataset contains eight weeks’ records of CGM blood glucose levels, insulin doses and self-reported life-event data for each of six patients with type 1 diabetes. To analyze this data, we divide these eight weeks into three hour intervals. The state variable $X_{i,t}$ is set to be a three-dimensional vector. Specifically, its first element $X_{i,t}^{(1)}$ is the average CGM blood glucose levels during the three hour interval $[t-1,t)$. The second covariate $X_{i,t}^{(2)}$ is constructed based on the $i$-patient’s self-reported time and the carbohydrate estimate for the meal. Suppose the patient has meals at time $t_{1},t_{2},\dots,t_{N}\in[t-1,t)$ with the carbohydrate estimates $\hbox{CE}_{1},\hbox{CE}_{2},\dots,\hbox{CE}_{N}$. Define

where $\gamma_{c}$ corresponds to the decay rate every five minutes. Here, we set $\gamma_{c}=0.5$. The third covariate $X_{i,t}^{(3)}$ is defined as an average of the basal rate during the three hour interval.

We discretize the action according to the amount of insulin injected in the three hour interval. Specifically, $A_{i,t}=1$ when the total amount of insulin delivered to the $i$-th patient is greater than one unit. Otherwise, we set $A_{i,t}=0$. The immediate reward $Y_{i,t}$ is defined according to the Index of Glycemic Control (Rodbard 2009, IGC,), which is a non-linear function of the blood glucose levels. Specifically, we set

A large IGC indicates the patient is in good health status. We set the discount factor $\gamma=0.5$, as in simulations.

For the $i$-th patient, we apply the double FQI algorithm to the data

to estimate a patient-specific optimal policy. Then we compute the estimator for the value function $V(\pi^{\tiny{opt}};X_{i,0})$ starting from the initial state variable $X_{i,0}$. In addition, we extend our methodology in Section 5.2 to construct the confidence interval for the value difference $V(\pi^{\tiny{opt}};X_{i,0})-V(b;X_{i,0})$ where $V(b;X_{i,0})$ corresponds to the value under the behavior policy. See Appendix B.2 for details. In Figure 4, we plot our proposed CI for the value difference, for each of the six patients, when the initial starting time is either 8:00 am or 2:00 pm in Day 1. It can be seen that the estimated value differences are strictly positive in all cases. This implies that the optimal value is strictly larger than the observed discounted cumulative reward. In some cases, the lower bound of our CI is larger than zero. The difference is thus significant. In Figure 5, we fix the starting time to 8:00 am and plot the CI of the value difference with different $\gamma$. Results show a similar qualitative pattern. This suggests applying reinforcement learning algorithms could potentially improve some patients’ health status.

Discussion

We discuss the advantages and limitations among the proposed method and DRL when inferring the value under a fixed decision rule. Generally speaking, the proposed method results in narrower CIs and would be preferred in cases where $(nT)^{-1/2}$-consistent estimation of the Q-function is feasible. This includes settings where the dimension of the state-vector is not large, as in our real data applications. In contrast, DRL would be preferred in ergodic environments with high-dimensional covariates where $(nT)^{-1/2}$-consistent estimation of the Q-function is infeasible.

Specifically, in Appendix C.2, we consider settings where both the behavior policy and the target policy are nondynamic. Under certain conditions, we prove that the variance of the DRL estimator is strictly larger than that of the proposed estimator. This in turn implies that our method yields a narrower CI in general.

In addition, we remark that the CI constructed by DRL requires the data to be generated from an ergodic environment. In the Cliff Walking example, the data are generated by a mixture of the optimal and random policy. Since the agent will be sent instantly to the starting point wherever it steps into the cliff or arrives at the destination, the Markov chain formed by the state-action pair is no longer ergodic. It can be seen from Table 1 where ECP of the CI is well below the nominal level in the Cliff Walking example. Although our procedure also requires the ergodicity assumption (see Condition (A3)(ii)), this assumption is not necessary. It can be seen from the proof of Theorem 1 that our CI is valid as long as the random matrix $\widehat{\bm{\Sigma}}_{\pi}$ stabilizes. This is consistent with the findings in Table 1 where our CI achieves nominal coverage in all cases.

However, to ensure the proposed CI is valid, we require the bias of our Q-estimator to decay at a rate of $o((nT)^{-1/2})$. Consequently, our estimator converges at a rate of $(nT)^{-1/2}$. This rate might not be achievable in high-dimensions. In contrast, DRL requires a weaker condition. The CI based on DRL is valid when both the Q-estimator and the estimated marginalized density ratio converge at a rate of $o((nT)^{-1/4})$.

Another potential limitation of our method is that in cases where $\widehat{\bm{\Sigma}}_{\pi}$ is close to a singular matrix, the resulting Q-estimator might suffer from over-fitting, leading to an unbounded outcome. In practice, we could add a ridge penalty to reduce over-fitting. We discuss in detail in Appendix D.4.

2 Number of basis functions

We outline a procedure to choose the number of basis function $L$ in this section. The idea is to simulate the model dynamics and select $L$ such that the resulting confidence interval achieves nominal coverage under the simulated model. Specifically, given the observed data $\{(X_{i,t},A_{i,t},Y_{i,t},X_{i,t+1})\}_{0\leq t<T_{i},1\leq i\leq n}$, we propose to learn the conditional density function of $(Y_{i,t},X_{i,t+1})$ given $(A_{i,t},X_{i,t})$. Following Janner et al. 2019, we recommend to use a Gaussian distribution to model the conditional density function in practice,

The conditional mean $\mu$ can be estimated via nonparametric regression (e.g., random forest). Let $\widehat{\mu}$ denote the corresponding estimator. Given the set of estimated residuals $\{(Y_{i,t},X_{i,t+1})^{\top}-\widehat{\mu}(A_{i,t},X_{i,t})\}_{i,t}$ the conditional covariance function $\Sigma$ can be estimated via regression as well. We remark that in addition to the Gaussian function, other density functions could be used to model the system dynamics as well.

The behavior policy $b$ can be similarly estimated via regression. Given an estimated behavior policy $\widehat{b}$ and $\widehat{\mu}$, $\widehat{\Sigma}$, we generate simulated trajectories to investigate the performance of the proposed CI with different choices of $L$.

Finally, we choose $L$ such that the resulting CI is the shortest among all CIs whose coverage probabilities are above certain level (e.g., 93%) under the simulated environment.

In Appendix D.1, we investigate the finite sample performance of such a method and find that it performs reasonably well. We remark that alternative to the aforementioned method, cross-validation could be applied to select $L$.

3 Sensitivity to the ordering of trajectories

The proposed sequential value evaluation procedure in Section 3.2 divides the data into blocks defined both by trajectories and by time. While there is a natural order in time, there does not appear to be a natural order in the trajectories. In Appendix D.2.3, we conduct additional simulation studies to investigate the sensitivity of our CI to the ordering of trajectories under Scenario (B). Results suggest that our CI is not overly sensitive under our simulation setting.

As suggested by one of the referees, we may aggregate CIs over multiple orderings in cases where the results depend strongly on the ordering of the trajectories. Dezeure et al. 2015 derived a CI for the regression coefficients in high-dimensional models by aggregating results over multiple sample splits using a quantile function. We can adopt their method to aggregate our CIs over multiple orderings. Alternatively, one may average the value estimates over sufficiently many orderings and apply similar methods developed in Wang et al. 2020; Shi et al. 2020a to derive the CI. However, these algorithms are much more time-consuming.

4 More on value-based method

In Section 3.1.1, we discuss a potential drawback of using nonparametric methods to directly model the value function. We remark that a kernel-type importance sampling estimator for the $V(\pi;x)$ will not suffer from this issue, since it does not directly model the value function, but uses an inverse propensity-score weighted method instead. Both IPW and regression type estimators have their own merits. In general, IPWEs might suffer from a large variance whereas regression-based estimators might suffer from a large bias. There exist methods that combine both for more robust off-policy evaluation (Kallus and Uehara 2019; Uehara et al. 2020; Tang et al. 2020; Shi et al. 2021, see e.g.,). However, as commented in Section 7.1, they might yield larger CIs compared to our method.

5 Rate of convergence of Q-learning type algorithms

Through authors’ communication, we found a recent independent work by Hu et al. 2021 that derived a nonasymptotic error bound on the value of the estimated optimal policy computed by Q-learning type algorithms under the margin condition. Their results are consistent with our theoretical findings in Theorems 3 and 4 that show the value of the estimated optimal policy converges to the optimal value at a faster rate than the estimated Q-function.

References

Appendix A Some technical conditions

where $\|\cdot\|_{TV}$ denotes the total variation norm.

A.2 Conditions A3*

(ii) The Markov chain $\{X_{0,t}\}_{t\geq 0}$ is geometrically ergodic.

We remark that Condition A3*(ii) is the same as A3(ii).

A.3 More on the margin condition

To better understand Condition A5, we consider a simple scenario where $\mathcal{A}=\{0,1\}$. Define $\tau(x)=Q^{\tiny{opt}}(x,1)-Q^{\tiny{opt}}(x,0)$. It follows that

As a result, (3.22) and (3.23) are equivalent to the followings:

for some $\alpha>0$. Then, with some calculations, we can show

Appendix B Additional details regarding the method

and $|\mathcal{I}|$ stands for the number of elements in $\mathcal{I}$.

B.2 Value difference between the target and behavior policy

In this section, we outline a method to evaluate the value difference function between the target and behavior policy. We first consider the scenario where the target policy is a fixed policy. We next consider the scenario where the target policy is an estimated optimal policy. To simplify the presentation, we assume $T_{1}=\cdots=T_{n}=T$. The proposed method can be similarly extended to on-policy settings.

Consider a data-independent policy $\pi$. We aim to evaluate the value difference function $\textrm{VD}(\pi;x)=V(\pi;x)-V(b;x)$ where $b$ is the unknown behavior policy. We first apply our method in Section 3.1.2 to compute an estimator value function $\widehat{V}(\pi;x)$ for $V(\pi;x)$.

The resulting estimates for $Q(b;x,a)$ can be derived as $\Phi_{L}^{\top}(x)\widehat{\beta}_{b,a}$. The corresponding estimator for $V(b,x)$ is given by $\widehat{V}(b,x)=\sum_{a}\widehat{b}(a|x)\Phi_{L}^{\top}(x)\widehat{\beta}_{b,a}$ where $\widehat{b}(a|x)$ denotes the sieve estimator $\Phi_{L}^{\top}(x)\widehat{\alpha}_{a}$ for $b(a|x)$ where

This yields the estimator for the value difference $\widehat{\textrm{VD}}(\pi;x)=\widehat{V}(\pi;x)-\widehat{V}(b;x)$.

We next derive a confidence interval for VD$(\pi;x)$ based on $\widehat{\textrm{VD}}(\pi;x)$. Similar to the proof of Theorem 1, we can show $\sqrt{nT}\{\widehat{\textrm{VD}}(\pi;x)-\textrm{VD}(\pi;x)\}$ is equivalent to

where $\varepsilon_{i,t}^{*}$ denotes the temporal difference error $Y_{i,t}+\gamma Q(b;X_{i,t+1},A_{i,t+1})-Q(b;X_{i,t},A_{i,t})$ and $\bm{\Psi}$ denotes the population limit of $\widehat{\Psi}$. Note that the RHS can be rewritten as $(nT)^{-1/2}\sum_{i=1}^{n}\sum_{t=0}^{T-1}\psi_{i,t}$ that corresponds to a sum of martingale difference. Its variance can be consistently estimated by $\widehat{\sigma}^{*2}(\pi;x)=(nT)^{-1}\sum_{i=1}^{n}\sum_{t=0}^{T-1}\widehat{\psi}_{i,t}^{2}$ where $\widehat{\psi}_{i,t}$ denotes some consistent estimator for $\psi_{i,t}$ based on $\widehat{Q}(\pi;\cdot,\cdot)$, $\widehat{Q}(b;\cdot,\cdot)$ and $\widehat{b}$. The confidence interval for VD$(\pi;x)$ is given by

B.2.2 Inference of the value difference under an estimated optimal policy

We begin by dividing the data into $K$ non-overlapping subsets $\cup_{k=1}^{K}\mathcal{I}_{k}$. Similar to Section 3.2.2, we construct the value difference estimator by

where $\widehat{\textrm{VD}}_{\mathcal{I}_{k}}$ and $\widehat{\sigma}^{*}_{\mathcal{I}_{k}}$ denote the versions of VD and $\widehat{\sigma}^{*}$ based on samples in $\mathcal{I}_{k}$ only. The corresponding confidence interval is given by

where $\widetilde{\sigma}^{*}(x)=(K-1)\{\sum_{k=1}^{K-1}\widehat{\sigma}_{k}^{*-1}(\widehat{\pi}_{\bar{\mathcal{I}}_{k}};x)\}^{-1}$.

Finally, we remark that such a confidence interval might not be valid in the extreme case where the behavior policy is equal to a deterministic optimal policy. To elaborate, notice that when the behavior policy is deterministic, the second line of (B.32) equal zero. In addition, when $\pi=b=\pi^{\tiny{opt}}$ for some optimal policy $\pi^{\tiny{opt}}$, the first line equals zero as well. In that case, $\widehat{\textrm{VD}}(\pi;x)$ would have a degenerate distribution. Suppose the estimated optimal policy is consistent for $\pi^{\tiny{opt}}$. Then $\widetilde{\textrm{VD}}(x)$ might not have a tractable limiting distribution, leading to an invalid confidence interval.

To address this concern, we could redefine the inverse weights $\widehat{\sigma}^{*}_{\mathcal{I}_{k+1}}(\widehat{\pi}_{\bar{\mathcal{I}}_{k}};x)$ by $\widehat{\sigma}^{*}_{\mathcal{I}_{k+1}}(\widehat{\pi}_{\bar{\mathcal{I}}_{k}};x,\delta)=\max\{\widehat{\sigma}^{*}_{\mathcal{I}_{k+1}}(\widehat{\pi}_{\bar{\mathcal{I}}_{k}};x),\delta\}$ for some $\delta>0$, as in Luedtke and van der Laan 2017. This guarantees that these inverse weights are strictly greater than zero. A similar approach is employed by Shi et al. 2020b for testing the overall qualitative treatment effects in single-stage decision making. In addition, one could allow $\delta$ to depend on $N$ and $T$. The resulting confidence interval would be valid as long as $\delta\gg(NT)^{-1/6}$(Shi et al. 2020b, see e.g., Theorem 3.1 of). However, a potential limitation is that it would yield a conservative confidence interval when the truncation is active, as discussed in Luedtke and van der Laan 2017.

B.3 Double fitted QQ-iteration

In this section, we introduce our algorithm for computing the estimated optimal policy in our numerical studies. The proposed algorithm is based on FQI that recursively updates the estimated optimal Q-function by some supervised learning method (see Example 2 in Section 3.2.3). In FQI, at each iteration, a maximization over estimated Q-function is used as an estimate of the maximum of the true Q-function. This can lead to a significant positive bias (Sutton and Barto 2018). Hasselt 2010 proposed a double Q-learning method to reduce the maximization bias. Here, we apply similar ideas to FQI to compute the estimated optimal policy. We use a pseudocode to summarize our algorithm below.

In Algorithm 1, we can apply any non-parametric models $Q(\cdot,\texttt{;}\theta)$ indexed by $\theta$ to model the optimal Q-function. In our implementation, we set $Q(\cdot,\cdot;\cdot)$ to be a linear combination of tensor product B-spline basis functions.

Appendix C Additional technical details

is positive semidefinite. It follows that

When $\pi$ is a deterministic policy, $\sum_{a\in\mathcal{A}}\bm{\xi}(x,a)\bm{\xi}^{\top}(x,a)b(a|x)-\gamma^{2}\bm{U}_{\pi}(x)\bm{U}_{\pi}^{\top}(x)$ is a block diagonal matrix. To show A4(i) holds, it suffices to show

Suppose $b$ is the $\epsilon$-greedy policy with respect to $\pi$, i.e, $b(a|x)=\epsilon m^{-1}+(1-\epsilon)\pi(a|x)$, for any $a\in\{1,\dots,m\}$ and $\epsilon$ satisfies $\epsilon\leq 1-\gamma^{2}$, we have

The condition in (C.33) is automatically satisfied when A2 holds (Burman and Chen 1989; Chen and Christensen 2015, see, e.g.,).

C.2 Additional details on the variance comparison

We consider a randomized study where $b(a|x)$ is a constant function of $x$. In addition, we assume the target policy is nondynamic, i.e., $\pi(a^{*}|x)=1$ for some $1\leq a^{*}\leq m$ and any $x$. We impose the following conditions.

(C1) The process $\{X_{0,t}\}_{t\geq 0}$ is stationary.

(C2) The temporal difference error $\varepsilon_{0,t}$ is independent of $(X_{0,t},A_{0,t})$.

We make some remarks. First, Condition (C1) is imposed to simplify the presentation. The same results hold as long as $\{X_{0,t}\}_{t}$ will converge to its stationary distribution. Second, the variances of our estimator and DRL are very difficult to analyse in general. Conditions (C2)-(C3) are imposed to simplify the calculation. Even when these conditions are violated, we expect the variance of the proposed estimator will be smaller in general, as reflected in our numerical study.

We next sketch a few lines to prove Theorem 5. Based on Theorem 1, the asymptotic variance of our estimator is given by

Using similar arguments in the proof of Theorem 16 of Kallus and Uehara 2019, we can show that the asymptotic variance of the DRL estimator equals

Under the given conditions, the above variance is equal to $\{nT(1-\gamma)^{2}p_{a^{*}}\}^{-1}\sigma_{*}^{2}\{1+{\mbox{Var}}(\omega(X_{0,0}))\}$ where $p_{a^{*}}={\mbox{Pr}}(A_{0,t}=a^{*})$. Consequently, it suffices to show

By the definition of $\bm{\Sigma}_{\pi}$ and $\bm{\xi}_{0,t}$, the left-hand-side of (C.35) is equal to

As such, the left-hand-side of (C.35) is upper bounded by

C.3 Additional details on on-policy evaluation

In this section, we show our proposed CI in Section 4 achieves nominal converge. To simplify the analysis, we focus on the setting where $K$ is finite, $T(1)=\cdots=T(K)=T$ and $L(1)=\cdots=L(K)=L$. When $K$ diverges, the sequences $\{T(k)\}_{k\geq 1}$ and $\{L(k)\}_{k\geq 1}$ shall be properly chosen to reduce the bias of the value estimates. We leave this for future research.

Similar to Appendix A.2, we assume the estimated policy $\widehat{\pi}_{\mathcal{I}}\in\Pi$ with probability $1$, for any $\mathcal{I}$. In on-policy settings, the behavior policy $b_{\pi}$ is a function of the estimated policy $\pi\in\Pi$. For instance, when an $\epsilon$-greedy policy is used to determine the behavior policy, then we have $b_{\pi}=(1-\epsilon)\pi+\epsilon\pi^{*}$ where $\pi^{*}$ denotes a uniform random policy. Let $\mathcal{B}=\{b_{\pi}:\pi\in\Pi\}$.

(A2’.) Assume $\nu_{0}$ and $q$ are uniformly bounded away from $0$ and $\infty$ on their supports.

(iii) There exists some constant $\bar{c}>0$ such that

Proof of Theorem 6 is omitted for brevity.

Appendix D Additional numerical results

We apply the proposed method detailed in Section 7.2 to Scenario (B) where the treatment assignment mechanism depends on the observed state, to investigate the finite sample performance of the resulting CI. Specifically, we apply the random forest algorithm to learn the conditional mean function $\mu$ and the behavior policy $b$. We assume $\Sigma(a,x)$ is a constant function of $(a,x)$ and estimated it by

We use the tensor product B-spline basis for $\Phi_{L}$, as in Section 5. Note that the state is a two-dimensional vector, $L$ is selected among the set $\{4^{2},5^{2},6^{2},7^{2},8^{2},9^{2}\}$. Specifically, we choose $L$ such that the resulting CI is the shortest among all CIs whose coverage probabilities are above 93%. If no such CI exists, we select the CI with the highest coverage probability.

We report the ECP and AL of the resulting CI in the left and middle panels of Figure 7. It can been seen that ECP is close to the nominal level in all cases and AL decays as either $n$ or $T$ increases. In the right panel of Figure 7, we report the number of basis functions that is being selected most by the proposed method as a function of $n$ and $T$ (denote by $L(n,T)$). It is clear from Figure 7 that $L(n,T)$ increases with the total number of observations $nT$. This is consistent with the following intuition: as $nT$ increases, more basis functions are needed to reduce the approximation error and guarantee the nominal coverage of the resulting CI.

D.2 Sensitivity analysis

In this section, we conduct the sensitivity test for the parameter $\eta$ in the number of basis $L=\lfloor(nT)^{\eta}\rfloor$. We consider the simulation of the off-policy evaluation with a fixed target policy in Section 5.1. For scenario (A), (B) and (C), we set $n=100$, $T=100$ and the different $\eta$’s are chosen from $(0.25,0.30,0.35,0.40,3/7,0.45)$. The result of the ECPs are plotted in Figure 8 where all the ECPs are close to the nominal coverage rate 0.95. It shows that the results of the coverage are not sensitive to the different choices of $\eta$.

D.2.2 Sensitivity test for γ\gamma

In Figure 9, we report the ECP and AL of the proposed CI and the MSE of our value estimate under Scenario B where the target policy is fixed, with $\gamma=0.3$ and $0.7$. In Figure 10, we report the ECP and AL of the proposed CI and the MSE of our value estimate under Scenario B where the target policy is an estimated optimal policy, with $\gamma=0.3$ and $0.7$. It can be seen that findings are very similar to those with $\gamma=0.5$.

In Tables 4 and 5, we report the ECP, AL and MSE of the proposed method and DRL under Scenario (C) where the target policy is fixed, with $\gamma=0.3$ and $0.7$. It can be seen that the proposed CI achieves nominal coverage in all cases. When $\gamma=0.3$, ECP of the DRL method is well below the nominal level in all cases. When $\gamma=0.7$, the AL and MSE of the proposed CI are much smaller than those based on DRL.

D.2.3 Sensitivity test for the ordering of trajectories

We focus on Scenario (B), detailed in Section 5.1, to examine the sensitivity of the proposed CI to the ordering of trajectories. Specifically, we first randomly permute all trajectories with some fixed random seed. We next apply our SAVE procedure to construct the CI. We use three random seeds to generate different random permutations and depict the corresponding results in Figure 11. It can be seen that our method is not overly sensitive to the ordering of trajectories.

D.3 Additional settings

D.4 Additional real data results

We use our real data example to discuss the issue of over-fitting in this section. Specifically, we apply the proposed method in Section 3.2.2 to evaluate the optimal value $V(\pi^{\tiny{opt}};X_{i,0})$ starting from the initial state variable $X_{i,0}$, for $i=1,2,\cdots,6$. When the initial starting time is 8:00 am in Day 1, CIs for Patient 5 and Patient 6 are $[-6.288,3.287]$ and $[-6.313,10.644]$, respectively. Both upper bounds are positive. However, according to our definition, the immediate reward is nonpositive. As such, the value and Q-function shall be nonpositive as well. This reflects one of the drawback of the proposed method. The resulting Q-estimator might suffer from over-fitting, leading to an unbounded outcome.

Specifically, it is due to that the matrix $\widehat{\bm{\Sigma}}_{\pi}$ is close to singular. Note that the regression coefficients $\widehat{\bm{\beta}}_{\pi}$ are computed by solving the linear equation

In our data example, the number of basis function equals 12. As such, $\widehat{\bm{\Sigma}}_{\pi}$ is a $12$ by $12$ matrix. When it is close to singular, the resulting Q-estimator might be unbounded.

To avoid offer-fitting, we note that in theory, $\widehat{\bm{\Sigma}}_{\pi}$ is a positive definite matrix under Condition (A3)(i). This motivates us to compute $\widehat{\bm{\beta}}_{\pi}$ by solving

where $I$ denotes the identity matrix. As long as $\lambda$ satisfies $\lambda=O(N^{-1}T^{-1})$, the proposed CI remains valid. In our real data example, we set $\lambda=1\times 10^{-9}$. The resulting CIs for Patient 5 and Patient 6 are $[-8.034,-5.667]$ and $[-9.866,-7.104]$. Both upper bounds are strictly negative.

Appendix E Technical proofs

The rest of the section is organized as follows. We first present the proof sketches for Theorems 1-4. We next present the detailed technical proofs.

We provide an outline for the proof in this section. The detailed proof can be found in Section E.5 of the supplementary article. We break the proof into three steps. In the first step, we show the estimator $\widehat{\bm{\beta}}_{\pi}$ satisfies

In the second step, we show the linear representation in (3.14) holds. The proof of (3.14) relies on the convergence rate of $\widehat{\bm{\beta}}$ established in the first step and some additional random matrix nequalities in Lemma 4.

In the last step, we show the leading term on the RHS of (3.14) is asymptotically normal, based on the martingale central limit theorem. The completes the proof of Theorem 1.

E.2 A sketch for the proof of Theorem 2

Similar to (3.19), for each $1\leq k\leq K$, we have

where $R_{k}^{(1)}$ denotes the remainder term and

for some remainder term $R_{k}^{(2)}$. Suppose $R_{k}^{(1)}$ and $R_{k}^{(2)}$ satisfy certain convergence rates. Combining (E.38) together with (E.37) yields,

due to our use of the inverse weighting trick. Theorem 2 thus follows from the martingale central limit theorem.

E.3 A sketch for the proofs of Theorems 3 and 4

For $j=0,1,2,\dots$, define a time-dependent policy $\widehat{\pi}_{\mathcal{I}}^{(j)}$ that executes $\widehat{\pi}_{\mathcal{I}}$ at the first $j$ time points and then follows $\pi^{\tiny{opt}}$. By definition, we have $\pi^{\tiny{opt}}=\widehat{\pi}_{\mathcal{I}}^{(0)}$ and $\widehat{\pi}_{\mathcal{I}}=\widehat{\pi}_{\mathcal{I}}^{(\infty)}$. Notice that

Let $q_{X}^{(j)}(\cdot|x)$ be the density function of $X_{0,j}$ conditional on $X_{0,0}=x$, following the estimated policy $\widehat{\pi}_{\mathcal{I}}$ at the first $j$ time points, we have

By A1, we have $\sup_{x,x^{\prime},a}q(x^{\prime}|x,a)\leq c$. Under the Markov assumption,

In addition, $\sum_{a\in\mathcal{A}}Q^{\tiny{opt}}(x^{\prime},a)\{\pi^{\tiny{opt}}(a|x^{\prime})-\widehat{\pi}_{\mathcal{I}}(a|x^{\prime})\}\geq 0$ for any $x^{\prime}$, by the definition of $\pi^{\tiny{opt}}$. Therefore, we obtain

E.4 Proof of Lemma 1

Since $r(\cdot,a)$ is $p$-smooth for any $a\in\mathcal{A}$, it suffices to show

is $p$-smooth for any $a\in\mathcal{A}$ and any policy $\pi$.

as $j\to\infty$. By the mean value theorem, we have

where $0\leq\theta_{x}\leq 1$ for all $x$. When $1<p\leq 2$, we have $\lfloor p\rfloor=1$. It follows from Condition A1 that

When $p>2$, $|\partial_{j}\partial_{j}q(x^{\prime}|x,a)|$ exists and is bounded by $c$ for any $x^{\prime}$, $x$ and $a$. It follows from the mean value theorem that

In addition, it follows from A1 and (E.44) that

where $\lambda(\cdot)$ denotes the Lebesgue measure. Using the same arguments, we can show for any $d$-tuple $\alpha=(\alpha_{1},\dots,\alpha_{d})^{\top}$ of nonnegative integers that satisfies $\|\alpha\|_{1}\leq\lfloor p\rfloor$,

Moreover, by A1, (E.44) and (E.45), we have for any $d$-tuple $\alpha$ with $\|\alpha\|_{1}=\lfloor p\rfloor$ that

E.5 Proof of Theorem 1

There exists some constant $c^{*}\geq 1$ such that

Suppose the conditions in Theorem 1 hold. We have as either $n\to\infty$ or $T\to\infty$ that $\|\bm{\Sigma}^{-1}\|_{2}\leq 3\bar{c}^{-1}$, $\|\bm{\Sigma}\|_{2}=O(1)$, $\|\widehat{\bm{\Sigma}}-\bm{\Sigma}\|_{2}=O_{p}\{L^{1/2}(nT)^{-1/2}\log(nT)\}$, $\|\widehat{\bm{\Sigma}}^{-1}-\bm{\Sigma}^{-1}\|_{2}=O_{p}\{L^{1/2}(nT)^{-1/2}\log(nT)\}$ and $\|\widehat{\bm{\Sigma}}^{-1}\|\leq 6\bar{c}^{-1}$ wpa1.

for some constant $C>0$. Let $\bm{\beta}^{*}=(\beta_{1}^{*T},\dots,\beta_{m}^{*T})^{\top}$, and

The condition ${\mbox{Pr}}(\max_{0\leq t\leq T-1}|Y_{i,t}|\leq c_{0})=1$ implies that $|Y_{i,t}|\leq c_{0},\forall i,t$, almost surely. By Lemma 1 and the definition of $p$-smooth functions, we obtain that $|Q(\pi;x,a)|\leq c^{\prime}$ for any $\pi,x,a$. It follows that

almost surely. In addition, it follows from (E.48) that

In the following, we show $\zeta_{2}=O_{p}\{L(nT)^{-1}\log(nT)\}$ and $\zeta_{3}=O_{p}(L^{-p/d})$ as either $n\to\infty$, or $T\to\infty$.

Error bound for $\|\zeta_{2}\|_{2}$: Let $\mathcal{F}_{i,t}$ denote the sub-dataset $\{X_{i,t},A_{i,t}\}\cup\{(Y_{i,j},A_{i,j},X_{i,j})\}_{1\leq j<t}$. By the Bellman equation in (3.9), MA and CMIA, we have

Notice that $\bm{\xi}_{i,t}$ is a function of $X_{i,t}$ and $A_{i,t}$ only, we have for any $0\leq t_{1}<t_{2}\leq T-1$ that

By Markov’s inequality, we obtain $(nT)^{-1}\sum_{i=1}^{n}\sum_{t=0}^{T-1}\bm{\xi}_{i,t}\varepsilon_{i,t}=O_{p}\{\sqrt{L/(nT)}\}$. Combining this together with Lemma 3 yields that $\zeta_{2}=O_{p}\{\sqrt{L/nT}\log(nT)\}O_{p}\{\sqrt{L/(nT)}\}=O_{p}\{L(nT)^{-1}\log(nT)\}$.

By Lemma 3, we have $\|\widehat{\bm{\Sigma}}^{-1}\|_{2}=O_{p}(1)$. Combining this together with (E.51) yields that $\zeta_{3}=O_{p}(L^{-p/d})$.

This completes the first step of the proof.

Step 2: Using similar arguments in bounding $\|\zeta_{2}\|_{2}$ in Step 1, we can show that $\|\zeta_{1}\|_{2}=O_{p}\{L^{1/2}(nT)^{-1/2}\}$ and thus

Since $\inf_{x,a}\omega(x,a)\geq c_{0}^{-1}$, it follows from Lemma 4 that

By Lemma 3, we have $\|\bm{\Sigma}\|_{2}=O(1)$, or equivalently, $\lambda_{\max}(\bm{\Sigma}^{\top}\bm{\Sigma})=O(1)$. This implies that $\lambda_{\min}\{\bm{\Sigma}^{-1}(\bm{\Sigma}^{\top})^{-1}\}\geq\bar{C}$ for some constant $\bar{C}>0$ and hence

by (E.56). Combining (E.57) together with (E.55) yields that

This completes the second step of the proof.

Let $\mathcal{F}^{(0)}=\{X_{1,0},A_{1,0}\}$. Then we iteratively define $\{\mathcal{F}^{(g)}\}_{1\leq g\leq nT}$ as follows:

Let $\bm{\xi}^{(g)}=\bm{\xi}_{i(g),t(g)}$ and $\varepsilon^{(g)}=\varepsilon_{i(g),t(g)}$. It follows that

By MA, CMIA and the Bellman equation in (3.9), we obtain that

Hence, the RHS of (E.59) forms a martingale with respect to the filtration $\{\sigma(\mathcal{F}^{(g)})\}_{g\geq 0}$, where $\sigma(\mathcal{F}^{(g)})$ stands for the $\sigma$-algebra generated by $\mathcal{F}^{(g)}$. To show the asymptotic normality, we use a martingale central limit theorem for triangular arrays (McLeish 1974, Corollary 2.8 of). This requires to verify the following two conditions:

where the first inequality follows from Cauchy-Schwarz inequality, the second inequality is due to (E.49), the third inequality is due to Lemma 2 and the fact that $\|\bm{\xi}^{(g)}\|_{2}\leq\sup_{x}\|\Phi_{L}(x)\|_{2}$, and the last inequality follows from (E.56). Since $L\ll\sqrt{nT}/\log(nT)$, (a) is proven. To verify (b), notice that

This can be proven using similar arguments in bounding $\|\widehat{\bm{\Sigma}}-\bm{\Sigma}\|_{2}$ in the proof of Lemma 3. In view of (E.58) and (E.59), we have by Slutsky’s theorem that

and hence $\|\bm{\Omega}\|_{2}=O(1)$. This together with Lemma 3 and the condition $L\ll\sqrt{nT}/\log(nT)$ yields that

Thus, it remains to show $\|\widehat{\bm{\Sigma}}^{-1}\widehat{\bm{\Omega}}(\widehat{\bm{\Sigma}}^{\top})^{-1}-\widehat{\bm{\Sigma}}^{-1}\bm{\Omega}(\widehat{\bm{\Sigma}}^{\top})^{-1}\|_{2}=o_{p}(1)$, or $\|\widehat{\bm{\Omega}}-\bm{\Omega}\|_{2}=o_{p}(1)$, by Lemma 3. In view of (E.60), it suffices to show $\|(nT)^{-1}\sum_{g=1}^{nT}(\varepsilon^{(g)})^{2}\bm{\xi}^{(g)}(\bm{\xi}^{(g)})^{\top}-\widehat{\bm{\Omega}}\|_{2}=o_{p}(1)$, or equivalently,