# Unbiased Test Error Estimation in the Poisson Means Problem via Coupled
Bootstrap Techniques

Natalia L. Oliveira1,2 Jing Lei1 Ryan J. Tibshirani1,2

(1Department of Statistics and Data Science, 2Machine Learning Department
Carnegie Mellon University)

###### Abstract

We propose a _coupled bootstrap_ estimator for the test error of an arbitrary
algorithm that estimates the mean in a Poisson sequence, often called the
Poisson means problem. The idea behind our method is to generate two
carefully-designed data vectors from the original data vector, by using
synthetic binomial noise. One such vector acts as the training sample and the
second acts as the test sample. To stabilize the test error estimate, we
average this over multiple draws of the synthetic noise. A key property of our
coupled bootstrap estimator is that it is unbiased for the test error in a
problem where the original mean has been shrunken by a small factor, driven by
the success probability $p$ in the binomial noise. Further, in the limit as
$p\to 0$, we show that the proposed estimator recovers a known unbiased
estimator for the test error, under no assumptions on the algorithm at hand
(in particular, no smoothness assumptions). Our methodology applies to two
central loss functions that can be used to a test error metric: Poisson
deviance and squared loss. Through a bias-variance decomposition, for each
loss function, we analyze the effects of the binomial success probability and
the number of bootstrap samples and on the accuracy of the estimator. We also
investigate our method empirically across a variety of settings, using
simulated as well as real data.

##  1 Introduction

We study the problem of estimating the test error of an algorithm in the
Poisson many means problem, also called the Poisson compound decision problem.
Given a data vector $Y=(Y_{1},\dots Y_{n})\in\mathbb{Z}_{+}^{n}$ (where we
write $\mathbb{Z}_{+}=\\{0,1,2,\dots\\}$ for the nonnegative integers)
distributed according to

$\text{$Y_{i}\sim\mathrm{Pois}(\mu_{i})$, \; independently, \; for
$i=1,\dots,n$},$ (1)

we seek to estimate the mean vector
$\mu=(\mu_{1},\dots,\mu_{n})\in\mathbb{R}^{n}_{+}$ (where we write
$\mathbb{R}_{+}=\\{x\in\mathbb{R}:x\geq 0\\}$ for the nonnegative reals). Let
$g:\mathbb{Z}_{+}^{n}\to\mathbb{R}_{+}^{n}$ be a measurable function that
estimates $\mu$ from the data $Y$, so that we can write $\hat{\mu}=g(Y)$. We
will often refer to $g$ as an algorithm, in the context of estimating $\mu$ in
the Poisson many means problem.

To evaluate the performance of $g$, we can use various metrics. One class of
metrics evaluate what we call _test error_ , based on a loss function
$L:\mathbb{Z}_{+}^{n}\times\mathbb{R}_{+}^{n}\to\mathbb{R}$,

$\text{$\mathrm{Err}(g)=\mathbb{E}[L({\tilde{Y}},g(Y))]$, \; where
${\tilde{Y}}$ is drawn from \eqref{eq:data_model}, independently of $Y$},$ (2)

which measures how well $g$ tracks an independent copy ${\tilde{Y}}$ of the
data. A second class of metrics evaluate what we call _risk_ , again based on
a loss function $L$,

$\mathrm{Risk}(g)=\mathbb{E}[L(\mu,g(Y))],$ (3)

which measures how well $g$ tracks the mean $\mu=\mathbb{E}[Y]$ of the data.
Admittedly, many authors use the terms “test error” and “risk”
interchangeably, but in this paper we are careful to use terminology that
distinguishes the two, for reasons that we will become apparent in the next
subsection.

###  1.1 Test error versus risk

In the classical normal means problem, where instead of (1) we observe
$Y_{i}\sim N(\mu_{i},\sigma^{2})$ independently, it is straightforward to show
that under a squared loss $L$, the test error (2) and risk (3) differ only by
the noise level $\sigma^{2}$. In the Poisson means problem, there is no direct
analogy for typical loss functions of interest, and the difference between (2)
and (3) will generally depend on $\mu$. This means that an estimator of one
metric (test error or risk) does not as easily translate into an estimator of
the other, since $\mu$ is of course unknown, and the primary estimand of
interest.

Thankfully, as we show here, when $L$ is a Bregman divergence the difference
between test error and risk does not depend on $g$. A Bregman divergence is a
loss function of the form $L(a,b)=D_{\phi}(a,b)$, where

$D_{\phi}(a,b)=\phi(a)-\phi(b)-\langle\nabla\phi(b),a-b\rangle,$ (4)

for a convex, differentiable function $\phi:\mathbb{R}^{n}\to\mathbb{R}$,
where here an subsequently we use $\langle u,v\rangle=u^{\mathsf{T}}v$ for
vectors $u,v$. In this case it is straightforward to see that

$\displaystyle\mathrm{Err}(g)-\mathrm{Risk}(g)$
$\displaystyle=\mathbb{E}[D_{\phi}({\tilde{Y}},g(Y))]-\mathbb{E}[D_{\phi}(\mu,g(Y))]$
$\displaystyle=\mathbb{E}[\phi({\tilde{Y}})]-\mathbb{E}[\phi(g(Y))]-\mathbb{E}[\langle\nabla\phi(g(Y)),{\tilde{Y}}-g(Y)\rangle]$
$\displaystyle\qquad-\phi(\mu)+\mathbb{E}[\phi(g(Y))]+\mathbb{E}[\langle\nabla\phi(g(Y)),\mu-g(Y)\rangle]$
$\displaystyle=\mathbb{E}[\phi(Y)]-\phi(\mu),$ (5)

where the cancellation of terms in the third line holds because
$Y,{\tilde{Y}}$ are i.i.d., and thus
$\mathbb{E}[\langle\nabla\phi(g(Y)),{\tilde{Y}}\rangle]=\langle\mathbb{E}[\nabla\phi(g(Y))],\mu\rangle$.
Observe that (5) is the gap in Jensen’s inequality. Therefore it is always
nonnegative, and $\mathrm{Err}(g)\geq\mathrm{Risk}(g)$.

In this paper, we will focus on estimating the test error (2) in the Poisson
means problem (1), for two special instances of a Bregman divergence: squared
loss and Poisson deviance, as will be discussed in the next subsection. Since
$\mathrm{Err}(g)-\mathrm{Risk}(g)=\mathbb{E}[\phi(Y)]-\phi(\mu)$ depends on
$\mu$, it will not be the case that we can automatically translate an
estimator of the test error of $g$ into an estimator of its risk. However, we
can still unbiasedly estimate the _difference_ in risk between two models $g$
and $h$, as discussed next.

#### Model comparisons.

The gap (5) does not depend on $g$. Thus for a comparison between two
algorithms $g$ and $h$, we always have (provided we use Bregman divergence to
define the test error and risk metrics):

$\mathrm{Err}(g)-\mathrm{Err}(h)=\mathrm{Risk}(g)-\mathrm{Risk}(h),$

To be clear, this means that if $\widehat{\mathrm{Err}}(g)$ is an unbiased
estimator of $\mathrm{Err}(g)$ for any $g$, just as we will produce in this
paper, then

$\text{$\widehat{\mathrm{Err}}(g)-\widehat{\mathrm{Err}}(h)$ \; is unbiased
for \; $\mathrm{Risk}(g)-\mathrm{Risk}(h)$, \; for any $g,h$}.$

As such, we can still use the tools developed in this paper to perform model
comparisons, or more broadly, model tuning (where $g_{s}$ is indexed by a
tuning parameter $s\in S$, and we select $s$ to minimize an unbiased estimate
of test error, or equivalently, risk).

#### Fixed-X Poisson regression.

A special case of our problem setting to keep in mind is _fixed-X_ Poisson
regression. Here we view $Y\in\mathbb{R}^{n}$ as a response vector and we have
an associated feature matrix $X\in\mathbb{R}^{n\times p}$. The algorithm $g$
typically performs a kind of Poisson regression of $Y$ on $X$. As long as we
consider $X$ to be fixed (nonrandom), we can still interpret this as a problem
of the form (1), with $\mu=\mu(X)$. In this setting, the test error metric (2)
translates to what is called fixed-X prediction error, where we evaluate
predictions at the same feature vectors (rows of $X$), but against new
responses (elements of ${\tilde{Y}}$).

While fixed-X analyses are more typical in classical statistics, the _random-
X_ perspective is great interest in modern prediction problems. Here the
feature vectors at which we make predictions are random, giving rise to
random-X prediction error as the metric of concern. Estimating random-X
prediction error is _not_ in general equivalent to estimating fixed-X
prediction error and the two can behave quite differently (see, e.g., Rosset
and Tibshirani (2020) for an extended discussion). The random-X perspective
eludes the framework of the current paper, but is an important topic for
future work.

###  1.2 Squared loss versus Poisson deviance

When $\phi(x)=\|x\|_{2}^{2}$, it is easy to check that

$D_{\phi}(a,b)=\|a-b\|_{2}^{2},$

which is the squared loss. Meanwhile, when $\phi(x)=2\sum_{i=1}^{n}x_{i}(\log
x_{i}-1)$, it follows that

$D_{\phi}(a,b)=2\sum_{i=1}^{n}\bigg{(}a_{i}\log\frac{a_{i}}{b_{i}}+b_{i}-a_{i}\bigg{)},$

which is known as _Poisson deviance_. We will take these to be the two loss
functions of primary interest in our work. Accordingly, we introduce the
notation for test error under squared loss and Poisson deviance:

$\displaystyle\mathrm{Err}^{\mathrm{sqr}}(g)$
$\displaystyle=\mathbb{E}\|{\tilde{Y}}-g(Y)\|_{2}^{2},$ (6)
$\displaystyle\mathrm{Err}^{\mathrm{dev}}(g)$
$\displaystyle=2\mathbb{E}\bigg{[}\sum_{i=1}^{n}\bigg{(}{\tilde{Y}}_{i}\log\frac{{\tilde{Y}}_{i}}{g_{i}(Y)}+g_{i}(Y)-{\tilde{Y}}_{i}\bigg{)}\bigg{]}.$
(7)

Squared loss is a standard choice in many estimation and prediction problems
and does not really need further motivation. Poisson deviance can be motivated
from different perspectives; one nice perspective is that, if we parametrize
$g_{i}(Y)=\exp(\theta_{i})$ for $i=1,\dots,n$, then fitting $g$ to minimize
Poisson deviance on the given data is equivalent to maximum likelihood in the
Poisson model,

$\mathop{\mathrm{minimize}}_{g}\;2\bigg{[}\sum_{i=1}^{n}\bigg{(}Y_{i}\log\frac{Y_{i}}{g_{i}(Y)}+g_{i}(Y)-Y_{i}\bigg{)}\bigg{]}\iff\mathop{\mathrm{minimize}}_{\theta}\;\sum_{i=1}^{n}\Big{(}-Y_{i}\theta_{i}+\exp(\theta_{i})\Big{)}.$

In the same vein, evaluating $g$ by Poisson deviance on ${\tilde{Y}}$ is
equivalent to evaluating $g$ by Poisson likelihood on an independent copy of
the training sample.

In our view, squared loss and Poisson deviance are each important loss
functions, and are each deserving of study. This is only strengthened by the
fact that they can have very different behaviors in certain problem settings.
As a simple example, suppose $n=1$, and we have two scenarios: in the first
${\tilde{Y}}=1$ and $g(Y)=2$, while in the second ${\tilde{Y}}=500$ and
$g(Y)=501$. The squared loss in each scenario is 1. However, the Poisson
deviance in first scenario is $\approx 0.307$, and in the second scenario it
is $\approx 0.001$. The difference here is driven by the fact that in the
Poisson model the variance scales with the mean. Hence according to Poisson
deviance (equivalent to Poisson likelihood), a prediction of 502 when the
predictand is 501 is not nearly as bad as a prediction of 2 when the
predictand is 1.

In Sections 5 and 6, we will present and discuss several examples that expose
differences in the behavior of squared loss and Poisson deviance in different
settings. That said, our primary focus is on estimating test error defined
with respect to these loss functions, and not on comparing them. A
comprehensive analysis of their differences is beyond the scope of the current
paper.

###  1.3 Hudson’s lemma

A fundamental result in this area is _Hudson’s lemma_ , due to Hudson (1978).
Hudson actually derived two identities, one each for continuous and discrete
exponential families. These can be viewed as extensions of Stein’s celebrated
identity (Stein, 1981) for the Gaussian family.111Stein’s work was actually
completed as a technical report in 1973, and was a motivation for Hudson’s
work, even though the publication dates of their papers do not reflect this.
According to Hudson, Stein already knew of the result in (8). For
concreteness, we state Hudson’s result for the Poisson case.

######  Lemma 1 (Hudson 1978).

Let $Y_{i}\sim\mathrm{Pois}(\mu_{i})$, independently, for $i=1,\dots,n$. Let
$g:\mathbb{Z}_{+}^{n}\to\mathbb{R}^{n}$ be such that
$\mathbb{E}|g_{i}(Y)|<\infty$, $i=1,\dots,n$. Then, denoting by
$e_{i}\in\mathbb{R}^{n}$ the vector whose $i^{\text{th}}$ entry is 1, with all
others 0,

$\mu_{i}\mathbb{E}[g_{i}(Y)]=\mathbb{E}[Y_{i}g_{i}(Y-e_{i})],\quad
i=1,\dots,n,$ (8)

