Simultaneous Statistical Inference for Off-Policy Evaluation in Reinforcement Learning

We sincerely appreciate the time and effort you dedicated to reviewing our manuscript and providing insightful feedback. Below, we address each of your comments in detail.

Q1&W2. Regarding the offline RL, challenges, and the goodness of dataset.

We will include [1]-[5] as literature on offline RL and highlight the challenges in offline RL. The following argument will also apear in the revision. Unlike online environments (e.g. Go, MiniGrid), collecting data with online interaction in many real-world applications (e.g. healthcare, autonomous driving) can be expensive or even dangerous. This limitation prevents the widespread adoption of traditional online RL methods. As an alternative, offline RL has been proposed to effectively utilize pre-collected experience data without requiring any online interaction.

However, this offline approach faces several challenges that make it different from online RL. For learning, the key challenges are:

• Distributional shift: The data generated by the behavior policy (used to collect the dataset) might mismatch the target policy (which is learned), which leads to poor performance or overfitting.
• Sample inefficiency: The learning must rely entirely on the static dataset, which means 'exploration' (as in online RL’s exploration-exploitation trade-off) is not possible during the learning process.

For evaluation, offline RL faces the difficulty of estimating policy performance without online deployment, necessitating off-policy evaluation (OPE)—the focus of our work. One of the core challenges in OPE is the gap between estimated off-policy evaluation and real-world performance. To address this, recent attempts, including pointwise CIs (as we concluded in the related work section), have been proposed to quantify the gap between offline assessments and real-world deployment through probabilistic measures. Our work strengthens the statistical tools from pointwise to global scale, providing more reliable guarantees for decision-making.

For the goodness of the dataset, the quality of an offline dataset depends critically on several factors:

• State-space coverage: Breadth of visited states
• Reward signal diversity: Range of observed returns
• Action variety: Exploration of possible actions
• Algorithm compatibility: Match between data distribution and RL method

As noted in [2], formalizing a non-trivial "sufficiency" condition for dataset quality remains an open problem. In our framework, we distill these complex considerations into Assumptions (A1)-(A3) (detailed in response Q5&W1). Under these assumptions, the goodness of the dataset can be roughly quantified by the sample size $N$. This simplification enables tractable theoretical analysis.

Q2. Regarding the $L(s)$ in (3.1).

$L(s)$ is a user-chosen scaling factor which shapes our SCR for different states $s\in\mathcal{S}$ ($L(s)>0$), especially the width. CIs in OPE consider $L(s)=\sqrt{U_\pi(s)^\top \Lambda_\pi U_\pi(s)}$ since the CLT theorem gaurantees $\sqrt{N}(\hat{V}(\pi;s)-V(\pi;s))/L(s)\to_d N(0,1)$ for fixed $s$.

In our framework, the choice of $L(s)$ can be very flexible ($L(s)=\sqrt{U_\pi(s)^\top \Lambda_\pi U_\pi(s)}$ or simply $L(s)=1$). This is achieved by applying a stronger convex Gaussian approximation theorem (Theorem 3.1) instead of the simple CLT (the detailed explanation is put in response Q6&W3).

Q3. Regarding the accessibility of the offline dataset for reliable results.

Our theorems ensure that the SCRs are reliable when the sample size $N$ is sufficiently large. A large offline dataset is accessible in many real-world settings since

• Offline RL leverages existing datasets (e.g., possibly historical logs, human demonstrations, or prior experiments).
• Many domains (healthcare, education) have abundant offline data.

Additionally, our theoretical reliability is built upon our Assumptions, which are mild for most cases and no stronger than those in serial theoretical work on offline RL (detailed in response Q5&W1). Thank you again for your comment, and the discussion on the access to offline data will be included in our revised paper.

Q4. Regarding the derivation of (3.13).

On the one hand, we write the probability in the first line of (3.13) by a more specific formulation:

On the other hand, the second line in (3.13) can be formulated as

Note that $L(s)>0$ and one can verify

thus, the derivation in (3.13) is correct.

Q5&W1. Regarding the intuition and rationality for Assumptions A1-A3.

