# A General Framework for Empirical Bayes Estimation in Discrete Linear
Exponential Family

Trambak Banerjee<EMAIL_ADDRESS>
Data Sciences and Operations
University of Southern California
Los Angeles, CA 90089, USA Qiang Liu<EMAIL_ADDRESS>
Computer Science
University of Texas at Austin
Austin, Texas 78712, USA Gourab Mukherjee<EMAIL_ADDRESS>
Data Sciences and Operations
University of Southern California
Los Angeles, CA 90089, USA Wenguang Sun<EMAIL_ADDRESS>
Data Sciences and Operations
University of Southern California
Los Angeles, CA 90089, USA

###### Abstract

We develop a Nonparametric Empirical Bayes (NEB) framework for compound
estimation in the discrete linear exponential family, which includes a wide
class of discrete distributions frequently arising from modern big data
applications. We propose to directly estimate the Bayes shrinkage factor in
the generalized Robbins’ formula via solving a scalable convex program, which
is carefully developed based on a RKHS representation of the Stein’s
discrepancy measure. The new NEB estimation framework is flexible for
incorporating various structural constraints into the data driven rule, and
provides a unified approach to compound estimation with both regular and
scaled squared error losses. We develop theory to show that the class of NEB
estimators enjoys strong asymptotic properties. Comprehensive simulation
studies as well as analyses of real data examples are carried out to
demonstrate the superiority of the NEB estimator over competing methods.

Keywords: Asymptotic Optimality; Empirical Bayes; Power Series Distributions;
Shrinkage estimation; Stein’s discrepancy

##  1 Introduction

Shrinkage methods, exemplified by the seminal work of James and Stein (1961),
have received renewed attention in modern large-scale inference problems
(Efron, 2012; Fourdrinier et al., 2018). Under this setting, the classical
Normal means problem has been extensively studied (Brown, 2008; Jiang et al.,
2009; Brown and Greenshtein, 2009; Efron, 2011; Xie et al., 2012; Weinstein et
al., 2018). However, in a variety of applications, the observed data are often
discrete. For instance, in the News Popularity study discussed in Section 5,
the goal is to estimate the popularity of a large number of news items based
on their frequencies of being shared in social media platforms such as
Facebook and LinkedIn. Another important application scenario arises from
genomics research, where estimating the expected number of mutations across a
large number of genomic locations can help identify key drivers or inhibitors
of a given phenotype of interest.

We mention two main limitations of existing shrinkage estimation methods.
First, the methodology and theory developed for continuous variables, in
particular for Normal means problem, may not be directly applicable to
discrete models. Second, existing methods have focused on the squared error
loss. However, the scaled loss (Clevenson and Zidek, 1975), which effectively
reflects the asymmetries in decision making [cf. Equation (3)], is a more
desirable choice for many discrete models such as Poisson, where the scaled
loss corresponds to the local Kulback-Leibler distance. The scaled loss also
provides a more desirable criterion in a range of sparse settings, for
example, when the goal is to estimate the rates of rare outcomes in Binomial
distributions (Fourdrinier and Robert, 1995). Much research is needed for
discrete estimation problems under various loss functions. This article
develops a general framework for empirical Bayes estimation for the discrete
linear exponential (DLE) family, also known as the family of discrete power
series distributions (Noack, 1950), under both regular and scaled error
losses.

The DLE family includes a wide class of popular members such as the Poisson,
Binomial, negative Binomial and Geometric distributions. Let $Y$ be a non-
negative integer valued random variable. Then $Y$ is said to belong to a DLE
family if its probability mass function (pmf) is of the form

$p(y|\theta)=\dfrac{a_{y}\theta^{y}}{g(\theta)},~{}~{}y\in\\{0,1,2,\cdots\\},$
(1)

where $a_{y}$ and $g(\theta)$ are known functions such that $a_{y}\geq 0$ is
independent of $\theta$ and $g(\theta)$ is a normalizing factor that is
differentiable at every $\theta$. Special cases of DLE include the
$\mbox{Poisson}(\lambda)$ distribution with $a_{y}=(y!)^{-1}$,
$\theta=\lambda$ and $g(\theta)=\exp{(\theta)}$, and the
$\mbox{Binomial}(m,q)$ distribution with $a_{y}=\binom{m}{y}$,
$\theta=q/(1-q)$ and $g(\theta)=(1+\theta)^{m}$.

Suppose $Y_{1},\ldots,Y_{n}$ obey the following hierarchical model

$Y_{i}~{}|~{}\theta_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\texttt{DLE}(\theta_{i}),~{}~{}\theta_{i}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}G(\cdot),$ (2)

where $G(\cdot)$ is an unspecified prior distributusuaion on $\theta_{i}$. The
problem of interest is to estimate
$\bm{\theta}=(\theta_{1},\ldots,\theta_{n})$ based on
$\bm{Y}=(Y_{1},\ldots,Y_{n})$. Empirical Bayes approaches to this compound
decision problem date back to the famous Robbins’ formula (Robbins, 1956)
under the Poisson model. Important recent progresses by Brown et al. (2013),
Koenker and Mizera (2014) and Koenker and Gu (2017) show that Robbins’
estimator can be vastly improved by incorporating smoothness and monotonicity
adjustments. The main idea of existing works is to approximate the shrinkage
factor in the Bayes estimator as smooth functionals of the unknown marginal
pmf $p(y)$. The pmf can be estimated in various ways including the observed
empirical frequencies (Robbins, 1956), the smoothness-adjusted estimator
(Brown et al., 2013) or the shape-constrained NPMLE approach (Koenker and
Mizera, 2014; Koenker and Gu, 2017).

This article develops a general non-parametric empirical Bayes (NEB) framework
for compound estimation in discrete models. We first derive generalized
Robbins’ formula (GRF) for the DLE Model (2), and then implement GRF via
solving a scalable convex program. The powerful convex program, which is
carefully developed based on a reproducing kernel Hilbert space (RKHS)
representation of Stein’s discrepancy measure, leads to a class of efficient
NEB shrinkage estimators. We develop theories to show that the NEB estimator
is $\sqrt{n}$ consistent up to certain logarithmic factors and enjoys superior
risk properties. Simulation studies are conducted to illustrate the
superiority of the proposed NEB estimator when compared to existing state-of-
the-art approaches such as Brown et al. (2013), Koenker and Mizera (2014) and
Koenker and Gu (2017). We show that the NEB estimator has smaller risk in all
comparisons and the efficiency gain is substantial in many settings.

There are several advantages of the proposed NEB estimation framework. First,
in contrast with existing methods such as the smoothness-adjusted Poisson
estimator in Brown et al. (2013), our methodology covers a much wider range of
distributions and presents a unified approach to compound estimation in
discrete models. Second, our proposed convex program is fast and scalable. It
directly produces stable estimates of optimal Bayes shrinkage factors and can
easily incorporate various structural constraints into the decision rule. By
contrast, the three-step estimator in Brown et al. (2013), which involves
smoothing, Rao-Blackwellization and monotonicity adjustments, is complicated,
computationally intensive and sometimes unstable (as the numerator and
denominator of the ratio are computed separately). Third, the RKHS
representation of Stein’s discrepancy measure provides a new analytical tool
for developing theories such as asymptotic optimality and convergence rates.
Finally, the NEB estimation framework is robust to departures from the true
model due to its utilization of a generic quadratic program that does not rely
on the specific form of a particular DLE family. Our numerical results in
Section 4 demonstrate that the NEB estimator has significantly better risk
performance than competitive approaches of Efron (2011), Brown et al. (2013)
and Koenker and Gu (2017) under a mis-specified Poisson model.

An alternative approach to compound estimation in discrete models, as
suggested and investigated by Brown et al. (2013), is to employ variance
stabilizing transformations, which converts the discrete problem to a
classical normal means problem. This allows estimation via Tweedie’s formula
for normal variables (Efron, 2011), where the marginal density can be
estimated using NPMLE (Jiang et al., 2009; Koenker and Mizera, 2014) or
through kernel density methods (Brown and Greenshtein, 2009). However, there
are several drawbacks of this approach compared to our NEB framework. First,
Tweedie’s formula is not applicable to scaled error loss whereas our
methodology is built upon the generalized Robbins’ formula, which covers both
regular and scaled squared error losses. Second, there can be information loss
in conventional data processing steps such as standardization, transformation
and continuity approximation. While investigating the impact of information
loss on compound estimation is of great interest, it is desirable to develop
methodologies directly based on generalized Robbins’ formula that is
specifically derived and tailored for discrete variables. Finally, our NEB
framework provides a convenient tool for developing asymptotic theories. By
contrast, convergence rates are yet to be developed for normality inducing
transformations, which can be highly non-trivial.

The rest of the paper is organized as follows. In Section 2, we introduce our
estimation framework while Section 3 presents a theoretical analysis of the
NEB estimator. The numerical performance of our method is investigated using
both simulated and real data in Sections 4 and 5 respectively. Additional
technical details and proofs are relegated to the Appendices.

##  2 A General Framework for Compound Estimation in DLE Family

This section describes the proposed NEB framework for compound estimation in
discrete models. We first introduce in Section 2.1 the generalized Robbins’
formula for the DLE family (2), then propose in Section 2.2 a convex
optimization approach for its practical implementation. Details for tuning
parameter selection are discussed in Section 2.3.

###  2.1 Generalized Robbins’ formula for DLE models

Denote $\delta_{i}$ an estimator of $\theta_{i}$. Consider a class of loss
functions

$\ell^{(k)}(\theta_{i},\delta_{i})=\theta_{i}^{-k}(\theta_{i}-\delta_{i})^{2}$
(3)

for $k\in\\{0,1\\}$, where $\ell^{(0)}(\theta_{i},\delta_{i})$ is the usual
squared error loss, and
$\ell^{(1)}(\theta_{i},\delta_{i})=\theta_{i}^{-1}(\delta_{i}-\theta_{i})^{2}$
corresponds to the scaled squared error loss (Clevenson and Zidek, 1975;
Fourdrinier and Robert, 1995). In compound estimation, one is concerned with
the average loss

$\mathcal{L}_{n}^{(k)}(\bm{\theta},\bm{\delta})=n^{-1}\sum_{i=1}^{n}\ell^{(k)}(\theta_{i},\delta_{i}).$

The associated risk is denoted
$R_{n}^{(k)}(\bm{\theta},\bm{\delta})=\mathbb{E_{\boldsymbol{Y}|\boldsymbol{\theta}}}\mathcal{L}_{n}^{(k)}(\bm{\theta},\bm{\delta})$.
Let $\boldsymbol{G}(\boldsymbol{\theta})$ denote the joint distribution of
$(\theta_{1},\cdots,\theta_{n})$. The Bayes estimator
$\bm{\delta}^{\pi}_{(k)}$ that minimizes the Bayes risk
$B_{n}^{(k)}(\bm{\theta})=\int
R_{n}^{(k)}(\bm{\theta},\bm{\delta})d\boldsymbol{G}(\bm{\theta})$ is given by
Lemma 1.

######  Lemma 1 (Generalized Robbins’ formula).

Consider the DLE Model (2). Let $p(\cdot)=\int p(\cdot|\theta)dG(\theta)$ be
the marginal pmf of $Y$. Define for $k\in\\{0,1\\}$,

$w_{p}^{(k)}(y_{i})=\dfrac{p(y_{i}-k)}{p(y_{i}+1-k)},\text{ for
}y_{i}=k,k+1,\cdots.$

Then the Bayes estimator that minimizes the risk $B_{n}^{(k)}(\bm{\theta})$ is
given by $\bm{\delta}^{\pi}_{(k)}=\\{\delta_{(k),i}^{\pi}(y_{i}):1\leq i\leq
n\\}$, where

$\delta_{(k),i}^{\pi}(y_{i})=\begin{cases}\dfrac{a_{y_{i}-k}/a_{y_{i}+1-k}}{w_{p}^{(k)}(y_{i})},\text{
for }y_{i}=k,k+1,\cdots\\\ 0,\quad\quad\quad\quad\quad\quad\text{ for
}y_{i}<k\end{cases}.$ (4)

######  Remark 1.

Under the squared error loss $(k=0)$ with
$Y_{i}~{}|~{}\theta_{i}\sim\text{Poi}(\theta_{i})$ and
$a_{y_{i}}=(y_{i}!)^{-1}$, Lemma 1 yields

$\delta_{(0),i}^{\pi}(y_{i})=(y_{i}+1)\dfrac{p(y_{i}+1)}{p(y_{i})},$ (5)

which recovers the classical Robbins’ formula (Robbins, 1956). In contrast,
under the scaled loss, we have

$\delta_{(1),i}^{\pi}(y_{i})=y_{i}\dfrac{p(y_{i})}{p(y_{i}-1)}\mbox{ for
}y_{i}>0\mbox{ and }\delta_{(1),i}^{\pi}(y_{i})=0\mbox{ otherwise}.$ (6)

Under scaled error loss the estimator (5) can be much outperformed by (6) (and
vice versa under the regular loss). We develop parallel results for the two
types of loss functions.

Next we discuss related works for implementing Robbins’ formula under the
empirical Bayes (EB) estimation framework. Inspecting (4) and (5), we can view
${a_{y_{i}-k}/a_{y_{i}+1-k}}$ as a naive and known estimator of $\theta_{i}$.
The ratio functional $w_{p}^{(k)}(y_{i})$, which is unknown in practice,
represents the optimal shrinkage factor that depends on $p(\cdot)$. Hence, a
simple EB approach, as done in the classical Robbins’ formula, is to estimate
$w^{(k)}_{p}(y)$ by plugging-in empirical frequencies:
$\hat{w}_{n}^{(0)}(y)=\hat{p}_{n}(y)/\hat{p}_{n}(y+1)$, where
$\hat{p}_{n}(y)=n^{-1}\sum_{i=1}^{n}\mathbb{I}(y_{i}=y)$. It is noted by Brown
et al. (2013) that this plug-in estimator can be highly inefficient especially
when $\theta_{i}$ are small. Moreover, the numerator and denominator in
$w_{p}^{(0)}(y)$ are estimated separately, which may lead to unstable ratios.
Brown et al. (2013) showed that Robbins’ formula can be dramatically improved
by imposing additional smoothness and monotonicity adjustments. An alternative
approach is to estimate $p(y)$ using NPMLE (Jiang et al., 2009) under
appropriate shape constraints (Koenker and Mizera, 2014). However, efficient
estimation of $p(y)$ may not directly translate into an efficient estimation
of the underlying ratio $p(y+1)/p(y)$. We recast the compound estimation
problems as a convex program, which directly produces consistent estimates of
the ratio functionals

$\bm{w}_{p}^{(k)}=\left\\{w_{p}^{(k)}(y_{1}),\ldots,w_{p}^{(k)}(y_{n})\right\\}$

from data. The estimators are shown to enjoy superior numerical and
theoretical properties. Unlike existing works that are limited to regular loss
and specific members in the DLE family, our method can handle a wide range of
distributions and various types of loss functions in a unified framework.

###  2.2 Shrinkage estimation by convex optimization

This section focuses on the scaled squared error loss $(k=1)$. Methodologies
and theories for the case with usual squared error loss $(k=0)$ can be derived
similarly; details are provided in Appendix A.1. We first introduce some
notations and then present the NEB estimator in Definition 1A.

Suppose $Y$ is a non-negative integer-valued random variable with pmf
$p(\cdot)$. Define

$h_{0}^{(1)}(y)=\begin{cases}1~{},\quad\quad\quad\quad\quad\text{if~{}}y=0\\\
1-w_{p}^{(1)}(y),\quad\text{if~{}}y\in\\{1,2,\ldots\\}.\end{cases}$ (7)

Let
$\mathcal{K}_{\lambda}(y,y^{\prime})=\exp\\{-\frac{1}{2\lambda}(y-y^{\prime})^{2}\\}$
be the positive definite Radial Basis Function (RBF) kernel with bandwidth
parameter $\lambda\in\Lambda$ where $\Lambda$ is a compact subset of
$\mathbb{R}^{+}$ bounded away from $0$. Given observations
$\bm{y}=(y_{1},\ldots,y_{n})$ from Model (2), let
$\bm{h}_{0}^{(1)}=\left\\{h_{0}^{(1)}(y_{1}),\ldots,h_{0}^{(1)}(y_{n})\right\\}$.
Define operators
$\Delta_{y}\mathcal{K}_{\lambda}(y,y^{\prime})=\mathcal{K}_{\lambda}(y+1,y^{\prime})-\mathcal{K}_{\lambda}(y,y^{\prime})$
and

$\Delta_{y,y^{\prime}}\mathcal{K}_{\lambda}(y,y^{\prime})=\Delta_{y^{\prime}}\Delta_{y}\mathcal{K}_{\lambda}(y,y^{\prime})=\Delta_{y}\Delta_{y^{\prime}}\mathcal{K}_{\lambda}(y,y^{\prime}).$

Consider the following $n\times n$ matrices, which are needed in the
definition of the NEB estimator:

$\bm{K}_{\lambda}={n^{-2}}[\mathcal{K}_{\lambda}(y_{i},y_{j})]_{ij},\quad{\Delta\bm{K}}_{\lambda}={n^{-2}}[\Delta_{y_{i}}\mathcal{K}_{\lambda}(y_{i},y_{j})]_{ij},\quad{\Delta_{2}\bm{K}}_{\lambda}={n^{-2}}[\Delta_{y_{i},y_{j}}\mathcal{K}_{\lambda}(y_{i},y_{j})]_{ij}.$

######  Definition 1A (NEB estimator).

Consider the DLE Model (2) with loss $\ell^{(1)}(\theta_{i},\delta_{i})$. For
any fixed $\lambda\in\Lambda$, let
$\hat{\bm{h}}_{n}^{(1)}(\lambda)=\left\\{\hat{h}_{1}^{(1)}(\lambda),\ldots,\hat{h}_{n}^{(1)}(\lambda)\right\\}$
be the solution to the following quadratic optimization problem:

$\displaystyle\min_{\bm{h}\in\bm{H}_{n}}~{}\hat{\mathbb{M}}_{\lambda,n}(\bm{h})=\bm{h}^{T}\bm{K}_{\lambda}\bm{h}+2\bm{h}^{T}\Delta\bm{K}_{\lambda}\bm{1}+\bm{1}^{T}\Delta_{2}\bm{K}_{\lambda}\bm{1},$
(8)

where
$\bm{H}_{n}=\\{\bm{h}=(h_{1},\ldots,h_{n}):\mathcal{A}\bm{h}\preceq\bm{b},~{}\mathcal{C}\bm{h}=\bm{d}\\}$
is a convex set and $\mathcal{A},\mathcal{C},\bm{b}$ and $\bm{d}$ are known
real matrices and vectors that enforce linear constraints on the components of
$\bm{h}$. Define $\hat{w}_{i}^{(1)}(\lambda)=1-\hat{h}_{i}^{(1)}(\lambda)$.
Then the NEB estimator is given by ${\bm{\delta}}^{\sf
neb}_{(1)}(\lambda)=\left\\{{\delta}_{(1),i}^{\sf neb}(\lambda):1\leq i\leq
n\right\\}$, where

${\delta}_{(1),i}^{\sf
neb}(\lambda)=\dfrac{a_{y_{i}-1}/a_{y_{i}}}{\hat{w}_{i}^{(1)}(\lambda)},\quad\mbox{
if }y_{i}\in\\{1,2,\ldots\\},$

and ${\delta}_{(1),i}^{\sf neb}(\lambda)=0$ if $y_{i}=0$.

Next we provide some insights on why the optimization criterion (8) works;
theories are developed in Section 3 to establish the properties of the NEB
estimator rigorously. Denote $h_{0}^{(1)}$ and $\tilde{h}^{(1)}$ as the ratio
functionals corresponding to pmfs $p$ and $\tilde{p}$, respectively. Suppose
$Y_{i}$ are i.i.d. samples obeying $p(y)$. Theorem 1 shows that

$\hat{\mathbb{M}}_{\lambda,n}(\tilde{\bm{h}})={\mathbb{M}}_{\lambda}(\tilde{\bm{h}})+O_{p}\Big{(}\dfrac{\log^{2}n}{n^{1/2}}\Big{)},$

where $\hat{\mathbb{M}}_{\lambda,n}(\tilde{\bm{h}})$ is the objective function
in (8) and ${\mathbb{M}}_{\lambda}(\tilde{\bm{h}})$, also denoted
$\mathcal{S}_{\lambda}[\tilde{p}](p)$, is the kernelized Stein’s discrepancy
(KSD). Roughly speaking, the KSD measures how different one distribution $p$
is from another distribution $\tilde{p}$, with
$\mathcal{S}_{\lambda}[\tilde{p}](p)=0$ if and only if $\tilde{p}=p$. A key
feature of the KSD is that $\mathcal{S}_{\lambda}[\tilde{p}](p)$ can be
equivalently represented by the discrepancy between the corresponding ratio
functionals $h_{0}^{(1)}$ and $\tilde{h}^{(1)}$. Hence, optimizing (8) is
asymptotically equivalent to finding $\tilde{h}^{(1)}$ that is as close as
possible to the true underlying $h_{0}^{(1)}$, which corresponds to the
optimal shrinkage factor in the compound estimation problem. Theorems 2A and
3A demonstrate that (8) is an effective convex program in the sense that the
minimizer $\hat{\bm{h}}_{n}$ is $\sqrt{n}$ consistent with respect to
$\bm{h}_{0}^{(1)}$, and the resultant NEB estimator converges to the Bayes
estimator.

