Bounding quantiles of Wasserstein distance between true and empirical measure

TL;DR

This paper analyzes the asymptotic behavior of Wasserstein distance quantiles between true and empirical measures, identifying worst-case distributions and providing confidence regions, with extensions to higher dimensions.

Contribution

It characterizes the asymptotic maximizers of Wasserstein quantiles and derives explicit confidence regions for the true distribution, extending to higher dimensions with numerical support.

Findings

Asymptotic maximizers are convex combinations of two-point and uniform distributions.

Explicit confidence regions for the true measure are derived.

Numerical evidence supports extensions to higher dimensions.

Abstract

Consider the empirical measure, $\hat{\mathbb{P}}_N$, associated to $N$ i.i.d. samples of a given probability distribution $\mathbb{P}$ on the unit interval. For fixed $\mathbb{P}$ the Wasserstein distance between $\hat{\mathbb{P}}_N$ and $\mathbb{P}$ is a random variable on the sample space $[0,1]^N$. Our main result is that its normalised quantiles are asymptotically maximised when $\mathbb{P}$ is a convex combination between the uniform distribution supported on the two points $\{0,1\}$ and the uniform distribution on the unit interval $[0,1]$. This allows us to obtain explicit asymptotic confidence regions for the underlying measure $\mathbb{P}$. We also suggest extensions to higher dimensions with numerical evidence.

Equations

$$\mathcal{W}(\text{mathds}P,\widetilde{\text{mathds}P}):=\underset{\pi \in \Pi (\text{mathds}P,\widetilde{\text{mathds}P})}{\inf }\left({\int }_{\mathbb{R}}|x-y|\pi (dx,dy)\right),$$

$${\widehat{\text{mathds}P}}_N:=\frac{1}{N}\sum _{k=1}^N{\delta }_{X_k},$$

$$\text{mathds}{P}^{\otimes }\left(\mathcal{W}({\widehat{\text{mathds}P}}_N,\text{mathds}P)\ge t\right)\le C\exp (-cN{t}^2)\forall \text{mathds}P\in \mathfrak{P}[0,1]$$

$$F^{-1}_{\mathcal{W}(\hat{\mathds{P}}_{N},\mathds{P})}(\alpha)\mathrel{\mathop{:}}=\inf\big{\{}x\in[0,1]\ |\ \mathds{P}^{\otimes}(\mathcal{W}(\hat{\mathds{P}}_{N},\mathds{P})\leq x)\geq\alpha\big{\}}$$

$$\underset{\text{mathds}P\in \mathfrak{P}([0,1])}{\mathrm{argmax}}{F}_{\mathcal{W}({\widehat{\text{mathds}P}}_N,\text{mathds}P)}^{-1}(\alpha ).$$

$$\text{mathds}{P}^{\lambda }:=\lambda \frac{{\delta }_0+{\delta }_1}{2}+(1-\lambda )\mathcal{U}([0,1])$$

$$\underset{\text{mathds}P\in \mathfrak{P}[0,1]}{\sup }{F}_{\mathcal{W}({\widehat{\text{mathds}P}}_N,\text{mathds}P)}^{-1}(\alpha )-F_{\mathcal{W}({\widehat{\text{mathds}P}}_N,\text{mathds}{P}^{\lambda })}^{-1}(\alpha )\le \frac{\epsilon }{N^{1/2}}.$$

$$N^{1/2}\mathcal{W}({\widehat{\text{mathds}P}}_N,\text{mathds}P)\Rightarrow {\int }_0^1|B(F_{\text{mathds}P}(t))|dt$$

$$x_1=0<x_2=\frac{1}{n-1}<x_3=\frac{2}{n-1}<\cdots <x_n=1.$$

$$N^{1/2}\mathcal{W}({\widehat{\text{mathds}P}}_N,\text{mathds}P)\Rightarrow \left(\sum _{i=1}^{n-1}\left|B(q_i)\right|(x_{i+1}-x_i)\right)=\sum _{i=1}^{n-1}\frac{\left|B(q_i)\right|}{n-1}$$

$${\mathfrak{B}}_n:=\sum _{i=1}^{n-1}\frac{|B(q_i)|}{n-1}$$

$${\mathbf{\Sigma }}_n^{p_1,\dots ,p_n}:=\left(\begin{array}{ccccc}{q}_1{\bar{q}}_1& q_1{\bar{q}}_2& q_1{\bar{q}}_3& \dots & q_1{\bar{q}}_{n-1}\\ q_1{\bar{q}}_2& q_2{\bar{q}}_2& q_2{\bar{q}}_3& \dots & q_2{\bar{q}}_{n-1}\\ ⋮& ⋮& ⋮& \dots & ⋮\\ q_1{\bar{q}}_{n-1}& q_2{\bar{q}}_{n-1}& q_3{\bar{q}}_{n-1}& \dots & q_{n-1}{\bar{q}}_{n-1}\end{array}\right).$$

$$\text{mathds}{P}^{\lambda (\alpha )}:=\lambda (\alpha )\frac{{\delta }_0+{\delta }_1}{2}+(1-\lambda (\alpha ))\mathcal{U}(\{x_1,\dots ,x_n\})$$

$$\Xi (t):=\{(y_1,\dots ,y_{n-1})\in {\mathbb{R}}^{n-1}||y_1|+\cdots +|y_{n-1}|\le t(n-1)\}\text{for }t\ge 0,$$

$$F_{{\mathfrak{B}}_n}(t)$$

$$={\int }_{\Xi (t)}\frac{1}{\sqrt{2\pi p_1}}{e}^{\frac{-y_1^2}{2p_1}}\frac{1}{\sqrt{2\pi p_2}}{e}^{\frac{-(y_2-y_1{)}^2}{2p_2}}\cdots$$

$$\cdot \frac{1}{\sqrt{2\pi p_{n-1}}}{e}^{\frac{-(y_{n-1}-y_{n-2}{)}^2}{2p_{n-1}}}\frac{1}{\sqrt{p_n}}{e}^{\frac{-y_{n-1}^2}{2p_n}}d{y}_1\cdots d{y}_{n-1}.$$

$$C(y_1,\dots ,y_{n-1}):=\frac{1}{\sqrt{2\pi p_2}}{e}^{\frac{-(y_2-y_1{)}^2}{2p_2}}\dots \frac{1}{\sqrt{2\pi p_{n-1}}}{e}^{\frac{-(y_{n-1}-y_{n-2}{)}^2}{2p_{n-1}}},$$

$$\frac{\partial }{\partial p_1}\left({\int }_{\Xi (t)}\frac{1}{\sqrt{2\pi p_1(1-p_1-c)}}{e}^{-\left(\frac{y_1^2}{2p_1}+\frac{y_{n-1}^2}{2(1-p_1-c)}\right)}C(y_1,\dots ,y_{n-1})d{y}_1\dots d{y}_n\right)$$

$$={\int }_{\Xi (t)}\frac{1}{\sqrt{2\pi }}{e}^{-\left(\frac{y_1^2}{2p_1}+\frac{y_{n-1}^2}{2(1-p_1-c)}\right)}C(y_1,\dots ,y_{n-1})\cdot$$

