GROS: A General Robust Aggregation Strategy.

Alejandro Cholaquidis1, Emilien Joly2 and Leonardo Moreno3

1 Centro de Matemática, Facultad de Ciencias, Universidad de la República,
Uruguay

2 Centro de Investigación en Matemáticas, CIMAT, México.

3 Departamento de Métodos Cuantitativos, Facultad de Ciencias Económicas y de
Administración, Universidad de la República, Uruguay.

###### Abstract

A new, very general, robust procedure for combining estimators in metric
spaces is introduced GROS. The method is reminiscent of the well-known median
of means, as described in [18]. Initially, the sample is divided into $K$
groups. Subsequently, an estimator is computed for each group. Finally, these
$K$ estimators are combined using a robust procedure. We prove that this
estimator is sub-Gaussian and we get its break-down point, in the sense of
Donoho. The robust procedure involves a minimization problem on a general
metric space, but we show that the same (up to a constant) sub-Gaussianity is
obtained if the minimization is taken over the sample, making GROS feasible in
practice. The performance of GROS is evaluated through five simulation
studies: the first one focuses on classification using $k$-means, the second
one on the multi-armed bandit problem, the third one on the regression
problem. The fourth one is the set estimation problem under a noisy model.
Lastly, we apply GROS to get a robust persistent diagram.

Keywords: Bandits, Median of means, Robustness, Sub-Gaussian estimator,
Topological data analysis

##  1 Introduction

The problem of combining estimators has been extensively studied in
statistics. There are recent proposals that merge regression estimators (see,
for instance, [5]), classifiers (see [12]), and density estimators (see [11]),
among others. In these scenarios, the aim is to merge the estimators to
generate one that, at least asymptotically, surpasses the best of the group.
In other instances, the aim is to derive a robust estimator. This is the case
of the well-known median of means (MOM) estimator, which aims to derive a
robust estimator of the mean of a random variable. In the MOM the data,
$\aleph_{n}$ (an i.i.d. sample of a random variable $X$), is first randomly
partitioned into $K$ groups. Subsequently, the mean of each group is computed.
The MOM estimator is then the median of these $K$ means. If the variance of
the data is assumed to be finite, this estimator is sub-Gaussian. For further
details, we refer to [18], [24], and the references therein. For the case of
random vectors, the so-called median of means tournament is introduced in
[29], where is proved to be sub-Gaussian. In [36] it is proved that the median
of means tournament has break-down point $\lfloor(K-1)/2\rfloor/n$, where
$\lfloor x\rfloor$ denotes the floor of $x$ and $n$ is the sample size.

The use of these estimators has proven to be valuable in varied statistical
scenarios, such as in machine learning, see [28]. In these contexts, it is
advisable to consider estimators that, without removing outliers, do not
reduce their precision. Robust statistics point in this direction, see [30]
and [2].

Following this idea of dividing into $K$ groups, calculating the estimator in
each group, and then combining them, we will introduce a new way to combine
the estimators in order to obtain, in a very general framework, a new one
that, under the assumption that the estimators by group are independent, turns
out to be sub-Gaussian (see Theorem 2 below). This new strategy, in what
follows: GROS, has breakdown point $\lceil K/2\rceil/n$, where $\lceil
x\rceil$ denotes the ceiling of $x$, see Section 2. The only assumption we
make is that the original sample comes from a random variable with finite
variance and that the space where the group estimators take their values is a
separable and complete metric space.

While the combined estimator requires solving a minimization problem in a
metric space, we prove that if it is minimized on the sample of the group
estimators, an appropriate candidate is obtained. We also determine how much
is lost by this choice, see Section 3.

Due to the immense generality of GROS , it can be applied to various areas of
statistics where robustness plays a key role. Furthermore, the space in which
the estimators reside doesn’t need to be a metric space: a pseudometric
suffices. This permits the consideration of estimators in the space of bounded
subsets of $\mathbb{R}^{d}$ equipped with the Hausdorff distance or the
measure distance, which is the case when the object to be estimated is, for
instance, a set. This space will be used in our fourth simulation example, see
subsection 4.4.

We have chosen to present five problems to demonstrate its performance,
comparing it with techniques explicitly crafted for these specific issues.
Some of these techniques were already designed to yield robust estimators.
Specifically, we treat:

  * •

The traditional clustering problem;

  * •

The multi-armed bandit problem with heavy-tailed rewards;

  * •

Regression in the presence of noisy data;

  * •

The estimation of a convex set when dealing with a noisy sample;

  * •

An application to topological data analysis.

It is worth noting that the fourth problem cannot be successfully treated
using conventional methods such as convex hulls or $r$-convex hulls, as we
will see.

As expected, in all cases the performance of GROS is noticeably better than
the proposals that do not consider the presence of outliers, see Section 4.
Moreover, the performance is good compared to methods that do consider the
presence of outliers, even surpassing some in certain cases.

The structure of the paper is outlined as follows: Section 2 introduces GROS
within a broad context and examines its robustness properties. Section 3
presents a modification of GROS to simplify its computational aspects. In
Section 4, the application of GROS across five problems is discussed. The
paper concludes with Section 5, where the findings and implications of the
study are elaborated.

##  2 Robust aggregation of weakly convergent estimators

In this section, we define and study a new proposal of a robust estimator
based on the aggregation of estimators. We assume given a sample
$\aleph_{n}=\\{X_{1},\dots,X_{n}\\}$ of i.i.d. random elements with common
distribution $P$. Let $\mu$ be a certain characteristic of $P$. We assume that
$\mu$ belongs to a complete and separable metric space $\mathcal{M}$ endowed
with a metric $d$. As we mentioned in the introduction, all the results in
this paper remains true if $d$ is a pseudo metric. In the context of robust
estimation, one goal is to obtain sub-Gaussian type inequalities for the
deviation of an estimator $\widehat{\mu}$ from $\mu$. A common way of defining
a robust, distribution-free estimator is to make $K$ disjoint groups out of
$\aleph_{n}$, hence, to create a collection of $K$ independent estimators
$\mu_{1},\dots,\mu_{K}\in\mathcal{M}$. We define the robust aggregation of
$\mu_{1},\dots,\mu_{K}$, GROS, by

$\mu^{*}=\operatorname*{argmin}_{\nu\in\mathcal{M}}\min_{I:|I|>\frac{K}{2}}\max_{j\in
I}d(\mu_{j},\nu).$ (1)

The minimization is taken over all the possible subsets $I$ of
$\\{1,\dots,K\\}$ that contain at least $\lfloor K/2\rfloor+1$ indices.
However, it is easy to see that it is enough to minimize over all possible
subsets $I$ whose cardinality, $|I|$, fulfills $|I|=\lfloor K/2\rfloor+1$.
Observe that, for any $\nu\in\mathcal{M}$,

$\min_{I:|I|>\frac{K}{2}}\max_{j\in I}d(\mu_{j},\nu)=:d(\nu,\nu_{(\lfloor
K/2\rfloor+1)-\text{NN}}),$

