# Improved rates for Wasserstein deconvolution with ordinary smooth error in dimension one

Jérôme Dedecker, Aurélie Fischer, Bertrand Michel

► To cite this version:

Jérôme Dedecker, Aurélie Fischer, Bertrand Michel. Improved rates for Wasserstein deconvolution with ordinary smooth error in dimension one. Electronic Journal of Statistics, 2015, 9 (1), pp.234-265. 10.1214/15-EJS997. hal-00971316v3

HAL Id: hal-00971316

https://hal.sorbonne-universite.fr/hal-00971316v3

Submitted on 27 Jan 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Improved rates for Wasserstein deconvolution with ordinary smooth error in dimension one

Jérôme Dedecker*

Laboratoire MAP5 UMR CNRS 8145

Université Paris Descartes

Sorbonne Paris Cité

75006 Paris, France

e-mail: jerome.dedecker@parisdescartes.fr

Aurélie Fischer

LPMA UMR CNRS 7599

Université Paris Diderot

Sorbonne Paris Cité

75205 Paris Cedex 13 France

e-mail: aurelie.fischer@univ-paris-diderot.fr

and

Bertrand Michel

LSTA

UPMC Université Paris 06

Sorbonne Universités

75252 Paris Cedex 05 France

e-mail: bertrand.michel@upmc.fr

**Abstract:** This paper deals with the estimation of a probability measure on the real line from data observed with an additive noise. We are interested in rates of convergence for the Wasserstein metric of order $p \ge 1$. The distribution of the errors is assumed to be known and to belong to a class of supersmooth or ordinary smooth distributions. We obtain in the univariate situation an improved upper bound in the ordinary smooth case and less restrictive conditions for the existing bound in the supersmooth one. In the ordinary smooth case, a lower bound is also provided, and numerical experiments illustrating the rates of convergence are presented.

**MSC 2010 subject classifications:** Primary 62G05, 62C20.

**Keywords and phrases:** Deconvolution, Wasserstein metrics, minimax rates.

Received April 2014.

<footnote>*The authors would like to thank the Editor, the Associate Editor and the referees for their insightful comments and suggestions. The authors were supported by the ANR project TopData ANR-13-BS01-0008.</footnote>
Contents

<table><tr><td>1</td><td>Introduction</td><td>235</td></tr><tr><td>2</td><td>On the case without error</td><td>237</td></tr><tr><td>3</td><td>Upper bounds for W<sub>p</sub> in deconvolution</td><td>239</td></tr><tr><td>3.1</td><td>Construction of the estimator</td><td>239</td></tr><tr><td>3.2</td><td>First upper bounds for W<sub>p</sub><sup>p</sup>(μ̂<sub>n</sub>, μ)</td><td>241</td></tr><tr><td>3.3</td><td>Main results</td><td>242</td></tr><tr><td>4</td><td>Lower bound</td><td>245</td></tr><tr><td>5</td><td>Proof of Proposition 3.1</td><td>248</td></tr><tr><td>6</td><td>Numerical experiments</td><td>253</td></tr><tr><td>6.1</td><td>Implementation of the deconvolution estimators</td><td>253</td></tr><tr><td>6.2</td><td>Computation of Wasserstein risks for simulated experiments</td><td>256</td></tr><tr><td>6.3</td><td>Estimation of the rates of convergence</td><td>257</td></tr><tr><td>6.4</td><td>Cantor set experiment</td><td>258</td></tr><tr><td>6.5</td><td>About the bandwidth choice</td><td>261</td></tr><tr><td></td><td>References</td><td>264</td></tr></table>

1. Introduction

Consider the following convolution model: we observe $n$ real-valued random variables $Y_1, \dots, Y_n$ such that

$$Y_i = X_i + \varepsilon_i, \qquad (1)$$

where the $X_i$'s are independent and identically distributed according to an unknown probability $\mu$, which we want to estimate. The random variables $\varepsilon_i$, $i = 1, \dots, n$, are independent and identically distributed according to a known probability measure $\mu_\varepsilon$, not necessarily symmetric. Moreover we assume that $(X_1, \dots, X_n)$ is independent of $(\varepsilon_1, \dots, \varepsilon_n)$.

The purpose of the paper is to investigate rates of convergence for the estimation of the measure $\mu$ under Wasserstein metrics. For $p \in [1, \infty)$, the Wasserstein distance $W_p$ between $\mu$ and $\nu$ is given by

$$W_p(\mu, \nu) = \inf_{\pi \in \Pi(\mu, \nu)} \left( \int_{\mathbb{R}^2} |x - y|^p \pi(dx, dy) \right)^{\frac{1}{p}},$$

where $\Pi(\mu, \nu)$ is the set of probability measures on $\mathbb{R} \times \mathbb{R}$ with marginal distributions $\mu$ and $\nu$ (see [24] or [26]). The distances $W_p$ are natural metrics for comparing measures. For instance they can compare two singular measures, which is of course impossible with the functional metrics commonly used in density estimation. Convergence of measure under Wasserstein distances is an active domain of research in probability and statistics. For instance, the rate of convergence of the empirical measure under these metrics has been obtained recently by both [14] and [19] in $\mathbb{R}^d$ and also by [1] in the one-dimensional framework. Moreover, Wasserstein metrics are involved in many fields of mathematics and computer sciences. For instance, in the field of Topological Data Analysis
(TDA) [5], Wasserstein distances recently appeared to be natural metrics for controlling the estimation of geometric and topological features of the sampling measure and its support. Indeed, in [7], a distance function to measures is introduced to solve geometric inference problems in a probabilistic setting: if a known measure $\nu$ is close enough with respect to $W_2$ to a measure $\mu$ concentrated on a given shape, then the topological properties of the shape can be recovered by using the distance to $\nu$. More generally, the Wasserstein loss could be used as a guide for inferring the support (see the Cantor experiment in Section 6.4). Other results in TDA with stability results involving the Wasserstein distances can be found in [20] and [8]. In practice, the data can be observed with noise, which motivates in this framework the study of the Wasserstein deconvolution problem [4], in particular if the deconvolved measure and the “true measure” are singular.

Rates of convergence in deconvolution have mostly been considered in density estimation, for pointwise or global convergence. Minimax rates can be found for instance in [17, 2, 3] and in the monograph of [23]. In this paper, however, we shall not assume that $\mu$ has a density with respect to the Lebesgue measure. In this context, rates of convergence for the $W_2$ Wasserstein distance have first been studied for several noise distributions by [4]. Recently, [10] have obtained optimal rates of convergence in the minimax sense for a class of supersmooth error distributions, in any dimension, under any Wasserstein metric $W_p$. The result relies on the fact that lower bounds in any dimension can be deduced in this case from the lower bounds in dimension 1. Such a method cannot be used in the ordinary smooth case, where the rate of convergence depends on the dimension. As noticed by [17], establishing optimal rates of convergence in the ordinary smooth case is more difficult than in the supersmooth one, even for pointwise estimation.

A key fact in the univariate context is that Wasserstein metrics are linked to integrated risks between cumulative distribution functions (cdf), see the upper bound (5) below. In dimension 1, when estimating the density of $\mu$, optimal rates of convergence for integrated risks can be found in [16, 18]. When estimating the cdf $F$ of $\mu$, optimal rates for the pointwise and integrated quadratic risks are given in [21], where it is shown in particular that the rate $\sqrt{n}$ can be reached when the error distribution is ordinary smooth with a smoothness index less than 1/2. Concerning the pointwise estimation of $F(x_0)$, optimal rates for the quadratic risk are also given in [9], when the density of $\mu$ belongs to a Sobolev class.

The case $\beta = 0$ in the upper bound (3.9) of [21] corresponds to the case where no assumption (except a moment assumption) is made on the measure $\mu$ (in particular $\mu$ is not assumed to be absolutely continuous with respect to the Lebesgue measure). This is precisely the case which we want to consider in the present paper. However the results by [21] cannot be applied to the Wasserstein deconvolution problems for two reasons: firstly, the integrated quadratic risk for estimating a cdf is not linked to Wasserstein distances, and secondly, the estimator of the cdf of $\mu$ proposed in [21] is the cdf of a signed measure, and is not well defined as an estimator of $\mu$ for the Wasserstein metric.
In the present contribution, we propose in the univariate situation an improved upper bound for deconvolving $\mu$ under $W_p$, and a lower bound when the error is ordinary smooth. We recover the optimal rate of convergence in the supersmooth case with slightly weaker regularity conditions than in [10]. The estimator of the cdf $F$ of $\mu$ is built in two steps: firstly, as in [21], we define a preliminary estimator through a classical kernel deconvolution method, and secondly we take an appropriate isotone approximation of this estimator. For controlling the random term, we use a moment inequality on the cdfs, which is due to [15]. To be complete, we show in Section 4 that for $p > 1$, the Wasserstein deconvolution problem is different from the cdf deconvolution problem with loss $L_p$ associated to Ébralidze's inequality (see (14) for the definition).

The paper is organized as follows. In Section 2, some facts about the case without error are recalled and discussed. The upper bounds for Wasserstein deconvolution with supersmooth or ordinary smooth errors are given in Section 3, and Section 4 is about lower bounds. The upper bound is proved in Section 5. Section 6 presents the implementation of the method and some experimental results. In particular, observed rates of convergence are compared with the theoretical bounds for the Wasserstein metrics $W_1$ and $W_2$, and we study as an illustrative example the deconvolution of the uniform measure on the Cantor set.

## 2. On the case without error

We begin by considering the simple case when one observes directly $X_1, \dots, X_n$ with values in $\mathbb{R}$ without error. Let us recall some results for the quantities $W_p(\mu_n, \mu)$, where $\mu_n$ is the empirical measure, given by

$$ \mu_n = \frac{1}{n} \sum_{i=1}^{n} \delta_{X_i}. $$

Let $F$ be the cdf of $X_1$, $F_n$ the cdf of $\mu_n$, and let $F^{-1}$ and $F_n^{-1}$ be their usual cadlag inverses. Recall that, for any $p \ge 1$,

$$ W_p^p(\mu_n, \mu) = \int_0^1 |F_n^{-1}(u) - F^{-1}(u)|^p du, \qquad (2) $$

and if $p = 1$:

$$ W_1(\mu_n, \mu) = \int |F_n^{-1}(u) - F^{-1}(u)| du = \int |F_n(t) - F(t)| dt. $$

The case $p = 1$ has been well understood since the paper by [11]. The random variable $\sqrt{n}W_1(\mu_n, \mu)$ converges in distribution to $\int |B(F(t))| dt$, where $B$ is a standard Brownian bridge, if and only if

$$ \int_0^\infty \sqrt{P(|X| > x)} dx < \infty, \qquad (3) $$
or equivalently if