###  2.3 Structural constraints and bandwidth selection

In problem (8) the linear inequality $\mathcal{A}\bm{h}\preceq\bm{b}$ can be
used to impose structural constraints on the NEB rule ${\bm{\delta}}^{\sf
neb}_{(1)}(\lambda)$. The structural constraints, which may take the form of
monotonicity constraints as pursued in, for example, Brown et al. (2013) and
Koenker and Mizera (2014), have been shown to be effective for stabilizing the
estimator and hence improving the accuracy. For example, when
$Y_{i}~{}|~{}\theta_{i}~{}\sim~{}\text{Poi}(\theta_{i})$ then
$\delta_{(1),i}^{\pi}(y_{i})=(y_{i}+1)/w_{p}^{(1)}(y_{i})$ and $\mathcal{A}$,
$\bm{b}$ can be chosen such that

$h_{(i-1)}-\frac{y_{(i-1)}+1}{y_{(i)}+1}h_{(i)}\leq
1-\frac{y_{(i-1)}+1}{y_{(i)}+1},\text{ for }2\leq i\leq n$

and $y_{(1)}\leq y_{(2)}\leq\cdots\leq y_{(n)}$. Moreover, when $y_{i}=0$ we
set $\delta_{(1),i}^{\sf neb}(\lambda)=0$ by convention (see lemma 1). The
equality constraints $\mathcal{C}\bm{h}=\bm{d}$ accommodate such boundary
conditions along with instances of ties for which we require $h_{i}=h_{j}$
whenever $y_{i}=y_{j}$.

The implementation of the quadratic program in (8) requires the choice of a
tuning parameter $\lambda$ in the RBF kernel. For practical applications,
$\lambda$ must be determined in a data-driven fashion. For infinitely
divisible random variables (Klenke, 2014) such as Poisson variables, Brown et
al. (2013) proposed a modified cross validation (MCV) method for choosing the
tuning parameter. However, the MCV method cannot be applied to distributions
with bounded support, e.g. variables that are not infinitely divisible (Sato
and Ken-Iti, 1999) such as the Binomial distribution. To provide a unified
estimation framework for the DLE family, we develop an alternative method for
choosing $\lambda$. The key idea is to derive an asymptotic risk estimate
${\sf ARE}^{(1)}_{n}(\lambda)$ that serves as an approximation of the true
risk $\mathcal{R}_{n}^{(1)}(\bm{\theta},\cdot)$. Then the tuning parameter is
chosen to minimize ${\sf ARE}^{(1)}_{n}(\lambda)$.

The methodology based on ARE is illustrated below for Poisson and Binomial
models under the scaled loss (see Definitions 2A and 3A, respectively). The
ideas can be extended to other members in the DLE family. In Appendix A.2, we
provide relevant details for choosing $\lambda$ under the regular loss
$\mathcal{L}_{n}^{(0)}$.

######  Definition 2A (ARE of ${\bm{\delta}}_{(1)}^{\sf neb}(\lambda)$ in the
Poisson model).

Suppose $Y_{i}~{}|~{}\theta_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\texttt{Poi}(\theta_{i})$. Under the loss
$\ell^{(1)}(\theta_{i},\cdot)$, an ARE of the true risk of
${\bm{\delta}}_{(1)}^{\sf neb}(\lambda)$ is

${\sf
ARE}^{(1,\mathcal{P})}_{n}(\lambda)=\dfrac{1}{n}\Bigl{\\{}\sum_{i=1}^{n}y_{i}+\sum_{i=1}^{n}{\psi}_{\lambda}(y_{i})-2\sum_{i=1}^{n}{\delta}_{(1),i}^{\sf
neb}(\lambda)\Bigr{\\}},\mbox{ where}$
$\psi_{\lambda}(y_{i})=\\{{\delta}_{(1),j}^{\sf
neb}(\lambda)\\}^{2}/(y_{i}+1),~{}y_{i}=0,1,\ldots.$

with $j\in\\{1,\ldots,n\\}$ such that $y_{j}=y_{i}+1$.

For the Binomial model, we proceed along similar lines and consider the
following asymptotic risk estimate of the true risk.

######  Definition 3A (ARE of ${\bm{\delta}}_{(1)}^{\sf neb}(\lambda)$ in the
Binomial model).

Suppose $Y_{i}~{}|~{}q_{i}\sim\mathcal{B}in(m,q_{i})$. Hence in Equations (1)
and (2) we have $a_{y_{i}}={m\choose y_{i}}$ and $\theta_{i}=q_{i}/(1-q_{i})$.
Under the loss $\ell^{(1)}(\theta_{i},\cdot)$, an ARE of the true risk of
${\bm{\delta}}_{(1)}^{\sf neb}(\lambda)$ is

${\sf
ARE}^{(1,\mathcal{B})}_{n}(\lambda)=\dfrac{1}{n}\Bigl{\\{}\sum_{i=1}^{n}\dfrac{y_{i}}{m-y_{i}+1}+\sum_{i=1}^{n}(m-y_{i}){\psi}_{\lambda}(y_{i})-2\sum_{i=1}^{n}{\delta}_{(1),i}^{\sf
neb}(\lambda)\Bigr{\\}},\mbox{ where}$
$\psi_{\lambda}(y_{i})=\\{{\delta}_{(1),j}^{\sf
neb}(\lambda)\\}^{2}/(y_{i}+1),~{}y_{i}=0,\ldots,m.\\\ $

with $j\in\\{1,\ldots,n\\}$ such that $y_{j}=y_{i}+1$.

######  Remark 2.

Although the expression for $\psi_{\lambda}$ appears to be identical across
Definitions 2A and 3A, it differs with respect to ${\bm{\delta}}_{(1)}^{\sf
neb}(\lambda)$. Specifically, in definition 2A, ${\bm{\delta}}_{(1)}^{\sf
neb}(\lambda)$ is the NEB estimator of the Poisson means, whereas in
Definition 3A, ${\bm{\delta}}_{(1)}^{\sf neb}(\lambda)$ is the NEB estimator
of the Binomial odds.

######  Remark 3.

If for some $i$, $y_{i}+1$ is not available in the observed sample $\bm{y}$,
${\psi}_{\lambda}(y_{i})$ can be calculated using cubic splines, and a linear
interpolation can be used to tackle the boundary point of the observed sample
maxima.

We propose the following estimate of the tuning parameter $\lambda$ based on
the ARE:

$\hat{\lambda}=\begin{cases}\operatorname*{arg\,min}_{\lambda\in{\Lambda}}{\sf
ARE}^{(1,\mathcal{P})}_{n}(\lambda),~{}\text{if~{}}Y_{i}~{}|~{}\theta_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\texttt{Poi}(\theta_{i})\\\
\operatorname*{arg\,min}_{\lambda\in{\Lambda}}{\sf
ARE}^{(1,\mathcal{B})}_{n}(\lambda),~{}\text{if~{}}Y_{i}~{}|~{}q_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\texttt{Bin}(m,q_{i})\end{cases}.$ (9)

In practice we recommend using $\Lambda=[10,10^{2}]$, which works well in all
our simulations and real data analyses. In Section 3, we present Lemmas 2 and
3 to provide asymptotic justifications for selecting $\lambda$ using equation
(9).

##  3 Theory

This section studies the theoretical properties for the NEB estimator under
the Poisson and Binomial models. We first investigate the large-sample
behavior of the KSD measure (Section 3.1), then turn to the performance of the
estimated risk ratios $\hat{\boldsymbol{w}}_{n}$ (Section 3.2), and finally
establish the consistency and risk properties of the proposed estimator
$\bm{\delta}^{\sf neb}$ (Section 3.3). The accuracy of the ${\sf ARE}$
criteria, which are used in choosing tuning parameter $\lambda$, will also be
investigated.

###  3.1 Theoretical properties of the KSD measure

To provide motivation and theoretical support for Definition 1A, we introduce
the Kernelized Stein’s Discrepancy (KSD) (Liu et al., 2016; Chwialkowski et
al., 2016) and discuss its connection to the quadratic program (8). While the
KSD has been used in various contexts including goodness of fit tests (Liu et
al., 2016), variational inference (Liu and Wang, 2016) and Monte Carlo
integration (Oates et al., 2017), our theory on its connection to the compound
estimation problem and empirical Bayes methodology is novel.

Assume that $Y$ and $Y^{\prime}$ are i.i.d. copies from the marginal pmf $p$.
Consider $h_{0}$ defined in Equation (7)111In Section 3.1 we shall drop the
superscript from $\bm{h}_{0}$, which is used to indicate whether the loss is
scaled or regular. The simplification has no impact since the general idea
holds for both types of losses and the discussion in this section focuses on
the scaled loss.. Let $\tilde{p}$ denote a pmf on the support of $Y$, for
which we similarly define $\tilde{h}$. The KSD, which is formally defined as

$\mathcal{S}_{\lambda}[\tilde{p}](p)=\mathbb{E}_{p}\Big{[}\left\\{\tilde{h}(Y)-h_{0}(Y)\right\\}\mathcal{K}_{\lambda}(Y,Y^{\prime})\left\\{\tilde{h}(Y^{\prime})-h_{0}(Y^{\prime})\right\\}\Big{]},$
(10)

provides a discrepancy measure between $p$ and $\tilde{p}$ in the sense that
(a)

$\mbox{$\mathcal{S}_{\lambda}[\tilde{p}](p)\geq 0$ and
$\mathcal{S}_{\lambda}[\tilde{p}](p)=0$ if and only if $p=\tilde{p}$},$

and (b) informally, $\mathcal{S}_{\lambda}[\tilde{p}](p)$ tends to increase
when there is a bigger disparity between $h_{0}$ and $\tilde{h}$ (or
equivalently, between $p$ and $\tilde{p}$).

The direct evaluation of $\mathcal{S}_{\lambda}[\tilde{p}](p)$ via Equation
(10) is difficult because $h_{0}$ is unknown. Note that while the pmf $p$ can
be learned well from a random sample $\\{Y_{1},\ldots,Y_{n}\\}\sim p$, we
introduce an alternative representation of KSD, developed by Liu et al.
(2016), in a reproducing kernel Hilbert space (RKHS) that does not directly
involve unknown $h_{0}$. Concretely, consider a smooth positive definite
kernel function $\kappa_{\lambda}[\tilde{p}]$:

$\kappa_{\lambda}[\tilde{p}](u,v)=\tilde{h}(u)\tilde{h}(v)\mathcal{K}_{\lambda}(u,v)+\tilde{h}(u)\Delta_{v}\mathcal{K}_{\lambda}(u,v)+\tilde{h}(v)\Delta_{u}\mathcal{K}_{\lambda}(u,v)+\Delta_{u,v}\mathcal{K}_{\lambda}(u,v).$
(11)

For i.i.d. copies $(Y,Y^{\prime})$ from distribution $p$, it can be shown that

$\displaystyle\mathcal{S}_{\lambda}[\tilde{p}](p)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{(Y,Y^{\prime})\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}p}\Big{[}\kappa_{\lambda}[\tilde{p}](Y,Y^{\prime})\Big{]}$
$\displaystyle=$
$\displaystyle\dfrac{1}{n(n-1)}\mathbb{E}_{p}\Big{[}\sum_{1\leq i\neq j\leq
n}\kappa_{\lambda}[\tilde{h}(Y_{i}),\tilde{h}(Y_{j})](Y_{i},Y_{j})\Big{]}$
$\displaystyle\coloneqq$
$\displaystyle{{\mathbb{M}}_{\lambda}(\tilde{\bm{h}}),}$

where $\\{Y_{1},\ldots,Y_{n}\\}$ is a random sample from $p$. It can be
similarly shown that $\mathbb{M}_{\lambda}(\tilde{\bm{h}})=0$ if and only if
$\tilde{\bm{h}}=\bm{h}_{0}$. Substituting the empirical distribution
$\hat{p}_{n}$ in place of the pmf $p$ in (3.1), we obtain the following
empirical evaluation scheme for $\mathcal{S}_{\lambda}[\tilde{p}](p)$ that is
both intuitive and computationally efficient:

$\mathcal{S}_{\lambda}[\tilde{p}](\hat{p}_{n})=\dfrac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\kappa_{\lambda}[{\tilde{h}(y_{i}),\tilde{h}(y_{j})}](y_{i},y_{j})\coloneqq\hat{\mathbb{M}}_{\lambda,n}(\tilde{\bm{h}}).$
(13)

Note that (13) is exactly the objective function of the quadratic program (8).

The empirical representation of KSD (13) provides an extremely useful tool for
solving the discrete compound decision problem under the EB estimation
framework. A key observation is that the kernel function
$\kappa_{\lambda}[\tilde{p}](u,v)$ depends on $\tilde{p}$ only through
$\tilde{h}$. Meanwhile, the EB implementation of the generalized Robbins’
formula [cf. Equations (4) and (7)] essentially boils down to the estimation
of $h_{0}$. Hence, if $\mathcal{S}_{\lambda}[\tilde{p}](\hat{p}_{n})$ is
asymptotically equal to $\mathcal{S}_{\lambda}[\tilde{p}](p)$, then minimizing
$\mathcal{S}_{\lambda}[\tilde{p}](\hat{p}_{n})$ with respect to the unknowns
$\tilde{\bm{h}}=\left\\{\tilde{h}(y_{1}),\ldots,\tilde{h}(y_{n})\right\\}$ is
effectively the process of finding an $\tilde{h}$ that is as close as possible
to $h_{0}$, which yields an asymptotically optimal solution to the EB
estimation problem. Therefore our formulation of the NEB estimator
$\bm{\delta}^{\sf neb}(\lambda)$ would be justified as long as we can
establish the asymptotic consistency of the sample criterion
$\mathcal{S}_{\lambda}[\tilde{p}](\hat{p}_{n})$ around the population
criterion $\mathcal{S}_{\lambda}[\tilde{p}](p)$ uniformly over $\lambda$
(Theorem 1).

For a fixed mass function $\tilde{p}$ on the support of $Y$, we impose the
following regularity conditions that are needed in our technical analysis.

  * (A1)

$\mathbb{E}_{p}|\kappa_{\lambda}[\tilde{h}(U),\tilde{h}(V)](U,V)|^{2}<\infty$
for all $\lambda\in\Lambda$ where $\Lambda$ is a compact subset of
$\mathbb{R}^{+}$ bounded away from $0$.

  * (A2)

For some $\epsilon\in(0,1)$,
$\varlimsup_{n\to\infty}n^{-1}\sum_{i=1}^{n}\exp(\epsilon\theta_{i})<\infty$.

  * (A3)

For any function $g$ that satisfies $0<\|g\|_{2}^{2}<\infty$, there exists a
constant $c>0$ such that
$\sum_{i,j=1}^{n}g(y_{i})\mathcal{K}_{\lambda}(y_{i},y_{j})g(y_{j})>c\|g\|_{2}^{2}$
for every $\lambda\in\Lambda$.

  * (A4)

The feasible solutions $\bm{h}_{n}$ to equations (8) and (15) satisfy
$\sup_{\bm{h}_{n}\in\bm{H}_{n}}\|\bm{h}_{n}\|_{1}=O(n\log n)$.

######  Remark 4.

Assumption (A1) is a standard moment condition on the kernel function related
to V-statistics, see, for example, Serfling (2009). Assumption (A2) ensures
that with high probability $\text{max}(Y_{1},\ldots,Y_{n})\leq\log n$ as
$n\to\infty$. This idea is formalized by Lemma A in Appendix B. Assumption
(A3) is a standard condition which ensures that the KSD
$\mathcal{S}_{\lambda}[\tilde{p}](p)$ is a valid discrepancy measure (Liu et
al., 2016; Chwialkowski et al., 2016). Assumption (A4) provides a control on
the growth rate of the $\ell_{1}$ norm of the feasible solutions. In
particular, both Assumptions (A3) and (A4) play a critical role in
establishing point-wise Lipschitz stability of the optimal solution
$\hat{\bm{h}}_{n}(\lambda)$ under perturbations on the tuning parameter
$\lambda\in\Lambda$ (see Lemma B in Appendix B).

######  Theorem 1.

If $p$ and $\tilde{p}$ are probability mass functions on the support of $Y$
then, under Assumptions (A1) and (A2), we have

$\sup_{\lambda\in\Lambda}\Big{|}\hat{\mathbb{M}}_{\lambda,n}(\tilde{\bm{h}})-{\mathbb{M}}_{\lambda}(\tilde{\bm{h}})\Big{|}=O_{p}\Big{(}\dfrac{\log^{2}n}{\sqrt{n}}\Big{)}.$

In the context of our compound estimation framework, Theorem 1 is significant
because it guarantees that the empirical version of the KSD measure given by
$\hat{\mathbb{M}}_{\lambda,n}(\tilde{\bm{h}})$ is asymptotically close to its
population counterpart $\mathbb{M}_{\lambda}(\tilde{\bm{h}})$ uniformly in
$\lambda\in\Lambda$. Moreover, along with the fact that
$\mathbb{M}_{\lambda}(\bm{h}_{0})=0$, Theorem 1 establishes that
$\hat{\mathbb{M}}_{\lambda,n}({\bm{h}})$ is the appropriate criteria to
minimize with respect to $\bm{h}\equiv\tilde{\bm{h}}$. In Theorem 2A, we
further show that the resulting estimator of the ratio functionals
$\bm{w}_{p}^{(1)}$ from equation (8) are consistent.

###  3.2 Theoretical properties of $\hat{\boldsymbol{w}}_{n}$

The optimization problem in (8) is defined over a convex set
$\bm{H}_{n}\subset\mathbb{R}^{n}$. However, the dimension of $\bm{H}_{n}$,
denoted by $\texttt{dim}(\bm{H}_{n})$, is usually much smaller than $n$.
Consider the Binomial case where $Y_{i}|q_{i}\sim\texttt{Bin}(m_{i},q_{i})$
with $q_{i}\in(0,1)$, $m_{i}\leq m<\infty$ and $\theta_{i}=q_{i}/(1-q_{i})$.
Here $\texttt{dim}(\bm{H}_{n})$ is at most $m$ since
$\text{max}(Y_{1},\ldots,Y_{n})\leq m$. While the boundedness of the support
is not always available outside the Binomial case, in most practical
applications it is reasonable to assume that the distribution of $\theta_{i}$
has some finite moments, which ensures that $\texttt{dim}(\bm{H}_{n})$ grows
slower than $\log n$; see Assumption (A2). In Lemma A we make this precise.
The next theorem establishes the asymptotic consistency of
$\hat{\bm{w}}_{n}^{(1)}(\lambda)$.

######  Theorem 2A.

Let $\mathcal{K}_{\lambda}(\cdot,\cdot)$ be the positive definite RBF kernel
with bandwidth parameter $\lambda\in\Lambda$. If
$\lim_{n\to\infty}c_{n}n^{-1/2}\log^{2}n=0$ then, under Assumptions (A1) -
(A3), we have for any $\lambda\in\Lambda$,

$\lim_{n\to\infty}\mathbb{P}\left\\{\Big{\|}\hat{\bm{w}}_{n}^{(1)}(\lambda)-{\bm{w}}_{p}^{(1)}\Big{\|}_{2}\geq
c_{n}^{-1}\epsilon\right\\}=0,~{}\text{for any }\epsilon>0,$

where $\hat{\bm{w}}_{n}^{(1)}(\lambda)=1-\hat{\bm{h}}_{n}^{(1)}(\lambda)$.

Theorem 2A shows that under the scaled squared error loss,
$\hat{\bm{w}}_{n}^{(1)}(\lambda)$, the optimizer of quadratic form (8),
provides a consistent estimator of ${\bm{w}}_{p}^{(1)}$, the optimal shrinkage
factor in the Bayes rule (Lemma 1). Theorem 2A is proved in Appendix B.3,
where we also include relevant details for proving a companion result under
the regular squared error loss.

######  Remark 5.

The estimation framework in Definition 1A may be used for producing consistent
estimators for any member in the DLE family. This allows the corresponding NEB
estimator to cover a much wider class of discrete distributions than
previously proposed. Compared to existing methods (Efron, 2011; Brown et al.,
2013; Koenker and Mizera, 2014; Koenker and Gu, 2017), our proposed NEB
estimation framework is robust against departures from the true data
generating process. This is due to the fact that the quadratic optimization
problem in (8) does not rely on the specific form of the distribution of
$Y|\theta$, and that the shrinkage factors are estimated in a non-parametric
fashion. The robustness of the estimator is corroborated by our numerical
results in Section 4.

###  3.3 Properties of the NEB estimator

