# Off-Policy Risk Assessment in Markov Decision Processes

Audrey Huang  Department of Computer Science, University of Illinois Urbana-
Champaign  Liu Leqi  Machine Learning Department, Carnegie Mellon University
Zachary C. Lipton  Machine Learning Department, Carnegie Mellon University
Kamyar Azizzadenesheli  Department of Computer Science, Purdue University

###### Abstract

Addressing such diverse ends as safety alignment with human preferences, and
the efficiency of learning, a growing line of reinforcement learning research
focuses on risk functionals that depend on the entire distribution of returns.
Recent work on _off-policy risk assessment_ (OPRA) for contextual bandits
introduced consistent estimators for the target policy’s CDF of returns along
with finite sample guarantees that extend to (and hold simultaneously over)
all risk. In this paper, we lift OPRA to Markov decision processes (MDPs),
where importance sampling (IS) CDF estimators suffer high variance on longer
trajectories due to small effective sample size. To mitigate these problems,
we incorporate model-based estimation to develop the first doubly robust (DR)
estimator for the CDF of returns in MDPs. This estimator enjoys significantly
less variance and, when the model is well specified, achieves the Cramer-Rao
variance lower bound. Moreover, for many risk functionals, the downstream
estimates enjoy both lower bias and lower variance. Additionally, we derive
the first minimax lower bounds for off-policy CDF and risk estimation, which
match our error bounds up to a constant factor. Finally, we demonstrate the
precision of our DR CDF estimates experimentally on several different
environments.

††††Published at the 25th International Conference on Artificial Intelligence
and Statistics (AISTATS) 2022.

##  1 Introduction

In off-policy evaluation (OPE), we aim to estimate the risk associated with a
_target_ policy given only offline datasets collected from a deployed
_behavior_ policy. Broadly, a risk functional maps from _the distribution of
returns_ to a real value. While OPE research has historically focused on
estimating the expected return, a wave of recent works in reinforcement
learning (RL) focus on other risk functionals, motivated by such diverse
desiderata as risk aversion (and seeking) and alignment with human
preferences. Examples of such alternative risk functionals include the
expectation, variance, value at risk, conditional value at risk, and broader
families such as distortion risk functionals, coherent risk functionals, and
Lipschitz risk functionals, Notably, the Lipschitz risks, introduced by Huang
et al. (2021b) subsume most common risk functionals.

In a new line of work on _off-policy risk assessment_ (OPRA), researchers have
developed methods for simultaneously estimating arbitrarily many risk
functionals, providing guarantees that hold simultaneously over an entire
collection of risks (Huang et al., 2021b; Chandak et al., 2021a). These
results are achieved through a two-step process: (i) estimate the cumulative
distribution function (CDF) of returns under the target policy; and (ii)
calculate the risk functionals of interest on the estimated CDF. This approach
has two advantages: (a) guarantees on the CDF estimation error yield
corresponding guarantees on the risk estimates; and (b) because the risk
estimates all share the same underlying CDF, the corresponding guarantees on
any number of downstream risk estimates hold simultaneously without loss in
statistical power.

To our knowledge, for MDPs, importance sampling (IS) is the only established
method for estimating off-policy CDFs (Chandak et al., 2021a). However, IS
estimators are known to suffer from high variance, render them impractical
over long trajectories. This variance can be exponential in the horizon (Liu
et al., 2018), as weights vanish for most trajectories and explode for others,
yielding tiny effective sample sizes. Consider, for example, that
practitioners often evaluate a deterministic target policy given a dataset
collected by a stochastic behavioral policy. In such cases, the IS weight is
zero for all trajectories where the target policy places zero probability on
_any_ action that the behavioral policy takes, even if the policies agree on
most actions. This renders most trajectories unusable, shrinking the effective
size of an offline dataset (Sussex et al., 2018).

In off-policy evaluation addressing the standard _expected return_ risk
functional, doubly robust estimators and their variants have been shown to
mitigate these problems (Dudik et al., 2011; Dudík et al., 2014; Jiang and Li,
2016; Thomas and Brunskill, 2016; Wang et al., 2017; Su et al., 2020). In the
contextual bandit setting, the paper introducing OPRA also proposed and
analyzed a doubly robust estimator for the CDF (Huang et al., 2021b), but
model-based and doubly robust risk estimation in MDPs, where it may be most
impactful, remains an open problem. One primary challenge is that, in the
expected return setting, extending the model-based and doubly robust
estimators from contextual bandits to MDPs (Jiang and Li, 2016) relied on
stepwise or recursive formulations of the importance sampling estimator, which
lack clear analogs in CDF estimation.

In this paper, we define and analyze the first doubly robust CDF estimator for
off-policy risk assessment in MDPs. Specifically, we contribute the following:

  1. 1.

A recursive formulation of the CDF of returns in MDPs (§ 3.2), which enables
us to derive importance sampling (§ 4), and model-based CDF estimators (§ 5).

  2. 2.

Rate-matching upper and minimax lower bounds (§ 4.2) for the IS CDF estimator.

  3. 3.

Variance-reduced IS estimators leveraging unique properties of the return
distribution (§ 4.4).

  4. 4.

The first doubly robust off-policy estimator for the CDF of returns (§ 6), and
analysis demonstrating its potential for significant variance reduction.

  5. 5.

The first Cramer-Rao lower bound for off-policy CDF estimation and analysis
demonstrating that the doubly robust estimator achieves this lower bound in
certain cases (§ 6.3).

  6. 6.

A demonstration that the risk estimates inherit the reduced variance and error
in CDF estimation (§ 7).

  7. 7.

Experimental evidence that the doubly robust CDF and risk estimates exhibit
significantly lower error and variance than importance sampling estimators on
benchmark problems, with applications to safe policy evaluation and
improvement (§ 8).

##  2 Related Work

Several estimators have been studied for off-policy estimation of the expected
return. A seminal work by Dudik et al. (2011) proposed the doubly robust
estimator for the expected return in contextual bandits, which was extended to
the MDP setting by Jiang and Li (2016). The doubly robust estimator leverages
both importance sampling weights and models to reduce variance while retaining
unbiasedness. Researchers have proposed several variants of the doubly robust
estimator, aiming to improve its performance, with Thomas and Brunskill
(2016); Wang et al. (2017) optimizing the balance between importance sampling
and model-based estimators and Su et al. (2020) optimizing the weights used in
the estimator. A Cramer-Rao variance lower bound (Jiang and Li, 2016) and
minimax lower bounds (Li et al., 2015; Wang et al., 2017; Ma et al., 2021) for
off-policy estimation of the mean have been derived.

While the expected return has been the focus of most previous off-policy
evaluation research, other risk functionals are increasingly of interest. In
the RL literature, the variance (Sani et al., 2013; Tamar et al., 2016),
variance-constrained mean (Mannor and Tsitsiklis, 2011), and conditional
value-at-risk (Rockafellar et al., 2000; Keramati et al., 2020; Huang et al.,
2021a) are often employed. Other risk functionals, including exponential
utility (Denardo et al., 2007), distortion risk functionals (Dabney et al.,
2018), and cumulative prospect weighting (Gopalan et al., 2017; Prashanth et
al., 2016, 2020) have also been investigated.

Several recent works have aimed to estimate risk functionals from off-policy
datasets. Of these, Chandak et al. (2021b) estimates the variance, while more
recent works (Huang et al., 2021b; Chandak et al., 2021a) tackle the
estimation of more general risks and are the closest works of comparison. Both
Huang et al. (2021b); Chandak et al. (2021a) take a two-step approach of first
estimating the off-policy CDF of returns; and then estimating their risks via
a plug-in approach. Chandak et al. (2021a) proposed an IS CDF estimator in
both stationary and nonstationary MDPs and derives confidence intervals for a
number of risk estimates on a per-risk basis. Huang et al. (2021b) proposed
both IS and DR estimators for the CDF of returns in contextual bandits,
establishing the first finite sample guarantees for off-policy CDF estimation
and showing that the estimated CDFs converge uniformly, achieving
$O(1/\sqrt{n})$ rates. Additionally, they introduced the broad class of
Lipschitz risk functionals, where the errors in risk estimation are upper-
bounded by the error in CDF estimation multiplied by each functional’s
Lipschitz constant.

While Chandak et al. (2021a) proposes _only_ IS CDF estimates for MDPs, Huang
et al. (2021b) proposes variance-reduced DR CDF estimates in _only_ contextual
bandits. Thus _variance-reduced risk estimation in MDPs_ , where high variance
from IS weights is most problematic (Liu et al., 2018), is an open problem.
Towards this end, we (1) develop provably variance-reduced CDF and risk
estimators in MDPs, and, (2) demonstrate their optimality by deriving the
first matching upper and lower bounds.

##  3 Preliminaries

###  3.1 Problem Setting

An MDP is defined by $\mathcal{M}=(\mathcal{S},\mathcal{A},P,R,\gamma,\mu)$,
where $\mathcal{S}$ is the state space, $\mathcal{A}$ is the action space,
$P:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow\mathbb{R}$ is the
transition function, $R:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$ is
the reward function, $\gamma\in(0,1)$ is the discount factor, and $\mu$ is the
fixed starting state distribution. We assume that rewards are bounded on the
support $[0,D]$. A stationary policy
$\pi:\mathcal{S}\rightarrow\Delta(\mathcal{A})$ maps a state to a probability
distribution over actions.

In the off-policy evaluation problem, we are concerned with estimating risk
functionals of a policy $\pi$ given a dataset of trajectories
$\mathcal{D}=\\{\tau^{i}\\}_{i=1}^{n}$ collected by the behavioral policy
$\beta$ interacting with MDP $\mathcal{M}$. Each
$\tau=(S_{1},A_{1},R_{1},\ldots,S_{H},A_{H},R_{H},S_{H+1})$ is an $H$-step
trajectory, with $S_{1}\sim\mu$, $R_{h}\sim R(S_{h},A_{h})$, and $S_{h+1}\sim
P(\cdot|S_{h},A_{h})$. In general, we will refer to random variables with
capital letters, and for convenience, we denote
$\tau_{j}:=(S_{1},A_{1},R_{1},\ldots S_{j},A_{j},R_{j})$, or a sub-trajectory
of $\tau$ from horizon $1$ to $j$.

Further, let $w(A_{h},S_{h}):=\frac{\pi(A_{h}|S_{h})}{\beta(A_{h}|S_{h})}$
denote the importance weight, and let
$w_{max}:=\max_{h,A_{h},S_{h}}w(A_{h},S_{h})$ be the maximum importance weight
over all horizons, states, and actions. For convenience, by
$w_{j}:=\prod_{h=1}^{j}w(A_{h},S_{h})$, we denote the importance weight of
subtrajectory $\tau_{j}$. Thus, $w_{H}$ is the importance weight of the
trajectory $\tau$.

We denote by $\mathbb{P}$ the distribution of trajectories induced by $\pi$ on
$\mathcal{M}$, and by $\mathbb{P}_{\beta}$ the distribution induced by
$\beta$. The random variable of returns of an $H$-step trajectory induced by
$\pi$ in $\mathcal{M}$ is given by
$Z^{\pi}=\left(\sum_{h=1}^{H}\gamma^{h-1}R_{h}\Big{|}\pi\right)$. We also
write the random variable of returns starting from a state $s_{h}$ at horizon
$h$, summed until the end of horizon $H$, to be
$Z^{\pi}_{s_{h}}=\left(\sum_{k=h}^{H}\gamma^{k-h}R_{k}\Big{|}\pi,s_{h}\right)$.
The random variable of returns conditioned on a state and action at horizon
$h$, $Z^{\pi}_{s_{h},a_{h}}$, is defined similarly.

The cumulative distribution function (CDF) of $Z^{\pi}$ is
$F^{\pi}(t)=\mathbb{E}_{\mathbb{P}}\left[\mathbbm{1}\left\\{Z^{\pi}\leq
t\right\\}\right]$, where $\mathbbm{1}\left\\{\cdot\right\\}$ is the indicator
function. Similarly, the CDF of conditional returns $Z_{s_{h}}^{\pi}$ is
defined as
$F^{\pi}_{s_{h}}(t)=\mathbb{E}_{\mathbb{P}}\left[\mathbbm{1}\left\\{Z^{\pi}_{s_{h}}\leq
t\right\\}|s_{h}\right]$, and the state-action dependent distribution
$F_{s_{h},a_{h}}(t)$ is defined similarly. We also define the _complementary
cumulative distribution function_ (CCDF)
$S^{\pi}(t):=1-F^{\pi}(t)=\mathbb{E}_{\mathbb{P}}\left[\mathbbm{1}\left\\{Z^{\pi}>t\right\\}\right]$,
and $S_{s_{h}}^{\pi}(t):=1-F_{s_{h}}^{\pi}(t)$.

For short, we write
$\mathbb{E}_{h}^{\beta}[\cdot]:=\mathbb{E}_{\mathbb{P}_{\beta}}[\cdot|\tau_{h-1}]$
to be the expectation under behavioral policy $\beta$, conditioned on a
trajectory up until time $h$. Similarly, for $\pi$, we write the conditional
expectation
$\mathbb{E}_{h}[\cdot]:=\mathbb{E}_{\mathbb{P}}[\cdot|\tau_{h-1}]$. We define
the conditional variances $\mathbb{V}_{h}^{\beta},\mathbb{V}_{h}$ similarly,
e.g.,
$\mathbb{V}_{h}^{\beta}=\mathbb{V}_{\mathbb{P}_{\beta}}[\cdot|\tau_{h-1}]$.

###  3.2 CDF Bellman Equation

Our off-policy CDF estimators rely on a novel _CDF Bellman equation_ that
provides a recursive formulation of the CDF of returns $F$. In the off-policy
setting, for state $s_{h}$ at horizon $h$, we have:

$\displaystyle F_{s_{h}}(t)$
$\displaystyle:=\mathbb{E}_{\mathbb{P}_{\beta}}\left[w(A_{h},s_{h})F_{S_{h+1}}\left(\frac{t-R_{h}}{\gamma}\right)\Big{|}s_{h}\right]$
(1)

where $F_{s_{H+1}}=\mathbbm{1}\\{0\leq t\\}$ for all $s_{H+1}$ at the end of
the horizon. Note that the CDF of returns under $\pi$ can also be written as
$F^{\pi}(t)=\mathbb{E}_{s_{0}\sim\mu}[F_{s_{0}}^{\pi}(t)]$, i.e. the
expectation under the starting state distribution $s_{0}\sim\mu$ of the CDFs
$F^{\pi}_{s_{0}}$.

This recursion was first introduced in Sobel (1982) for deterministic rewards
in the value iteration setting, which we have extended to the more general
stochastic setting. Using the definition of equivalence in distribution
($\stackrel{{\scriptstyle D}}{{=}}$), it can be seen that the CDF Bellman
operator is equivalent to the distributional Bellman operator (Bellemare et
al., 2017; Dabney et al., 2018), $Z_{s_{h}}\stackrel{{\scriptstyle
D}}{{=}}R(s_{h},A_{h})+\gamma Z_{S_{h+1}}$.

###  3.3 Risk Functionals

Let $Z\in\mathcal{L}_{\infty}(\Omega,\mathcal{F}_{Z},\mathbb{P}_{Z})$ denote a
real-valued random variable that admits a CDF
$F_{Z}\in\mathcal{L}_{\infty}(\mathbb{R},\mathcal{B}(\mathbb{R}))$. A _risk
functional_
$\rho:\mathcal{L}_{\infty}(\Omega,\mathcal{F}_{Z},\mathbb{P}_{Z})\rightarrow\mathbb{R}$
is said to be law-invariant if $\rho(Z)$ depends only on the distribution of
$Z$ (Kusuoka, 2001), i.e., if for any pair of random variables $Z$ and
$Z^{\prime}$,
$F_{Z}=F_{Z^{\prime}}\;\Longrightarrow\;\rho(Z)=\rho(Z^{\prime})$. When clear
from the context, we write $\rho(F_{Z})$ in place of $\rho(Z)$. In this paper,
we restrict our focus to law invariant risk functionals, as it may not be
practical to estimate risk functionals that are not law invariant from data
(Balbás et al., 2009).

Examples of popularly studied law-invariant risk functionals include the
expected value, variance, mean-variance, coherent risks (including conditional
value-at-risk (CVaR) and proportional hazard), distortion risks, and
cumulative prospect theory risks, which we discuss more extensively in
Appendix A due to space constraints. We give particular attention to
distortion risk functionals due to their widespread usage and general form.
Distortion risks apply a distortion function $g:[0,1]\to[0,1]$ to the CDF
$F_{Z}$ and have the formula

$\rho(Z)=\int_{0}^{\infty}g(1-F_{Z}(t))dt=\int_{0}^{\infty}g(S(t))dt.$

Notably, expected value is a special case of distortion risk with $g(x)=x$,
while CVaR at level $\alpha$ ($\text{CVaR}_{\alpha}$) has a distortion
function $g(x)=\min\\{\frac{x}{1-\alpha},1\\}$. A large number of the
aforementioned risk functionals satisfy a notion of smoothness in the CDF,
which Huang et al. (2021b) define as the class of _Lipschitz risk functionals_
:

######  Definition 3.1 (Lipschitz Risk Functionals).

A law invariant risk functional $\rho$ is $L$-Lipschitz if there exists
$L\in[0,\infty)$ such that for any pair of CDFs $F_{Z}$ and $F_{Z^{\prime}}$,
it satisfies