$$\displaystyle\qquad\Bigg{(}-\frac{1-2p_{1}-c}{2(p_{1}(1-p_{1}-c))^{3/2}}+\frac{-\frac{y_{1}^{2}}{2p_{1}^{2}}+\frac{y_{n-1}^{2}}{2(1-p_{1}-c)^{2}}}{\sqrt{p_{1}(1-p_{1}-c)}}\Bigg{)}dy_{1}\dots dy_{n-1}.$$

$$\displaystyle\int_{\Xi(t)}\frac{1}{\sqrt{2\pi}}e^{-\left(\frac{y_{1}^{2}}{2(1-c)}+\frac{y_{n-1}^{2}}{2(1-c)}\right)}C(y_{1},\dots,y_{n-1})\Bigg{(}\frac{-\frac{y_{1}^{2}}{(1-c)^{2}/2}+\frac{y_{n-1}^{2}}{(1-c)^{2}/2}}{(1-c)/2}\Bigg{)}dy=0$$

$$\frac{\partial }{\partial p_1}{F}_{{\mathfrak{B}}_n}(t;p_1)>0$$

$$\frac{\partial }{\partial p_1}{F}_{{\mathfrak{B}}_n}(t;p_1)$$

$$\ge$$

$$\frac{\partial }{\partial p_1}$$

$$\text{mathds}{P}^{\lambda }=\lambda \frac{{\delta }_0+{\delta }_1}{2}+(1-\lambda )\text{leavevmode}\mathcal{U}(\{x_1,\dots ,x_n\})$$

$$p=\left(p_1,\frac{1-2p_1}{n-2},\dots ,\frac{1-2p_1}{n-2},p_1\right).$$

$$\displaystyle\frac{\partial}{\partial p_{1}}\Bigg{(}\int_{\Xi(t)}\frac{p_{1}^{-1}}{(2\pi)^{(n-1)/2}}\left(\frac{1-2p_{1}}{n-2}\right)^{-(n-2)/2}$$

$$\displaystyle\qquad\cdot\exp\bigg{(}-\frac{1}{2}\bigg{(}\frac{1}{p_{1}}(y_{1}^{2}+y_{n-1}^{2})+\frac{n-2}{1-2p_{1}}\sum_{i=2}^{n-1}(y_{i}-y_{i-1})^{2}\bigg{)}\bigg{)}dx\Bigg{)}$$

$$=\frac{1}{(2\pi {)}^{(n-1)/2}}\frac{1}{p_1{\left(\frac{1-2p_1}{n-2}\right)}^{(n-2)/2}}\cdot$$

$$\displaystyle\qquad\int_{\Xi(t)}\bigg{(}\frac{np_{1}-1}{(1-2p_{1})p_{1}}+\frac{y_{1}^{2}+y_{n-1}^{2}}{2p_{1}^{2}}-\frac{n-2}{(1-2p_{1})^{2}}\sum_{i=2}^{n-1}(y_{i}-y_{i-1})^{2}\bigg{)}\cdot$$

$$\tilde{C}(y_1,\dots ,y_{n-1})d{y}_1\dots d{y}_{n-1},$$

$$\begin{split}0&\leq\tilde{C}(y_{1},\dots,y_{n})\\ &\mathrel{\mathop{:}}=\exp\bigg{(}-\frac{1}{2}\left(\frac{1}{p_{1}}(y_{1}^{2}+y_{n-1}^{2})+\frac{n-2}{1-2p_{1}}\left((y_{2}-y_{1})^{2}+\dots+(y_{n-1}-y_{n-2})^{2}\right)\right)\bigg{)}\end{split}$$

$${\alpha }_1\ge {\alpha }_2\ge \cdots \ge {\alpha }_d.$$

Full text

We also suggest extensions to higher dimensions with numerical evidence.

MSC: 60F05, 62G15, 60E15

Keywords: Wasserstein distance, confidence interval, central limit theorem, stochastic dominance, empirical measure, concentration inequalities

1. Notation and main result

Many problems in statistics depend on understanding how close the empirical distribution of a sample is from its generating distribution. In this paper, we show that, for measures on $[0,1]$ with distances measured by the $1$-Wasserstein metric, the quantiles of this (random) distance are maximised by a specific family of sampling distributions. This yields a tight upper bound on the law of the Wasserstein distance between the true and empirical measures.

Let $\mathfrak{P}([0,1])$ denote the space of probability measures on $[0,1]$. For two measures $\mathds{P},\tilde{\mathds{P}}\in\mathfrak{P}([0,1])$ we denote the $1$-Wasserstein distance by

[TABLE]

where $\Pi(\mathds{P},\tilde{\mathds{P}})$ is the set of couplings $\pi\in\mathfrak{P}([0,1]^{2})$ with marginals $\mathds{P}$ and $\tilde{\mathds{P}}$. For independent samples $X_{1},\dots,X_{N}$ distributed according to $\mathds{P}$ we define the empirical measure

[TABLE]

where $\delta_{X_{k}}$ denotes the Dirac measure at $X_{k}$. As $\hat{\mathds{P}}_{N}$ depends on the realisations $X_{1},\dots,X_{N}$, the observed distance $\mathcal{W}(\hat{\mathds{P}}_{N},\mathds{P})$ is a random variable with domain $[0,1]^{N}$. We write $\mathds{P}^{\otimes}=\mathds{P}\otimes\mathds{P}\otimes\cdots$ for the product measure on $[0,1]^{\mathbb{N}}$, the space of all sequences of observations.

Denote weak convergence of a sequence of measures $(\mathds{P}_{N})_{N\in\mathbb{N}}$ to a measure $\mathds{P}\in\mathfrak{P}([0,1])$ by $\mathds{P}_{N}\Rightarrow\mathds{P}$, and similarly $Y_{N}\Rightarrow Y$ for convergence in law (under $\mathds{P}^{\otimes}$) of random variables. While the well-known Glivenko–Cantelli theorem gives the weak convergence $\hat{\mathds{P}}_{N}\Rightarrow\mathds{P}$ ($\mathds{P}^{\otimes}$-a.s.), the problem of controlling the distance between $\hat{\mathds{P}}_{N}$ and $\mathds{P}$ has a long history in probability theory, statistics and many related fields, see for example [9, 7, 17, 26, 1, 22, 4, 14, 19, 5, 11, 12] and the references therein for recent results in quantisation, particle systems and MCMC to name just a few. Building on the results of Dereich, Scheutzow and Schottstedt [8], non-asymptotic concentration inequalities of the form

[TABLE]

were recently obtained by Fournier and Guillin [11], however the constants $c,C$ are not explicit. In statistical applications, the values of these constants along with their (non-)optimality are of importance, for example, for the calculation of confidence intervals.

In this article we try to shed light on this matter by investigating the distribution of $\mathcal{W}(\hat{\mathds{P}}_{N},\mathds{P})$ for large $N\in\mathbb{N}$. Specifically, we aim to determine the probability measure $\mathds{P}$ that maximises the quantile