In this section we discuss the risk properties of the NEB estimator. We begin
with two lemmas showing that uniformly in $\lambda\in\Lambda$, the gap between
the estimated risk ${\sf ARE}^{(1)}_{n}(\lambda)$ and true risk is
asymptotically negligible. This justifies our proposed methodology for
choosing the tuning parameter $\lambda$ in Section 2.3. In the following two
lemmas, we let $c_{n}$ be a sequence satisfying
$\lim_{n\to\infty}c_{n}n^{-1/2}\log^{5/2}n=0$.

######  Lemma 2.

Under Assumptions (A3) and (A4) and the Binomial model, we have

$\displaystyle(a).$ $\displaystyle c_{n}\sup_{\lambda\in\Lambda}\Big{|}{\sf
ARE}^{(1,\mathcal{B})}_{n}(\lambda,\bm{Y})-\mathcal{R}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))\Big{|}=o_{p}(1);$ $\displaystyle(b).$ $\displaystyle
c_{n}\sup_{\lambda\in\Lambda}\Big{|}{\sf
ARE}^{(1,\mathcal{B})}_{n}(\lambda,\bm{Y})-\mathcal{L}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))\Big{|}=o_{p}(1).$

######  Lemma 3.

Under Assumptions (A2), (A3) and (A4) and the Poisson model, we have

$\displaystyle(a).$ $\displaystyle c_{n}\sup_{\lambda\in\Lambda}\Big{|}{\sf
ARE}^{(1,\mathcal{P})}_{n}(\lambda,\bm{Y})-\mathcal{R}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))\Big{|}=o_{p}(1);$ $\displaystyle(b).$ $\displaystyle
c_{n}\sup_{\lambda\in\Lambda}\Big{|}{\sf
ARE}^{(1,\mathcal{P})}_{n}(\lambda,\bm{Y})-\mathcal{L}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))\Big{|}=o_{p}(1).$

To analyze the quality of the data-driven bandwidth $\hat{\lambda}$ [cf.
Equation (9)], we consider an oracle loss estimator $\bm{\delta}^{\sf
or}_{(1)}\coloneqq{\bm{\delta}_{(1)}^{\sf neb}}(\lambda^{\sf orc}_{1})$, where

$\lambda^{\sf
orc}_{1}\coloneqq\operatorname*{arg\,min}_{\lambda\in\Lambda}\mathcal{L}_{n}^{(1)}\left\\{\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda)\right\\}.$

The oracle bandwidth $\lambda^{\sf orc}_{1}$ is not available in practice
since it requires the knowledge of unknown $\boldsymbol{\theta}$. However, it
provides a benchmark for assessing the effectiveness of the data-driven
bandwidth selection procedure in Section 2.3. The following lemma shows that
the loss of ${\bm{\delta}}^{\sf neb}_{(1)}(\hat{\lambda})$ converges in
probability to the loss of $\bm{\delta}^{\sf or}_{(1)}$.

######  Lemma 4.

Under Assumptions (A2) - (A4), if
$\lim_{n\to\infty}c_{n}n^{-1/2}\log^{5/2}n=0$, then for both the Poisson and
Binomial models, we have

$\lim_{n\to\infty}\mathbb{P}\Big{[}\mathcal{L}_{n}^{(1)}\left\\{\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\hat{\lambda})\right\\}\geq\mathcal{L}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
or})+c_{n}^{-1}\epsilon\Big{]}=0\text{ for any }\epsilon>0.$

Obviously, the estimator $\bm{\delta}_{(1)}^{\sf neb}(\lambda^{\sf orc}_{1})$
is lower bounded by the risk of the optimal solution $\bm{\delta}_{(1)}^{\pi}$
(Lemma 1). Next we study the _asymptotic optimality_ of $\bm{\delta}^{\sf
neb}_{(1)}$, which aims to provide decision theoretic guarantees on
$\bm{\delta}^{\sf neb}_{(1)}$ in relation to $\bm{\delta}_{(1)}^{\pi}$.
Theorem 3A establishes the optimality theory by showing that (a) the largest
coordinate-wise gap between ${\bm{\delta}}^{\sf neb}_{(1)}({\hat{\lambda}})$
and $\bm{\delta}_{(1)}^{\pi}$ is asymptotically small, and (b) the estimation
loss of the NEB estimator converges in probability to the loss of the
corresponding Bayes estimator as $n\to\infty$.

######  Theorem 3A.

Under the conditions of Theorem 2A, if
$\lim_{n\to\infty}c_{n}n^{-1/2}\log^{4}n=0$, then for both the Poisson and
Binomial models, we have

$c_{n}\Big{\|}{\bm{\delta}}_{(1)}^{\sf
neb}({\hat{\lambda}})-\bm{\delta}_{(1)}^{\pi}\Big{\|}_{\infty}=o_{p}(1).$

Furthermore, under the same conditions, we have for both the Poisson and
Binomial models,

$\lim_{n\to\infty}\mathbb{P}\Big{[}\mathcal{L}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\hat{\lambda}))\geq\mathcal{L}_{n}^{(1)}(\bm{\theta},\bm{\delta}_{(1)}^{\pi})+c_{n}^{-1}\epsilon\Big{]}=0\text{
for any }\epsilon>0.$

######  Remark 6.

The second part of the statement of Theorem 3A follows from the first part and
Lemma 4. In Appendix A.2, we discuss the counterpart to Theorem 3A under the
squared error loss $\mathcal{L}_{n}^{(0)}$.

##  4 Numerical Results

In this section we first discuss, in Section 4.1, the implementation details
of the convex program (8) and bandwidth selection process (9) (see also (19)
in Appendix A.2). Then we investigate the numerical performance of the NEB
estimator for Poisson and Binomial compound decision problems, respectively in
Sections 4.2 and 4.3. Both regular and scales losses will be considered. Our
numerical results demonstrate that the proposed the NEB estimator enjoys
superior numerical performance and the efficiency gain over competitive
methods is substantial in many settings.

We have developed an R package, npeb, to implement our proposed NEB estimator
in Definitions 1A (and Definition 1B in Appendix A.1), for the Poisson and
Binomial models under both regular and scaled losses. Moreover, the R code
that reproduces the numerical results in simulations can be downloaded from
the following link:

https://github.com/trambakbanerjee/DLE_paper.

###  4.1 Implementation Details

For a fixed $\lambda$ we use the R-package CVXR (Fu et al., 2017) to solve the
optimization problem in Equations (8) [and (15) in Appendix A.1]. As discussed
in Section 2.2, in the implementation under the scaled squared error loss
($k=1$) the linear inequality constraints, given by
$\mathcal{\bm{A}}\bm{h}\preceq\bm{b}$, ensure that the resulting decision rule
${\bm{\delta}}^{\sf neb}_{(1)}(\lambda)$ is monotonic, while the equality
constraints $\mathcal{\bm{C}}\bm{h}=\bm{d}$ handle boundary cases that involve
$y_{i}=0$ and ties. Moreover, since $w_{p}^{(1)}(y)>0$, the inequality
constraints also ensure that $h_{i}<1$ whenever $y_{i}>0$. Implementation
under the squared error loss ($k=0$) follows along similar lines and the
inequality constraints in this case ensure that $h_{i}+y_{i}>0$ whenever
$y_{i}\geq 0$.

A data-driven choice of the tuning parameter $\lambda$ is obtained by first
solving problems (8) and (15) over a grid of $\lambda$ values, i.e.
$\\{\lambda_{1},\ldots,\lambda_{s}\\}$, and then computing the corresponding
asymptotic risk estimate $\texttt{ARE}_{n}^{(k)}(\lambda_{j})$ for
$j=1,\ldots,s$. Then $\lambda$ is chosen according to

$\hat{\lambda}_{k}\coloneqq\operatorname*{arg\,min}_{\lambda\in{\lambda_{1},\ldots,\lambda_{s}}}\texttt{ARE}_{n}^{(k)}(\lambda),$

where $k\in\\{0,1\\}$. For all simulations and real data analyses considered
in this paper, we have fixed $s=10$ and employed an equi-spaced grid over
$[10,10^{2}]$.

###  4.2 Simulations: Poisson Distribution

In this section, we generate observations
$Y_{i}~{}|~{}\theta_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\texttt{Poi}(\theta_{i})$ for $i=1,\ldots,n$ and vary $n$ from
$500$ to $5000$ in increments of $500$. We consider four different scenarios
to simulate $\theta_{i}$. For each scenario, the following competing
estimators are considered:

  * •

the proposed estimator, denoted NEB;

  * •

the oracle NEB estimator $\bm{\delta}^{\sf or}\coloneqq{\bm{\delta}^{\sf
neb}}(\lambda^{\sf orc})$, denoted NEB OR;

  * •

the estimator of Poisson means based on Brown et al. (2013), denoted BGR;

  * •

the estimator of Poisson means based on Koenker and Gu (2017), denoted KM;

  * •

Tweedie’s formula based on Efron (2011) for the Poisson model, denoted TF OR;

  * •

Tweedie’s formula for the Normal means problem based on transformed data and
the convex optimization approach in Koenker and Mizera (2014), denoted TF
Gauss. The approach using transformation was suggested by Brown et al. (2013).

The risk performance of the TF OR method relies heavily on the choice of a
suitable bandwidth parameter $h>0$. We use the oracle loss estimate $h^{\sf
orc}$, which is obtained by minimizing the true loss $\mathcal{L}_{n}^{(0)}$.
The TF Gauss methodology is only applicable for the Normal means problem, and
uses a variance stabilization transformation on $Y_{i}$ to get
$Z_{i}=2\sqrt{Y_{i}+0.25}$. The $Z_{i}$ are then treated as approximate Normal
random variables with mean $\mu_{i}$ and variances $1$. To estimate normal
means $\mu_{i}$, we use the NPMLE approach proposed by Koenker and Mizera
(2014). Finally, $\theta_{i}$ is estimated as $0.25\hat{\mu_{i}}^{2}$.

It is important to note that the competitors to our NEB estimator only focus
on the regular loss $\mathcal{L}_{n}^{(0)}$. Nevertheless, in our simulation
we assess the performance of these estimators for estimating $\bm{\theta}$
under both $\mathcal{L}_{n}^{(0)}$ and $\mathcal{L}_{n}^{(1)}$. Consider the
following settings:

Scenario 1: We generate $\theta_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}\texttt{Unif}(0.5,15)$ for $i=1,\ldots,n$.

Scenario 2: We generate $\theta_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.75~{}\texttt{Gamma}(5,1)+0.25~{}\texttt{Gamma}(10,1)$ for
$i=1,\ldots,n$.

In the next two scenarios we assess the robustness of the five competing
estimators to departures from the Poisson model. Specifically consider the
Conway-Maxwell-Poisson distribution (Shmueli et al., 2005)
$\texttt{CMP}(\theta_{i},\nu)$. The CMP distribution is a generalization of
some well-known discrete distributions. With $\nu<1$, CMP represents a
discrete distribution that has longer tails than the Poisson distribution with
parameter $\theta_{i}$.

Scenario 3: We generate $\theta_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.5~{}\delta_{\\{10\\}}+0.5~{}\texttt{Gamma}(5,2)$ for each $i$
and let

$Y_{i}~{}|~{}\theta_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}0.8~{}\texttt{Poi}(\theta_{i})+0.2~{}\texttt{CMP}(\theta_{i},\nu),$

where we let $\nu=0.8$ for the CMP distribution.

Scenario 4:  We let $\bm{\theta}$ to be an equi-spaced vector of length $n$ in
$[1,5]$ and let $Y_{i}~{}|~{}\theta_{i}$ be distributed as the CMP
distribution with parameters $\theta_{i}$ and $\nu=0.8$.

(a) Scenario 1: Estimation of $\bm{\theta}$ under loss $\mathcal{L}_{n}^{(1)}$
where $\theta_{i}\stackrel{{\scriptstyle iid.}}{{\sim}}\text{Unif}(0.5,15)$.

(b) Scenario 2: Estimation of $\bm{\theta}$ under loss $\mathcal{L}_{n}^{(1)}$
where $\theta_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.75~{}\texttt{Gamma}(5,1)+0.25~{}\texttt{Gamma}(10,1)$.

(c) Scenario 3: Estimation of $\bm{\theta}$ under loss $\mathcal{L}_{n}^{(1)}$
where $\theta_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.5~{}\delta_{\\{10\\}}+0.5~{}\texttt{Gamma}(5,2)$.

(d) Scenario 4: Estimation of $\bm{\theta}$ under loss $\mathcal{L}_{n}^{(1)}$
where $\bm{\theta}$ is an equi-spaced vector of length $n$ in $[1,5]$ and
$Y_{i}|\theta_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\texttt{CMP}(\theta_{i},0.8)$ .

Figure 1: Poisson compound decision problem under scaled squared error loss:
Risk estimates of the various estimators for scenarios 1 to 4.

(a) Scenario 1: Estimation of $\bm{\theta}$ under loss $\mathcal{L}_{n}^{(0)}$
where $\theta_{i}\stackrel{{\scriptstyle iid.}}{{\sim}}\texttt{Unif}(0.5,15)$.

(b) Scenario 2: Estimation of $\bm{\theta}$ under loss $\mathcal{L}_{n}^{(0)}$
where $\theta_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.75~{}\texttt{Gamma}(5,1)+0.25~{}\texttt{Gamma}(10,1)$.

(c) Scenario 3: Estimation of $\bm{\theta}$ under loss $\mathcal{L}_{n}^{(0)}$
where $\theta_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.5~{}\delta_{\\{10\\}}+0.5~{}\texttt{Gamma}(5,2)$.

(d) Scenario 4: Estimation of $\bm{\theta}$ under loss $\mathcal{L}_{n}^{(0)}$
where $\bm{\theta}$ is an equi-spaced vector of length $n$ in $[1,5]$ and
$Y_{i}|\theta_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\texttt{CMP}(\theta_{i},0.8)$ .

Figure 2: Poisson compound decision problem under squared error loss: Risk
estimates of the various estimators for scenarios 1 to 4.

Table 1: Poisson compound decision problem under scaled squared error loss:
Risk ratios
$\mathcal{R}_{n}^{(1)}(\bm{\theta},\cdot)/\mathcal{R}_{n}^{(1)}(\bm{\theta},\bm{\delta}^{\sf
neb}_{k})$ at $n=5000$ for estimating $\bm{\theta}$.

| Scenario
---|---
Method | 1 | 2 | 3 | 4
KM | 1.10 | 1.04 | 1.17 | 1.37
TF Gauss | 1.03 | 1.01 | 1.10 | 1.25
TF OR | 1.00 | 1.02 | 1.14 | 1.22
BGR | 1.22 | 1.07 | 1.16 | 1.37
NEB | 1.00 | 1.00 | 1.00 | 1.00
NEB OR | 1.00 | 1.00 | 0.90 | 1.00

Table 2: Poisson compound decision problem under squared error loss: Risk
ratios
$\mathcal{R}_{n}^{(0)}(\bm{\theta},\cdot)/\mathcal{R}_{n}^{(0)}(\bm{\theta},\bm{\delta}^{\sf
neb}_{k})$ at $n=5000$ for estimating $\bm{\theta}$.

| Scenario
---|---
Method | 1 | 2 | 3 | 4
KM | 1.00 | 1.00 | 1.59 | 1.21
TF Gauss | 1.00 | 1.01 | 1.51 | 1.08
TF OR | 1.07 | 1.12 | 1.66 | 1.12
BGR | 1.01 | 1.01 | 1.55 | 1.15
NEB | 1.00 | 1.00 | 1.00 | 1.00
NEB OR | 1.00 | 1.00 | 0.90 | 1.00

The performances of these four estimators are presented in figures 1 and 2
wherein the risk $\mathcal{R}_{n}^{(k)}(\bm{\theta},\cdot)$ of the various
estimators is estimated using $50$ Monte Carlo repetitions for varying $n$.
Tables 2 and 2 report the risk ratios
$\mathcal{R}_{n}^{(k)}(\bm{\theta},\cdot)/\mathcal{R}_{n}^{(k)}(\bm{\theta},\bm{\delta}^{\sf
neb}_{k})$ at $n=5000$ and for $k=1,0$ respectively, where a risk ratio bigger
than $1$ demonstrates a smaller estimation risk for the NEB estimator. For BGR
the modified cross validation approach of choosing the bandwidth parameter was
extremely slow in our simulations and we therefore report its risk performance
only at $n=5000$. From figures 1, 2 and tables 2, 2, we note that the NEB
estimator demonstrates an overall competitive risk performance. In particular,
we see that when estimation is conducted under loss $\mathcal{L}_{n}^{(1)}$
the risk ratios of the competing estimators in table 2 reflect a relatively
better performance of the NEB estimator which is not surprising considering
the fact that KM, TF Gauss and TF OR are designed to estimate $\bm{\theta}$
under loss $\mathcal{L}_{n}^{(0)}$. We note that TF Gauss is highly
competitive against KM (Koenker and Mizera, 2014) and this observation was
also reported in Brown et al. (2013). Of particular interest are Scenarios 3
and 4, which reflect the relative performance of these estimators under
departures from the Poisson model. The NEB estimator has a significantly
better risk performance in these settings across both types of losses.

###  4.3 Simulations: Binomial Distribution

In this section, we generate $Y_{i}~{}|~{}q_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\mathcal{B}in(m_{i},~{}q_{i})$ for $i=1,\ldots,n$ and vary $n$
from $500$ to $5000$ in increments of $500$. We consider four different
scenarios to simulate $\theta_{i}=q_{i}/(1-q_{i})$ and for each scenario we
consider the following competing estimators:

  * •

the proposed estimator, denoted NEB;

  * •

the oracle NEB estimator $\bm{\delta}^{\sf or}\coloneqq{\bm{\delta}^{\sf
neb}}(\lambda^{\sf orc})$, denoted NEB OR;

  * •

the estimator of Binomial means based on Koenker and Gu (2017), denoted KM;

  * •

Tweedie’s formula for Binomial log odds based on Efron (2011) and Fu et al.
(2018), denoted TF OR;

  * •

Tweedie’s formula for the Normal means problem based on transformed data and
the convex optimization approach in Koenker and Mizera (2014), denoted TF
Gauss.

For TF OR, analogous to the Poisson case, we continue to use the oracle loss
estimate $h^{\sf orc}$ as a choice for the bandwidth parameter. Since the TF
Gauss methodology is only applicable for the Normal means problem, it uses a
variance stabilization transformation on $Y_{i}$ to get $Z_{i}={\sf
arcsin}\sqrt{(Y_{i}+0.25)/(m_{i}+0.5)}$. The $Z_{i}$ are then treated as
approximate Normal random variables with mean $\mu_{i}$, variances
$(4m_{i})^{-1}$, and estimate of the means $\mu_{i}$’s are obtained using the
NPMLE approach of Koenker and Mizera (2014). Finally, $q_{i}$ is estimated as
$\\{\sin(\hat{\mu}_{i})\\}^{2}$. We note that unlike the Poisson case
discussed earlier, the competitors to our NEB estimator do not directly
estimate the odds $\bm{\theta}$. For instance, under a squared error loss both
KM and TF Gauss estimate the success probabilities $\bm{q}$ while TF OR
estimates $\log\bm{\theta}$. Nevertheless, in this simulation experiment we
assess the performance of these estimators for estimating the odds under both
squared error loss and its scaled version.

The following settings are considered in our simulation:

Scenario 1: We generate $q_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.4~{}\delta_{\\{0.5\\}}+0.6~{}\texttt{Beta}(2,5)$ and fix
$m_{i}=5$ for $i=1,\ldots,n$.

Scenario 2: We let $\theta_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.8~{}\delta_{\\{0.5\\}}+0.2~{}\texttt{Gamma}(1,2)$ and fix
$m_{i}=10$ for $i=1,\ldots,n$. In this scenario we let the odds $\theta_{i}$
arise from a mixture model that has $80\%$ point mass at $0.5$.

Scenario 3: $\theta_{i}\stackrel{{\scriptstyle iid.}}{{\sim}}\chi^{2}_{2}$ and
fix $m_{i}=5$ for $i=1,\ldots,n$. This scenario is similar to scenario 2 where
we let the odds $\theta_{i}$ arise from a Chi-square distribution with 2
degrees of freedom.

Scenario 4: We generate $q_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.5~{}\texttt{Beta}(1,1)+0.5~{}\texttt{Beta}(1,3)$ and fix
$m_{i}=10$ for $i=1,\ldots,n$.

The simulation results are presented in Figures 3 and 4 wherein the risks of
various estimators are calculated by averaging over $50$ Monte Carlo
repetitions for varying $n$. Tables 4 and 4 report the risk ratios
$\mathcal{R}_{n}^{(k)}(\bm{\theta},\cdot)/\mathcal{R}_{n}^{(k)}(\bm{\theta},\bm{\delta}^{\sf
neb}_{k})$ at $n=5000$ and for $k=1,0$ respectively, where a risk ratio bigger
than $1$ demonstrates a smaller estimation risk for the NEB estimator.

(a) Scenario 1: Estimation of odds $\bm{\theta}$ under loss
$\mathcal{L}_{n}^{(1)}$ where $q_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.4~{}\delta_{\\{0.5\\}}+0.6~{}\texttt{Beta}(2,5)$ and
$m_{i}=5$.