$|\rho(F_{Z})-\rho(F_{Z^{\prime}})|\leq L\|F_{Z}-F_{Z^{\prime}}\|_{\infty}.$
(2)

Lipschitz risk functionals include distortion risk functionals with Lipschitz
distortion functions $g$. As special cases, expected value is $D$-Lipschitz
and CVaR at level $\alpha$ is $D/(1-\alpha)$-Lipschitz. Additionally, variance
is ($3D^{2}$-)Lipschitz and so is mean-variance (Huang et al., 2021b).

As we analyze the theoretical properties of CDF and risk estimators later in
our paper, we show that estimates of Lipschitz risk functionals retain many
desirable properties, including characterizable finite sample error and
consistency.

###  3.4 Off-Policy Risk Assessment

In the remainder of this paper, we will take the general approach of Algorithm
1 (OPRA) for off-policy risk assessment. OPRA takes as input a set of risk
functionals of interest, and outputs estimates of their values under the
target policy $\pi$. It first estimates the CDF of returns under $\pi$ from
off-policy data, then generates plug-in risk estimates on the CDF. In Sections
4, 5 and 6, we define and analyze CDF estimators. In Section 7, we turn to
plug-in risk estimation.

Input: Dataset $\mathcal{D}$, policy $\pi$, probability $\delta$, risk
functionals $\left\\{\rho_{p}\right\\}_{p=1}^{P}$ with Lipschitz constants
$\left\\{L_{p}\right\\}_{p=1}^{P}$.

1 Estimate the CDF $\widehat{F}$ with error $\epsilon$;

2 For $p=1...P$, estimate risk $\widehat{\rho}_{p}=\rho_{p}(\widehat{F})$;

Output: Estimates with errors $\\{\widehat{\rho}_{p}\pm
L_{p}\epsilon\\}_{p=1}^{P}$.

Algorithm 1 Off-Policy Risk Assessment (OPRA)

##  4 Importance Sampling CDF Estimation

Off-policy evaluation faces the unique challenge that only rewards for actions
taken by the behavioral policy are observed. Importance sampling (IS), which
reweights samples according to the ratio between target and behavioral policy
probabilities, can be used to account for this distribution shift. While the
IS estimator for off-policy estimation was previously introduced in Huang et
al. (2021b); Chandak et al. (2021a), finite sample convergence in MDPs and
lower bounds for off-policy CDF estimation (in general) remained open
questions. In this section, we derive the first minimax lower bound for off-
policy CDF estimation, and show that it matches the IS estimator upper bound
up to a constant factor. Further, we leverage unique properties of the return
distribution to develop new variance-reduced IS estimators.

###  4.1 Importance Sampling Estimator

The IS CDF estimator forms an empirical estimate of the recursion in (1):

$\displaystyle\widehat{F}_{\textsc{IS}}(t)=\frac{1}{n}\sum_{i=1}^{n}\widehat{F}_{s_{1}^{i}}(t),$
(3)
$\displaystyle\quad\text{where}\quad\widehat{F}_{s_{h}^{i}}(t)=w(a_{h}^{i},s_{h}^{i})\widehat{F}_{s_{h+1}^{i}}\left(\frac{t-r_{h}^{i}}{\gamma}\right)$

At each horizon $h$, the recursive IS estimator applies the importance weight
to the estimated CDF at the next state $s_{h+1}$ convolved with the reward
$r_{h}$ and discount factor $\gamma$. Unrolling the recursion recovers the
more intuitive form of the IS estimator, also introduced in Chandak et al.
(2021a), that reweights each sample by the importance weight of the entire
trajectory:

$\widehat{F}_{\textsc{IS}}(t)=\frac{1}{n}\sum_{i=1}^{n}w_{H}^{i}\mathbbm{1}\left\\{\sum_{h=1}^{H}\gamma^{h-1}r_{h}^{i}\leq
t\right\\}.$ (4)

######  Lemma 4.1.

The IS CDF estimator (3) is unbiased and its variance is recursively given by

$\displaystyle\mathbb{V}_{h}^{\beta}\left[\widehat{F}_{S_{h}}(t)\right]=$
$\displaystyle\mathbb{E}_{h}^{\beta}\left[\mathbb{V}_{h}^{\beta}\left[w(A_{h},S_{h})F_{S_{h},A_{h}}(t)\Big{|}S_{h}\right]\right]$
$\displaystyle+\mathbb{E}_{h}^{\beta}\left[w(A_{h},S_{h})^{2}\mathbb{V}_{h+1}\left[\widehat{F}_{S_{h+1}}\left(\frac{t-R_{h}}{\gamma}\right)\Big{|}S_{h},A_{h}\right]\right]$
$\displaystyle+\mathbb{V}_{h}\left[F_{S_{h}}(t)\right].$ (5)

The first term expresses variance from importance-weighted actions sampled as
$A_{h}\sim\beta(\cdot|S_{h})$, while the third term expresses variance from
transitions. The second term encompasses variance from random rewards and
transitions, and recurses to the next time step. At horizon $h$, each term in
the variance is multiplied by the product of importance weights from the
trajectory thus far, $w_{h-1}$.

###  4.2 Finite Sample Error Bound

While the IS estimator is known to be consistent, obtaining finite sample
error bounds in MDPs has remained an important open problem. Such bounds
characterize the convergence rate of the estimate towards the true $F$, and
provide confidence intervals around estimates in the finite sample regimes of
practical scenarios.

First, however, it is important to note that, given finite samples, the IS
estimator may not be a valid CDF due to the use of importance weighting, e.g.,
it can be greater than 1. To address this issue, we can use a weighted
importance sampling (WIS) estimator instead:
$\widehat{F}_{\textsc{WIS}}(t)=\frac{n}{\sum_{j=1}^{n}w_{H}^{j}}\widehat{F}_{\textsc{IS}}(t)$,
which normalizes the estimate using the sum of importance weights instead of
the sample size, guaranteeing that $\widehat{F}_{\textsc{WIS}}\in[0,1]$. As
Chandak et al. (2021a) demonstrates, the WIS estimator is biased but uniformly
consistent.

Alternatively, the IS estimator can be clipped to the unit range, a strategy
also employed by Huang et al. (2021b) for off-policy CDF estimation in
contextual bandits, shown below:

$\widehat{F}_{\textsc{IS-
clip}}(t)=\min\left\\{\widehat{F}_{\textsc{IS}}(t),1\right\\},$

The IS-clip estimator, too, is biased but uniformly consistent, and its error
is uniformly bounded in the following lemma:

######  Lemma 4.2.

With probability at least $1-\delta$, for universal constants $c_{1},c_{2}>0$

$\displaystyle\|\widehat{F}_{\textsc{IS-clip}}-F\|_{\infty}\leq
c_{1}\sqrt{\frac{\mathbb{E}_{\mathbb{P}_{\beta}}[w_{H}^{2}]}{n}}+c_{2}\frac{w_{max}}{n}.$
(6)

The IS estimator converges at a rate of
$O(\mathbb{E}_{\mathbb{P}_{\beta}}[w_{H}^{2}]/\sqrt{n})$, subsuming previous
results on convergence in contextual bandits (Huang et al., 2021b). For
readability throughout the paper, we leave constants implicit and give their
exact form in the proof (Appendix B).

How close to optimal is this convergence rate? For large enough sample size
$n$, the IS estimator is in fact minimax rate-optimal, as the following lower
bound demonstrates:

######  Theorem 4.1.

For a universal constant $c$, when $n\geq
c\frac{w_{max}^{2}}{\mathbb{E}_{\mathbb{P}_{\beta}}[w_{H}^{2}]}$, we have

$\displaystyle\inf_{\widehat{F}}\sup_{F\in\mathcal{F}}\mathbb{E}\left[\|\widehat{F}-F^{\pi}\|_{\infty}\right]\geq
c\sqrt{\frac{\mathbb{E}_{\mathbb{P}_{\beta}}[w_{H}^{2}]}{n}}$ (7)

###  4.3 Variance-Reduced IS Estimation

Unlike the IS estimator of the expected return (Dudik et al., 2011; Jiang and
Li, 2016), the IS estimator of the CDF has the unique property that its
variance differs across its support as a function of $t$. Examining the
variance for $\widehat{F}_{\textsc{IS}}$ (5) at a fixed $t$, the first term is
small if $F_{S_{h},A_{h}}(t)$ is close to 0, which can occur at small $t$.
However if $F_{S_{h},A_{h}}(t)$ is close to its maximal value of $1$, which
occurs when $t$ is close to the maximum return $D$, the first term is close to
its largest possible value.

####  4.3.1 Estimating the Complementary Cumulative Distribution Function
(CCDF)

Rather than estimating the CDF $F$, another possibility is to estimate the
CCDF $S=1-F$, which is directly used in the expression for many risk
expressions (e.g. distortion risk functionals, see § 3). IS estimators for $S$
take a similar form as those for $F$, but with the indicator function
inequality swapped to $\mathbbm{1}\\{\;\cdot\;>t\\}$. For example,

$\widehat{S}_{\textsc{IS}}(t)=\frac{1}{n}\sum_{i=1}^{n}w_{H}^{i}\mathbbm{1}\left\\{\sum_{h=1}^{H}\gamma^{h-1}r_{h}^{i}>t\right\\},$
(8)

Going forward, we refer to the estimator in (4) as
$\widehat{F}_{\textsc{F-IS}}$, and
$\widehat{F}_{\textsc{S-IS}}:=1-\widehat{S}_{\textsc{IS}}$ from (8) to avoid
confusion. Perhaps non-intuitively, $\widehat{F}_{\textsc{F-IS}}$ and
$\widehat{F}_{\textsc{S-IS}}$ have different theoretical properties, namely
its variance:

$\displaystyle\mathbb{V}_{h}^{\beta}\left[\widehat{F}_{\textsc{S-IS},\;S_{h}}(t)\right]=\;$
$\displaystyle\mathbb{V}_{h}^{\beta}\left[1-\widehat{S}_{S_{h}}(t)\right]=\mathbb{E}_{h}^{\beta}\left[\mathbb{V}_{h}^{\beta}\left[w(A_{h},S_{h})(1-F_{S_{h},A_{h}}(t))\Big{|}S_{h}\right]\right]$
$\displaystyle\quad+\mathbb{E}_{h}^{\beta}\left[w(A_{h},S_{h})^{2}\mathbb{V}_{h+1}\left[\widehat{S}_{S_{h+1}}\left(\frac{t-R_{h}}{\gamma}\right)\Big{|}S_{h},A_{h}\right]\right]$
$\displaystyle\quad+\mathbb{V}_{h}\left[S_{S_{h}}(t)\right]$

The first term of the variance expression for $\widehat{F}_{\textsc{S-IS}}$
takes the exact opposite trend as that of $\widehat{F}_{\textsc{F-IS}}$: it is
low where $F_{S_{h},A_{h}}(t)$ is high, and vice versa. For example, if
$\widehat{F}_{\textsc{F-IS}}$ has low variance in estimating the lower tails
but high variance in the upper tails, $\widehat{F}_{\textsc{S-IS}}$ will have
low variance in the upper tails but high variance in the lower tails. Lower
variance leads to lower error of the CDF estimate on certain portions of the
distribution. To formalize this argument, we present the following bounds on
the error of the CDF estimate in terms of its variance at each $t$. As we
later show, this will also have important consequences for downstream
estimation of different risks, which up-weight different parts of the
distribution.

######  Lemma 4.3.

If $\mathbb{V}[\widehat{F}_{\textsc{IS}}]$ and
$\mathbb{V}[\widehat{S}_{\textsc{IS}}]$ have bounded variation (see Assumption
B.1), there exists universal constants $c_{1},c_{2}$ such that with
probability at least $1-\delta$, $\forall t$,

$\displaystyle|F(t)-\widehat{F}_{\textsc{F-IS}}(t)|$
$\displaystyle\leq\frac{c_{1}\log(\sqrt{n}/\delta)}{n}+\sqrt{\frac{c_{2}\mathbb{V}[\widehat{F}_{\textsc{IS}}(t)]\log(\sqrt{n}/\delta)}{n}},$
$\displaystyle|F(t)-\widehat{F}_{\textsc{S-IS}}(t)|$
$\displaystyle\leq\frac{c_{1}\log(\sqrt{n}/\delta)}{n}+\sqrt{\frac{c_{2}\mathbb{V}[\widehat{S}_{\textsc{IS}}(t)]\log(\sqrt{n}/\delta)}{n}}$

This difference has important implications for downstream risk estimation. The
choice of whether to estimate $S$ or $F$ can reduce the variance of the plug-
in risk estimator, especially for risk functionals such as CVaR that heavily
up-weight one tail of the distribution. In general, under reward
distributions, $\widehat{F}_{\textsc{S-IS}}$ will lead to lower variance
estimators for risk-averse risk functionals that upweight the lower tails of
the distribution. $\widehat{F}_{\textsc{F-IS}}$ will lead to lower-variance
risk estimators for risk-seeking risk functionals that upweight the upper
tails.

The RHS of the error bound in Lemma 4.3 gives a uniform confidence band on the
estimated CDF that contains $F$ with high probability. However, the variance
of $\widehat{F}_{\textsc{F-IS}}$ or $\widehat{F}_{\textsc{S-IS}}$ cannot be
computed without knowledge of the underlying MDP. Instead, we can use the
empirical variance of $\widehat{F}_{\textsc{F-IS}}$ or
$\widehat{F}_{\textsc{S-IS}}$:

######  Lemma 4.4 (Empirical Bernstein Bound).

Let $\mathbb{V}_{n}$ denote the sample variance. For any choice of $M$ points
$\\{t_{j}\\}_{j=1}^{M}$ where $t_{j}\in[0,D]$, with probability at least
$1-\delta$ we have that $\forall j\in[M]$,

$\displaystyle|F(t_{j})-\widehat{F}_{\textsc{F-IS}}(t_{j})|$
$\displaystyle\leq\frac{c_{1}^{\prime}\log(2M/\delta)}{n}+\sqrt{\frac{c_{2}^{\prime}\mathbb{V}_{n}[\widehat{F}_{\textsc{IS}}(t_{j})]\log(2M/\delta)}{n}},$
$\displaystyle|F(t_{j})-\widehat{F}_{\textsc{S-IS}}(t_{j})|$
$\displaystyle\leq\frac{c_{1}^{\prime}\log(2M/\delta)}{n}+\sqrt{\frac{c\mathbb{V}_{n}[\widehat{S}_{\textsc{IS}}(t_{j})]\log(2M/\delta)}{n}}.$

Note that Lemma 4.4 is a bound for a finite number of points, while Lemma 4.3
is a uniform bound over all $t$. In practice, the sample variance of only a
finite number of points can be calculated.

###  4.4 Combining IS Estimators

Lemma 4.3 also indicates that the variance of the resulting estimator can be
further reduced by combining the $F$ and $S$ estimators pointwise, such that
the estimator with the lower variance is used at each $t$. Such a combined
estimator may be especially effective for reducing the variance of the plug-in
estimators for distortion risk functionals that place weight on both upper and
lower tails of the distribution, such as the expected return.

How should the choice between $\widehat{F}_{\textsc{F-IS}}(t)$ or
$\widehat{F}_{\textsc{S-IS}}(t)$ be made for each $t$? Lemma 4.4 suggests that
one effective method is to estimate the empirical variance of
$\widehat{F}_{\textsc{F-IS}}(t)$ and $\widehat{F}_{\textsc{S-IS}}(t)$ at each
$t$, and to simply use the estimate with lower empirical variance:

$\widehat{F}_{\textsc{C-IS}}(t)=\begin{cases}\widehat{F}_{\textsc{F-IS}}(t),&\;\text{if}\;\mathbb{V}_{n}[\widehat{F}_{\textsc{F-IS}}(t)]<\mathbb{V}_{n}[\widehat{F}_{\textsc{S-IS}}(t)],\\\
\widehat{F}_{\textsc{S-IS}}(t),&\;\text{otherwise.}\end{cases}$ (9)

The error bound of such a combined estimator is given below.

######  Proposition 4.1.

Let
$\mathbb{V}_{n}^{min}(t)=\min\\{\mathbb{V}_{n}[\widehat{F}_{\textsc{F-IS}}(t)],\mathbb{V}_{n}[\widehat{F}_{\textsc{S-IS}}(t)]\\}$.
For any choice of $M$ points $\\{t_{j}\\}_{j=1}^{M}$ for $t_{j}\in[0,D]$ and
$\delta\in(0,1)$, define

$\Delta(\mathbb{V}_{n}):=\sqrt{\frac{2\log(2M/\delta)}{n-1}}+\max_{j\in[M]}\Big{|}\mathbb{V}_{n}[\widehat{F}_{\textsc{F-IS}}](t_{j})-\mathbb{V}_{n}[\widehat{F}_{\textsc{S-IS}}](t_{j})\Big{|}.$

Then, for $n\geq\frac{8\log(2n/\delta)}{\Delta(\mathbb{V}_{n})^{2}}$, we have
with probability at least $1-\delta$ that $\forall j\in[M]$,

$\displaystyle|F(t_{j})-\widehat{F}_{\textsc{C-IS}}(t_{j})|$
$\displaystyle\leq\frac{c_{1}^{\prime\prime}\log(2M/\delta)}{n}+\sqrt{\frac{c_{2}^{\prime\prime}\mathbb{V}_{n}^{min}(t_{j})\log(2M/\delta)}{n}}.$