$$\int_{\mathbb{R}} dx < \infty.$$

More recently, [1] have shown that the rate of $\mathbb{E}W_1(\mu_n, \mu)$ can be characterized by the quantities

$$\int_{4nF(x)(1-F(x))\le 1} F(x)(1-F(x))dx$$

and

$$\int_{4nF(x)(1-F(x))>1} \sqrt{F(x)(1-F(x))} dx.$$

More precisely, the rate $1/\sqrt{n}$ is achieved if and only if (3) is satisfied. When this is not the case, $\mathbb{E}W_1(\mu_n, \mu)$ may decay at an arbitrary slow rate. See the Theorems 3.3 and 3.5 in their paper.

For $p > 1$, the situation is more complicated. Extra conditions are necessary to ensure that $W_p(\mu_n, \mu)$ is of order $1/\sqrt{n}$. If the random variables take their values in a compact interval $[a, b]$ and if the cdf $F$ is continuously differentiable on $[a, b]$ with strictly positive derivative $f$, then $n^{p/2}W_p^p(\mu_n, \mu)$ converges in distribution to $\int_0^1 |B(u)|^p/|f \circ F^{-1}(u)|^p du$ (see Lemma 3.9.23 in [25]). But in general, the rate can be much slower. The convergence in distribution for the case $p = 2$ has been studied in detail by [12]. Under additional conditions on $F$ (see condition (2.7) in [12], which requires in particular that $F$ is twice differentiable), the rate of convergence depends on the behavior of $F^{-1}$ in a neighborhood of 0 and 1. For instance, if

$$F(t) = \left(1 - \frac{1}{t^{\alpha-1}}\right) \mathbf{1}_{\{t>1\}},$$

where $\alpha > 3$, it follows from Theorem 4.7 in [12] that

$$n^{(\alpha-3)/(\alpha-1)} W_2^2(\mu_n, \mu) \qquad (4)$$

converges in distribution. The limiting distribution is explicitly given in [12].

The rates of decay of $\mathbb{E}W_p(\mu_n, \mu)$ and $[\mathbb{E}W_p^p(\mu_n, \mu)]^{1/p}$ have been studied more recently in [1]. They show that these quantities decay at the standard rate $1/\sqrt{n}$ if and only if

$$J_p(\mu) = \int_{\mathbb{R}} \frac{[F(x)(1-F(x))]^{p/2}}{f(x)^{p-1}} dx < \infty,$$

where $f$ is the density of the absolutely continuous component of $\mu$. In particular (see their Theorem 5.6), they show that

$$[\mathbb{E}W_p^p(\mu_n, \mu)]^{1/p} \le \frac{5p}{\sqrt{n+2}} J_p^{1/p}(\mu).$$

However, this approach cannot be applied when the measure $\mu$ and the Lebesgue measure are singular. An alternative approach to obtain the rate of decay of
$\mathbb{E}W_p^p(\mu_n, \mu)$ is to use the following inequality, due to [15] (see also Sections 7.4 and 7.5 in [1]): for any $p \ge 1$,

$$W_p^p(\mu, \nu) \le \kappa_p \int |x|^{p-1} |F_\mu - F_\nu|(x)dx, \quad (5)$$

where $\kappa_p = 2^{p-1}p$. Starting from (5), we get that

$$\begin{align*}
\mathbb{E}W_p^p(\mu_n, \mu) &\le \kappa_p \int |x|^{p-1} \mathbb{E}|F_n(x) - F(x)|dx \\
&\le \kappa_p \int |x|^{p-1} \sqrt{\mathbb{E}|F_n(x) - F(x)|^2}dx \\
&\le \frac{\kappa_p}{\sqrt{n}} \int |x|^{p-1} \sqrt{F(x)(1-F(x))}dx
\end{align*}$$

where $F_n$ is the empirical cdf. Now, it is easy to see that this last integral is finite if and only if

$$\int_0^\infty |x|^{p-1} \sqrt{P(|X| > x)} dx < \infty. \quad (6)$$

It follows that $\mathbb{E}W_p^p(\mu_n, \mu) \le Cn^{-1/2}$ as soon as (6) is satisfied. For instance, taking $p=2$, a tail satisfying $P(|X| > x) = O(\frac{1}{x^4 \log x^{2+\varepsilon}})$ gives the rate $1/\sqrt{n}$. Hence, we obtain the same rate as in (4) for $\alpha = 5$, with a slightly stronger tail condition (due to the fact that we control the expectation), but without additional assumptions on the cdf $F$.

Since we want to estimate singular measures, we shall follow this approach in the sequel.

# 3. Upper bounds for $W_p$ in deconvolution

## 3.1. Construction of the estimator

Let us start with some notations. For $\mu$ a probability measure and $\nu$ another probability measure, with density $g$, we denote by $\mu \star g$ the density of $\mu \star \nu$, given by

$$\mu \star g(x) = \int_{\mathbb{R}} g(x - y)\mu(dy).$$

We further denote by $\mu^*$ (respectively $f^*$) the Fourier transform of the probability measure $\mu$ (respectively of the integrable function $f$), that is

$$\mu^*(x) = \int_{\mathbb{R}} e^{iux} \mu(du) \quad \text{and} \quad f^*(x) = \int_{\mathbb{R}} e^{iux} f(u)du.$$

Finally, let $F$ be the cumulative distribution function of $\mu$.

The estimator $\tilde{\mu}_n$ of the measure $\mu$ is built in two steps:

1. A preliminary estimator of $F$. Let $[p]$ be the least integer greater than or equal to $p$. We first introduce a symmetric nonnegative kernel $k$ such that its Fourier transform $k^*$ is $[p]$ times differentiable with Lipschitz
[p]-th derivative and is supported on [-1, 1]. An example of such a kernel is given by

$$k(x) = C_p \left[ \frac{(2[p/2] + 2) \sin \frac{x}{2[p/2]+2}}{x} \right]^{2[p/2]+2}, \quad (7)$$

where $C_p$ is such that $\int k(x)dx = 1$.

We define now a preliminary estimator $\hat{F}_n$ of F:

$$\hat{F}_n(t) = \frac{1}{nh} \int_{-\infty}^{t} \sum_{k=1}^{n} \tilde{k}_h \left( \frac{u - Y_k}{h} \right) du \quad (8)$$

where

$$\tilde{k}_h(x) = \frac{1}{2\pi} \int \frac{e^{iux} k^*(u)}{\mu^*_e(-u/h)} du.$$

Let us first give some conditions under which these quantities are well defined. Clearly, $\tilde{k}_h(x)$ is well defined as soon as $\mu^*_e$ does not vanish, since in that case it is the Fourier transform of a continuous and compactly supported function (it can be easily checked that $\tilde{k}_h(x)$ is a real function). In the sequel, we shall always assume that $r_e = 1/\mu^*_e$ is at least two times continuously differentiable. In that case, the function $w(u) = \frac{k^*(u)}{\mu^*_e(-u/h)}$ is two times differentiable with bounded and compactly supported derivatives. An integration by parts yields

$$\tilde{k}_h(x) = -\frac{1}{2\pi i x} \int e^{iux} w'(u) du \quad \text{and} \quad \tilde{k}_h(x) = -\frac{1}{2\pi x^2} \int e^{iux} w''(u) du.$$

It follows that $\tilde{k}_h$ is a continuous function such that $\tilde{k}_h(x) = O(1/(1+x^2))$. Hence $\tilde{k}_h$ belongs to $L^1(dx)$ and $\hat{F}_n$ is well defined. Now the inverse Fourier formula gives that $\tilde{k}_h^*(x) = \frac{k^*(u)}{\mu^*_e(u/h)}$. Consequently $\tilde{k}_h^*(0) = 1$, proving that $\int \tilde{k}_h(x)dx = 1$ and that $\lim_{t \to \infty} \hat{F}_n(t) = 1$.

However, this estimator $\hat{F}_n$, based on the standard deconvolution kernel density estimator $\tilde{k}_h$ first introduced by [6], is not a cumulative distribution function since it is not necessarily non-decreasing.

2. **Isotone approximation.** We need to define an estimator $\tilde{F}_n$ of F which is a cumulative distribution function. We choose the estimator $\tilde{F}_n$ as an approximate minimizer over all distribution functions of the quantity $\int_{\mathbb{R}} |x|^{p-1}|\hat{F}_n - G|(x)dx$. More precisely, given $\rho > 0$, let $\tilde{F}_n$ be such that, for every distribution function G,

$$\int |x|^{p-1} |\hat{F}_n - \tilde{F}_n|(x) dx \le \int |x|^{p-1} |\hat{F}_n - G|(x) dx + \rho.$$

Here, $\rho$ may be chosen equal to 0 (best isotone approximation) but the condition $\rho = O(n^{-1/2})$ is the only condition required to get the rates of Section 3.3 below.
The estimator $\tilde{\mu}_n$ is then defined by:

$$ \tilde{\mu}_n \text{ is the probability measure with distribution function } \tilde{F}_n. \quad (9) $$

**Remark 3.1.** The second step is different from that of [10], who choose $\tilde{\mu}_n$ as the (normalized) positive part of $\mu_n$. As we shall see, the isotone approximation allows to get better rates of convergence in the ordinary smooth case. The superiority of the isotone estimator will also be clearly highlighted through the simulations (see Section 6.2). However, this approach works only in the one-dimensional case.

One may argue that the estimator $\tilde{F}_n$ is not explicit, and can be quite difficult to compute, because the minimization is done over an infinite dimensional set. In fact, this is not an issue, because powerful algorithms have been developed to deal with this situation. In Section 6, we shall use the function gpava from the R package isotonic [22] (see Section 6.1 for more details).

## 3.2. First upper bounds for $W_p^p(\tilde{\mu}_n, \mu)$

The control of $W_p^p(\tilde{\mu}_n, \mu)$ is done in three steps:

1. **A bias/random decomposition.** Denoting by $K_h$ the function $h^{-1}k(\cdot/h)$, we have that

$$ W_p^p(\tilde{\mu}_n, \mu) \le 2^{p-1} W_p^p(\mu \star K_h, \mu) + 2^{p-1} W_p^p(\tilde{\mu}_n, \mu \star K_h). \quad (10) $$

The non-random quantity $W_p^p(\mu \star K_h, \mu)$ is the *bias* of the estimator $\tilde{\mu}_n$.

2. **Control of the bias.** Let $V_h$ be a random variable with distribution $K_h$ and independent of $X_1$, in such a way that the distribution of $X_1 + V_h$ is $\mu \star K_h$. By definition of $W_p$, we have

$$ W_p^p(\mu \star K_h, \mu) \le \mathbb{E}|X_1 + V_h - X_1|^p = \mathbb{E}|V_h|^p = h^p \int |x|^p k(x) dx. \quad (11) $$