where $\nu_{(\lfloor K/2\rfloor+1)-\text{NN}}$ denotes the $(\lfloor
K/2\rfloor+1)$-nearest neighbor of $\nu$ in $\mu_{1},\dots,\mu_{K}$. This last
quantity is a measure of the depth of $\nu$ inside the set
$\mu_{1},\dots,\mu_{K}$. Then, $\mu^{*}$ is the point with the least depth
from all the candidates $\nu\in\mathcal{M}$. In full generality, the set of
minimizers of (1) may not be unique. In that case, we still denote by
$\mu^{*}$ one of the minimizers arbitrarily chosen. Since we assumed that
$\mathcal{M}$ is a Polish space, it is easy to see that
$\mu^{*}\neq\emptyset$.

######  Remark.

A natural generalization of $\mu^{*}$ is to define, for $q\in[1/2,1]$,

$\mu_{q}^{*}=\operatorname*{argmin}_{\nu\in\mathcal{M}}\min_{I:|I|>Kq}\max_{j\in
I}d(\mu_{j},\nu).$

All the results we present are for $q=1/2$ but they remain true for any
$q\in[1/2,1]$.

As we said in the Introduction, we aim to combine the estimators
$\mu_{1},\dots,\mu_{k}$ in a robust way. More precisely, let us recall the
definition of finite-sample breakdown point introduced by Donoho (see [19]).

######  Definition 1.

Let $\mathbf{x}=\\{x_{1},\dots,x_{n}\\}$ be a dataset, $\theta$ an unknown
parameter lying in a metric space $\Theta$, and
$\widehat{\theta}_{n}=\widehat{\theta}_{n}(\mathbf{x})$ an estimator based on
$\mathbf{x}$. Let $\mathcal{X}_{p}$ be the set of all datasets $\mathbf{y}$ of
size $n$ having $n-p$ elements in common with $\mathbf{x}$:

$\mathcal{X}_{p}=\\{\mathbf{y}:|\mathbf{y}|=n\text{ and
}|\mathbf{x}\cap\mathbf{y}|=n-p\\}.$

Then, the breakdown point of $\widehat{\theta}_{n}$ at $\mathbf{x}$ is
$\epsilon_{n}^{*}(\widehat{\theta}_{n},\mathbf{x})=p^{*}/n,$ where

$p^{*}=\max\\{p\geq
0:\forall\mathbf{y}\in\mathcal{X}_{p},\widehat{\theta}_{n}(\mathbf{y})\textrm{
is bounded and}\\\ \text{bounded away from the boundary}\
\partial\Theta,\textrm{ if }\partial\Theta\neq\emptyset\\}.$

From (1) it follows easily that the finite-sample breakdown point of $\mu^{*}$
is $\lceil K/2\rceil/n$, which is the same order obtained in [36] for the MOM
aggregation strategy mentioned in the Introduction.

The following lemma states that if there exists an $\eta$ for which at least
$K/2$ of the $\mu_{i}$ are at a distance at most $t$ from $\eta$, then any
minimum in (1) is at a distance at most $2t$ from $\eta$. This, as we will
see, implies the robustness and sub-Gaussianity of the estimator (1). Let us
write $[K]=\\{1,\dots,K\\}$.

######  Lemma 1.

Assume that there exist an $\eta\in\mathcal{M}$ and an $I\subset[K]$ with
$|I|>K/2$ such that for all $j\in I$, $d(\mu_{j},\eta)\leq t$. Then,
$d(\mu^{*},\eta)\leq 2t$.

###### Proof.

By hypothesis, there exists a set $I$ of cardinality at least $K/2$ such that
$\max_{j\in I}d(\mu_{j},\eta)\leq t$. Since $\mu^{*}$ is a minimizer of (1),
there exists a set $I_{0}$ (a priori different from $I$) with $|I_{0}|>K/2$
such that $\max_{j\in I_{0}}d(\mu_{j},\mu^{*})\leq t.$ Now, note that
$|I|+|I_{0}|>K$ and so there exists $j\in I\cap I_{0}$ such that
$d(\mu^{*},\eta)\leq d(\mu^{*},\mu_{j})+d(\mu_{j},\eta)\leq 2t$, which
concludes the proof. ∎

Lemma 1 can be applied when $\eta=\mu$, in which case if a group of at least
$K/2$ estimators is reasonably close to the objective $\mu$, then $\mu^{*}$
itself is reasonably close. Such an estimator is robust to outliers since this
effect will not be altered by the bad behavior of up to $K/2-1$ estimators.
This lemma is a technical fact that will allow us to use the so called
binomial argument. Indeed, assume that $\mu$ is such that for any $0<p<1/2$,
there exists $t=t(n,K)$ such that for all $k\in 1,\dots,K$,
$\mathbb{P}\left(d(\mu_{k},\mu)>t\right)\leq p$. Since the estimators
$\mu_{1},\dots,\mu_{K}$ are independent,

$\displaystyle\mathbb{P}\left(d(\mu^{*},\mu)>2t\right)$
$\displaystyle\leq\mathbb{P}(\exists I:|I|\geq\lfloor K/2\rfloor\ \text{ and
}\forall i\in I,d(\mu_{i},\mu)>t)$
$\displaystyle=\mathbb{P}\Bigg{(}\sum_{k=1}^{K}\mathbb{I}_{\\{d(\mu_{k},\mu)>t\\}}\geq\lfloor
K/2\rfloor\Bigg{)}$ $\displaystyle\leq\mathbb{P}\left(B_{K,p}\geq\lfloor
K/2\rfloor\right)$ $\displaystyle\leq e^{\frac{-2(\lfloor K/2\rfloor-
Kp)^{2}}{K}},$ (2)

where $B_{K,p}$ denotes a random variable with binomial distribution, with
parameters $K$ and $p$. Let us assume that the $\mu_{i}$ are identically
distributed. Then, if we choose $p=1/4$ and $K=\lceil
8\log(\delta^{-1})\rceil$, using the fact that $\lfloor\lceil
x\rceil/2\rfloor\leq\lceil x\rceil/2$, we get from (2) together with Markov’s
inequality

$\mathbb{P}\left(d(\mu^{*},\mu)>4\sqrt{\mathbb{E}d^{2}(\mu_{1},\mu)}\right)\leq\delta.$

Lastly, we have proved the following theorem:

######  Theorem 2.

Let $\aleph_{n}=\\{X_{1},\dots,X_{n}\\}$ be an i.i.d. sample from a
distribution $P$. Let $\mu\in\mathcal{M}$ a certain characteristic of $P$
where $(\mathcal{M},d)$ is a metric space. Let $\delta>0$ and $K=\lceil
8\log(\delta^{-1})\rceil$. We split $\aleph_{n}$ into $K$ disjoint groups
(assume that $n$ is large enough to guarantee that $n/K=\ell\in\mathbb{N}$),
and create $K$ estimators $\mu_{1},\dots,\mu_{K}$ of $\mu$. Let $\mu^{*}$ be
the aggregation defined by (1). Then,