where by convention we set $g_{i}(-1)=0$, $i=1,\dots,n$.

Compared to Stein’s identity, which requires that $g$ is weakly
differentiable, Hudson’s identity (8) holds without any smoothness assumptions
on $g$ (of course, even formulating precisely what smoothness would mean over
a discrete domain like $\mathbb{Z}_{n}^{+}$ would be tricky, but the lack of
assumptions needed for Lemma 1 are remarkable nonetheless). Hudson’s main
interest was in developing inadmissibility results for estimators of the
location parameter in an exponential family distribution. The identities he
established were used as tools in his analysis, which parallels Stein’s use of
his own identity in Stein (1981).

Moreover, analogous to what can be done with Stein’s lemma, Hudson’s lemma can
be used to derived unbiased estimators for various risk metrics in exponential
families. An important contribution in this area is Eldar (2009), and further
contributions (along with a comprehensive summary of available tools and
results from the literature) are given in Deledalle (2017).

###  1.4 Unbiased estimation

Our focus in this paper is slightly unique, since we consider test error (2)
as the primary target and not risk (3), as considered by Eldar (2009);
Deledalle (2017), and most other authors in the literature. Nonetheless, the
estimators developed by these authors have natural analogues for test error.
In fact, the story is for test error is simpler, and an unbiased estimator can
be obtained for any Bregman divergence loss function.

To see this, we first recall a general decomposition of test error for Bregman
divergence losses known as _Efron’s optimism theorem_ , due to Efron (1975,
1986, 2004): this shows that for any Bregman divergence $D_{\phi}$ in (4) and
any algorithm $g$, this difference in test error and training error satisfies

$\displaystyle\mathbb{E}[D_{\phi}({\tilde{Y}},g(Y))]-\mathbb{E}[D_{\phi}(Y,g(Y))]$
$\displaystyle=\mathbb{E}[\phi({\tilde{Y}})]-\mathbb{E}[\phi(g(Y))]-\mathbb{E}[\langle\nabla\phi(g(Y)),{\tilde{Y}}-g(Y)\rangle]$
$\displaystyle\qquad-\mathbb{E}[\phi(Y)]+\mathbb{E}[\phi(g(Y))]+\mathbb{E}[\langle\nabla\phi(g(Y)),Y-g(Y)\rangle]$
$\displaystyle=\mathbb{E}[\langle\nabla\phi(g(Y)),Y\rangle]-\mathbb{E}[\langle\nabla\phi(g(Y)),\mu\rangle].$
(9)

The second line follows from the fact that
$\mathbb{E}[\phi({\tilde{Y}})]=\mathbb{E}[\phi(Y)]$.222This exposes the reason
why the analogous decomposition for risk can be more complex: when we replace
${\tilde{Y}}$ with $\mu$ in the calculation that led to (9), we are left with
an extra term $\phi(\mu)-\mathbb{E}[\phi(Y)]$ that does not cancel and must be
estimated. Simply rewriting the above, we see that if we are able to construct
an unbiased estimator $V(g)$ of
$\mathbb{E}[\langle\nabla\phi(g(Y)),\mu\rangle]$ then

$D_{\phi}(Y,g(Y))+\langle\nabla\phi(g(Y)),Y\rangle-V(g)$

will be an unbiased estimator for $\mathbb{E}[D_{\phi}({\tilde{Y}},g(Y))]$. In
the Poisson case, Hudson’s identity (8) precisely gives the unbiased estimator
$V(g)$ that we require, which leads to the following result.

######  Proposition 1.

Let $Y_{i}\sim\mathrm{Pois}(\mu_{i})$, independently, for $i=1,\dots,n$. Let
$g:\mathbb{Z}_{+}^{n}\to\mathbb{R}^{n}$ be any algorithm and $D_{\phi}$ be any
Bregman divergence loss function (indexed by a convex, differentiable function
$\phi:\mathbb{R}^{n}\to\mathbb{R}$) such that $\mathbb{E}|\phi(g(Y))|<\infty$
and $\mathbb{E}|\nabla_{i}\phi(g(Y))|<\infty$, $i=1,\dots,n$. Then

$\mathrm{UE}(g)=D_{\phi}(Y,g(Y))+\langle\nabla\phi(g(Y)),Y\rangle-\langle\nabla\phi(g_{-}(Y)),Y\rangle$
(10)

is unbiased for $\mathrm{Err}(g)=\mathbb{E}[D_{\phi}({\tilde{Y}},g_{i}(Y))]$,
where we abbreviate $g_{-}(Y)=(g_{1}(Y-e_{1}),\dots,g_{n}(Y-e_{n}))$, and as
usual, ${\tilde{Y}}$ denotes an independent copy of $Y$.

As a consequence, we have the following unbiased estimators for squared loss
and Poisson deviance:

$\displaystyle\mathrm{UE}^{\mathrm{sqr}}(Y)$
$\displaystyle=\|Y\|_{2}^{2}+\|g(Y)\|_{2}^{2}-2\langle g_{-}(Y),Y\rangle,$
(11) $\displaystyle\mathrm{UE}^{\mathrm{dev}}(Y)$
$\displaystyle=2\sum_{i=1}^{n}\Big{(}Y_{i}\log Y_{i}-Y_{i}\log
g_{i}(Y-e_{i})+g_{i}(Y)-Y_{i}\Big{)}.$ (12)

These are altogether highly similar to the unbiased risk estimators in Eldar
(2009); Deledalle (2017), and to be clear, we do not consider (11), (12) to be
major (or even original) contributions of our work. That said, we have not yet
seen the general unbiased estimator for Bregman divergence (10) noted in the
literature, thus we believe it may be useful to record it (along with the
observation that estimation of test error can be easier than estimation of
risk).

The estimators in (11), (12) have a clear strength: they are unbiased _for any
algorithm $g$_. This is a strong property; recall that by comparison, in the
Gaussian model, the analogous estimator is Stein’s unbiased risk estimator
(SURE), which requires $g$ to be weakly differentiable. The estimators in
(11), (12) also have a clear downside: they require the algorithm $g$ to be
run $n+1$ times, once to obtain the original fit $g(Y)$, and then $n$ more
times to obtain $g_{-}(Y)$, which recall has entries $g_{i}(Y-e_{i})$,
$i=1,\dots,n$. Thus we can liken (11), (12) to leave-one-out cross-validation,
in terms of computational cost.

This draws a clear line of motivation to the main contribution of our paper:
in what follows, we develop an unbiased estimator of test error, for any
Bregman divergence loss function $D_{\phi}$ and any algorithm $g$, using a
carefully-crafted parametric bootstrap scheme. The computational cost (number
of runs of $g$) here is tied to a user-controlled parameter $B$, the number of
bootstrap samples. In general, increasing $B$ decreases the variance of the
estimator, but any choice of $B\geq 1$ yields an estimator that is unbiased
for the test error in a _mean-shrunken_ Poisson problem, with mean $(1-p)\mu$,
where $p>0$ is another user-controlled parameter.

###  1.5 Summary of contributions

The following gives a summary of our main contributions and an outline for
this paper.

  * •

In Section 2, we introduce the coupled bootstrap (CB) estimator, and prove
that it is unbiased for the test error in a mean-shrunken Poisson problem.

  * •

In Section 3, we analyze the behavior of the CB estimator as $B\to\infty$ and
$p\to 0$, and prove that the limiting CB estimator is exactly the unbiased
estimator (10) from Hudson’s lemma.

  * •

In Section 4, we study the bias and variance of the CB estimator and quantify
how they depend on $B,p$ and other problem parameters.

  * •

In Section 5, we compare the CB and the unbiased estimator on various
simulated data sets, and find that the performance of the CB estimator is
favorable, especially when the algorithm $g$ is unstable.

  * •

In Section 6, we examine the use of the CB estimator for model
tuning—selecting from a family $g_{s}$, $s\in S$ of algorithms—in two
applications: image denoising and density estimation. We find that using
Poisson deviance (to define the test error metric) consistently delivers more
regularized models than using squared loss.

  * •

In Section 7, we conclude with a brief discussion and ideas for future work.

###  1.6 Related work

Estimating risk and test error is of central importance in statistics and
machine learning. In the random-X prediction setting (which recall does not
fit in the framework of our work) the most ubiquitous estimator is arguably
cross-validation, which itself carries a long line of literature. We do not
describe this literature here, but highlight Bates et al. (2021) as a nice
recent paper that carefully reexamines this classic estimator, and also
provides a nice overview of literature on cross-validation.

In the fixed-X prediction setting—or in general, parametric many means
problems—there has also been a long history of work in statistics, with Akaike
(1973); Mallows (1973); Efron (1975); Stein (1981); Efron (1986) marking early
important contributions. This has been particularly well-studied in the
Gaussian means problem, and in this area, we draw attention to Breiman (1992);
Ye (1998); Efron (2004), and particularly to Oliveira et al. (2021), as
motivation for our current work in the Poisson means problem. These papers use
auxiliary noise—they inject synthetic (Gaussian) noise into the data at
hand—in order to estimate the risk or fixed-X prediction error of an arbitrary
algorithm $g$. Our previous work, Oliveira et al. (2021), proposes a coupled
bootstrap (CB) scheme for doing so that has a simple, intuitive target of
estimation for any auxiliary noise level. In particular, for any auxiliary
noise level $\alpha>0$ (a user-controlled parameter), the CB method produces
an unbiased estimator for the risk in a Gaussian means problem that has an
inflated noise variance $(1+\alpha)\sigma^{2}$ (where $\sigma^{2}$ denotes the
original noise level). The current paper builds off this idea, and develops a
coupled boostrap scheme in the Poisson model that enjoys analogous properties.

Relative to the Gaussian case, risk and test error estimation in the Poisson
means model has been less well-studied. However, there has still certainly
been important and influential work in the area. This includes Hudson (1978),
as already described in the introduction, and also Shen et al. (2004); Eldar
(2009); Deledalle (2017). Meanwhile, the literature on _mean estimation_ —the
role played by what we are calling the algorithm $g$—in the Poisson many means
problem is vast. Quite a lot of work on this topic has been done in the signal
processing community, where it is often called Poisson denoising; see, e.g.,
Harmany et al. (2009); Luisier et al. (2010); Raginsky et al. (2010); Harmany
et al. (2012); Salmon et al. (2014); Cao and Xie (2016). Therefore we believe
that the techniques we develop for estimating test error in the Poisson means
model should have widespread practical applications, in signal processing and
elsewhere.

Lastly, we mention a concurrent, related line of work on auxiliary
randomization approaches that allow for rigorous post-selection inference in
parametric many means models, including the Poisson means model. We highlight
Leiner et al. (2021); Neufeld et al. (2023) as two nice recent papers in the
area. In particular, a core piece of the auxiliary randomization procedure in
our work was directly inspired by the former paper, and their use of binomial
auxiliary noise in the Poisson model.

##  2 Coupled boostrap estimator

In this section, we introduce the CB estimator in the Poisson means model, and
investigate some of its basic properties.

###  2.1 Proposed estimator

The following is a simple but key “three-point” formula for expected Bregman
divergence loss from Oliveira et al. (2021) that will drive our main proposal
in this paper.

######  Proposition 2.

Let $U,V,W\in\mathbb{R}^{n}$ be independent random vectors. For any $g$, and
Bregman divergence $D_{\phi}$,

$\mathbb{E}[D_{\phi}(V,g(U))]-\mathbb{E}[D_{\phi}(W,g(U))]=\mathbb{E}[\phi(V)]-\mathbb{E}[\phi(W)]+\langle\mathbb{E}[\nabla\phi(U)],\mathbb{E}[W]-\mathbb{E}[V]\rangle,$
(13)

assuming all expectations exist and are finite. In particular, if $U,V$ are
i.i.d. and $\mathbb{E}[U]=\mathbb{E}[W]$, then

$\mathbb{E}[D_{\phi}(V,g(U))]=\mathbb{E}[D_{\phi}(W,g(U))]+\mathbb{E}[\phi(U)]-\mathbb{E}[\phi(W)].$
(14)

###### Proof.

The first statement (13) follows from the definition of Bregman divergence
(4), and the independence of $U,V,W$. The second result (14) follows from the
first, by noting that if $U,V$ are i.i.d. and $\mathbb{E}[U]=\mathbb{E}[W]$
then $\mathbb{E}[V]=\mathbb{E}[W]$, thus the last term on the right-hand side
in (13) is zero, and the first term is $\mathbb{E}[\phi(U)]$. ∎

While simple to state and prove, the results in Proposition 2 are useful
observations. To map them onto to the problem of estimating of test error in
the Poisson model, consider the following. Given a Poisson data vector $Y$
from (1), suppose that we can generate a pair of vectors
$(U,W)=(Y^{*},Y^{\dagger})$ that are independent of each other and have the
same mean. Then (14) says that

$\text{$D_{\phi}(Y^{\dagger},g(Y^{*}))+\phi(Y^{*})-\phi(Y^{\dagger})$ \; is
unbiased for \; $\mathbb{E}[D_{\phi}({\tilde{Y}}^{*},g(Y^{*}))]$},$ (15)

where ${\tilde{Y}}^{*}$ is an independent copy of $Y^{*}$. In other words, the
above constructs an unbiased esitmator for the test error in a problem in
which the original data vector was $Y^{*}$, rather than $Y$. Thus if _$Y^{*}$
was “close” in distribution to $Y$_, then this estimator would be meaningful.
(Ideally, we would like $Y^{*}$ to be be identical in distribution to $Y$, but
that will not be generically possible without knowledge of $\mu$.)