Assumptions A1-A3 are very mild assumptions and serve as the minimal requirement for the goodness of the offline dataset to support feasible evaluation.

• Assumption (A1) requires the density $\mu(s)$ for the invariant distribution to be uniformly bounded away from $0$ and $\infty$ so that the Markov chain would not be trapped in a small subset of the entire space. The density $\nu_0(s)$ for the initial state should also follow the boundary condition so that every state is possible to be the initial state.
• Assumption (A2) relaxes the condition on sample size. Previous work (e.g., [5]) requires the number of trajectories $n\to\infty$. (A2) additionally allows fixed $n$, but length $T\to \infty$ when the action variety is sufficiently large on each chain. The geometrical ergodic for the Markov chain is technical for theoretical deduction and commonly assumed for deriving limit theory.
• Assumption (A3) requires the reward signal diversity. The $\omega_\pi(s,a)\geq c_0^{-1}$ requires that the reward random variable is nondegenerate (not always the same). $P(\max_{0\leq t\leq T-1}|R_{0,t}|\leq c_0)=1$ means the rewards are bounded.

These assumptions are the same as those needed to establish pointwise CIs in OPE as discussed in Remark 3.1.

Q6&W3. Regarding the discussion of Theorem 3.1.

Theorem 3.1 is our main theorem, which provides a higher convex Gaussian approximation towards $\hat{\theta}\_\pi=\hat{\beta}\_\pi-\hat{\beta}\_\pi\^*$ in the OPE estimation error $\hat{V}(\pi;s)-V(\pi;s)=\Phi(s)^\top \sum_{a\in\mathcal{A}}\pi(a|s)\hat{\theta}\_\pi$. This is a different approach from the previous pointwise CI results obtained by the simple CLT. The CLT focuses on deriving a limiting theorem on the distribution function of inner product value $\Phi(s)^\top \sum_{a\in\mathcal{A}}\pi(a|s)\hat{\theta}\_\pi$ for fixed $s\in\mathcal{S}$. The generalization from fixed $s\in\mathcal{S}$ to arbitrary $s\in\mathcal{S}$ is not straightforward, since one needs to control

where $\mathbf Z\_\pi$ is a random Gaussian vector. Note that for both scenarios of finite collection of $s$ and all $s\in \mathcal S$, the above can be rewritten as $\Delta\_{\mathbb O}=\left|\mathrm{P}\left(\sqrt{N} \hat{\theta}\_\pi \in \mathbb{O}\right)-\mathrm{P}\left(\mathbf{Z}\_\pi \in \mathbb{O}\right)\right|$ for some covex region $\mathbb O$ (see Remark 3.2 for details). When considering finite sample properties, we provide a bound concerning sample size $N$ and basis number $K$ from the proof of Theorem 3.1 (details in (S.37), (S.38), and (S.39) in the supplement),

In the proof of Theorem 3.1, we also verify that the bound can converge to $0$ as $N\to \infty$ when $K$ is smaller than $N^{2/7}$. In this way we address the space-time issues in the Gaussian process in Theorem 3.1. Our convex Gaussian approximation focuses on the probability of the random vector $\hat{\theta}_\pi$ on all convex sets in $\mathbb{R}^{mK}$. We will add a new remark in revision for the explanation.

W3, part II. Regarding the experimental demonstrations.

Due to the page limit, please kindly see our Q2bW2 response to Reviewer MvGP for our demonstration on OhioT1DM experiment.

Minor. Regarding Assumption on Eq (2.2).

Eq (2.2) assumes the conditional mean of the reward $R_{0,t}$ is independent of $(S\_{0,k},R\_{0,k},A\_{0,k})\_{k<t}$, i.e.

for some reward function $r$.

[1] Fujimoto S, Gu S S. A minimalist approach to offline reinforcement learning.

[2] Levine S, Kumar A, Tucker G, et al. Offline reinforcement learning: Tutorial, review, and perspectives on open problems.

[3] Prudencio R F, Maximo M R O A, Colombini E L. A survey on offline reinforcement learning: Taxonomy, review, and open problems.

[4] Shi, C. (2025). Statistical inference in reinforcement learning: A selective survey.