3. **Control of the random term.** Note that

$$ \mathbb{E}[\hat{F}_n(t)] = \int_{-\infty}^{t} \mu \star K_h(x) dx $$

is the cdf of $\mu \star K_h$. Applying Ébrailidze's inequality (5), we obtain that

$$ W_p^p(\tilde{\mu}_n, \mu \star K_h) \le \kappa_p \int |x|^{p-1} |\tilde{F}_n - \mathbb{E}[\hat{F}_n]|(x) dx. $$

Now, by the triangle inequality and the definition of $\tilde{F}_n$,

$$ \begin{aligned} W_p^p(\tilde{\mu}_n, \mu \star K_h) & \le \kappa_p \left( \int |x|^{p-1} |\tilde{F}_n - \hat{F}_n|(x) dx + \int |x|^{p-1} |\hat{F}_n - \mathbb{E}[\hat{F}_n]|(x) dx \right) \\ & \le \rho + 2\kappa_p \int |x|^{p-1} |\hat{F}_n - \mathbb{E}[\hat{F}_n]|(x) dx. \end{aligned} \quad (12) $$
From (10), (11) and (12), to get explicit rates of convergence for $\mathbb{E}[W_p^p(\tilde{\mu}_n, \mu)]$, it remains to control the term

$$ \mathbb{E} \left( \int |x|^{p-1} |\hat{F}_n - \mathbb{E}[\hat{F}_n]|(x) dx \right). $$

**Remark 3.2.** Another main difference between the present paper and [10] is the use of Ébralidze's inequality (5) to control the random term. In [10] the term $W_p^p(\tilde{\mu}_n, \mu \star K_h)$ (for another choice of $\tilde{\mu}_n$) is bounded by a term involving the variation norm between $\tilde{\mu}_n$ and $\mu$. In our case, this upper bound would give a worse rate of convergence.

Note that Inequality (5) is used here to control the random term only. A possible alternative approach is to use (5) directly, as in the case without error (see Section 2). This would give the upper bound

$$ W_p^p(\tilde{\mu}_n, \mu) \le \kappa_p \int |x|^{p-1} |\tilde{F}_n - F|(x) dx. \quad (13) $$

In that case, the bias term would be

$$ \int |x|^{p-1} |\mathbb{E}[\hat{F}_n] - F|(x) dx. $$

However, without extra regularity assumptions on $\mu$, this would give a bias term of order $h$, and then the same rate of convergence as in the case $p = 1$, that is $n^{1/(2\beta+1)}$ under the assumptions of Theorem 3.1 (Item 2) of the next section. But this rate is always too slow for $p > 1$, see again Theorem 3.1. Moreover, there is no hope to obtain a better rate from (13) because $n^{1/(2\beta+1)}$ is also the minimax rate to estimate $F$ with the loss function

$$ L_p(G, F) = \int |x|^{p-1} |G - F|(x) dx, \quad (14) $$

This last assertion comes from the lower bound stated in Theorem 4.2 of Section 4.

### 3.3. Main results

Let $r_\varepsilon = 1/\mu_\varepsilon^*$, and let $r_\varepsilon^{(\ell)}$ be the $\ell$-th derivative of $r_\varepsilon$. Let $m_0$ denote the least integer strictly greater than $p + \frac{1}{2}$, and $m_1$ be the least integer strictly greater than $p - \frac{1}{2}$.

Our first result is a general proposition which gives an upper bound for $\mathbb{E}W_p^p(\tilde{\mu}_n, \mu)$ involving a tail condition on $Y$ and the regularity of $r_\varepsilon$.

**Proposition 3.1.** Let $\rho \le n^{-1/2}$, and let $\tilde{\mu}_n$ be the estimator defined in (9). Assume that $r_\varepsilon$ is $m_0$ times continuously differentiable. For any $h \le 1$, we have

$$ \mathbb{E}W_p^p(\tilde{\mu}_n, \mu) \le \frac{1}{\sqrt{n}} + h^p 2^{p-1} \int |x|^p k(x) dx + \frac{C}{\sqrt{n}} (A_1 + A_2 + A_3 + A_4) $$
where

$$A_1 = \left( \sup_{t \in [-2,2]} \sum_{\ell=0}^{1} |r_{\varepsilon}^{(\ell)}(t)| \right) \int_{0}^{\infty} |x|^{p-1} \sqrt{P(|Y| \ge x)} dx$$

$$A_2 = \sup_{t \in [-2,2]} \sum_{\ell=0}^{m_0} |r_{\varepsilon}^{(\ell)}(t)|$$

$$A_3 = \left[ \mathbb{E}|Y|^{2p-\frac{1}{2}} \int_{-1/h}^{1/h} \frac{|r_{\varepsilon}(x)|^2}{|x|^2} \mathbf{1}_{[-1,1]^c}(x) dx \right]^{1/2}$$

$$A_4 = \left[ \sum_{\ell=0}^{m_1} \int_{-1/h}^{1/h} \frac{|r_{\varepsilon}^{(\ell)}(x)|^2}{|x|^2} \mathbf{1}_{[-1,1]^c}(x) dx \right]^{1/2}.$$

For the sake of readability, the proof of Proposition 3.1 is postponed to Section 5.

We are now in a position to give the rates of convergence for the Wasserstein deconvolution, for a class of supersmooth error distributions, and for a class of ordinary smooth error distributions.

**Theorem 3.1.** Let $\rho \le n^{-1/2}$, and let $\tilde{\mu}_n$ be the estimator defined in (9). Assume that

$$\int_{0}^{\infty} |x|^{p-1} \sqrt{P(|Y| \ge x)} dx < \infty \quad \text{and} \quad \sup_{t \in [-2,2]} |r_{\varepsilon}^{(m_0)}(t)| < \infty. \quad (15)$$

1. Assume that there exist $\beta > 0$, $\tilde{\beta} \ge 0$, $\gamma > 0$ and $c > 0$, such that for every $\ell \in \{0, 1, \dots, m_1\}$ and every $t \in \mathbb{R}$,

$$|r_{\varepsilon}^{(\ell)}(t)| \le c(1 + |t|)^{\tilde{\beta}} \exp(|t|^{\beta}/\gamma). \quad (16)$$

Then, taking $h = (4/(\gamma \log n))^{1/\beta}$, there exists a positive constant $C$ such that

$$\mathbb{E}W_p^p(\tilde{\mu}_n, \mu) \le C(\log n)^{-p/\beta}.$$

2. Assume that there exist $\beta > 0$ and $c > 0$, such that for every $\ell \in \{0, 1, \dots, m_1\}$ and every $t \in \mathbb{R}$,

$$|r_{\varepsilon}^{(\ell)}(t)| \le c(1 + |t|)^{\beta}. \quad (17)$$

Then, taking $h = n^{-1/(2p+(2\beta-1))}$, there exists a positive constant $C$ such that

$$\mathbb{E}W_p^p(\tilde{\mu}_n, \mu) \le C\psi_n, \quad (18)$$

where

$$\psi_n = \begin{cases} n^{-\frac{p}{2p+(2\beta-1)}} & \text{if } \beta > \frac{1}{2} \\ \sqrt{\frac{\log n}{n}} & \text{if } \beta = \frac{1}{2} \\ \frac{1}{\sqrt{n}} & \text{if } \beta < \frac{1}{2}. \end{cases}$$
This result requires several comments.

**Remark 3.3.** In the ordinary smooth case, when $\beta < 1/2$, any bandwidth $h = O(n^{-1/2p})$ leads to the rate $n^{-1/2}$. The fact that there are three different situations according as $\beta > 1/2$, $\beta = 1/2$ or $\beta < 1/2$ has already been pointed out in Theorem 3.2 of [21] and in Theorem 2.1 of [9] for the estimation of the cdf $F$. Note that the estimator $\hat{F}_n$ of [21] is exactly the estimator defined in (8) (with possibly a slightly different kernel). Hence it is not always non-decreasing and cannot be used directly to estimate $\mu$ with respect to Wasserstein metrics.

For instance, for a Laplace error distribution, the estimator $\hat{F}_n$ of [21] is such that

$$ \left( \mathbb{E} \left[ \int |\hat{F}_n(t) - F(t)|^2 dt \right] \right)^{1/2} \leq C n^{-1/8}, $$

while the rate of convergence of our estimator for $W_1$ is

$$ \mathbb{E}W_1(\tilde{\mu}_n, \mu) = \mathbb{E} \left[ \int |\tilde{F}_n(t) - F(t)| dt \right] \leq C n^{-1/5}. $$

In both cases, there are no assumptions on $\mu$, except moment assumptions; in particular, $\mu$ needs not be absolutely continuous with respect to the Lebesgue measure. It is then a different context than that considered by [9] for the point-wise estimation of $F(x_0)$. In this paper, the authors always assume that $\mu$ is absolutely continuous with respect to the Lebesgue measure, with a density $f$ belonging to a Sobolev space of order $\alpha > -1/2$.

Note that the two rates described in this remark are minimax (see Section 4 for our estimator).

**Remark 3.4.** Since the function $H_Y(x) = P(|Y| \geq x)$ is non-increasing, the tail condition

$$ \int_{0}^{\infty} |x|^{p-1} \sqrt{P(|Y| \geq x)} dx < \infty \qquad (19) $$

in Assumption (15) implies that $H_Y(x) = O(1/|x|^{2p})$. Hence $|Y|$ has a weak moment of order $2p$, which implies a strong moment of order $q$ for any $q < 2p$. Note that (19) is the same as the tail condition (6) obtained in Section 2 to get the rate $\mathbb{E}W_p^p(\mu_n, \mu) \leq Cn^{-1/2}$ in the case without noise. Recall that, in the case without noise when $p = 1$, this condition is necessary and sufficient for the weak convergence of $\sqrt{n}W_1(\mu_n, \mu)$. Note also that

$$ (19) \text{ holds iff (6) holds and } \int_{0}^{\infty} |x|^{p-1} \sqrt{P(|\varepsilon| \geq x)} dx < \infty. $$

The “if” part follows easily from the simple inequality $P(|X+\varepsilon| > x) \leq P(|X| > x/2) + P(|\varepsilon| > x/2)$. To prove the “only if” part, note that, since $X$ and $\varepsilon$ are independent, (19) can be written

$$ \int_{0}^{\infty} \int_{0}^{\infty} |x|^{p-1} \sqrt{P(|X+y| \geq x)} dx \, d\mu_{\varepsilon}(dy) < \infty. \qquad (20) $$
But this implies that

$$ \int_0^\infty |x|^{p-1} \sqrt{P(|X+y| \ge x)} dx < \infty \quad (21) $$