[TABLE]

for a given confidence level $\alpha\in[0,1]$, i.e.

[TABLE]

Our main result is the following:

Theorem 1.1**.**

Let $\alpha\in[0,1]$ and let $\ \mathcal{U}([0,1])$ denote the uniform distribution on $[0,1]$. Take $\varepsilon>0$. Then there exists $\lambda=\lambda(\alpha,\varepsilon)\in[0,1]$ such that

[TABLE]

asymptotically maximises the normalised quantile $N^{1/2}F^{-1}_{\mathcal{W}(\hat{\mathds{P}}_{N},\mathds{P})}(\alpha)$, i.e. there exists $N_{0}\in\mathbb{N}$ such that for all $N\geq N_{0}$

[TABLE]

Furthermore $\lim_{\alpha\downarrow 0}\lambda(\alpha,\varepsilon)=0$ and $\lim_{\alpha\uparrow 1}\lambda(\alpha,\varepsilon)=1$ for all $\varepsilon>0$.

Remark 1.2*.*

We note that the interval $[0,1]$ in the statement above can be replaced by any closed and bounded interval. Furthermore Theorem 1.1 can be extended to the case where we only consider $\mathds{P}\in\mathfrak{P}([0,1])$ with a density bounded from below. We then restrict to $\lambda\in[0,b]$, where $b\in[0,1]$ is chosen such that $\mathds{P}^{\lambda}$ fulfils this constraint.

Our method of proof also gives a version of this result when measures are supported on finitely many (equally spaced) points, in which case, $\mathcal{U}([0,1])$ is replaced by a uniform measure over these points.

The remainder of this article is structured as follows: Section 2 is concerned with the proof of Theorem 1.1 for finitely supported measures. This is achieved via the use of a functional central limit theorem, approximating $N^{1/2}\mathcal{W}(\hat{\mathds{P}}_{N},\mathds{P})$ in law by a certain integral over a Brownian bridge. As it turns out, we can explicitly find the measure $\mathds{P}$ maximising the quantile of this integral. We then extend this result to general measures $\mathds{P}$ on the unit interval $[0,1]$ in Section 3 using again an approximation argument and a particularly powerful concentration inequality by Talagrand [23]. We also show some numerical results in Section 4 and discuss further applications and possible extensions of Theorem 1.1 in Section 5.

2. Brownian bridge approximation and quantile order for measures with finite support

Even though it is explicitly known, the (rescaled) distribution of $\mathcal{W}(\hat{\mathds{P}}_{N},\mathds{P})$ is hard to control for moderate to large $N\in\mathbb{N}$. We thus resort to asymptotic analysis. The following result has a long history in statistical theory and is crucial for our results:

Lemma 2.1** (del Bario, Giné , Matrán,[6, Theorem 2.1]).**

For $\mathds{P}\in\mathfrak{P}([0,1])$ we have

[TABLE]

for $N\to\infty$, where $(B(q))_{0\leq q\leq 1}$ is a standard Brownian bridge, $F_{\mathds{P}}$ the distribution function associated with $\mathds{P}$ and $\Rightarrow$ indicates convergence in law under the measure $\mathds{P}^{\otimes}$.

Instead of finding a solution to (1.1) for general measures $\mathds{P}\in\mathfrak{P}([0,1])$, we first restrict to finitely supported measures: To this end, let us fix $n\in\mathbb{N}$ and take equidistant points

[TABLE]

The set of all measures $\mathds{P}\in\mathfrak{P}([0,1])$ supported on $\{x_{1},\dots,x_{n}\}$ is then characterised by the probabilities $p_{i}\mathrel{\mathop{:}}=\mathds{P}(\{x_{i}\})$, where $p_{1},\dots,p_{n}\in[0,1]$ with $\sum_{i=1}^{n}p_{i}=1$. For such a measure $\mathds{P}=\sum_{i=1}^{n}p_{i}\delta_{x_{i}}$, Lemma 2.1 reads

[TABLE]

where $q_{1}=p_{1}$ and $q_{i}=p_{1}+\dots+p_{i}$. For notational simplicity, we write

[TABLE]

in this case. Furthermore we set $\bar{q}_{i}=1-q_{i}$ and define the covariance matrix

[TABLE]

It is well-known that the discretised Brownian bridge $(B(q_{i}))_{1\leq i\leq n-1}$ has covariance matrix $\boldsymbol{\Sigma}_{n}^{p_{1},\dots,p_{n}}$ and, as its mean is zero, $\boldsymbol{\Sigma}_{n}^{p_{1},\dots,p_{n}}$ completely determines its distribution. By the peculiar structure of the covariance matrix $\boldsymbol{\Sigma}_{n}^{p_{1},\dots,p_{n}}$ there is in fact a one-to-one correspondence between $(p_{1},\dots,p_{n})$ and $\boldsymbol{\Sigma}_{n}^{p_{1},\dots,p_{n}}$ (up to symmetry). It is thus sufficient to investigate the tractable properties of $\boldsymbol{\Sigma}_{n}^{p_{1},\dots,p_{n}}$ in order to determine the dependence of $F^{-1}_{\mathfrak{B}_{n}}(\alpha)$ on $(p_{1},\dots,p_{n})$.

Note that, although $\mathcal{W}(\hat{\mathds{P}}_{N},\mathds{P})$ is bounded by one, the support of its normalised version $N^{1/2}\mathcal{W}(\hat{\mathds{P}}_{N},\mathds{P})$ is unbounded for $N\to\infty$; the same is true for the distribution of $\mathfrak{B}_{n}$.

The proof of Theorem 1.1 for finitely supported measures is carried out in three steps: We first establish (1.2) in Proposition 2.2 and then prove the limiting behaviour of $\lambda(\alpha)$ for $\alpha\downarrow 0$ and $\alpha\uparrow 1$ in Propositions 2.3 and 2.6.

Proposition 2.2**.**

Let $\alpha\in[0,1]$ and let $\mathcal{U}(\{x_{1},\dots,x_{n}\})$ denote the uniform distribution on the equidistant points $\{x_{1},\dots,x_{n}\}$. Then there exists $\lambda(\alpha)\in[0,1]$ such that

[TABLE]

maximises the quantile $F^{-1}_{\mathfrak{B}_{n}}(\alpha)$.

Proof.

We note that maximising the quantile $F^{-1}_{\mathfrak{B}_{n}}(\alpha)$ for a given $\alpha\in[0,1]$ is equivalent to minimising the distribution function $F_{\mathfrak{B}_{n}}(t)$ for a given $t\in\mathbb{R}_{+}$. Given a vector $(p_{1},\dots,p_{n})$ of probabilities, as described above, the density of $\mathfrak{B}_{n}$ is fully specified: indeed define the ball

[TABLE]

then

[TABLE]