$\mathbb{P}\left(d(\mu^{*},\mu)>4\sqrt{\mathbb{E}d^{2}(\mu_{1},\mu)}\right)\leq\delta.$
(3)

######  Remark.

In Theorem 2, the restriction that $n$ is large enough to guarantee
$n/K=\ell\in\mathbb{N}$ is purely technical, to ensure that the
$\mu_{1},\dots,\mu_{K}$ are identically distributed and then to get the clean
expression (3). If this is not the case, that is, the groups are unbalanced,
the estimators $\mu_{1},\dots,\mu_{K}$ are independent but not necessarily
identically distributed, and the obtained bound is

$\mathbb{P}\left(d(\mu^{*},\mu)>4\sqrt{\max_{i=1,\dots,K}\mathbb{E}d^{2}(\mu_{i},\mu)}\right)\leq\delta.$
(4)

###  2.1 Mis-specification of the set $\mathcal{M}$

Lemma 1 uses the fact that $\eta$ is inside the set $\mathcal{M}$. Now, assume
given a imperfect set $\widetilde{\mathcal{M}}$ that is mis-specified in the
sense that
$d(\eta,\widetilde{\mathcal{M}})=\inf_{\nu\in\widetilde{\mathcal{M}}}d(\eta,\nu)=\epsilon>0$.
The true parameter of interest does not belong to the set of features
$\widetilde{\mathcal{M}}$. The minimization is then given by

$\widetilde{\mu}=\operatorname*{argmin}_{\nu\in\widetilde{\mathcal{M}}}\min_{I:|I|>\frac{K}{2}}\max_{j\in
I}d(\mu_{j},\nu).$ (5)

######  Lemma 3.

Assume that there exist an $\eta\in\mathcal{M}$ and an $I\subset[K]$ with
$|I|>K/2$ such that for all $j\in I$, $d(\mu_{j},\eta)\leq t$. Then,
$d(\widetilde{\mu},\eta)<2t+\epsilon$.

###### Proof.

By the triangle inequality $\forall\delta>0$, there exists an
$\eta_{\delta}\in\widetilde{\mathcal{M}}$ such that

$\max_{j\in I}d(\mu_{j},\eta_{\delta})\leq t+\epsilon+\delta.$

So there exists a set $I_{0}$ of cardinality $|I_{0}|>K/2$, such that
$\max_{j\in I_{0}}d(\mu_{j},\widetilde{\mu})\leq t+\epsilon+\delta.$

Since $|I|+|I_{0}|>K$, there exists $j_{0}\in I\cap I_{0}$. Then,
$d(\widetilde{\mu},\eta)\leq
d(\widetilde{\mu},\mu_{j_{0}})+d(\mu_{j_{0}},\eta)\leq 2t+\epsilon+\delta$.
Since this holds for all $\delta>0$, it follows that
$d(\widetilde{\mu},\eta)<2t+\epsilon$. ∎

The following corollary is a direct consequence of Lemma 3.

######  Corollary 2.1.

Let $\aleph_{n}=\\{X_{1},\dots,X_{n}\\}$ be an i.i.d. sample from a
distribution $P$. Let $\mu\in\mathcal{M}$ be a certain characteristic of $P$,
where $(\mathcal{M},d)$ is a metric space. We assume given a set
$\widetilde{\mathcal{M}}\subset\mathcal{M}$ (possibly random) such that
$d(\eta,\widetilde{\mathcal{M}})=\epsilon>0$. Let $\delta>0$ and $K=\lceil
8\log(\delta^{-1})\rceil$. We construct the $K$ disjoint groups and $K$
estimators $\mu_{1},\dots,\mu_{K}$ of $\mu$ as in Theorem 2. Then,

$\mathbb{P}\left(d(\widetilde{\mu},\eta)>4\sqrt{\mathbb{E}d^{2}(\mu_{1},\mu)}+\epsilon\right)\leq\delta.$
(6)

##  3 Computational aspects

Equation (1) supposes that one is able to find minimizers of a complex
functional on the metric space $\mathcal{M}$, which is often an unfeasible
problem. To simplify that task, one can restrict the minimization to the set
of estimators $\mu_{1},\dots,\mu_{K}$. That is, we find the index $j^{*}$ such
that

$j^{*}=\operatorname*{argmin}_{j=1,\dots,K}\min_{I:|I|>K/2}\max_{i\in
I}d(\mu_{i},\mu_{j}).$ (7)

The next lemma and theorem state that $\mu_{j^{*}}$ has the same sub-Gaussian
type bound (up to a constant) as $\mu^{*}$.

######  Lemma 4.

Assume that there exists an $I\subset[K]$ such that $|I|>K/2$, and for all
$j\in I$, $d(\mu,\mu_{j})\leq t$. Then $d(\mu_{j^{*}},\mu)\leq 3t$.

###### Proof.

Let $I\subset[K]$ be such that $|I|>K/2$, and for all $j\in I$,
$d(\mu,\mu_{j})\leq t$, we have $d(\mu_{i},\mu_{j})\leq
d(\mu_{i},\mu)+d(\mu,\mu_{j})\leq 2t$ for all $i,j\in I$. Then, there exists
an $I_{0}$ with cardinality at least $K/2$ such that
$d(\mu_{j^{*}},\mu_{i})\leq 2t$ for all $i\in I_{0}$. Since $|I|+|I_{0}|>K$,
there exists $j_{0}\in I\cap I_{0}$. Lastly, $d(\mu_{j^{*}},\mu)\leq
d(\mu_{j^{*}},\mu_{j_{0}})+d(\mu_{j_{0}},\mu)\leq 3t$. ∎

By means of Lemma 4, it is possible to give a practical version of Theorem 2.

######  Theorem 5.

Assume the hypotheses of Theorem 2. Let
$\mu_{j^{*}}\in\\{\mu_{1},\dots,\mu_{K}\\}$ defined by the optimization (7).
Then,

$\mathbb{P}\left(d(\mu_{j^{*}},\mu)>6\sqrt{\mathbb{E}d^{2}(\mu_{1},\mu)}\right)\leq\delta.$
(8)

Note that Equation (8) is the same as Equation (3) except that it has the
constant 6 in the right-hand side. This shows that the practical version of
the estimator, $\mu_{j^{*}}$, has essentially the same rate of convergence as
$\mu^{*}$, but it can be deteriorated by a constant factor.

##  4 Some applications of GROS

###  4.1 Classification by $k$-means

One of the most popular procedures for determining clusters in a dataset is
the $k$-means method. Although the precursors of this algorithm were MacQueen
in 1967, see [31], and Hartigan in 1978, see [23]. Pollard in 1981 [34] proved
the strong consistency of the method and in 1982, in [35], determined its
asymptotic distribution.