##  5 Model-Based CDF Estimation

In the expected value OPE literature, model-based estimation can avoid the
high variance associated with importance sampling by directly learning the
expected return under $\pi$ using a model of the underlying MDP learned from
data (Jiang and Li, 2016). Taking inspiration from such approaches, we now
develop the first model-based off-policy CDF estimator. First, a model of the
MDP
$\overline{\mathcal{M}}=(\mathcal{S},\mathcal{A},\overline{P},\overline{R},\gamma,\overline{\mu})$
is formed from data, upon which a model $\overline{F}_{s_{h},a_{h}}$ of the
return distribution under $\pi$ can be directly computed using the recursion
in (1) for all horizons $h$, states $s_{h}$, and $a_{h}$. Formally, for
horizons $h=H,\ldots,1$, $\overline{F}$ can be computed recursively:

$\displaystyle\overline{F}_{s_{h},a_{h}}(t)=\mathbb{E}_{\overline{P},\overline{R}}\left[\overline{F}_{S_{h+1}}\left(\frac{t-R_{h}}{\gamma}\right)\right]$
$\displaystyle\overline{F}_{s_{h}}(t)=\mathbb{E}_{\pi}\left[\overline{F}_{s_{h},A_{h}}(t)\right]$

with $\overline{F}^{0}_{s_{H+1}}(t)=\mathbbm{1}\\{0\leq t\\}$ for all
$s_{H+1}$ at the end of the horizon. The models $\overline{F}$ can then be
used for the _direct method_ (DM) of estimation, which simply averages
$\overline{F}$ for the starting states observed in the data:

$\widehat{F}_{\textsc{DM}}(t)=\frac{1}{n}\sum_{i=1}^{n}\overline{F}_{s_{1}^{i}}(t)$
(10)

Because this direct method of estimation directly deploys the target policy
$\pi$ in the estimated model, the error of the estimators
$\widehat{F}_{\textsc{DM}}$ and $\overline{F}$ are directly related to the
misspecification of the MDP model $\overline{\mathcal{M}}$. The model
$\overline{\mathcal{M}}$ and thus the direct estimator is frequently prone to
uncharacterizable and possibly high bias, e.g., when the dataset coverage is
inadequate or the state space is high dimensional.

##  6 Doubly Robust CDF Estimation

We now define a new _doubly robust_ (DR) estimator for the distribution of
returns in MDPs (§ 6.1), which takes advantage of both importance sampling (§
4) and model-based (§ 5) estimators to retain the unbiasedness of the IS CDF
estimator, with potentially significant reduction in variance. In § 6.3, we
derive the first Cramer-Rao lower bound on off-policy CDF estimation, and show
that the DR estimator achieves this lower bound when the model $\overline{F}$
is equal to the true distribution.

###  6.1 Doubly Robust (DR) Estimator

The recursive formulation in (3) is key for defining the doubly robust CDF
estimator. At each time step $h$ in the recursive formulation, we use the
model $\overline{F}$ as a baseline and apply an IS-weighted data-dependent
correction to obtain the doubly robust estimator. Specifically, the DR
estimator has the following recursive form:

$\displaystyle\quad\quad\quad\widehat{F}_{\textsc{DR}}(t)=\frac{1}{n}\sum_{i=1}^{n}\widehat{F}_{s_{1}^{i}}(t),$
where
$\displaystyle\widehat{F}_{s_{h}^{i}}(t)=\overline{F}_{s_{h}^{i}}(t)+w(a_{h}^{i},s_{h}^{i})\left(\widehat{F}_{s_{h+1}^{i}}\left(\frac{t-r_{h}^{i}}{\gamma}\right)-\overline{F}_{s_{h}^{i},a_{h}^{i}}(t)\right)$
(11)

The model $\overline{F}$ can be obtained using the methods described in § 5,
and as before, $\widehat{F}^{0}_{s_{H+1}}(t)=\mathbbm{1}\left\\{0\leq
t\right\\}$ for all $s_{H+1}$. Note that when $H=1$, (11) recovers the DR
estimator for contextual bandits previously derived in Huang et al. (2021b).
Letting $z_{h}^{i}=\sum_{k=1}^{h}\gamma^{k-1}r_{k}^{i}$ be the return of
trajectory $i$ up until step $h$, the recursion in (11) can be unrolled over
the $H$ timesteps of the trajectory to obtain

$\displaystyle\widehat{F}_{\textsc{DR}}(t)=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\bigg{(}w_{H}^{i}\mathbbm{1}\\{z_{H}^{i}\leq
t\\}+\sum_{h=1}^{H}w_{h-1}^{i}\overline{F}_{s^{i}_{h}}\left(\frac{t-z_{h}^{i}}{\gamma^{h}}\right)-\sum_{h=1}^{H}w_{h}^{i}\overline{F}_{s^{i}_{h},a^{i}_{h}}\left(\frac{t-z_{h}^{i}}{\gamma^{h}}\right)\Bigg{)}.$

The bias and variance of the DR estimator is derived below.

######  Lemma 6.1.

The doubly robust estimator is unbiased, and letting $\widehat{F}_{s_{h}}$ be
given as in (11), its variance is:

$\displaystyle\mathbb{V}_{h}^{\beta}\left[\widehat{F}_{S_{h}}(t)\right]=$
$\displaystyle\mathbb{E}_{h}^{\beta}\left[\mathbb{V}_{h}^{\beta}\left[w(A_{h},S_{h})\Delta_{S_{h},A_{h}}(t)\Big{|}S_{h}\right]\right]$
$\displaystyle+\mathbb{E}_{h}^{\beta}\left[w(A_{h},S_{h})^{2}\mathbb{V}_{h+1}\left[\widehat{F}_{S_{h+1}}\left(\frac{t-R_{h}}{\gamma}\right)\Big{|}S_{h},A_{h}\right]\right]$
$\displaystyle+\mathbb{V}_{h}\left[F_{S_{h}}(t)\right]$

where $\Delta_{s,a}(t)=\overline{F}_{s,a}(t)-F_{s,a}(t)$.

Comparing the variance expression above with that of the IS estimator (Lemma
4.1), we see that the first term in the variance of the DR estimator contains
a difference $\Delta_{s,a}=\overline{F}_{s,a}-F_{s,a}$, while the variance of
the IS estimator contains only $F$. The variance in the first term is over
stochastic actions given the current state.

When $\Delta_{s,a}<F_{s,a}$, which is often the case in practice, the variance
contributed by the first term can be significantly reduced in the DR
estimator. Consequently, the DR estimator can be seen to reduce variance
contributed by importance sampling for stochastic actions, but not variance
from rewards or transitions, which arise from the last two terms in Lemma 6.1.

To further reduce variance, taking inspiration from Thomas et al. (2015), we
can combine weighted importance sampling with the DR estimator (11). Letting
$\widetilde{w}_{h}=\sum_{i=1}^{n}w_{h}^{i}$ and
$z_{h}^{i}=\sum_{k=1}^{h}\gamma^{k-1}r^{i}_{k}$,

$\displaystyle\widehat{F}_{\textsc{WDR}}(t)=$
$\displaystyle\sum_{i=1}^{n}\Bigg{(}\frac{w^{i}_{H}}{\widetilde{w}_{H}}\mathbbm{1}\left\\{z_{H}^{i}\leq
t\right\\}+\sum_{h=1}^{H}\frac{w^{i}_{h-1}}{\widetilde{w}_{h-1}}\overline{F}_{s^{i}_{h}}\left(\frac{t-z_{h}^{i}}{\gamma^{h}}\right)-\sum_{h=1}^{H}\frac{w^{i}_{h-1}}{\widetilde{w}_{h-1}}w(a^{i}_{h},s^{i}_{h})\overline{F}_{s^{i}_{h},a^{i}_{h}}\left(\frac{t-z_{h}^{i}}{\gamma^{h}}\right)\Bigg{)}.$

Neither the DR nor WDR estimators are guaranteed to valid CDFs; they may be
outside the $[0,1]$ interval or lose monotonicity due to the subtracted terms
(particularly at lower sample sizes). In practice, we transform the DR (and
WDR) estimators to valid CDFs using:

$\widehat{F}_{\textsc{M-DR}}(t)=\text{Clip}\left(\max_{t^{\prime}\leq
t}\widehat{F}_{\textsc{DR}}(t^{\prime}),0,1\right),$ (12)

which simply does not allow the CDF to decrease, and clips the value to the
unit interval.

###  6.2 Error Bound and Consistency

The WDR estimator is difficult to analyze, and we empirically demonstrate its
efficacy in § 8. The advantage of the M-DR estimator in (12), however, is that
it obeys the following finite sample error bound.

######  Lemma 6.2.

With probability at least $1-\delta$,

$\displaystyle\left\|\widehat{F}_{\textsc{M-DR}}-F\right\|_{\infty}\leq
w_{max}^{H}\sqrt{\frac{72}{n}\log\frac{8n^{1/2}}{\delta}},$

and thus $\widehat{F}_{\textsc{M-DR}}$ is uniformly consistent.

We achieve a convergence rate of $O(w_{max}^{H}/\sqrt{n})$. However, we note
that the above bound is conservative in the sense that it assumes worst-case
model misspecification, under which the DR estimator can perform worse than
IS. This contributes looseness to the result, and a tighter bound remains an
open problem.

###  6.3 Cramer-Rao Lower Bound

Though we have demonstrated the variance reduction potential of the DR
estimator, a natural question to ask is whether there exists an unbiased
estimator of the CDF that can even further reduce variance. Inspired by Jiang
and Li (2016); Huang and Jiang (2020), we derive the first Cramer-Rao lower
bound on the variance of off-policy CDF estimation in MDPs, and show that the
DR estimator can achieve this lower bound.

######  Theorem 6.1.

For unique directed acyclic graphs (Definition C.1), the variance of any
unbiased off-policy CDF estimator is pointwise lower bounded by

$\sum_{h=1}^{H}\mathbb{E}_{\mathbb{P}_{\beta}}\left[w_{h-1}^{2}\mathbb{V}_{h}\left[F_{S_{h}}\left(t-\sum_{k=1}^{h}R_{k}\right)\right]\right].$

This lower bound on the variance shows that difficulties introduced by the
intrinsic state transition stochasticity of the underlying MDP cannot be
eliminated. In fact, the variance of the DR estimator in unique directed
acyclic graphs is equal to the Cramer Rao lower bound when $\overline{F}$ is a
perfect model, that is, $\overline{F}_{s_{h},a_{h}}(t)=F_{s_{h},a_{h}}(t)$ for
all $s_{h},a_{h}$ and $t$. Further, as the error between $\overline{F}$ and
$F$ decreases, the variance approaches the Cramer-Rao lower bound.

##  7 Risk Estimation

We have developed several methods for CDF estimation and their error bounds,
which we now translate to off-policy risk assessment. Returning to Algorithm
1, we have the CDF estimator $\widehat{F}$ and the CDF estimation error
$\epsilon$ (Line 2). We form the plug-in estimate of a risk functional $\rho$
as:

$\widehat{\rho}=\rho(\widehat{F})$ (13)

In this section, we characterize the theoretical properties of the risk
estimate $\widehat{\rho}$ and show that the properties of the CDF estimators
translate directly to the risk estimators built upon them. Notably, the
reduced variance of the DR CDF estimator results in risk estimators with
reduced bias and variance, and for Lipschitz risk functionals (Definition
3.1), the error of $\widehat{\rho}$ is proportional $\epsilon$. For uniformly
consistent CDF estimators, such as the previously introduced IS (4) and DR
(11) estimators, this implies that their downstream risk estimates will be
consistent as well.

###  7.1 Bias, Variance, and Error

In general, the risk estimator $\widehat{\rho}$ (13) can be biased even if the
CDF estimator is unbiased. Intuitively, however, $\widehat{\rho}$ may have
lower variance (and thus error) when it is estimated from a CDF estimate
$\widehat{F}$ with lower variance, because it is computed directly from
$\widehat{F}$ itself.

For Lipschitz risk functionals, we formalize this intuition as a finite-sample
convergence guarantee (Corollary 7.1), which demonstrates that the error of
any Lipschitz risk estimator is upper bounded by its Lipschitz constant (from
Definition 3.1), and the error (upper bound) of the CDF estimate it utilizes.

######  Corollary 7.1.

For all Lipschitz risk functionals simultaneously, given a CDF estimator
$\widehat{F}$ with error $\epsilon$, we have with probability at least
$1-\delta$ that

$|\rho(F)-\rho_{p}(\widehat{F})|\leq L\epsilon,$ (14)

where $L$ is the Lipschitz constant of $\rho$. Further, if $\widehat{F}$ is
uniformly consistent, then $\rho(\widehat{F})\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\rho(F)$.

As we have seen in Sections 4 and 6, the error guarantees for $\widehat{F}$
are functions of its variance. Thus, Corollary 7.1 demonstrates that we can
expect to have faster convergence for risk estimates estimated on lower-
variance CDFs. As an additional consequence, downstream risk estimators built
from the IS or DR CDF estimators, which have error that scale with
$O(1/\sqrt{n})$, are guaranteed to be consistent.

###  7.2 Lower Bound

We have previously demonstrated that the IS CDF estimator is minimax rate-
optimal, and demonstrated that the DRCDF estimator can achieve the variance
lower bound for off-policy CDF estimation. However, the natural question to
ask is: are their downstream risk estimators also rate-optimal? For general
risk functionals this remains an open problem, but we show that for the
Lipschitz risk functionals of $\text{CVaR}_{\alpha}$ and expected return
(which is equivalent to $\text{CVaR}_{0}$), their off-policy risk estimator
$\widehat{\rho}$ has the following lower bound.

######  Theorem 7.1.

Let $\mathcal{F}$ be the family of CDFs with bounded support and
$\rho=\text{CVaR}_{\alpha}$. For $\alpha\in[0,1)$ and $n\geq
c\frac{w_{max}^{2}}{\mathbb{E}_{\mathbb{P}_{\beta}}[w_{H}^{2}]}$ for a
universal constant $c>0$,

$\displaystyle\inf_{\widehat{\rho}}\sup_{F^{\pi}\in\mathcal{F}}\mathbb{E}\left[|\widehat{\rho}-\rho(F^{\pi})|\right]\geq\frac{cD}{1-\alpha}\sqrt{\frac{\mathbb{E}_{\mathbb{P}_{\beta}}[w_{H}^{2}]}{n}}.$

For large samples sizes, also matches existing results on off-policy mean
estimation (Li et al., 2015; Wang et al., 2017; Ma et al., 2021). As the
Lipschitz constants for $\text{CVaR}_{\alpha}$ and expected return are
$\frac{D}{1-\alpha}$ and $D$, respectively, the convergence rate of Corollary
7.1 with IS estimation is indeed minimax rate-optimal.

Figure 1: Sup-norm error of the CDF estimate (row 1), and mean squared error
of the expected return (row 2) and CVaR 0.25 (row 3) plug-in estimates for
different $\lambda$ (column) in Cliffwalk and Simglucose.

##  8 Experiments

We now provide empirical support for our theoretical results, and compare the
performance of the CDF and plug-in risk estimators in Figure 1. We compare the
DR and WDR estimators against the importance sampling estimators— F-IS (4),
S-IS (8), C-IS (9), WIS ($\S$ 4.2)—and model-based DM estimator. Our primary
baselines of comparison are the F-IS and WIS estimators, which are the only
two off-policy CDF/risk estimators in existing work (Chandak et al., 2021a).
As we will see, the DR and WDR estimators outperform the IS-based estimators
in all studied cases.

##### Environment

We utilize two experimental domains: a tabular Cliffwalk (Sutton and Barto,
2018), and continuous diabetes treatment simulator (Xie, 2018). The horizons
of both are $H=200$, and the behavioral policy $\beta$ is a mixture between
the target $\pi$ and a uniform policy UNIF, i.e.,
$\beta=\lambda\pi+(1-\lambda)\text{UNIF}$. Appendix E includes the full setup
and additional experiments. Both environments have risk-related
interpretations: in Cliffwalk the agent must take the shortest path while
avoiding a slippery, high-cost cliff, while in Simglucose the agent must
maintain a patient’s glucose levels within a healthy range.

##### Risk Functionals

As such, we estimate one each of risk-neutral, risk-averse, and risk-seeking
risk functionals using OPRA: the expected return, $\text{CVaR}_{\alpha}$, and
a risk functional we call $\text{CCaR}_{\alpha}$ (conditional cost-at-risk),
respectively. While $\text{CVaR}_{\alpha}$ is the expectation of the worst-
case $\alpha$ tail of the distribution, CCaR is the expectation of the best-
case $\alpha$ tail.

##### Results

In both environments, the DR and WDR estimators outperform the IS estimators
and their variants in CDF and risk estimation error (Figure 1). These
environments present challenges for IS estimators because of their long
horizons, but DR and WDR effectively overcome them by incorporating model
information. At very low sample sizes, however, the DM estimator can
outperform the other estimators.