What remains is a precise scheme in the Poisson setting to generate the pair
$(Y^{*},Y^{\dagger})$ from $Y$ such that $Y^{*},Y^{\dagger}$ are independent,
$\mathbb{E}[Y^{*}]=\mathbb{E}[Y^{\dagger}]$, and $Y^{*},Y$ are “close” in
distribution. The next lemma does the trick and fulfills these three
properties precisely. It was brought to our attention by Leiner et al. (2021)
who used it in a distinct but generally related post-selection inference
context. For completeness, we provide a proof in Appendix A.1. Here and
henceforth, we use the following abbreviations: we write
$Y\sim\mathrm{Pois}(\mu)$ to mean that we draw
$Y_{i}\sim\mathrm{Pois}(\mu_{i})$, independently, for $i=1,\dots,n$, and
similarly $Z\sim\mathrm{Binom}(N,p)$ to mean that we draw
$Z_{i}\sim\mathrm{Binom}(N_{i},p)$, independently, for $i=1,\dots,n$.

######  Lemma 2.

Given $Y\sim\mathrm{Pois}(\mu)$, fix any $0<p<1$ and let
$\omega\,|\,Y\sim\mathrm{Binom}(Y,p)$. Then, defining $Y^{*}=Y-\omega$ and
$Y^{\dagger}=(1-p)/p\cdot\omega$, it holds that:

  1. (i)

$Y^{*},Y^{\dagger}$ are independent;

  2. (ii)

$\mathbb{E}[Y^{*}]=\mathbb{E}[Y^{\dagger}]$; and

  3. (iii)

$Y^{*}\sim\mathrm{Pois}((1-p)\mu)$.

Lemma 2, combined with the observation in (15), forms the basis for the CB
test error estimator in the Poisson many means problem. To stabilize the
estimator, we can simply repeat the draws of binomial noise from Lemma 2 over
independent repetitions $b=1,\dots,B$. To be concrete, this leads to the
following method: we first generate samples according to

$\begin{gathered}\text{$\omega^{b}\,|\,Y\sim\mathrm{Binom}(Y,p)$, \;
independently, \; for $b=1,\dots,B$},\\\
Y^{*b}=Y-\omega^{b},\quad,Y^{\dagger{b}}=\frac{1-p}{p}\omega^{b},\quad\text{for
$b=1,\dots,B$},\end{gathered}$ (16)

for an arbitrary binomial success probability $0<p<1$, and a number of
bootstrap draws $B\geq 1$; then we define the _coupled bootstrap_ (CB)
estimator, for test error under Bregman divergence loss $D_{\phi}$, by:

$\mathrm{CB}_{p}(g)=\frac{1}{B}\sum_{b=1}^{B}\Big{(}D_{\phi}(Y^{\dagger{b}},g(Y^{*b}))+\phi(Y^{*b})-\phi(Y^{\dagger{b}})\Big{)}.$
(17)

We can view each $Y^{*b}$ as a synthetic training set for $g$, and each
$Y^{\dagger{b}}$ as a synthetic test set. The correction term
$\phi(Y^{*b})-\phi(Y^{\dagger{b}})$ accounts for the fact that
$Y^{*b},Y^{\dagger{b}}$ do not have the same distribution (though recall they
do have the same mean, by construction).

For the two loss functions of primary interest, squared loss and Poisson
deviance, the CB estimator in (17) becomes, respectively:

$\displaystyle\mathrm{CB}_{p}^{\mathrm{sqr}}$
$\displaystyle=\frac{1}{B}\sum_{b=1}^{B}\Big{(}\|Y^{\dagger{b}}-g(Y^{*b})\|_{2}^{2}+\|Y^{*b}\|^{2}_{2}-\|Y^{\dagger{b}}\|_{2}^{2})\Big{)},$
(18) $\displaystyle\mathrm{CB}_{p}^{\mathrm{dev}}$
$\displaystyle=\frac{2}{B}\sum_{b=1}^{B}\sum_{i=1}^{n}\Big{(}Y^{*b}_{i}\log
Y^{*b}_{i}-Y^{\dagger{b}}_{i}\log
g_{i}(Y^{*b})+g_{i}(Y^{*b})-Y^{\dagger{b}}_{i}\Big{)}.$ (19)

#### Interlude: special care with deviance estimators.

We take a brief but practically important detour to note that special care
must be taken with test error estimators with respect to Poisson deviance
loss. In this case, each of the unbiased (12) and the coupled bootstrap (19)
estimators can diverge if the coordinate functions of $g$ can output zero. For
the unbiased estimator this occurs when $Y_{i}\not=0$ and $g_{i}(Y-e_{i})=0$;
for the coupled bootstrap estimator this occurs when
$Y^{\dagger{b}}_{i}\not=0$ and $g_{i}(Y^{*b})=0$. As a safety mechanism, we
can simply pad the output of $g$ so that zero is never in the range of its
coordinate functions: say, we can define a modified algorithm

$\tilde{g}_{i}(y)=g_{i}(y)1\\{g_{i}(y)\not=0\\}+c1\\{g_{i}(y)=0\\},\quad
i=1,\dots,n,$

for a small constant $c>0$. This is reasonable because even the population
Poisson deviance (7) can itself diverge when the coordinate functions of $g$
can output zero. With a modified rule like the one above, we may ask how
frequently the padding is actually in effect in the computation of the
estimators (12) and (19). We study this in Appendix A.2 and show that, in a
sense, it is typically in effect less frequently in the CB estimator (19) than
in the unbiased one (12).

###  2.2 Unbiasedness for mean-shrunken target

The next result is immediate from Lemma 2 and (15).

######  Corollary 1.

Let $Y\sim\mathrm{Pois}(\mu)$. Let $g:\mathbb{Z}_{+}^{n}\to\mathbb{R}^{n}$ be
any algorithm, let $D_{\phi}$ be any Bregman divergence loss function, and let
$0<p<1$ and $B\geq 1$ be arbitrary. Then the CB estimator $\mathrm{CB}_{p}(g)$
in (17) is unbiased for $\mathrm{Err}_{p}(g)$ (assuming all terms in (17) have
finite expectations), where $\mathrm{Err}_{p}(g)$ is the test error of $g$
with respect to a mean-shrunken Poisson problem:

$\mathrm{Err}_{p}(g)=\mathbb{E}[D_{\phi}({\tilde{Y}}_{p},g(Y_{p}))]$, where
$Y_{p},{\tilde{Y}}_{p}\sim\mathrm{Pois}((1-p)\mu)$, and
$Y_{p},{\tilde{Y}}_{p}$ are independent. (20)

The strength of Corollary 1 rests on the fact that the estimand in (17) of
$\mathrm{CB}_{p}(g)$ for any choice of $p>0$ is highly intuitive: it is
$\mathrm{Err}_{p}(g)$ in (20), which is the test error that we would encounter
in a slightly harder version of our original problem, where the mean $\mu$ has
been replaced by $(1-p)\mu$.

Why is this important? It means that we do _not_ have to send $p\to 0$ in
order to be able to interpret the estimand of the CB estimator, and thus
justify its use. Any nonzero (noninfinitesimal) $p$ will still result in a
target that has a clear, intuitive meaning. This is good news for the CB
estimator, because when $p$ is away from zero, we can generally choose a
reasonably small number of bootstrap draws $B$ in order to stabilize the
variance of the estimator, which presents a computational advantage over the
unbiased estimator in (10). We will learn more about the behavior of the CB
estimator, as we vary $p$ and $B$, in Sections 4 and 5 (where we formally
analyze the bias and variance, and carry out empirical comparisons,
respectively).

###  2.3 Smoothness of mean-shrunken target

Now that we have shown that $\mathrm{CB}_{p}(g)$ is unbiased for
$\mathrm{Err}_{p}(g)$, it is natural to ask whether $\mathrm{Err}_{p}(g)$ will
be close to $\mathrm{Err}(g)$ for small $p$. Our next result gives a partial
answer by proving that if $g$ satisfies some mild moment conditions, then the
map $p\mapsto\mathrm{Err}_{p}(g)$ will be continuous (and in fact, it can be
continuously differentiable, depending on the number of moments assumed) in an
interval containing $p=0$. Later on, in Section 4, we will derive results that
give a more quantitative sense of how close $\mathrm{Err}_{p}(g)$ can be to
$\mathrm{Err}(g)$.

######  Proposition 3.

For $0\leq p<1$, let $\mathrm{Err}_{p}(g)$ be as defined in (20). If for some
integer $k\geq 0$,

$\mathbb{E}\big{[}D_{\phi}({\tilde{Y}},g(Y))\langle{\tilde{Y}}+Y,1_{n}\rangle^{m}\big{]}<\infty,\quad
m=0,\dots,k,$

where in the above above $Y,{\tilde{Y}}\sim\mathrm{Pois}(\mu)$ are
independent, and $1_{n}\in\mathbb{R}^{n}$ denotes the vector of all 1s, then
the map $p\mapsto\mathrm{Err}_{p}(g)$ has $k$ continuous derivatives on
$[0,1)$.

The proof of this result is not difficult but a bit technical and deferred to
Appendix A.3. We remark that when $k=0$, the assumption in Proposition 3 is
simply $\mathrm{Err}(g)=\mathbb{E}[D_{\phi}({\tilde{Y}},g(Y))]<\infty$ (i.e.,
the original test error is finite), which is extremely weak, and even in this
case we get that $\mathrm{Err}_{p}(g)\to\mathrm{Err}(g)$ as $p\to 0$.

##  3 Noiseless limit

In this section, we consider the _infinite-bootstrap_ version of the CB
estimator, $\mathrm{CB}_{p}^{\infty}(g)=\lim_{B\to\infty}\mathrm{CB}_{p}(g)$.
By the law of large numbers, this estimator is equivalent to taking the
expectation over the binomial noise,

$\mathrm{CB}_{p}^{\infty}(g)=\mathbb{E}[\mathrm{CB}_{p}(g)\,|\,Y]=\mathbb{E}\bigg{[}D_{\phi}\bigg{(}\frac{1-p}{p}\omega,g(Y-\omega)\bigg{)}+\phi(Y-\omega)-\phi\bigg{(}\frac{1-p}{p}\omega\bigg{)}\bigg{]},$
(21)

where $\omega\,|\,Y\sim\mathrm{Binom}(Y,p)$, as in (16).

The next result considers the _noiseless limit_ of the infinite-bootstrap
version of the CB estimator (21), where $p\to 0$. Its proof is deferred until
Appendix B.

######  Theorem 1.

Let $Y\sim\mathrm{Pois}(\mu)$. Let $g:\mathbb{Z}_{+}^{n}\to\mathbb{R}^{n}$ be
any algorithm, let $D_{\phi}$ be any Bregman divergence loss function, and
assume that $|\phi_{i}(g(Y))|<\infty$, $|\nabla_{i}\phi(g(Y))|<\infty$, and
$|\nabla_{i}\phi(g(Y-e_{i}))|<\infty$ almost surely, for each $i=1,\dots,n$.
Then

$\lim_{p\to 0}\mathrm{CB}_{p}^{\infty}(g)=\mathrm{UE}(g),\quad\text{almost
surely},$ (22)

where $\mathrm{UE}(g)$ is the unbiased estimator defined in (10). Thus as a
consequence, the noiseless limit of $\mathrm{CB}_{p}^{\infty}(g)$ is unbiased
for $\mathrm{Err}(g)$.

That the limiting CB estimator (17) recovers the unbiased estimator (10) based
on Hudson’s lemma, as $B\to\infty$ and $p\to 0$, is certainly an encouraging
property for the former. We recall that in the Gaussian many means problem,
the analogous result was derived in Oliveira et al. (2021): there, the CB
estimator recovers the unbiased estimator based on Stein’s lemma, in the
noiseless limit. However, this Gaussian result requires $g$ to be weakly
differentiable (which is the condition required for Stein’s unbiased estimator
to be valid in the first place). In the current Poisson many means problem,
note that Theorem 1 requires no such restrictions on $g$ (and indeed, recall,
the unbiased estimator does not either, from Proposition 1).

Figure 1 illustrates this difference via a simple simulation; see the figure
caption for details.

Figure 1: Density of the CB and unbiased estimators of test error in Gaussian
and Poisson settings, where the CB estimators effectively take large $B$ and a
small amount of auxiliary noise. This is based on a simulation with $n=30$,
where we generate Gaussian or Poisson data with constant mean, and consider
two estimators: soft-thresholding and hard-thresholding (the latter violating
weak differentiability). The black vertical line in each panel marks the true
test error. For small auxiliary noise in the Gaussian hard-thresholding
setting, it is clear that CB and the unbiased estimator are separated (and the
latter is far from unbiased).

##  4 Bias and variance

In this section, we study a bias-variance decomposition of the estimator
$\mathrm{CB}_{p}(g)$ in (17), when targeting the _original_ error
$\mathrm{Err}(g)$. We will consider an arbitrary Bregman divergence loss (used
to define $\mathrm{Err}(g)$), and use the decomposition