Given $k$, the $k$-means clustering procedure partitions a set
$\\{x_{1},\ldots,x_{n}\\}\subset\mathbb{R}^{d}$ into $k$ groups as follows:
first $k$ cluster centres $a_{j}$ are chosen, in such a way as to minimise

$W_{n}=\frac{1}{n}\sum_{l=1}^{n}\min_{1\leq j\leq k}\|x_{l}-a_{j}\|^{2}.$

Then it assigns each $x_{l}$ to its nearest cluster centre. In this way, each
centre acquires a subset $C_{l}$ as its associated cluster. The mean of the
points in $C_{l}$ must equal $a_{l}$, otherwise $W_{n}$ could be decreased by
replacing $a_{l}$ with the cluster mean in the first instance, and then
reassigning some of the $x$’s to their new centres. This criterion is then
equivalent to that of minimising the sum of squares between clusters.

The standard $k$-means algorithm starts from a set
$a_{1}^{(1)},\ldots,a_{k}^{(1)},$ and alternates, up to a stopping criterion,
two steps: Assignment step: Assigns each observation $x_{j}$ to the cluster
whose centre is the closest one. Update step: Recalculate means (centroids)
for the observations assigned to each cluster, by averaging the observations
in each cluster.

These $k$-centres are not robust to the presence of outliers, nor to the
distribution’s possession of heavy tails. There are several proposals in the
literature that seek to make the $k$-means algorithm more robust. An example
is the $k$-medoids algorithm (called PAM), see [25, 26]. In this algorithm, in
the second step, one chooses that point in the cluster which minimises the sum
of the distances to the remaining observations. These points are called the
medoids.

Another proposal for a robust version of $k$-means, TClust, is developed in
[16]. It is based on an $\alpha>0$ trimming of the data. This trimming is
self-determined by the data and aims to mitigate the impact of extreme data.

We propose a simple modification to the second step of the $k$-means
algorithm, which will be referred to as RobustkM. The idea is very simple:
instead of calculating the centroids by taking the arithmetic mean in each
group, the centroids are determined using (7).

####  4.1.1 Simulations

To evaluate the performance of this proposal, we run a small simulation study.
In all cases, the data consist of an i.i.d. sample $X_{1},\ldots,X_{n}$, whose
common distribution, $F_{X}$, is given by the mixture of three bi-variate
Student distributions. More precisely,

$F_{X}(x)=0.45T\left(x,\mu_{1},\nu,\Sigma_{1}\right)+0.45T\left(x,\mu_{2},\nu,\Sigma_{2}\right)+0.1T\left(x,\mu_{3},\nu,\Sigma_{3}\right),$
(9)

where $T\left(x,\mu,\nu,\Sigma\right)$ denotes the bi-variate cumulative
Student distribution function with mean $\mu\in\mathbb{R}^{2}$, variance and
covariance matrix $\Sigma$, and $\nu=2$ degrees of freedom.

In the examples, we chose $\mu_{1}=(6,0)$, $\mu_{2}=(-6,0)$, $\mu_{3}=(0,6)$,
$\Sigma_{1}=\Sigma_{2}=\bigl{(}\begin{smallmatrix}3&0\\\
0&3\end{smallmatrix}\bigl{)}$ and
$\Sigma_{3}=\bigl{(}\begin{smallmatrix}4&1\\\ 1&9\end{smallmatrix}\bigl{)}$.

Figure 1 shows a simulation with $n=1000$ points. It can be seen that the
dispersion of the third group makes the clustering problem more difficult.

(a)  Figure 1: Simulation of 1000 observations of the multivariate Student
mixture (9). Observations are colored according to the component of the
mixture which the data comes from.

Let $\pi$ be a permutation of the set of labels $\\{1,2,3\\}$. Denote by
$C(x)\in\\{1,2,3\\}$ the true (unknown) label, and
$\widehat{C}(x)\in\\{1,2,3\\}$ the label assigned by the algorithm to
observation $x$. Then the classification error is given by

$\min_{\pi}\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{\Big{\\{}C(x_{i})\neq\pi(\widehat{C}(x_{i}))\Big{\\}}}.$
(10)

Figure 2 shows the performance of RobustkM (with $K=10$), $k$-means, PAM, and
TClust (with $\alpha=0.01$), over $1000$ replications. In TClust, the trimmed
data (at the end of the algorithm) are assigned to the nearest centres. This
toy example shows that the proposed algorithm is a competitive alternative to
other methods that “robustify” $k$-means.

(a)  Figure 2: Box plot of classification errors, according to (10), of
$K$-means, TClust, PAM and RobustKM over $1000$ replicates.

###  4.2 Bandits

The bandit problem was first proposed by Thompson in 1933, and has been
recently been gaining increasing attention. There have been several
adaptations to various frameworks, as evidenced by comprehensive surveys, such
as [7] and [9], which cover a broad spectrum of applications.

According to [27], a bandit problem can be described as a sequential game
played between a learner and an environment, spanning over a number of rounds,
$T$, which is a positive integer. In each round, denoted by $t$, the learner
first chooses an action (denoted by $\mathcal{A}_{t}$), from a defined set of
actions (denoted by $\mathcal{A}$). Following this, the environment unveils a
reward $X_{t}$, which is a real number. Given the action $\mathcal{A}_{t}$,
the reward $X_{t}$ is assumed to be independent of the past.

In the stochastic $L$-armed bandit problem, there are $L$ possible actions.
The rewards are derived from a set of distributions denoted by $P_{1}$ through
$P_{L}$ (not depending on $t$). The aim of the learner is to identify an arm,
denoted by $i^{*}$, whose distribution yields the highest mean, denoted by
$\mathbb{E}(P_{i^{*}})$. This is equivalent to minimizing the regret

$R_{T}=T\mathbb{E}(P_{i^{*}})-\sum_{t=1}^{T}\mathbb{E}(P_{\mathcal{A}_{t}}).$

A classic example of such an algorithm is the well-known “upper confidence
bands” (UCB). Let us recall that, in the UCB, the first choice
$\mathcal{A}_{1}$ is made at random, and for $t\in\\{1,\dots,T-1\\}$,

$\mathcal{A}_{t+1}=\underset{j\in\\{1,\dots,L\\}}{\text{argmax}}\ \
\frac{X_{1}\mathbb{I}_{\\{\mathcal{A}_{1}=j\\}}+\dots+X_{t}\mathbb{I}_{\\{\mathcal{A}_{t}=j\\}}}{N_{t,j}}+\sqrt{\frac{\log(t)}{N_{t,j}}}$
(11)

where $X_{i}\mathbb{I}_{\\{\mathcal{A}_{1}=j\\}}$ is the reward of arm $j$ at
time $i$ if this arm is chosen at that time and
$N_{t,j}=\sum_{s=1}^{t}\mathbb{I}_{\\{\mathcal{A}_{s}=j\\}}$ is the number of
times that the arm $j$ was chosen.

####  4.2.1 Heavy-tailed bandits