Our results also demonstrate the variability of the CDF estimator performance
for different risk functionals; the gap between F-IS and S-IS is especially
obvious. Under the Cliffwalk cost distribution, S-IS far outperforms F-IS for
both the risk-neutral expected return and risk-sensitive
$\text{CVaR}_{\alpha}$, For the risk-seeking $\text{CCaR}_{\alpha}$, however,
F-IS outperforms S-IS. The combined estimator C-IS never does worse than
either F-IS or S-IS, and as $n$ increases can do better than both, as
Proposition 4.1 implied. The DM estimator has low error under
$\text{CCaR}_{\alpha}$ because the model estimates the lower tail well, but
the upper tail poorly. The DR and WDR estimator straddle best of both worlds;
they perform well under all evaluated risk functionals.

##  9 Discussion

In this work, we have introduced and analyzed variance-reduced off-policy CDF
and risk estimators for MDPs, including the first doubly robust estimator for
return CDF estimation. Several future research directions are of interest.
First, it is possible that estimators and bounds can be tailored and tightened
for individual risk functionals, and a thorough comparison of different risk
estimators is an important avenue of future work. Second, our results have
also highlighted the risk estimation performance disparity between different
types of CDF estimators. This disparity is more severe for risk-averse or
risk-seeking risk functionals than expected value. As such, adaptive estimator
selection (based on risk functionals of interest) may prove to be of both
theoretical and practical interest. Lastly, in off-policy evaluation of the
mean, marginalized importance sampling (Xie et al., 2019) has been shown to
achieve the Cramer-Rao lower bound, and model-free regression (Duan et al.,
2020) has been proven to be minimax optimal. One important future direction is
to investigate whether these methods can be extended to off-policy CDF or risk
estimation, which, unlike the mean, are nonlinear functions of the return.

## References

  * Artzner et al. [1999] Philippe Artzner, Freddy Delbaen, Jean-Marc Eber, and David Heath.  Coherent measures of risk.  _Mathematical finance_ , 9(3):203–228, 1999.
  * Balbás et al. [2009] Alejandro Balbás, José Garrido, and Silvia Mayoral.  Properties of distortion risk measures.  _Methodology and Computing in Applied Probability_ , 11(3):385–399, 2009.
  * Bellemare et al. [2017] Marc G Bellemare, Will Dabney, and Rémi Munos.  A distributional perspective on reinforcement learning.  In _International Conference on Machine Learning_ , pages 449–458. PMLR, 2017.
  * Chandak et al. [2021a] Yash Chandak, Scott Niekum, Bruno Castro da Silva, Erik Learned-Miller, Emma Brunskill, and Philip S Thomas.  Universal off-policy evaluation.  _arXiv preprint arXiv:2104.12820_ , 2021a.
  * Chandak et al. [2021b] Yash Chandak, Shiv Shankar, and Philip S. Thomas.  High-confidence off-policy (or counterfactual) variance estimation, 2021b.
  * Chang et al. [2020] Joel QL Chang, Qiuyu Zhu, and Vincent YF Tan.  Risk-constrained thompson sampling for cvar bandits.  _arXiv preprint arXiv:2011.08046_ , 2020.
  * Dabney et al. [2018] Will Dabney, Georg Ostrovski, David Silver, and Rémi Munos.  Implicit quantile networks for distributional reinforcement learning.  In _International conference on machine learning_ , pages 1096–1105. PMLR, 2018.
  * Delbaen [2002] Freddy Delbaen.  Coherent risk measures on general probability spaces.  In _Advances in finance and stochastics_ , pages 1–37. Springer, 2002\.
  * Denardo et al. [2007] Eric V Denardo, Haechurl Park, and Uriel G Rothblum.  Risk-sensitive and risk-neutral multiarmed bandits.  _Mathematics of Operations Research_ , 32(2):374–394, 2007.
  * Denneberg [1990] Dieter Denneberg.  Distorted probabilities and insurance premiums.  _Methods of Operations Research_ , 63(3):3–5, 1990.
  * Duan et al. [2020] Yaqi Duan, Zeyu Jia, and Mengdi Wang.  Minimax-optimal off-policy evaluation with linear function approximation.  In _International Conference on Machine Learning_ , pages 2701–2709. PMLR, 2020.
  * Dudik et al. [2011] Miroslav Dudik, John Langford, and Lihong Li.  Doubly robust policy evaluation and learning, 2011.
  * Dudík et al. [2014] Miroslav Dudík, Dumitru Erhan, John Langford, Lihong Li, et al.  Doubly robust policy evaluation and optimization.  _Statistical Science_ , 29(4):485–511, 2014.
  * Gopalan et al. [2017] Aditya Gopalan, LA Prashanth, Michael Fu, and Steve Marcus.  Weighted bandits or: How bandits learn distorted values that are not expected.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 31, 2017.
  * Huang et al. [2021a] Audrey Huang, Liu Leqi, Zachary C. Lipton, and Kamyar Azizzadenesheli.  On the convergence and optimality of policy gradient for markov coherent risk, 2021a.
  * Huang et al. [2021b] Audrey Huang, Liu Leqi, Zachary C Lipton, and Kamyar Azizzadenesheli.  Off-policy risk assessment in contextual bandits.  _arXiv preprint arXiv:2104.08977_ , 2021b.
  * Huang and Jiang [2020] Jiawei Huang and Nan Jiang.  From importance sampling to doubly robust policy gradient.  In _International Conference on Machine Learning_ , pages 4434–4443. PMLR, 2020.
  * Jiang and Li [2016] Nan Jiang and Lihong Li.  Doubly robust off-policy value evaluation for reinforcement learning.  In _International Conference on Machine Learning_ , pages 652–661. PMLR, 2016.
  * Keramati et al. [2020] Ramtin Keramati, Christoph Dann, Alex Tamkin, and Emma Brunskill.  Being optimistic to be conservative: Quickly learning a cvar policy.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pages 4436–4443, 2020.
  * Kusuoka [2001] Shigeo Kusuoka.  On law invariant coherent risk measures.  In _Advances in mathematical economics_ , pages 83–95. Springer, 2001\.
  * Li et al. [2015] Lihong Li, Rémi Munos, and Csaba Szepesvári.  Toward minimax off-policy value estimation.  In _Artificial Intelligence and Statistics_ , pages 608–616. PMLR, 2015.
  * Liu et al. [2018] Qiang Liu, Lihong Li, Ziyang Tang, and Dengyong Zhou.  Breaking the curse of horizon: Infinite-horizon off-policy estimation.  _arXiv preprint arXiv:1810.12429_ , 2018.
  * Ma et al. [2021] Cong Ma, Banghua Zhu, Jiantao Jiao, and Martin J Wainwright.  Minimax off-policy evaluation for multi-armed bandits.  _arXiv preprint arXiv:2101.07781_ , 2021.
  * Mannor and Tsitsiklis [2011] Shie Mannor and John Tsitsiklis.  Mean-variance optimization in markov decision processes.  _arXiv preprint arXiv:1104.5601_ , 2011.
  * Maurer and Pontil [2009] Andreas Maurer and Massimiliano Pontil.  Empirical bernstein bounds and sample variance penalization.  _arXiv preprint arXiv:0907.3740_ , 2009.
  * Prashanth et al. [2016] LA Prashanth, Cheng Jie, Michael Fu, Steve Marcus, and Csaba Szepesvári.  Cumulative prospect theory meets reinforcement learning: Prediction and control.  In _International Conference on Machine Learning_ , pages 1406–1415. PMLR, 2016.
  * Prashanth et al. [2020] LA Prashanth, Krishna Jagannathan, and Ravi Kolla.  Concentration bounds for cvar estimation: The cases of light-tailed and heavy-tailed distributions.  In _International Conference on Machine Learning_ , pages 5577–5586. PMLR, 2020.
  * Rockafellar et al. [2000] R Tyrrell Rockafellar, Stanislav Uryasev, et al.  Optimization of conditional value-at-risk.  _Journal of risk_ , 2:21–42, 2000.
  * Sani et al. [2013] Amir Sani, Alessandro Lazaric, and Rémi Munos.  Risk-aversion in multi-armed bandits.  _arXiv preprint arXiv:1301.1936_ , 2013.
  * Shapiro et al. [2014] Alexander Shapiro, Darinka Dentcheva, and Andrzej Ruszczyński.  _Lectures on stochastic programming: modeling and theory_.  SIAM, 2014.
  * Sobel [1982] Matthew J Sobel.  The variance of discounted markov decision processes.  _Journal of Applied Probability_ , pages 794–802, 1982.
  * Su et al. [2020] Yi Su, Maria Dimakopoulou, Akshay Krishnamurthy, and Miroslav Dudík.  Doubly robust off-policy evaluation with shrinkage.  In _International Conference on Machine Learning_ , pages 9167–9176. PMLR, 2020.
  * Sussex et al. [2018] Scott Sussex, Omer Gottesman, Yao Liu, Susan Murphy, Emma Brunskill, and Finale Doshi-Velez.  Stitched trajectories for off-policy learning.  2018\.
  * Sutton and Barto [2018] Richard S Sutton and Andrew G Barto.  _Reinforcement learning: An introduction_.  MIT press, 2018.
  * Tamar et al. [2015] Aviv Tamar, Yinlam Chow, Mohammad Ghavamzadeh, and Shie Mannor.  Policy gradient for coherent risk measures.  In _Proceedings of the 28th International Conference on Neural Information Processing Systems-Volume 1_ , pages 1468–1476, 2015.
  * Tamar et al. [2016] Aviv Tamar, Dotan Di Castro, and Shie Mannor.  Learning the variance of the reward-to-go.  _The Journal of Machine Learning Research_ , 17(1):361–396, 2016.
  * Tamkin et al. [2019] Alex Tamkin, Ramtin Keramati, Christoph Dann, and Emma Brunskill.  Distributionally-aware exploration for cvar bandits.  In _NeurIPS 2019 Workshop on Safety and Robustness on Decision Making_ , 2019.
  * Thomas and Brunskill [2016] Philip Thomas and Emma Brunskill.  Data-efficient off-policy policy evaluation for reinforcement learning.  In _International Conference on Machine Learning_ , pages 2139–2148. PMLR, 2016.
  * Thomas et al. [2015] Philip Thomas, Georgios Theocharous, and Mohammad Ghavamzadeh.  High-confidence off-policy evaluation.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 29, 2015.
  * Wang [1996] Shaun Wang.  Premium calculation by transforming the layer premium density.  _ASTIN Bulletin: The Journal of the IAA_ , 26(1):71–92, 1996.
  * Wang et al. [1997] Shaun S Wang, Virginia R Young, and Harry H Panjer.  Axiomatic characterization of insurance prices.  _Insurance: Mathematics and economics_ , 21(2):173–183, 1997.
  * Wang et al. [2017] Yu-Xiang Wang, Alekh Agarwal, and Miroslav Dudık.  Optimal and adaptive off-policy evaluation in contextual bandits.  In _International Conference on Machine Learning_ , pages 3589–3597. PMLR, 2017.
  * Wirch and Hardy [2001] Julia L Wirch and Mary R Hardy.  Distortion risk measures: Coherence and stochastic dominance.  In _International congress on insurance: Mathematics and economics_ , pages 15–17, 2001.
  * Xie [2018] Jinyu Xie.  Simglucose v0. 2.1 (2018).  _URL https://github. com/jxx123/simglucose_ , 2018.
  * Xie et al. [2019] Tengyang Xie, Yifei Ma, and Yu-Xiang Wang.  Towards optimal off-policy evaluation for reinforcement learning with marginalized importance sampling.  _Advances in Neural Information Processing Systems_ , 32, 2019.

##  Appendix A Preliminaries: Risk Overview

Let $Z$ be a random variable. A _risk functional_ $\rho$ is a mapping from a
space of random variables to the space of real numbers
$\rho:\mathcal{L}_{\infty}(\Omega,\mathcal{F}_{Z},\mathbb{P}_{Z})\rightarrow\mathbb{R}$.
Risk functionals are law-invariant if they depend only on the distribution of
$Z$ [Kusuoka, 2001]. Formally $\rho(Z)$ is law invariant if for any pair of
random variables $Z$ and $Z^{\prime}$,
$F_{Z}=F_{Z^{\prime}}\;\Longrightarrow\;\rho(Z)=\rho(Z^{\prime}).$

Lipschitz risk functionals are a subset of law-invariatn risk functionals, and
we now delineate popular classes of risks and their associated Lipschitz
constants, where possible, adapted from Huang et al. [2021a].

##### Distorted Risk Functionals

For $Z\geq 0$, distortion risk functionals are defined to be Denneberg [1990],
Wang [1996], Wang et al. [1997], Balbás et al. [2009]

$\displaystyle\rho(F_{Z})=\int_{0}^{\infty}g(1-F_{Z}(t))dt=\int_{0}^{\infty}g(S_{Z}(t))dt$

where the distortion $g:[0,1]\rightarrow[0,1]$ has $g(0)=0$ and $g(1)=1$, and
is increasing. When $g(x)=x$, the expected value is recovered. When
$g(x)=\min\\{\frac{x}{1-\alpha},1\\}$ for $\alpha\in(0,1)$, CVaR at level
$\alpha$ is recovered.

When $g(x)=F(F^{-1}(x)-F^{-1}(\alpha))$, the Wang risk functional at level
$\alpha$ [Wang, 1996] is recovered, and the proportional hazard risk
functional can by obtained by setting $g(x)=x^{\alpha}$ for $\alpha<1$.
Distorted risk functionals are coherent if and only if $g$ is concave [Wirch
and Hardy, 2001]. Not all distorted risk functionals are coherent. For
example, setting $g(s)=\mathbbm{1}_{\\{s\geq 1-\alpha\\}}$ recovers the value-
at-risk (VaR), which is not coherent.

##### Coherent Risk Functionals

Coherent risk functionals satisfy properties called monotonicity,
subadditivity, translation invariance, and positive homogeneity Artzner et al.
[1999], Delbaen [2002]. Examples include expected value, conditional value-at-
risk (CVaR), entropic value-at-risk, and mean semideviation, proportional
hazard, and Wang transform [Chang et al., 2020, Tamkin et al., 2019, Tamar et
al., 2015, Shapiro et al., 2014, Wirch and Hardy, 2001, Wang, 1996].

######  Lemma A.1 (Lipschitzness of Coherent and Distorted Risk Functionals
Huang et al. [2021b]).

On the space of random variables with support in $[0,D]$, the distorted risk
functional of any $\frac{L}{D}$-Lipschitz distortion function is a
$L$-Lipschitz risk functional.

The expected value risk functional is $D$-Lipschitz and $\text{CVaR}_{\alpha}$
is $\frac{D}{\alpha}$-Lipschitz.

##### Cumulative Prospect Theory (CPT) Risk Functionals

CPT risks Prashanth et al. [2016] are defined as:

$\displaystyle\rho(F_{Z})=\int_{0}^{+\infty}g^{+}\left(1-F_{u^{+}(Z)}(t)\right)dt-\int_{0}^{+\infty}g^{-}\left(1-F_{u^{-}(Z)}(t)\right)dt,$

where $g^{+},g^{-}$: $[0,1]\rightarrow[0,1]$ are continuous with
$g^{+/-}(0)=0$ and $g^{+/-}(1)=1$. The functions
$u^{+},u^{-}:\mathbb{R}\rightarrow\mathbb{R}_{+}$ are continuous, with
$u^{+}(z)=0$ when $z\geq c$ and $u^{-}(z)=0$ when $z<c$ for some constant
$c\in\mathbb{R}$. The CPT functional handles gains and losses separately,
applying distortion $g^{+}$ to the part of the distribution corresponding to
gains (controlled by $u^{+}$), and $g^{-}$ to the part corresponding to losses
(controlled by $u^{-}$).

######  Lemma A.2 (Lipschitzness of CPT Functional).

If the CPT distortion functions $g^{+}$ and $g^{-}$ are both
$\frac{L}{D}$-Lipschitz, then the CPT risk functional is $L$-Lipschitz.

##### Other Risk Functionals

Many other risk functionals, though not members of a specific class, can be
shown to be Lipschitz. This includes the variance and mean-variance. The
variance is defined as
$\rho(F_{Z})=2\int_{0}^{\infty}t(1-F_{Z}(t))dt-\left(\int_{0}^{\infty}(1-F_{Z}(t))dt\right)^{2}.$

######  Lemma A.3 (Lipschitzness of Variance).

On the space of random variables with support in $[0,D]$, variance is a
$3D^{2}$-Lipschitz risk functional. For mean-variance given by
$\rho(Z)=\mathbb{E}[Z]+\lambda\mathbb{V}(Z)$ for some $\lambda>0$, it is
$(1+3\lambda D^{2})$-Lipschitz.

##  Appendix B Proofs for Importance Sampling CDF Estimation (Section 4)

###  B.1 Bias and Variance (Proof of Lemma 4.1)

As the IS estimator can be seen as a special case of the DR estimator with
$\overline{F}=0$, we defer the proof of bias and variance to Appendix C.1.

###  B.2 Error Upper Bound (Proof of Lemma 4.2)

The proof of the upper bound closely follows the proof of Theorem 5.1 from
Huang et al. [2021b].

For any $\lambda>0$, we have

$\displaystyle\lambda\mathbb{E}\left[\|\widehat{F}-F\|_{\infty}\right]=$
$\displaystyle\mathbb{E}\left[\log\exp\left(\lambda\|\widehat{F}-F\|_{\infty}\right)\right]$
$\displaystyle\leq\log\mathbb{E}\left[\exp\left(\lambda\|\widehat{F}-F\|_{\infty}\right)\right]$