For fixed $t\geq 0$, this expression is symmetric in $p_{1}$ and $p_{n}$, and in the middle terms for $p_{2},\dots,p_{n-1}$ respectively. We now argue that it is optimal to choose $p_{1}=p_{n}$ and $p_{2}=\dots=p_{n-1}$. Indeed, let us fix $p_{2},\dots,p_{n-1}$, $t\geq 0$ and define the constant $c\mathrel{\mathop{:}}=p_{2}+\dots+p_{n-1}\geq 0$ as well as the function

[TABLE]

which is independent of $p_{1}$. Differentiating (2) with respect to $p_{1}$ yields

[TABLE]

Plugging in $p_{1}=(1-c)/2$ we obtain

[TABLE]

by symmetry. We now show that $p_{1}\mapsto(\partial/\partial p_{1})F_{\mathfrak{B}_{n}}(t;p_{1})$ has a global minimum in $p_{1}$. For this we first show that

[TABLE]

for all $p_{1}\in((1-c)/2,1-c]$, which follows from elementary calculus: indeed we first note that the functions $p_{1}\mapsto-(1-2p_{1}-c)$, $p_{1}\mapsto-y_{1}^{2}/(2p_{1}^{2})$ and $p_{1}\mapsto y_{n-1}^{2}/(2(1-p_{1}-c)^{2})$ are strictly increasing for every $y_{1},y_{n-1}\in\mathbb{R}\setminus\{0\}$. As $p_{1}\mapsto 1/(p_{1}(1-p_{1}-c))$ is strictly increasing for $p_{1}>(1-c)/2$ and $C(y_{1},\dots,y_{n-1})>0$ we thus conclude by (2) that for some strictly positive function $p_{1}\mapsto\tilde{C}(p_{1})$

[TABLE]

for all $p_{1}\in((1-c)/2,1]$. Next considering $p_{1}\in[0,(1-c)/2)$, we see that (up to exchanging the roles of $y_{1}$ and $y_{n-1}$, which does not affect the sign) the function $p_{1}\mapsto(\partial/\partial p_{1})F_{\mathfrak{B}_{n}}(t;p_{1})$ is antisymmetric about $((1-c)/2,0)$. We conclude that $p_{1}\mapsto F_{\mathfrak{B}_{n}}(t)$ has a global minimum at $p_{1}=(1-c)/2$.

By a similar argument we obtain $p_{2}=\dots=p_{n-1}$. Finally $p_{1}\geq p_{2}$ follows as, again by symmetry considerations, the derivative

[TABLE]

evaluated at $p_{1}=(1-c)/2$ is negative: indeed, this can be seen by the same calculation as in (2) replacing $y_{n-1}$ with $(y_{2}-y_{1})$ and noting the asymmetry of the terms involving $y_{1}$ and $y_{2}$ and the fact that the terms $y_{1}\mapsto\exp(-y^{2}_{1}/(2p_{1}))$ and $y_{2}-y_{1}\mapsto\exp(-(y_{2}-y^{2}_{1})^{2}/(2(1-p_{1}-c)))$ are radially symmetric and strictly decreasing as a function of the radius $y_{1}^{2}$ and $(y_{2}-y_{1})^{2}$ respectively. ∎

In order to finish the proof of Theorem 1.1 for finitely supported measures, we now investigate the behaviour of $\lambda(\alpha)$ when $\alpha$ approaches one of the two boundary points $\{0,1\}$. This is trivial whenever $n\leq 2$. We start with the case $\alpha\downarrow 0$:

Proposition 2.3**.**

Let $n\geq 3$. We have $\lim_{\alpha\downarrow 0}\lambda(\alpha)=0$ in Proposition 2.2.

Proof.

By Proposition 2.2 it is sufficient to restrict to measures

[TABLE]

for $\lambda\in[0,1]$. Note that any such distribution can be identified with a probability vector

[TABLE]

We aim to find $p_{1}$ which minimises $\mathds{P}((B(q_{1}),\dots,B(q_{n-1}))\in\Xi(t))$ for small $t>0$. Substituting and taking the derivative in $p_{1}$ in (2) we have

[TABLE]

where

[TABLE]

and $p_{1}\in(1/n,1/2)$. Carefully examining (2) we conclude that for $t=t(p_{1})>0$ small enough the derivative $(\partial/\partial{p_{1}})F_{\mathfrak{B}_{n}}(t)$ is positive. The claim follows. ∎

Next we address the case $\lim_{\alpha\uparrow 1}\lambda(\alpha)$.

Lemma 2.4**.**

The entries $(q_{i}\bar{q}_{j})_{1\leq i\leq j\leq n}$ of the covariance matrix $\boldsymbol{\Sigma}_{n}^{p_{1},\dots,p_{n}}$ are maximised componentwise by the choice $p_{1}=p_{n}=1/2$ and $p_{i}=0$ for $i=2,\dots,n-1$. For this choice $q_{i}\bar{q}_{j}=1/4$ for all $1\leq i\leq j\leq n$.

Proof.

Evidently $q_{i}\bar{q}_{j}\leq q_{i}\bar{q}_{i}\leq 1/4$ for $j\geq i$ and the claim follows. ∎

Before stating the next lemma let us recall the standard notation $f\sim g$ iff $\lim_{t\to\infty}f(t)/g(t)=1$ and $f\lesssim g$ iff $\limsup_{t\to\infty}f(t)/g(t)\leq 1$ for functions $f,g:\mathbb{R}\to\mathbb{R}$.

Lemma 2.5**.**

Consider $\boldsymbol{X}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$, where $\boldsymbol{\Sigma}\in\mathbb{R}^{d\times d}$. Let $\{\alpha_{j}\ |\ j=1,\dots,d\}$ be the eigenvalues of $\boldsymbol{\Sigma}$ such that

[TABLE]

Let $m$ be the multiplicity of $\alpha_{1}$ and denote by $\Gamma(s)=\int_{0}^{\infty}t^{s-1}e^{-t}dt$ the Gamma function. Then, as $t\to\infty$,

[TABLE]

Furthermore, for fixed $\varepsilon>0$, (2.5) holds uniformly in $\{(\alpha_{1},\dots,\alpha_{d})\ |\ 1/\varepsilon\geq\alpha_{1}\geq\alpha_{1}-\varepsilon\geq\alpha_{m+1}\geq\cdots\geq\alpha_{d}\}$.

Proof.

The first claim is stated in [13, Theorem 1] and [18, Theorem 3.6 & Remark 3.8]. In particular, [13, Proof of Theorem 1] proceeds as follows: Noting that

[TABLE]

where $U_{i}$ are i.i.d. normal random variables, approximation results for the tails of $\chi^{2}$-distributions can be manipulated to show the claim. A careful examination of the estimates in the proof show that these approximations are indeed uniform in $\{(\alpha_{1},\dots,\alpha_{d})\ |\ 1/\varepsilon\geq\alpha_{1}\geq\alpha_{1}-\varepsilon\geq\alpha_{m+1}\geq\cdots\geq\alpha_{d}\}$. ∎

Proposition 2.6**.**

Let $n\geq 3$. We have $\lim_{\alpha\uparrow 1}\lambda(\alpha)=1$ in Proposition 2.2.