To get a logarithmic regret for the UCB (i.e., $R_{T}/\log(T)\to C$), it is
usually assumed that the distributions $P_{i}$ are sub-Gaussian. This can be
weakened to require the distributions to have finite moment generating
function, see [3]. To overcome this limitation, instead of just estimating the
mean at each step $t$ of the UCB, in [8] it is proposed to use a robust
estimator, $\widehat{\mu}_{i,t}$, of $\mathbb{E}(P_{i})$, that fulfills the
following assumption:

Assumption 1 Let $\epsilon\in(0,1]$ be a positive parameter, and let $c,v$ be
a positive constant. Let $X_{1},\dots,X_{n}$ be i.i.d. random variables with
finite mean $\mu$. Suppose that for all $\delta\in(0,1)$ there exists an
estimator $\widehat{\mu}=\widehat{\mu}(n,\delta)$ such that, with probability
at least $1-\delta$,

$\widehat{\mu}\leq\mu+v^{1/(1+\epsilon)}\left(\frac{c\log(1/\delta)}{n}\right)^{\epsilon/(1+\epsilon)}$
(12)

and also, with probability at least $1-\delta$,

$\mu\leq\widehat{\mu}+v^{1/(1+\epsilon)}\left(\frac{c\log(1/\delta)}{n}\right)^{\epsilon/(1+\epsilon)}.$
(13)

In that case we say that $\widehat{\mu}$ fulfills Assumption 1.

######  Remark.

From (3), our robust proposal $\mu^{*}$, applied on each arm, and based on
$K=\lceil 8\log(1/\delta)\rceil$ groups, fulfills Assumption 1, for
$\epsilon=1$, $v=\mathbb{V}(P_{i})$, ($\mathbb{V}(P_{i})$ being the variance
of the distribution $P_{i}$), and $c=16$. In that case, for an arm
$i=1,\dots,L$, $n=N_{t,i}/K$.

[8] proposes the following robust variant of the UCB: given
$\epsilon\in(0,1]$, for arm $i$, define $\widehat{\mu}_{i,s,t}$ as the
estimator $\widehat{\mu}(s,t^{-2})$ based on the first $s$ observed values
$X_{i,1},\dots,X_{i,s}$ of the reward of arm $i$. Define the index

$B_{i,s,t}=\widehat{\mu}_{i,s,t}+v^{1/(1+\epsilon)}\left(\frac{c\log(t^{2})}{s}\right)^{\epsilon/(1+\epsilon)}$

for $s,t\geq 1$ and $B_{i,0,t}=+\infty$. Then, at time $t$ draw an arm
maximizing $B_{i,N_{i,t-1},t}$.

We propose the following algorithm. First we choose the arms at random from
$t=1,\dots,t_{0}$, where $t_{0}$ guarantees that for all $i=1,\dots,L$,
$N_{t_{0},i}/\lceil 8\log(t_{0}^{2})\rceil\geq 1$. We compute, for each arm,
the estimator (1), denoted by $\mu^{*}_{t_{0},i}$, where we split the
$N_{t_{0},i}$ observations into $K=\lceil 8\log(t_{0}^{2})\rceil$ groups, and
compute the mean of each group. Define the index

$\mathcal{B}_{i,N_{t_{0},i},t_{0}}=\mu^{*}_{t_{0},i}+4\sqrt{\widehat{\mathbb{V}}(P_{i})}\left(\frac{\log(t_{0}^{2})}{N_{t_{0},i}}\right)^{1/2},$

where $\widehat{\mathbb{V}}(P_{i})$ is any consistent estimation of
$\mathbb{V}(P_{i})$, or an upper bound. At time $t_{0}+1$, we choose the arm
that maximize $\mathcal{B}_{i,N_{t_{0},i},t_{0}}$, at time $t_{0}+2$ we choose
the arm that maximize $\mathcal{B}_{i,N_{t_{0}+1,i},t_{0}+1}$, and so on.

Proposition 1 in [8] proves that this algorithm attains logaritmic regret.
More precisely,

$R_{T}\leq\sum_{i:\Delta_{i}>0}\left(32\left(\frac{\mathbb{V}(P_{i})}{\Delta_{i}}\right)\log(T)+5\Delta_{i}\right),$

where $\Delta_{i}=\mathbb{E}(P_{i^{*}})-\mathbb{E}(P_{i})$.

####  4.2.2 Simulations

As a toy example, let us consider the classical two-armed bandit problem and
rewards given by $X_{t}|\mathcal{A}_{t}=j\sim\mu_{j}+S(3),$ for $j=1,2$, where
$S(3)$ is a random variable following Student’s distribution with $3$ degrees
of freedom, $\mu_{1}=7$, and $\mu_{2}=8$. In our algorithm, indicated by RUCB
in Figure 3, the first $t_{0}=40$ arms are chosen at random.

The results are shown in Figure 3, where the red dotted horizontal line
$(y=8)$ is the maximum expected gain. The dashed lines corresponds to the mean
rewards of the UCB (orange) and the RUCB (blue) respectively. As can be seen,
it takes the RUCB algorithm about 120 steps to outperform the UCB algorithm,
and the difference grows larger as the number of steps increases.

(a)  Figure 3: Cumulative gains over 500 replications, for $t=1,\ldots,750$.
The red dotted horizontal line ($y=8$) is the maximum expected gain. The black
dotted vertical line ($x=40$) indicates the number of random warm-up runs in
the RUCB algorithm. The dashed lines depict the mean reward of the UCB
(orange) and RUCB (blue) algorithms.

###  4.3 Robust regression

A current topic of interest in statistics is that of regression models in the
presence of noise. It is known that a small fraction of outliers can cause
severe biases in classical regression estimators. These classical models
generally assume additive residuals in the model with finite second moment
(e.g., Gaussian). An extensive review of robust outlier estimation methods for
nonparametric regression models is provided in [38].

We tackle the problem of estimating the function $m:\mathcal{X}\to\mathbb{R}$,
from an i.i.d. sample of a random element
$(X,Y)\in\mathcal{X}\times\mathbb{R}$, that satisfies the model

$Y=m(X)+\epsilon$

where $\epsilon$ is a noise term such that
$\mathbb{E}\left[\epsilon|X\right]=0$. To keep the following discussion fairly
simple, we will only cover the case $\mathcal{X}=\mathbb{R}$. In the context
where very little information is known about $m$, a general estimator is often
given by the kernel estimator

$\widehat{m}(x)=\frac{\sum_{i=1}^{n}K_{h}(X_{i}-x)Y_{i}}{\sum_{i=1}^{n}K_{h}(X_{i}-x)},$
(14)