$\mathbb{E}[\mathrm{CB}_{p}(g)-\mathrm{Err}_{p}(g)]^{2}=\underbrace{\big{[}\mathrm{Err}_{p}(g)-\mathrm{Err}(g)\big{]}^{2}}_{\mathrm{Bias}^{2}(\mathrm{CB}_{p}(g))}\,+\,\underbrace{\mathbb{E}\big{[}\mathrm{Var}(\mathrm{CB}_{p}(g)\,|\,Y)\big{]}}_{\mathrm{RVar}(\mathrm{CB}_{p}(g))}\,+\,\underbrace{\mathrm{Var}\big{(}\mathbb{E}[\mathrm{CB}_{p}(g)\,|\,Y]\big{)}}_{\mathrm{IVar}(\mathrm{CB}_{p}(g))}.$
(23)

This is the usual bias-variance decomposition of squared error loss, where we
have used $\mathbb{E}[\mathrm{CB}_{p}(g)]=\mathrm{Err}_{p}(g)$ in the bias
term, and we have further expanded the usual variance term (using the law of
total variance) into two components which we call the _reducible_ and
_irreducible_ variance, respectively, as in Oliveira et al. (2021). We note
that as the number of bootstrap draws $B$ grows, the reducible variance
shrinks, but the irreducible variance does not; the latter does not depend on
$B$ at all, and in fact, it can be viewed as the variance of the infinite-
bootstrap version of the estimator,
$\mathrm{CB}_{p}^{\infty}(g)=\mathbb{E}[\mathrm{CB}_{p}(g)\,|\,Y]$.

In what follows, we will analyze each of the three terms in (23) to understand
their behavior as functions of $p$ and $B$, with a focus on small $p$ and
large $B$. As usual, we assume throughout that $Y\sim\mathrm{Pois}(\mu)$,
where $Y_{p}\sim\mathrm{Pois}((1-p)\mu)$ for $p\geq 0$, and we denote by
${\tilde{Y}},{\tilde{Y}}_{p}$ independent copies of $Y,Y_{p}$, respectively.
Lastly, $D_{\phi}$ represents an arbitrary Bregman divergence loss.

###  4.1 Bias

First we give an exact expression for the bias,
$\mathrm{Bias}(\mathrm{CB}_{p}(g))=\mathrm{Err}_{p}(g)-\mathrm{Err}(g)$, and
an upper bound on its magnitude for small $p$, under an assumption of monotone
variance. The proof is given in Appendix C.1.

######  Proposition 4.

Assume that $\mathbb{E}[D_{\phi}({\tilde{Y}}_{p},g_{p}(Y))\langle
Y_{p}+{\tilde{Y}}_{p},1_{n}\rangle]<\infty$. Then for all $p\in[0,1)$,

$\mathrm{Err}_{p}(g)-\mathrm{Err}(g)=-\sqrt{2\sum_{i=1}^{n}\mu_{i}}\int_{0}^{p}\frac{1}{\sqrt{1-t}}\mathrm{Cor}\Big{(}D_{\phi}({\tilde{Y}}_{t},g(Y_{t})),\langle{\tilde{Y}}_{t}+Y_{t},1_{n}\rangle\Big{)}\sqrt{\mathrm{Var}\big{[}D_{\phi}({\tilde{Y}}_{t},g(Y_{t}))\big{]}}\,dt.$
(24)

Further, if $\mathrm{Var}[D_{\phi}({\tilde{Y}}_{p},g(Y_{p}))]$ is decreasing
in $p$ on $[0,1/2]$, then for any $p$ in this range,

$|\mathrm{Err}_{p}(g)-\mathrm{Err}(g)|\leq\frac{5p}{3}\sqrt{\mathrm{Var}\big{[}D_{\phi}({\tilde{Y}},g(Y))\big{]}\sum_{i=1}^{n}\mu_{i}}$
(25)

We remark that the assumption of decreasing variance of the loss is fairly
natural (because the variance of each component of $Y_{p}$ decreases
monotonically to $0$ as $p$ increases to $1$). We can also drop this
condition, and replace the variance term in the bound (25) by
$\sup_{t\in(0,p)}\mathrm{Var}[D_{\phi}({\tilde{Y}}_{t},g(Y_{t}))]$.

###  4.2 Reducible variance

Next we bound the reducible variance, $\mathrm{RVar}(\mathrm{CB}_{p}(g))$. We
focus on the dependence on $p$ and $B$, for small $p$ and large $B$. The
notation $O(\cdot)$ is to be interpreted in this regime (small $p$, large
$B$), and hides factors that may depend on the mean $\mu$, which may in turn
depend on the dimensionality $n$. The proof is given in Appendix C.2.

######  Proposition 5.

Assume the variables $f(Y_{p})$, $f^{2}(Y_{p})$, $f(Y_{p})\langle
Y_{p},1_{n}\rangle$, $f^{2}(Y_{p})\langle Y_{p},1_{n}\rangle$ all have finite
$L^{1}$ norm, uniformly bounded over all functions $f\in\mathcal{F}$ and all
$p\in[0,q)$, for some $q>0$, where

$\mathcal{F}=\big{\\{}y\mapsto
D_{\phi}(y,g(y))\big{\\}}\cup\big{\\{}y\mapsto\nabla_{i}\phi(g(y)):i=1,\dots,n\big{\\}}.$

Then for all $p\in[0,q)$,

$\mathrm{RVar}(\mathrm{CB}_{p}(g))\leq\frac{2}{B}\mathrm{Var}\Big{[}D_{\phi}(Y,g(Y)))+\langle
Y,\nabla\phi(g(Y))\rangle\Big{]}+\frac{2}{Bp}\sum_{i=1}^{n}\mu_{i}\mathbb{E}\big{[}\nabla_{i}\phi(g(Y_{p}))^{2}\big{]}+{}\\\
\frac{2}{B}\sum_{i=1}^{n}\mu_{i}^{2}\mathrm{Var}\big{[}\nabla_{i}\phi(g(Y_{p}))\big{]}+O\bigg{(}\frac{p}{B}\bigg{)}.$
(26)

A simple simulation, whose results are presented in Figure 2, shows that the
reducible variance bound (26) appears to have the right dependence on $\mu$
and $B$. See the figure caption for details.

(a) Squared loss

(b) Deviance loss

Figure 2: Comparison of the true reducible variance (approximated by Monte
Carlo) and the bound given in (26), for squared and deviance loss, in a
simulation with $n=100$ and $p=0.1$. The data vector $Y$ has Poisson entries,
and $\mu$ denotes the common mean of each component; we use a simple linear
shrinkage estimator $g$. We can see that the behavior for varying $B,\mu$
looks qualitatively similar across the true reducible variance heatmap and the
bound heatmap, for each loss function.

###  4.3 Irreducible variance

Last we analyze the irreducible variance, $\mathrm{IVar}(\mathrm{CB}_{p}(g))$.
The proof is given in Appendix C.3.

######  Proposition 6.

Assume that $\mathbb{E}[D^{2}_{\phi}(Y,g(Y))]<\infty$ and
$\mathbb{E}[\langle\nabla(g(Y)),Y\rangle^{2}]<\infty$. Define $\Phi_{g}$ to
have component functions

$\Phi_{g,i}(y)=\sup_{0\leq z\leq y}\,\big{|}\nabla_{i}\phi(g(y))\big{|},\quad
i=1,\dots,n.$

(Here when we write $0\leq z\leq y$, all inequalities are to be interpreted
componentwise.) Then,

$\lim_{p\to 0}\mathrm{IVar}(\mathrm{CB}_{p}(g))\leq
2\mathrm{Var}\Big{[}D_{\phi}(Y,g(Y))+\langle\nabla\phi(g(Y)),Y\rangle\Big{]}+2\mathbb{E}\big{[}\langle\Phi_{g}(Y),Y\rangle^{2}\big{]}.$
(27)

###  4.4 Discussion of bias and variance results

We discuss interpretation of the results above. The bias bound (25) decreases
linearly with $p$, which suggests that we should take $p$ to be as small as
possible in order to decrease the bias. The irreducible variance bound (27)
provides no resistance to this idea, as it has a stable noiseless limit, as we
send $p\to 0$. The behavior of the reducible variance bound (26), however, is
more intricate. The second term on the right-hand side in (26) diverges as
$p\to 0$, but this can be offset by sending $B\to\infty$.

How large do we need to take $B$? Altogether, there are really only two
quantities on the right-hand side in (26) that $B$ needs to balance out, which
are the second and third terms. This suggests that we should take (up to
constants):

$B\geq\max\bigg{\\{}\frac{1}{p}\sum_{i=1}^{n}\mu_{i},\,\sum_{i=1}^{n}\mu_{i}^{2}\bigg{\\}}.$
(28)

Table 1 shows this “minimum value” of $B$ for $n=1$, and various combinations
of $p$ and $\mu$. For fixed $\mu$, the value of $B$ in (28) grows with $1/p$;
for fixed $p$, it grows with $\mu^{2}$.

Table 1: Minimum value for $B$ as in (28) for different combinations of $\mu$ and $p$, when $n=1$.  $p=0.1$ | $\mu$ | 0.5 | 1 | 5 | 10 | 30 | 100
---|---|---|---|---|---|---|---
$B$ | 5 | 10 | 50 | 100 | 900 | 10000
$\mu=5$ | $p$ | 0.001 | 0.005 | 0.01 | 0.1 | 0.3 | 0.5
$B$ | 5000 | 1000 | 500 | 50 | 25 | 25

Of course, in practice, the mean vector $\mu\in\mathbb{R}^{n}$ is unknown, so
our recommendation is to use a plug-in estimate, simply replacing $\mu$ with
the observed data vector $Y\in\mathbb{R}^{n}$, in (28).

##  5 Simulated experiments

In this section, we run two sets of simulations. The first compares the
unbiased estimator (UE) in (10) and the CB estimator in (17), across different
data models, loss functions, algorithms $g$, and binomial noise levels $p$. We
find that CB performs favorably overall, delivering similar error estimates to
UE for small values of $p$, and importantly, CB can have much smaller variance
than UE when $g$ is unstable. The second set fixes the data model and
algorithm $g$, and focuses on the performance of CB and UE after stratifying
the two to have roughly equal computational cost—accomplished by sampling
summands in (11) or (12). We find that CB generally performs favorably unless
the underlying mean has very large components.

###  5.1 Comparing CB and UE

Here we compare the CB and UE estimators, for squared and deviance loss
functions: see (18), (19) for CB and (11), (12) for UE. Throughout, we set
$n=100$, and use $B=100$ bootstrap samples for CB. We consider the following
combinations of different data models for $Y$, and algorithms $g$:

  * •

_Denoising._ We draw $Y_{i}\sim\mathrm{Pois}(\mu_{i})$, independently,
$i=1,\dots,n$, where $\mu_{i}=10$ for $i\leq 10$ and $\mu_{i}=0.5$ for $i>10$.
In this setting, we take $g$ to be a 1-step improvement on an empirical Bayes
(EB) estimator as described in Brown et al. (2013), with tuning parameter
fixed at $h=0.85$.

  * •

_Low-dimensional regression._ We set $p=10$, draw features $X_{i}\sim
N(\theta,I_{p})$, independently, $i=1,\dots,n$, where each $\theta_{j}=3$ and
$I_{p}$ denotes the $p\times p$ identity matrix; then we draw responses
$Y_{i}\sim\mathrm{Pois}(X_{i}^{\mathsf{T}}\beta)$, independently,
$i=1,\dots,n$, where each $\beta_{j}=0.05$. This corresponds to a signal-to-
noise ratio (SNR) of approximately 2\. We examine two choices for $g$: Poisson
regression and a regression tree.

  * •

_High-dimensional regression._ We set $p=200$, and use a similar setup to the
above, except with features $X_{i}\sim N(0,\sigma^{2}I_{p})$, where
$\sigma^{2}=1.5$, and responses
$Y_{i}\sim\mathrm{Pois}(X_{i}^{\mathsf{T}}\beta)$, where each
$\beta_{j}=0.13$. This was done to maintain an SNR of roughly 2. In this
setting, we take $g$ to be a lasso Poisson regression, with the tuning
parameter $\lambda$ chosen by 5-fold cross-validation (CV).

In each setting, we perform 100 repetitions (i.e., we draw the data vector
$Y\in\mathbb{R}^{n}$ 100 times from the specified data model); but we note
that in the regression settings, the features are drawn once and fixed
throughout. We consider a range of noise levels $p$ for the CB estimator:
0.05, 0.1, 0.3, 0.5, and 0.7. All error metrics and error estimators, here and
throughout all empirical examples, are scaled by $1/n$.

(a) Low-dim regression, $g=$ Poisson regression

(b) Low-dim regression, $g=$ regression tree

(c) Denoising, $g=$ EB 1-step estimator

(d) High-dim regression, $g=$ 5-fold CV lasso

Figure 3: Comparison of CB and UE across different data models, algorithms,
and loss functions.

The results are displayed in Figure 3, with each panel (a)–(d) displaying a
different combination of data model and algorithm $g$. In each panel, the
average test error estimate is displayed for each method (CB or UE), as well
as standard errors measured over the 100 repetitions. Furthermore, the black
points denote the true estimands (computed via Monte Carlo): $\mathrm{Err}(g)$
for UE, and $\mathrm{Err}_{p}(g)$ for CB. As expected, all estimators are seen
to be roughly unbiased for their targets: $\mathrm{Err}(g)$ or
$\mathrm{Err}(p)$. Interestingly, we also see that for squared loss, the
target $\mathrm{Err}_{p}(g)$ clearly decreases as $p$ grows, but for deviance
loss, the behavior of $\mathrm{Err}_{p}(g)$ tends to be more robust to growing
$p$ (and can increase or decrease, depending on the setting).