for $\mu_\varepsilon$ almost every $y$. Now if (21) holds for one $y$, then it holds for every $y$, proving that (6) holds (and the same is true for $\varepsilon$ by interchanging $X$ and $\varepsilon$ in (20)). As we have seen, the tail condition on $\varepsilon$ implies that $|\varepsilon|$ has a moment of order $k$ for any integer $k$ strictly less than $2p$, hence $\mu_\varepsilon^*$ is at least $k$ times continuously differentiable.

**Remark 3.5.** The rate $\mathbb{E}W_p^p(\tilde{\mu}_n, \mu) \le C(\log n)^{-p/\beta}$ in the supersmooth case has already been given in Theorem 4 of [10] and is valid in any dimension. However the condition on the regularity of $r_\varepsilon$ is more restrictive in the paper by [10], since it is assumed there that Condition (16) is true for $\ell \in \{0, 1, \dots, \lfloor p \rfloor + 1\}$. Note that this rate is minimax, as stated in Theorem 2 of [10].

**Remark 3.6.** Applying Proposition 1 in [10], if Condition (17) is true for $\ell \in \{0, 1, \dots, \lfloor p \rfloor + 1\}$, one can build an explicit estimator $\tilde{\mu}_n$ such that $\mathbb{E}W_p^p(\tilde{\mu}_n, \mu) \le Cn^{-p/(2p+2\beta+1)}$, which is worse than (18). The estimator $\tilde{\mu}_n$ is the “naive” estimator defined in Section 6.1. However, the procedure given in [10] works also when the observations $Y_i$ are $\mathbb{R}^d$-valued, whereas the estimator $\tilde{\mu}_n$ defined in (9) is well defined for $d = 1$ only. Hence, a reasonable question is: can we improve on Proposition 1 of [10] in any dimension?

*Proof.* We first prove Item 1. From Proposition 3.1 and Assumptions (15) and (16), we obtain the upper bound

$$ \mathbb{E}W_p^p(\tilde{\mu}_n, \mu) \le C \left( h^p + \frac{1}{\sqrt{n}} \frac{1}{h^{\beta}} e^{1/h^{\beta}\gamma} \right). $$

Taking $h = (4/(\gamma \log(n)))^{1/\beta}$ gives the result.

We now prove Item 2. From Proposition 3.1 and Assumptions (15) and (17), we obtain

$$ \mathbb{E}W_p^p(\tilde{\mu}_n, \mu) \le \begin{cases} C \left(h^p + \frac{1}{\sqrt{n}} \frac{1}{h^{\beta-1/2}}\right) & \text{if } \beta > \frac{1}{2} \\ C \left(h^p + \frac{1}{\sqrt{n}} \sqrt{\log(\frac{1}{h})}\right) & \text{if } \beta = \frac{1}{2} \\ C \left(h^p + \frac{1}{\sqrt{n}}\right) & \text{if } \beta < \frac{1}{2}. \end{cases} $$

Taking $h = n^{-\frac{1}{2p+(2\beta-1)+}}$ gives the result. $\square$

## 4. Lower bound

For some $M > 0$ and $q \ge 1$, we denote by $\mathcal{D}(M, q)$ the set of measures $\mu$ on $\mathbb{R}$ such that $\int |x|^q d\mu(x) \le M$.
**Theorem 4.1.** Let $M > 0$ and $q \ge 1$. Assume that there exist $\beta > 0$ and $c > 0$, such that for every $\ell \in \{0, 1, 2\}$ and every $t \in \mathbb{R}$,

$$|\mu_{\varepsilon}^{*(\ell)}(t)| \le c(1 + |t|)^{-\beta}. \qquad (22)$$

Then, there exists a constant $C > 0$ such that, for any estimator $\hat{\mu}$,

$$\liminf_{n \to \infty} n^{\frac{p}{2\beta+1}} \sup_{\mu \in \mathcal{D}(M,q)} \mathbb{E} W_p^p(\hat{\mu}, \mu) > C.$$

**Remark 4.1.** For $W_1$, this lower bound matches the upper bound given in Theorem 3.1 for $\beta \ge 1/2$. For $W_p$ ($p > 1$), we conjecture that the upper bounds given by Theorem 3.1 are appropriate under the assumed tail conditions. Getting better rates of convergence for $W_p$ ($p > 1$) is an open question. From Section 2, it seems reasonable to think that better rates can be obtained when $\mu$ has an absolutely continuous component with respect to the Lebesgue measure which is strictly positive on the support of $\mu$ (and also that this should be a necessary condition condition to reach the lower bound when $\beta > 1/2$).

We also give a lower bound for the cdf deconvolution problem with loss $L_p$ defined in (14).

**Theorem 4.2.** Let $M > 0$ and $q \ge 1$. Assume that there exist $\beta > 0$ and $c > 0$, such that (22) is satisfied for every $\ell \in \{0, 1, 2\}$ and every $t \in \mathbb{R}$. Then, there exists a constant $C > 0$ such that, for any estimator $\hat{F}$ of $F$:

$$\liminf_{n \to \infty} n^{\frac{1}{2\beta+1}} \sup_{\mu \in \mathcal{D}(M,q)} \mathbb{E} \int_{\mathbb{R}} |x|^{p-1} |\hat{F}(x) - F(x)| dx > C.$$

We give below the proof of Theorem 4.1 for the Wasserstein metric. The proof of Theorem 4.2 is similar, it can be easily adapted from the proofs of Theorem 4.1 and of Theorem 3 in [10].

*Proof.* Let $M > 0$ and $q \ge 1$. The proof is similar to the proof of Theorem 3 in [10] and thus we only give here a sketch of the proof. We first define a finite family in $\mathcal{D}(M, q)$ using the densities

$$f_{0,r}(t) := C_r(1 + t^2)^{-r} \qquad (23)$$

with some $r > (1 + q)/2$. Next, let $b_n$ be the sequence

$$b_n := \left[ n^{\frac{1}{2\beta+1}} \right] \vee 1, \qquad (24)$$

where $[\cdot]$ is the integer part. For any $\theta \in \{0, 1\}^{b_n}$, let

$$f_{\theta}(t) = f_{0,r}(t) + C \sum_{s=1}^{b_n} \theta_s H(b_n(t - t_{s,n})), \quad t \in \mathbb{R}, \qquad (25)$$

where $C$ is a positive constant and $t_{s,n} = (s - 1)/b_n$. The function $H$ is a bounded function whose integral on the line is 0. Moreover, we may choose a function $H$ such that (see for instance [17] or [18]):
$$ (A1) \quad \int_{-\infty}^{+\infty} H(t) dt = 0 \text{ and } \int_{0}^{1} |H^{(-1)}(t)| dt > 0, $$

$$ (A2) \quad |H(t)| \le c(1+t^2)^{-r_0} \text{ where } r_0 > \max(3/2, (1+q)/2), $$

$$ (A3) \quad H^*(z) = 0 \text{ outside } [1, 2], $$

where $H^{(-1)}(t) := \int_{-\infty}^{t} H(u) du$ is a primitive of $H$. Note that by replacing $H$ by $H/C$ in the following, we finally can take $C = 1$ in (25). Let $\mu_{\theta}$ be the measure of density $f_{\theta}$ with respect to the Lebesgue measure. Then we can find some $M$ large enough such that for all $\theta \in \{0, 1\}^{b_n}$, $\mu_{\theta} \in \mathcal{D}(M, q)$. Moreover, under these assumptions the first two derivatives of $H^*$ are continuous and bounded.

For $\theta \in \{0, 1\}^{b_n}$ and $s \in \{1, \dots, b_n\}$, let us define the probability measures $\mu_{\theta,s,0}$ and $\mu_{\theta,s,1}$ with densities

$$ f_{\theta,s,0} := f(\theta_1, \dots, \theta_{s-1}, 0, \theta_{s+1}, \dots, \theta_{b_n}) \quad \text{and} \quad f_{\theta,s,1} := f(\theta_1, \dots, \theta_{s-1}, 1, \theta_{s+1}, \dots, \theta_{b_n}). $$

We also consider the densities $h_{\theta,s,u} = f_{\theta,s,u} \star \mu_{\varepsilon}$ for $u = 0$ or $1$. Since $W_1$ is dominated by $W_p$, and using Jensen's inequality, it follows that

$$ \sup_{\mu \in \mathcal{D}(M,q)} \mathbb{E}_{(\mu \star \mu_{\varepsilon})^{\otimes n}} W_p^p (\mu, \tilde{\mu}_n) \ge \sup_{\mu \in \mathcal{D}(M,q)} \mathbb{E}_{(\mu \star \mu_{\varepsilon})^{\otimes n}} W_1^p (\mu, \tilde{\mu}_n) \\ \ge \left( \sup_{\mu \in \mathcal{D}(M,q)} \mathbb{E}_{(\mu \star \mu_{\varepsilon})^{\otimes n}} W_1 (\mu, \tilde{\mu}_n) \right)^p. \quad (26) $$

Using a standard randomization argument (see for the instance the proof of Theorem 3 in [10] for the multivariate case), it can be shown that there exists a constant $C > 0$ such that

$$ \sup_{\mu \in \mathcal{D}(M,q)} \mathbb{E}_{(\mu \star \mu_{\varepsilon})^{\otimes n}} W_1 (\mu, \tilde{\mu}_n) \ge \frac{C}{b_n} \int_0^1 |H^{(-1)}(u)| du \quad (27) $$

as soon as there exists a constant $c > 0$ such that, for any $\theta \in \{0, 1\}^{b_n}$,

$$ \chi^2(h_{\theta,s,0}, h_{\theta,s,1}) \le \frac{c}{n}, \quad (28) $$

where the $\chi^2$ distance between two densities $h_1$ and $h_2$ on $\mathbb{R}$ is defined by

$$ \chi^2(h_1, h_2) = \int \frac{\{(h_1(x) - h_2(x))^2\}}{h_1(x)} dx. $$

If (28) is satisfied, we take $b_n$ as in (24) and the theorem is thus proved according to (26), (27) and (A1).

It remains to prove (28). Using (A2), we can find a constant $C > 0$ such that for any $t \in \mathbb{R}$ and any $s \in \{1, \dots, b_n\}$,

$$ \chi^2(h_{\theta,s,0}, h_{\theta,s,1}) \le C b_n^{-1} \int \frac{\{\int H(v-y) \mu_{\varepsilon}(dy/b_n)\}^2}{f_{0,r} \star \mu_{\varepsilon}(v/b_n)} dv. \quad (29) $$

The right side of (29) is typically the kind of $\chi^2$ divergence that is upper bounded in the proofs of Theorems 4 and 5 in [17] for computing pointwise rates of
convergence: under Assumption (22), it gives that there exists a constant C such that