[5] Jiang, N. Li, L. (2016) Doubly robust off-policy value evaluation for reinforcement learning.

Thank you for acknowledging our clarifications. We respectfully emphasize that offline RL was specifically created to take advantage of large historical datasets. This isn't simply an added feature, but rather represents a core objective of the methodology. The method exists precisely for situations where we can't do online interaction (like in healthcare due to safety concerns) but do have existing data to learn from.

Our empirical results demonstrate reliable performance of our large sample theory with dataset sizes around N=2000, and outperform modified classic methods (such as SAVE) when the sample size is larger. This scale is actually quite modest compared to typical offline RL applications. For example, the HiRID [1] and MIMIC-III [2] datasets both contain extensive ICU records collected over multiple years, with $N>10^8$. These datasets are widely used in medical AI research and are notable for their large volume and richness in temporal resolution. In the field of computer science, there are also many large-scale datasets available, see for example, D4RL (Datasets for Deep Data-Driven Reinforcement Learning), where most datasets have $N > 10^8$. Notice that the above are examples of publicly available data sets. Therefore, our approach can be reliably applied to evaluate policies across a wide range of scenarios.

From a theoretical perspective, we understand your concern about sample size may stem from the fact that traditional asymptotic frameworks for evaluating maximum deviation lead to a Gumbel distribution (from extreme value theory), which exhibits slow convergence rates. Therefore, we stress that our theoretical results are built upon non-asymptotic results, which do not involve any convergence from extreme value theory in statistics. In literature, even if the explicit critical value can be obtained from the limiting distribution derived from extreme value theory, one might still prefer bootstrap inference for better finite sample performance. Similarly, in our case, the bootstrap procedure for widths of the proposed SCR in our methods results in good performance under moderate sample size, as supported by simulations.

We hope this evidence reassures you that our dataset requirements are not only realistic but also conservative compared to common practice.

[1] HiRID: Faltys, M., Zimmermann, M., Lyu, X., Hüser, M., Hyland, S., Rätsch, G., & Merz, T. (2021). HiRID, a high time-resolution ICU dataset (version 1.1.1). PhysioNet. RRID:SCR_007345.

[2] MIMIC-III: Johnson, A., Pollard, T., & Mark, R. (2016). MIMIC-III Clinical Database (version 1.4). PhysioNet. RRID:SCR_007345.

We sincerely appreciate your thorough review and insightful feedback on our manuscript. Thank you for recognizing the significance of our work in advancing uncertainty quantification for OPE in RL. We are particularly grateful for your acknowledgment of our framework's contributions regarding its efficiency, practical applicability, and flexibility. We have carefully considered all of your comments. In what follows, we address each of the questions you raised in detail.

W1. Regarding the constraint on $K$.

The condition $K=O(N^{2/7}\log^{-1}N)$ is imposed to ensure the convex Gaussian approximation error in Theorem 3.1 can attain $o(1)$. In contrast to the pointwise CLT ($K=O(N^{1/2}\log^{-1}N)$), this condition narrows the range of $K$ but does not lead to conflict in most cases.

Take a simple case where the state space $\mathcal{S}$ is one-dimensional$(d=1)$ and the true Q-function lies in a H"older space with smoothness $p=2$. On the one hand, the mathematical approximation theorem ensures approximation accuracy for the Q-function is of rate $\epsilon_K=O(K^{-2})$ (details included in Section D.1 of the appendix). On the other hand, our Proposition 3.1 shows our SCR width is of rate $O(N^{-1/2}\log^{1/2} N)$. Therefore, we wish the $K$ is sufficiently large so that $\epsilon_K$ is negligible compared to our SCR, i.e., $\epsilon_K\ll N^{-1/2}\log^{1/2} N$. This yields $K\gg N^{1/4}\log^{-1/4} N$ and does not conflict with $K=O(N^{2/7}\log^{-1}N)$ so that we can find a feasible $K=N^{\alpha}$ for some $\alpha\in (1/4,2/7)$ to support valid asymptotic inference.