To upper bound the RHS, define the following function class:

$\displaystyle\mathbb{F}(n):=\left\\{f(g):=\varrho\frac{1}{n}\mathbbm{1}_{\left\\{g\leq
t\right\\}}:\forall t\in\mathbb{R};\forall
g\in\mathbb{Q},\varrho\in\left\\{-1,+1\right\\}\right\\}$

Note that this is a countable set. Denote
$w(\tau)=\prod_{h=1}^{H}w(a_{h},s_{h})$, and
$g^{i}=\sum_{h=1}^{H}\gamma^{h}r_{h}^{i}$ the return of trajectory $i$ for
short. Using this definition, we have

$\displaystyle\sup_{t\in\mathbb{R}}\left|\widehat{F}_{\textsc{IS}}(t)-F(t)\right|=\sup_{f\in\mathbb{F}(n)}\left|\frac{1}{n}\sum_{i=1}^{n}w(\tau^{i})f(g^{i})-\mathbb{E}_{\mathbb{P}_{\beta}}\left[\frac{1}{n}\sum_{i=1}^{n}w(\tau^{i})f(g^{i})\right]\right|$

Then adapting the proof of Theorem 5.1 in Huang et al. [2021b], we have for
any $\lambda\in(0,\frac{n}{2w_{max}})$ that,

$\displaystyle\mathbb{E}\left[\exp\left(\lambda\|\widehat{F}-F\|_{\infty}\right)\right]$
$\displaystyle\leq
4\exp\left(\frac{n\lambda^{2}\frac{4\mathbb{E}_{\mathbb{P}_{\beta}}[w(\tau)^{2}])}{n^{2}}}{2\left(1-\lambda\frac{2}{n}w_{max}\right)}\right)$

Then

$\displaystyle\mathbb{E}\left[\|\widehat{F}-F\|_{\infty}\right]$
$\displaystyle\leq\frac{1}{\lambda}\log\left(4\exp\left(\frac{n\lambda^{2}\frac{4\mathbb{E}_{\mathbb{P}_{\beta}}[w(\tau)^{2}])}{n^{2}}}{2\left(1-\lambda\frac{2}{n}w_{max}\right)}\right)\right)$
$\displaystyle\leq\frac{\log
4}{\lambda}+\frac{2\lambda\mathbb{E}_{\mathbb{P}_{\beta}}[w(\tau)^{2}]}{n(1-\lambda\frac{2}{n}w_{max})}$
(15)

Since this holds for any $\lambda\in(0,\frac{n}{2w_{max}})$, it also holds for
$\min_{\lambda}$ of the RHS. To find the minimizer, we solve

$\displaystyle\frac{d}{d\lambda}\left(\frac{\log
4}{\lambda}+\frac{2\lambda\mathbb{E}_{\mathbb{P}_{\beta}}[w(\tau)^{2}]}{n-2\lambda
w_{max}}\right)=0$

which WolframAlpha tells us is solved by

$\displaystyle\lambda^{*}=\frac{n}{\sqrt{\mathbb{E}_{\mathbb{P}_{\beta}}[w(\tau)^{2}]n\log
2}+2w_{max}}$

if $\mathbb{E}_{\mathbb{P}_{\beta}}[w(\tau)^{2}]n\neq w_{max}^{2}\log 16$. We
observe that $\lambda^{*}\in(0,\frac{n}{2w_{max}})$ satisfies the constraints
on $\lambda$. Plugging $\lambda^{*}$ back into (15) we have,

$\displaystyle\mathbb{E}\left[\|\widehat{F}-F\|_{\infty}\right]$
$\displaystyle\leq\frac{\log
4}{\frac{n}{\sqrt{\mathbb{E}_{\mathbb{P}_{\beta}}[w(\tau)^{2}]n\log
2}+2w_{max}}}+\frac{\frac{2n\mathbb{E}_{\mathbb{P}_{\beta}}[w(\tau)^{2}]}{\sqrt{\mathbb{E}_{\mathbb{P}_{\beta}}[w(\tau)^{2}]n\log
2}+2w_{max}}}{n-\frac{2w_{max}n}{\sqrt{\mathbb{E}_{\mathbb{P}_{\beta}}[w(\tau)^{2}]n\log
2}+2w_{max}}}$ $\displaystyle=\log 4\sqrt{\log
2}\sqrt{\frac{\mathbb{E}_{\mathbb{P}_{\beta}}[w(\tau)^{2}]}{n}}+\frac{2\log
4w_{max}}{n}+2\sqrt{\frac{\mathbb{E}_{\mathbb{P}_{\beta}}[w(\tau)^{2}]}{n\log
2}}$ $\displaystyle\leq
c_{1}\sqrt{\frac{\mathbb{E}_{\mathbb{P}_{\beta}}[w(\tau)^{2}]}{n}}+c_{2}\frac{w_{max}}{n}$

###  B.3 CDF Estimation Lower Bound (Proof of Theorem 4.1)

We begin by deriving a CDF estimation lower bound in multi-armed bandits
(MABs). We can reduce the estimation problem in the MDP case to the MAB case
by treating entire trajectories $\tau$ as “actions".

Let a multi-armed bandit be defined by $K$ arms, and suppose the behavioral
policy is given by $\beta(a)$ for $a\in[K]$, and the target policy is $\pi(a)$
for $a\in[K]$, with importance weight $w(a):=\beta(a)/\pi(a)$. Let
$w_{max}=\max_{a}w(a)$.

######  Lemma B.1 (CDF Lower Bound in MABs).

When $n\geq c\frac{w_{max}^{2}}{\sum_{a}\beta(a)w^{2}(a)}$ for a universal
constant $c$, the minimax lower bound on off-policy CDF estimation in MABs is

$\inf_{\widehat{F}}\sup_{F\in\mathcal{F}}\mathbb{E}\left[\|\widehat{F}-F^{\pi}\|_{\infty}\right]\geq
c\sqrt{\frac{\sum_{a\in[K]}\beta(a)w^{2}(a)}{n}}.$ (16)

We now lift this result to the MDP case. For a horizon of length $H$, the
space of possible actions of size
$|\tau|=|\mathcal{S}|^{H+1}|\mathcal{A}|^{H}$. Recall that in MDPs,
$w_{max}=\max_{\tau}w(\tau)$, and
$w(\tau)=\frac{\pi(\tau)}{\beta(\tau)}=\prod_{h=1}^{H}\frac{\pi(a_{h}|s_{h})}{\beta(a_{h}|s_{h})}$.
Thus using Lemma B.1, when $n\geq
c\frac{w_{max}^{2}}{\sum_{\tau}\mathbb{P}_{\beta}(\tau)w^{2}(\tau)},$ we have

$\displaystyle\inf_{\widehat{F}}\sup_{F\in\mathcal{F}}\mathbb{E}\left[\|\widehat{F}-F^{\pi}\|_{\infty}\right]\geq
c\sqrt{\frac{\sum_{\tau}\mathbb{P}_{\beta}(\tau)w^{2}(\tau)}{n}}.$ (17)

This matches our Bernstein-style bound on the error of IS estimation up to a
constant:

$\displaystyle\mathbb{E}\left[\|\widehat{F}-F\|_{\infty}\right]\leq
c_{1}\sqrt{\frac{\sum_{\tau}\mathbb{P}_{\beta}(\tau)w^{2}(\tau)}{n}}+c_{2}\frac{w_{max}}{n}.$
(18)

###### Proof of Lemma B.1.

We use Le Cam’s method to derive the lower bound. We first need to select two
problem MAB problem instances, by defining their reward distributions $R_{1}$,
$R_{2}$. Then define distributions $P_{1}=\beta\circ R_{1}$, and
$P_{2}=\beta\circ R_{2}$.

Intuitively, we want to select them so that the difference
$\|F_{1}^{\pi}-F_{2}^{\pi}\|_{\infty}$ is large enough that an estimator can
distinguish between them, but with the distance
$\text{D}_{\text{KL}}(P_{1}\|P_{2})$ to be small enough that the testing
problem is not too easy.

For each action $a\in[K]$, let $R_{2}(\cdot|a)$ be a Bernoulli distribution on
the set $\\{0,1\\}$ with probability $p=1/2$, and let $R_{1}(\cdot|a)$ be a
Bernoulli distribution with $p=1/2+\delta(a)$ for some $\delta(a)\in[0,1/2]$.
Then

$\displaystyle\|F_{1}^{\pi}-F_{2}^{\pi}\|_{\infty}$
$\displaystyle=\max_{t\in\left\\{0,1\right\\}}\left|F_{1}^{\pi}(t)-F_{2}^{\pi}(t)\right|$
$\displaystyle=\left|F_{1}^{\pi}(t=0)-F_{2}^{\pi}(t=0)\right|$
$\displaystyle=\left|\sum_{a}\pi(a)\delta(a)\right|$

At the same time,

$\displaystyle\text{D}_{\text{KL}}(P_{1}\|P_{2})$
$\displaystyle=\sum_{a}\beta(a)\text{D}_{\text{KL}}(R_{1}\|R_{2})\leq
4\sum_{a}\beta(a)\delta^{2}(a)$

Since we want the probability of error in our test to be a constant, we want
to upper bound the above by $\frac{1}{2n}$. Thus to obtain as tight a lower
bound as possible, we want to solve the optimization problem

$\max_{\delta\in\mathbb{R}^{K}}\sum_{a}\pi(a)\delta(a)\quad\text{s.t.}\quad
4\sum_{a}\beta(a)\delta^{2}(a)\leq\frac{1}{2n},\;\delta\in[0,1/2]$ (19)

Using Lemma 3 from Ma et al. [2021], the Fenchel dual of the optimization
problem in the above panel is

$\min_{\nu\in\mathbb{R}^{K}}\sqrt{\frac{1}{8n}\sum_{a}\frac{(\pi(a)-\nu(a))^{2}}{\beta(a)}}+\frac{1}{2}\sum_{a}|\nu(a)|$
(20)

Let $\nu^{*}$ be the solution to (20) and define the set of actions

$S^{*}=\\{a\in[K]|\nu^{*}(a)>0\\}$

Then Lemma 2 from Ma et al. [2021] states that

$\displaystyle\min_{\nu\in\mathbb{R}^{K}}\sqrt{\frac{1}{8n}\sum_{a}\frac{(\pi(a)-\nu(a))^{2}}{\beta(a)}}+\frac{1}{2}\sum_{a}|\nu(a)|\asymp\pi(S^{*})+\sqrt{\frac{\sum_{a\notin
S^{*}}\beta(a)w^{2}(a)}{n}}$

Thus we have that for some universal constant $c$,

$\inf_{\widehat{F}}\sup_{F\in\mathcal{F}}\mathbb{E}\left[\|\widehat{F}-F^{\pi}\|_{\infty}\right]\geq
c\left\\{\pi(S^{*})+\sqrt{\frac{\sum_{a\notin
S^{*}}\beta(a)w^{2}(a)}{n}}\right\\}$ (21)

When we have sample size $n$ such that $n\geq
c\frac{w_{max}^{2}}{\sum_{a}\beta(a)w^{2}(a)}$, $S^{*}=\emptyset$ and the
lemma statement follows. ∎

###  B.4 Proofs for Section 4.3

####  B.4.1 Bernstein Bounds for $S$ and $F$ (Proof of Lemma 4.3)

######  Definition B.1 (Bounded Variation).

A function $f$ has bounded variation if there exists $C<\infty$ such that

$\sup_{N,\\{t_{i}\\}_{i=1}^{N}}\sum_{i=1}^{N}|f(t_{i+1})-f(t_{i})|\leq C$

######  Assumption B.1.

Suppose that the random variable of returns takes values on a compact set.
Assume the variance of the IS estimator $\mathbb{V}^{\widehat{F}}(t)$ has
bounded variation for all $t\in\mathbb{R}$.

######  Lemma B.2.

If a function $f$ has bounded variation, then there exists monotone functions
$f^{+}$, $f^{-}$ with bounded variation such that $\forall t$,

$f(t)=f^{+}(t)-f^{-}(t)$

###### Proof.

Using Lemma B.2, let $\mathbb{V}^{\widehat{F}}=\mathbb{V}^{+}-\mathbb{V}^{-}$,
and note that $\mathbb{V}^{+},\mathbb{V}^{-}\in[0,d_{\infty}]$ and are
monotone functions. As they are monotone functions, let
$\\{s^{+}_{j}\\}_{j=1}^{M}$ and $\\{s^{-}_{j}\\}_{j=1}^{M}$ be the respective
stepping points for their stepwise approximations. Further, let
$\\{s^{F}_{k}\\}_{k=1}^{M^{\prime}}$ be the stepping points for the CDF $F$.
Finally, overriding notation slightly, let $\\{s_{j}\\}_{j=1}^{M^{\prime}+2M}$
be the union of these three sets of stepping points. We know that between
adjacent stepping points, $F$, $\mathbb{V}^{+}$, and $\mathbb{V}^{-}$ do not
change much. Let $\zeta^{F},\zeta^{+},\zeta^{-}$ be the step functions
constructed using the stepping points for each function, respectively, and let
$\zeta^{\mathbb{V}^{\widehat{F}}}=\zeta^{+}-\zeta^{-}$. Formally, this means
that

$\displaystyle\sup_{t\in\mathbb{R}}|F(t)-\zeta^{F}(t)|$ $\displaystyle\leq
1/2M^{\prime}$
$\displaystyle\sup_{t\in\mathbb{R}}|\mathbb{V}^{+}(t)-\zeta^{+}(t)|$
$\displaystyle\leq 1/2M$
$\displaystyle\sup_{t\in\mathbb{R}}|\mathbb{V}^{-}(t)-\zeta^{-}(t)|$
$\displaystyle\leq 1/2M$
$\displaystyle\sup_{t\in\mathbb{R}}|\mathbb{V}^{\widehat{F}}(t)-\zeta^{\mathbb{V}^{\widehat{F}}}(t)|$
$\displaystyle\leq 1/M$

At a given $s_{j}$, Bernstein’s inequality gives us

$\displaystyle|F(s_{j})-\widehat{F}(s_{j})|$
$\displaystyle\leq\frac{c_{1}\log(1/\delta)}{n}+\sqrt{\frac{c_{2}\mathbb{V}^{\widehat{F}}(s_{j})\log(1/\delta)}{n}}$
$\displaystyle=\frac{c_{1}\log(1/\delta)}{n}+\sqrt{\frac{c_{2}\left(\mathbb{V}^{+}(s_{j})-\mathbb{V}^{-}(s_{j})\right)\log(1/\delta)}{n}}$

Taking a union bound, this holds for all $\\{s_{j}\\}_{j=1}^{M^{\prime}+2M}$
with probability at least $1-\delta(M^{\prime}+2M)$.

For each $j$ and $t\in[s_{j},s_{j+1}]$, if $F(t)\geq\widehat{F}(t)$,

$\displaystyle F(t)-\widehat{F}(t)$ $\displaystyle\leq
F(s_{j+1})-\widehat{F}(s_{j})$
$\displaystyle\leq\frac{1}{M^{\prime}}+F(s_{j})-\widehat{F}(s_{j})$
$\displaystyle\leq\frac{1}{M^{\prime}}+\frac{c_{1}\log(1/\delta)}{n}+\sqrt{\frac{c_{2}\mathbb{V}^{\widehat{F}}(s_{j})\log(1/\delta)}{n}}$
$\displaystyle\leq\frac{1}{M^{\prime}}+\frac{c_{1}\log(1/\delta)}{n}+\sqrt{\frac{c_{2}\mathbb{V}^{\widehat{F}}(t)\log(1/\delta)}{n}}+\sqrt{\frac{c_{2}\log(1/\delta)}{Mn}}$

If $\widehat{F}(t)>F(t)$,

$\displaystyle\widehat{F}(t)-F(t)$
$\displaystyle\leq\widehat{F}(s_{j+1})-F(s_{j})$
$\displaystyle\leq\frac{1}{M^{\prime}}+|\widehat{F}(s_{j+1})-F(s_{j+1})|$
$\displaystyle\leq\frac{1}{M^{\prime}}+\frac{c_{1}\log(1/\delta)}{n}+\sqrt{\frac{c_{2}\mathbb{V}^{\widehat{F}}(s_{j+1})\log(1/\delta)}{n}}$
$\displaystyle\leq\frac{1}{M^{\prime}}+\frac{c_{1}\log(1/\delta)}{n}+\sqrt{\frac{c_{2}\mathbb{V}^{\widehat{F}}(t)\log(1/\delta)}{n}}+\sqrt{\frac{c_{2}\log(1/\delta)}{Mn}}$

Choosing $M,M^{\prime}\propto\sqrt{n}$ gives the resulting bound.

Similarly, for $\widehat{F}_{\textsc{S-IS}}$, we have

$\displaystyle
F-\widehat{F}_{\textsc{S-IS}}=1-S-(1-\widehat{S}_{\textsc{IS}})=\widehat{S}_{\textsc{IS}}-S$

The remainder of the proof occurs in the same manner, but with
$\mathbb{V}[\widehat{S}_{\textsc{IS}}]$ on the RHS. ∎

####  B.4.2 Empirical Bernstein Bound (Proof of Lemma 4.4)

###### Proof.

