The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
36,264 characters
Robustifying Empirical Bayes
\title{Robustifying Empirical Bayes}
\author{Roger Koenker and Jiaying Gu}
\thanks{Version: \today . The authors wish to express their appreciation to Pat Kline, Toru
Kitagawa and Peter Bickel for comments on a previous draft. Details on all the computations for the
figures and tables are available from \url{https://rkoenker.github.io/www/roger/research/ebayes/ebayes.html}.}
\begin{abstract}
Two strategies are explored for robustifying classical denoising procedures for the
Gaussian sequence model. First, the Hodges and Lehmann (1952)
restricted Bayes approach is used to reduce sensitivity to the specification
of the initial prior distribution. Second, alternatives to the Gaussian
noise assumption are explored. In both cases proposals of Huber (1964)
and Mallows (1978) play a crucial role.
\end{abstract}
\maketitle
\section{Introduction}
The Gaussian sequence model can be viewed as a compound decision problem with
observed $X_i \sim \mathcal{N} (\theta_i, 1), \; i = 1, \dots , n$. The objective is to estimate the
$\theta \in \mathbb{R}^n$ subject to quadratic loss. We will denote the standard Gaussian
density and cumulative by $\varphi$ and $\Phi$ respectively. Observations are assumed to be
exchangeable, so their marginal density is given by,
\[
f_G (x) = \int \varphi (x| \theta) dG(\theta),
\]
for some mixing distribution $G$.
Were $G$ known the optimal (Bayes) decision rule is given by Tweedie's formula, \cite{efron.11}
\[
\hat \theta_i = \delta^B (x_i) = x_i + f_G'(x_i)/f_G (x_i).
\]
When $G$ is unknown various shrinkage procedures have been proposed, initiated by the
fundamental papers of \cite{stein56} and \cite{robbins56}. The extensive literature on
Stein shrinkage has offered a rich assortment of practical frequentist and Bayesian procedures
for improving upon the naive maximum likelihood estimator, $\delta (x_i) = x_i$ in terms
of quadratic loss. Among these procedures more recently the
nonparametric maximum likelihood estimator (NPMLE) of \cite{kw},
\[
\hat G = \mbox{argmax}_{G \in \mathcal{G}} \big \{ \sum_{i=1}^n \log (f_G (x_i) \big \},
\]
has been proposed as a plug-in estimator for $G$, \cite{jz}, \cite{km} and \cite{soloff}.
This $G$-modeling strategy -- in the terminology of \cite{e19} -- performs well in simulations,
e.g. \cite{km}, \cite{gk16}, \cite{kg26}, relative to alternatives that attempt to
estimate $f_G$ directly or
that make a priori assumptions about the form of $G$. However, it is obviously subject
to the criticism that the Gaussian assumption on the likelihood is quite strong,
and priors are never terribly convincing.
In what follows we consider two basic strategies for robustifying empirical Bayes
procedures. The first, following a proposal of \cite{hodgeslehmann}, seeks protection
from excessive confidence in our initial prior on $G$ by bounding pointwise risk
thereby offering a compromise between minimax and Bayes decision rules. The second,
acknowledges scepticism about the strictly Gaussian form of $\varphi$, the distribution
of the model noise.
We find that in accordance with familiar robustness lore that modest modifications of
of our initial prior or the Gaussian noise assumption can significantly improve performance
of empirical Bayes decision rules while sacrificing only modest performance in the
event that the initial prior or the Gaussian noise assumptions are valid.
\section{Bayes Risk and Brown's Identity}
Bayes risk in the Gaussian sequence model is,
\[
r(G, \delta) = \int R(\delta, \theta) dG(\theta),
\]
with
\[
R(\delta, \theta) =\mathbb{E}_\theta[(\delta(X) - \theta)^2] =
\int ( \delta(x)-\theta)^2 \varphi(x - \theta) dx.
\]
Plugging the optimal Bayes rule back into $r(G, \delta)$, we obtain
Brown's identity, \cite{Brown71}:
\begin{align*}
r(G, \delta^{B}) & =
\int (x-\theta + \frac{f_G'(x)}{f_G(x)})^2 \varphi(x - \theta) dx dG(\theta)\\
& = 1 + 2 \mathbb{E}\Big [\Big (\frac{f_G'(X)}{f_G(X)}\Big )'\Big ] +
\mathbb{E}\Big [\Big (\frac{f_G'(X)}{f_G(X)}\Big )^2 \Big ]\\
& = 1 +2 \mathbb{E}\Big [ \frac{f_G^{''}(X)}{f_G(x)} - \Big( \frac{f_G'(x)}{f_G(x)}\Big)^2 \Big ]
+ \mathbb{E}\Big [\Big (\frac{f_G'(X)}{f_G(X)}\Big )^2 \Big ]\\
& = 1- \mathbb{E}\Big[ \Big (\frac{f_G'(X)}{f_G(X)}\Big )^2 \Big ]\\
& = 1 - I(\Phi * G).
\end{align*}
The second equality follows from Stein's lemma, the fourth from the fact that
$\int f_G^{''} (x) dx = 0$,
and the last equality from the definition of Fisher
information for distributions with absolutely continuous densities.
It may seem curious that Bayes risk of the optimal empirical Bayes rule for the Gaussian
sequence model reduces to Fisher information for a scalar location parameter of the convolution
distribution $\Phi * G$. We will exploit the latter connection in the next section to consider least favorable
alternative priors for an initial prior in which we lack complete confidence. Such modified priors
offer some compromise between strictly Bayesian and minimax procedures.
Choosing contamination models by minimizing Fisher information
is a classical strategy for choosing alternatives for the univariate Gaussian location
and regression problems following \cite{huber64}. The convolution form of the Fisher
information for Bayes risk leads us back to an alternative proposal of \cite{mallows78}
as well.
A natural objection to many empirical Bayes procedures is that they place unjustified reliance on an initial
prior. While such procedures may still perform well with respect to ensemble risk, they may also fail spectacularly
for some subpopulations or individuals. This concern underlies the limited translation
proposal of \cite{efronmorris}. We will see that bounding minimax risk by modifying an initial
prior can often serve to soften the impact of these failings.
\section{Restricted Bayes Solutions}
In an effort to balance minimax and Bayes solutions, \cite{hodgeslehmann} proposed solving,
\[
\min_\delta r(G_0,\delta) \; s.t. \; \max_\theta R(\delta, \theta) \leq 1 + t,
\]
for an initial prior $G_0$ and some $t > 0$. Since $\max_\theta R(X,\theta)=1$ corresponding to
the worst pointwise risk among all decision rules, achieved by the MLE estimator $\delta(X) = X$, the
Hodges and Lehmann modified decision rule is thereby constrained to do uniformly
well over the entire parameter space with risk bounded by $1+t$, while minimizing
the Bayes risk under prior $G_0$, hence called the restricted Bayes rule. They show that this is equivalent
to solving,
\[
\min_\delta \; \max_{G \in \mathcal{G}_\epsilon(G_0)} r(G, \delta),
\]
with $\mathcal{G}_\epsilon(G_0) = \{ G = (1 - \epsilon) G_0 + \epsilon H \}$ for some $\epsilon \in (0,1)$ depending
upon $t$, where $H$ is an arbitrary distribution for $\theta$. Due to the minimax theorem we can switch the
minimization and the maximization, and the resulting optimal rule $\delta^*$ is the posterior mean
of $\theta$ under the least favorable prior from the class $\mathcal{G}_{\epsilon}(G_0)$. \cite{berger85} comments that
``it is \emph{very} difficult to determine such $\delta^*$; furthermore, this 'optimal' $\delta^*$ is
usually extremely messy and difficult to work with.'' On the contrary, with the aid of
modern convex optimization techniques we find them quite tractable and elegant.
\cite{bickel83} considers the case with $G_0$ having point mass one at zero. Then, by the
Brown identity, the Hodges and Lehmann problem is equivalent to solving,
\[
\max_{G \in \mathcal{G}_\epsilon(\delta_0)} (1 - I (\Phi * G)) = \min_{G \in \mathcal{G}_\epsilon(\delta_0)} I (\Phi * G)
\]
This is the problem posed by \cite{mallows78} motivated by robustness considerations for
time-series problems with additive outliers.
Mallows conjectured that the least favorable $G$ would be discrete,
supported on the integers with mass declining exponentially.
\cite{bickel83} reports a modified conjecture of Donoho that relaxes the spacing
of the Mallows mass points, but is otherwise similar. Neither conjecture seems to be strictly
correct, but numerical computations confirm the nearly exponential decay of the mass.
\cite{bickel1983minimizing} provide a detailed discussion of the discrete nature of the Mallows
solutions based on the analyticity of the objective function. See also \cite{johnstone94}.
\subsection{Computing Mallows's Least Favorable Distribution} \label{sec: computation}
As noted by \cite{bickel1983minimizing} and \cite{marazzi} the \cite{mallows78} problem is convex.
\cite{marazzi} suggests a gridding strategy that imposes an exponentially declining mass condition.
This produces a remarkably accurate solution for an initial
Gaussian prior employing generic optimization software. Using modern convex optimization software,
we show that the problem can be efficiently solved numerically for any prior distribution.
Our implementation embodied in the function \texttt{HodgesLehmann} in the our REBayes package
for the R language employs the Mosek \cite{mosek} optimizer and provides a general interface for
computing either the Huber or Mallows solutions for an arbitrary initial prior, $G_0$.
We now describe our procedure for the simplest (Dirac) initial prior, $G_0 = \delta_0$.
Our objective is to solve
\[
\underset{f \in \mathcal{K} }{\min} \int \frac{f'(x)^2}{f(x)} dx
\]
with
\[
\mathcal{K}_\epsilon =\Big \{f = (1 - \epsilon) \int \varphi(x- \theta)d \delta_0(\theta) +
\epsilon \int \varphi(x- \theta) dH(\theta)\Big \}
\]
This can be solved, on a grid of $x$, $\{x_1< x_2 < \dots <x_M\}$
and a grid of $\theta$ as $\{ \theta_1 < \theta_2 < \dots < \theta_L\}$
as the rotated quadratic cone convex optimization problem:
\[
\min \sum_{i = 1}^M w_i
\]
subject to
\begin{align*}
u_i &= f_{i+1}-f_i, \\
v_i &= \frac{1}{2}( f_{i+1}+f_i),\\
u_i^2 &\leq 2v_i w_i\\
f_i &= (1 - \epsilon) \varphi(x_i) + \epsilon \sum_{j=1}^L \varphi(x_i - \theta_j) h_j\\
h & \in \mathcal{S} \equiv \{h \in \mathbb{R}^L : \sum_{j=1}^L h_j = 1, \; h_j \geq 0, \; j = 1, \dots , L \}
\end{align*}
Provided that the grids are sufficiently finely spaced interior point optimization in Mosek
is capable of producing very accurate solutions very efficiently.
In Figure \ref{fig.Mallows} we illustrate a Mallows marginal density, $f^M$, its corresponding
mixing distribution, $P^M$, and plot the log mass of the discrete mass points of $P^M$ at their
respective locations. At first glance, it seems that the mass points are approximately equally spaced and
have mass that declines exponentially. However, on closer examination the spacing of
the mass points in the right tail are estimated to be: $\{ 1.96, 1.80, 1.70, 1.61, 1.52,
1.39, 1.29, 1.37, 1.55, 1.70, 1.91\}$, which seems sufficiently non-uniform to call the
uniform spacing conjecture into question. \cite{djm} suggest an alternative computational
strategy for the Mallows problem using a parametric model that assumes equal spacing of the
mass points. The grid for evaluation of $f$ is equally spaced
from -30 to 30 with 500 points of evaluation. The grid for evaluation of $P^M$ is also equally
spaced from -20 to 20 with 4003 points of evaluation. For purposes of illustration the mass
at $\theta = 0$ is taken to be 0.2.
\begin{figure}
\begin{center}
\resizebox{ \textwidth}{!}{{\includegraphics{figs/Mallows.pdf}}}
\end{center}
\caption{Mallows least favorable marginal density, probability mass function of the Mallows
mixing distribution and log mass of the
mixing distribution as a function of location indicating the approximate exponentiality
of the mixing distribution as conjectured by Mallows.}
\label{fig.Mallows}
\end{figure}
\begin{figure}
\begin{center}
\resizebox{ .6\textwidth}{!}{{\includegraphics{figs/HLmass.pdf}}}
\end{center}
\caption{Mass points of the Mallows contamination distribution $H^*$ with
$G_0 = 0.5 \delta_{-2} + 0.5 \delta_2$ for $\epsilon = 0.2$. The solid black
curve depicts the pointwise risk function of the decision rule $\delta^*$, constructed using Mallows's least favorable prior $G^*$. The vertical green lines depict the location of the mass points of
the solution $H^*$, while their length represents the amount of mass assigned to each. In
accordance with the Tukey ``hanging rootogram'' principle these lengths are rescaled as the square
root of the respective masses. The dashed horizontal line represents the bound on
the pointwise risk imposed by the Hodges-Lehmann constraint in this case approximately
1.67 induced by the choice of $\epsilon = 0.2$. The constraint is binding at mass points of $H^*$.}
\label{fig.HLmass}
\end{figure}
In the previous example we have taken the initial prior, $G_0$ as Dirac, but
there is no obstacle to starting from any other initial prior. To provide some additional
intuition about the nature of the Mallows solution for general $G_0$, we illustrate in Figure \ref{fig.HLmass}
a plot of the pointwise risk function $R(\delta^* , \theta) := \mathbb{E}_\theta[(\delta^*(X)-\theta)^2]$ of the decision rule $\delta^*(\cdot)$, constructed as the posterior mean of $\theta$ using the least favorable Mallows prior
\[
G^* (\theta) = (1 - \epsilon ) G_0 (\theta) + \epsilon H^* (\theta)
\]
where $G_0 (\theta)$ is taken to be an equally weighted mixture of two point masses at -2 and 2.
The tangencies in this plot with the horizontal dotted line marking out $\sup_\theta R(\delta^*, \theta) = 1 + t$ coincide
with the location of the mass points of the solution of $H^*$ indicated in the plot by the
vertical green lines. The bound, $1 + t$, is the pointwise
risk bound chosen to constrain the Hodges-Lehmann restricted Bayes rule, which can be constructed using Mallows's least favorable prior $G^*$.
Since the initial prior $G_0$ places all its mass on the two points $\{-2,2\}$ the Mallows modification
hedges this bet by placing a considerable mass at zero and exponentially declining mass at a few points
below -2 and above +2. This figure is strongly reminiscent of Figure 5.5 of \cite{lindsay}
illustrating the location of mass points of the NPMLE.
\subsection{Some Examples}
We now consider several special cases of the Hodges and Lehmann restricted Bayes approach.
In each case we consider not only the Mallows equivalent form of the Hodges-Lehmann modification, but
also a Huber alternative that relaxes the Mallows objective of minimizing the Fisher information over convolutions by
minimizing over the entire class of contamination distributions for $X$. Taking the least favorable density and plug into the Tweedie formula gives rise the Huber procedure to estimate $\theta$ for each value of $x$.
The Mallows rule has the obvious advantage that it yields a Bayes decision rule while
the corresponding Huber procedure does not. This is particularly evident in the third
example of Casella and Strawderman where the unrestricted Bayes rule is minimax, so the Hodges-Lehman's
restriction on point-wise risk is unbinding and the restricted Bayes rule coincides with the unrestricted,
but the Huber procedure is inadmissible. On the other hand there is something attractive about the Huber
rules that it can be shown that they are necessarily monotone. \cite{donohoreeves} propose
an alternative based on the \cite{huber74} spline that minimizes Fisher information over a Kolmogorov
neighborhood of the marginal density of $X$ specified by a finite number of evaluations of its quantile function.
They then apply the Tweedie formula with the resulting least favorable density.
An implementation of this procedure is also included in the REBayes package with the function
\texttt{HuberSpline}, although we do not pursue it further in this paper.
\vspace{5mm}
\begin{description}
\item[Dirac $G_0$] This is the case considered by \cite{bickel83} and anticipated by \cite{mallows78}.
Bickel first simplifies the problem of minimizing $I(\Phi * G)$ by considering the relaxed
\cite{huber64} problem of minimizing $I(F)$ over $\mathcal{F} = \{ F = (1 - \epsilon) \Phi + \epsilon W\}$ where $W$ is an arbitrary distribution for $X$.
The least favorable distribution $F$ has the well known score function,
\[
-f^\prime (x)/f(x) =
\begin{cases} x & |x| \leq k\\
k \; \text{sign}(x) & |x| > k.
\end{cases}
\]
Applying Tweedie's formula gives rise to the hard thresholding Huber rule,
\[
\delta(x) = x + f^\prime (x)/f(x) =
\begin{cases} 0 & |x| \leq k\\
x - k \; \text{sign}(x) & |x| > k
\end{cases}
\]
We contrast this with the numerical solution of the corresponding
Mallows problem in Figure \ref{fig.dirac}. The Mallows rule offers a soft thresholding
alternative to the piecewise linear Huber rule that oscillates around the Huber rule
in the tails.
\begin{figure}
\begin{center}
\resizebox{ 0.6\textwidth}{!}{{\includegraphics{figs/dirac.pdf}}}
\end{center}
\caption{The figure contrasts the restricted Hodges-Lehmann decision rules based on an initial Dirac
prior with mass one at zero: the piecewise linear Huber rule imposes hard thresholding near zero
while the Mallows rule allows soft thresholding near zero and oscillates around the Huber rule in
the tails. Here $\epsilon = 0.4$.}
\label{fig.dirac}
\end{figure}
\vspace{5mm}
\item[Gaussian $G_0$] This is the case considered by \cite{efronmorris}. The initial prior is
Gaussian, $G_0 \sim \mathcal{N}(0, A)$, so the score function of the least favorable Huber density, which minimizes Fisher information of distributions in the contamination class $\mathcal{F}_\epsilon = \{ F = (1-\epsilon) N(0,A+1) + \epsilon W\}$, is
\[
-f^\prime (x)/f(x) =
\begin{cases} \frac{1}{A+1}x & |x| \leq k(A+1)\\
k \; \text{sign}(x) & |x| > k(A+1)
\end{cases}
\]
and Tweedie's formula yields the piecewise linear ``limited translation'' rule,
\[
\delta(x) = x + f^\prime (x)/f(x) =
\begin{cases} \frac{A}{A+1} x & |x| \leq k(A+1)\\
x - k \; \text{sign} (x) & |x| > k(A+1),
\end{cases}
\]
proposed by Efron and Morris.
Figure \ref{fig.gauss} contrasts the Huber and Mallows forms of the Hodges-Lehmann restricted Bayes rule
with the unrestricted Bayes rule. We stick to linear shrinkage for values of $x$ in the middle but refrain
from shrinking at the two tails. Again, the Mallows rule smooths the Huber rule in the center
and oscillates around the Huber rule in the tails. Remarkably, \cite{efronmorris} in
their Appendix Figure B' already illustrates this behavior of the restricted Mallows rule.
See also Figure 1 in \cite{marazzi}. Extension of this example to James-Stein forms for
$G_0$ is straightforward although computation of the associated risk is not as simple as shown by
\cite{efronmorrisii}.
\begin{figure}
\begin{center}
\resizebox{ 0.6\textwidth}{!}{{\includegraphics{figs/gauss.pdf}}}
\end{center}
\caption{The figure contrasts the Huber and Mallows forms of the restricted
Hodges-Lehmann decision rules based on an initial Gaussian
prior with variance one: the piecewise linear Huber rule is linear in the center
while the Mallows rule is smoother near zero and oscillates around the Huber rule in
the tails. Again, $\epsilon = 0.4$.}
\label{fig.gauss}
\end{figure}
\vspace{5mm}
\item[Two Point $G_0$ I] This is the case considered by \cite{casella1981estimating}, the
initial prior is $G_0 = 0.5 \delta_{-1} + 0.5 \delta_{1}$ and we restrict the domain
of the parameter $\theta$ to the interval $[-1,1]$. Casella and Strawderman show
that the Bayes rule,
\[
\delta^B (x) = \tanh (x)
\]
attains maximal pointwise risk of
\[
R(\delta^B, 1) = \mathbb{E}_{X \sim \mathcal{N}(1,1)} [ \tanh(X) - \theta)^2] \approx 0.45 < 1,
\]
at $\theta = \pm 1$. Consequently, $\delta^B$ is minimax and $G_0$ is least favorable for
all $t(\epsilon) \geq 0$ as asserted in Theorem 3.1 of \cite{casella1981estimating}.
In contrast, the Huber modification of the Bayes rule,
\[
\delta^H(x) =
\begin{cases}
x + k & x < -b\\
\delta^B(x) & |x| < b\\
x -k & x > b
\end{cases}
\]
with $k = -(\log f_0)'(b)$ and $b$ solves $\int_{-b}^b f_0(x) dx + \frac{2 f_0(b)}{k} = (1-\epsilon)^{-1}$, has strictly greater risk for all $\theta \in \Theta = [-1,1]$.
To see this, let
$f_0(x) = \frac{1}{2} \varphi(x-1) + \frac{1}{2} \varphi(x+1)$,
$g^H(x) = \delta^H(x) - x$ and $g^B(x) = \delta^B(x) - x = f_0'(x)/f_0(x)$.
By Stein's lemma,
\begin{align*}
R(\delta^H ,\theta) & - R( \delta^B , \theta)\\
& = 2 \mathbb{E}_\theta[ (g^H(X))' - (g^B(X))'] + \mathbb{E}_\theta [ (g^H(X))^2 - (g^B(X))^2]\\
& = -2 \mathbb{E}_\theta[ 1\{|X|>b\} (\log f_0(X))'' ]\\
& \quad + \mathbb{E}_\theta [ 1\{|X|>b\} (k^2 - ((\log f_0(X))')^2]
\end{align*}
The first term is positive because $f_0(x)$ is log-concave.
The second term is also positive because $k = -(\log f_0)'(b)$ and
$k^2 >( (\log f_0(x))')^2$ for $|x|>b$, hence we conclude $\delta^H$ is inadmissible.
Figure \ref{fig.2ptrisk} illustrates the pointwise risk functions of the Bayes, Mallows
and Huber decision rules on $[-1,1]$. Since the Bayes rule and its Mallows modification
are identical the two are indistinguishable in the figure.
\begin{figure}
\begin{center}
\resizebox{.6\textwidth}{!}{{\includegraphics{figs/2ptrisk.pdf}}}
\caption{Pointwise risk for the Bayes rule and its Huber and Mallows
restricted modifications. The Huber risk function assumes $\epsilon = 0.4$.
The Bayes rule and its Mallows modification
are indistinguishable so the Mallows risk has been artificially increased
by 0.01 to make them both almost distinguishable.} \label{fig.2ptrisk}
\end{center}
\end{figure}
\vspace{5mm}
\item[Two point $G_0$ II] When the two points of support of $G_0$ are more widely separated, for example,
$G_0 = \frac{1}{2} \delta_2 + \frac{1}{2} \delta_{-2}$, and the support of $\theta$ is the whole
real line, the marginal $ f_0(x) = \int \varphi(x-\theta) dG_0(\theta)$ is bimodal, the Huber
least favaroable density for $X$ is described in the following proposition whose proof appears in
Appendix \ref{app.A}.
\begin{proposition} \label{prop:Huber}
For $\epsilon$ sufficiently large,\footnote{For very small $\epsilon$ no modification
in the center of the distribution is required; only the tail behavior is modified.
In our example this threshold is about $\epsilon_0 < 0.0001$.}
the solution to,
\[
\min_{F \in \mathcal{F}_\epsilon} I(F)
\]
with $\mathcal{F}_\epsilon = \{F: F= (1-\epsilon) F_0 + \epsilon W\}$ where $F_0(x) = \frac{1}{2} \Phi(x+2) + \frac{1}{2} \Phi(x-2)$ and $W$ is arbitrary distribution of $X$. The least favorable distribution has its density of
the form:
\[
f^*(x) = \begin{cases}
(1-\epsilon) f_0(x) e^{k(x+b)} & x \leq -b\\
(1-\epsilon) f_0(x) & |x| \in [c,b]\\
A^2 \cosh^2(kx/2) & x \in [-c,c]\\
(1-\epsilon) f_0(x) e^{-k(x-b)} & x \geq b
\end{cases}
\]
where for a given $k$, and $(b,c,A)$,
\begin{align*}
k &= -(\log f_0)'(b)\\
k \cdot \tanh(kc/2) & = (\log f_0)'(c)\\
A^2 \cosh^2(kc/2) &= (1-\epsilon) f_0(c).
\end{align*}
The constant $k$ is defined implicitly by,
\[
\int_{-c}^c A^2 \cosh^2(kx/2) dx + 2 \int_c^b (1-\epsilon) f_0(x)dx + \frac{2(1-\epsilon) f_0(b)}{k} = 1.
\]
The associated Huber decision rule is then,
\[
\delta^*(x) =
\begin{cases}
x + k \cdot \tanh(kx/2) & x\in [-c,c]\\
x + f_0(x)'/f_0(x) & |x| \in [c,b]\\
x-k & x \geq b\\
x+k & x \leq -b
\end{cases}
\]
\end{proposition}
In Figure \ref{fig.2pt2} we compare this shrinkage rule with the initial Bayes rule
and the Hodges-Lehmann rule (labeled as Mallows) in the left panel. Again we
see that the Mallows's rule oscillates around the Huber rule in the tails, while being
somewhat smoother in the around zero. In the right panel of Figure \ref{fig.2pt2} we
depict the pointwise risk of the three procedures for the choice, $\epsilon = 0.2$.
Worst case risk for the Huber and Mallows rules is about 1.67 while the worst case
risk of the Bayes rule is unbounded. The Huber rule is now
clearly admissible, but still not as attractive as the Mallows rule, in the sense that
its Bayes risk is strictly larger than the Mallows rule.
\begin{figure}[h!]
\centering
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{figs/rule_twopointG}
\end{subfigure}
\begin{subfigure}[b]{0.5\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/risk_twopointG}
\end{subfigure}
\caption{Decision rules and associated pointwise risk functions for three procedures.
The initial prior is: $G_0 = \frac{1}{2} \delta_2 + \frac{1}{2} \delta_{-2}$,
with contamination level $\epsilon = 0.2$. The worst case risk for
Mallows and Huber rules is about 1.67. }\label{fig.2pt2}
\end{figure}
\end{description}
\subsection{Restricted Empirical Bayes rules}
Having examined several examples of the Hodges and Lehmann restricted Bayes rules for some simple initial
prior distributions, we now consider an empirical Bayes counterpart with $G_0$ estimated by maximum
likelihood as proposed by \cite{kw} and anticipated by \cite{r50}.
Given a sample from the compound decision problem posed in the introduction, we consider
the nonparametric maximum likelihood estimator $\hat G$ as the initial $G_0$ and then
proceed to construct a modified prior according to the principles laid out by Hodges and
Lehmann.
\begin{figure}
\begin{center}
\resizebox{\textwidth}{!}{{\includegraphics{figs/npmle.pdf}}}
\end{center}
\caption{The left panel of the figure contrasts the NPMLE prior with the modified
Mallows prior. The right panel contrasts the Huber and Mallows forms of the restricted
Hodges-Lehmann decision rules based on the initial Kiefer-Wolfowitz NPMLE
prior: the Huber rule is almost linear in the center
while the Mallows rule is smoother near zero and oscillates around the Huber rule in
the tails. In this example we take $\epsilon = 0.1$ on the presumption that the initial
empirical prior is more reliable than in the prior examples.}
\label{fig.npmle}
\end{figure}
We illustrate the consequences of this in Figure \ref{fig.npmle}. Data is generated from the
standard Gaussian sequence model with $G_0 \sim U[0,3]$. Heavy black vertical lines indicate the
original NPMLE $\hat G$ while the red vertical lines indicate the mass points of the modified
prior. While some alteration of the mass in the center of the estimated mixing distribution can
be seen, the main change is the new mass points in the tails which decline exponentially in
accordance with the Mallows's conjecture. This feature resembles the \cite{efronmorris} limited
translation estimator that imposed linear shrinkage in the center of the distribution, but eschewed
shrinkage in the tails. In baseball terms: a few extremely good hitters deserve their exalted averages.
Note that the restricted and unrestricted prior decision rules illustrated in the right panel of the
figure agree quite closely on the support of the true $\theta$'s, but diverge sharply beyond this support.
In the following theorem, we establish that the excess Bayes risk, the difference between an
oracle Mallows rule $\delta^M$ with known $G_0$ and an empirical Mallows rule $\hat \delta^M$ with
estimated (NPMLE) $\hat G$ vanishes asymptotically. The proof makes use of machinery from variational
analysis and the important feature that the score function of the oracle and EB Mallows
least favorable density is uniformly bounded.
The proof appears in Appendix \ref{app.B}.
\begin{thm}\label{thm: EBMallows}
Provided the NPMLE estimator of the marginal density $f_{\hat G}$ is Hellinger consistent
for the true marginal density $f_{G_0}$, then as $n \to \infty$,
\[
r(G_0, \hat \delta^M) - r(G_0, \delta^M) \to 0
\]
\end{thm}
\subsection{Some simulation experience} \label{sec: simu}
To evaluate the cost of imposing restrictions on the prior of the Hodges-Lehmann type we consider
three examples in this section:
\begin{description}
\item[Gaussian $G_0$] $G_0 \sim \mathcal{N} (0,1)$.
\item[Uniform $G_0$] $G_0 \sim U[-2,2]$.
\item[Twopoint $G_0$] $G_0 \sim 0.5(\delta_{-2} + \delta_2)$.
\end{description}
For each of these settings we compute mean squared error (MSE) for each of the following
decision rules:
\begin{description}
\item[MLE] Minimax Rule.
\item[$\delta^L$] Best Linear Rule.
\item[$\delta^B$] Bayes Rule.
\item[$\hat \delta_{\hat G}^B$] Bayes Rule with NPMLE $\hat G$.
\item[$\delta^H$] Huber Modified Bayes Rule.
\item[$\delta^M$] Mallows Modified Bayes Rule.
\item[$\hat \delta_{\hat G}^H$] Huber Modified Bayes Rule with NPMLE $\hat G$.
\item[$\hat \delta_{\hat G}^M$] Mallows Modified Bayes Rule with NPMLE $\hat G$.
\end{description}
The last four rules are evaluated for four distinct values of $\epsilon \in \{ 0.05, 0.1, 0.2, 0.4\}$.
All the simulations are based on 500 replications, for each compound decision problem.
The linear rules work reasonably well for the Gaussian and Uniform settings, however they
perform poorly in the two-point setting. The cost of the Hodges-Lehmann restricted priors
is modest for small $\epsilon$, but not surprisingly grows substantially when $\epsilon$ is
larger. With only $n = 100$ observations, the NPMLE rule, $\delta_{\hat G}^B$, is too variable,
but for the larger sample sizes it is nearly competitive with the (oracle) Bayes rules.
For each $G_0$ we can evaluate for any $\epsilon$ the corresponding worst case pointwise risk of each rule.
The minimax rule achieves a worst case pointwise risk of one, while the Bayes rule typically has unbounded
pointwise risk. For the Mallows rule, this can be evaluated by $\mathbb{E}_{\theta^*}[(\delta^M(X)-\theta^*)^2]$
where $\theta^*$ is any mass point of $H^*$ as discussed in Section \ref{sec: computation}.
For the Huber rule, we can show that for all the $G_0$ we considered,
$\sup_\theta R(\delta^H, \theta) = 1+k_\epsilon^2$ with $k_\epsilon = \sup_{x} | (\log f_\epsilon^H(X))'|$
in which $f_\epsilon^H$ is the least favorable Huber density for the given $G_0$ and $\epsilon$.
Details along with some simulation results appear in Appendix \ref{app.C}.
\input tabs1/sim1a
\input tabs1/sim2a
\input tabs1/sim3a
\section{Robustified Gaussian Likelihoods}
Rather than robustifying the prior an alternative strategy is to robustify the
likelihood.
We will consider two variants of this: the first
following \cite{huber64} and the second following \cite{mallows78}.
The classical procedure of Huber for estimating a location parameter is easily adapted
to the Gaussian sequence compound decision problem. In place of the Gaussian likelihood
in the NPMLE problem we simply insert the Huber log likelihood with density,
\[
\varphi (u) =
\begin{cases}
(1-\epsilon) \varphi (k) \exp(-k(u-k)) & u > k\\
(1-\epsilon)\varphi(u) & |u| \leq k\\
(1-\epsilon) \varphi (k)\exp(k(u+k)) & u < -k
\end{cases}
\]
where $\epsilon$ and $k$ are linked by $2 \varphi(k)/k - 2 \Phi(-k) = \epsilon/(1-\epsilon)$.
This density is least favorable, that is has minimal Fisher information for location,
in the contamination model,
\[
\Psi_\epsilon= \{ \Psi = (1 - \epsilon) \Phi + \epsilon H \},
\]
over all symmetric distributions $H$. When $\epsilon = 1/2$ the least favorable
Huber distribution is Laplace, or double exponential, and can be viewed as least favorable
against asymmetric noise as well as symmetric.
\cite{mallows78} proposes to consider minimizing $I(\Phi * G)$ over $\mathcal{G}$, the set of all
distributions with mass $1 - \epsilon$ at zero, provides an alternative to the Huber
contamination model. Rather than assuming iid innovations each arising from the Huber
mixture model, Mallows considers an additive outlier model in which with probability $\epsilon$
innovations are standard Gaussian, but occasionally are generated by the convolution $\Phi * H$.
\begin{figure}
\begin{center}
\resizebox{0.6 \textwidth}{!}{{\includegraphics{figs/rules.pdf}}}
\end{center}
\caption{The Huber and Mallows decision rules are contrasted with the Gaussian rule
in a setting with $G = U[0,3]$. The heavier tail behavior of the Huber and Mallows
base distribution, $\varphi$ results in more aggressive shrinkage with extreme
observations discounted as the consequence of noise rather than signal. The Huber
and Mallows rules both set $\epsilon = 0.1$ for this figure.}
\label{fig.rules}
\end{figure}
In Figure \ref{fig.rules} we contrast the Huber and Mallows decision rules with the traditional
Gaussian rule. The true mixing distribution, $G$, is chosen to be $U[0,3]$ so the Gaussian rule is
itself somewhat curved, not the linear rule we would expect were $G$ itself Gaussian. In contrast
the Huber and Mallows rules impose a more aggressive form of shrinkage. With Gaussian $\varphi$
extreme observations can be confidently attributed to signal, while the heavier tailed $\varphi$
of the Huber and Mallows rules tend to attribute such observations to noise. As we have seen
previously, the Mallows rule oscillates around the Huber rule in the tails, but otherwise their
behavior is quite similar.
When the usual Gaussian $\varphi$ is replaced by either the Huber or Mallows alternative in the
nonparametric maximum likelihood estimation of $G$ solutions are accordingly more concentrated
with fewer extreme mass points. This effect accentuates the more aggressive shrinkage effect
observed in Figure \ref{fig.rules}.
Both the Mallows and Huber least favorable contamination models offer principled alternatives to the
strictly Gaussian noise model. They preserve the convexity of the underlying NPMLE
problem and therefore can be easily implemented in software. In the next section we
compare performance of several variants of these procedures for a few simulated compound decision
settings.
\subsection{Some Simulation Experience}
We consider the compound decision problem with observations generated from,
\[
Y_i = \theta_i + U_i, \quad i = 1, \dots, n,
\]
with $\theta_i$ and $U_i$ independent and each generated iidly from $G$ and $\Psi$ respectively.
There are two choices of $G$: Either $G \sim U[0,3]$ or $G \sim 0.9 \delta_0 + 0.1 \delta_3$.
And three choices of $\Psi$: $\Psi \sim 0.8 \Phi + 0.2 \Phi(\cdot / 3)$, $\Psi \sim \text{Laplace}$
and $\Psi = \Phi$, which we label Tukey, Laplace and Gauss respectively.
In Table \ref{tab.sim} we compare mean squared error performance of ten options with an infeasible oracle
procedure that ``knows'' both the $\Psi$ and $G$ distributions. The experiment has 500 replications each
with sample size $n = 500$. The competing feasible decision rules are: GLmix, the Gaussian NPMLE;
Laplace, the Laplacian NPMLE; HLmix$(\epsilon)$, the Huber NPMLE; MLmix$(\epsilon)$, the Mallows NPMLE,
with $\epsilon \in \{ 0.20, 0.10, 0.05, 0.025 \}$.
\input tabs2/sima.tex
It is evident from the table that the Gaussian NPMLE, GLmix, performs best when
the noise distribution is actually Gaussian; however, when $\Psi \neq \Phi$
it pays to consider one of the alternatives. The Mallows NPMLE procedures seem
to perform slightly better than the corresponding Huber methods, while the Laplace
NPMLE, LLmix, which can be regarded as a ``median-type'' estimator performs surprisingly
well over all the experimental settings.
\section{conclusion}
We have considered two distinct strategies for robustifying empirical Bayes decision rules for the
Gaussian sequence model. In the first motivated by the seminal paper of Hodges and Lehmann we would
like protection against deviations from an initial Bayes prior. In the second we seek
protection against non-Gaussian behavior in the noise distribution. Both strategies rely on
the classical robustness proposals of \cite{huber64} and \cite{mallows78}. Some combination of
the two strategies is obviously possible, but choice of tuning parameters remains a delicate issue.