$$ \int \frac{\left\{ \int H(v-y) \mu_{\varepsilon}(dy/b_n) \right\}^2}{f_{0,r} \star \mu_{\varepsilon}(v/b_n)} dv \le C b_n^{-2\beta} $$

and (28) is proved. □

## 5. Proof of Proposition 3.1

Throughout, C will denote a positive constant depending on p which may change from line to line.

We start from the basic inequality (10). Inequality (11) yields the bias term

$$ h^p 2^{p-1} \int |x|^p k(x) dx, $$

and it remains to control the term $\mathbb{E}W_p^p(\tilde{\mu}_n, \mu \star K_h)$.

By (12), we have

$$ \begin{aligned} \mathbb{E}W_p^p(\tilde{\mu}_n, \mu \star K_h) &\le C \int |x|^{p-1} \mathbb{E}|\hat{F}_n - \mathbb{E}[\hat{F}_n]|(x) dx + \rho \\ &\le C \int |x|^{p-1} \sqrt{\mathrm{Var}(\hat{F}_n)(x)} dx + \rho. \end{aligned} \quad (30) $$

Now, let $\phi$ denote a symmetric function, $[p]+1$ times continuously differentiable, equal to 1 on the interval $[-1, 1]$ and to 0 outside $[-2, 2]$. Our preliminary estimator $\hat{F}_n$ may be written

$$ \begin{aligned} \hat{F}_n(t) &= \frac{1}{nh} \int_{-\infty}^{t} \sum_{k=1}^{n} \tilde{k}_h \left( \frac{u - Y_k}{h} \right) du \\ &= \frac{1}{n} \sum_{k=1}^{n} G_{1,h} \left( \frac{t - Y_k}{h} \right) + \frac{1}{n} \sum_{k=1}^{n} G_{2,h} \left( \frac{t - Y_k}{h} \right) \\ &:= \hat{F}_{1,n} + \hat{F}_{2,n}, \end{aligned} $$

where

$$ G_{1,h}(x) = \int_{-\infty}^{x} \tilde{k}_{1,h}(u) du \quad \text{and} \quad G_{2,h}(x) = \int_{-\infty}^{x} \tilde{k}_{2,h}(u) du. $$

Here,

$$ \tilde{k}_{1,h}(u) = \frac{1}{2\pi} \int \frac{e^{itu} k^*(t) \phi(t/h)}{\mu_{\varepsilon}^*(-t/h)} dt, \quad \tilde{k}_{2,h}(u) = \frac{1}{2\pi} \int \frac{e^{itu} k^*(t)(1 - \phi(t/h))}{\mu_{\varepsilon}^*(-t/h)} dt. $$

From (30), we infer that

$$ \mathbb{E}W_p^p(\tilde{\mu}_n, \mu \star K_h) \le C(I + J) + \rho, \quad (31) $$
where

$$I = \int |x|^{p-1} \sqrt{\mathrm{Var}(\hat{F}_{1,n})(x)} dx \quad \text{and} \quad J = \int |x|^{p-1} \sqrt{\mathrm{Var}(\hat{F}_{2,n})(x)} dx.$$

To prove Proposition 3.1, we shall give some upper bounds for the terms I and J.

**Control of I** We first split the integral into two parts:

$$I = \int_{-\infty}^{0} |x|^{p-1} \sqrt{\mathrm{Var}(\hat{F}_{1,n})(x)} dx + \int_{0}^{\infty} |x|^{p-1} \sqrt{\mathrm{Var}(\hat{F}_{1,n})(x)} dx := I^{-} + I^{+}.$$

Now,

$$\begin{aligned} I^{-} &= \int_{-\infty}^{0} |x|^{p-1} \sqrt{\mathrm{Var}(\hat{F}_{1,n})(x)} dx \\ &\le \frac{C}{\sqrt{n}} \int_{-\infty}^{0} |x|^{p-1} \sqrt{\mathrm{E}\left[G_{1,h}\left(\frac{x-Y}{h}\right)\right]^2} dx \\ &\le \frac{C}{\sqrt{n}} \int_{-\infty}^{0} |x|^{p-1} \sqrt{\mathrm{E}\left[\int \tilde{k}_{1,h}(u) \mathbf{1}_{\{u \le \frac{x-Y}{h}\}} du\right]^2} dx. \end{aligned}$$

Then, letting $z = uh$ and applying Cauchy-Schwarz's inequality we obtain, for any $a \in ]0, 1[$,

$$\begin{aligned} I^{-} &\le \frac{C}{\sqrt{n}} \int_{-\infty}^{0} |x|^{p-1} \sqrt{\mathrm{E}\left[\int \frac{\tilde{k}_{1,h}(z/h)}{h} \mathbf{1}_{\{Y+z \le x\}} dz\right]^2} dx \\ &\le \frac{C}{\sqrt{n}} \int_{-\infty}^{0} |x|^{p-1} \sqrt{\mathrm{E}\left[\int (1+|z|^{1+a}) \left(\frac{\tilde{k}_{1,h}(z/h)}{h}\right)^2 \mathbf{1}_{\{Y+z \le x\}} dz\right]} dx. \end{aligned}$$

Noticing that $\mathbf{1}_{\{Y+z \le x\}} \le \mathbf{1}_{\{Y \le \frac{x}{2}\}} + \mathbf{1}_{\{z \le \frac{x}{2}\}}$, we obtain that $I^{-} \le I_{1}^{-} + I_{2}^{-}$, where

$$\begin{aligned} I_{1}^{-} &= \frac{C}{\sqrt{n}} \int_{-\infty}^{0} |x|^{p-1} \sqrt{\mathrm{E}\left[\int (1 + |z|^{1+a}) \left(\frac{\tilde{k}_{1,h}(z/h)}{h}\right)^2 \mathbf{1}_{\{Y \le \frac{x}{2}\}} dz\right]} dx \\ I_{2}^{-} &= \frac{C}{\sqrt{n}} \int_{-\infty}^{0} |x|^{p-1} \sqrt{\int (1 + |z|^{1+a}) \left(\frac{\tilde{k}_{1,h}(z/h)}{h}\right)^2 \mathbf{1}_{\{z \le \frac{x}{2}\}} dz} dx. \end{aligned}$$

To control the term $I_1^-$, note that

$$I_1^- \le \frac{C}{\sqrt{n}} \sqrt{\int (1 + |z|^2) \left(\frac{\tilde{k}_{1,h}(z/h)}{h}\right)^2 dz} \int_{-\infty}^{0} |x|^{p-1} \sqrt{P(Y \le \frac{x}{2})} dx.$$

Here we shall use the following lemma.
**Lemma 5.1.** For any nonnegative integer $k$ and any $h \le 1$ we have

$$ \int |z|^{2k} \left( \frac{\tilde{k}_{1,h}(z/h)}{h} \right)^2 dz \le C \left( \sup_{t \in [-2,2]} \sum_{\ell=0}^{k} |r_{\varepsilon}^{(\ell)}(t)| \right)^2 . $$

*Proof.* By definition of $\tilde{k}_{1,h}$,

$$ \int |z|^{2k} \left( \frac{\tilde{k}_{1,h}(z/h)}{h} \right)^2 dz \le \frac{1}{4\pi^2} \int |z|^{2k} \left| \int \frac{e^{iuz} k^*(uh)\phi(u)}{\mu_{\varepsilon}^*(-u)} du \right|^2 dz. $$

Now, by Parseval-Plancherel's identity,

$$ \int |z|^{2k} \left| \int \frac{e^{iuz} k^*(uh)\phi(u)}{\mu_{\varepsilon}^*(-u)} du \right|^2 dz = 2\pi \int \left| \left( \frac{k^*(th)\phi(t)}{\mu_{\varepsilon}^*(-t)} \right)^{(k)} \right|^2 dt. $$

It can be checked that, for $h \le 1$,

$$ \left| \left( \frac{k^*(th)\phi(t)}{\mu_{\varepsilon}^*(-t)} \right)^{(k)} \right| \le C \sum_{\ell=0}^{k} |r_{\varepsilon}^{(\ell)}(t)| \mathbf{1}_{[-2,2]}(t), $$

which concludes the proof of the Lemma.

<figure><img src="image_2.png" /></figure>

Applying Lemma 5.1 with $k=1$, we obtain that

$$ I_1^- \le \frac{C}{\sqrt{n}} \left( \sup_{t \in [-2,2]} \sum_{\ell=0}^{1} |r_{\varepsilon}^{(\ell)}(t)| \right) \int_{-\infty}^{0} |x|^{p-1} \sqrt{P\left(Y \le \frac{x}{2}\right)} dx. \quad (32) $$