We begin with the Bernstein bounds for $\widehat{F}_{\textsc{F-IS}}(t_{j})$,
$\widehat{F}_{\textsc{S-IS}}(t_{j})$ at a given $t_{j}$. That is,

$\displaystyle|F(t_{j})-\widehat{F}_{\textsc{F-IS}}(t_{j})|$
$\displaystyle\leq\frac{c_{1}\log(1/\delta)}{n}+\sqrt{\frac{c_{2}\mathbb{V}[\widehat{F}_{\textsc{IS}}(t_{j})]\log(1/\delta)}{n}}$

and a similar statement holds for $\widehat{F}_{\textsc{S-IS}}$. From Maurer
and Pontil [2009], we have the following bound on estimation of the variance:

######  Lemma B.3 (Theorem 10 from Maurer and Pontil [2009]).

With probability at least $1-\delta$, for random variables $X\in[0,1]$,

$\left|\sqrt{\mathbb{V}_{n}}-\sqrt{\mathbb{V}}\right|\leq\sqrt{\frac{2\ln
1/\delta}{n-1}}$ (22)

Taking a union bound between Lemma B.3 and Lemma 4.3 at all $M$ points, we
have that with probability at least $1-\delta$ that for all $t_{j}$,

$\displaystyle\left|F(t_{j})-\widehat{F}_{\textsc{F-IS}}(t_{j})\right|\leq$
$\displaystyle\frac{c_{1}\log(2M/\delta)}{n}+\sqrt{\frac{c_{2}\mathbb{V}_{n}[\widehat{F}_{\textsc{IS}}(t)]\log(2M/\delta)}{n}}$
$\displaystyle+\sqrt{\frac{2c_{2}\log^{2}(2M/\delta)}{n(n-1)}}+\sqrt{\frac{c_{2}\log(2M/\delta)}{Mn}}$

Similar steps give the bound for $\widehat{F}_{\textsc{S-IS}}$, with
$\mathbb{V}_{n}[\widehat{S}_{\textsc{IS}}(t)]$ on the RHS.

∎

####  B.4.3 Empirical Bernstein Bound for Combined Estimator (Proof of
Proposition 4.1)

In order to correctly choose either $\widehat{F}_{\textsc{F-IS}}$ or
$\widehat{F}_{\textsc{S-IS}}$ for each $t_{j}$, we need that for all $t_{j}$,
the empirical variance estimates $\mathbb{V}_{n}$ are separated by at least
$\epsilon_{\mathbb{V}}(t)=\frac{1}{2}|\mathbb{V}[\widehat{F}_{\textsc{F-IS}}(t)]-\mathbb{V}[\widehat{F}_{\textsc{S-IS}}(t)]|\leq\frac{1}{2}|\mathbb{V}_{n}[\widehat{F}_{\textsc{F-IS}}(t)]-\mathbb{V}_{n}[\widehat{F}_{\textsc{S-IS}}(t)]|+\sqrt{\frac{\ln
1/\delta^{\prime}}{2(n-1)}}$. For this to be the case, we set

$\displaystyle\delta(t)=\frac{1}{2}\exp{\left(-\frac{1}{2}(n-1)\epsilon_{\mathbb{V}}(t)^{2}\right)},$

and $\delta=\max_{j}\delta(t_{j})$. Thus if

$\displaystyle
n\geq\frac{8\log(2n/\delta)}{\max_{j}\epsilon_{\mathbb{V}}(t_{j})^{2}},$

then we correctly choose $\widehat{F}_{\textsc{F-IS}}$ or
$\widehat{F}_{\textsc{S-IS}}$ for each $t_{j}$, which gives the result.

##  Appendix C Proofs for Doubly Robust CDF Estimation (Section 6)

###  C.1 Bias and Variance (Proof of Lemma 6.1)

We will prove the unbiasedness of the DR estimator by induction. As the base
case, $F_{s_{H+1}}^{0}(t)=0$ for all states $s$ and for all $t$, which is
unbiased by definition.

As the induction assumption, we assume that the DR estimator at the step $h+1$
is unbiased, that is
$\mathbb{E}_{h+1}[\widehat{F}_{S_{h+1}}(t)]=\mathbb{E}_{h+1}[F_{S_{h+1}}](t)$.

We need to prove that it is also unbiased for the $h$th step. Recall for now
we assume that rewards are deterministic. Then for any $t$:

$\displaystyle\mathbb{E}^{\beta}_{h}[\widehat{F}_{S_{h}}(t)]$
$\displaystyle=\mathbb{E}_{\mathbb{P}_{\beta}}\left[\overline{F}^{H+1-h}_{S_{h}}(t)+w(A_{h},S_{h})\left(\widehat{F}^{H-h}_{S_{h+1}}\left(\frac{t-R(S_{h},A_{h})}{\gamma}\right)-\overline{F}^{H+1-h}_{S_{h},A_{h}}(t)\right)\right]$
$\displaystyle=\mathbb{E}^{\beta}_{h}\left[\overline{F}^{H+1-h}_{S_{h}}(t)-w(A_{h},S_{h})\overline{F}^{H+1-h}_{S_{h},A_{h}}(t)\right]+\mathbb{E}^{\beta}_{h}\left[w(A_{h},S_{h})\widehat{F}^{H-h}_{S_{h+1}}\left(\frac{t-R(S_{h},A_{h})}{\gamma}\right)\right]$
$\displaystyle=\mathbb{E}_{h}\left[\overline{F}^{H+1-h}_{Z(S_{h})}(t)-\mathbb{E}_{h}\left[\overline{F}^{H+1-h}_{S_{h},A_{h}}(t)\Big{|}S_{h}\right]\right]$
$\displaystyle\quad+\mathbb{E}_{h}\left[\widehat{F}^{H-h}_{S_{h+1}}\left(\frac{t-R(S_{h},A_{h})}{\gamma}\right)\right]$
$\displaystyle=\mathbb{E}_{h}\left[\overline{F}^{H+1-h}_{Z(s_{h})}(t)-\overline{F}^{H+1-h}_{Z(s_{h})}(t)\right]+\mathbb{E}_{h}\left[\mathbb{E}_{h+1}\left[\widehat{F}^{H-h+1}_{S_{h}+1}\left(\frac{t-R(S_{h},A_{h})}{\gamma}\right)\right]\right]$
$\displaystyle=\mathbb{E}_{h}\left[\mathbb{E}_{h+1}\left[F^{H-h+1}_{S_{h}+1}\left(\frac{t-R(S_{h},A_{h})}{\gamma}\right)\right]\right]$
$\displaystyle=\mathbb{E}_{h}\left[F_{S_{h}}(t)\right]$

where the third equality incorporates the importance sampling weight $w$ into
a change of measure from $\beta$ to $\pi$, the foruth equality uses the
conditional expectation. The second to last equality uses the induction
assumption and the last uses the recursive identity.

Next, we derive the variance of the DR estimator.

$\displaystyle\mathbb{V}_{h}^{\beta}\left[\widehat{F}_{S_{h}}(t)\right]=$
$\displaystyle\mathbb{E}_{h}^{\beta}\left[\left(\widehat{F}_{S_{h}}(t)\right)^{2}\right]-\mathbb{E}_{h}^{\beta}\left[\widehat{F}_{S_{h}}(t)\right]^{2}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{h}^{\beta}\left[\left(\overline{F}^{H+1-h}_{S_{h}}(t)+w(A_{h},S_{h})\left(\widehat{F}^{H-h}_{S_{h+1}}\left(\frac{t-R(S_{h},A_{h})}{\gamma}\right)-\overline{F}^{H+1-h}_{S_{h},A_{h}}(t)\right)\right)^{2}-F_{S_{h}}(t)^{2}\right]$
$\displaystyle+\mathbb{E}_{h}^{\beta}\left[F_{S_{h}}(t)^{2}\right]-\mathbb{E}_{h}^{\beta}\left[F_{S_{h}}(t)\right]^{2}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{h}^{\beta}\left[\left(\overline{F}^{H+1-h}_{S_{h}}(t)+w(A_{h},S_{h})\left(\widehat{F}^{H-h}_{S_{h+1}}\left(\frac{t-R(S_{h},A_{h})}{\gamma}\right)-\overline{F}^{H+1-h}_{S_{h},A_{h}}(t)\right)\right)^{2}-F_{S_{h}}(t)^{2}\right]$
$\displaystyle+\mathbb{V}_{h}\left[F_{S_{h}}(t)\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}_{h}^{\beta}\bigg{[}\bigg{(}\overline{F}^{H+1-h}_{S_{h}}(t)+w(A_{h},S_{h})F_{S_{h},A_{h}}(t)-w(A_{h},S_{h})\overline{F}^{H+1-h}_{S_{h},A_{h}}$
$\displaystyle+w(A_{h},S_{h})\left(\widehat{F}_{S_{h+1}}\left(\frac{t-R(S_{h},A_{h})}{\gamma}\right)-F_{S_{h},A_{h}}(t)\right)\bigg{)}^{2}-F_{S_{h}}(t)^{2}\bigg{]}+\mathbb{V}_{h}\left[F_{S_{h}}(t)\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{h}^{\beta}\bigg{[}\bigg{(}\overline{F}^{H+1-h}_{S_{h}}(t)-w(A_{h},S_{h})\Delta_{S_{h},A_{h}}(t)$
$\displaystyle+w(A_{h},S_{h})\left(\widehat{F}_{S_{h+1}}\left(\frac{t-R(S_{h},A_{h})}{\gamma}\right)-\mathbb{E}_{h+1}\left[\widehat{F}_{S_{h+1}}\left(\frac{t-R(S_{h},A_{h})}{\gamma}\right)\right]\right)\bigg{)}^{2}-F_{S_{h}}(t)^{2}\bigg{]}$
$\displaystyle+\mathbb{V}_{h}\left[F_{S_{h}}(t)\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}_{h}^{\beta}\left[\left(\overline{F}^{H+1-h}_{S_{h}}(t)-w(A_{h},S_{h})\Delta_{S_{h},A_{h}}(t)\right)^{2}-F_{S_{h}}(t)^{2}\right]$
$\displaystyle+\mathbb{E}_{h}\left[\left(w(A_{h},S_{h})\left(\widehat{F}_{S_{h+1}}\left(\frac{t-R(S_{h},A_{h})}{\gamma}\right)-\mathbb{E}_{h+1}\left[\widehat{F}_{S_{h+1}}\left(\frac{t-R(S_{h},A_{h})}{\gamma}\right)\right]\right)\right)^{2}\right]$
$\displaystyle+\mathbb{V}_{h}\left[F_{S_{h}}(t)\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}_{h}^{\beta}\left[\mathbb{V}_{h}^{\beta}\left[w(A_{h},S_{h})\Delta_{S_{h},A_{h}}(t)\Big{|}S_{h}\right]\right]+\mathbb{E}_{h}^{\beta}\left[w(A_{h},S_{h})^{2}\mathbb{V}_{h+1}\left[\widehat{F}_{S_{h+1}}\left(\frac{t-R(S_{h},A_{h})}{\gamma}\right)\Big{|}S_{h},A_{h}\right]\right]$
$\displaystyle+\mathbb{V}_{h}\left[F_{S_{h}}(t)\right]$

where $\Delta_{s,a}(t)=\overline{F}_{s,a}(t)-F_{s,a}(t)$.

###  C.2 Error Bound on Doubly Robust Estimator (Proof of Lemma 6.2)

We separate the bound into two parts using the definition of
$\widehat{F}_{\textsc{DR}}$:

$\displaystyle\left\|F-\widehat{F}_{\textsc{DR}}\right\|_{\infty}\leq$
$\displaystyle\left\|F-\frac{1}{n}\sum_{i=1}^{n}w(\tau^{i})\mathbbm{1}\left\\{\sum_{h=1}^{H}\gamma^{h}r_{h}\leq
t\right\\}\right\|_{\infty}$
$\displaystyle+\left\|\frac{1}{n}\sum_{i=1}^{n}w(\tau_{1:h-1}^{i})\left(w(s^{i}_{h},a^{i}_{h})\overline{F}_{S_{h},A_{h}}\left(\frac{t-\sum_{k=1}^{h}R_{k}^{i}}{\gamma^{h}}\right)-\overline{F}_{S_{h}}\left(\frac{t-\sum_{k=1}^{h}R_{k}^{i}}{\gamma^{h}}\right)\right)\right\|_{\infty}$
(23)

Lemma 4.2 upper bounds the first term, and we are left to bound the second.

Let
$\tau=(s^{0},a^{0},r^{0},\ldots,s^{h},a^{h},r^{h},\ldots,s^{H},a^{H},r^{H},s^{H+1})$
be a trajectory, and define $w(\tau_{1:h})=\prod_{k=1}^{h}w(s_{h},a_{h})$.
Define the function class

$\displaystyle\mathbb{F}(m,w)=\left\\{f\left(\\{\\{x_{j}^{h}\\}_{j=1}^{M}\\}_{h=1}^{H},\tau\right):=\frac{\varrho}{m}\sum_{h=1}^{H}w(\tau_{1:h-1})\sum_{j=1}^{m}w(a^{h},s^{h})\mathbbm{1}_{\\{x_{j}^{h}\leq
t\\}}:\forall
t\in\mathbb{R},\varrho\in\\{-1,+1\\},\\{x_{j}^{h}\\}_{j,h=1}^{M,H}\in\mathbb{Q}^{M\times
H}\right\\}$

Note that (23) is

$\displaystyle\leq\sup_{\zeta\in\mathbb{F}(m,w)}\left|\frac{1}{n}\sum_{i=1}^{n}\zeta\left(\\{x^{j}_{S_{i}^{h},A_{i}^{h}}\\},\tau_{i}\right)-\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{\mathbb{P}_{\beta}}\left[\zeta\left(\\{x^{j}_{S_{i}^{h},A_{i}^{h}}\\},\tau_{i}\right)\Big{|}\\{S_{i}^{h}\\}_{h=1}^{H}\right]\right|+\frac{1}{m}$

Going forward, we refer to
$\zeta\left(\\{x^{j}_{S_{i}^{h},A_{i}^{h}}\\},\tau\right)$ as
$\zeta(\tau_{i})$ for short. For $\lambda>0$ we have

$\displaystyle\mathbb{E}_{\mathbb{P}_{\beta}}\left[\exp\left(\lambda\sup_{\zeta\in\mathbb{F}(m,w)}\left(\frac{1}{n}\sum_{i=1}^{n}\zeta(\tau_{i})-\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{\mathbb{P}_{\beta}}\left[\zeta(\tau_{i})\Big{|}\\{S^{h}_{i}\\}_{h=1}^{H}\right]\right)\right)\right]$
$\displaystyle=\mathbb{E}_{\mathbb{P}_{\beta}}\left[\exp\left(\lambda\sup_{\zeta\in\mathbb{F}(m,w)}\left(\frac{1}{n}\sum_{i=1}^{n}\zeta(\tau_{i})-\frac{1}{n}\sum_{i=1}^{n}\zeta(\tau^{\prime}_{i})\right)\right)\right]$

where $\tau^{\prime}_{i}=\\{s_{i}^{h},a^{{}^{\prime}h}_{i}\\}_{h=1}^{H}$.
Using symmetrization, this is

$\displaystyle\leq\mathbb{E}_{\mathbb{P}_{\beta},\Re}\left[\exp\left(2\lambda\sup_{\zeta\in\mathbb{F}(m,w)}\frac{1}{n}\sum_{i=1}^{n}\xi_{i}\zeta(\tau_{i})\right)\right]$
$\displaystyle=\mathbb{E}_{\mathbb{P}_{\beta},\Re}\left[\sup_{\zeta\in\mathbb{F}(m,w)}\exp\left(2\lambda\frac{1}{n}\sum_{i=1}^{n}\xi_{i}\zeta(\tau_{i})\right)\right]$
$\displaystyle=\mathbb{E}_{\mathbb{P}_{\beta},\Re}\left[\sup_{\varrho,t}\exp\left(2\lambda\frac{\varrho}{nm}\sum_{i=1}^{n}\xi_{i}\sum_{h=1}^{H}w(\tau_{i}^{1:h-1})\sum_{j=1}^{m}w(A_{i}^{h},S_{i}^{h})\mathbbm{1}_{\\{x^{j}_{S_{i}^{h},A_{i}^{h}}\leq
t\\}}\right)\right]$

Now for each $h,j$, permute the indices of $i$ with a new indexing $(i)$ such
that

$x^{j}_{S_{(1)}^{h},A_{(1)}^{h}}\leq\ldots\leq
x^{j}_{S_{(i)}^{h},A_{p(i)}^{h}}\leq\ldots\leq
x^{j}_{S_{(n)}^{h},A_{(n)}^{h}}.$

Then we have

$\displaystyle\mathbb{E}_{\mathbb{P}_{\beta},\Re}\left[\sup_{\varrho,t}\exp\left(2\lambda\frac{\varrho}{nm}\sum_{i=1}^{n}\xi_{i}\sum_{h=1}^{H}w(\tau_{i}^{1:h-1})\sum_{j=1}^{m}w(A_{i}^{h},S_{i}^{h})\mathbbm{1}_{\\{x^{j}_{S_{i}^{h},A_{i}^{h}}\leq
t\\}}\right)\right]$
$\displaystyle=\mathbb{E}_{\mathbb{P}_{\beta},\Re}\left[\max_{\varrho,k}\exp\left(2\lambda\frac{\varrho}{nm}\sum_{h=1}^{H}\sum_{j=1}^{m}\sum_{i=1}^{k}\xi_{(i)}w(\tau_{(i)}^{1:h-1})w(A_{(i)}^{h},S_{(i)}^{h})\right)\right]$
$\displaystyle\leq\mathbb{E}_{\mathbb{P}_{\beta},\Re}\left[\max_{\varrho,j,k}\exp\left(2\lambda\frac{\varrho}{n}\sum_{i=1}^{k}\xi_{(i)}\sum_{h=1}^{H}w(\tau_{(i)}^{1:h-1})w(A_{(i)}^{h},S_{(i)}^{h})\right)\right]$
$\displaystyle\leq
2\sum_{j=1}^{m}\mathbb{E}_{\mathbb{P}_{\beta},\Re}\left[\exp\left(\frac{2\lambda}{n}\sum_{i=1}^{n}\xi_{j(i)}\sum_{h=1}^{H}w(\tau_{j(i)}^{1:h})\right)\mathbbm{1}_{\\{\sum_{i=1}^{n}\xi_{j(i)}\sum_{h=1}^{H}w(\tau_{j(i)}^{1:h})\geq
0\\}}\right]$ $\displaystyle\leq
2m\mathbb{E}_{\mathbb{P}_{\beta},\Re}\left[\exp\left(\frac{2\lambda}{n}\sum_{i=1}^{n}\xi_{j(i)}\sum_{h=1}^{H}w(\tau_{j(i)}^{1:h})\right)\right]$
$\displaystyle\leq
2m\exp\left(\frac{2\lambda^{2}w_{max}(w_{max}^{H}-1)}{n(w_{max}-1)}\right)$

And we have the result,

$\displaystyle\mathbb{P}_{\beta}\left(\cdot\geq\epsilon+\frac{1}{m}\right)\leq
4m\exp\left(\frac{2\lambda^{2}w_{max}(w_{max}^{H}-1)}{n(w_{max}-1)}-\lambda\epsilon\right)$
$\displaystyle\leq
4m\min_{\lambda>0}\exp\left(\frac{2\lambda^{2}w_{max}(w_{max}^{H}-1)}{n(w_{max}-1)}-\lambda\epsilon\right)$
$\displaystyle=4m\exp\left(\frac{-n\epsilon(w_{max}-1)^{2}}{w_{max}^{2}(w_{max}^{H}-1)^{2}}\right)$

where the second line results because the inequality holds for all
$\lambda>0$, so we can minimize over $\lambda$.

Then we have

$\displaystyle\mathbb{P}_{\beta}\left(\cdot\geq\sqrt{\frac{w_{max}^{2}(w_{max}^{H}-1)^{2}}{n(w_{max}-1)^{2}}\log\frac{4m}{\delta}}+\frac{1}{m}\right)\leq\delta$

Combining this with Lemma 4.2 to bound the first term gives the result and
choice of $m\propto\sqrt{n}$ gives the result.

###  C.3 Cramer-Rao Variance Lower Bound of Off-Policy CDF Estimation (Proof
of Theorem 6.1)