In panel (a), which is the low-dimensional regression setting, with Poisson
regression as the algorithm $g$, we can see that UE has a noticeably lower
standard error than CB at the lowest binomial noise level of $p=0.05$,
particularly for deviance loss. This is the only setting in which this
happens. In all others, the CB estimator at the lowest noise level has either
comparable or smaller variability than UE. In fact, in panel (d), which is the
high-dimensional regression setting, with CV-tuned lasso as the algorithm $g$,
we see that UE has a dramatically higher standard error than CB at any level
of noise $p$. The algorithm $g$ is inherently unstable here, because CV
(operating in high-dimensions, and at a moderate SNR) can choose very
different tuning parameter values across different data instances. Despite
this, CB is able to deliver estimates of reasonably low variance, since it
averages across draws of auxiliary binomial noise, which acts as a method of
smoothing (like bagging). We note that the analogous phenomenon also occurs in
the Gaussian setting, as observed by Oliveira et al. (2021).

###  5.2 Sampling summands

The unbiased estimator in (10) requires $n+1$ runs of the algorithm $g$,
making it computationally expensive for large sample sizes. In contrast, the
number of runs of $g$ required by the CB estimator in (17) is $B$, which is a
user-chosen parameter (recall that any choice of $B$ results in an unbiased
estimator $\mathrm{CB}_{p}(g)$ for $\mathrm{Err}_{p}(g)$, whereas larger $B$
reduces the variance of the estimator). In the last subsection, we fixed
$n=B=100$. In this one, we consider larger much sample sizes, with $n$ ranging
from $10^{3}$ to $10^{5}$. We maintain $B=100$, but we equate computational
costs between UE and CB by sampling $m=100$ summands uniformly at random (and
without replacement) from (11) or (12), and then scaling up the resulting sum
by $n/m$. We denote the estimator resulting from this “sampling summands”
approach by $\mathrm{UE}_{\mathrm{ss}}(g)$, which is unbiased for
$\mathrm{Err}(g)$.

For the data model, we draw $Y_{i}\sim\mathrm{Pois}(\mu)$, independently,
$i=1,\dots,n$, where $\mu$ ranges from 0.5 to 30. For the algorithm, we use a
simple linear shrinkage estimator:
$g(y)=0.8y+0.2\bar{y}+0.011\\{\bar{y}=0\\}$. We note that this choice is
generally favorable to UE, and more unstable algorithms $g$ would only create
more variability for UE relative to CB, and thus look more favorable to CB, as
observed in the last subsection. For the binomial noise parameter, we set

$p=\min\bigg{\\{}0.1,\sum_{i=1}^{n}\mu_{i}/\sum_{i=1}^{n}\mu_{i}^{2}\bigg{\\}},$

which roughly controls the leading terms in the reducible variance upper bound
(26).

(a) Deviance loss

(b) Squared loss

Figure 4: Variability of CB, $\mathrm{UE}_{\mathrm{ss}}$, and UE as functions
of $n$ and $\mu$, for a simple linear shrinkage estimator.

The results are displayed in Figure 4. For deviance loss, the results are
overall quite favorable for CB: it has a lower variance than
$\mathrm{UE}_{\mathrm{ss}}$ (darker shade of blue) in all but the top left
corner, which corresponds to small $n$ and large $\mu$. In fact, the
variability of CB is quite similar (across all $n,\mu$) to that of UE for
deviance loss, even though the former is considerably cheaper ($B=100$ runs of
the algorithm $g$, versus $n$ runs). For squared loss, there is more of a
clear tradeoff: for large values of $\mu$ (i.e., roughly $\log\mu>2$, or
$\mu>7.38$), we see that CB has greater variability than
$\mathrm{UE}_{\mathrm{ss}}$; for moderate values of $\mu$ (roughly $\log\mu$
between 0 and 2, or $\mu$ between 1 and 7.38), CB has comparable variability
for small $n$ and smaller variability for large $n$; while for small values of
$\mu$ (roughly $\log\mu<0$, or $\mu<1$), CB has smaller variability than
$\mathrm{UE}_{\mathrm{ss}}$.

##  6 Applications

###  6.1 Image denoising

We consider the following Poisson image denoising framework from Harmany et
al. (2012) (motivated by the study of Poisson noise or shot noise in areas
such as microscopy and astrophotography). We observe data
$Y_{i}\sim\mathrm{Pois}(f^{*}_{i})$, independently, $i=1,\dots,n$, where
$f^{*}\in\mathbb{R}^{n}_{+}$ is an unknown signal of interest, which we assume
has the structure of an $N\times N$ image, where $n=\sqrt{N}$. We consider an
estimator $\hat{f}=g(Y)$ for $f^{*}$ given by solving the optimization
problem:

$\mathop{\mathrm{minimize}}_{f\geq
0}\;\sum_{i=1}^{n}\big{(}-Y_{i}\log(f_{i}+\rho)+f_{i}+\rho\big{)}+\tau\sum_{i\sim
j}|f_{i}-f_{j}|,$ (29)

where $\rho$ is a small positive constant to avoid the singularity $f=0$, and
we write $i\sim j$ to indicate that indices $i,j$ are adjacent to each other
in the ordering determined by the underlying image. This estimator is a form
of total variation (TV) regularized Poisson image denoising; we note that the
loss term is equivalent (up to constants) to the Poisson deviance between $f$
and $Y$, and the penalty term encourages the estimated image $\hat{f}$ to be
piecewise constant, with the tuning parameter $\tau\geq 0$ determining the
strength of regularization.

As an example, we consider the well-known synthetic phantom image $f^{*}$, of
resolution $128\times 128$ (so that $n=16384$). We use CB to estimate the test
error of $\hat{f}$, the Poisson image denoising estimator defined in (29),
over a range of values of the tuning parameter $\tau$. We consider both
squared (18) and deviance (19) loss, and set $B=50$ and $p=0.1$. We did not
consider the unbiased estimator in (11) or (12) in this experiment, due to its
prohibitive computational cost (it requires $n=16384$ refits of the TV
denoising estimator (29)).

Figure 5: Comparison of CB and true test error curves, as functions of the
tuning parameter $\tau$, for a Poisson image denoising estimator. Figure 6:
Original, noisy, and denoised estimates at the CB-optimal value of $\tau$, for
squared and deviance loss.

Figure 5 displays the CB error curve and true error curve (approximated by
Monte Carlo), as functions of $\tau$, with separate panels for squared and
deviance loss. There are two points worth noting. First, despite the gap
between the CB and true test error curves (unsurprising, because $p=0.1$),
their curvature is similar; in particular, the value of $\tau$ minimizing the
CB curve is close to the value minimizing test error. This is the case for
both squared and deviance loss, and it shows that the CB estimator can be
useful for model tuning, even when $p$ is not small. Second, the value of
$\tau$ minimizing the CB curve is larger for deviance loss than it is for
squared loss, marked by the dotted lines in each panel. This translates into a
greater degree of regularization, as can be seen clearly in Figure 6, which
plots the denoised estimates themselves at the CB-optimal values of $\tau$,
for squared and deviance loss.

###  6.2 Density estimation

We study density estimation, which can be turned into a Poisson regression via
Lindsey’s method (Lindsey, 1974; Lindsey and Mersch, 1992; Efron and
Tibshirani, 1996). The basic idea is to discretize the domain into bins, and
model the count in each bin as a Poisson random variable, with the mean
parameter constrained or regularized to be smoothly varying across the bins.
We can then estimate the mean parameter by (regularized) maximum likelihood,
which in turn gives a discretized density estimate.

In particular, we consider an example from Phillips et al. (2006) on the
distribution of _Bradypus variegatus_ , a lowland species of sloth found
across Central and South America. Each data point consists of latitude and
longitude pair, representing a site where a sloth was seen, and the data set
contains 116 total sightings. To form a 2d density estimate, we apply
Lindsey’s method, with 200 equally-spaced bins along the latitude and
longitude axes, and we use a P-spline to model the Poisson mean parameter.
P-splines were first proposed by Eilers and Marx (1996); details for the 2d
case can be found in Eilers and Marx (2003); Eilers et al. (2006); Eilers and
Marx (2021). We use a 2d cubic B-spline parametrization for the mean function
with 30 knots in each dimension, and we use a penalty on the sum of squared
second-order differences across adjacent B-spline parameters along each
dimension. Moreover, we consider two versions of this penalty: an
_anisotropic_ version, which decouples the regularization strength along each
dimension, and has two tuning parameters $\lambda_{1},\lambda_{2}\geq 0$; and
an _isotropic_ version, which ties together the regularization strength over
the dimensions, and has a single tuning parameter $\lambda\geq 0$.

Figure 8 shows the results of using the CB method, with $B=100$ and $p=0.1$,
to estimate both squared and deviance loss, across a range of tuning parameter
values. For either the anisotropic or isotropic penalty, it is clear that
minimizing CB-estimated deviance loss leads to larger tuning parameter
values—and hence more regularized density estimates—than minimizing CB-
estimated squared loss. As we can see in Figure 8, this leads to more
plausible looking density estimates (top row). The estimates obtained
optimizing CB-estimated squared loss (bottom row) appear too concentrated
around the observations themselves.

(a) Anisotropic penalty, deviance loss

(b) Anisotropic penalty, squared loss

(c) Isotropic penalty

Figure 7: CB curves for anisotropic and isotropic penalties as a function of
the tuning parameter(s).

Figure 8: Density estimates at the CB-optimized values of the tuning
parameters for anisotropic and isotropic penalties, and squared and deviance
loss.

##  7 Discussion

We proposed and studied a coupled bootstrap (CB) method for test error
estimation in the Poisson means problem, with a focus on squared and Poisson
deviance loss functions. The CB estimator, for any choice of the binomial
noise parameter $p>0$, is unbiased for an intuitive target:
$\mathrm{Err}_{p}(g)$, the test error of the given algorithm $g$, when the
mean vector in the Poisson model has been shrunk from $\mu$ to $(1-p)\mu$.
Importantly, this unbiasedness requires no assumptions on $g$ whatsoever.
Furthermore, we proved that in the noiseless limit $p\to 0$, the CB estimator
(with infinite bootstrap iterations) reduces to the natural unbiased estimator
(UE) for test error that comes from an application of Hudson’s lemma. However,
CB has two key advantages over UE. First, it requires running the algorithm
$g$ in question $B$ times (which is a user-controlled parameter in CB), versus
$n+1$ times (which comes directly from the form of UE). Second, as we show in
our experiments, CB can often have smaller variance than UE, particularly when
the underlying algorithm $g$ is unstable.

We finish by emphasizing that it would be interesting to extend the CB method
to other data models—beyond Gaussian, as in Oliveira et al. (2021), and
Poisson, as in the current paper. For this, we would need a construction
analogous to Lemma 2, above. Given a data vector $Y$ from a particular model,
we need to be able to construct $Y^{*},Y^{\dagger}$ that are independent, and
satisfy $\mathbb{E}[Y^{*}]=\mathbb{E}[Y^{\dagger}]$; and furthermore, the
distribution of $Y^{*}$ needs to be “close enough” to that of $Y$ such that
modeling $Y^{*}$ would be considered an “interesting” proxy problem. In the
Gaussian problem, our construction yields $Y^{*}$ with an inflated variance,
and in the Poisson problem, it yields $Y^{*}$ with a shrunken mean. It remains
to find appropriate auxiliary randomization schemes for other distributions,
with binary classification (the binomial distribution) remaining a case that
is of central interest but also highly challenging.

## Acknowledgements

The authors thank Aaditya Ramdas and Boyan Duan for helpful and inspiring
discussions. NLO was supported by an Amazon Fellowship.