Proof.

As in Lemma 2.5 we denote the eigenvalues of

[TABLE]

by $\{\alpha_{i}(\lambda)\ |\ i=1,\dots,n-1\}$, where

[TABLE]

Lemma 2.5 implies that for $t\to\infty$

[TABLE]

where $\mathds{P}^{\lambda}$ are the measures defined in Proposition 2.2, and

[TABLE]

By the Perron–Frobenius theorem, $\alpha_{1}(\lambda)$ has multiplicity one and is less than the maximum over the sum of its rows, in particular by Lemma 2.4 it is strictly less than $(n-1)/4$ for $\lambda<1$. Using this fact, together with continuity of the map

[TABLE]

and thus also continuity of the eigenvalues $\{\alpha_{i}(\lambda)\ |\ i=1,\dots,n-1\}$ as a function of $\lambda$, and compactness of $[0,1]$, we conclude that the term

[TABLE]

is uniformly bounded in $\lambda\in[0,1]$.

Let us write $F_{\mathfrak{B}_{n}}(t;\lambda)$ for the corresponding distribution functions of $\mathfrak{B}_{n}=\sum_{i}|B(q_{i})|/(n-1)$, where $q$ is chosen as in Lemma 2.1, and $B$ is a Brownian bridge under $\mathds{P}^{\lambda}$. By Hölder’s inequality we have $\|x\|_{\ell_{1}}\leq\sqrt{n-1}\|x\|_{\ell_{2}}$ for $x\in\mathbb{R}^{n-1}$. We then compute

[TABLE]

As $\mathfrak{B}_{n}$ is supported on $\mathbb{R}_{+}$, writing

[TABLE]

for the minimiser at a point $t\in[0,\infty)$, in order to show that

[TABLE]

it is sufficient to show that every sequence $(t_{k})_{k\in\mathbb{N}}$ converging to infinity has a subsequence $(t_{k_{l}})_{l\in\mathbb{N}}$ such that $\lim_{l\to\infty}\bar{\lambda}(t_{k_{l}})=1$. But if $\lim_{l\to\infty}\bar{\lambda}(t_{k_{l}})<1$ then there exists $\varepsilon>0$ and $l_{0}\in\mathbb{N}$ such that $\bar{\lambda}(t_{k_{l}})\leq 1-\varepsilon$ for all $l\geq l_{0}$. By compactness of the interval $[0,1-\varepsilon]$ there furthermore exists $\delta>0$ such that

[TABLE]

for all $\bar{\lambda}\in[0,1-\varepsilon]$. By (2) above and Lemma 2.5, we now conclude that, for $l\to\infty$,

[TABLE]

where the last line follows from noting that

[TABLE]

and $C(\bar{\lambda})$ is bounded in $\bar{\lambda}\in[0,1]$. This leads to a contradiction and so we conclude $\bar{\lambda}(t)\to 1$. ∎

Remark 2.7*.*

By rearranging (2), as $F_{\mathfrak{B}_{2}}(t;\lambda)\equiv F_{\mathfrak{B}_{2}}(t;1)$ we see that as $t\to\infty$

[TABLE]

in particular there exists a continuous function $\bar{C}$ with $\lim_{t\to\infty}\bar{C}(t)=1$ such that

[TABLE]

As a consequence of the results above, we can identify the maximal choice of $F_{\mathfrak{B}_{n}}(t)$ in terms of second order stochastic dominance. For this we recall the definition of a directionally convex function $f:\mathbb{R}^{d}\to\mathbb{R}$:

Definition 2.8**.**

Define the difference operator

[TABLE]

where $e_{i}$ is the $i^{\text{th}}$ unit vector and $\varepsilon>0$. A function $f:\mathbb{R}^{d}\to\mathbb{R}$ is called directionally convex if $\Delta_{i}^{\varepsilon}\Delta_{j}^{\delta}f(x)\geq 0$ for all $x\in\mathbb{R}^{n}$, $1\leq i,j\leq n$ and all $\varepsilon,\delta>0$.

We note that directional convexity neither implies nor is implied by conventional convexity for $d\geq 2$.

Remark 2.9*.*

Even for $n=3$, the determinant of $\boldsymbol{\Sigma}_{3}^{1/2,0,1/2}-\boldsymbol{\Sigma}_{3}^{1/3,1/3,1/3}$ is negative, so there is no hope of using results connecting positive semi-definiteness of matrices and stochastic order. In particular Bergmann [3, Proposition 2.5] does not apply.

We use the following result:

Lemma 2.10** (Müller, [15, Theorem 12]).**

Take vectors $\mathbf{X}\sim\mathcal{N}(\boldsymbol{\mu},\mathbf{\Sigma})$ and $\mathbf{Y}\sim\mathcal{N}(\boldsymbol{\mu}^{\prime},\mathbf{\Sigma}^{\prime})$. Then $\mathbb{E}(f(\boldsymbol{X}))\leq\mathbb{E}(f(\boldsymbol{Y}))$ for all directionally convex functions $f$ if and only if $\boldsymbol{\mu}=\boldsymbol{\mu}^{\prime}$ and $\sigma_{i,j}\leq\sigma_{i,j}^{\prime}$ for all $i,j\in\{1,\dots,n\}$.

Corollary 2.11**.**

For every $(p_{1},\dots,p_{n})$ and all $K\geq 0$ we have

[TABLE]

where $\mathfrak{B}_{n}=\sum_{i=1}^{n-1}\frac{\left|B(q_{i})\right|}{n-1}$ and $\{B(q_{i})\}_{1\leq i\leq n-1}$ has covariance matrix $\boldsymbol{\Sigma}^{p_{1},\dots,p_{n}}_{n}$ under $\mathds{P}^{(p_{1},\dots,p_{n})}$.

Proof.

From Lemma 2.4 and 2.10 we conclude that the distribution of $(B(q_{i}))_{1\leq i\leq n-1}$ is dominated by $\mathds{P}^{(1/2,0,\dots,0,1/2)}$ in the sense that $\mathbb{E}[f(\boldsymbol{X})]\leq\mathbb{E}[f(\boldsymbol{Y})]$ for all directionally convex functions $f$, where $\mathbf{X}\sim\mathcal{N}(\boldsymbol{0},\mathbf{\Sigma}_{n}^{p_{1},\dots,p_{n}})$ and $\mathbf{Y}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}_{n}^{1/2,0,\dots,0,1/2})$. Furthermore we note that, for all $\lambda\in[0,1]$ and all $K\geq 0$,

[TABLE]

As $f(x_{1},\dots,x_{n})=((|x_{1}|+\dots+|x_{n}|)-K)^{+}$ is directionally convex, we conclude from Lemma 2.10 that (2.7) holds for all $K\geq 0$. ∎

3. Proof of Theorem 1.1

Given the results of Section 2, the proof of Theorem 1.1 is mainly concerned with a probabilistic extension of the results obtained in Section 2 to general measures $\mathds{P}\in\mathfrak{P}([0,1])$. This is obtained from weak compactness of the set $\mathfrak{P}([0,1])$ and the following powerful concentration inequality for empirical measures:

Lemma 3.1** (Talagrand, [23], Theorem 1.4).**

Let $\mathbf{X}_{1},\dots,\mathbf{X}_{N}$ be $\mathds{P}^{\otimes}$-independent random variables with values in a measurable space $(S,\mathcal{S})$, let $\mathcal{F}$ be a countable class of measurable functions on $S$ and let

[TABLE]

Then there exists a universal constant $K$ (independent of $\mathbf{X}_{i},N,S$ and $\mathcal{F}$) such that

[TABLE]

We now prove Theorem 1.1 for general measures $\mathds{P}\in\mathfrak{P}([0,1])$, by extending our result for finitely supported measures.

Proof of Theorem 1.1.

For any given measure $\mathds{P}\in\mathfrak{P}([0,1])$ consider a sequence $(\mathds{P}^{n})_{n\in\mathbb{N}}$ of measures supported on the equally spaced grid

[TABLE]

such that we have the weak convergence $\mathds{P}^{n}\Rightarrow\mathds{P}$. Consequently, we have the convergence in law (under $\mathds{P}^{\otimes}$)

[TABLE]

as $n\to\infty$. This suggests that it is sufficient to consider finitely supported measures.

Now fix $\alpha\in[0,1]$. For every $n\in\mathbb{N}$ we denote by $\mathfrak{P}(\mathbf{x}_{n})$ the measures supported on $\mathbf{x}_{n}\subset[0,1]$:

[TABLE]

We conclude by Proposition 2.2 that the quantile $F^{-1}_{\mathfrak{B}_{n}}(\alpha)$ is maximised for measures in $\mathfrak{P}(\mathbf{x}_{n})$ by

[TABLE]

where $\lambda^{n}=\lambda^{n}(\alpha)\in[0,1]$ and $\mathcal{U}(\mathbf{x}_{n})$ is the uniform measure over $\mathbf{x}_{n}$. After taking a subsequence, $\lambda^{n}$ converges to some $\lambda=\lambda(\alpha)\in[0,1]$ and thus

[TABLE]

We wish to show that $\mathds{P}^{\lambda}$ is the measure maximising the desired quantile in $\mathfrak{P}([0,1])$.

Next we need to establish that a general measure can be approximated sufficiently well by a measure on $\mathbf{x}_{n}$, in particular, our intermediate goal is to establish

[TABLE]

where $F^{-1}$ is the quantile under the corresponding measure $\mathds{P}^{\otimes}$. To do this, we follow a similar outline to that in del Barrio, Giné and Matrán [6, Proofs of Theorems 4.2 & 5.1], which in turn heavily relies on Lemma 3.1. Given our observation random variables $X_{i}$, define the quantised versions

[TABLE]

We denote by $\mathds{P}^{\lfloor n\rfloor}$ and $\hat{\mathds{P}}^{\lfloor n\rfloor}_{N}$ the corresponding true and empirical distributions of our quantised observations. In order to obtain (3.2), we need to show weak convergence to zero of

[TABLE]

for $n\to\infty$, uniformly in $N$. Writing

[TABLE]

we are interested in the convergence of $N^{-1/2}|Z_{N}-Z_{N}^{\lfloor n\rfloor}|$. Define the random variables

[TABLE]

Since $L_{1}(\mathbb{R})$ is separable, by an application of the Hahn–Banach theorem there exists a countable subset of the unit ball of $L_{\infty}(\mathbb{R})$, which we denote by $\mathcal{F}$, such that

[TABLE]

By construction,

[TABLE]

For the family of functions $\mathcal{F}$ and random variables $\{\mathbf{X}_{i}\}_{i=1}^{N}$, we next estimate the constants $U$ and $V$ appearing in Lemma 3.1. Clearly $U\leq 1$ by definition of $\mathcal{F}$ and

[TABLE]

Similarly, we compute

[TABLE]

Combining with (3), we have

[TABLE]

We conclude, using Lemma 3.1, that there exists a universal constant $K$ such that

[TABLE]

for all $n\in\mathbb{N}$ and $t>0$. Evaluating (3.4) at $t=\sqrt{N/\log(n)}$ yields

[TABLE]

where the third line follows from the fact that $x\mapsto x\log(1+c/x)$ is increasing for all $c>0$. In particular, the bound in (3) is independent of $N$ and the distribution $\mathds{P}\in\mathfrak{P}([0,1])$, so we have the convergence in probability (under $\mathds{P}^{\otimes}$) of $N^{-1/2}|Z_{N}^{\lfloor n\rfloor}-Z_{N}|\to 0$ as $n\to\infty$, uniformly in $N\in\mathbb{N}$ and $\mathds{P}\in\mathfrak{P}([0,1])$.

We now use this convergence to establish (3.2). Fix $\varepsilon>0,\alpha\in[0,1]$ and take $n\in\mathbb{N}$ large enough that,

[TABLE]

For any $x\in\mathbb{R}$ such that

[TABLE]

by (3) we have

[TABLE]

It follows that

[TABLE]

In particular, this implies

[TABLE]

from which (3.2) follows.

Now that we have the desired uniform approximation result (3.2), we can combine (3) with our result for finitely supported measures (Lemma 2.1), to conclude

[TABLE]

For fixed $\varepsilon>0$ it is thus sufficient to consider $\mathds{P}^{\lambda^{n}}$ instead of $\mathds{P}^{\lambda}$ for large $n\in\mathbb{N}$ in the statement of Theorem 1.1. As the asymptotic relations $\lim_{\alpha\downarrow 0}\lambda^{n}(\alpha)=0$ and $\lim_{\alpha\uparrow 1}\lambda^{n}(\alpha)=1$ follow from Propositions 2.3 and 2.6, this concludes the proof. ∎

4. Numerical results for measures with finite support

We illustrate the implications of Theorem 1.1 for measures with finite support, i.e. Proposition 2.2, by a numerical experiment. We also begin to consider measures in higher dimension, optimising the empirical quantiles associated with some measures in $\mathfrak{P}([0,1]^{2})$.

4.1. Illustration of Proposition 2.2

Here we look at the distribution of the discretised Brownian bridge $(B(q_{i}))_{1\leq i\leq n-1}$ in a ball $\Xi(t)$ as a function of the radius $t$. To this end, we estimate the integral in (2) via Monte Carlo

[TABLE]

where $\textbf{B}^{(i)}$, $i=1,\dots,M$ are i.i.d. draws from $\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma}_{n}^{p_{1},\dots,p_{n}})$. Figure 1 shows the result with a Monte Carlo sample of size $M=10^{5}$ for $n=10$ (left panel) and $n=1000$ (right panel). In both panels we plot the extremal cases from $\mathds{P}_{\delta}=\frac{1}{2}(\delta_{0}+\delta_{1})$ in black and from $\mathds{P}_{\mathcal{U}}=\mathcal{U}(\{x_{1},\dots,x_{n}\})$ in blue. We also include the distribution function from convex combinations $\mathds{P}^{\lambda}=\lambda\mathds{P}_{\delta}+(1-\lambda)\mathds{P}_{\mathcal{U}}$, $\lambda\in(0,1)$ in red and from randomly generated measures in green.