We now control the term $I_2^-$. Let $b \in ]0,1[$. Applying Cauchy-Schwarz's inequality

$$ I_2^- \le \frac{C}{\sqrt{n}} \sqrt{ \int_{-\infty}^{0} |x|^{2p-2} (1+|x|^{1+b}) \int_{-\infty}^{\frac{x}{2}} (1+|z|^{1+a}) \left( \frac{\tilde{k}_{1,h}(z/h)}{h} \right)^2 dz dx. } $$

Consequently, by Fubini's Theorem

$$ \begin{aligned} I_2^- & \le \frac{C}{\sqrt{n}} \sqrt{ \int_{-\infty}^{0} (1 + |z|^{1+a}) \left( \frac{\tilde{k}_{1,h}(z/h)}{h} \right)^2 \int_{2z}^{0} |x|^{2p-2} (1 + |x|^{1+b}) dx dz } \\ & \le \frac{C}{\sqrt{n}} \sqrt{ \int (1 + |z|^{2p+1+a+b}) \left( \frac{\tilde{k}_{1,h}(z/h)}{h} \right)^2 dz } \end{aligned} $$

Let $m_0$ be the least integer strictly greater than $p+1/2$. Taking $a$ and $b$ close enough to 0, it follows that

$$ I_2^- \le \frac{C}{\sqrt{n}} \sqrt{ \int (1 + |z|^{2m_0}) \left( \frac{\tilde{k}_{1,h}(z/h)}{h} \right)^2 dz } $$
Applying Lemma 5.1 with $k = m_0$, it follows that

$$I_2^{-} \le \frac{C}{\sqrt{n}} \left( \sup_{t \in [-2,2]} \sum_{\ell=0}^{m_0} |r_{\varepsilon}^{(\ell)}(t)| \right). \quad (33)$$

In the same way, we have

$$ \begin{aligned} I^+ &= \int_0^\infty |x|^{p-1} \sqrt{\mathrm{Var}(1 - \hat{F}_{1,n})(x)} dx \\ &\le \frac{C}{\sqrt{n}} \int_0^\infty |x|^{p-1} \sqrt{\mathbb{E}\left[1 - G_{1,h}\left(\frac{x-Y}{h}\right)\right]^2} dx \\ &\le \frac{C}{\sqrt{n}} \int_0^\infty |x|^{p-1} \sqrt{\mathbb{E}\left[\int \tilde{k}_{1,h}(u) \mathbf{1}_{\{u \ge \frac{x-Y}{h}\}} du\right]^2} dx. \end{aligned} $$

Using the same arguments as for $I^-$, we obtain,

$$ \begin{aligned} I^+ \le \frac{C}{\sqrt{n}} \left( \sup_{t \in [-2,2]} \sum_{\ell=0}^{1} |r_{\varepsilon}^{(\ell)}(t)| \right) \int_0^\infty |x|^{p-1} \sqrt{P\left(Y \ge \frac{x}{2}\right)} dx \\ + \frac{C}{\sqrt{n}} \left( \sup_{t \in [-2,2]} \sum_{\ell=0}^{m_0} |r_{\varepsilon}^{(\ell)}(t)| \right). \end{aligned} \quad (34) $$

Consequently, gathering (32), (33) and (34) we obtain that

$$ \begin{aligned} I \le \frac{C}{\sqrt{n}} \left( \sup_{t \in [-2,2]} \sum_{\ell=0}^{1} |r_{\varepsilon}^{(\ell)}(t)| \right) \int_0^\infty |x|^{p-1} \sqrt{P(|Y| \ge x)} dx \\ + \frac{C}{\sqrt{n}} \left( \sup_{t \in [-2,2]} \sum_{\ell=0}^{m_0} |r_{\varepsilon}^{(\ell)}(t)| \right). \end{aligned} \quad (35) $$

**Control of J** Let $a \in ]0, 1/2[$. By definition of the term $J$, and applying Cauchy-Schwarz's inequality,

$$ \begin{aligned} J &\le \frac{C}{\sqrt{n}} \int |x|^{p-1} \sqrt{\mathbb{E}\left[G_{2,h}\left(\frac{x-Y}{h}\right)\right]^2} dx \\ &\le \frac{C}{\sqrt{n}} \sqrt{\int |x|^{2p-2} (1 + |x|^{1+a}) \mathbb{E}\left[G_{2,h}\left(\frac{x-Y}{h}\right)\right]^2} dx. \end{aligned} $$

Let us write

$$ G_{2,h}(x) = \lim_{T \to -\infty} \int_T^x \tilde{k}_{2,h}(u) du = \lim_{T \to -\infty} \int \mathbf{1}_{[T,x]}(u) \tilde{k}_{2,h}(u) du. $$

Using Parseval-Plancherel's identity, we get

$$ G_{2,h}(x) = \lim_{T \to -\infty} \frac{1}{2\pi} \int \overline{\mathbf{1}_{[T,x]}(u)} \tilde{k}_{2,h}^*(u) du $$
$$ = -\frac{1}{2\pi i} \left[ \int \frac{e^{-itx} k^*(t)(1 - \phi(t/h))}{t \mu_{\varepsilon}^*(t/h)} dt - \lim_{T \to -\infty} \int \frac{e^{-itT} k^*(t)(1 - \phi(t/h))}{t \mu_{\varepsilon}^*(t/h)} dt \right] $$

Since the function $t \mapsto \frac{k^*(t)(1-\phi(t/h))}{t\mu_{\varepsilon}^*(t/h)}$ is integrable, the Riemann-Lebesgue Lemma ensures that

$$ \lim_{T \to -\infty} \int \frac{e^{-itT} k^*(t)(1 - \phi(t/h))}{t \mu_{\varepsilon}^*(t/h)} dt = 0, $$

so that

$$ G_{2,h}(x) = -\frac{1}{2\pi i} \int \frac{e^{-itx} k^*(t)(1 - \phi(t/h))}{t \mu_{\varepsilon}^*(t/h)} dt. $$

Consequently,

$$ J \le \frac{C}{\sqrt{n}} \sqrt{\mathbb{E} \left( \left( \int (1 + |x|^{2p-1+a}) \left( -\frac{1}{2\pi i} \int \frac{e^{-it \frac{x-Y}{h}} k^*(t)(1-\phi(t/h))}{t \mu_{\varepsilon}^*(t/h)} dt \right)^2 dx \right)^{1/2} \right)}. $$

Setting $u = t/h$ and using the fact that $|x|^q \le 2^{q-1}|x - Y|^q + 2^{q-1}|Y|^q$ for any $q \ge 1$, we obtain that

$$ J \le \frac{C}{\sqrt{n}} \left[ \mathbb{E} \left( \int |x - Y|^{2p-1+a} \left( -\frac{1}{2\pi i} \int \frac{e^{-iu(x-Y)} k^*(uh)(1-\phi(u))}{u\mu_{\varepsilon}^*(u)} du \right)^2 dx \right) + \mathbb{E} \left( (1 + |Y|^{2p-\frac{1}{2}}) \int \left( -\frac{1}{2\pi i} \int \frac{e^{-iu(x-Y)} k^*(uh)(1-\phi(u))}{u\mu_{\varepsilon}^*(u)} du \right)^2 dx \right) \right]^{1/2}. $$

Thus,

$$ J \le \frac{C}{\sqrt{n}} \left[ \int (1 + |x|^{2p-1+a}) \left( -\frac{1}{2\pi i} \int \frac{e^{-iu x} k^*(uh)(1-\phi(u))}{u\mu_{\varepsilon}^*(u)} du \right)^2 dx + \mathbb{E} |Y|^{2p-\frac{1}{2}} \int \left( -\frac{1}{2\pi i} \int \frac{e^{-iu x} k^*(uh)(1-\phi(u))}{u\mu_{\varepsilon}^*(u)} du \right)^2 dx \right]^{1/2}. $$

Let $m_1$ be the least integer strictly greater than $p - \frac{1}{2}$. Taking $a$ close enough to zero, it follows that

$$ J \le \frac{C}{\sqrt{n}} \left[ \int (1 + |x|^{2m_1}) \left| \frac{1}{2\pi} \int \frac{e^{-iu x} k^*(uh)(1 - \phi(u))}{u\mu_{\varepsilon}^*(u)} du \right|^2 dx + \mathbb{E} |Y|^{2p-\frac{1}{2}} \int \left| \frac{1}{2\pi} \int \frac{e^{-iu x} k^*(uh)(1 - \phi(u))}{u\mu_{\varepsilon}^*(u)} du \right|^2 dx \right]^{1/2}. $$
By Parseval-Plancherel's identity,

$$ \int \left| \int \frac{e^{-iux} k^*(uh)(1-\phi(u))}{u\mu_{\varepsilon}^*(u)} du \right|^2 dx = 2\pi \int \left| \frac{k^*(th)(1-\phi(t))}{t\mu_{\varepsilon}^*(-t)} \right|^2 dt, $$

and

$$ \int |x|^{2m_1} \left| \int \frac{e^{-iux} k^*(uh)(1-\phi(u))}{u\mu_{\varepsilon}^*(u)} du \right|^2 dx \\ = 2\pi \int \left| \left( \frac{k^*(th)(1-\phi(t))}{t\mu_{\varepsilon}^*(-t)} \right)^{(m_1)} \right|^2 dt. $$

Now, for $h \le 1$,

$$ \left| \left( \frac{k^*(th)(1-\phi(t))}{t\mu_{\varepsilon}^*(-t)} \right)^{(m_1)} \right| \le C \sum_{j=0}^{m_1} \sum_{\ell=0}^{j} \frac{|r_{\varepsilon}^{(\ell)}(-t)|}{|t|^{j-\ell+1}} \mathbf{1}_{[-1,1]^c}(t) \\ \le C \sum_{\ell=0}^{m_1} \frac{|r_{\varepsilon}^{(\ell)}(-t)|}{|t|} \mathbf{1}_{[-1,1]^c}(t). $$

Finally,

$$ J \le \frac{C}{\sqrt{n}} \left[ \mathbb{E}|Y|^{2p-\frac{1}{2}} \int_{-1/h}^{1/h} \frac{|r_{\varepsilon}(x)|^2}{|x|^2} \mathbf{1}_{[-1,1]^c}(x) dx \right. \\ \left. + \sum_{\ell=0}^{m_1} \int_{-1/h}^{1/h} \frac{|r_{\varepsilon}^{(\ell)}(x)|^2}{|x|^2} \mathbf{1}_{[-1,1]^c}(x) dx \right]^{1/2} . \quad (36) $$

Starting from (10) and gathering the upper bounds (11), (31), (35) and (36), the proof of Proposition 3.1 is complete.

# 6. Numerical experiments

This section is devoted to the implementation of the deconvolution estimators. We continue the experiments of [4] about Wasserstein deconvolution in the ordinary smooth case. In particular, we study the $W_1$ and $W_2$ univariate deconvolution problems and we compare our numerical results with the upper and lower bounds given in the previous sections. We also apply our procedure to the deconvolution of the uniform measure on the Cantor set. The deconvolution method is implemented in R.

## 6.1. Implementation of the deconvolution estimators

For all the experiments we use the kernel

$$ k(x) = \frac{3}{16\pi} \left( \frac{8 \sin(x/8)}{x} \right)^4 $$
which corresponds to the kernel given by (7) with $p = 2$ and a Fourier support over $[-1/2, 1/2]$. Computing the deconvolution estimators requires to evaluate many times the function

$$ \tilde{k}_h : x \mapsto \frac{1}{2\pi} \int \frac{e^{iux} k^*(u)}{\mu_{\varepsilon}^*(-u/h)} du $$

which is the Fourier transform of

$$ \psi_h : u \mapsto \frac{1}{2\pi} \frac{k^*(u)}{\mu_{\varepsilon}^*(-u/h)}. $$

The Fourier decomposition of $\psi_h$ is given by $\psi_h(u) = \sum_{k \in \mathbb{Z}} a_{k,h} e^{2i\pi k u}$ where $a_{k,h} = \int_{-1/2}^{1/2} \psi_h(u) e^{-2i\pi k u} du$. In this section we consider symmetric distributions for $\mu_{\varepsilon}$. Thus $k^*$ and $\mu_{\varepsilon}^*$ are even functions, and the $a_{k,h}$'s are real coefficients. Next,

$$ \begin{aligned} \tilde{k}_h(x) &= \int_{-1/2}^{1/2} \psi_h(u) e^{ixu} du \\ &= \sum_{k \in \mathbb{Z}} a_{k,h} \int_{-1/2}^{1/2} e^{i(2\pi k + x)u} du \\ &= \sum_{k \in \mathbb{Z}} a_{k,h} \text{sinc} \left( \frac{2\pi k + x}{2} \right). \end{aligned} $$

For large $N$, the coefficient $a_{k,h}$ can be approximated by the $k$-th coefficient of a discrete Fourier transform taken at $(\psi_h(0), \psi_h(1/N), \dots, \psi_h(1-1/N))$, denoted $\hat{a}_{k,h,N}$ in the sequel. Of course we use the Fast Fourier Transform algorithm to compute these quantities. For some large $K$, we evaluate $\tilde{k}_h$ at some point $x$ by

$$ \hat{\tilde{k}}_h(x) \approx \sum_{|k| \le K} \hat{a}_{k,h,N} \text{sinc} \left( \frac{2\pi k + x}{2} \right). \quad (37) $$

For intensive simulation, it may be relevant to preliminary compute $\hat{\tilde{k}}_h$ on a grid of high resolution rather than calling this function each time.

We first define a discrete approximation of the function

$$ \hat{\mu}_{n,h} : u \mapsto \frac{1}{nh} \sum_{k=1}^{n} \tilde{k}_h \left( \frac{u - Y_k}{h} \right). $$

Let $\mathcal{P} = \{t_1 < \dots < t_q\}$ be a finite regular grid of points in $\mathbb{R}$ with resolution $\eta$. A discrete approximation $\hat{\mu}_{n,h}^d$ of $\hat{\mu}_{n,h}$ is defined on $\mathcal{P}$ by

$$ \hat{\mu}_{n,h}^d = \eta \sum_{j=1}^{q} \hat{\mu}_{n,h}(t_j) \delta_{t_j}, $$
where $\delta_x$ is the Dirac distribution at $x$. Since $\hat{\mu}_{n,h}(t_j)$ can be negative, the first method for estimating $\mu$ consists in taking the positive part of $\hat{\mu}_{n,h}^d$:

$$ \hat{\mu}_{n,h}^{\text{naive}} := \frac{\sum_{j=1}^{q} (\hat{\mu}_{n,h}^{d}(t_j))^{+} \delta_{t_j}}{\sum_{j=1}^{q} (\hat{\mu}_{n,h}^{d}(t_j))^{+}} $$

This first estimator is called the “naive” deconvolution estimator henceforth. Note that it was studied in [4] and [10]. For implementing the alternative estimator $\tilde{\mu}_{n,h}$ proposed in this paper, we first need to find some probability distribution $\tilde{F}_{n,h}$ on $\mathbb{R}$ such that

$$ \int_{\mathbb{R}} |x|^{p-1} |\hat{F}_{n,h} - \tilde{F}_{n,h}|(x)dx \\ \approx \inf \left\{ \int_{\mathbb{R}} |x|^{p-1} |\hat{F}_{n,h} - G|(x)dx \mid G \text{ probability distribution on } \mathbb{R} \right\}. \quad (38) $$

In practice, this corresponds to finding a distribution function close to the step function

$$ \hat{F}_{n,h}^{d} : t \mapsto \sum_{j=1}^{q} \hat{\mu}_{n,h}^{d}(t_j) \mathbf{1}_{\{t_j \le t\}}. $$

Since $\hat{F}_{n,h}^d$ may take its values outside $[0, 1]$, we can also look for a distribution function close to $t \mapsto \hat{F}_{n,h}^d(t) \mathbf{1}_{\hat{F}_{n,h}^d(t) \in [0, 1]}$. In other terms, we compute the isotone regression of $t \mapsto \hat{F}_{n,h}^d(t) \mathbf{1}_{\hat{F}_{n,h}^d(t) \in [0, 1]}$ with weights $t_j^{p-1}$:

$$ \hat{F}_{n,h}^{\text{isot},p} := \operatorname{argmin} \left\{ \sum_{j=1}^{q} |t_j|^{p-1} \left| G(t_j) - \hat{F}_{n,h}^{d}(t_j) \mathbf{1}_{\hat{F}_{n,h}^{d}(t_j) \in [0,1]} \right|^{p} \mid G \text{ non-decreasing} \right\}. $$

We compute $\hat{F}_{n,h}^{\text{isot},p}$ thanks to the function `gpava` from the R package `isotonic` [22]. The measure $\mu$ is finally estimated by the absolutely continuous measure $\hat{\mu}_{n,h}^{\text{isot},p}$ whose distribution function is $\hat{F}_{n,h}^{\text{isot},p}$. We call this estimator the isotone deconvolution estimator for the metric $W_p$.

The construction of $\hat{\mu}_{n,h}^{\text{isot},p}$ depends on many parameters, for instance $K$, $h$, $N$ and $\eta$. Tuning all these parameters is a tricky issue. For this paper we only tune these quantity by hand. The bandwidth choice is discussed in Section 6.5. Note that one crucial point is the length $N$ of the vector we use for computing the the $a_{k,h,N}$'s with the FFT. For ordinary smooth distributions, we observe that $\tilde{k}_h$ decreases slowly for small $\beta$ for the range of bandwidths $h$ giving minimum Wasserstein risks. Consequently, a small $\beta$ requires many terms in the expansion (37), and hence a large $N$. For $\beta$ smaller than 0.5, it was necessary to take $N \approx 10^4$.
<figure><img src="image_3.png" /><figcaption>FIG 1. Estimated Wasserstein risks for the Dirac experiment. The noise distribution is the symmetrized Gamma distribution with $\beta = 2$. The twenty curves correspond to samples of length $n$ taken between 100 and 2000.</figcaption></figure>

## 6.2. Computation of Wasserstein risks for simulated experiments

For fixed distributions $\mu$ and $\mu_\varepsilon$, we simulate $Y_1, \dots, Y_n$ according to the convolution model (1). For a given bandwidth $h$ and $p \ge 1$, we can compute $W_p^p(\hat{\mu}_n^{\text{naive}}, \mu)$ and $W_p^p(\hat{\mu}_{n,h}^{\text{isot},p}, \mu)$ using the quantile functions of the measures, thanks to the relation (2). The Wasserstein risks $\mathcal{R}^{\text{naive}}(n, h) := \mathbb{E}W_p^p(\hat{\mu}_{n,h}^{\text{naive}}, \mu)$ and $\mathcal{R}^{\text{isot}}(n, h) := \mathbb{E}W_p^p(\hat{\mu}_{n,h}^{\text{isot},p}, \mu)$ can be estimated by an elementary Monte Carlo method by repeating the simulation of the $Y_i$'s and averaging the Wasserstein distances. Let $\tilde{r}_p^{\text{isot}}(n, h)$ and $\tilde{r}_p^{\text{naive}}(n, h)$ be the estimated risks obtained this way (see Figure 1 for an illustration of such curves for the Dirac experi-
<table><caption>TABLE 1
Ordinary smooth distributions used for the error</caption><thead><tr><th>Distribution</th><th>μ<sub>ε</sub><sup>*</sup></th><th>β</th></tr></thead><tbody><tr><td>Symmetrized Gamma</td><td>t ↦ (1 + t<sup>2</sup>)<sup>-β/2</sup></td><td>0.3, 0.5, 1.2, 2,3,4</td></tr><tr><td>Laplace</td><td>t ↦ (1 + t<sup>2</sup>)<sup>-1</sup></td><td>2</td></tr><tr><td>Symmetrized χ<sup>2</sup></td><td>t ↦ (1 + 4t<sup>2</sup>)<sup>(-1/2)</sup></td><td>1</td></tr></tbody></table>

ment). For each $n$, an approximation of the minimal risks over the bandwidths is proposed by