## References

  * Akaike (1973) Hirotogu Akaike.  Information theory and an extension of the maximum likelihood principle.  In _Second International Symposium on Information Theory_ , pages 267–281, 1973.
  * Bates et al. (2021) Stephen Bates, Trevor Hastie, and Robert Tibshirani.  Cross-validation: What does it estimate and how well does it do it?  arXiv: 2104.00673, 2021.
  * Breiman (1992) Leo Breiman.  The little bootstrap and other methobds for dimensionality selection in regression: X-fixed prediction error.  _Journal of the American Statistical Association_ , 87(419):738–754, 1992.
  * Brown et al. (2013) Lawrence D. Brown, Eitan Greenshtein, and Ya’acov Ritov.  The Poisson compound decision problem revisited.  _Journal of the American Statistical Association_ , 108(502):741–749, 2013.
  * Cao and Xie (2016) Yang Cao and Yao Xie.  Poisson matrix recovery and completion.  _IEEE Transactions on Signal Processing_ , 64(6):1609–1620, 2016.
  * Deledalle (2017) Charles-Alban Deledalle.  Estimation of Kullback-Leibler losses for noisy recovery problems within the exponential family.  _Electronic Journal of Statistics_ , 11(2):3141–3164, 2017.
  * Efron (1975) Bradley Efron.  Defining the curvature of a statistical problem (with applications to second order efficiency).  _Annals of Statistics_ , 3(6):1189–1242, 1975\.
  * Efron (1986) Bradley Efron.  How biased is the apparent error rate of a prediction rule?  _Journal of the American Statistical Association_ , 81(394):461–470, 1986.
  * Efron (2004) Bradley Efron.  The estimation of prediction error: Covariance penalties and cross-validation.  _Journal of the American Statistical Association_ , 99(467):619–632, 2004.
  * Efron and Tibshirani (1996) Bradley Efron and Robert Tibshirani.  Using specially designed exponential families for density estimation.  _Annals of Statistics_ , 24(6):2431–2461, 1996\.
  * Eilers and Marx (1996) Paul H. C. Eilers and Brian D. Marx.  Flexible smoothing with B-splines and penalties.  _Statistical Science_ , 11(2):89–121, 1996.
  * Eilers and Marx (2003) Paul H. C. Eilers and Brian D. Marx.  Multivariate calibration with temperature interaction using two-dimensional penalized signal regression.  _Chemometrics and Intelligent Laboratory Systems_ , 66:159–174, 2003.
  * Eilers and Marx (2021) Paul H. C. Eilers and Brian D. Marx.  _Practical Smoothing: The Joys of P-splines_.  Cambridge University Press, 2021.
  * Eilers et al. (2006) Paul H. C. Eilers, Iain D. Currie, and Maria Durbán.  Fast and compact smoothing on large multidimensional grids.  _Computational Statistics & Data Analysis_, 50(1):61–76, 2006.
  * Eldar (2009) Yonina C. Eldar.  Generalized SURE for exponential families: Applications to regularization.  _IEEE Transactions on Signal Processing_ , 57(2):471–481, 2009.
  * Harmany et al. (2009) Zachary T. Harmany, Roummel F. Marcia, and Rebecca M. Willett.  Sparse Poisson intensity reconstruction algorithms.  In _IEEE Workshop on Statistical Signal Processing_ , 2009.
  * Harmany et al. (2012) Zachary T. Harmany, Roummel F. Marcia, and Rebecca M. Willett.  This is SPIRAL-TAP: Sparse Poisson intensity reconstruction algorithms–theory and practice.  _IEEE Transactions on Image Processing_ , 21(3):1084–1096, 2012.
  * Hudson (1978) H. Malcolm Hudson.  A natural identity for exponential families with applications in multiparameter estimation.  _Annals of Statistics_ , 6(3):473–484, 1978.
  * Leiner et al. (2021) James Leiner, Boyan Duan, Larry Wasserman, and Aaditya Ramdas.  Data fission: Splitting a single data point.  arXiv: 2112.11079, 2021.
  * Lindsey (1974) J. K. Lindsey.  Comparison of probability distributions.  _Journal of the Royal Statistical Society: Series B_ , 36(1):38–47, 1974.
  * Lindsey and Mersch (1992) J. K. Lindsey and G. Mersch.  Fitting and comparing probability distributions with log linear models.  _Computational Statistics & Data Analysis_, 13(4):373–384, 1992.
  * Luisier et al. (2010) Florian Luisier, Cédric Vonesch, Thierry Blu, and Michael Unser.  Fast interscale wavelet denoising of Poisson-corrupted images.  _Signal Processing_ , 90(2):415–427, 2010.
  * Mallows (1973) Colin Mallows.  Some comments on $C_{p}$.  _Technometrics_ , 15(4):661–675, 1973.
  * Neufeld et al. (2023) Anna Neufeld, Ameer Dharamshi, Lucy L. Gao, and Daniela Witten.  Data thinning for convolution-closed distributions.  arXiv: 2301.07276, 2023.
  * Oliveira et al. (2021) Natalia L. Oliveira, Jing Lei, and Ryan J. Tibshirani.  Unbiased risk estimation in the normal means problem via coupled bootstrap techniques.  arXiv: 2111.09447, 2021.
  * Phillips et al. (2006) Steven J. Phillips, Robert P. Anderson, and Robert E. Schapire.  Maximum entropy modeling of species geographic distributions.  _Ecological Modelling_ , 190(3):231–259, 2006\.
  * Raginsky et al. (2010) Maxim Raginsky, Rebecca M. Willett, Zachary T. Harmany, and Roummel F. Marcia.  Compressed sensing performance bounds under Poisson noise.  _IEEE Transactions on Signal Processing_ , 58(8):3990–4002, 2010.
  * Rosset and Tibshirani (2020) Saharon Rosset and Ryan J. Tibshirani.  From fixed-X to random-X regression: Bias-variance decompositions, covariance penalties, and prediction error estimation.  _Journal of the American Statistical Association_ , 15(529):138–151, 2020.
  * Salmon et al. (2014) Joseph Salmon, Zachary Harmany, Charles-Alban Deledalle, and Rebecca Willett.  Poisson noise reduction with non-local PCA.  _Journal of Mathematical Imaging and Vision_ , 48(2):279–294, 2014.
  * Shen et al. (2004) Xiaotong Shen, Hsin-Cheng Huang, and Jimmy Ye.  Adaptive model selection and assessment for exponential family distributions.  _Technometrics_ , 46(3):306–317, 2004.
  * Stein (1981) Charles Stein.  Estimation of the mean of a multivariate normal distribution.  _Annals of Statistics_ , 9(6):1135–1151, 1981\.
  * Ye (1998) Jianming Ye.  On measuring and correcting the effects of data mining and model selection.  _Journal of the American Statistical Association_ , 93(441):120–131, 1998.

##  Appendix A Proofs and additional details for Section 2

###  A.1 Proof of Lemma 2

The joint probability mass function of $(Y,\omega)$ is given by

$\displaystyle\mathbb{P}(Y=y,\omega=w)$
$\displaystyle=\mathbb{P}(\omega=k\,|\,Y=y)\mathbb{P}(Y=y)$
$\displaystyle=\binom{y}{k}p^{k}(1-p)^{y-k}\frac{\mu^{y}e^{-\mu}}{y!},$

for all $y\geq 0$ and $k\in\\{0,\dots,y\\}$. The probability mass function of
$(U,V)=(Y-\omega,\omega)$ is thus

$\displaystyle\mathbb{P}(U=u,V=v)$ $\displaystyle=\mathbb{P}(Y=u+v,\omega=v)$
$\displaystyle=\binom{u+v}{v}p^{v}(1-p)^{u}\frac{\mu^{u+v}e^{-\mu}}{(u+v)!}$
$\displaystyle=\frac{1}{u!v!}p^{v}(1-p)^{u}\mu^{u}\mu^{v}e^{-(p+1-p)\mu}$
$\displaystyle=\frac{((1-p)\mu)^{u}e^{-(1-p)\mu}}{u!}\frac{(p\mu)^{v}e^{-p\mu}}{v!}.$

This shows that $U$ and $V$ are independent $\mathrm{Pois}((1-p)\mu)$ and
$\mathrm{Pois}(p\mu)$ random variables, respectively, which proves the desired
result.

###  A.2 Divergence deviance terms in UE versus CB

For deviance loss, a summand in the unbiased estimator (12) is undefined if
$Y_{i}\not=0$ and $g_{i}(Y-e_{i})=0$, and a summand in the CB estimator (19)
is undefined if $Y^{\dagger}_{i}\not=0$ and $g_{i}(Y^{*})=0$, where
$Y^{*},Y^{\dagger}$ are a sample from (16) (we have hidden the dependence on
$b$). Suppose for simplicity that $n=1$ and $g(0)=0$. We can then compute,
under $Y\sim\mathrm{Pois}(\mu)$, the probability with which this happens for
the unbiased and CB estimators. For the unbiased estimator, this is:

$\mathbb{P}(Y=1)=e^{-\mu}\mu.$

For the CB estimator, this is:

$\displaystyle\sum_{y=1}^{\infty}\mathbb{P}(Y^{*}=0\,|\,Y=y)\mathbb{P}(Y=y)$
$\displaystyle=\sum_{y=1}^{\infty}\mathbb{P}(\omega=Y\,|\,Y=y)\mathbb{P}(Y=y)$
$\displaystyle=\sum_{y=1}^{\infty}p^{y}e^{-\mu}\mu^{y}/y!$
$\displaystyle=e^{-(1-p)\mu}\sum_{y=1}^{\infty}e^{-p\mu}(p\mu)^{y}/y!$
$\displaystyle=e^{-(1-p)\mu}(1-e^{-p\mu})$
$\displaystyle=e^{-\mu}(e^{p\mu}-1).$

Figure 9 plots these two probabilities, a functions of $\mu$, for $p$ ranging
over 0.01, 0.1, 0.3, 0.5. For small $p$, it is clear that the CB estimator has
much lower probability of being ill-defined than the unbiased estimator.

Figure 9: Comparison of the probabilities of an individual summand from the
unbiased and CB estimators being ill-defined, as functions of the mean $\mu$.
The four panels show different values of $p$.

###  A.3 Proof of Proposition 3

In this proof we use $f(Y_{p},{\tilde{Y}}_{p})$ to denote
$D_{\phi}({\tilde{Y}}_{p},g(Y_{p}))$. First, we show that the map is
continuous. For any $p\in[0,1)$,

$\displaystyle\lim_{t\to p}\mathbb{E}[f(Y_{t},{\tilde{Y}}_{t})]$
$\displaystyle=\lim_{t\to
p}\sum_{y_{1},\dots,y_{n}=0}^{\infty}\sum_{{\tilde{y}}_{1},\dots,{\tilde{y}}_{n}=0}^{\infty}f(y,{\tilde{y}})\frac{e^{-2(1-t)\sum_{i=1}^{n}\mu_{i}}(\prod_{i=1}^{n}\mu_{i}^{y_{i}+{\tilde{y}}_{i}})(1-t)^{\sum_{i=1}^{n}y_{i}+{\tilde{y}}_{i}}}{\prod_{i=1}^{n}y_{i}!{\tilde{y}}_{i}!}$
$\displaystyle=\sum_{y_{1},\dots,y_{n}=0}^{\infty}\sum_{{\tilde{y}}_{1},\dots,{\tilde{y}}_{n}=0}^{\infty}f(y,{\tilde{y}})\frac{e^{-2(1-p)\sum_{i=1}^{n}\mu_{i}}(\prod_{i=1}^{n}\mu_{i}^{y_{i}+{\tilde{y}}_{i}})(1-p)^{\sum_{i=1}^{n}y_{i}+{\tilde{y}}_{i}}}{\prod_{i=1}^{n}y_{i}!{\tilde{y}}_{i}!}$
$\displaystyle=\mathbb{E}[f(Y_{p},{\tilde{Y}}_{p})].$

To switch the infinite sum and the limit, we used the dominated convergence
theorem (DCT). The dominating function is given by

$h({\tilde{y}},y)=f(y,{\tilde{y}})\frac{\prod_{i=1}^{n}\mu_{i}^{y_{i}+{\tilde{y}}_{i}}}{\prod_{i=1}^{n}y_{i}!{\tilde{y}}_{i}!},$

which is integrable by assumption:

$\displaystyle\sum_{y_{1},\dots,y_{n}=0}^{\infty}\sum_{{\tilde{y}}_{1},\dots,{\tilde{y}}_{n}=0}^{\infty}h({\tilde{y}},y)$
$\displaystyle=\sum_{y_{1},\dots,y_{n}=0}^{\infty}\sum_{{\tilde{y}}_{1},\dots,{\tilde{y}}_{n}=0}^{\infty}f(y,{\tilde{y}})\frac{\prod_{i=1}^{n}\mu_{i}^{y_{i}+{\tilde{y}}_{i}}}{\prod_{i=1}^{n}y_{i}!{\tilde{y}}_{i}!}e^{-2\sum_{i=1}^{n}\mu_{i}}e^{2\sum_{i=1}^{n}\mu_{i}}$
$\displaystyle=e^{2\sum_{i=1}^{n}\mu_{i}}\mathbb{E}[f(Y_{0},{\tilde{Y}}_{0})]<\infty.$

Next, for the first derivative, note that

$\displaystyle\frac{\partial}{\partial p}\mathbb{E}[f(Y_{p},{\tilde{Y}}_{p})]$
$\displaystyle=\frac{\partial}{\partial
p}\sum_{y_{1},\dots,y_{n}=0}^{\infty}\sum_{{\tilde{y}}_{1},\dots,{\tilde{y}}_{n}=0}^{\infty}f(y,{\tilde{y}})\frac{e^{-2(1-p)\sum_{i=1}^{n}\mu_{i}}\prod_{i=1}^{n}\mu_{i}^{y_{i}+{\tilde{y}}_{i}}(1-p)^{\sum_{i=1}^{n}y_{i}+{\tilde{y}}_{i}}}{\prod_{i=1}^{n}y_{i}!{\tilde{y}}_{i}!}$
$\displaystyle=\sum_{y_{1},\dots,y_{n}=0}^{\infty}\sum_{{\tilde{y}}_{1},\dots,{\tilde{y}}_{n}=0}^{\infty}\frac{f(y,{\tilde{y}})\prod_{i=1}^{n}\mu_{i}^{y_{i}+{\tilde{y}}_{i}}}{\prod_{i=1}^{n}y_{i}!{\tilde{y}}_{i}!}\frac{\partial}{\partial
p}e^{-2(1-p)\sum_{i=1}^{n}\mu_{i}}(1-p)^{\sum_{i=1}^{n}y_{i}+{\tilde{y}}_{i}}$
$\displaystyle=\sum_{y_{1},\dots,y_{n}=0}^{\infty}\sum_{{\tilde{y}}_{1},\dots,{\tilde{y}}_{n}=0}^{\infty}\frac{f(y,{\tilde{y}})\prod_{i=1}^{n}\mu_{i}^{y_{i}+{\tilde{y}}_{i}}}{\prod_{i=1}^{n}y_{i}!{\tilde{y}}_{i}!}\bigg{(}2\sum_{i=1}^{n}\mu_{i}e^{-2(1-p)\sum_{i=1}^{n}\mu_{i}}(1-p)^{\sum_{i=1}^{n}y_{i}+{\tilde{y}}_{i}}$
$\displaystyle\qquad-\sum_{i=1}^{n}(y_{i}+{\tilde{y}}_{i})e^{-2(1-p)\sum_{i=1}^{n}\mu_{i}}(1-p)^{\sum_{i=1}^{n}y_{i}+{\tilde{y}}_{i}-1}\bigg{)}$
$\displaystyle=2\sum_{i=1}^{n}\mu_{i}\mathbb{E}[f(Y_{p},{\tilde{Y}}_{p})]-\frac{1}{1-p}\mathbb{E}\big{[}f(Y_{p},{\tilde{Y}}_{p})\langle{\tilde{Y}}_{p}+Y_{p},1_{n}\rangle\big{]},$