(b) Scenario 2: Estimation of odds $\bm{\theta}$ under loss
$\mathcal{L}_{n}^{(1)}$ where $\theta_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.8~{}\delta_{\\{0.5\\}}+0.2~{}\texttt{Gamma}(1,2)$ and fix
$m_{i}=10$.

(c) Scenario 3: Estimation of odds $\bm{\theta}$ under loss
$\mathcal{L}_{n}^{(1)}$ where $\theta_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}\chi^{2}_{2}$ and $m_{i}=5$.

(d) Scenario 4: Estimation of odds $\bm{\theta}$ under loss
$\mathcal{L}_{n}^{(1)}$ where $q_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.5~{}\mathcal{B}eta(1,1)+0.5~{}\texttt{Beta}(1,3)$ and
$m_{i}=10$.

Figure 3: Binomial compound decision problem under scaled squared error loss:
Risk estimates of the various estimators for Scenarios 1 to 4.

(a) Scenario 1: Estimation of odds $\bm{\theta}$ under loss
$\mathcal{L}_{n}^{(0)}$ where $q_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.4~{}\delta_{\\{0.5\\}}+0.6~{}\texttt{Beta}(2,5)$ and
$m_{i}=5$.

(b) Scenario 2: Estimation of odds $\bm{\theta}$ under loss
$\mathcal{L}_{n}^{(0)}$ where $\theta_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.8~{}\delta_{\\{0.5\\}}+0.2~{}\texttt{Gamma}(1,2)$ and
$m_{i}=10$.

(c) Scenario 3: Estimation of odds $\bm{\theta}$ under loss
$\mathcal{L}_{n}^{(0)}$ where $\theta_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}\chi^{2}_{2}$ and $m_{i}=5$.

(d) Scenario 4: Estimation of odds $\bm{\theta}$ under loss
$\mathcal{L}_{n}^{(0)}$ where $q_{i}\stackrel{{\scriptstyle
iid.}}{{\sim}}0.5~{}\mathcal{B}eta(1,1)+0.5~{}\texttt{Beta}(1,3)$ and
$m_{i}=10$.

Figure 4: Binomial compound decision problem under squared error loss: Risk
estimates of the various estimators for Scenarios 1 to 4.

Table 3: The Binomial compound decision problem under scaled squared error
loss: Risk ratios
$\mathcal{R}_{n}^{(1)}(\bm{\theta},\cdot)/\mathcal{R}_{n}^{(1)}(\bm{\theta},\bm{\delta}^{\sf
neb}_{k})$ at $n=5000$ for estimating $\bm{\theta}$.

| Scenario
---|---
Method | 1 | 2 | 3 | 4
KM | 1.22 | 1.38 | 1.26 | 1.03
TF Gauss | 1.23 | 1.51 | 1.34 | 1.08
TF OR | $>10$ | $>10$ | $>10$ | $>10$
NEB | 1.00 | 1.00 | 1.00 | 1.00
NEB OR | 1.00 | 1.00 | 1.00 | 1.00

Table 4: The Binomial compound decision problem under the usual squared error
loss: Risk ratios
$\mathcal{R}_{n}^{(0)}(\bm{\theta},\cdot)/\mathcal{R}_{n}^{(0)}(\bm{\theta},\bm{\delta}^{\sf
neb}_{k})$ at $n=5000$ for estimating $\bm{\theta}$

| Scenario
---|---
Method | 1 | 2 | 3 | 4
KM | 1.01 | 1.08 | 1.08 | 1.06
TF Gauss | 1.06 | 1.06 | 1.09 | 1.17
TF OR | $>10$ | $>10$ | $>10$ | $>10$
NEB | 1.00 | 1.00 | 1.00 | 1.00
NEB OR | 1.00 | 0.99 | 0.99 | 1.00

We can see from the simulation results that the NEB estimator demonstrates an
overall superior risk performance than its competitors. In particular, we see
that when estimation is conducted under loss $\mathcal{L}_{n}^{(1)}$ the risk
ratios of the competing estimators in Table 4 reflect a significantly better
performance of the NEB estimator. This is not surprising because KM, TF Gauss
and TF OR are designed to estimate $\bm{q}$ under loss
$\mathcal{L}_{n}^{(0)}$. This also explains the relatively improved
performance of these estimators as seen through their risk ratios in table 4
wherein the estimation is conducted under the usual squared error loss
$\mathcal{L}_{n}^{(0)}$. Across the four scenarios, TF OR exhibits the poorest
performance and appears to suffer from the fragmented approach of estimating
the gradient of the log density $\log p(\bm{y})$ wherein $p(\bm{y})$ and its
first derivative with respect to $\bm{y}$ are estimated separately using a
Gaussian kernel with common bandwidth $h^{\sf orc}$. The approach of using a
variance stabilizing transformation to convert the data to approximate
normality renders TF Gauss highly competitive to KM (Koenker and Mizera,
2014). A similar phenomenon was also reported in Brown et al. (2013) in the
context of the Poisson model. However, under the Binomial model, when the
primary goal is to estimate the odds $\bm{\theta}$, the risk ratios reported
in Tables 4 and 4 suggest that the proposed NEB estimator is by far the best
amongst these competitors under both types of losses.

##  5 Real Data Analyses

This section illustrates the proposed method for estimating the Juvenile
Delinquency rates from Poisson models and news popularity in Binomial models.

###  5.1 Estimation of Juvenile Delinquency rates

In this section we consider an application for analysis of the Uniform Crime
Reporting Program (UCRP) Database (US Department of Justice and Federal Bureau
of Investigation, 2014) that holds county-level counts of arrests and offenses
ranging from robbery to weapons violations in 2012. The database is maintained
by the National Archive of Criminal Justice Data (NACJD) and is one of the
most widely used database for research related to factors that affect juvenile
delinquency (JD) rates across the United States; see for example Aizer and
Doyle Jr (2015) and Damm and Dustmann (2014); Koski et al. (2018). A
preliminary and important goal in these analyses is to estimate the JD rates
based on the observed arrest data and determine the counties that are amongst
the worst or least affected. However with almost 3,000 counties being
evaluated the JD rates are susceptible to selection bias, wherein some of the
data points are in the extremes merely by chance and traditional estimators
may underestimate or overestimate the corresponding means, especially in
counties with fewer total number of arrests across all age groups.

Figure 5: Observed Juvenile Delinquency rates in 2012. The top $500$ and
bottom $500$ counties are plotted. The data on Florida arrests is not
available in the US Department of Justice and Federal Bureau of Investigation
(2014) database.

Figure 6: Estimated Juvenile Delinquency rates of the $1000$ counties exhibited in figure 6. Left: Estimation under squared error loss $(k=0)$. Right: Estimation under scaled squared error loss $(k=1)$. Table 5: Loss ratios of the competing methods for estimating $\bm{\theta}$. $(n=2,803)$ | Loss ratios
---|---
Method | $k=1$ | $k=0$
NEB | 1.00 | 1.00
BGR | 1.18 | 1.03
KM | 1.19 | 1.03
TF Gauss | 1.12 | 1.01
TF OR | 1.11 | 1.05

For the purpose of our analyses, we use the 2012 UCRP data that spans
$n=3,178$ counties in the U.S. and consider estimating the mean JD rate
$\bm{\theta}$ as a vector of Poisson means. The observed data for county $i$
is denoted $y_{i}$, which represents the number of juvenile arrests expressed
as a percentage of total arrests in that county in Year 2012. We assume that
$Y_{i}~{}|~{}\theta_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\texttt{Poi}(\theta_{i})$ for $i=1,\ldots,n$. Figure 6 plots the
observed data for the top $500$ and the bottom $500$ counties that have at
least $1$ juvenile arrest. Campbell county in South Dakota, followed by Fulton
county in New York, exhibits the highest observed JD rates in Year 2012.

As discussed in Section 4.2, we consider the following 5 estimators of
$\bm{\theta}$: NEB, BGR (Brown et al., 2013), KM (Koenker and Gu, 2017), TF OR
(Efron, 2011) and TF Gauss (Koenker and Mizera, 2014; Brown et al., 2013). We
use the 2014 UCRP data (US Department of Justice and Federal Bureau of
Investigation, 2017) to compare their estimation accuracies under both
$\mathcal{L}_{n}^{(0)}$ and $\mathcal{L}_{n}^{(1)}$ losses. The data were
cleaned prior to any analyses which ensured that all counties in the year 2012
had at least one arrest (juvenile or not). This resulted in $n=2803$ counties
where all methods are applied to.

In Figure 6 we visualize the shrinkage estimates of JD rates for those $1000$
counties considered in Figure 6. The left plot presents the estimates under
the squared error loss while right plot presents the results under the scaled
squared error loss. Notably, the scaled error loss exhibits a larger magnitude
of shrinkage for the bigger observations than the squared error loss. Table 5
reports the loss ratios
$\mathcal{L}_{n}^{(k)}(\bm{\theta},\bm{\delta})/\mathcal{L}_{n}^{(k)}(\bm{\theta},{\bm{\delta}}^{\sf
neb}_{(k)})$ where for any estimator $\bm{\delta}$ of $\bm{\theta}$, a ratio
bigger than $1$ indicates a smaller estimation loss for ${\bm{\delta}}^{\sf
neb}_{(k)}$. We can see that for estimating $\bm{\theta}$, all four
competitors exhibit loss ratios bigger than $1$ under the scaled squared error
loss ($k=1$). This is not surprising since these competitors are designed to
estimate $\bm{\theta}$ under the regular squared error loss ($k=0$).
Interestingly, even under the regular loss, the NEB estimator continues to
provide a better estimation accuracy than TF OR, BGR and KM, and demonstrates
a competitive performance against TF Gauss.

###  5.2 News popularity in social media platforms

Journalists and editors often face the critical task of assessing the
popularity of various news items and determining which articles are likely to
become popular; hence existing content generation resources can be efficiently
managed and optimally allocated to avenues with maximum potential. Due to the
dynamic nature of the news articles, popularity is usually measured by how
quickly the article propagates (frequency) and the number of readers that the
article can reach (severity) through social media platforms like Twitter,
Youtube, Facebook and LinkedIn. As such predicting these two aspects of
popularity based on early trends is extremely valuable to journalists and
content generators (Bandari et al., 2012).

Table 6: Loss ratios of the competing methods for estimating $\bm{\theta}$.
News article genre: Economy and social media: Facebook

$(n=3,972)$ | Loss Ratios
---|---
Method | $k=1$ | $k=0$
NEB | 1.00 | 1.00
KM | 12.25 | 13.25
TF Gauss | 4.06 | 3.25
TF OR | 81.57 | 36.24

Table 7: Loss ratios of the competing methods for estimating $\bm{\theta}$.
News article genre: Microsoft and social media: LinkedIn

$(n=3,850)$ | Loss Ratios
---|---
Method | $k=1$ | $k=0$
NEB | 1.00 | 1.00
KM | 41.06 | 59.64
TF Gauss | 9.31 | 7.44
TF OR | 36.14 | 12.76

In this section, we assess the popularity of several news items based on their
frequency of propagation and analyze a dataset from Moniz and Torgo (2018)
that holds $48$ hours worth of social media feedback data on a large
collection of news articles since the time of first publication. For the
purposes of our analysis, we consider two popular genres of news from this
data set: Economy and Microsoft, and examine how frequently these articles
were shared in Facebook and LinkedIn, respectively, over a period of $48$
hours from the time of their first publication. Each news article in the data
has a unique identifier and $16$ consecutive time intervals, each of length
$180$ minutes, to detect whether the article was shared at least once in that
time interval. Let $Z_{ij}=1$ if article $i$ was shared in time interval $j$
and $0$ otherwise, where $i=1,\ldots,n$ and $j=1,\ldots,16$. Suppose
$q_{ij}\in[0,1]$ denote the probability that news article $i$ is shared in
interval $j$. We let $q_{ij}=q_{i}$ for all $j=1,\ldots,16$ and assume that
for each $i$, $Z_{ij}$ are independent realizations from
$\texttt{Ber}(q_{i})$. It follows that
$Y_{i}=\sum_{j=1}^{8}Z_{ij}\stackrel{{\scriptstyle
ind.}}{{\sim}}\texttt{Bin}(8,q_{i})$.

To assess the popularity of article $i$, we estimate its odds of sharing given
by $\theta_{i}=q_{i}/(1-q_{i})$ and consider the following 4 estimators of
$\bm{\theta}$: NEB, KM (Koenker and Gu, 2017), TF Gauss (Koenker and Mizera,
2014) and TF OR (Efron, 2011; Fu et al., 2018). We use the data on time points
$j=9,\ldots,16$ to compare the estimation accuracy of these estimators under
both $\mathcal{L}_{n}^{(0)}$ and $\mathcal{L}_{n}^{(1)}$ losses. Tables 7 and
7 report the loss ratios
$\mathcal{L}_{n}^{(k)}(\bm{\theta},\bm{\delta})/\mathcal{L}_{n}^{(k)}(\bm{\theta},{\bm{\delta}}^{\sf
neb}_{(k)})$ where for any estimator $\bm{\delta}$ of $\bm{\theta}$, a ratio
bigger than $1$ indicates a smaller estimation loss for ${\bm{\delta}}^{\sf
neb}_{(k)}$. We observe that the three competitors to the NEB estimator
exhibit loss ratios substantially bigger than $1$ under both the losses. This
is not surprising since these competitors are designed to estimate $\bm{q}$
and $\log\bm{\theta}$ under a squared error loss ($k=0$). However, when the
primary goal is to estimate the odds, the proposed NEB estimator is by far the
best amongst these competitors under both losses.

## Acknowledgments

G. Mukherjee’s work was partly supported by the NSF grant DMS-1811866. W.
Sun’s work was partly supported by the NSF grant DMS-1712983. Q. Liu’s work is
supported in part by NSF CRII 1830161 and NSF CAREER 1846421.

Supplementary Material for “EB Estimation in Discrete Linear Exponential
Family”

In this supplement, we first present in Appendix A the results for the NEB
estimator under the regular squared loss, then provide in Appendix B the
proofs and technical details of all theories in the main text and Appendix A.

##  Appendix A Results Under the Squared Error Loss

###  A.1 The NEB estimator

In this section we discuss the estimation of $\bm{w}_{p}^{(0)}$ that appear in
lemma 1 under the usual squared error loss $(k=0)$. Let $Y$ be a non-negative
integer-valued random variable with probability mass function (pmf) $p$ and
define

$h_{0}^{(0)}(y)=\dfrac{y+1}{w_{p}^{(0)}(y)}-y~{},y\in\\{0\\}\cup\mathbb{N}$
(14)

Suppose
$\mathcal{K}_{\lambda}(y,y^{\prime})=\exp\\{-0.5\lambda^{-2}(y-y^{\prime})^{2}\\}$
be the positive definite RBF kernel with bandwidth parameter
$\lambda\in\Lambda$ where $\Lambda$ is a compact subset of $\mathbb{R}^{+}$
bounded away from $0$. Given observations $\bm{y}=(y_{1},\ldots,y_{n})$ from
model (2), let
$\bm{h}_{0}^{(0)}=(h_{0}^{(0)}(y_{1}),\ldots,h_{0}^{(0)}(y_{n}))$ and define
the following $n\times n$ matrices:
$n^{2}\bm{K}_{\lambda}=[\mathcal{K}_{\lambda}(y_{i},y_{j})]_{ij}$,
$n^{2}{\Delta\bm{K}}_{\lambda}=[\Delta_{y_{i}}\mathcal{K}_{\lambda}(y_{i},y_{j}+1)]_{ij}$
and
$n^{2}{\Delta_{2}\bm{K}}_{\lambda}=[\Delta_{y_{i},y_{j}}\mathcal{K}_{\lambda}(y_{i},y_{j})]_{ij}$
where
$\Delta_{y}\mathcal{K}_{\lambda}(y,y^{\prime})=\mathcal{K}_{\lambda}(y+1,y^{\prime})-\mathcal{K}_{\lambda}(y,y^{\prime})$
and
$\Delta_{y,y^{\prime}}\mathcal{K}_{\lambda}(y,y^{\prime})=\Delta_{y^{\prime}}\Delta_{y}\mathcal{K}_{\lambda}(y,y^{\prime})=\Delta_{y}\Delta_{y^{\prime}}\mathcal{K}_{\lambda}(y,y^{\prime})$.

######  Definition 1B (NEB estimator of $\theta_{i}$).

Consider the DLE Model (2) with loss $\ell^{(0)}(\theta_{i},\delta_{i})$. For
a fixed $\lambda\in\Lambda$, let
$\hat{w}_{i}^{(0)}(\lambda)=(y_{i}+1)/(y_{i}+\hat{h}_{i}^{(0)}(\lambda))$ and
$\hat{\bm{h}}_{n}^{(0)}(\lambda)=\left\\{\hat{h}_{1}^{(0)}(\lambda),\ldots,\hat{h}_{n}^{(0)}(\lambda)\right\\}$
be the solution to the following quadratic optimization problem:

$\displaystyle\min_{\bm{h}\in\bm{H}_{n}}~{}\hat{\mathbb{M}}_{\lambda,n}(\bm{h})=\bm{h}^{T}\bm{K}_{\lambda}\bm{h}+2\bm{h}^{T}\Delta\bm{K}_{\lambda}\bm{y}+\bm{y}^{T}\Delta_{2}\bm{K}_{\lambda}\bm{y},$
(15)

where
$\bm{H}_{n}=\\{\bm{h}=(h_{1},\ldots,h_{n}):\mathcal{A}\bm{h}\preceq\bm{b},~{}\mathcal{C}\bm{h}=\bm{d}\\}$
is a convex set and $\mathcal{A},\mathcal{C},\bm{b},\bm{d}$ are known real
matrices and vectors that enforce linear constraints on the components of
$\bm{h}$. Then the NEB estimator for a fixed $\lambda$ is given by
${\bm{\delta}}^{\sf neb}_{(0)}(\lambda)=\left\\{{\delta}_{(0),i}^{\sf
neb}(\lambda):1\leq i\leq n\right\\}$, where

${\delta}_{(0),i}^{\sf
neb}(\lambda)=\dfrac{a_{y_{i}}/a_{y_{i}+1}}{\hat{w}_{i}^{(0)}(\lambda)},\quad\mbox{
if }y_{i}\in\\{0,1,2,\ldots\\}$

######  Theorem 2B.

Let $\mathcal{K}_{\lambda}(\cdot,\cdot)$ be the positive definite RBF kernel
with bandwidth parameter $\lambda\in\Lambda$. If
$\lim_{n\to\infty}c_{n}n^{-1/2}\log n=0$ then, under assumptions $A1-A3$, we
have for any $\lambda\in\Lambda$,

$\lim_{n\to\infty}\mathbb{P}\Big{[}\Big{\|}\hat{\bm{w}}_{n}^{(0)}(\lambda)-{\bm{w}}_{p}^{(0)}\Big{\|}_{2}\geq
c_{n}^{-1}\epsilon\Big{]}=0,~{}\text{for any }\epsilon>0$

where
$\hat{\bm{w}}_{n}^{(0)}(\lambda)=[(Y_{i}+1)/(\hat{h}_{i}^{(0)}(\lambda)+Y_{i})]_{i}$.

We now provide some motivation behind the minimization problem in definition
1B for estimating the ratio functionals ${\bm{w}}_{p}^{(0)}$. Suppose
$\tilde{p}$ be a probability mass function on the support of $Y$ and define

$\mathcal{S}_{\lambda}[\tilde{p}](p)=\mathbb{E}_{p}\Big{[}(\tilde{h}^{(0)}(Y)-h_{0}^{(0)}(Y))\mathcal{K}_{\lambda}(Y+1,Y^{\prime}+1)(\tilde{h}^{(0)}(Y^{\prime})-h_{0}^{(0)}(Y^{\prime}))\Big{]}$
(16)

where $h_{0}^{(0)},\tilde{h}^{(0)}$ are as defined in equation (14) and
$Y,Y^{\prime}$ are i.i.d copies from the marginal distribution that has mass
function $p$. $\mathcal{S}_{\lambda}[\tilde{p}](p)$ in equation (16) is the
Kernelized Stein’s Discrepancy (KSD) measure that can be used to distinguish
between two distributions with mass functions $p,\tilde{p}$ such that
$\mathcal{S}_{\lambda}[\tilde{p}](p)\geq 0$ and
$\mathcal{S}_{\lambda}[\tilde{p}](p)=0$ if and only if $p=\tilde{p}$ (Liu et
al., 2016; Chwialkowski et al., 2016). Moreover
$\mathcal{S}_{\lambda}[\tilde{p}](p)=\mathbb{E}_{(U,V)\stackrel{{\scriptstyle
i.i.d}}{{\sim}}}p\Big{[}\kappa_{\lambda}[\tilde{p}](U,V)\Big{]}$ where
$\kappa_{\lambda}[\tilde{p}](u,v)$ is

$\tilde{h}^{(0)}(u)\tilde{h}^{(0)}(v)\mathcal{K}_{\lambda}(u,v)+\tilde{h}^{(0)}(u)v\Delta_{v}\mathcal{K}_{\lambda}(u+1,v)+\tilde{h}^{(0)}(v)u\Delta_{u}\mathcal{K}_{\lambda}(u,v+1)+uv\Delta_{u,v}\mathcal{K}_{\lambda}(u,v)$
(17)