For the high-dimensional setting where $d\to \infty$, the curse of dimensionality is involved in the asymptotic rate. To the best of our knowledge, constructing confidence intervals/regions for the high-dimensional data in off-policy evaluation is still an open problem, which could be a future work.

W2. Regarding the geometric properties of the chosen basis functions.

We appreciate your insightful observation. Indeed, the validity of the logarithmic rate for SCR width relies on the geometric properties specified in Condition (3.15). While these properties can be rigorously verified for commonly used basis functions (e.g., Legendre polynomials, splines, and Fourier series in Section D of the supplement), extending the verification to arbitrary or highly complex bases may require additional substantial theoretical work.

In practice, most applications employ such standard bases due to their well-characterized properties and computational tractability. Therefore, for the current analysis, we focus on establishing guarantees for those widely adopted bases that possess the geometric properties, and leave the generalizations to broader function classes as future work.

W3. Regarding the exploration of finite-sample behaviors under various data characteristics.

Thank you for your question. Taking Scenario 2 (i.e., Table 1) as an illustrative example, we conducted additional simulations for both smaller and larger values of $N$. We kindly refer you to our second response to Reviewer MvGP for the results. Moreover, for high noise and non-Gaussian errors, we modified the state transition rule, specifically the $z_{0,t}$ in line 268, to be an i.i.d. two-dimensional $t(8)$ random variable, while keeping the rest unchanged. The results are shown in Table R2. These results demonstrate that for noise with heavier tails, our coverage and length remain robust.

Table R2: Results for Scenario 2 with $t(8)$ noises across 500 repetitions.. Format: ECP(AL).

W4. Regarding comparisons against baseline other than Bonferroni correction.

Thank you for your question. In addition to the Bonferroni correction, we also compared our method with SAVE using Sidak correction [1], and the results are presented in Table R1 (we also paste it below for your convenience), and Table R2 in the previous question. The results show that, while ensuring coverage (ECP), our method achieves the shortest averaged length (AL).

Table R1: Results for Scenario 2 across 500 repetitions.. Format: ECP(AL).

[1] Abdi, Hervé. "Bonferroni and Šidák corrections for multiple comparisons." Encyclopedia of measurement and statistics 3.01 (2007): 2007.

Thank the authors for their detailed response. The additional experiment and other response did help clarify the paper. Therefore, I would like to keep my score.

We sincerely appreciate your careful review of our manuscript and your thoughtful feedback on our work. Thank you for remarking on the clarity of our presentation, theoretical rigor, and experimental validation, as well as for your kind recognition of the novelty of our method in addressing a key gap in the literature. We are also grateful for your attention to detail in identifying typos. Your acknowledgment of the soundness of our methodology and the relevance of our empirical results is greatly encouraging.

W1. Regarding offline RL limitation.

We acknowledge your observation that our method is currently limited to offline RL. This focus aligns with our goal of addressing critical challenges in off-policy evaluation, particularly in safety-sensitive domains where online exploration is costly or impractical (e.g., healthcare, autonomous systems).

W2. Regarding code for implementation.

Thank you for your very nice question! We have set up a github repository to make the code public, however, due to the rebuttal regulations, we cannot give out links now. We will put a link to the code in the revised version. The datasets used in our experiments are publicly available and cited in the manuscript. We will also incorporate the data in our code to ensure good reproducibility.

W3. Regarding computational cost and technical constraints.

Thank you for your question. In terms of computation, our procedure performs fast and can be completed on a personal laptop. For example, all simulations in this paper were run on an Apple M1 machine with 16GB of RAM and macOS Sonoma. From a theoretical perspective, the term $\hat{\beta}_{\pi}$ in (3.6) is similar to a least squares estimation, which can be efficiently computed. Additionally, Steps 1 and 4 of the bootstrap procedure are both linear, as we use a linear approximation method. Therefore, the computations are efficient. When higher precision is required, we can choose a larger $B$ (bootstrap size) and employ parallel computing. Specifically, the time complexity of our method is $O(NK^2+BK)$ and the space complexity is $O(N+BK)$. This also addresses Question 1 and will be included in our revised paper.

Minor. Regarding typos

Thank you and we have corrected all the typos you mentioned in the revision.