where we used DCT to switch the sums and derivative, using a similar
dominating function as above and recognizing that the summand is Lipschitz in
$p$ with Lipschitz constant depending on $\mu,y,{\tilde{y}}$. Now to prove
continuity of the first derivative, we apply the above continuity result with
$f(Y_{p},{\tilde{Y}}_{p})\langle{\tilde{Y}}_{p}+Y_{p},1_{n}\rangle$ in place
of $f(Y_{p},{\tilde{Y}}_{p})$. For $k^{\text{th}}$ derivatives, the argument
follows from sequential applications of the same continuity result and similar
derivative calculations.

##  Appendix B Proof of Theorem 1

We start with a lemma that contains two key results to be used in the proof of
Theorem 1.

######  Lemma 3.

Let $h:\mathbb{Z}^{n}_{+}\to\mathbb{R}^{n}_{+}$, and set $h(z)=0$ for
$z\notin\mathbb{Z}^{n}_{+}$. Fix any $y\in\mathbb{Z}_{+}^{n}$, and draw
$\omega_{i}\sim\mathrm{Binom}(y_{i},p)$, independently, for $i=1,\dots,n$,
where $p\in[0,1)$. Then, for each $i=1,\dots,n$,

  1. (a)

$\lim_{p\to 0}\mathbb{E}[h_{i}(y-\omega)]=h_{i}(y)$;

  2. (b)

$\lim_{p\to
0}\frac{1-p}{p}\cdot\mathbb{E}[\omega_{i}h_{i}(y-\omega)]=y_{i}h_{i}(y-e_{i})$.

Recall, we use $e_{i}\in\mathbb{R}^{n}$ to denote the vector whose
$i^{\text{th}}$ entry is 1, with all others 0.

###### Proof.

Define the following sets:

$\displaystyle\Omega$
$\displaystyle=\\{(\omega_{1},\dots,\omega_{n}):\omega_{j}\in\\{0,\dots,y_{j}\\},j=1,\dots,n\\},$
$\displaystyle\Omega_{\setminus 0}$ $\displaystyle=\Omega\setminus\\{0\\},$
$\displaystyle\Omega_{i0}$
$\displaystyle=\\{(\omega_{1},\dots,\omega_{n}):\omega_{i}=0,\omega_{j}\in\\{0,\dots,y_{j}\\},j\not=i\\},$
$\displaystyle\Omega_{i10}$
$\displaystyle=\\{(\omega_{1},\dots,\omega_{n}):\omega_{i}=1,\omega_{j}=0,j\not=i\\},$
$\displaystyle\Omega_{i1\bar{0}}$
$\displaystyle=\\{(\omega_{1},\dots,\omega_{n}):\omega_{i}\geq
1,\omega_{j}\in\\{0,\dots,y_{j}\\},j\not=i\\}\setminus\Omega_{i10}.$

For the first result (a), we have

$\displaystyle\lim_{p\to 0}\mathbb{E}[h_{i}(y-\omega)]$
$\displaystyle=\lim_{p\to
0}\sum_{\omega\in\Omega}h_{i}(y-\omega)p^{\sum_{k=1}^{n}\omega_{k}}(1-p)^{\sum_{k=1}^{n}y_{k}-\omega_{k}}\prod_{j=1}^{n}\binom{y_{j}}{\omega_{j}}$
$\displaystyle=\lim_{p\to 0}h_{i}(y)(1-p)^{\sum_{k=1}^{n}y_{k}}$
$\displaystyle\qquad+\sum_{\omega\in\Omega_{\setminus 0}}\lim_{p\to
0}h_{i}(y-\omega)p^{\sum_{k=1}^{n}\omega_{k}}(1-p)^{\sum_{k=1}^{n}y_{k}-\omega_{k}}\prod_{j=1}^{n}\binom{y_{j}}{\omega_{j}}$
$\displaystyle=h_{i}(y),$

since for $\omega\in\Omega_{\setminus 0}$, we have that
$\sum_{k=1}^{n}\omega_{k}>0$. For the second result (b),

$\displaystyle\lim_{p\to 0}\frac{1-p}{p}\mathbb{E}[\omega_{i}h_{i}(y-\omega)]$
$\displaystyle=\lim_{p\to
0}\frac{1-p}{p}\sum_{\omega\in\Omega}\omega_{i}h_{i}(y-\omega)p^{\sum_{k=1}^{n}\omega_{k}}(1-p)^{\sum_{k=1}^{n}y_{k}-\omega_{k}}\prod_{j=1}^{n}\binom{y_{j}}{\omega_{j}}$
$\displaystyle=\lim_{p\to
0}\sum_{\omega\in\Omega}\omega_{i}h_{i}(y-\omega)p^{\sum_{k=1}^{n}\omega_{k}-1}(1-p)^{1+\sum_{k=1}^{n}y_{k}-\omega_{k}}\prod_{j=1}^{n}\binom{y_{j}}{\omega_{j}}$
$\displaystyle=\lim_{p\to
0}\sum_{\omega\in\Omega_{i10}}\omega_{i}h_{i}(y-\omega)p^{\sum_{k=1}^{n}\omega_{k}-1}(1-p)^{1+\sum_{k=1}^{n}Y_{k}-\omega_{k}}\prod_{j=1}^{n}\binom{y_{j}}{\omega_{j}}$
$\displaystyle=\lim_{p\to
0}h_{i}(y-e_{i})p^{0}(1-p)^{1+\sum_{k=1}^{n}y_{k}-\omega_{k}}\binom{y_{i}}{1}\prod_{j\neq
i,j=1}^{n}\binom{y_{j}}{0}$ $\displaystyle=y_{i}h_{i}(y-e_{i}),$

where we use the fact that
$\Omega=\Omega_{i0}\cup\Omega_{i10}\cup\Omega_{i1\bar{0}}$ and

$\lim_{p\to
0}\sum_{\omega\in\Omega_{i0}}\omega_{i}h_{i}(y-\omega)p^{\sum_{k=1}^{n}\omega_{k}-1}(1-p)^{1+\sum_{k=1}^{n}y_{k}-\omega_{k}}\prod_{j=1}^{n}\binom{y_{j}}{\omega_{j}}\\\
=\lim_{p\to
0}\sum_{\omega\in\Omega_{i0}}0h_{i}(y-\omega)p^{\sum_{k=1}^{n}\omega_{k}-1}(1-p)^{1+\sum_{k=1}^{n}y_{k}-\omega_{k}}\prod_{j=1}^{n}\binom{y_{j}}{\omega_{j}}=0$

as well as

$\lim_{p\to
0}\sum_{\omega\in\Omega_{i1\bar{0}}}\omega_{i}h_{i}(y-\omega)p^{\sum_{k=1}^{n}\omega_{k}-1}(1-p)^{1+\sum_{k=1}^{n}y_{k}-\omega_{k}}\prod_{j=1}^{n}\binom{y_{j}}{\omega_{j}}\\\
=\sum_{\omega\in\Omega_{i1\bar{0}}}\omega_{i}h_{i}(y-\omega)0^{\sum_{k=1}^{n}\omega_{k}-1}\prod_{j=1}^{n}\binom{y_{j}}{\omega_{j}}=0,$

since for $\omega\in\Omega_{i1\bar{0}}$, we have that
$\sum_{k=1}^{n}\omega_{k}-1>0$. ∎

Now we are ready to prove Theorem 1. We start by expanding the infinite
bootstrap estimator

$\displaystyle\mathbb{E}[\mathrm{CB}_{p}(g)\,|\,Y]$
$\displaystyle=\sum_{i=1}^{n}\mathbb{E}\Big{[}D_{\phi}(Y^{\dagger}_{i},g_{i}(Y^{*}))+\phi(Y^{*}_{i})-\phi(Y^{\dagger}_{i})\,\big{|}\,Y\Big{]}$
$\displaystyle=\sum_{i=1}^{n}\mathbb{E}\Big{[}\phi(Y^{*}_{i})-\phi(g_{i}(Y^{*}))-\nabla_{i}\phi(g(Y^{*}))(Y^{\dagger}_{i}-Y^{*}_{i})\,\big{|}\,Y\Big{]}$
$\displaystyle=\sum_{i=1}^{n}\mathbb{E}\bigg{[}\phi(Y_{i}-\omega_{i})-\phi(g_{i}(Y-\omega))-\nabla_{i}\phi(g(Y-\omega))\bigg{(}\frac{1-p}{p}\omega_{i}-(Y_{i}-\omega_{i})\bigg{)}\,\Big{|}\,Y\bigg{]}$
$\displaystyle=\sum_{i=1}^{n}\mathbb{E}\bigg{[}\phi(Y_{i}-\omega_{i})-\phi(g_{i}(Y-\omega))-\frac{1-p}{p}\nabla_{i}\phi(g(Y-\omega))\omega_{i}+\nabla_{i}\phi(g(Y-\omega))(Y_{i}-\omega_{i})\,\Big{|}\,Y\bigg{]}.$

Then, taking the limit in the last line,

$\displaystyle\sum_{i=1}^{n}\lim_{p\to 0}\,$
$\displaystyle\mathbb{E}\bigg{[}\phi(Y_{i}-\omega_{i})-\phi(g_{i}(Y-\omega))-\frac{1-p}{p}\nabla_{i}\phi(g(Y-\omega))\omega_{i}+\nabla_{i}\phi(g(Y-\omega))(Y_{i}-\omega_{i})\,\Big{|}\,Y\bigg{]}$
$\displaystyle=\sum_{i=1}^{n}\lim_{p\to
0}\mathbb{E}\Big{[}\phi(Y_{i}-\omega_{i})-\phi(g_{i}(Y-\omega))+\nabla_{i}\phi(g(Y-\omega))(Y_{i}-\omega_{i})\,\big{|}\,Y\Big{]}-Y_{i}\nabla_{i}\phi(g(Y-e_{i}))$
$\displaystyle=\sum_{i=1}^{n}\phi(Y_{i})-\phi(g_{i}(Y))+\nabla_{i}\phi(g(Y))(Y_{i})-Y_{i}\nabla_{i}\phi(g(Y-e_{i}))$
$\displaystyle=\mathrm{UE}(Y),$

where in the second-to-last and last lines we used Lemma 3 parts (b) and (a),
respectively. This completes the proof.

##  Appendix C Proofs for Section 4

###  C.1 Proof of Proposition 4

From Proposition 3, the mapping $p\mapsto\mathrm{Err}_{p}(g)$ has a continuous
derivative for $p\in[0,1)$. By an application of the fundamental theorem of
calculus, we can write the bias as