An empirical evaluation scheme for $\mathcal{S}_{\lambda}[\tilde{p}](p)$ is
given by $\mathcal{S}_{\lambda}[\tilde{p}](\hat{p}_{n})$ where

$\mathcal{S}_{\lambda}[\tilde{p}](\hat{p}_{n})=\dfrac{1}{n^{2}}\sum_{i,j=1}^{n}\kappa_{\lambda}[\tilde{p}](y_{i},y_{j})$
(18)

and $\bm{Y}=(Y_{1},\ldots,Y_{n})$ is a random sample from the marginal
distribution with mass function $p$ with empirical CDF $\hat{p}_{n}$. Note
that $\kappa_{\lambda}[\tilde{p}](u,v)$ in equation (17) involves $\tilde{p}$
only through $\tilde{h}^{(0)}$ and may analogously be denoted by
$\kappa_{\lambda}[\tilde{h}(u),\tilde{h}(v)](u,v)\coloneqq\kappa_{\lambda}[\tilde{p}](u,v)$
where we have dropped the superscript from $\tilde{h}$ that indicates that the
loss in question is the regular squared error loss. This slight abuse of
notation is harmless as the discussion in this section is geared towards the
squared error loss only.

Under the empirical Bayes compound estimation framework of model (2), our goal
is to estimate $h_{0}^{(0)}$. To do that we minimize
$\mathcal{S}_{\lambda}[\tilde{p}](\hat{p}_{n})$ in equation (18) with respect
to the unknowns $\bm{\tilde{}}{h}=(\tilde{h}(y_{1}),\ldots,\tilde{h}(y_{n}))$
and the sample criteria

$\hat{\mathbb{M}}_{\lambda,n}(\tilde{\bm{h}})=\dfrac{1}{n^{2}}\sum_{i,j=1}^{n}\kappa_{\lambda}[\bm{\tilde{}}{h}(y_{i}),\tilde{h}(y_{j})](y_{i},y_{j})$

is thus the objective function of the optimisation problem in equation (15)
with optimisation variables $h_{i}\equiv\tilde{h}(y_{i})$ for $i=1,\ldots,n$.

Note that $\hat{\mathbb{M}}_{\lambda,n}(\tilde{\bm{h}})$ above is a
V-statistic and a biased estimator of the population criteria
$\mathbb{M}_{\lambda}(\tilde{\bm{h}})$ that is defined in equation (3.1) with
$\kappa_{\lambda}[\tilde{h}(u),\tilde{h}(v)](u,v)$ given by equation (17).

###  A.2 Bandwidth choice and asymptotic properties

We propose the following asymptotic risk estimate of the true risk of
${\bm{\delta}}_{(0)}^{\sf neb}(\lambda)$ in the Poisson and Binomial model.

######  Definition 2B (ARE of ${\bm{\delta}}_{(0)}^{\sf neb}(\lambda)$ in the
Poisson model).

Suppose $Y_{i}~{}|~{}\theta_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\texttt{Pois}(\theta_{i})$. Under the loss
$\ell^{(0)}(\theta_{i},\cdot)$ an ARE of the true risk of
${\bm{\delta}}_{(0)}^{\sf neb}(\lambda)$ is

${\sf
ARE}^{(0,\mathcal{P})}_{n}(\lambda)=\dfrac{1}{n}\Bigl{\\{}\sum_{i=1}^{n}y_{i}(y_{i}-1)-2\sum_{i=1}^{n}y_{i}{\psi}_{\lambda}(y_{i})+\sum_{i=1}^{n}\\{{\delta}_{(0),i}^{\sf
neb}(\lambda)\\}^{2}\Bigr{\\}}$

where

$\psi_{\lambda}(y_{i})={\delta}_{(0),j}^{\sf neb}(\lambda),~{}y_{i}\geq 1$

with $j\in\\{1,\ldots,n\\}$ such that $y_{j}=y_{i}-1$.

######  Definition 3B (ARE of ${\bm{\delta}}_{(0)}^{\sf neb}(\lambda)$ in the
Binomial model).

Suppose $Y_{i}~{}|~{}q_{i}\sim\texttt{Bin}(m,q_{i})$ so that in equation (1)
$a_{y_{i}}={m\choose y_{i}}$ and $\theta_{i}=q_{i}/(1-q_{i})$ in equation (2).
Under the loss $\ell^{(0)}(\theta_{i},\cdot)$ an ARE of the true risk of
${\bm{\delta}}_{(0)}^{\sf neb}(\lambda)$ is

${\sf
ARE}^{(0,\mathcal{B})}_{n}(\lambda)=\dfrac{1}{n}\Bigl{\\{}\sum_{i=1}^{n}\dfrac{y_{i}(y_{i}-1)}{(m-y_{i}+2)(m-y_{i}+1)}-2\sum_{i=1}^{n}y_{i}{\psi}_{\lambda}(y_{i})+\sum_{i=1}^{n}\\{{\delta}_{(0),i}^{\sf
neb}(\lambda)\\}^{2}\Bigr{\\}}$

where

$\psi_{\lambda}(y_{i})={\delta}_{(0),j}^{\sf
neb}(\lambda)/(m-y_{i}+1),~{}y_{i}\geq 1$

with $j\in\\{1,\ldots,n\\}$ such that $y_{j}=y_{i}-1$.

Note that if for some index $i$, $y_{i}-1$ is not available in the observed
sample $\bm{y}$, ${\psi}_{\lambda}(y_{i})$ can be calculated using cubic
splines. We propose the following estimate of the tuning parameter $\lambda$
based on the ARE:

$\hat{\lambda}=\begin{cases}\operatorname*{arg\,min}_{\lambda\in{\Lambda}}{\sf
ARE}^{(0,\mathcal{P})}_{n}(\lambda),~{}\text{if~{}}Y_{i}~{}|~{}\theta_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\texttt{Pois}(\theta_{i})\\\
\operatorname*{arg\,min}_{\lambda\in{\Lambda}}{\sf
ARE}^{(0,\mathcal{B})}_{n}(\lambda),~{}\text{if~{}}Y_{i}~{}|~{}q_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\texttt{Bin}(m,q_{i})\end{cases}$ (19)

where a choice of $\Lambda=[10,10^{2}]$ work well in the simulations and real
data analyses of sections 4 and 5. Lemmata 3 and 2 continue to provide the
large-sample properties of the proposed ${\sf ARE}^{(0,\mathcal{P})}_{n}$,
${\sf ARE}^{(0,\mathcal{B})}_{n}$ criteria.

To analyze the quality of the estimates $\hat{\lambda}$ obtained from equation
(19), we consider an oracle loss estimator $\bm{\delta}^{\sf
or}_{(0)}\coloneqq{\bm{\delta}}_{(0)}^{\sf neb}(\lambda^{\sf orc}_{0})$ where

$\lambda^{\sf
orc}_{0}=\operatorname*{arg\,min}_{\lambda\in\Lambda}\mathcal{L}_{n}^{(0)}(\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}(\lambda))$

and Lemma 4 establishes the asymptotic optimality of $\hat{\lambda}$ obtained
from equation (19). In theorem 3B below we provide decision theoretic
guarantees on the NEB estimator and show that the largest coordinate-wise gap
between ${\bm{\delta}}^{\sf neb}_{(0)}({\hat{\lambda}})$ and
$\bm{\delta}_{(0)}^{\pi}$ is asymptotically small.

######  Theorem 3B.

Under the conditions of Theorem 2B, if
$\lim_{n\to\infty}c_{n}n^{-1/2}\log^{3}n=0$ then, for the Poisson and the
Binomial model,

$c_{n}\Big{\|}{\bm{\delta}}_{(0)}^{\sf
neb}({\hat{\lambda}})-\bm{\delta}_{(0)}^{\pi}\Big{\|}_{\infty}=o_{p}(1).$

Furthermore, under the same conditions, we have for the Poisson and the
Binomial model,

$\lim_{n\to\infty}\mathbb{P}\Big{[}\mathcal{L}_{n}^{(0)}\Bigl{\\{}\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}(\hat{\lambda})\Bigr{\\}}\geq\mathcal{L}_{n}^{(0)}(\bm{\theta},\bm{\delta}_{(0)}^{\pi})+c_{n}^{-1}\epsilon\Big{]}=0\text{
for any }\epsilon>0.$

##  Appendix B Technical Details and Proofs

We will begin this section with some notations and then state two lemmas that
will be used in proving the statements discussed in Section 3.

Let $c_{0},c_{1},\ldots$ denote some generic positive constants which may vary
in different statements. Let $\mathcal{D}_{n}=\\{0,1,2,\ldots,C\log n\\}$ and
given a random sample $(Y_{1},\ldots,Y_{n})$ from model (2) denote $B_{n}$ to
be the event $\\{\max_{1\leq i\leq n}Y_{i}\leq C\log n\\}$ where $C$ is the
constant given by lemma A below under assumption $(A2)$.

######  Lemma A.

Assumption $(A2)$ implies that with probability tending to $1$ as
$n\to\infty$,

$\max(Y_{1},\ldots,Y_{n})\leq C\log n$

where $C>0$ is a constant depending on $\epsilon$.

Our next lemma below is a statement on the pointwise Lipschitz stability of
the optimal solution $\hat{\bm{h}}_{n}^{(k)}(\lambda)$ under perturbations on
the parameter $\lambda\in\Lambda$. See, for example, Bonnans and Shapiro
(2013) for general results on the stability and sensitivity of parametrized
optimization problems.

######  Lemma B.

Let $\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})$ be the solution to problems (8) and
(15), respectively, for $k\in\\{0,1\\}$ and for some $\lambda_{0}\in\Lambda$.
Then, under Assumption $(A3)$, there exists a constant $L>0$ such that for any
$\lambda\in\Lambda$ the solution $\hat{\bm{h}}_{n}^{(k)}(\lambda)$ to problems
(8) and (15) satisfies

$\Big{\|}\hat{\bm{h}}_{n}^{(k)}(\lambda)-\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\Big{\|}_{2}\leq
L|\lambda-\lambda_{0}|$

###  B.1 Proof of Lemma 1

First note that for any coordinate $i$, the integrated Bayes risk of an
estimator $\delta_{(k),i}$ of $\theta_{i}$ is $\sum_{y_{i}}\int
p(y_{i}|\theta_{i})\ell_{n}^{(k)}(\theta_{i},\delta_{(k),i})dG(\theta_{i})$
which is minimized with respect to $\delta_{(k),i}$ if for each $y_{i}$,
$\delta_{(k),i}(y_{i})$ is defined as

$\delta_{(k),i}^{\pi}(y_{i})=\operatorname*{arg\,min}_{\delta_{(k),i}}\int
p(y_{i}|\theta_{i})\ell_{n}^{(k)}(\theta_{i},\delta_{(k),i})dG(\theta_{i})$

However, $\int
p(y_{i}|\theta_{i})\ell_{n}^{(k)}(\theta_{i},\delta_{(k),i})dG(\theta_{i})$ is
a minimum with respect to $\delta_{(k),i}$ when

$\delta_{(k),i}^{\pi}(y_{i})=\dfrac{\int
p(y_{i}|\theta_{i})\theta_{i}^{1-k}dG(\theta_{i})}{\int
p(y_{i}|\theta_{i})\theta_{i}^{-k}dG(\theta_{i})}$

The result then follows by noting that $p(y_{i}-k)=\int
a_{y_{i}-k}\theta_{i}^{y_{i}-k}/g(\theta_{i})dG(\theta_{i})$, and
$p(y_{i}+1-k)=\int
a_{y_{i}+1-k}\theta_{i}^{y_{i}+1-k}/g(\theta_{i})dG(\theta_{i})$ for
$~{}y_{i}=k,k+1,\ldots$.

###  B.2 Proof of Theorem 1

Define
$\tilde{\mathbb{M}}_{\lambda}(\bm{h})=\sum_{i,j\in\mathcal{D}_{n}}\kappa_{\lambda}[h(i),h(j)](i,j)\mathbb{P}(Y=i)\mathbb{P}(Y=j)$
and re-write $\hat{\mathbb{M}}_{\lambda,n}(\bm{h})$ as

$\hat{\mathbb{M}}_{\lambda,n}(\bm{h})=\dfrac{1}{n^{2}}\sum_{i,j\in\mathcal{D}_{n}}\kappa_{\lambda}[h(i),h(j)](i,j)\mathcal{C}_{ij},$

where $\mathcal{C}_{ij}$ is the number of pairs $(Y_{r},Y_{s})$ in the sample
that has $Y_{r}=i,Y_{s}=j$. Now, we have

$\sup_{\lambda\in\Lambda}\Big{|}\hat{\mathbb{M}}_{\lambda,n}(\bm{h})-{\mathbb{M}}_{\lambda}(\bm{h})\Big{|}\leq\sup_{\lambda\in\Lambda}\Big{|}\hat{\mathbb{M}}_{\lambda,n}(\bm{h})-\tilde{\mathbb{M}}_{\lambda}(\bm{h})\Big{|}+\sup_{\lambda\in\Lambda}\Big{|}{\mathbb{M}}_{\lambda}(\bm{h})-\tilde{\mathbb{M}}_{\lambda}(\bm{h})\Big{|}$
(20)

Consider the first term on the right hand side of the inequality in equation
(20) above and with $P_{i}\coloneqq\mathbb{P}(Y=i)$ note that assumption $A2$
and lemma A imply

$\displaystyle\mathbb{E}\sup_{\lambda\in\Lambda}\Big{|}\hat{\mathbb{M}}_{\lambda,n}(\bm{h})-\tilde{\mathbb{M}}_{\lambda}(\bm{h})\Big{|}\leq\sum_{i,j\in\mathcal{D}_{n}}\mathbb{E}\Big{[}\sup_{\lambda\in\Lambda}\Big{|}\kappa_{\lambda}[h(i),h(j)](i,j)\Big{|}\Big{|}\dfrac{\mathcal{C}_{ij}}{n^{2}}-P_{i}P_{j}\Big{|}\Big{]}\Bigl{\\{}1+o(1)\Bigr{\\}}$
$\displaystyle\leq\sum_{i,j\in\mathcal{D}_{n}}\Bigl{\\{}\mathbb{E}\Big{[}\sup_{\lambda\in\Lambda}\Big{|}\kappa_{\lambda}[h(i),h(j)](i,j)\Big{|}\Big{]}^{2}\mathbb{E}\Big{|}\dfrac{\mathcal{C}_{ij}}{n^{2}}-P_{i}P_{j}\Big{|}^{2}\Bigr{\\}}^{1/2}\Bigl{\\{}1+o(1)\Bigr{\\}}$
(21)

In equation (B.2) above, $\mathbb{E}|n^{-2}{\mathcal{C}_{ij}}-P_{i}P_{j}|^{2}$
is $O(1/n)$ and assumption $A1$ together with the compactness of $\Lambda$ and
the continuity of $\kappa_{\lambda}[h(i),h(j)](i,j)$ with respect to $\lambda$
imply that
$\mathbb{E}[\sup_{\lambda\in\Lambda}|\kappa_{\lambda}[h(i),h(j)](i,j)|]^{2}<\infty$.
Thus
$\mathbb{E}\sup_{\lambda\in\Lambda}|\hat{\mathbb{M}}_{\lambda,n}(\bm{h})-\tilde{\mathbb{M}}_{\lambda}(\bm{h})|$
is $O(\log^{2}n/n^{1/2})$. Now consider the second term on the right hand side
of the inequality in equation (20) and note that it is bounded above by the
following tail sums

$2\sum_{i\in\mathcal{D}_{n},j\notin\mathcal{D}_{n}}\sup_{\lambda\in\Lambda}|\kappa_{\lambda}[h(i),h(j)](i,j)|P_{i}P_{j}+\sum_{i,j\notin\mathcal{D}_{n}}\sup_{\lambda\in\Lambda}|\kappa_{\lambda}[h(i),h(j)](i,j)|P_{i}P_{j}.$

But from Assumption (A1),
$\mathbb{E}_{p}\sup_{\lambda\in\Lambda}|\kappa_{\lambda}[h(U),h(V)](U,V)|<\infty$
and together with assumption (A2) and proof of lemma A, it follows that the
terms in the display above are $O(n^{-\nu})$ for some $\nu>1/2$.

Now fix an $\epsilon>0$ and let $c_{n}=n^{1/2}/\log^{2}n$. Since
$\mathbb{E}_{p}\sup_{\lambda\in\Lambda}|\hat{\mathbb{M}}_{\lambda,n}(\bm{h})-\tilde{\mathbb{M}}_{\lambda}(\bm{h})|$
is $O(\log^{2}n/n^{1/2})$ there exists a finite constant $M>0$ and an $N_{1}$
such that
$c_{n}\mathbb{E}\sup_{\lambda\in\Lambda}|\hat{\mathbb{M}}_{\lambda,n}(\bm{h})-\tilde{\mathbb{M}}_{\lambda}(\bm{h})|\leq
M$ for all $n\geq N_{1}$. Moreover since
$\sup_{\lambda\in\Lambda}|{\mathbb{M}}_{\lambda}(\bm{h})-\tilde{\mathbb{M}}_{\lambda}(\bm{h})|\to
0$ as $n\to\infty$, there exists an $N_{2}$ such that
$\sup_{\lambda\in\Lambda}|{\mathbb{M}}_{\lambda}(\bm{h})-\tilde{\mathbb{M}}_{\lambda}(\bm{h})|\leq
M/c_{n}$ for all $n\geq N_{2}$. Thus with $t=4M/\epsilon$, we have
$\mathbb{P}(c_{n}\sup_{\lambda\in\Lambda}|\hat{\mathbb{M}}_{\lambda,n}(\bm{h})-{\mathbb{M}}_{\lambda}(\bm{h})|>t)<\epsilon$
for all $n\geq\max(N_{1},N_{2})$ which suffices to prove the desired result.

###  B.3 Proofs of Theorems 2A and 2B

We will first prove Theorem 2A. Note that from equation (7),

$\Big{\|}\hat{\bm{w}}_{n}^{(1)}(\lambda)-\bm{w}_{p}^{(1)}\Big{\|}_{2}=\Big{\|}\hat{\bm{h}}_{n}^{(1)}(\lambda)-\bm{h}_{0}^{(1)}\Big{\|}_{2}.$

Now from assumption $A3$ and for any $\epsilon>0$, there exists a $\delta>0$
such that for any $\lambda\in\Lambda$,

$\mathbb{P}\Big{[}c_{n}\Big{\|}\hat{\bm{h}}_{n}^{(1)}(\lambda)-\bm{h}_{0}^{(1)}\Big{\|}_{2}\geq\epsilon\Big{]}\leq\mathbb{P}\Big{[}c_{n}\Bigl{\\{}\mathbb{M}_{\lambda}(\hat{\bm{h}}_{n}^{(1)})-\mathbb{M}_{\lambda}(\bm{h}_{0}^{(1)})\Bigr{\\}}\geq\delta\Big{]}.$

But the right hand side is upper bounded by the sum of
$\mathbb{P}\Big{[}c_{n}\Bigl{\\{}\mathbb{M}_{\lambda}(\hat{\bm{h}}_{n}^{(1)})-\hat{\mathbb{M}}_{\lambda,n}(\hat{\bm{h}}_{n}^{(1)})\Bigr{\\}}\geq\delta/3\Big{]}$,
$\mathbb{P}\Big{[}c_{n}\Bigl{\\{}\hat{\mathbb{M}}_{\lambda,n}(\hat{\bm{h}}_{n}^{(1)})-\hat{\mathbb{M}}_{\lambda,n}({\bm{h}}_{0}^{(1)})\Bigr{\\}}\geq\delta/3\Big{]}$
and
$\mathbb{P}\Big{[}c_{n}\Bigl{\\{}\hat{\mathbb{M}}_{\lambda,n}({\bm{h}}_{0}^{(1)})-{\mathbb{M}}_{\lambda}({\bm{h}}_{0}^{(1)})\Bigr{\\}}\geq\delta/3\Big{]}$.
From theorem 1, the first and third terms go to zero as $n\to\infty$ while the
second term is zero since
$\hat{\mathbb{M}}_{\lambda,n}(\hat{\bm{h}}_{n}^{(1)})\leq\hat{\mathbb{M}}_{\lambda,n}({\bm{h}}_{0}^{(1)})$.
This proves the statement of theorem 2A. To prove theorem 2B first note that
from equation (14),

$\Big{\|}\hat{\bm{w}}_{n}^{(0)}(\lambda)-\bm{w}_{p}^{(0)}\Big{\|}_{2}^{2}\leq\sum_{i=1}^{n}\Bigl{\\{}\dfrac{Y_{i}+1}{\hat{h}_{n,i}^{(0)}(\lambda)h_{0,i}^{(0)}}\Bigr{\\}}^{2}\Bigl{\\{}\hat{h}_{n,i}^{(0)}(\lambda)-h_{0,i}^{(0)}\Bigr{\\}}^{2}.$