Note that, for $\alpha\approx 0.4$, intermediate distributions with $\lambda\in(0,1)$ yield the most extreme quantiles. This is highlighted in Figure 2 for the case $n=10$, where we plot the $\lambda$ which generates the dominating measure $\mathds{P}^{\lambda}$ as a function of $\alpha\in(0,1)$. (The ‘step function’ behaviour of the graph is due to numerical noise.) Note the behaviour $\lambda(\alpha)\uparrow 1$ as $\alpha\uparrow 1$ and $\lambda(\alpha)\downarrow 0$ as $\alpha\downarrow 0$, that is, dominance by the two extreme cases $\mathds{P}_{\delta}$ and $\mathds{P}_{\mathcal{U}}$ respectively. This is as expected from Theorem 1.1.

4.2. Empirical quantiles for measures on the unit square

Consider measures in $\mathfrak{P}([0,1]^{2})$ with finite support $\textbf{x}\subset[0,1]^{2}$. For simplicity, we take x to be a Cartesian product

[TABLE]

which we index as a $n\times m$ matrix with $(x,y)$-elements

[TABLE]

A measure $\mathds{P}\in\mathfrak{P}([0,1]^{2})$ with support x can then be identified with a $n\times m$ matrix p with elements $p_{ij}\in[0,1]$ summing up to one

[TABLE]

We denote by $\mathfrak{P}(\textbf{x})$ the set of all such probability matrices.

The Wasserstein distance with $\ell_{1}$-norm between $\mathds{P}$ and a measure $\tilde{\mathds{P}}$ with support $\tilde{\textbf{x}}$ of size $\tilde{n}\times\tilde{m}$ and probabilities $\tilde{\textbf{p}}\in\mathfrak{P}(\tilde{\textbf{x}})$ is given by

[TABLE]

where $\Pi(\mathds{P},\tilde{\mathds{P}})$ is the set of couplings in $\mathfrak{P}([0,1]^{4})$ that can be identified with tensors $\boldsymbol{\pi}$ of dimension $n\times m\times\tilde{n}\times\tilde{m}$ such that

[TABLE]

Note that the objective in (4.2) with constraints (4.3) is a linear problem which can be solved numerically by linear programming.

To investigate empirically the distribution of $\mathcal{W}(\mathds{P},\hat{\mathds{P}}_{N})$ we fix a reference measure $\mathds{P}=(\textbf{x},\textbf{p})$ with support being a regular Cartesian product of equidistant points $x_{i},y_{i}=(i-1)/(n-1)$, $i=1,\dots,n$, that is, of size $n\times n$. Figure 3 shows an example measure for $n=3$.

Next, we sample points $X_{1},\dots,X_{N}$ from to $\mathds{P}$ and associate the empirical measure $\hat{\mathds{P}}_{N}=(\textbf{x},\hat{\textbf{p}}_{N})$, where $\hat{\textbf{p}}_{N}$ has elements

[TABLE]

Using this, we can compute a sample point of $\mathcal{W}(\mathds{P},\hat{\mathds{P}}_{N})=:\mathcal{W}(\textbf{p},\hat{\textbf{p}}_{N})$. If we repeat this procedure $M$ times we obtain an $M$-sample of the distance from which we can compute an empirical distribution function $\hat{F}_{\mathcal{W}(\textbf{p},\hat{\textbf{p}}_{N})}$. Figure 4 shows the result with the measure in Figure 3 as reference for $N=100$ and $M=500$ generated outcomes of the Wasserstein distance with $\ell_{1}$-norm. To compute the latter, we use the interior-point method of Andersen and Andersen [2] to solve the linear program.

For a fixed confidence level $\alpha\in[0,1]$, we can now approach optimising the empirical quantile function

[TABLE]

with a numerical optimiser.

Note, however, that the global problem (4.4) is non-trivial. The objective is high dimensional, expensive to evaluate (each sample-point of the distance is computed with a numerical linear-program solver) and with estimation error—in all making a case for noisy, gradient-free optimisation. To this end, we employ Bayesian optimisation, see e.g. Osborne, Garnett and Roberts [16] or Shahriari et al [21] for a review. These methods use Bayesian theory to frame the optimisation as a sequential decision problem through a probabilistic surrogate of the “black-box” objective function. The surrogate — in our case a noisy Gaussian process model (see Williams and Rasmussen [27]) of the empirical quantile function — is used to carefully select a set of unseen arguments as candidates for an optimum, where the selection is with a Bayesian expected-loss criterion. The true objective function is then evaluated at selected locations, where a new potential optimum is recorded. The surrogate model is also updated and re-estimated with acquired data; in our case the Gaussian process posterior and its hyper-parameters. The sequential procedure continues until a stopping criterion kicks in, here a given number of function evaluations.

Figure 5 shows a maximising measure of $\hat{F}^{-1}_{\mathcal{W}(\mathds{P},\hat{\mathds{P}}_{100})}(0.95)$ where the distribution function is estimated from $M=100$ draws of the distance. Note that the quantile is maximised by an extremal measure with non-zero and equal mass at two points only, where the points are placed at two opposite corners. This is the case for high confidence levels: in Figure 6 we include the results for $\alpha\in\{0.75,0.80,0.85,0.95\}$. Note in particular $\alpha=0.85$, where the non-zero point-masses are placed at the other pair of opposite corners.

Similarly, for low confidence levels $\alpha\in\{0.05,0.10\}$ the result is shown in Figure 7. Here the estimated quantiles are maximised by measures with mass at all points of the support. Although not completely uniform, the apparent deviation from such a measure is likely due to estimation- and numerical errors: either in the empirical distribution, or in the solution of the optimisation problem (or in both).

In any case, from the numerical results in $\mathbb{R}^{2}$ and comparison with Theorem 1.1 we derive the following conjecture:

Conjecture 4.1**.**

Let $d\geq 1$, $K\subseteq\mathbb{R}^{d}$ be a convex compact set such that $0<\mu(K)$, where $\mu$ denotes the Lebesgue measure on $\mathbb{R}^{d}$. Furthermore let $x_{0},x_{1}\in\mathbb{R}^{d}$ attain

[TABLE]

Let $\alpha\in[0,1]$ and let $\ \mathcal{U}(K)$ denote the uniform distribution on $K$. Take $\varepsilon>0$. Then there exists $\lambda=\lambda(\alpha,\varepsilon)\in[0,1]$ such that

[TABLE]

asymptotically maximises the normalised quantile $N^{1/2}F^{-1}_{\mathcal{W}(\hat{\mathds{P}}_{N},\mathds{P})}(\alpha)$, i.e. there exists $N_{0}\in\mathbb{N}$ such that for all $N\geq N_{0}$

[TABLE]