$\displaystyle\mathrm{Err}_{p}(g)$
$\displaystyle-\mathrm{Err}(g)=\int_{0}^{p}\frac{\partial}{\partial
t}\mathrm{Err}_{t}(g)\,dt$
$\displaystyle=\int_{0}^{p}\bigg{\\{}2\sum_{i=1}^{n}\mu_{i}\mathbb{E}[D_{\phi}({\tilde{Y}}_{t},g(Y_{t}))]-\frac{1}{(1-t)}\mathbb{E}\big{[}D_{\phi}({\tilde{Y}}_{t},g(Y_{t}))\langle{\tilde{Y}}_{t}+Y_{t},1_{n}\rangle\big{]}\bigg{\\}}\,dt$
$\displaystyle=\int_{0}^{p}\bigg{\\{}2\sum_{i=1}^{n}\mu_{i}\mathrm{Err}_{t}(g)-\frac{1}{(1-t)}\mathrm{Cov}\Big{(}D_{\phi}({\tilde{Y}}_{t},g(Y_{t})),\langle{\tilde{Y}}_{t}+Y_{t},1_{n}\rangle\Big{)}-\frac{1}{(1-t)}2\mathrm{Err}_{t}(g)(1-t)\sum_{i=1}^{n}\mu_{i}\bigg{\\}}\,dt$
$\displaystyle=-\int_{0}^{p}\frac{1}{(1-t)}\mathrm{Cov}\Big{(}D_{\phi}({\tilde{Y}}_{t},g(Y_{t})),\langle{\tilde{Y}}_{t}+Y_{t},1_{n}\rangle\Big{)}\,dt$
$\displaystyle=-\int_{0}^{p}\frac{1}{(1-t)}\mathrm{Cor}\Big{(}D_{\phi}({\tilde{Y}}_{t},g(Y_{t})),\langle{\tilde{Y}}_{t}+Y_{t},1_{n}\rangle\Big{)}\sqrt{\mathrm{Var}\big{[}D_{\phi}({\tilde{Y}}_{t},g(Y_{t}))\big{]}}\sqrt{\mathrm{Var}\big{[}\langle{\tilde{Y}}_{t}+Y_{t},1_{n}\rangle\big{]}}\,dt$
$\displaystyle=-\int_{0}^{p}\frac{1}{(1-t)}\mathrm{Cor}\Big{(}D_{\phi}({\tilde{Y}}_{t},g(Y_{t})),\langle{\tilde{Y}}_{t}+Y_{t},1_{n}\rangle\Big{)}\sqrt{\mathrm{Var}\big{[}D_{\phi}({\tilde{Y}}_{t},g(Y_{t}))\big{]}}\sqrt{(2(1-t)\sum_{i=1}^{n}\mu_{i}}\,dt$
$\displaystyle=-\int_{0}^{p}\frac{1}{\sqrt{1-t}}\mathrm{Cor}\big{[}D_{\phi}({\tilde{Y}}_{t},g(Y_{t})),\langle{\tilde{Y}}_{t}+Y_{t},1_{n}\rangle\Big{)}\sqrt{\mathrm{Var}\big{[}D_{\phi}({\tilde{Y}}_{t},g(Y_{t}))\big{]}}\,dt\cdot\sqrt{2\sum_{i=1}^{n}\mu_{i}}.$

which proves (24). Upper bounding the correlation by $1$, and using the
monotone variance assumption for $p\in[0,1/2]$, we have

$\displaystyle|\mathrm{Err}_{p}(g)-\mathrm{Err}(g)|$
$\displaystyle\leq\int_{0}^{p}\frac{1}{\sqrt{1-t}}\sqrt{\mathrm{Var}\big{[}D_{\phi}({\tilde{Y}}_{t},g(Y_{t}))\Big{)}}\,dt\cdot\sqrt{2\sum_{i=1}^{n}\mu_{i}}$
$\displaystyle\leq\int_{0}^{p}\frac{1}{\sqrt{1-t}}\sqrt{\mathrm{Var}\big{[}D_{\phi}({\tilde{Y}},g(Y))\Big{)}}\,dt\cdot\sqrt{2\sum_{i=1}^{n}\mu_{i}}$
$\displaystyle=\sqrt{\mathrm{Var}\big{[}D_{\phi}({\tilde{Y}},g(Y))\Big{)}}\sqrt{2\sum_{i=1}^{n}\mu_{i}}\int_{0}^{p}\frac{1}{\sqrt{1-t}}\,dt$
$\displaystyle=\sqrt{\mathrm{Var}\big{[}D_{\phi}({\tilde{Y}},g(Y))\Big{)}}\sqrt{2\sum_{i=1}^{n}\mu_{i}}\frac{2p}{1+\sqrt{1-p}}$
$\displaystyle=\sqrt{\mathrm{Var}\big{[}D_{\phi}({\tilde{Y}},g(Y))\Big{)}}\sqrt{2\sum_{i=1}^{n}\mu_{i}}\frac{5p}{3},$

which proves (25).

###  C.2 Proof of Proposition 5

We start from the fact that

$\mathrm{Var}[\mathrm{CB}_{p}\,|\,Y]=\frac{1}{B}\mathrm{Var}\Big{[}D_{\phi}(Y^{\dagger},g(Y^{*}))+\phi(Y^{*})-\phi(Y^{\dagger})\,\big{|}\,Y\Big{]}.$

Then

$\displaystyle\mathbb{E}\big{[}\mathrm{Var}[\mathrm{CB}_{p}\,|\,Y]\big{]}$
$\displaystyle=\frac{1}{B}\mathbb{E}\bigg{[}\mathrm{Var}\Big{[}D_{\phi}(Y^{\dagger},g(Y^{*}))+\phi(Y^{*})-\phi(Y^{\dagger})\,\big{|}\,Y\Big{]}\bigg{]}$
$\displaystyle=\frac{1}{B}\mathbb{E}\bigg{[}\mathrm{Var}\Big{[}\phi(Y^{*})-\phi(g(Y^{*}))-\langle\nabla\phi(g(Y^{*})),Y^{\dagger}-g(Y^{*})\rangle\,\big{|}\,Y\Big{]}\bigg{]}$
$\displaystyle=\frac{1}{B}\mathbb{E}\bigg{[}\mathrm{Var}\Big{[}\phi(Y^{*})-\phi(g(Y^{*}))-\langle\nabla\phi(g(Y^{*})),Y^{\dagger}-Y^{*}+Y^{*}-g(Y^{*})\rangle\,\big{|}\,Y\Big{]}\bigg{]}$
$\displaystyle=\frac{1}{B}\mathbb{E}\bigg{[}\mathrm{Var}\Big{[}D_{\phi}(Y^{*},g(Y^{*}))-\langle\nabla\phi(g(Y^{*})),Y^{\dagger}-Y^{*}\rangle\,\big{|}\,Y\Big{]}\bigg{]}$
$\displaystyle=\frac{1}{B}\mathbb{E}\bigg{[}\mathrm{Var}\Big{[}D_{\phi}(Y^{*},g(Y^{*}))-\frac{1}{p}\langle\omega,\nabla\phi(g(Y^{*}))\rangle+\langle\nabla\phi(g(Y^{*})),Y\rangle\,\big{|}\,Y\Big{]}\bigg{]}$
$\displaystyle\leq\frac{2}{B}\mathbb{E}\bigg{[}\mathrm{Var}\Big{[}D_{\phi}(Y^{*},g(Y^{*}))+\langle\nabla\phi(g(Y^{*})),Y\rangle\,\big{|}\,Y\Big{]}\bigg{]}+\frac{2}{Bp^{2}}\mathbb{E}\bigg{[}\mathrm{Var}\Big{[}\langle\omega,\nabla\phi(g(Y^{*}))\rangle\,\big{|}\,Y\Big{]}\bigg{]}$
$\displaystyle\leq\frac{2}{B}\mathrm{Var}\Big{[}D_{\phi}(Y^{*},g(Y^{*}))+\langle
Y^{*},\nabla\phi(g(Y^{*}))\rangle\Big{]}+\frac{2}{Bp^{2}}\mathrm{Var}\Big{[}\langle\omega,\nabla\phi(g(Y^{*}))\rangle\Big{]}$
$\displaystyle=\frac{2}{B}\mathrm{Var}\Big{[}D_{\phi}(Y,g(Y))+\langle
Y,\nabla\phi(g(Y))\rangle\Big{]}+\frac{2}{B}O(p)+\frac{2}{Bp^{2}}\mathrm{Var}\big{[}\langle\omega,\nabla\phi(g(Y^{*}))\rangle\big{]},$

where the second-to-last line uses the law of total variance, and the last
line uses a Taylor expansion which follows from the continuity in $p$ of
expected value of functions of $Y_{p}$ as assumed in the proof of Proposition
3. For the last term in the above display, note that $\omega$ and $Y^{*}$ are
independent. Therefore, we can use that fact that if $X_{1}$ and $X_{2}$ are
independent, then
$\mathrm{Var}[X_{1}X_{2}]=\mathrm{Var}[X_{1}]\mathrm{Var}[X_{2}]+\mathrm{Var}[X_{1}]\mathbb{E}^{2}[X_{2}]+\mathrm{Var}[X_{2}]\mathbb{E}^{2}[X_{1}]$,
which translates to

$\displaystyle\frac{2}{Bp^{2}}\mathrm{Var}\big{[}$
$\displaystyle\langle\omega,\nabla\phi(g(Y^{*}))\rangle\big{]}$
$\displaystyle=\frac{2}{Bp^{2}}\sum_{i=1}^{n}\mathrm{Var}\big{[}\omega_{i}\nabla_{i}\phi(g(Y^{*}))\big{]}$
$\displaystyle=\frac{2}{Bp^{2}}\sum_{i=1}^{n}\Big{(}\mathrm{Var}[\omega_{i}]\mathrm{Var}\big{[}\nabla_{i}\phi(g(Y^{*}))\big{]}+\mathrm{Var}[\omega_{i}]\mathbb{E}^{2}\big{[}\nabla_{i}\phi(g(Y^{*}))\big{]}+\mathrm{Var}\big{[}\nabla_{i}\phi(g(Y^{*}))\big{]}\mathbb{E}^{2}[\omega_{i}]\Big{)}$
$\displaystyle=\frac{2}{Bp^{2}}\sum_{i=1}^{n}\Big{(}p\mu_{i}\mathrm{Var}\big{[}\nabla_{i}\phi(g(Y^{*}))\big{]}+p\mu_{i}\mathbb{E}^{2}\big{[}\nabla_{i}\phi(g(Y^{*}))\big{]}+\mathrm{Var}\big{[}\nabla_{i}\phi(g(Y^{*}))\big{]}p^{2}\mu_{i}^{2}\Big{)}$
$\displaystyle=\frac{2}{Bp}\sum_{i=1}^{n}\mu_{i}\mathbb{E}\big{[}\nabla_{i}\phi(g(Y^{*}))^{2}\big{]}+\frac{2}{B}\sum_{i=1}^{n}\mu_{i}^{2}\mathrm{Var}\big{[}\nabla_{i}\phi(g(Y^{*}))\big{]}$
$\displaystyle=\frac{2}{Bp}\sum_{i=1}^{n}\mu_{i}\mathbb{E}\big{[}\nabla_{i}\phi(g(Y^{*}))^{2}\big{]}.$

Putting it all together gives the desired result (26).

###  C.3 Proof of Proposition 6

Note that the irreducible variance
$\mathrm{Var}(\mathbb{E}[\mathrm{CB}_{p}(g)\,|\,Y])$ does not depend on $B$
(because the inner expectation does not), so we assume without a loss of
generality that $B=1$ henceforth. Observe that

$\displaystyle\mathrm{Var}\big{(}$
$\displaystyle\mathbb{E}[\mathrm{CB}_{p}(g)\,|\,Y]\big{)}$
$\displaystyle=\mathrm{Var}\Big{(}\mathbb{E}\Big{[}D_{\phi}(Y^{*},g(Y^{*}))+\langle\nabla\phi(g(Y^{*})),Y^{*}-Y^{\dagger}\rangle\,\big{|}\,Y\Big{]}\Big{)}$
$\displaystyle=\mathrm{Var}\bigg{(}\mathbb{E}\bigg{[}D_{\phi}(Y^{*},g(Y^{*}))+\langle\nabla\phi(g(Y^{*})),Y^{*}\rangle-\frac{1-p}{p}\langle\nabla\phi(g(Y^{*})),\omega\rangle\,\Big{|}\,Y\bigg{]}\bigg{)}$
$\displaystyle\leq
2\mathrm{Var}\Big{(}\mathbb{E}\Big{[}D_{\phi}(Y^{*},g(Y^{*}))+\langle\nabla\phi(g(Y^{*})),Y^{*}\rangle\big{|}\,Y\Big{]}\Big{)}+2\mathrm{Var}\bigg{(}\mathbb{E}\bigg{[}\frac{1-p}{p}\langle\nabla\phi(g(Y^{*})),\omega\rangle\Big{|}\,Y\bigg{]}\bigg{)}.$
(30)

For the first term in (30), we apply the law of total variance and

$\displaystyle\mathrm{Var}\Big{(}\mathbb{E}\Big{[}D_{\phi}(Y^{*},g(Y^{*}))+\langle\nabla\phi(g(Y^{*})),Y^{*}\rangle$
$\displaystyle\leq\mathrm{Var}\Big{[}D_{\phi}(Y^{*},g(Y^{*}))+\langle\nabla\phi(g(Y^{*})),Y^{*}\rangle\Big{]}$
$\displaystyle\to\mathrm{Var}\Big{[}D_{\phi}(Y,g(Y))+\langle\nabla\phi(g(Y)),Y\rangle\Big{]},\quad\text{as
$p\to 0$},$

where the last convergence is guaranteed provided that $D^{2}_{\phi}(Y,g(Y))$
and $\langle\nabla\phi(g(Y)),Y\rangle^{2}$ have finite expectation, according
to the first step in the proof of Proposition 3.

For the second term in (30) recalling that $\Phi_{g}$ as defined in the
statement of the proposition,

$\displaystyle\lim_{p\to
0}\mathrm{Var}\bigg{(}\mathbb{E}\bigg{[}\frac{1-p}{p}\langle\nabla\phi(g(Y^{*})),\omega\rangle\,\Big{|}\,Y\bigg{]}\bigg{)}$
$\displaystyle\leq\lim_{p\to
0}\mathbb{E}\bigg{(}\mathbb{E}^{2}\bigg{[}\frac{1-p}{p}\langle\nabla\phi(g(Y^{*})),\omega\rangle\,\Big{|}\,Y\bigg{]}\bigg{)}$
$\displaystyle\leq\lim_{p\to
0}\mathbb{E}\bigg{(}\mathbb{E}^{2}\bigg{[}\frac{1-p}{p}\langle\Phi_{g}(Y),\omega\,\Big{|}\,Y\bigg{]}\bigg{)}$
$\displaystyle=\lim_{p\to
0}\,(1-p)^{2}\mathbb{E}\Big{[}\langle\Phi_{g}(Y),\mathbb{E}[\omega/p\,|\,Y]\rangle^{2}\Big{]}$
$\displaystyle=\lim_{p\to
0}\,(1-p)^{2}\mathbb{E}\big{[}\langle\Phi_{g}(Y),Y\rangle^{2}\big{]}$
$\displaystyle=\mathbb{E}\big{[}\langle\Phi_{g}(Y),Y\rangle^{2}\big{]}.$

Putting it all together gives the desired result (27).