From assumption $A2$ and Lemma A, there exists a constant $c_{0}>0$ such that
for large $n$, $\max_{1\leq i\leq n}(Y_{i}+1)\leq c_{0}\log n$ with high
probability. Moreover for $i=1,\cdots,n$, since
$\hat{w}_{n,i}^{(0)}(\lambda)>0$ for every $\lambda\in\Lambda$ and
$w_{p,i}^{(0)}>0$, lemma A and equation (14) together imply
$\hat{h}_{n,i}^{(0)}(\lambda)+c_{0}\log n>0$ and ${h}_{0,i}^{(0)}+c_{0}\log
n>0$. Thus, conditional on the event $\\{\max_{1\leq i\leq n}(Y_{i}+1)\leq
c_{0}\log n\\}$ and for any $\epsilon>0$,

$\mathbb{P}\Big{[}c_{n}\log
n\Big{\|}\hat{\bm{w}}_{n}^{(0)}(\lambda)-\bm{w}_{p}^{(0)}\Big{\|}_{2}\geq\epsilon\Big{]}\leq\mathbb{P}\Big{[}c_{1}c_{n}\Big{\|}\hat{\bm{h}}_{n}^{(0)}(\lambda)-\bm{h}_{0}^{(0)}\Big{\|}_{2}\geq\epsilon\Big{]}$

for some constant $c_{1}>0$. The proof of the statement of theorem 2B thus
follows from the proof of theorem 2A above and lemma A.

###  B.4 Proof of Lemma 2

We will first prove the two statements of lemma 2 under the scaled squared
error loss. The proof for the squared error loss will follow from similar
arguments and we will highlight only the important steps. Throughout the
proof, we will denote $d_{1}\coloneqq\inf_{\lambda\in\Lambda}\inf_{1\leq i\leq
n}(1-\hat{h}_{n,i}^{(1)}(\lambda))>0$ and
$d_{2}\coloneqq\inf_{\lambda\in\Lambda}\inf_{1\leq i\leq
n}\hat{w}_{n,i}^{(0)}(\lambda)>0$.

Proof of statement 1 for the scaled squared error loss $(k=1)$

Note that by triangle inequality $\sup_{\lambda\in\Lambda}\Big{|}{\sf
ARE}^{(1,\mathcal{B})}_{n}(\lambda,\bm{Y})-\mathcal{R}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))\Big{|}$ is upper bounded by the following sum

$\sup_{\lambda\in\Lambda}\Big{|}{\sf
ARE}^{(1,\mathcal{B})}_{n}(\lambda,\bm{Y})-\mathbb{E}{\sf
ARE}^{(1,\mathcal{B})}_{n}(\lambda,\bm{Y})\Big{|}+\sup_{\lambda\in\Lambda}\Big{|}\mathbb{E}{\sf
ARE}^{(1,\mathcal{B})}_{n}(\lambda,\bm{Y})-\mathcal{R}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))\Big{|}.$

Consider the first term. Using definition 3A, this term is upper bounded by

$\displaystyle\Big{|}\dfrac{1}{n}\sum_{i=1}^{n}U_{i}\Big{|}$ $\displaystyle+$
$\displaystyle\sup_{\lambda\in\Lambda}\Big{|}\dfrac{2}{n}\sum_{i=1}^{n}\Bigl{\\{}{\delta}_{(1),i}^{\sf
neb}(\lambda)-\mathbb{E}{\delta}_{(1),i}^{\sf neb}(\lambda)\Bigr{\\}}\Big{|}$
(22) $\displaystyle+$
$\displaystyle\sup_{\lambda\in\Lambda}\Big{|}\dfrac{1}{n}\sum_{i=1}^{n}\Bigl{\\{}(m-Y_{i})\psi_{\lambda}(Y_{i})-\mathbb{E}(m-Y_{i})\psi_{\lambda}(Y_{i})\Bigr{\\}}\Big{|}$

where $U_{i}=Y_{i}/(m-Y_{i}+1)-\theta_{i}$, $\mathbb{E}U_{i}=0$ and
$\mathbb{E}U_{i}^{2}<\infty$ since $Y_{i}/(m-Y_{i}+1)\leq m<\infty$ for all
$i=1,\cdots,n$. So, $\Big{|}n^{-1}\sum_{i=1}^{n}U_{i}\Big{|}=O_{p}(n^{-1/2})$.

Now consider the second term in equation (22) above and define
$Z_{n}(\lambda)=n^{-1}\sum_{i=1}^{n}V_{i}(\lambda)$ where
$V_{i}(\lambda)={\delta}_{(1),i}^{\sf
neb}(\lambda)-\mathbb{E}{\delta}_{(1),i}^{\sf neb}(\lambda)$. Recall that for
every $i$ such that $Y_{i}>0$, $1-\hat{h}_{n,i}^{(1)}(\lambda)>0$ for all
$\lambda\in\Lambda$. Moreover, from definition 1A for the Binomial model,

${\delta}_{(1),i}^{\sf
neb}(\lambda)=[Y_{i}/(m-Y_{i}+1)]\left\\{1-\hat{h}_{n,i}^{(1)}(\lambda)\right\\}^{-1}.$

Thus, ${\delta}_{(1),i}^{\sf neb}(\lambda)\leq m/d_{1}$ and
$|V_{i}(\lambda)|\leq 2m/d_{1}$ for all $i=1,\ldots,n$. Along with the fact
that $\mathbb{E}V_{i}(\lambda)=0$ Hoeffding’s inequality gives, for a fixed
$\lambda$ and (for now) arbitrary $r_{n}>1$

$\mathbb{P}\Bigl{\\{}|Z_{n}(\lambda)|>\dfrac{r_{n}}{\sqrt{n}}\Bigr{\\}}\leq
2\exp\Bigr{\\{}-\dfrac{r_{n}^{2}d_{1}^{2}}{8m^{2}}\Bigl{\\}}$ (23)

Next for a perturbation $\lambda^{\prime}$ of $\lambda$ such that
$(\lambda,\lambda^{\prime})\in\Lambda\coloneqq[\lambda_{l},\lambda_{u}]$, we
will bound the increments
$\Big{|}Z_{n}(\lambda)-Z_{n}(\lambda^{\prime})\Big{|}$. To that effect, note
that

$n|Z_{n}(\lambda)-Z_{n}(\lambda^{\prime})|\leq\|{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda)-{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda^{\prime})\|_{1}+\mathbb{E}\|{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda)-{\bm{\delta}}_{(1)}^{\sf neb}(\lambda^{\prime})\|_{1}$

and $d_{1}^{2}\|{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda)-{\bm{\delta}}_{(1)}^{\sf neb}(\lambda^{\prime})\|_{1}\leq
m\|\hat{\bm{h}}^{(1)}_{n}(\lambda)-\hat{\bm{h}}^{(1)}_{n}(\lambda^{\prime})\|_{1}$.
Now from lemma B we know that

$\|\hat{\bm{h}}^{(1)}_{n}(\lambda)-\hat{\bm{h}}^{(1)}_{n}(\lambda^{\prime})\|_{1}\leq
n^{1/2}c^{-1}|\lambda-\lambda^{\prime}|\sup_{{\bm{h}}\in
N_{\delta}(\hat{\bm{h}}_{n}^{(1)}(\lambda))}\|\nabla^{2}_{\bm{h}_{n},\lambda}\hat{\mathbb{M}}_{\lambda,n}({\bm{h}})+o(1)\|_{2}.$

Moreover, the Binomial model with $m<\infty$ and assumption $A4$ imply that
the supremum in the display above is $O(\log n)$. Thus so long as
$|\lambda-\lambda^{\prime}|\leq\epsilon_{n}$,

$\Big{|}Z_{n}(\lambda)-Z_{n}(\lambda^{\prime})\Big{|}\leq
c_{0}\epsilon_{n}\dfrac{\log n}{\sqrt{n}}.$

Now choose $\lambda_{j}=j\epsilon_{n}\in\Lambda$ and note that
$A_{n}=\\{\sup_{\lambda\in\Lambda}|Z_{n}(\lambda)|>3r_{n}/\sqrt{n}\\}\subset
D_{n}\cup E_{n}$ where

$D_{n}=\\{\sup_{j}|Z_{n}(\lambda_{j})|>r_{n}/\sqrt{n}\\}\text{ and
}E_{n}=\\{\sup_{j}\sup_{|\lambda-\lambda_{j}|\leq\epsilon_{n}}|Z_{n}(\lambda)-Z_{n}(\lambda_{j})|>2r_{n}/\sqrt{n}\\}.$

Choose $\epsilon_{n}$ so that $\epsilon_{n}n^{-1/2}\log n=o(n^{-1/2})$ and
note that $\mathbb{P}(E_{n})\leq\exp\\{-2r_{n}^{2}\\}$ and
$\mathbb{P}(D_{n})\leq
2(\lambda_{u}/\epsilon_{n})\exp\\{-r_{n}^{2}d_{1}^{2}/8m^{2}\\}$ from equation
(23) and the cardinality of $\lambda_{j}$. Thus,

$P(A_{n})\leq
2(\lambda_{u}/\epsilon_{n})\exp\\{-r_{n}^{2}d_{1}^{2}/8m^{2}\\}+\exp\\{-2r_{n}^{2}\\}$

Set $r_{n}=(s\log n)^{1/2}(m\sqrt{8}/d_{1})=O(\sqrt{\log n})$. Then
$\mathbb{P}(A_{n})\leq n^{-s}3(\lambda_{u}/\epsilon_{n})$ and thus the second
term in equation (22) is $O_{p}(\sqrt{\log n/n})$.

We will now consider the third term in equation (22) and analyze it in a
similar manner to the second term of equation (22). Here we will assume that
the set $\mathcal{I}_{i}=\\{j:Y_{j}=Y_{i}+1\\}$ is non-empty for every
$i=1,\ldots,n$. Recall from definition 3A that
$(m-Y_{i})\psi_{\lambda}(Y_{i})\leq m/d_{1}^{2}$ for all $Y_{i}=0,1,\ldots,m$.
Define $Z_{n}(\lambda)=n^{-1}\sum_{i=1}^{n}V_{i}(\lambda)$ where
$V_{i}(\lambda)=(m-Y_{i})\psi_{\lambda}(Y_{i})-\mathbb{E}(m-Y_{i})\psi_{\lambda}(Y_{i})$
with $\mathbb{E}V_{i}(\lambda)=0$ and $|V_{i}(\lambda)|\leq 2m/d_{1}^{2}$.
Hoeffding’s inequality gives, for a fixed $\lambda$ and (for now) arbitrary
$r_{n}>1$

$\mathbb{P}\Bigl{\\{}|Z_{n}(\lambda)|>\dfrac{r_{n}}{\sqrt{n}}\Bigr{\\}}\leq
2\exp\Bigr{\\{}-\dfrac{r_{n}^{2}d_{1}^{4}}{8m^{2}}\Bigl{\\}}$

Moreover, for a perturbation $\lambda^{\prime}$ of $\lambda$ such that
$(\lambda,\lambda^{\prime})\in\Lambda$, the increments
$n|Z_{n}(\lambda)-Z_{n}(\lambda^{\prime})|$ are bounded above by
$m\|\bm{\psi}_{\lambda}-\bm{\psi}_{\lambda^{\prime}}\|_{1}+m\mathbb{E}\|\bm{\psi}_{\lambda}-\bm{\psi}_{\lambda^{\prime}}\|_{1}$
where
$\bm{\psi}_{\lambda}=(\psi_{\lambda}(Y_{1}),\ldots,\psi_{\lambda}(Y_{n}))$ for
any $\lambda\in\Lambda$. Now, note that from definition 3A and for the
Binomial model,

$|\psi_{\lambda}(Y_{i})-\psi_{\lambda^{\prime}}(Y_{i})|\leq
2(m^{2}/d_{1}^{3})|\hat{h}_{n,j}(\lambda)-\hat{h}_{n,j}(\lambda^{\prime})|,$

where $j\in\mathcal{I}_{i}.$ Therefore

$\displaystyle|Z_{n}(\lambda)-Z_{n}(\lambda^{\prime})|\leq
2\dfrac{m^{3}}{nd_{1}^{3}}\Big{[}\|\hat{\bm{h}}^{(1)}_{n}(\lambda)-\hat{\bm{h}}^{(1)}_{n}(\lambda^{\prime})\|_{1}+\mathbb{E}\|\hat{\bm{h}}^{(1)}_{n}(\lambda)-\hat{\bm{h}}^{(1)}_{n}(\lambda^{\prime})\|_{1}\Big{]}.$

Now lemma B and assumption $A4$ imply that the right hand side of the
inequality above is $|\lambda-\lambda^{\prime}|O(n^{-1/2}\log n)$. So as long
as $|\lambda-\lambda^{\prime}|\leq\epsilon_{n}$, choose
$\lambda_{j}=j\epsilon_{n}\in\Lambda$. In a manner similar to the second term
of equation (22) define the events $A_{n},~{}D_{n}$ and $E_{n}$, and set
$r_{n}=(s\log n)^{1/2}(m\sqrt{8}/d_{1}^{2})$ to conclude that the third term
of equation (22) continues to be is $O_{p}(\sqrt{\log n/n})$ which suffices to
prove the statement of the result.

Proof of statement 2 for the scaled squared error loss $(k=1)$

Note that $\sup_{\lambda\in\Lambda}|{\sf
ARE}^{(1,\mathcal{B})}_{n}(\lambda,\bm{Y})-\mathcal{L}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))|$ is bounded above by the sum of:
$\sup_{\lambda\in\Lambda}|{\sf
ARE}^{(1,\mathcal{B})}_{n}(\lambda,\bm{Y})-\mathcal{R}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))|$ and
$\sup_{\lambda\in\Lambda}|\mathcal{L}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))-\mathbb{E}\mathcal{L}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))|$. The first term is $O_{p}(\sqrt{\log n/n})$ from statement 1.
The second term is bounded above by

$\dfrac{2}{n}\sup_{\lambda\in\Lambda}\Big{|}\sum_{i=1}^{n}\Bigl{\\{}{\delta}_{(1),i}^{\sf
neb}(\lambda)-\mathbb{E}{\delta}_{(1),i}^{\sf
neb}(\lambda)\Bigr{\\}}\Big{|}+\dfrac{1}{n}\sup_{\lambda\in\Lambda}\Big{|}\sum_{i=1}^{n}\theta_{i}^{-1}\Bigl{\\{}\Big{[}{\delta}_{(1),i}^{\sf
neb}(\lambda)\Big{]}^{2}-\mathbb{E}\Big{[}{\delta}_{(1),i}^{\sf
neb}(\lambda)\Big{]}^{2}\Bigr{\\}}\Big{|}$ (24)

where the first term in equation (24) is $O_{p}(\sqrt{\log n/n})$ from the
proof of statement 1. Now consider the second term in equation (24) and define
$Z_{n}(\lambda)=n^{-1}\sum_{i=1}^{n}V_{i}(\lambda)$ where
$\theta_{i}V_{i}(\lambda)=[{\delta}_{(1),i}^{\sf
neb}(\lambda)]^{2}-\mathbb{E}[{\delta}_{(1),i}^{\sf neb}(\lambda)]^{2}$.
Recall that from definition 1A for the Binomial model,

${\delta}_{(1),i}^{\sf
neb}(\lambda)=[Y_{i}/(m-Y_{i}+1)]\left\\{1-\hat{h}_{n,i}^{(1)}(\lambda)\right\\}^{-1}.$

Thus, $[{\delta}_{(1),i}^{\sf neb}(\lambda)]^{2}\leq m^{2}/d_{1}^{2}$ and
$|V_{i}(\lambda)|\leq 2\theta_{i}^{-1}m^{2}/d_{1}^{2}$ for all $i=1,\ldots,n$.
For an arbitrary $r_{n}>1$ and $\lambda$ fixed, we have from Hoeffding’s
inequality,

$\mathbb{P}\Bigl{\\{}|Z_{n}(\lambda)|>\dfrac{r_{n}}{\sqrt{n}}\Bigr{\\}}\leq
2\exp\Bigr{\\{}-\dfrac{r_{n}^{2}d_{1}^{4}n}{8m^{2}\sum_{i=1}^{n}\theta_{i}^{-1}}\Bigl{\\}}$
(25)

Moreover for a perturbation $\lambda^{\prime}$ of $\lambda$ such that
$(\lambda,\lambda^{\prime})\in\Lambda$,

$\sum_{i=1}^{n}\theta_{i}^{-1}\Big{|}[{\delta}_{(1),i}^{\sf
neb}(\lambda)]^{2}-[{\delta}_{(1),i}^{\sf
neb}(\lambda^{\prime})]^{2}\Big{|}\leq\dfrac{2m^{2}}{d_{1}^{2}}\Big{\|}\hat{\bm{h}}^{(1)}_{n}(\lambda)-\hat{\bm{h}}^{(1)}_{n}(\lambda^{\prime})\Big{\|}_{2}\sum_{i=1}^{n}\theta_{i}^{-1}.$

Thus lemma B, assumption $A4$ and the display above together imply that
$n|Z_{n}(\lambda)-Z_{n}(\lambda^{\prime})|$ is bounded above by
$c_{1}|\lambda-\lambda^{\prime}|\log n\sum_{i=1}^{n}\theta_{i}^{-1}$. Now as
long as $|\lambda-\lambda^{\prime}|\leq\epsilon_{n}$, choose $\epsilon_{n}$ so
that $\epsilon_{n}n^{-1}\log n\sum_{i=1}^{n}\theta_{i}^{-1}=o(n^{-1/2})$ and
along with equation (25), follow the steps outlined in the proof of the second
term in equation (22) to conclude that the second term in equation (24) is
$O_{p}(\sqrt{\log n/n})$ from which the desired result follows.

Proof of statement 1 for the squared error loss $(k=0)$

The proof of this statement is very similar to the proof of statement 1 under
the scaled squared error loss and therefore we highlight the important steps
here. To prove statement 1, we will only look at the term
$\sup_{\lambda\in\Lambda}\Big{|}{\sf
ARE}^{(0,\mathcal{B})}_{n}(\lambda,\bm{Y})-\mathbb{E}{\sf
ARE}^{(0,\mathcal{B})}_{n}(\lambda,\bm{Y})\Big{|}$ because under the Binomial
model, it can be verified using definition 3B that $\mathbb{E}{\sf
ARE}^{(0,\mathcal{B})}_{n}(\lambda)=\mathcal{R}_{n}^{(0)}(\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}(\lambda))$. Now note that $\sup_{\lambda\in\Lambda}\Big{|}{\sf
ARE}^{(0,\mathcal{B})}_{n}(\lambda,\bm{Y})-\mathbb{E}{\sf
ARE}^{(0,\mathcal{B})}_{n}(\lambda,\bm{Y})\Big{|}$ is bounded above by

$\displaystyle\Big{|}\dfrac{1}{n}\sum_{i=1}^{n}U_{i}\Big{|}$ $\displaystyle+$
$\displaystyle\sup_{\lambda\in\Lambda}\Big{|}\dfrac{2}{n}\sum_{i=1}^{n}\Bigl{\\{}Y_{i}\psi_{\lambda}(Y_{i})-\mathbb{E}Y_{i}\psi_{\lambda}(Y_{i})\\}\Big{|}$
(26) $\displaystyle+$
$\displaystyle\dfrac{1}{n}\sup_{\lambda\in\Lambda}\Big{|}\sum_{i=1}^{n}\Bigl{\\{}\Big{[}{\delta}_{(0),i}^{\sf
neb}(\lambda)\Big{]}^{2}-\mathbb{E}\Big{[}{\delta}_{(0),i}^{\sf
neb}(\lambda)\Big{]}^{2}\Bigr{\\}}\Big{|}$

where $U_{i}=Y_{i}(Y_{i}-1)/[(m-Y_{i}+2)(m-Y_{i}+1)]-\theta_{i}^{2}$,
$\mathbb{E}U_{i}=0$ and $\mathbb{E}U_{i}^{2}<\infty$ since $|U_{i}|\leq
m^{2}<\infty$ for all $i=1,\cdots,n$. So,
$\Big{|}n^{-1}\sum_{i=1}^{n}U_{i}\Big{|}=O_{p}(n^{-1/2})$. For the second term
in equation (26) note that from definition 3B, $Y_{i}\psi_{\lambda}(Y_{i})\leq
m^{2}/d_{2}$ and for a perturbation $\lambda^{\prime}$ of $\lambda$ such that
$(\lambda,\lambda^{\prime})\in\Lambda$,

$\sum_{i=1}^{n}Y_{i}\Big{|}\psi_{\lambda}(Y_{i})-\psi_{\lambda^{\prime}}(Y_{i})\Big{|}\leq
c_{2}|\lambda-\lambda^{\prime}|\sqrt{n}\log n\\{1+o(1)\\}.$