$$ \bar{r}_{p,*}^{\text{isot}}(n) := \min_{h \in \mathcal{H}} \bar{r}_{p}^{\text{isot}}(n, h) $$

and

$$ \bar{r}_{p,*}^{\text{naive}}(n) := \min_{h \in \mathcal{H}} \bar{r}_{p}^{\text{naive}}(n, h) $$

where $\mathcal{H}$ is a grid of bandwidth values.

### 6.3. Estimation of the rates of convergence

In this experiment we study the rates of convergence of the estimators for the deconvolution of three distributions:

* Dirac distribution at 0,

* Uniform distribution on [-0.5, 0.5],

* Mixture of the Dirac distribution at 0 and the uniform distribution on [-0.5, 0].

We take for $\mu_{\epsilon}$ the ordinary smooth distributions summarized in Table 1. Recall that the coefficient $\beta$ of a symmetrized Gamma distribution is twice the shape parameter of the distribution. For each error distribution and for $n$ chosen between 100 and 2000, we simulate 200 times a sample of length $n$ from which we compute the estimated minimal risks $\bar{r}_{p,*}^{\text{isot}}(n)$ and $\bar{r}_{p,*}^{\text{naive}}(n)$. We study the Wasserstein risks $W_1$ and $W_2$. We obtain some estimation of the exponent of the rate of convergence for each deconvolution problem by computing the linear regression of $\log \bar{r}_{p,*}(n)$ by $\log n$. See Figure 2 for an illustration and Figures 7 and 8 at the end of the paper for the complete outputs of the Dirac case. A linear trend can be observed in all cases. As expected, the risks are smaller for the isotone estimators than for the naive ones.

The estimated exponents of the convergences rates are plotted in Figure 3 as functions of $\beta$. These estimated rates can be compared with the upper and lower bounds obtained in the paper. Of course the rates of convergence of the isotone estimator have no reason to match exactly the lower bounds. However it can be checked that the estimated rates we obtain are consistent with the theoretic bounds proved before. In particular we see that the parametric rate is reached for values of $\beta$ close to 0, at least in the Dirac case. These results also suggest that the correct minimax rate for $W_2$ probably corresponds to the upper bound given in Theorem 3.1 (that is, when no further assumption is made on the unknown distribution $\mu$).
<figure><img src="image_4.png" /><figcaption>FIG 2. Estimated rates of convergence to zero of the W₁-risk for the naive method and the isotone method for μ being a Dirac distribution at 0. The noise distribution is the symmetrized Gamma distribution with β = 2.</figcaption></figure>

6.4. Cantor set experiment

We now illustrate the deconvolution method with a more original experiment. We take for μ the uniform distribution on the Cantor set $\mathbb{C}$. Remember that the Cantor set can be defined by repeatedly deleting the open middle thirds of a set of line segments:

$$ \mathbb{C} = \bigcap_{m \ge 1} F_n $$

where $F_0 = [0, 1]$ and $F_{m+1}$ is obtained by cutting out the middle thirds of all the intervals of $F_m$: $F_1 = [0, \frac{1}{3}] \cup [\frac{2}{3}, 1]$ and $F_2 = [0, \frac{1}{9}] \cup [\frac{2}{9}, \frac{1}{3}] \cup [\frac{2}{3}, \frac{7}{9}] \cup [\frac{8}{9}, 1]$, etc... The uniform measure $\mu_{\mathbb{C}}$ on $\mathbb{C}$ can be defined as the distribution of the random variable $X := 2 \sum_{k \ge 1} 3^{-k} B_k$ where $(B_k)_{k \ge 1}$ is a sequence of independent random variables with Bernoulli distribution of parameter 1/2. Note that the Lebesgue measure of $\mathbb{C}$ is zero and thus the Lebesgue measure and $\mu_{\mathbb{C}}$ are singular. The deconvolution estimators being densities for the Lebesgue measure, the Wasserstein distances are relevant metrics for comparing these with $\mu_{\mathbb{C}}$.