Furthermore $\lim_{\alpha\downarrow 0}\lambda=0$ and $\lim_{\alpha\uparrow 1}\lambda=1$.

A proof of such a multidimensional extension of Theorem 1.1 in the spirit of Conjecture 4.1 does not seem straightforward to us, at least given the tools we use: While Section 3 is mainly concerned with controlling the quantisation error and should have a multidimensional analogue, we rely on a representation of $\mathcal{W}(\hat{\mathds{P}}_{N},\mathds{P})$ through the integral over cumulative distribution functions $F_{\hat{\mathds{P}}_{N}}$ and $F_{\mathds{P}}$ in the course of the whole paper. This goes back to the theory of optimal transport, and seems to be a fundamentally one-dimensional result.

For the same reasons our method of proof only covers the case of the $1$-Wasserstein distance, while we expect a similar result to hold for the $p$-Wasserstein distance for $p>1$ also. We refer to Johnson & Samworth [20, Cor. 2.7(b), p.8] for a convergence result corresponding to Lemma 2.1 for $p>1$ for finitely supported measures.

5. Applications & outlook

For $\mathds{P}=(\delta_{0}+\delta_{1})/2$ clearly $\int_{0}^{1}|B(F_{\mathds{P}}(t))|dt=|B(1/2)|$ and in particular

[TABLE]

where $\Phi$ denotes the cumulative distribution function of a standard normal random variable. Thus, intuitively, Theorem 1.1 shows that for any measure $\mathds{P}$, we have the approximate relation

[TABLE]

for large $t\geq 0$ and large $N\in\mathbb{N}$.

On the other hand, for small $t\geq 0$ we are interested in $\mathds{P}=\mathcal{U}([0,1])$ and so we need to compute the distribution of $\int_{0}^{1}|B(F_{\mathds{P}}(t))|dt=\int_{0}^{1}|B(t)|dt$. We refer to Tolmatz [24, 25], where the distribution of this integral has been studied.

These observations allow us to compute an explicit non-parametric confidence region for a sampling distribution $\mathds{P}$, at least in the case $\alpha\approx 1$.

Theorem 5.1**.**

Let $X_{1},\dots,X_{N}$ be i.i.d. samples from an unknown measure $\mathds{P}\in\mathfrak{P}([0,1])$, with empirical distribution $\hat{\mathds{P}}_{N}(\omega)$. For confidence levels $\alpha\to 1$, the random set of measures

[TABLE]

where

[TABLE]

is an asymptotic confidence region for $\mathds{P}$; in the sense that

[TABLE]

where $g:[0,1]\to[0,1]$ is a continuous function with $g(\alpha)\sim 1-\alpha$ as $\alpha\to 1$.

Proof.

Let us first fix $n\in\mathbb{N}$ and let us recall the finitely supported measures $\mathds{P}^{\bar{\lambda}(t,n)}$ from the proof of Theorem 1.1. From Remark 2.7, there is a continuous function $\bar{C}\geq 1$ with $\bar{C}(t)\to 1$ as $t\to\infty$ such that, for all $t\geq 0$ and all $\alpha\in[0,1]$,

[TABLE]

As described above, we also know the distribution of $\mathcal{W}(\hat{\mathds{P}}_{N}(\omega),\mathds{P}^{1})$ under $(\mathds{P}^{1})^{\otimes}$, that is,

[TABLE]

Define

[TABLE]

where $k$ is as in the statement of the theorem. Observe that $g$ satisfies $g(\alpha)\sim 1-\alpha$ as $\alpha\to 1$. For large $N$, we know from (5.1) and (5.2) that $k$ satisfies

[TABLE]

From the approximation results in Section 3, for any $\mathds{P}$, it follows that

[TABLE]

as desired.

∎

It is worth considering how these intervals compare with classical parametric estimation.

Example 5.2*.*

Let $X$ and $X_{1},...,X_{N}$ be iid random variables from some distribution on $[0,1]$. For a function $f\in\mathrm{Lip}_{1}$ (that is $|f(x)-f(y)|\leq|x-y|$), define $Y=f(X)$ (By scaling, similar results hold for $K$-Lipschitz functions for any $K>0$). Observe that $\mathrm{Var}(Y)\leq 1/4$ (from Popoviciu’s inequality for variances), and that this bound is attainable for appropriate choices111In particular, Popoviciu’s inequality is attained when $X$ has distribution of the form $(\delta_{0}+\delta_{1})/2$ and $f(x)=x+y$ for some $y\in\mathbb{R}$. This is the same distribution achieving our extreme quantiles in Theorem 1.1 for large $\alpha$. of $\mathds{P}$ and $f$.

A classical confidence bound for $\mathbb{E}[Y]$ (using this worst-case bound on the variance and the central limit theorem) given observations $X_{1},...,X_{N}$ is

[TABLE]

Using an approach based on the Wasserstein distance $\mathcal{W}$, the Kantorovich–Rubinstein duality formula states that, for any $\mathds{P}^{\prime}\in\mathfrak{P}([0,1])$,

[TABLE]

So, for our specific choice of $f$, we obtain the corresponding bounds

[TABLE]

Theorem 5.1 states that we have approximately $\alpha$-confidence that the true sampling distribution $\mathds{P}$ lies within the region $\hat{C}=\{{\mathds{P}^{\prime}}:\mathcal{W}({\mathds{P}^{\prime}},\hat{\mathds{P}}_{N})\leq k/\sqrt{N}\}$. Substituting in the value of $k$ and evaluating $\mathbb{E}_{\hat{\mathds{P}}_{N}}[f(X)]$, we obtain the confidence bounds

[TABLE]

These non-parametric bounds are identical to those obtained through classical parametric approximations for every choice of $f$. On the one hand, by deriving them using a non-parametric approach, we obtain a bound which is uniform in $f\in\mathrm{Lip}_{1}$, that is, if $\hat{C}_{f}$ denotes the interval given by (5.3), we have

[TABLE]

whereas the classical approach only gives us an interval for a specific choice of $f$. It is somewhat surprising that the interval computed in this non-parametric manner, for a bound which is uniform in $f$, is no wider than that obtained earlier for a specific $f$. In particular, no Bonferroni-type correction of the interval width is required (cf. for example, the classical methods of Dunn [10]).

On the other hand, the Wasserstein distance derived bound is only asymptotic for large $\alpha$ (although our numerical results suggest this approximation is very good for $\alpha>0.7$) and large $N$ (this is also true for those described in the parametric case, as they are based on the central limit theorem). Tighter bounds can be obtained in the parametric case by estimating the variance of $Y$ from observations.

Acknowledgements

Samuel Cohen and Martin Tegnér thank the Oxford–Man Institute for research support. Samuel Cohen acknowledges the support of The Alan Turing Institute under the Engineering and Physical Sciences Research Council grant EP/N510129/1. Johannes Wiesel and Samuel Cohen gratefully acknowledge support from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement no. 335421. Johannes Wiesel furthermore acknowledges support of the German Academic Scholarship foundation and of the Erwin Schrödinger Institute for Mathematics and Physics.