[1] Abdi, Hervé. "Bonferroni and Šidák corrections for multiple comparisons." Encyclopedia of measurement and statistics 3.01 (2007): 2007.

We sincerely appreciate your review of our manuscript and your insightful feedback. We are grateful for your recognition of the interest and soundness of our work. In the following, we will address the questions you raised.

Q1. Regarding consolidating assumptions in Section 3

We appreciate your helpful suggestion for improving our paper’s readability. We agree that presenting a consolidated list of assumptions at the beginning of Section 3, which includes policy stationarity and other key conditions, would significantly improve the paper's clarity and help contextualize the theoretical results. While we cannot modify the current NeurIPS submission due to venue constraints, we will incorporate this change in the camera-ready version if accepted, and in any extended/publicly archived versions of the work.

Q2a&W1,W3. Regarding more comprehensive evaluation on synthetic settings.

Thank you for your question. Following your comments in weakness 1 and 3, We have conducted a broader range of additional simulations based on the settings in Table 1, and the results are shown in Table below. In addition to the original SAVE method with Bonferroni correction compared in the paper, we also present results using SAVE with Sidak correction [1]. The results indicate that the coverage of our method remains around 0.95. moreover, as $n \times T$ increases, the averaged length (AL) generally shows a decreasing trend.

Table R1: Results for Scenario 2 across 500 repetitions. Format: ECP(AL).

Q2b&W2. Regarding the real data application on OhioT1DM dataset.

Thank you for your question. To clarify, we do not pre-assume that double FQI policy ($\pi^{opt}$) outperforms random forest policy ($b$). Instead, our key contribution lies in proposing a simultaneous framework to compare $\pi^{opt}$ and $b$ across all different initial states, whereas prior works (e.g., [3]) could only assess them pointwisely. Figure 2 in our paper illustrates that, when using $\pi^{opt}$, the number of states where the value function is significantly non-zero, which indicates that glucose levels outside the normal range (either hypoglycemic or hyperglycemic) are fewer, suggesting that $\pi^{opt}$ performs better. This conclusion cannot be drawn from pointwise comparison.

We sincerely thank you for your insightful suggestion to compare the lengths in real data analysis. The true values of the value function under both policies are unknown; however, we can report the difference in average length of the confidence band at the same significance level. Taking $V(\pi^{opt}, S_0)$ as an example, our method yields an averaged length of 27.0, while SAVE with Bonferroni correction produces an AL of 32.4, which is 20% larger than ours. We will include these results in the revised version of the paper. This reply also addresses the weakness 2.

Q2c. Regarding more baselines other than SAVE.

Thank you for your question. We would like to clarify that SAVE is a direct method, and as such, we focus on comparing our method with IS. To illustrate, we performed additional simulations by modifying Scenario 2. Specifically, we set $S_0 = (-2 + 0.2i, -2 + 0.2j)^\top$, where $0 \leq i, j \leq 20$. For each combination of $i$ and $j$, we generated 10 trajectories, each of length 10. All other settings remained unchanged. For bootstrapping IS, we used the bootstrap algorithm from [2], which provides confidence intervals for each $V(\pi, S_0)$. We then applied the Bonferroni correction to obtain the SCB. We present the results in figure (similar to Figure 1 in the main article), and it will be included in the revised version. The results show that the ECP of IS is mostly conservative, although some coverage rates fall below 0.95. As for AL, it is on average longer than our method. These results demonstrate that our approach offers certain advantages over IS. This reply also addresses the weakness 3.

[1] Abdi, Hervé. Bonferroni and Šidák corrections for multiple comparisons. Encyclopedia of measurement and statistics 3.01 (2007): 2007.

[2] Hanna, Josiah, Peter Stone, and Scott Niekum. Bootstrapping with models: Confidence intervals for off-policy evaluation. Proceedings of the AAAI Conference on Artificial Intelligence. Vol. 31. No. 1. 2017.

[3] Shi, C., Zhang, S., Lu, W., & Song, R. (2022). Statistical inference of the value function for reinforcement learning in infinite-horizon settings. Journal of the Royal Statistical Society Series B: Statistical Methodology, 84(3), 765-793.