We prove a Cramer-Rao lower bound for UDAGs, defined below.

######  Definition C.1 (Unique Directed Acyclic Graph).

An MDP is unique directed acyclic graph (UDAG) if the state and action spaces
are finite. A state only occurs at a particular horizon $h$, that is, for any
$s\in\mathcal{S}$, there exists a unique $h$ such that
$\max_{\pi}P(s_{h}=s|\pi)>0$. In addition, each state has a unique reward
function. For simplicity, assume $\gamma=1$.

###### Proof.

First, we transform the MDP of Definition C.1 into one with an augmented state
space, including the history of rewards encountered. That is, a trajectory is
now defined as

$\tau=(s_{0},g_{0},a_{0},...s_{H},g_{H},a_{H},s_{H+1},g_{H+1}),$

where $g_{h}=\sum_{k=0}^{h-1}r_{h}$. Note that
$P(g_{h+1}|s_{h},g_{h},a_{h})=P(r_{h}|s_{h},a_{h})$ and
$P(s_{h+1}|s_{h},g_{h},a_{h})=P(s_{h+1}|s_{h},a_{h})$.

For brevity, denote $\widetilde{s}=(s,g)$ to be the augmented state space. We
now derive a pointwise Constrained Cramer-Rao Bound (CCRB) for $F(t)$ in the
augmented MDP,

$KU(UIU)^{-1}U^{\top}K^{\top}$

where $K$ is the Jacobian of $F(t)$ and $I$ is the Fisher information matrix.
$U$ is a matrix whose column vectors consist of an orthonormal basis for the
null space of a matrix $N$, which is a block-diagonal matrix with
$N_{(\widetilde{s},a),(\widetilde{s}^{\prime},a^{\prime})}=\mathbbm{1}\left\\{\widetilde{s},a=\widetilde{s}^{\prime},\widetilde{a}^{\prime}\right\\}$.
Further, define
$\eta_{\widetilde{s},a,\widetilde{s}^{\prime},g^{\prime}}=P(\widetilde{s}^{\prime}|\widetilde{s},a)$,
and note that $N\eta=\mathbf{1}$.

To compute the Fisher information matrix $I$,

$I=\mathbb{E}\left[\left(\frac{\partial\log
P_{\beta}(\tau)}{\partial\eta}\right)\left(\frac{\partial\log
P_{\beta}(\tau)}{\partial\eta}\right)^{\top}\right]$

where $P_{\beta}(\tau)$ is the probability of a trajectory $\tau$ under the
behavioral policy $\beta$,

$\displaystyle P_{\beta}(\tau)$
$\displaystyle=\mu(s_{1},g_{1})\beta(a_{1}|s_{1})P(s_{2},|s_{1},a_{1})P(g_{2}|s_{1},a_{1},g_{1})\ldots
P(s_{H},g_{H}|s_{H-1},g_{H-1}a_{H-1})\beta(a_{H}|s_{H})P(s_{H+1},g_{H+1}|s_{H},g_{H},a_{H})$
$\displaystyle=\mu(\widetilde{s}_{1})\beta(a_{1}|\widetilde{s}_{1})P(\widetilde{s}_{2},|\widetilde{s}_{1},a_{1})P(\widetilde{s}_{2}|\widetilde{s}_{1},a_{1})\ldots
P(\widetilde{s}_{H}|\widetilde{s}_{H-1}a_{H-1})\beta(a_{H}|\widetilde{s}_{H})P(\widetilde{s}_{H+1}|\widetilde{s}_{H},a_{H})$

Define a new indicator vector $g(\tau)$ with
$g(\tau)_{\widetilde{s}_{h},a_{h},\widetilde{s}_{h+1}}=1$ if
$(\widetilde{s}_{h},a_{h},\widetilde{s}_{h+1})\in\tau$. Then, letting $\circ$
denote an elementwise operation, we have

$\frac{\partial\log P_{\beta}(\tau)}{\partial\eta}=\eta^{\circ-1}\circ
g(\tau)$

which makes

$I=\mathbb{E}\left[[1/\eta_{i}\eta_{j}]_{ij}\circ\left(g(\tau)g(\tau)^{\top}\right)\right]=[1/\eta_{i}\eta_{j}]_{ij}\circ\mathbb{E}\left[g(\tau)g(\tau)^{\top}\right]$

As a result, the elements of $I$ are as follows, where $P_{M}$ refers to the
marginal probability:

  * •

diagonal:
$\frac{P_{M}(\widetilde{s}_{h},a_{h})}{P(\widetilde{s}_{h+1}|\widetilde{s}_{h},a_{h})}$

  * •

row $=(\widetilde{s}_{h},a_{h},\widetilde{s}_{h+1})$, column
$=(\widetilde{s}^{\prime}_{h},a^{\prime}_{h},\widetilde{s}^{\prime}_{h+1})$:
$P_{M}(\widetilde{s}_{h}^{\prime},a_{h}^{\prime})P_{M}(\widetilde{s}_{h},a_{h}|\widetilde{s}_{h+1}^{\prime})$

  * •

otherwise: 0

We know calculate the term $(U^{\top}IU)^{-1}$, and we use a similar strategy
to Jiang and Li [2016], Huang and Jiang [2020] to avoid taking its inverse.
Note that

$U^{\top}IU=U^{\top}\left(N^{\top}X^{\top}+I+XN\right)U.$

where $X$ is arbitrary, $N$ was defined previously. Let
$D=N^{\top}X^{\top}+I+XN$. Our goal is to define $X$ such that $D$ is a
diagonal matrix with diagonal identical to the diagonal of $I$. Note that $XN$
and $N^{\top}X^{\top}$ are symmetric, so we can design the former to eliminate
the upper triangle of $I$, and the latter to eliminate the lower triangle. To
verify we can do this, set

  * •

row $=(\widetilde{s}_{h},a_{h},\widetilde{s}_{h+1})$, column
$=(\widetilde{s}_{h},a_{h},2)$ : $-P_{M}(\widetilde{s}_{h},a_{h})$

  * •

row $=(\widetilde{s}_{h},a_{h},\widetilde{s}_{h+1})$, column
$=(\widetilde{s}_{h}^{\prime},a_{h}^{\prime},\\{1\;\text{or}\;2\\})$ :
$-P_{M}(\widetilde{s}_{h},a_{h})P_{M}(\widetilde{s}_{h}^{\prime},a_{h}^{\prime},\widetilde{s}_{h+1})\mathbbm{1}\left\\{h<h^{\prime}\right\\}$

  * •

otherwise: 0

Then with the proper choice of $U$,

$U(U^{\top}IU)^{-1}U^{\top}=\text{Diag}(\left\\{B(\widetilde{s}_{h},a_{h})\right\\}_{0}^{H})$

where $\text{Diag}(\left\\{\cdot\right\\})$ is a block-diagonal matrix
consisting of matrices in the set $\left\\{\cdot\right\\}$, and

$B(\widetilde{s}_{h},a_{h})=\frac{\text{Diag}(P(\cdot|\widetilde{s}_{h},a_{h}))-P(\cdot|\widetilde{s}_{h},a_{h})P(\cdot|\widetilde{s}_{h},a_{h})^{\top}}{P_{M}(\widetilde{s}_{h},a_{h})}$

where $P(\cdot|\widetilde{s}_{h},a_{h})$ is the transition vector.

We now need to calculate the Jacobian $K$. Recall that the estimation
objective is

$\displaystyle F(t)$
$\displaystyle=\sum_{\widetilde{s}_{1}}\mu(\widetilde{s}_{1})\sum_{a_{1}}\pi(a_{1}|\widetilde{s}_{1})\ldots\sum_{s_{H}}P(\widetilde{s}_{H}|\widetilde{s}_{H-1},a_{H-1})\mathbbm{1}\left\\{g_{H}\leq
t\right\\}$

Then $K(\widetilde{s}_{h},a_{h},\widetilde{s}_{h+1})$ is

$\displaystyle\frac{\partial F(t)}{\partial P(r_{h},s_{h+1}|s_{h},a_{h})}$
$\displaystyle=\frac{\partial F(t)}{\partial
P(s_{h+1},g_{h+1}|s_{h},g_{h},a_{h})}=\mathbb{P}(s_{h},g_{h},a_{h})P(r_{h}|s_{h},a_{h})P(s_{h+1}|s_{h},a_{h})F_{s_{h+1}}(t-g_{h}-r_{h})$

Denote by $K_{(s_{h},a_{h},\cdot)}$ the vector fragment of $K$ whose index
tuple starts with $s_{h},a_{h}$. Then we have that

$\displaystyle KU(UIU)^{-1}U^{\top}K^{\top}$
$\displaystyle=\sum_{h=1}^{H+1}\sum_{\widetilde{s}_{h},a_{h}}\left(K_{(\widetilde{s}_{h},a_{h},\cdot)}\right)^{\top}B(\widetilde{s}_{h},a_{h})K_{(\widetilde{s}_{h},a_{h},\cdot)}$
(24)

We have

$\displaystyle(K_{(s_{h},g_{h},a_{h},\cdot)})^{\top}B(s_{h},g_{h},a_{h})K_{(s_{h},g_{h},a_{h},\cdot)}$
$\displaystyle=\frac{\mathbb{P}(s_{h},g_{h},a_{h})^{2}}{\mathbb{P}_{\beta}(s_{h},g_{h},a_{h})}\Bigg{(}\sum_{r_{h}}P(r_{h}|s_{h},a_{h})\sum_{s_{h+1}}P(s_{h+1}|s_{h},a_{h})F_{s_{h+1}}(t-g_{h}-r_{h})^{2}$
$\displaystyle\quad-\left(\sum_{r_{h}}P(r_{h}|s_{h},a_{h})\sum_{s_{h+1}}P(s_{h+1}|s_{h},a_{h})F_{s_{h+1}}(t-g_{h}-r_{h})\right)^{2}\Bigg{)}$
$\displaystyle=\frac{\mathbb{P}(s_{h},g_{h},a_{h})^{2}}{\mathbb{P}_{\beta}(s_{h},g_{h},a_{h})}\mathbb{V}_{r_{h},s_{h+1}|s_{h},a_{h}}\left[F_{s_{h+1}}(t-g_{h+1})\right]$

Under Definition C.1, we have that $\mathbb{P}(g_{h})=\mathbb{P}(\tau_{h})$ so
the lower bound is

$\displaystyle\sum_{h=1}^{H}\mathbb{E}_{\mathbb{P}_{\beta}}\left[w_{1:h-1}^{2}\mathbb{V}_{r_{h},s_{h+1}|s_{h},a_{h}}\left[F_{s_{h+1}}\left(t-\sum_{k=1}^{h}r_{k}\right)\right]\right].$

∎

##### Comparison with DR Variance.

When expanded over $H$ horizons, the DR variance (6.1) in MDP of Definition
C.1 is

$\displaystyle\sum_{h=1}^{H}\mathbb{E}_{\mathbb{P}_{\beta}}\left[w_{1:h-1}^{2}\mathbb{V}_{h}\left[w(A_{h},S_{h})\Delta_{S_{h},A_{h}}\left(t-\sum_{k=1}^{h}R_{k}\right)\Big{|}S_{h}\right]\right]+\sum_{h=1}^{H}\mathbb{E}_{\mathbb{P}_{\beta}}\left[w_{1:h-1}^{2}\mathbb{V}_{h}\left[F_{S_{h}}\left(t-\sum_{k=1}^{h}R_{k}\right)\right]\right]$

where $\Delta_{s,a}(t)=\overline{F}_{s,a}(t)-F(t)$. Thus when
$\overline{F}=F$, the first term of the above expression goes to 0 and the DR
variance becomes

$\displaystyle\sum_{h=1}^{H}\mathbb{E}_{\mathbb{P}_{\beta}}\left[w_{1:h-1}^{2}\mathbb{V}_{h}\left[F_{S_{h}}\left(t-\sum_{k=1}^{h}R_{k}\right)\right]\right]$

which attains the Cramer-Rao lower bound.

##  Appendix D Proofs for Risk Estimation (Section 7)

###  D.1 CVaR Risk Lower Bound (Proof of Theorem 7.1)

To prove the lower bound for CVaR estimation, we make an adaptation of the
proof in B.3. From Le Cam’s method, we have

$\inf_{\widehat{F}}\sup_{P_{1},P_{2}\in\mathcal{P}}\mathbb{E}_{P}\left[|\widehat{\rho}-\rho_{P}^{\pi}|\right]\geq\frac{1}{8}|\rho_{1}^{\pi}-\rho_{2}^{\pi}|e^{-n\text{D}_{\text{KL}}(P_{1}\|P_{2})}$
(25)

For a distortion risk functional $\rho$,

$\displaystyle\left|\rho_{1}^{\pi}-\rho_{2}^{\pi}\right|$
$\displaystyle=\left|\int_{0}^{D}g(1-F_{1}^{\pi}(t))-g(1-F_{2}^{\pi}(t))dt\right|$

In the case of $\text{CVaR}_{\alpha}$, with $\alpha\in[0,1)$, we have
$g(x)=\min\\{\frac{x}{1-\alpha},1\\}$ so

$\displaystyle\left|\text{CVaR}_{\alpha}(F_{1}^{\pi})-\text{CVaR}_{\alpha}(F_{2}^{\pi})\right|$
$\displaystyle=\left|\int_{0}^{D}\min\left\\{\frac{1-F_{1}^{\pi}(t)}{1-\alpha},1\right\\}-\min\left\\{\frac{1-F_{2}^{\pi}(t)}{1-\alpha},1\right\\}dt\right|$

Note that, under this definition, the expected return is CVaR with level
$\alpha=0$.

Recall the Bernoulli problem instances we previously used to prove the lower
bound for MABs. For the MAB $P_{2}$, we set the arm reward distribution to be
$R_{2}(0|a)=c$ for all $a\in\mathcal{A}$. For the MAB $P_{1}$, we set the arm
reward distribution to be $R_{1}(0|a)=c+\delta(a)$ for some
$\delta(a)\in[0,1-c]$. Then