Let $\mu_{\mathbb{C},K}$ be the distribution of the random variable defined by the partial sum $\tilde{X} := 2 \sum_{k=1}^{K} 3^{-k} B_k$ where the $B_k$'s are defined as before. The distribution $\mu_{\mathbb{C},K}$ is an approximation of $\mu_{\mathbb{C}}$ which can be computed in practice. We simulate a sample of $n = 10^4$ observations from $\mu_{\mathbb{C},K}$ with $K = 100$. These observations are contaminated by random variables with symmetrized Gamma distribution (the shape parameter is equal to 1/4 (so that $\beta = 0.5$) and the scale parameter is equal to 1/2).
<figure><img src="image_5.png" /><figcaption>FIG 3. Estimated exponents of the convergence rates of the naive and the isotone deconvolution estimator for three distributions $\mu$. The exponents are given as functions of the ordinary smooth coefficient $\beta$. Regarding the noise distribution, the $\chi^2$ and the Laplace distributions are indicated directly on the graph, the others experiments have been done with the symmetrized Gamma distribution. The top graph corresponds to the $W_1$ deconvolution and the bottom one to the $W_2$ deconvolution.</figcaption></figure>

In Figure 4, the isotone estimators for $W_1$ and $W_2$ and the naive estimator are plotted on the first four levels $F_m$ of the Cantor set. The bandwidths are chosen by minimizing the Wasserstein risks over a grid, as in Section 6.3. This requires to approximate the quantile functions for the isotone deconvolution estimator
<figure><img src="image_6.png" /><figcaption>FIG 4. Deconvolution of the uniform measure on the Cantor set.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG 5. Approximation of the Devil's staircase and distributions functions of the $W_1$ and $W_2$ isotone deconvolution estimators.</figcaption></figure>

and for the $\mu_{\mathbb{C}}$. Regarding the quantile function of $\mu_{\mathbb{C}}$, we simulate a large sample according to $\mu_{\mathbb{C},100}$ and we compute the corresponding empirical distribution function. This last cdf is an approximation of the so called “Devil’s staircase” (see Figure 5). For the naive deconvolution estimator we find $h = 0.011$ for $W_1$ and $h = 0.018$ for $W_2$. For the $W_1$-isotone deconvolution estimator we find
<figure><img src="image_8.png" /><figcaption>FIG 6. Comparison of the locations of the minima of the W₂-risks for five distribution measures µ. The noise distribution µε is the symmetrized Gamma distribution with β = 0.75. Each risk curve has been normalized by its minimum value for facilitating the comparison.</figcaption></figure>

$h = 0.002$ and $h = 0.01$ for the W₂-isotone estimator. Note that these values are consistent with the fact that the bandwidth increases with the parameter $p$ of the Wasserstein metric, as shown by Theorem 3.1. On Figure 4, the W₁-isotone deconvolution estimator is able to “see” the first three levels of the Cantor set and the three other deconvolution methods recover the first two levels. A kernel density estimator (with no deconvolution) only recovers the first level.

6.5. About the bandwidth choice

In practice, we need to choose a bandwidth $h$ for the deconvolution estimators. As was explained in [4] (see Remark 3 in this paper), it seems that the influence of the measure $\mu$ is weak. We now propose a simple experiment to check this principle. We choose for $\mu_{\epsilon}$ the symmetrized gamma distribution with a shape parameter equal to 0.375 ($\beta = 0.75$) and we simulate contaminated observations from the following various distributions:

* Truncated standard Gaussian distribution on [-1, 1],

* Uniform distribution on [-0.5, 0.5],

* Uniform distribution on the Cantor set,

* Mixture of the Dirac distribution at 0 and the uniform distribution on [-0.5, 0],
<figure><img src="image_9.png" /><figcaption>FIG 7. Deconvolution of the Dirac distribution at zero observed with one of the noise distributions listed in Table 1: log-log plots of the estimated $W_1$-risks for the naive method and the isotone method.</figcaption></figure>

• Mixture of Dirac distributions at -0.5, -0.2 and 0.3 with proportions 1/4, 1/4 and 1/2,

• Dirac distribution at 0.

We focus here on the study of the $W_2$-isotone deconvolution estimator. Figure 6 compares the locations of the minimums of the five risk curves $h \mapsto \bar{r}_{2,h}^{\text{isot}}$ by
<figure><img src="image_10.png" /><figcaption>FIG 8. Deconvolution of the Dirac distribution at zero observed with one of the noise distributions listed in Table 1: log-log plots of the estimated W₂-risks for the naive method and the isotone method.</figcaption></figure>

averaging over 200 samples of 1000 contaminated observations. For this experiment, the sensitivity of the minimum risk location to the distribution $\mu$ is not very large.

On another hand, from Figure 3, it seems that the rates for the mixture Dirac Uniform are quite slow (in particular, they are close to the minimax rates for $W_1$).
From these remarks, it seems that the bandwidth minimizing the risk computed for the mixture Dirac Uniform should be a reasonable choice for deconvolving other distributions. Of course, this is in some sense a “minimax choice”, and it will not give the appropriate rate for measures which are easier to estimate (for instance measures with smooth densities).

A bootstrap method in the spirit of [13] may give a more satisfactory answer to this problem. However, note that the use of the Wasserstein metric makes difficult the asymptotical analysis of the risk. This interesting problem is out of the scope of this paper, we intend to investigate it in a future work.

## References

[1] BOBKOV, S. and LEDOUX, M., One-dimensional empirical measures, order statistics and Kantorovich transport distances. Preprint, 2014. MR2797936

[2] BUTUCEA, C. and TSYBAKOV, B., Sharp optimality in density deconvolution with dominating bias. I. *Theory Probab. Appl.*, 52:24–39, 2008a. MR2354572

[3] BUTUCEA, C. and TSYBAKOV, B., Sharp optimality in density deconvolution with dominating bias. II. *Theory Probab. Appl.*, 52:237–249, 2008b. MR2742504

[4] CAILLERIE, C., CHAZAL, F., DEDECKER, J., and MICHEL, B., Deconvolution for the Wasserstein metric and geometric inference. *Electron. J. Stat.*, 5:1394–1423, 2011. MR2851684

[5] CARLSSON, G., Topology and data. *Bull. Amer. Math. Soc.*, 46:255–308, 2009. MR2476414

[6] CARROLL, R.J. and HALL, P., Optimal rates of convergence for deconvolving a density. *J. Amer. Statist. Assoc.*, 83:1184–1186, 1988. MR0997599

[7] CHAZAL, F., COHEN-STEINER, D., and MÉRIGOT, Q., Geometric inference for probability measures. *Found. Comput. Math.*, 11:733–751, 2011. MR2859954

[8] CHAZAL, F., FASY, B.T., LECCI, F., MICHEL, B., RINALDO, A., and WASSERMAN, L., Subsampling methods for persistent homology. *arXiv:1406.1901*, 2014.

[9] DATTNER, I., GOLDENSHLUGER, A., and JUDITSKY, A., On deconvolution of distribution functions. *Ann. Statist.*, 39:2477–2501, 2011. MR2906875

[10] DEDECKER, J. and MICHEL, B., Minimax rates of convergence for Wasserstein deconvolution with supersmooth errors in any dimension. *J. Multivar. Anal.*, 122:278–291, 2013. MR3189324

[11] DEL BARRIO, E., GINÉ, E., and MATRÁN, C., The central limit theorem for the Wasserstein distance between the empirical and the true distributions. *Ann. Probab.*, 27:1009–1971, 1999. MR1698999

[12] DEL BARRIO, E., GINÉ, E., and UTZET, F., Asymptotics for $L_2$ functionals of the empirical quantile process, with applications to tests of fit based on weighted Wasserstein distances. *Bernoulli*, 11:131–189, 2005. MR2121458
[13] DELAIGLE, A. and GIJBELS, I., Bootstrap bandwidth selection in kernel density estimation from a contaminated sample. *Ann. I. Stat. Math.*, 56(1):19–47, 2004. MR2053727

[14] DEREICH, S., SCHEUTZOW, M., and SCHOTTSTEDT, R., Construc- tive quantization: Approximation by empirical measures. *Ann. Inst. H. Poincaré Probab. Statist.*, 49:1183–1203, 2013. MR3127919

[15] ĒBRALIDZE, Š.S., Inequalities for the probabilities of large deviations in terms of pseudomoments. *Teor. Verojatnost. i Primenen.*, 16:760–765, 1971. MR0295417

[16] FAN, J., Global behavior of deconvolution kernel estimates. *Statist. Sinica*, 2:541–551, 1991a. MR1130132

[17] FAN, J., On the optimal rates of convergence for nonparametric deconvolution problems. *Ann. Stat.*, 19:1257–1272, 1991b. MR1126324

[18] FAN, J., Adaptively local one-dimensional subproblems with application to a deconvolution problem. *Ann. Stat.*, 21:600–610, 1993. MR1232507

[19] FOURNIER, N. and GUILLIN, A., On the rate of convergence in Wasserstein distance of the empirical measure. *To appear in Probability Theory and Related Fields*, 2014.

[20] GUIBAS, L., MOROZOV, D., and MÉRIGOT, Q., Witnessed k-distance. *Discrete Comput. Geom.*, 49:22–45, 2013. MR3010216

[21] HALL, P. and LAHIRI, S.N., Estimation of distributions, moments and quantiles in deconvolution problems. *Ann. Statist*, 36:2110–2134, 2008. MR2458181

[22] MAIR, P., HORNIK, K., and DE LEEUW, J., Isotone optimization in R: pool-adjacent-violators algorithm (PAVA) and active set methods. *J. Stat. Softw.*, 32(5):1–24, 2009.

[23] MEISTER, A., *Deconvolution Problems in Nonparametric Statistics*. Lecture Notes in Statistics. Springer, 2009. MR2768576

[24] RACHEV, S.T. and RÜSCHENDORF, L., *Mass Transportation Problems*, volume II of *Probability and Its Applications*. Springer-Verlag, 1998. MR1619171

[25] VAN DER VAART, A.W. and WELLNER, J.A., *Weak Convergence and Empirical Processes*. Springer Series in Statistics. Springer, 1996. MR1385671

[26] VILLANI, C., *Optimal Transport: Old and New*. Grundlehren Der Mathematischen Wissenschaften. Springer-Verlag, 2008. MR2459454