The last inequality in the display above follows from definition 3B, lemma B
and assumption $A4$. Thus the upper bound on the second term of equation (26)
and the corresponding upper bound on its increments over
$(\lambda,\lambda^{\prime})\in\Lambda$ suffice to show that this term is
$O_{p}(\sqrt{\log n/n})$. Finally, the third term in equation (26) is bounded
above by $4m^{2}/d_{2}^{2}$ and $\|{\bm{\delta}}_{(0)}^{\sf
neb}(\lambda)-{\bm{\delta}}_{(0)}^{\sf neb}(\lambda^{\prime})\|_{1}$ is
$|\lambda-\lambda^{\prime}|O(\sqrt{n}\log n)$ for
$(\lambda,\lambda^{\prime})\in\Lambda$ from which the desired follows that
$\sup_{\lambda\in\Lambda}\Big{|}{\sf
ARE}^{(0,\mathcal{B})}_{n}(\lambda,\bm{Y})-\mathbb{E}{\sf
ARE}^{(0,\mathcal{B})}_{n}(\lambda,\bm{Y})\Big{|}$ is $O_{p}(\sqrt{\log
n/n})$.

Proof of statement 2 for the squared error loss $(k=0)$

We will only look at the term
$\sup_{\lambda\in\Lambda}|\mathcal{L}_{n}^{(0)}(\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}(\lambda))-\mathbb{E}\mathcal{L}_{n}^{(0)}(\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}(\lambda))|$ and show that it is $O_{p}(\sqrt{\log n/n})$. Note that

$\dfrac{2}{n}\sup_{\lambda\in\Lambda}\Big{|}\sum_{i=1}^{n}\theta_{i}\Bigl{\\{}{\delta}_{(0),i}^{\sf
neb}(\lambda)-\mathbb{E}{\delta}_{(0),i}^{\sf
neb}(\lambda)\Bigr{\\}}\Big{|}+\dfrac{1}{n}\sup_{\lambda\in\Lambda}\Big{|}\sum_{i=1}^{n}\Bigl{\\{}\Big{[}{\delta}_{(0),i}^{\sf
neb}(\lambda)\Big{]}^{2}-\mathbb{E}\Big{[}{\delta}_{(0),i}^{\sf
neb}(\lambda)\Big{]}^{2}\Bigr{\\}}\Big{|}$ (27)

is an upper bound on
$\sup_{\lambda\in\Lambda}|\mathcal{L}_{n}^{(0)}(\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}(\lambda))-\mathbb{E}\mathcal{L}_{n}^{(0)}(\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}(\lambda))|$. Analogous to the preceding proofs under lemma 2, we will
provide the upper bounds and bounds on the increments with respect to
perturbations on $\lambda$ for the two terms in equation (27). The rest of the
proof will then follow from the proof of statement 1 in lemma 2.

From definition 1B, we know that ${\delta}_{(0),i}^{\sf neb}(\lambda)\leq
2m/d_{2}$. Thus

$\sum_{i=1}^{n}|\theta_{i}{\delta}_{(0),i}^{\sf neb}(\lambda)|\leq
2(m/d_{2})\sum_{i=1}^{n}\theta_{i}.$

Moreover, from Lemma B and Assumption (A4),

$\sum_{i=1}^{n}\theta_{i}\Big{|}{\delta}_{(0),i}^{\sf
neb}(\lambda)-{\delta}_{(0),i}^{\sf neb}(\lambda^{\prime})\Big{|}\leq
c_{3}|\lambda-\lambda^{\prime}|\sum_{i=1}^{n}\theta_{i}\log n\\{1+o(1)\\}$

and

$\sum_{i=1}^{n}\Big{|}\Big{[}{\delta}_{(0),i}^{\sf
neb}(\lambda)\Big{]}^{2}-[{\delta}_{(0),i}^{\sf
neb}(\lambda^{\prime})\Big{]}^{2}\Big{|}\leq
c_{4}|\lambda-\lambda^{\prime}|\sqrt{n}\log n\\{1+o(1)\\}$

for $(\lambda,\lambda^{\prime})\in\Lambda$. Hoeffding’s inequality and the
steps outlined in the proof of statement 1 ($k=1$) of lemma 2 will then show
that both the terms in equation (27) are $O_{p}(\sqrt{\log n/n})$ which proves
the desired result.

###  B.5 Proof of Lemma 3

We will first prove the two statements of lemma 3 under the scaled squared
error loss. The proof for the squared error loss will follow from similar
arguments and we will highlight only the important steps. Throughout the
proof, we will denote $d_{1}\coloneqq\inf_{\lambda\in\Lambda}\inf_{1\leq i\leq
n}(1-\hat{h}_{n,i}^{(1)}(\lambda))>0$ and
$d_{2}\coloneqq\inf_{\lambda\in\Lambda}\inf_{1\leq i\leq
n}\hat{w}_{n,i}^{(0)}(\lambda)>0$.

Proof of statement 1 for the scaled squared error loss $(k=1)$

Note that by triangle inequality $\sup_{\lambda\in\Lambda}\Big{|}{\sf
ARE}^{(1,\mathcal{P})}_{n}(\lambda,\bm{Y})-\mathcal{R}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))\Big{|}$ is upper bounded by the following sum:

$\sup_{\lambda\in\Lambda}\Big{|}{\sf
ARE}^{(1,\mathcal{P})}_{n}(\lambda,\bm{Y})-\mathbb{E}{\sf
ARE}^{(1,\mathcal{P})}_{n}(\lambda,\bm{Y})\Big{|}+\sup_{\lambda\in\Lambda}\Big{|}\mathbb{E}{\sf
ARE}^{(1,\mathcal{P})}_{n}(\lambda,\bm{Y})-\mathcal{R}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))\Big{|}.$

Under the Poisson model and definition 2A it can be shown that

$\mathbb{E}\left\\{{\sf
ARE}^{(1,\mathcal{P})}_{n}(\lambda,\bm{Y})\right\\}=\mathcal{R}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda)).$

So the second term in the display above is zero. Now consider the first term
and note that it is bounded above by

$\displaystyle\Big{|}\dfrac{1}{n}\sum_{i=1}^{n}U_{i}\Big{|}+\sup_{\lambda\in\Lambda}\Big{|}\dfrac{2}{n}\sum_{i=1}^{n}\Bigl{\\{}{\delta}_{(1),i}^{\sf
neb}(\lambda)-\mathbb{E}{\delta}_{(1),i}^{\sf
neb}(\lambda)\Bigr{\\}}\Big{|}+\sup_{\lambda\in\Lambda}\Big{|}\dfrac{1}{n}\sum_{i=1}^{n}\Bigl{\\{}\psi_{\lambda}(Y_{i})-\mathbb{E}\psi_{\lambda}(Y_{i})\Bigr{\\}}\Big{|}$
(28)

where $U_{i}=Y_{i}-\theta_{i}$ and $\mathbb{E}U_{i}=0$. Then using Hoeffding’s
inequality on $\Big{|}n^{-1}\sum_{i=1}^{n}U_{i}\Big{|}$ along with assumption
$A3$ and lemma A gives that the first term in equation (28) is
$O_{p}(\log^{3/2}n/\sqrt{n})$. Now consider the second term in equation (28)
and define $V_{i}(\lambda)={\delta}_{(1),i}^{\sf
neb}(\lambda)-\mathbb{E}{\delta}_{(1),i}^{\sf neb}(\lambda)$ with
$Z_{n}(\lambda)=n^{-1}\sum_{i=1}^{n}V_{i}(\lambda)$. Conditional on the event
$B_{n}$, we have $|V_{i}(\lambda)|\leq 2C\log n/d_{1}$. Moreover conditional
on $B_{n}$, assumptions $A2-A4$ and lemma B give

$\Big{|}Z_{n}(\lambda)-Z_{n}(\lambda^{\prime})\Big{|}\leq
c_{0}\epsilon_{n}\dfrac{\log^{4}n}{\sqrt{n}}\\{1+o(1)\\}$

whenever $(\lambda,\lambda^{\prime})\in\Lambda$ and
$|\lambda-\lambda^{\prime}|\leq\epsilon_{n}$. The proof of statement 1 for
lemma 2 and the bounds in the display above along with those developed for
$V_{i}(\lambda)$ establish that the second term in equation (28) is
$O_{p}(\log^{3/2}n/\sqrt{n})$. Now for the third term in equation (28), define
$V_{i}(\lambda)=\psi_{\lambda}(Y_{i})-\mathbb{E}\psi_{\lambda}(Y_{i})$. From
assumption $A3$ and lemma A, there exists a constant $C^{\prime}>0$ such that
for large $n$, $\max_{1\leq i\leq n}(Y_{i}+1)\leq C^{\prime}\log n$ with high
probability which gives, conditional on the event
$B_{n}^{\prime}=\\{\max_{1\leq i\leq n}(Y_{i}+1)\leq C^{\prime}\log n\\}$,
$|V_{i}(\lambda)|\leq c_{1}\log^{2}n$. Moreover with
$Z_{n}(\lambda)=n^{-1}\sum_{i=1}^{n}V_{i}(\lambda)$ and conditional on the
event $B_{n}^{\prime}$, assumptions $A2-A4$ and lemma B give

$\Big{|}Z_{n}(\lambda)-Z_{n}(\lambda^{\prime})\Big{|}\leq
c_{2}\epsilon_{n}\dfrac{\log^{5}n}{\sqrt{n}}\\{1+o(1)\\}.$

Now we mimic the proof of statement 1 for lemma 2 to establish that the third
term in equation (28) is $O_{p}(\log^{5/2}n/\sqrt{n})$ which proves the
statement of the lemma.

Proof of statement 2 for the scaled squared error loss $(k=1)$

For the proof of this statement, we will show that
$\sup_{\lambda\in\Lambda}|\mathcal{L}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))-\mathbb{E}\mathcal{L}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}(\lambda))|$ is $O_{p}(\log^{5/2}n/\sqrt{n})$. This term is bounded above
by

$\dfrac{2}{n}\sup_{\lambda\in\Lambda}\Big{|}\sum_{i=1}^{n}\Bigl{\\{}{\delta}_{(1),i}^{\sf
neb}(\lambda)-\mathbb{E}{\delta}_{(1),i}^{\sf
neb}(\lambda)\Bigr{\\}}\Big{|}+\dfrac{1}{n}\sup_{\lambda\in\Lambda}\Big{|}\sum_{i=1}^{n}\theta_{i}^{-1}\Bigl{\\{}\Big{[}{\delta}_{(1),i}^{\sf
neb}(\lambda)\Big{]}^{2}-\mathbb{E}\Big{[}{\delta}_{(1),i}^{\sf
neb}(\lambda)\Big{]}^{2}\Bigr{\\}}\Big{|}$ (29)

where the first term in equation (29) is $O_{p}(\log^{3/2}n/\sqrt{n})$ from
the proof of statement 1. Now consider the second term in equation (29) and
define $Z_{n}(\lambda)=n^{-1}\sum_{i=1}^{n}V_{i}(\lambda)$ where
$\theta_{i}V_{i}(\lambda)=[{\delta}_{(1),i}^{\sf
neb}(\lambda)]^{2}-\mathbb{E}[{\delta}_{(1),i}^{\sf neb}(\lambda)]^{2}$.
Conditional on the event $B_{n}$, we have $|V_{i}(\lambda)|\leq
2\theta_{i}^{-1}(C/d_{1})^{2}\log^{2}n$. Moreover for a perturbation
$\lambda^{\prime}$ of $\lambda$ such that
$(\lambda,\lambda^{\prime})\in\Lambda$,

$\sum_{i=1}^{n}\theta_{i}^{-1}\Big{|}[{\delta}_{(1),i}^{\sf
neb}(\lambda)]^{2}-[{\delta}_{(1),i}^{\sf
neb}(\lambda^{\prime})]^{2}\Big{|}\leq\dfrac{2C^{2}\log^{2}n}{d_{1}^{3}}\Big{\|}\hat{\bm{h}}^{(1)}_{n}(\lambda)-\hat{\bm{h}}^{(1)}_{n}(\lambda^{\prime})\Big{\|}_{2}\sum_{i=1}^{n}\theta_{i}^{-1}$

conditional on $B_{n}$. Thus assumptions $A2-A4$, lemma B and the above
display together imply that $n|Z_{n}(\lambda)-Z_{n}(\lambda^{\prime})|$ is
bounded above by
$c_{0}\epsilon_{n}\log^{5}n\sum_{i=1}^{n}\theta_{i}^{-1}\\{1+o(1)\\}$ whenever
$|\lambda-\lambda^{\prime}|\leq\epsilon_{n}$. Now we follow the steps outlined
in the proof of statement 1 for lemma 3 to conclude that the second term in
equation (29) is $O_{p}(\log^{5/2}n/\sqrt{n})$ from which the desired result
follows.

Proof of statement 1 for the squared error loss $(k=0)$

The proof of this statement is very similar to the proof of statement 1 under
the scaled squared error loss and therefore we highlight the important steps
here. To prove statement 1, we will only show that the term
$\sup_{\lambda\in\Lambda}\Big{|}{\sf
ARE}^{(0,\mathcal{P})}_{n}(\lambda,\bm{Y})-\mathbb{E}{\sf
ARE}^{(0,\mathcal{P})}_{n}(\lambda,\bm{Y})\Big{|}$ is $O_{p}(\log
n^{5/2}/\sqrt{n})$ because under the Poisson model, it can be verified using
definition 2B that $\mathbb{E}{\sf
ARE}^{(0,\mathcal{P})}_{n}(\lambda)=\mathcal{R}_{n}^{(0)}(\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}(\lambda))$. Now note that $\sup_{\lambda\in\Lambda}\Big{|}{\sf
ARE}^{(0,\mathcal{P})}_{n}(\lambda,\bm{Y})-\mathbb{E}{\sf
ARE}^{(0,\mathcal{P})}_{n}(\lambda,\bm{Y})\Big{|}$ is bounded above by

$\displaystyle\Big{|}\dfrac{1}{n}\sum_{i=1}^{n}U_{i}\Big{|}$ $\displaystyle+$
$\displaystyle\sup_{\lambda\in\Lambda}\Big{|}\dfrac{2}{n}\sum_{i=1}^{n}\Bigl{\\{}Y_{i}\psi_{\lambda}(Y_{i})-\mathbb{E}Y_{i}\psi_{\lambda}(Y_{i})\\}\Big{|}$
(30) $\displaystyle+$
$\displaystyle\dfrac{1}{n}\sup_{\lambda\in\Lambda}\Big{|}\sum_{i=1}^{n}\Bigl{\\{}\Big{[}{\delta}_{(0),i}^{\sf
neb}(\lambda)\Big{]}^{2}-\mathbb{E}\Big{[}{\delta}_{(0),i}^{\sf
neb}(\lambda)\Big{]}^{2}\Bigr{\\}}\Big{|}$

where $U_{i}=Y_{i}(Y_{i}-1)-\theta_{i}^{2}$ and $\mathbb{E}U_{i}=0$. Then
using Hoeffding’s inequality on $\Big{|}n^{-1}\sum_{i=1}^{n}U_{i}\Big{|}$,
along with assumption $A3$ and lemma A, gives that the first term in equation
(30) is $O_{p}(\log^{5/2}n/\sqrt{n})$. Now for the second term in equation
(30), define
$V_{i}(\lambda)=Y_{i}\psi_{\lambda}(Y_{i})-\mathbb{E}Y_{i}\psi_{\lambda}(Y_{i})$.
Conditional on the event $B_{n}$, $|V_{i}(\lambda)|\leq c_{0}\log^{2}n$.
Moreover with $Z_{n}(\lambda)=n^{-1}\sum_{i=1}^{n}V_{i}(\lambda)$ and
conditional on the event $B_{n}$, assumptions $A2-A4$ and lemma B give

$|Z_{n}(\lambda)-Z_{n}(\lambda^{\prime})|\leq
c_{1}\epsilon_{n}{n^{-1/2}\log^{5}n}\\{1+o(1)\\}$

whenever $|\lambda-\lambda^{\prime}|\leq\epsilon_{n}$ for
$(\lambda,\lambda^{\prime})\in\Lambda$. The proof of statement 1 ($k=1$ case)
for lemma 3 and the bounds in the display above along with those developed for
$V_{i}(\lambda)$ establish that the second term in equation (30) is
$O_{p}(\log^{5/2}n/\sqrt{n})$. For the third term in equation (30), we proceed
in a similar manner and define $V_{i}(\lambda)=[{\delta}_{(0),i}^{\sf
neb}(\lambda)]^{2}-\mathbb{E}[{\delta}_{(0),i}^{\sf neb}(\lambda)]^{2}$ and
$Z_{n}(\lambda)=n^{-1}\sum_{i=1}^{n}V_{i}(\lambda)$. From Assumption (A3) and
Lemma A, there exists a constant $C^{\prime}>0$ such that for large $n$,
$\max_{1\leq i\leq n}(Y_{i}+1)\leq C^{\prime}\log n$ with high probability
which gives, conditional on the event $B_{n}^{\prime}=\\{\max_{1\leq i\leq
n}(Y_{i}+1)\leq C^{\prime}\log n\\}$, (i) $|V_{i}(\lambda)|\leq
c_{2}\log^{2}n$, and (ii) under Assumptions (A2)-(A4) and Lemma B,

$|Z_{n}(\lambda)-Z_{n}(\lambda^{\prime})|\leq
c_{3}\epsilon_{n}n^{-1/2}\log^{4}n\\{1+o(1)\\}$

whenever $|\lambda-\lambda^{\prime}|\leq\epsilon_{n}$ for
$(\lambda,\lambda^{\prime})\in\Lambda$. Thus the third term in equation (30)
is $O_{p}(\log^{5/2}/\sqrt{n})$ which follows from the proof of statement 1
($k=1$ case) for lemma 3 together with the preceding bounds developed for
$V_{i}(\lambda)$ and $|Z_{n}(\lambda)-Z_{n}(\lambda^{\prime})|$, and suffices
to prove the desired result.

Proof of statement 2 for the squared error loss $(k=0)$

We will only look at the term
$\sup_{\lambda\in\Lambda}|\mathcal{L}_{n}^{(0)}(\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}(\lambda))-\mathbb{E}\mathcal{L}_{n}^{(0)}(\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}(\lambda))|$ and show that it is $O_{p}(\log^{5/2}n/\sqrt{n})$. Note that

$\dfrac{2}{n}\sup_{\lambda\in\Lambda}\Big{|}\sum_{i=1}^{n}\theta_{i}\Bigl{\\{}{\delta}_{(0),i}^{\sf
neb}(\lambda)-\mathbb{E}{\delta}_{(0),i}^{\sf
neb}(\lambda)\Bigr{\\}}\Big{|}+\dfrac{1}{n}\sup_{\lambda\in\Lambda}\Big{|}\sum_{i=1}^{n}\Bigl{\\{}\Big{[}{\delta}_{(0),i}^{\sf
neb}(\lambda)\Big{]}^{2}-\mathbb{E}\Big{[}{\delta}_{(0),i}^{\sf
neb}(\lambda)\Big{]}^{2}\Bigr{\\}}\Big{|}$ (31)

is an upper bound on
$\sup_{\lambda\in\Lambda}|\mathcal{L}_{n}^{(0)}(\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}(\lambda))-\mathbb{E}\mathcal{L}_{n}^{(0)}(\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}(\lambda))|$. Analogous to the preceding proofs under lemma 3, we will
provide the upper bounds and bounds on the increments with respect to
perturbations on $\lambda$ for the two terms in equation (31). The rest of the
proof will then follow from the proof of statement 1 in lemma 3.

From definition 1B, we know that under the Poisson model and conditional on
the event $B_{n}^{\prime}$, ${\delta}_{(0),i}^{\sf neb}(\lambda)\leq
C^{\prime}\log n/d_{2}$. Moreover, from lemma B,
$\sum_{i=1}^{n}\theta_{i}|{\delta}_{(0),i}^{\sf
neb}(\lambda)-{\delta}_{(0),i}^{\sf neb}(\lambda^{\prime})|$ is bounded above
by
$c_{0}|\lambda-\lambda^{\prime}|\sum_{i=1}^{n}\theta_{i}\log^{3}n\\{1+o(1)\\}$
and $\sum_{i=1}^{n}|[{\delta}_{(0),i}^{\sf
neb}(\lambda)]^{2}-[{\delta}_{(0),i}^{\sf neb}(\lambda^{\prime})]^{2}|$ is
bounded above by $c_{1}|\lambda-\lambda^{\prime}|n^{1/2}\log^{4}n\\{1+o(1)\\}$
for $(\lambda,\lambda^{\prime})\in\Lambda$. Hoeffding’s inequality and the
steps outlined in the proof of statement 1 ($k=1$) of lemma 3 will then show
that the first term in equation (31) is $O_{p}(\log^{3/2}n/\sqrt{n})$and the
second term is $O_{p}(\log^{5/2}n/\sqrt{n})$ which proves the desired result.

###  B.6 Proof of Lemma 4

The statement of this lemma follows from part (2) of Lemmata 3 and 2. We will
prove this lemma for the Poisson case first. Note that for any $\epsilon>0$
and $k\in\\{0,1\\}$, the probability
$\mathbb{P}\Big{[}\mathcal{L}_{n}^{(k)}(\bm{\theta},\bm{\delta}_{(k)}^{\sf
neb}(\hat{\lambda}))\geq\mathcal{L}_{n}^{(k)}(\bm{\theta},\bm{\delta}_{(k)}^{\sf
or})+c_{n}^{-1}\epsilon\Big{]}$ is bounded above by