$\displaystyle\left|\text{CVaR}_{\alpha}(F_{1}^{\pi})-\text{CVaR}_{\alpha}(F_{2}^{\pi})\right|$
$\displaystyle=\left|\int_{0}^{1}\min\left\\{\frac{1-F_{1}^{\pi}(t)}{1-\alpha},1\right\\}-\min\left\\{\frac{1-F_{2}^{\pi}(t)}{1-\alpha},1\right\\}dt\right|$
$\displaystyle=\left|\min\left\\{\frac{1-F_{1}^{\pi}(t=0)}{1-\alpha},1\right\\}-\min\left\\{\frac{1-F_{2}^{\pi}(t=0)}{1-\alpha},1\right\\}\right|$
$\displaystyle=\left|\min\left\\{\frac{1-c-\sum_{a}\pi(a)\delta(a)}{1-\alpha},1\right\\}-\min\left\\{\frac{1-c}{1-\alpha},1\right\\}\right|$

Then note that the difference above has three possibilities:

$=\begin{cases}0,&\text{if}\;1-\alpha\leq 1-c-\sum_{a}\pi(a)\delta(a)\\\
\frac{|\alpha-c|}{1-\alpha},&\text{if}\;1-c-\sum_{a}\pi(a)\delta(a)\leq
1-\alpha\leq 1-c\\\
\frac{\sum_{a}\pi(a)\delta(a)}{1-\alpha},&\text{if}\;1-c\leq
1-\alpha\end{cases}$

By setting $c=\alpha$ and $\delta(a)\in[0,1-\alpha]$, it can be seen that the
first scenario (a vacuous lower bound) is impossible, and that the third
scenario encompasses the second scenario. That is, if $\delta(a)=0\;\forall
a$, then both the second and third cases become 0.

Then since

$\displaystyle\text{D}_{\text{KL}}(P_{1}||P_{2})=\sum_{a}\beta(a)\text{D}_{\text{KL}}(R_{1}(a)||R_{2}(a))\leq
4\sum_{a}\beta(a)\delta(a)^{2},$

Our optimization problem boils down to

$\max_{\delta\in\mathbb{R}^{\mathcal{A}}}\frac{\sum_{a}\pi(a)\delta(a)}{1-\alpha}\quad\text{s.t.}\quad\sum_{a}\beta(a)\delta(a)^{2}\leq\frac{1}{8n},\;0\leq\delta\leq
1-\alpha$

Thus we have the Lagrangian

$\displaystyle L(\delta,\lambda,\nu)$
$\displaystyle=-\frac{\sum_{a}\pi(a)\delta(a)}{1-\alpha}+\lambda\left(\sum_{a}\beta(a)\delta(a)^{2}-\frac{1}{8n}\right)+\sum_{a}\nu(a)\left(\delta(a)-1+\alpha\right)$
$\displaystyle=-\frac{\lambda}{8n}+\sum_{a}\left\\{\lambda\beta(a)\delta^{2}(a)+\delta(a)\left(\nu(a)-\frac{\pi(a)}{1-\alpha}\right)-(1-\alpha)\nu(a)\right\\}$

The optimal solution is

$\displaystyle\delta^{*}(a)=\frac{\frac{\pi(a)}{1-\alpha}-\nu(a)}{2\lambda\beta(a)}$

When we plug this in, we have

$\displaystyle-\frac{\sum_{a}\pi(a)\delta^{*}(a)}{1-\alpha}$
$\displaystyle=\max_{\lambda\geq 0,\nu\geq 0}g(\lambda,\nu)=\max_{\nu\geq
0}\left\\{-\sqrt{\frac{1}{8n}\sum_{a}\frac{\left(\frac{\pi(a)}{1-\alpha}-\nu(a)\right)^{2}}{\beta(a)}}-(1-\alpha)\sum_{a}\nu(a)\right\\}$
$\displaystyle=-\min_{\nu\geq
0}\left\\{\sqrt{\frac{1}{8n}\sum_{a}\frac{\left(\frac{\pi(a)}{1-\alpha}-\nu(a)\right)^{2}}{\beta(a)}}+(1-\alpha)\sum_{a}\nu(a)\right\\}$
(26)

Finally, we have from Lemma D.1 that

$\displaystyle\min_{\nu\geq
0}\left\\{\sqrt{\frac{1}{8n}\sum_{a}\frac{\left(\frac{\pi(a)}{1-\alpha}-\nu(a)\right)^{2}}{\beta(a)}}+(1-\alpha)\sum_{a}\nu(a)\right\\}\asymp\frac{1}{1-\alpha}\left\\{\sum_{a\in
S^{*}}\pi(a)+\sqrt{\frac{\sum_{a\notin S}\beta(a)w(a)^{2}}{n}}\right\\}$

######  Lemma D.1.

$\displaystyle\min_{\nu\geq
0}\left\\{\sqrt{\frac{1}{8n}\sum_{a}\frac{\left(\frac{\pi(a)}{1-\alpha}-\nu(a)\right)^{2}}{\beta(a)}}+(1-\alpha)\sum_{a}\nu(a)\right\\}\asymp\sum_{a\in
S^{*}}\pi(a)+\frac{1}{1-\alpha}\sqrt{\frac{\sum_{a\notin
S}\beta(a)w(a)^{2}}{n}}$ (27)

###### Proof.

First, note that $0\leq\nu^{*}\leq\pi(a)$ from the optimization problem (26).
Further, WLOG we can assume that $\pi(a)>0$ and $\beta(a)>0$. This is because
if $\pi(a)=0$, then $w(a)=0$, which contributes 0 to both sides of the
equation. Similarly, if $\beta(a)=0$, then $\nu^{*}(a)=\pi(a)$, which means
$a\in S^{*}$, and $\pi(a)$ is contributed to both sides of the equation.

We separate the proof into three cases: $v^{*}=\frac{\pi}{1-\alpha}$,
$v^{*}=0$, $0\neq v^{*}\neq\frac{\pi}{1-\alpha}$.

#####  $\mathbf{v^{*}=\pi}$.

In this case, $S^{*}=[k]$, and the LHS of (27) is

$\displaystyle(1-\alpha)\sum_{a}\frac{\pi(a)}{1-\alpha}=1$

while the RHS is

$\displaystyle\sum_{a}\pi(a)=1.$

#####  $\mathbf{v^{*}=0}$.

In this case, $S^{*}=\\{\\}$, and the LHS of (27) is

$\displaystyle\sqrt{\frac{1}{8n}\sum_{a}\frac{\frac{\pi(a)}{1-\alpha}^{2}}{\beta(a)}}$
$\displaystyle=\frac{1}{1-\alpha}\sqrt{\frac{1}{8n}\sum_{a}\beta(a)w(a)^{2}},$

while the RHS is

$\displaystyle\frac{1}{1-\alpha}\sqrt{\frac{\sum_{a}\beta(a)w(a)^{2}}{n}}.$

##### $\mathbf{0\neq v^{*}\neq\pi}.$

First, we write the optimality conditions of (26):

$\displaystyle\left(\frac{w(a)}{1-\alpha}-\frac{\nu^{*}(a)}{\beta(a)}\right)^{2}=8n(1-\alpha)^{2}\sum_{a}\frac{\left(\frac{\pi(a)}{1-\alpha}-\nu^{*}(a)\right)^{2}}{\beta(a)}\quad\quad$
$\displaystyle\text{for}\;a\in S^{*}$ (28) $\displaystyle w(a)^{2}\leq
8n(1-\alpha)^{4}\sum_{a}\frac{\left(\frac{\pi(a)}{1-\alpha}-\nu^{*}(a)\right)^{2}}{\beta(a)}\quad\quad$
$\displaystyle\text{for}\;a\notin S^{*}$ (29)

Denote

$\displaystyle T_{1}$ $\displaystyle:=\sum_{a\in
S^{*}}\frac{\left(\frac{\pi(a)}{1-\alpha}-\nu^{*}(a)\right)^{2}}{\beta(a)}$
$\displaystyle T_{2}$ $\displaystyle:=\sum_{a\notin
S^{*}}\frac{\left(\frac{\pi(a)}{1-\alpha}-\nu^{*}(a)\right)^{2}}{\beta(a)}=\frac{1}{(1-\alpha)^{2}}\sum_{a\notin
S^{*}}\beta(a)w(a)^{2}$

Next, we will show that $T_{1}=\frac{1-\epsilon}{\epsilon}T_{2}$ for some
$\epsilon\in(0,1)$ such that $\beta(S^{*})=(1-\epsilon)/8(1-\alpha)^{2}n$.

Summing (28) over $a\in S^{*}$, we have

$\displaystyle T_{1}$ $\displaystyle=\sum_{a\in
S^{*}}\beta(a)\left(\frac{w(a)}{1-\alpha}-\frac{\nu^{*}(a)}{\beta(a)}\right)^{2}$
$\displaystyle=8n(1-\alpha)^{2}\beta(S^{*})\sum_{a}\frac{\left(\frac{\pi(a)}{1-\alpha}-\nu^{*}(a)\right)^{2}}{\beta(a)}$
$\displaystyle=8n(1-\alpha)^{2}\beta(S^{*})(T_{1}+T_{2})$

which implies that $S^{*}\neq[k]$. Then since $w(a)>0$ and $\nu^{*}\neq 0$,
both $T_{2}>0$ and $\beta(S^{*})>0$. Then we have

$\displaystyle T_{1}>8n(1-\alpha)^{2}\beta(S^{*})T_{1}$

which implies that

$\frac{1}{8n(1-\alpha)^{2}}>\beta(S^{*})>0.$

This means that for some $\epsilon\in(0,1)$,

$\displaystyle\beta(S^{*})=\frac{1-\epsilon}{8n(1-\alpha)^{2}}.$

Then combining this with the previous relationship,

$\displaystyle T_{1}$
$\displaystyle=8n(1-\alpha)^{2}\frac{1-\epsilon}{8n(1-\alpha)^{2}}(T_{1}+T_{2})$
$\displaystyle=\frac{1-\epsilon}{\epsilon}T_{2}.$

Next, we show that (26) has optimal value

$\displaystyle\pi(S^{*})+\sqrt{\frac{T_{2}}{8n}}\epsilon.$

From (28), for $a\in S^{*}$ we have that

$\displaystyle\nu^{*}(a)$
$\displaystyle=\frac{\pi(a)}{1-\alpha}-\beta(a)\sqrt{8(1-\alpha)^{2}n(T_{1}+T_{2})}.$

Then

$\displaystyle\sqrt{\frac{1}{8n}\sum_{a}\frac{\left(\frac{\pi(a)}{1-\alpha}-\nu^{*}(a)\right)^{2}}{\beta(a)}}+(1-\alpha)\sum_{a}|\nu^{*}(a)|=$
$\displaystyle\sqrt{\frac{1}{8n}(T_{1}+T_{2})}$
$\displaystyle+(1-\alpha)\sum_{a\in
S^{*}}\left(\frac{\pi(a)}{1-\alpha}-\beta(a)\sqrt{8(1-\alpha)^{2}n(T_{1}+T_{2})}\right)$
$\displaystyle=$
$\displaystyle\pi(S^{*})+\left(1-\beta(S^{*})8n(1-\alpha)^{2}\right)\sqrt{\frac{1}{8n}(T_{1}+T_{2})}$
$\displaystyle=$
$\displaystyle\pi(S^{*})+\epsilon\sqrt{\frac{1}{8n}(T_{1}+T_{2})}$
$\displaystyle=$ $\displaystyle\pi(S^{*})+\sqrt{\frac{\epsilon}{8n}T_{2}}$

Going back to our original problem, note that for $\epsilon\in[1/2,1]$ we
achieve the desired equivalence. Thus we focus on the case where
$\epsilon\in(0,1/2)$, which means that
$\frac{1}{16(1-\alpha)^{2}n}<\beta(S^{*})<\frac{1}{8(1-\alpha)^{2}n}$. Assume
WLOG that the actions are ordered according to likelihood ratios
$w(1)\leq\ldots\leq w(k)$. From (28) and (29), $S^{*}=\\{t+1,\ldots,k\\}$ for
some $t\in[k]$. Then applying (28), we have that for action $a=t+1$,

$\displaystyle\frac{w(t+1)^{2}}{(1-\alpha)^{2}}\geq\left(\frac{w(t+1)}{1-\alpha}-\frac{\nu^{*}(t+1)}{\beta(t+1)}\right)^{2}=8n(1-\alpha)^{2}(T_{1}+T_{2})=\frac{8n(1-\alpha)^{2}T_{2}}{\epsilon}$

Using the fact that $\frac{\pi(S^{*})}{\beta(S^{*})}\geq w(t+1)$, we have that

$\displaystyle\frac{\pi(S^{*})}{\beta(S^{*})}\geq(1-\alpha)^{2}\sqrt{\frac{8nT_{2}}{\epsilon}}$

Then since $\beta(S^{*})>\frac{1}{16(1-\alpha)^{2}n}$,

$\displaystyle\pi(S^{*})\geq\sqrt{\frac{T_{2}}{8n\epsilon}}\geq\sqrt{\frac{T_{2}}{4n}\epsilon}$

Then

$\displaystyle\sqrt{\frac{1}{8n}\sum_{a}\frac{\left(\frac{\pi(a)}{1-\alpha}-\nu^{*}(a)\right)^{2}}{\beta(a)}}+(1-\alpha)\sum_{a}|\nu^{*}(a)|$
$\displaystyle=\pi(S^{*})+\sqrt{\frac{\epsilon}{8n}T_{2}}\asymp\pi(S^{*})\asymp\pi(S^{*})+\sqrt{\frac{T_{2}}{n}}$

which achieves the target equivalence using the definition of $T_{2}$.

Finally, examining the optimality condition for $a\notin S^{*}$ (29) shows
that $\nu^{*}=0$ when

$\displaystyle
n\geq\frac{w_{max}^{2}}{8(1-\alpha)^{2}\sum_{a}\beta(a)w^{2}(a)},$

which means that the lower bound is

$\displaystyle\frac{cD}{1-\alpha}\sqrt{\frac{\sum_{a}\beta(a)w(a)^{2}}{n}}$

for some universal constant $c$, giving the theorem statement. ∎

##  Appendix E Experiments

###  E.1 Implementation

We describe each environment and their modeling methods below. For both
environments we use a maximum horizon $H=200$ and $\gamma=1$.

##### Cliffwalk.

Cliffwalk Sutton and Barto [2018] is a $4\times 12$ tabular environment where
the agent must travel from the start state (lower left corner) to the goal
state (lower right corner). The agent can take one of four cardinal actions
(left, down up, right). The bottom row of cells represents a cliff, and moving
into the cliff incurs a cost of 100. Each other state incurs a cost of 1, thus
incentivizing the agent to take the shortest path between start a goal. As the
original Cliffwalk is deterministic, we introduce a random transition towards
the cliff with probability $p=0.05$ in each state.

The target policy $\pi$ is learned via Q-learning, and the behavioral policy
is a mixture between the target and a uniform policy, that is for a constant
$\lambda\in[0,1]$, $\beta=\lambda\pi+(1-\lambda)\text{UNIF}$. For estimators
(DM, DR, WDR) involving models $\overline{F}$, half of the dataset is used to
construct the model and the other half is used for estimation. This process is
then repeated with the halves switched, and the resulting estimate is an
average of the two. First, an estimate of the MDP $\overline{M}$ is
constructed using empirical averages for the transitions and rewards.
Following Thomas et al. [2015], we assume that the model is given imperfect
information of the horizon as $H=201$. $\overline{F}$ is then computed
recursively via the relations from Section 5:

$\displaystyle\overline{F}_{s_{h},a_{h}}^{H+1-h}(t)=\mathbb{E}_{\overline{P},\overline{R}}\left[\overline{F}_{S_{h+1}}^{H-h}\left(\frac{t-R_{h}}{\gamma}\right)\right]$
$\displaystyle\overline{F}_{s_{h}}^{H+1-h}(t)=\mathbb{E}_{\pi}\left[\overline{F}_{s_{h},A_{h}}^{H+1-h}(t)\right]$

The Cliffwalk results shown in Figure 1 were averaged over 1000 repetitions.

##### Simglucose.

In Simglucose, the agent must control insulin bolus injections to a patient
with type 1 diabetes. The state is a continuous vector with the patient’s
blood glucose levels and the carbohydrate intake from the last meal. We
discretize the space of possible bolus injections into 6 actions. The agent
receives a reward according to whether the patient’s blood glucose levels are
within acceptable limits or not. If the patient’s blood glucose exceeds 180, a
condition called hyperglycemia, the agent receives a reward of -1. If it falls
under 70, called hypoglycemia, the agent receives a reward of -2. Otherwise,
the agent receives a reward of +1.

The target policy is the built-in controller for the patient, and the
behavioral policy is defined as the same way for Cliffwalk. As the state space
is continuous, the model was built by first discretizing the state space, then
using empirical averages to estimate the MDP $\overline{M}$ in the discretized
space, upon which $\overline{F}$ is calculated.

The Simglucose results shown in Figure 1 were averaged over 100 repetitions.

###  E.2 Mean vs Plug-in + CDF

As an additional experiment, we compare the plug-in mean estimate on the
estimated CDFs with direct mean estimation for existing
$\textsc{IS},\textsc{WIS}$ and DR Jiang and Li [2016] estimates in the
Cliffwalk environment. For WIS, at low sample sizes the direct mean estimates
can outperform the CDF estimates, but reach the same error as $n$ increases.
For the IS and DR estimates, however, the CDF version of the estimates perform
just as well, if not better, than direct mean estimation.

Figure 2: Normalized MSE for different $\lambda$ in the Cliffwalk environment.
Dashed lines show mean analogues of CDF estimates.