where $K_{h}(x)=h^{-1}\phi(h^{-1}x)$ and $h>0$ is some bandwidth parameter.
The function $\phi$ is non-negative and such that
$\int_{\mathbb{R}}\phi(x)dx=1$. In general, the continuity of the function $m$
is enough to get the consistency of the kernel estimator $\widehat{m}$ and a
light tailed behavior for $\epsilon$ gives sub-exponential deviation bounds
for the estimator around its mean value $m$. Nevertheless, in various concrete
settings, one can face heavy tailed distributions for $\epsilon$ in such a way
that the estimator proposed in (14) becomes highly unstable. Indeed, the
presence of (virtually) one outlier is enough to drive the estimator towards
values very distant from $m$. One natural way to measure the quality of the
estimator $\widehat{m}$ is to consider the $L_{2}$ distance

$d_{2}(\widehat{m},m)=\mathbb{E}\left[|\widehat{m}(X)-m(X)|^{2}\right]^{1/2}.$

In this context, it is possible to introduce the robust version of the
estimator by considering GROS with the distance $d_{2}$. In practice, the
distance $d_{2}$ is actually intractable. To overcome this difficulty, it is
common to use a discretization on a mesh for the space $\mathcal{X}$ and
approximate the integral by a Riemann sum. This estimator verifies a bound as
in Theorem 2 for the associated distance. The estimator $m^{*}$ is, in the
sense of Definition (1), the best candidate in the class of the estimator
$m_{1},\dots,m_{K}$ constructed on the disjoint groups $G_{1},\dots,G_{K}$ of
data points. These estimators $m_{1},\dots,m_{K}$ are defined by

$m_{j}(x)=\frac{\sum_{i\in G_{j}}K_{h}(X_{i}-x)Y_{i}}{\sum_{i\in
G_{j}}K_{h}(X_{i}-x)},$

so that these estimators are independent. Following subsection 2.1, the
following estimator for $\widehat{m}$ can be proposed:

$\widehat{m}=\operatorname*{argmin}_{m\in\widetilde{\mathcal{M}}}\min_{I:|I|>\frac{K}{2}}\max_{j\in
I}d_{2}(m_{j},m),$

where the $(m_{j})_{j=1}^{K}$ are the kernel estimates of $m$ on each of the
$K$ groups, and the set $\widetilde{\mathcal{M}}$ denotes the set of functions
that are piece-wise equal to one of the $m_{i}$. The difference from the naive
definition is that the set of minimization does not end with an estimator of
the form $m_{i^{*}}$, which allows avoiding choosing a function that may be a
good fit in some regions of the space but sensitive to outliers in other
parts.

In order to fully use the result of Equation (4), we cite a result that gives
an upper bound for the mean squared error for any of the estimators $m_{j}$
previously defined. In chapter 5 of [22, Theorem 5.2], we get that if $m$ is
$\lambda$-Lipschitz, and that $\text{Var}(Y|X=x)\leq\sigma^{2}$ for all
$x\in\mathcal{X}$, then

$\mathbb{E}\left[\|m_{1}-m\|^{2}\right]\leq
c\left(\frac{\sigma^{2}+\sup_{x\in\mathcal{X}}m(x)^{2}}{(n/K)h^{d}}\right)+\lambda^{2}h^{2}.$

This induces a choice of $h$ of the form

$h=c^{\prime}\left(\frac{K}{n}\right)^{\frac{1}{d+2}}.$

The upper bound then takes the form of

$\mathbb{E}\left[\|m_{1}-m\|^{2}\right]\leq
c^{\prime\prime}\left(\sigma^{2}+\sup_{x\in\mathcal{X}}m(x)^{2}\right)^{2/d+2}\times\left(\frac{n}{K}\right)^{-2/d+2},$

which can be directly plugged into the main bound of the theorem.

####  4.3.1 Simulations

We compare, by means of simulations, the performance of GROS with some
classical and robust regression alternatives proposed in the literature.

In this example we consider a sample $X_{1},\ldots,X_{1000}$ with uniform
distribution on $[0,5]$. Let $m(x)=4\sin(3x)$ and suppose that the centered
noise follow the skew-normal Student distribution, see [21, 4], whose density
is defined as follows: denote by $t(x,\kappa,\nu)$ the density of the non-
standardized Student’s distribution with location $\kappa$ and and $\nu$
degrees of freedom, and with cumulative distribution function
$T(x,\kappa,\nu)$. Then, the density of the skew-normal Student is

$f(x;\kappa,\nu,\sigma,\xi)=\frac{1}{\sigma}t(x/\sigma,\kappa,\nu)T\left(\frac{\xi
x}{\sigma},\kappa,\nu\right),$

where $\sigma>0$ denotes the scale parameter. The slant (or skewness
parameter) is $\xi$.

We will write NW for the non-parametric Nadaraya–Watson kernel regression
estimator (14), see [32, 41]. As robust estimators, we consider the proposals
developed in [33] (ONL estimator) and [6] (SBMB estimator). Both estimators
are implemented in the R language: the first in the fields library and the
other in the RBF library. Our estimator in this context will be called RANW
(Robust Aggregation for Naradaya–Watson). For the latter, $K=12$ was
considered. Figure 4 shows the estimated regressions in each scenario and in
red the true function $f$. It is clear that the NW estimator performs poorly
in all cases.

The performance of each estimator is measured by the average distance
$d_{2}(\widehat{m},m)$ over $1000$ replicates, considering $\kappa=0$, $\nu=3$
(note that the noises have heavy tails) and varying the parameters $\sigma$
and $\xi$. For the bandwidth parameter $h$, we choose $0.2$ in all cases.

Figure 4 shows box plots of the errors for four choices of the parameters. As
can be seen, when the noise is asymmetric (i.e., $\xi\neq 1$), the best-
performing estimators are RANW and ONL, depending on the value of the scale
parameter. Note that when the distribution of the noise is symmetric (i.e.,
$\xi=1$), SBMB and RANW perform the best.

(a) $\sigma=9,\nu=3,\xi=1$

(b) $\sigma=9,\nu=3,\xi=9$

(c) $\sigma=16,\nu=3,\xi=1$

(d) $\sigma=16,\nu=3,\xi=9$

Figure 4: Box plot of classification errors (according to L2 distance) in
$1000$ replicates. The different scenarios are obtained in the skew-normal
Student distribution with $\sigma\in\\{9,16\\}$ and $\xi\in\\{1,9\\}$, fixed
$\nu=3$ and $\kappa=0$.

(a) $\sigma=9,\nu=3,\xi=1$

(b) $\sigma=9,\nu=3,\xi=9$

(c) $\sigma=16,\nu=3,\xi=1$

(d) $\sigma=16,\nu=3,\xi=9$

Figure 5: Regression functions estimated with the RANW (orange), NW (black),
ONL (light blue) and SBMB (blue) in one replicate. The true function is shown
in red. The different scenarios are obtained in the skew-normal Student
distribution with $\sigma\in\\{9,16\\}$ and $\xi\in\\{1,9\\}$, fixed $\mu=0$
and $\nu=3$.

###  4.4 Robust set estimation