$\mathbb{P}\Big{[}\mathcal{L}_{n}^{(k)}(\bm{\theta},\bm{\delta}_{(k)}^{\sf
neb}(\hat{\lambda}))-{\sf
ARE}^{(1,\mathcal{P})}_{n}(\hat{\lambda},\bm{Y})\geq\mathcal{L}_{n}^{(k)}(\bm{\theta},\bm{\delta}_{(k)}^{\sf
or})-{\sf ARE}^{(1,\mathcal{P})}_{n}(\lambda^{\sf
orc},\bm{Y})+c_{n}^{-1}\epsilon\Big{]},$

which converges to 0 by part (2) of Lemma 3. For the Binomial case, similar
arguments using part (2) of Lemma 2 suffice.

###  B.7 Proofs of Theorems 3A, 3B

We will first prove Theorem 3A. Note that $||{\bm{\delta}}_{(1)}^{\sf
neb}({\hat{\lambda}})-\bm{\delta}_{(1)}^{\pi}||_{\infty}^{2}\leq||{\bm{\delta}}_{(1)}^{\sf
neb}({\hat{\lambda}})-\bm{\delta}_{(1)}^{\pi}||_{2}^{2}$ and

$\Big{\|}{\bm{\delta}}_{(1)}^{\sf
neb}({\hat{\lambda}})-\bm{\delta}_{(1)}^{\pi}\Big{\|}_{2}^{2}=\sum_{i=1}^{n}\Big{[}\dfrac{a_{Y_{i}-1}/a_{Y_{i}}}{\hat{w}^{(1)}_{n,i}(\hat{\lambda})w^{(1)}_{p,i}}\Big{]}^{2}\Big{[}\hat{w}^{(1)}_{n,i}(\hat{\lambda})-w^{(1)}_{p,i}\Big{]}^{2}$

Now, $\hat{w}_{n,i}^{(1)}(\lambda)>0$ for every $\lambda\in\Lambda$ and
$w_{p,i}^{(1)}>0$. This fact along with assumption $A2$ and lemma A imply that
there exists a constant $c_{0}>0$ such that $||{\bm{\delta}}_{(1)}^{\sf
neb}({\hat{\lambda}})-\bm{\delta}_{(1)}^{\pi}||_{2}\leq c_{0}\log
n\|\hat{\bm{w}}_{n}^{(1)}(\hat{\lambda})-\bm{w}_{p}^{(1)}\|_{2}$. The first
result thus follows from the above inequality and Theorem 2A. To prove the
second part of the theorem, note that
$n^{-1}|\mathcal{L}_{n}^{(1)}(\bm{\theta},\bm{\delta}_{(1)}^{\pi})-\mathcal{L}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}({\hat{\lambda}}))|$ is upper bounded by

$\Big{|}\sum_{i=1}^{n}\theta_{i}^{-1}\Bigl{\\{}\Big{[}{\delta}_{(1),i}^{\sf
neb}(\hat{\lambda})\Big{]}^{2}-\Big{[}{\delta}_{(1),i}^{\pi}\Big{]}^{2}\Bigr{\\}}\Big{|}+2\Big{|}\sum_{i=1}^{n}\Big{[}{\delta}_{(1),i}^{\sf
neb}(\hat{\lambda})-{\delta}_{(1),i}^{\pi}\Big{]}\Big{|}$ (32)

Now use the fact that $\hat{w}_{n,i}^{(1)}(\lambda)>0$ for every
$\lambda\in\Lambda$ and $w_{p,i}^{(1)}>0$ to deduce, from assumption $A2$ and
Lemma A, that $|{\delta}_{(1),i}^{\sf
neb}(\hat{\lambda})+{\delta}_{(1),i}^{\pi}|\leq c_{1}\log n$ for some constant
$c_{1}>0$. Using this inequality we can now upper bound the display in
equation (32) by:

$\sum_{i=1}^{n}\Big{|}{\delta}_{(1),i}^{\sf
neb}(\hat{\lambda})-{\delta}_{(1),i}^{\pi}\Big{|}\Bigl{\\{}2+\dfrac{c_{1}\log
n}{\theta_{i}}\Bigr{\\}}\leq\Big{\|}{\bm{\delta}}_{(1)}^{\sf
neb}({\hat{\lambda}})-\bm{\delta}_{(1)}^{\pi}\Big{\|}_{\infty}\Bigl{\\{}2n+c_{1}\log
n\sum_{i=1}^{n}\theta_{i}^{-1}\Bigr{\\}}$

Thus,

$\Big{|}\mathcal{L}_{n}^{(1)}(\bm{\theta},\bm{\delta}_{(1)}^{\pi})-\mathcal{L}_{n}^{(1)}(\bm{\theta},{\bm{\delta}}_{(1)}^{\sf
neb}({\hat{\lambda}}))\Big{|}\leq\Big{\|}{\bm{\delta}}_{(1)}^{\sf
neb}({\hat{\lambda}})-\bm{\delta}_{(1)}^{\pi}\Big{\|}_{\infty}\Bigl{\\{}2+c_{1}\dfrac{\log
n}{n}\sum_{i=1}^{n}\theta_{i}^{-1}\Bigr{\\}}$

Finally, the result follows from the above display and the first part of this
theorem after noting that $\theta_{i}>0$ for all $i=1,2,\ldots$.

We will now prove Theorem 3B. Using theorem 2B, the first part of theorem 3B
follows along similar lines as the first part of theorem 3A. To prove the
second part of the theorem, note that
$n^{-1}\Big{|}\mathcal{L}_{n}^{(0)}(\bm{\theta},\bm{\delta}_{(0)}^{\pi})-\mathcal{L}_{n}^{(0)}(\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}({\hat{\lambda}}))\Big{|}$ is upper bounded by

$\Big{|}\sum_{i=1}^{n}\Bigl{\\{}\Big{[}{\delta}_{(0),i}^{\sf
neb}(\hat{\lambda})\Big{]}^{2}-\Big{[}{\delta}_{(0),i}^{\pi}\Big{]}^{2}\Bigr{\\}}\Big{|}+2\Big{|}\sum_{i=1}^{n}\theta_{i}\Big{[}{\delta}_{(0),i}^{\sf
neb}(\hat{\lambda})-{\delta}_{(0),i}^{\pi}\Big{]}\Big{|}$ (33)

and the display in equation (33) is less than or equal to

$\|{\bm{\delta}}_{(0)}^{\sf
neb}({\hat{\lambda}})-\bm{\delta}_{(0)}^{\pi}\|_{\infty}\\{\sum_{i=1}^{n}2\theta_{i}+c_{2}n\log
n\\},$

where we have used the fact that $\hat{w}_{n,i}^{(0)}(\lambda)>0$ for every
$\lambda\in\Lambda$, $w_{p,i}^{(0)}>0$ and along with assumption $A2$ and
Lemma A, $|{\delta}_{(0),i}^{\sf
neb}(\hat{\lambda})+{\delta}_{(0),i}^{\pi}|\leq c_{2}\log n$ for some constant
$c_{2}>0$. Thus for $n$ large,

$\Big{|}\mathcal{L}_{n}^{(0)}(\bm{\theta},\bm{\delta}_{(0)}^{\pi})-\mathcal{L}_{n}^{(0)}(\bm{\theta},{\bm{\delta}}_{(0)}^{\sf
neb}({\hat{\lambda}}))\Big{|}\leq c_{3}\log n\Big{\|}{\bm{\delta}}_{(0)}^{\sf
neb}({\hat{\lambda}})-\bm{\delta}_{(0)}^{\pi}\Big{\|}_{\infty}$

from which the desired result follows.

###  B.8 Proofs of Lemmata A and B

Proof of Lemma A

First note that from assumption $A4$ if $N(\delta,n)$ denotes the cardinality
of the set $\\{i:\theta_{i}\geq\epsilon^{-(1+\delta)}\log n\\}$ for some
$\delta>0$, then $N(\delta,n)\to 0$ as $n\to\infty$. We will now prove the
statement of lemma A for the case when
$Y_{i}|\theta_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\texttt{Poi}(\theta_{i})$. For distributions with bounded
support, like the Binomial model, the lemma follows trivially.

Under the Poisson model, we have
$\mathbb{P}(Y_{i}\geq\theta_{i}+t)\leq\exp\\{-0.5t^{2}/(\theta_{i}+t)\\}$ for
any $t>0$. The above inequality follows from an application of Bennett
inequality to the Poisson MGF (see Pollard (2015)). Now consider
$\mathbb{P}(\max_{i=1,\ldots,n}Y_{i}\leq\theta_{i}+t)$ and note that since
$Y_{i}$ are all independent, this probability is given by
$\prod_{i=1}^{n}[1-\exp\\{-0.5t^{2}/(\theta_{i}+t)\\}]$. Take $t=s\log n$
where $s^{2}/\\{s+\epsilon^{-(1+\delta)}\\}>4$. Then with
$\theta_{i}\leq\epsilon^{-(1+\delta)}\log n$, the above probability is bounded
below by $a_{n}=\\{1-n^{-(1+\nu)}\\}^{n}$ for some $\nu>0$. As $n\to\infty$,
$a_{n}\to 1$ which proves the statement of the lemma.

Proof of Lemma B

We begin with some remarks on the optimization problems (8) and (15). Note
that the feasible set $\bm{H}_{n}$ in equation (8) (and (15)) is compact and
independent of $\lambda$. Moreover, the optimization problem in definitions 1A
and 1B is convex. Consequently, (i) for all $\lambda\in\Lambda$, the
optimization takes place in a compact set, and (ii) the optimal solution set
corresponding to any $\lambda\in\Lambda$ is a singleton,
$\\{\hat{\bm{h}}_{n}^{(k)}(\lambda)\\}$. Now fix an $\epsilon>0$. Then for any
$\lambda\in N_{\epsilon}(\lambda_{0})\cap\Lambda$ there exists a $\delta>0$
such that the optimal solution
${\bm{h}}_{n}\coloneqq\hat{\bm{h}}_{n}^{(k)}(\lambda)\in
N_{\delta}(\hat{\bm{h}}_{n}^{(k)}(\lambda_{0}))$ and
$\hat{\mathbb{M}}_{\lambda,n}\\{{\bm{h}}_{n}\\}-\hat{\mathbb{M}}_{\lambda,n}\\{\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\\}\leq
0$. Moreover, we can re-write
$\hat{\mathbb{M}}_{\lambda_{0},n}\\{{\bm{h}}_{n}\\}-\hat{\mathbb{M}}_{\lambda_{0},n}\\{\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\\}$
as

$\displaystyle\hat{\mathbb{M}}_{\lambda_{0},n}\\{{\bm{h}}_{n}\\}-\hat{\mathbb{M}}_{\lambda,n}\\{{\bm{h}}_{n}\\}-\hat{\mathbb{M}}_{\lambda_{0},n}\\{\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\\}+\hat{\mathbb{M}}_{\lambda,n}\\{\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\\}+\hat{\mathbb{M}}_{\lambda,n}\\{{\bm{h}}_{n}\\}-\hat{\mathbb{M}}_{\lambda,n}\\{\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\\}$

The last term in the display above is negative and thus we can upper bound
$\hat{\mathbb{M}}_{\lambda_{0},n}\\{{\bm{h}}_{n}\\}-\hat{\mathbb{M}}_{\lambda_{0},n}\\{\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\\}$
by

$\hat{\mathbb{M}}_{\lambda_{0},n}\\{{\bm{h}}_{n}\\}-\hat{\mathbb{M}}_{\lambda,n}\\{{\bm{h}}_{n}\\}-\hat{\mathbb{M}}_{\lambda_{0},n}\\{\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\\}+\hat{\mathbb{M}}_{\lambda,n}\\{\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\\}$

Now apply the mean value theorem with respect to $\bm{h}_{n}$ to the function
$\hat{\mathbb{M}}_{\lambda_{0},n}\\{{\bm{h}}_{n}\\}-\hat{\mathbb{M}}_{\lambda,n}\\{{\bm{h}}_{n}\\}$
in the display above and notice that
$\hat{\mathbb{M}}_{\lambda_{0},n}\\{{\bm{h}}_{n}^{(k)}\\}-\hat{\mathbb{M}}_{\lambda_{0},n}\\{{\bm{h}}_{n}^{(k)}(\lambda_{0})\\}$
is bounded above by

$\Big{[}\nabla_{\bm{h}_{n}}\Bigl{\\{}\hat{\mathbb{M}}_{\lambda_{0},n}(\bar{\bm{h}}_{n})-\hat{\mathbb{M}}_{\lambda,n}(\bar{\bm{h}}_{n})\Bigr{\\}}\Big{]}^{T}\Big{[}\bm{h}_{n}-{\bm{h}}_{n}^{(k)}(\lambda_{0})\Big{]}$

where
$\bar{\bm{h}}_{n}=\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})+\tau\\{\bm{h}_{n}-\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\\}$
for some $\tau\in(0,1)$ and
$\nabla_{\bm{h}}\hat{\mathbb{M}}_{\lambda,n}({\bm{h}})$ is the partial
derivative of $\hat{\mathbb{M}}_{\lambda,n}({\bm{h}})$ with respect to
$\bm{h}$. Using
$\nabla_{\bm{h}_{n}}[\hat{\mathbb{M}}_{\lambda_{0},n}({\bm{h}}_{n})-\hat{\mathbb{M}}_{\lambda,n}({\bm{h}}_{n})]=\nabla^{2}_{\bm{h}_{n},\lambda}\hat{\mathbb{M}}_{\lambda_{0},n}({\bm{h}}_{n})(\lambda-\lambda_{0})+o(|\lambda-\lambda_{0}|)$
we get

$\hat{\mathbb{M}}_{\lambda_{0},n}\\{{\bm{h}}_{n}\\}-\hat{\mathbb{M}}_{\lambda_{0},n}\\{\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\\}\leq\sup_{{\bm{h}}\in
N_{\delta}(\hat{\bm{h}}_{n}^{(k)}(\lambda_{0}))}\Big{[}\Big{\|}\nabla^{2}_{\bm{h}_{n},\lambda}\hat{\mathbb{M}}_{\lambda_{0},n}({\bm{h}})+o(1)\Big{\|}_{2}\Big{]}\Big{|}\lambda-\lambda_{0}\Big{|}\Big{\|}{\bm{h}}_{n}-\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\Big{\|}_{2}$

Moreover assumption $A3$ implies that

$\hat{\mathbb{M}}_{\lambda_{0},n}\\{{\bm{h}}_{n}^{(k)}\\}-\hat{\mathbb{M}}_{\lambda_{0},n}\\{\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\\}\geq
c\Big{\|}{\bm{h}}_{n}^{(k)}-\hat{\bm{h}}_{n}^{(k)}(\lambda_{0})\Big{\|}_{2}^{2}$

The desired result thus follows from the above two displays with

$L=\sup_{{\bm{h}}\in
N_{\delta}(\hat{\bm{h}}_{n}^{(k)}(\lambda_{0}))}\|\nabla^{2}_{\bm{h}_{n},\lambda}\hat{\mathbb{M}}_{\lambda_{0},n}({\bm{h}})+o(1)\|_{2}/c.$

## References

  * Aizer and Doyle Jr (2015) Anna Aizer and Joseph J Doyle Jr.  Juvenile incarceration, human capital, and future crime: Evidence from randomly assigned judges.  _The Quarterly Journal of Economics_ , 130(2):759–803, 2015.
  * Bandari et al. (2012) Roja Bandari, Sitaram Asur, and Bernardo A Huberman.  The pulse of news in social media: Forecasting popularity.  _ICWSM_ , 12:26–33, 2012.
  * Bonnans and Shapiro (2013) J Frédéric Bonnans and Alexander Shapiro.  _Perturbation analysis of optimization problems_.  Springer Science & Business Media, 2013.
  * Brown (2008) Lawrence D Brown.  In-season prediction of batting averages: A field test of empirical bayes and bayes methodologies.  _The Annals of Applied Statistics_ , pages 113–152, 2008.
  * Brown and Greenshtein (2009) Lawrence D Brown and Eitan Greenshtein.  Nonparametric empirical bayes and compound decision approaches to estimation of a high-dimensional vector of normal means.  _The Annals of Statistics_ , pages 1685–1704, 2009.
  * 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.
  * Chwialkowski et al. (2016) Kacper Chwialkowski, Heiko Strathmann, and Arthur Gretton.  A kernel test of goodness of fit.  JMLR: Workshop and Conference Proceedings, 2016.
  * Clevenson and Zidek (1975) M Lawrence Clevenson and James V Zidek.  Simultaneous estimation of the means of independent poisson laws.  _Journal of the American Statistical Association_ , 70(351a):698–705, 1975.
  * Damm and Dustmann (2014) Anna Piil Damm and Christian Dustmann.  Does growing up in a high crime neighborhood affect youth criminal behavior?  _American Economic Review_ , 104(6):1806–32, 2014\.
  * Efron (2011) Bradley Efron.  Tweedie’s formula and selection bias.  _Journal of the American Statistical Association_ , 106(496):1602–1614, 2011.
  * Efron (2012) Bradley Efron.  _Large-scale inference: empirical Bayes methods for estimation, testing, and prediction_ , volume 1.  Cambridge University Press, 2012.
  * Fourdrinier and Robert (1995) Dominique Fourdrinier and Christian P Robert.  Intrinsic losses for empirical bayes estimation: A note on normal and poisson cases.  _Statistics & probability letters_, 23(1):35–44, 1995.
  * Fourdrinier et al. (2018) Dominique Fourdrinier, William E Strawderman, and Martin T Wells.  _Shrinkage Estimation_.  Springer, 2018.
  * Fu et al. (2017) Anqi Fu, Balasubramanian Narasimhan, and Stephen Boyd.  Cvxr: An r package for disciplined convex optimization.  _arXiv preprint arXiv:1711.07582_ , 2017.
  * Fu et al. (2018) Luella Fu, Gareth James, and Wenguang Sun.  Nonparametric empirical bayes estimation on heterogeneous data.  2018\.
  * James and Stein (1961) William James and Charles Stein.  Estimation with quadratic loss.  In _Proceedings of the fourth Berkeley symposium on mathematical statistics and probability_ , volume 1, pages 361–379, 1961.
  * Jiang et al. (2009) Wenhua Jiang, Cun-Hui Zhang, et al.  General maximum likelihood empirical bayes estimation of normal means.  _The Annals of Statistics_ , 37(4):1647–1684, 2009.
  * Klenke (2014) Achim Klenke.  _Probability Theory: A Comprehensive Course_ , pages 331–349.  Springer London, 2014.
  * Koenker and Gu (2017) Roger Koenker and Jiaying Gu.  Rebayes: An r package for empirical bayes mixture methods.  _Journal of Statistical Software_ , 82(1):1–26, 2017.
  * Koenker and Mizera (2014) Roger Koenker and Ivan Mizera.  Convex optimization, shape constraints, compound decisions, and empirical bayes rules.  _Journal of the American Statistical Association_ , 109(506):674–685, 2014.
  * Koski et al. (2018) Susan V Koski, David Bowers, and SE Costanza.  State and institutional correlates of reported victimization and consensual sexual activity in juvenile correctional facilities.  _Child and Adolescent Social Work Journal_ , 35(3):243–255, 2018.
  * Liu and Wang (2016) Qiang Liu and Dilin Wang.  Stein variational gradient descent: A general purpose bayesian inference algorithm.  In _Advances In Neural Information Processing Systems_ , pages 2378–2386, 2016.
  * Liu et al. (2016) Qiang Liu, Jason D Lee, and Michael I Jordan.  A kernelized stein discrepancy for goodness-of-fit tests.  In _Proceedings of the International Conference on Machine Learning (ICML)_ , 2016.
  * Moniz and Torgo (2018) Nuno Moniz and Luís Torgo.  Multi-source social feedback of online news feeds.  _arXiv preprint arXiv:1801.07055_ , 2018.
  * Noack (1950) Albert Noack.  A class of random variables with discrete distributions.  _The Annals of Mathematical Statistics_ , 21(1):127–132, 1950.
  * Oates et al. (2017) Chris J Oates, Mark Girolami, and Nicolas Chopin.  Control functionals for monte carlo integration.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 79(3):695–718, 2017.
  * Pollard (2015) David Pollard.  A few good inequalities.  _available at:http://www.stat.yale.edu/~pollard/Books/Mini/Basic.pdf_, 2015.
  * Robbins (1956) Herbert Robbins.  An empirical bayes approach to statistics.  Technical report, COLUMBIA UNIVERSITY New York City United States, 1956\.
  * Sato and Ken-Iti (1999) Ken-iti Sato and Sato Ken-Iti.  _Lévy processes and infinitely divisible distributions_.  Cambridge university press, 1999.
  * Serfling (2009) Robert J Serfling.  _Approximation theorems of mathematical statistics_ , volume 162.  John Wiley & Sons, 2009.
  * Shmueli et al. (2005) Galit Shmueli, Thomas P Minka, Joseph B Kadane, Sharad Borle, and Peter Boatwright.  A useful distribution for fitting discrete data: revival of the conway–maxwell–poisson distribution.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 54(1):127–142, 2005.