Set estimation involves determining a set, or a function of that set, based on
a random sample of points. This set could represent various things, such as
the support of a probability distribution (see, for instance, [37, 10]), its
boundary (see, for instance [17], the surface area of the boundary (see, for
instance [1]), to name a few.

Within this framework, shape constraints are usually imposed, convexity being
the most well-known. However, it can be restrictive for some applications,
such as, for instance, if the set is the home-range of a species (see [13],
[14], [15] and references therein). Usually, the available data is an i.i.d.
sample of a random vector whose support is the unknown set. For classic
estimators such as the convex hull, $r$-convex hull, or cone-convex hull, any
noise in the sample, no matter how small, drastically changes the estimators.
This becomes even more pronounced in the case of convexity.

####  4.4.1 Simulations

We show the performance of our robust proposal under a noisy model, where the
aim is to estimate a convex set. More precisely, let $D(r,R)$ be the uniform
distribution on the ring in $\mathbb{R}^{2}$ with inner radius $r$ and outer
radius $R$. We simulated $2000$ i.i.d. observations of the mixture

$(1-\lambda)D(0,1)+\lambda D(1,1.25).$

The aim is to estimate $D(0,1)$ from this sample. We have chosen
$\lambda=0.01$ as the proportion of noise. The estimator, referred to as
RChull, is the one proposed in Section 3, and we considered the Hausdorff
distance between compact sets to measure the discrepancy between $D(0,1)$ and
the estimator. To build our estimator, we first split the original sample at
random into $K=20$ disjoint groups of points of size $100$, and compute the
convex hull on each group. Then, we select from the $K$ hulls
$H_{1},\dots,H_{20}$, the hull $H_{j^{*}}$, with $j^{*}$ as in (7).

To gain an insight into the improvement resulting from using this robust
procedure, we show on the left of Figure 6 the set to be estimated, the
classical estimator (the convex hull (Chull) of the whole sample) and the 20
hulls of the subsamples. The right panel shows the convex hull of the whole
sample, together with $H_{j^{*}}$.

(a)

(b)

Figure 6: In the blue dotted line the boundary of the ball of radius 1. The
sample is shown as solid black points. Outside this ball the sample generated
from $D(1,1.25)$. On the left there are shown the 20 convex hulls of the
selected subsamples (of size 100). On the right the convex hullf of the whole
sample (Chull) and the robust estimator based on the $20$-convex hulls,
(RChull).

To evaluate the performance of GROS , we run 100 replicates and estimate the
set by both methods (the classical convex hull and our robust proposal) in
each replicate. We calculate the Hausdorff distance of the estimated sets
RChull and Chull from the set D(0,1). In Figure 7 box plots of these distances
are shown. The RChull estimator outperforms the classical Chull estimator, but
as can be seen, it has larger variability.

(a)  Figure 7:

###  4.5 Robust persistent diagram

The robustness of persistent homology to perturbations in data measured by the
Hausdorff distance is well established. However, it has a high sensitivity to
outliers, as discussed in [39, 40].

In this section, we introduce a robust persistence diagram, which we call the
Robust Wasserstein Estimator, using (7).

The measure of dissimilarity between two persistence diagrams $P_{1}$ and
$P_{2}$ is quantified by the 1-Wasserstein distance $W_{1}(P_{1},P_{2})$. This
distance quantifies the cost associated with achieving the optimal alignment
of points between the two diagrams, as detailed in [20, p. 202].

####  4.5.1 Simulations

The example examines uniformly simulated data on $S^{1}$ consisting of 600
points (baseline sample). It explores two potential scenarios of sample
distortion, as depicted in Figure 8.

(a) Baseline sample (b) Scenario 1: Local perturbation (c) Scenario 2: Groups
of outliers

Figure 8:  (a) Baseline sample of 600 points uniformly distributed on $S^{1}$.
(b) Locally perturbed sample as described in Scenario 1. (c) Sample perturbed
in accordance with Scenario 2.

  * •

Scenario 1:  Local perturbation. The original sample is perturbed using
Gaussian noise centered at each sample point, with a standard deviation matrix
$0.05\times Id$.

  * •

Scenario 2: Group of Outliers: We randomly selected $90\%$ of the perturbed
sample as defined in Scenario 1. The remaining $10\%$ are derived from a
Matérn cluster process within the square region $[-0.5,0.5]^{2}$. This process
is characterized by an intensity of $3$ for the Poisson process of cluster
centers, a scale of $0.25$, and an average of $20$ points per cluster.

The persistence diagrams for the previously described tree samples were
computed. Figure 9 displays these diagrams, where $Dgm$, $Dgm_{1}$, and
$Dgm_{2}$ represent the persistence diagrams of the baseline sample, the
locally perturbed sample (Scenario 1), and the sample with outliers (Scenario
2), respectively. It is evident that the persistence diagram is more distorted
in scenarios with groups of outliers than it is with those with local
perturbations.

(a) $Dgm$ (b) $Dgm_{1}$ (c) $Dgm_{2}$

Figure 9: (a): Persistence diagram of the baseline sample. (b) and (c):
Persistence diagrams for the samples perturbed according to Scenarios 1 and 2,
respectively.

To evaluate the performance of our robust proposal (7), the sample is divided
into $K=6$ distinct groups. The robust persistence diagrams for Scenarios 1
and 2 are labeled as $RDgm_{1}$ and $RDgm_{2}$, respectively. We computed
$W_{1}(Dgm_{1},Dgm)$, $W_{1}(Dgm_{2},Dgm)$, $W_{1}(RDgm_{1},Dgm)$, and
$W_{1}(RDgm_{2},Dgm)$.

A total of 100 independent iterations of this experiment were conducted, and
box plots of the respective distances are shown in Figure 10. The results
indicate that in both scenarios, the robust estimates of the persistence
diagrams show improved performance, being closer to the baseline diagram in
terms of the Wasserstein distance.

(a)  Figure 10: Box plots illustrating the distances between the persistence
diagrams of perturbed samples and the baseline sample diagram $Dgm$, as well
as the distances between the robust persistence diagrams and the baseline
sample diagram.

##  5 Concluding remarks

We have demonstrated through simulations that GROS , which is applicable to a
broad range of problems, significantly improves upon the non-robust, problem-
specific solutions for each of the five examples treated. It is also
competitive with robust solutions designed for each specific case, even
showing some improvements. GROS proposal’s flexibility makes it applicable to
a wide variety of problems, including those already presented, as well as any
other scenario where robustness plays a crucial role. Furthermore, more
examples using only pseudo-distances may be of interest for future research.

## References

  * [1] Catherine Aaron, A. Cholaquidis, and R. Fraiman. Estimation of surface area. Electron. J. Statist. 16(2): 3751-3788, 2022.
  * [2] Catherine Aaron, Alejandro Cholaquidis, Ricardo Fraiman, and Badih Ghattas. Multivariate and functional robust fusion methods for structured big data. Journal of Multivariate Analysis, 170:149–161, 2019.
  * [3] Rajeev Agrawal. Sample mean based index policies by $o$ (log $n$) regret for the multi-armed bandit problem. Advances in Applied Probability, 27(4):1054–1078, 1995.
  * [4] Adelchi Azzalini. The Skew-Normal and Related Families. Institute of Mathematical Statistics Monographs. Cambridge University Press, 2013.
  * [5] Gérard Biau, Aurélie Fischer, Benjamin Guedj, and James D. Malley. COBRA: A combined regression strategy. Journal of Multivariate Analysis, 146:18–28, 2016.
  * [6] Graciela Boente, Alejandra Martínez, and Matías Salibián-Barrera. Robust estimators for additive models using backfitting. Journal of Nonparametric Statistics, 29(4):744–767, 2017.
  * [7] Etienne Boursier and Vianney Perchet. A survey on multi-player bandits. arXiv preprint arXiv:2211.16275, 2022.
  * [8] Sébastien Bubeck, Nicolo Cesa-Bianchi, and Gábor Lugosi. Bandits with heavy tail. IEEE Transactions on Information Theory, 59(11):7711–7717, 2013\.
  * [9] Giuseppe Burtini, Jason Loeppky, and Ramon Lawrence. A survey of online experiment design with the stochastic multi-armed bandit. arXiv preprint arXiv:1510.00757, 2015.
  * [10] A. Cholaquidis, A. Cuevas, and R. Fraiman. On Poincaré cone property. The Annals of Statistics, 42(1):255–284, 2014.
  * [11] Alejandro Cholaquidis, Ricardo Fraiman, Badih Ghattas, and Juan Kalemkerian. A combined strategy for multivariate density estimation. Journal of Nonparametric Statistics, 33(1):39–59, 2021.
  * [12] Alejandro Cholaquidis, Ricardo Fraiman, Juan Kalemkerian, and Pamela Llop. A nonlinear aggregation type classifier. Journal of Multivariate Analysis, 146:269–281, 2016.
  * [13] Alejandro Cholaquidis, Ricardo Fraiman, Gábor Lugosi, and Beatriz Pateiro-López. Set Estimation from Reflected Brownian Motion. Journal of the Royal Statistical Society Series B: Statistical Methodology, 78(5):1057–1078, 01 2016.
  * [14] Alejandro Cholaquidis, Ricardo Fraiman, Ernesto Mordecki, and Cecilia Papalardo. Level set and drift estimation from a reflected Brownian motion with drift. Statistica Sinica, 31:29–51, 2021.
  * [15] Alejandro Cholaquidis, Manuel Hernandez, and Ricardo Fraiman. Home range estimation under a restricted sampling scheme. Journal of Nonparametric Statistics, 2023. To appear.
  * [16] Juan Antonio Cuesta-Albertos, Alfonso Gordaliza, and Carlos Matrán. Trimmed $k$-means: an attempt to robustify quantizers. Ann. Statist., 25(2):553–576, 1997.
  * [17] A. Cuevas and A. Rodríguez-Casal. On boundary estimation. Advances in Applied Probability, 36(2):340–354, 2004.
  * [18] Luc Devroye, Matthieu Lerasle, Gabor Lugosi, and Roberto I. Oliveira. Sub-Gaussian mean estimators. Ann. Statist., 44(6):2695–2725, 2016.
  * [19] David L. Donoho. Breakdown properties of multivariate location estimators. Technical report, Harvard University, Boston, 1982.
  * [20] Herbert Edelsbrunner and John L. Harer. Computational Topology: An Introduction. American Mathematical Society, 2022.
  * [21] Carmen Fernández and Mark F.J. Steel. On Bayesian modeling of fat tails and skewness. Journal of the American Statistical Association, 93(441):359–371, 1998.
  * [22] László Györfi, Michael Kohler, Adam Krzyzak, Harro Walk, et al. A Distribution-free Theory of Nonparametric Regression, Volume 1. Springer, 2002.
  * [23] J.A. Hartigan. Asymptotic distributions for clustering criteria. Ann. Statist., 6:117–131, 1978.
  * [24] Emilien Joly, Gábor Lugosi, and Roberto Imbuzeiro Oliveira. On the estimation of the mean of a random vector. Electron. J. Statist., 11(1):440–451, 2017.
  * [25] Leonard Kaufman. Partitioning Around Medoids. Finding Groups in Data, 344:68–125, 1990.
  * [26] Leonard Kaufman and Peter J. Rousseeuw. Finding Groups in Data: An Introduction to Cluster Analysis. Wiley, 2009.
  * [27] Tor Lattimore and Csaba Szepesvári. Bandit Algorithms. Cambridge University Press, 2020.
  * [28] Guillaume Lecué and Matthieu Lerasle. Robust machine learning by median-of-means: Theory and practice. Ann. Statist., 48(2):906–931, 2020.
  * [29] Gábor Lugosi and Shahar Mendelson. Sub-Gaussian estimators of the mean of a random vector. Ann. Statist., 47(2):783–704, 2019.
  * [30] Ricardo A. Maronna, R. Douglas Martin, Victor J. Yohai, and Matías Salibián-Barrera. Robust Statistics: Theory and Methods (with R). Wiley, 2019.
  * [31] James, B. McQueen. Some methods for classification and analysis of multivariate observations. In: Proceedings of the 5th Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Statistics, University of California Press, Berkeley, 281-297. 1967.
  * [32] Elizbar A. Nadaraya. On estimating regression. Theory of Probability & Its Applications, 9(1):141–142, 1964.
  * [33] Hee-Seok Oh, Douglas W. Nychka, and Thomas C.M. Lee. The role of pseudo data for robust smoothing with application to wavelet regression. Biometrika, 94(4):893–904, 2007.
  * [34] David Pollard. Strong consistency of $k$-means clustering. Annals of Probability, 9(1):135–140, 1981.
  * [35] David Pollard. A central limit theorem for $k$-means clustering. Annals of Probability, 10(4):919–926, 1982.
  * [36] Daniela Rodriguez and Marina Valdora. The breakdown point of the median of means tournament. Statistics & Probability Letters, 153:108–112, 2019.
  * [37] A. Rodríguez Casal. Set estimation under convexity type assumptions. Annales de l’IHP Probabilités et statistiques, 43(6):763–774, 2007.
  * [38] Matias Salibian-Barrera. Robust nonparametric regression: Review and practical considerations. Econometrics and Statistics, 2023.
  * [39] Siddharth Vishwanath, Kenji Fukumizu, Satoshi Kuriki, and Bharath K. Sriperumbudur. Robust persistence diagrams using reproducing kernels. Advances in Neural Information Processing Systems, 33:21900–21911, 2020.
  * [40] Siddharth Vishwanath, Bharath K. Sriperumbudur, Kenji Fukumizu, and Satoshi Kuriki. Robust topological inference in the presence of outliers. arXiv preprint arXiv:2206.01795, 2022.
  * [41] Geoffrey S. Watson. Smooth regression analysis. Sankhyā: The Indian Journal of Statistics, Series A, 26(4) 359–372, 1964.