# On the use of the M-quantiles for outlier detection in multivariate data

Sajal<EMAIL_ADDRESS>[ Ram
<EMAIL_ADDRESS>[ A. Alexandre
<EMAIL_ADDRESS>[ Department of Mathematics and
Statistics, Texas Tech University, USA presep=, ]e1,e2,e3

###### Abstract

Defining a successful notion of a multivariate quantile has been an open
problem for more than half a century, motivating a plethora of possible
solutions. Of these, the approach of [8] and [25] leading to M-quantiles, is
very appealing for its mathematical elegance – combining elements of convex
analysis and probability theory. The key idea is the description of a convex
function (the K-function) whose gradient (the K-transform) is in one-to-one
correspondence between all of $\mathbb{R}^{d}$ and the unit ball in
$\mathbb{R}^{d}$. By analogy with the $d=1$ case where the K-transform is a
cumulative distribution function-like object (an M-distribution), the fact
that its inverse is guaranteed to exist lends itself naturally to providing
the basis for the definition of a quantile function for all $d\geq 1$. Over
the past twenty years the resulting M-quantiles have seen applications in a
variety of fields, primarily for the purpose of detecting outliers in
multidimensional spaces. In this article we prove that for odd $d\geq 3$, it
is not the gradient but a poly-Laplacian of the K-function that is (almost
everywhere) proportional to the density function. For $d$ even one cannot
establish a differential equation connecting the K-function with the density.
These results show that usage of the K-transform for outlier detection in
higher odd-dimensions is in principle flawed, as the K-transform does not
originate from inversion of a true M-distribution. We demonstrate these
conclusions in two dimensions through examples from non-standard asymmetric
distributions. Our examples illustrate a feature of the K-transform whereby
regions in the domain with higher density map to larger volumes in the co-
domain, thereby producing a magnification effect that moves inliers closer to
the boundary of the co-domain than outliers. This feature obviously disrupts
any outlier detection mechanism that relies on the inverse K-transform.

multivariate distribution function,

multivariate quantile,

geometric quantile,

outlier detection,

###### keywords:

and and

###### Contents

  1. 1 Introduction
  2. 2 The K-transform: Definition and Outlier Detection
  3. 3 The K-transform in Relation to the Distribution
  4. 4 The zoom-in effect for $\mathbb{R}^{2}$
  5. 5 Conclusion
  6. A Algorithm for Computation of Geometric Quantiles
  7. B Proof of Theorem 2.1
  8. C Proof of Lemma 3.2
  9. D Proof of Lemma 3.3
  10. E Proof of Lemma 3.4
  11. F Proof of Lemma 3.5
  12. G Proof of Theorem 3.6
  13. H Proof of Corollary 1

##  1 Introduction

The subject of outlier detection has a long history in the statistical and
data-mining literature. This harks at the difficulty of the enterprise, which
reflects the fact that a precise definition is ellusive. Thus, while there is
general agreement that an outlier is an observation that deviates so much from
the others as to arouse suspicion, or appears to be inconsistent with the
remainder of the data, there is much debate on how to measure the actual
“deviation” or “inconsistency” [3].

For univariate data that can plausibly be modeled as a realization of a random
sample (independent and identically distributed, or i.i.d.), most detection
methods rely on the quantile function of the data generating mechanism. In the
taxonomy of the field one can distinguish between parametric methods which
classify outliers as those lying on the tails of the distribution, and
nonparametric methods that rely on inter-quantile ranges [47]. For
multivariate data, the inherent difficulty in defining quantiles in
multidimensional spaces has shifted the focus to measuring the distance of
individual points from the center of the data cloud. Thus one finds a
substantial body of literature on such _distance functions_ , or their
inverses, _depth functions_ ; see for example [34] for a survey, and [49] for
formal definitions and theoretical development.

Quantiles are calculated by inverting the cumulative distribution function
(cdf). Mathematically, for any random variable $X$ defined on $\mathbb{R}$,
the most generally accepted definition of the quantile function at probability
level $s$, is the mapping $Q:[0,1]\mapsto\mathbb{R}$, given by

$Q(s)=\inf\\{x:F(x)\geq s\\},$ (1)

where $F(x)=P(X\leq x)$ denotes the cdf of $X$, corresponding to a left-to-
right accumulation of mass. In one dimension, the only other essentially
different mass accumulation would be center-outward (with “center” taken to
be, e.g., the median). We term $\overline{F}(x)=2\,F(x)-1$ the _signed_ cdf.
It is the so-called _M-distribution_ arrived at by [25] when attempting to
extend the notion of a quantile. [28] notes that goodness-of-fit tests like
Kolmogorov-Smirnov based on mass accumulation dictated by $|\overline{F}(x)|$,
which corresponds to a center-outward ordering, are more powerful at detecting
scale differences. Such a center-outward ordering also features prominently in
[23], who define a multivariate cdf (and associated quantile function) based
on the intriguing concept of transportation-based ranks and signs. One may
therefore define an alternate quantile function based on the mapping
$\overline{Q}:[-1,1]\mapsto\mathbb{R}$, given by

$\overline{Q}(s)=\inf\\{x:\overline{F}(x)\geq s\\}.$ (2)

As there is no natural ordering in multidimensional spaces (dimension $d$ in
general), it is difficult to extend the notion of a quantile. Early
contributions to this field of study were made by [22], [5], [36], [4], and
[11], among others. [42] provides a comprehensive overview and an exhaustive
comparison of the different methodologies. Although it is not our intent to
give an updated survey here, we highlight two approaches that serve to give
some idea of the difficulties involved. [4] define an M-quantile through the
notion of an M-estimator and its associated loss and influence functions.
Although this leads to sensible definitions for the median, e.g., under
radially symmetric loss we obtain the geometric median (discussed in detail
below), the quantile extension is cumbersome and suffers from a high
computational burden. Another definition for finite as well as for infinite-
dimensional spaces (e.g., functional data) was proposed by [18]. Their
approach however, which builds on earlier ideas by [26] and others that rely
on directional projections of the probability measure onto the unit sphere,
leads to quantiles that are not affinely equivariant. (Instead of plotting
quantile contours, they proposed computing principal quantile directions in
order to detect outliers.)

A common thread in these various attempts to generalize the notion of a
quantile, is the imposition of the requirement that the specialization to the
univariate case corresponds to the quantile function $Q$ or $\overline{Q}$.
Any such attempt must also by necessity be grounded on the simpler problem of
leading to the special case of a well-defined median. In this regard, it has
long been noted that $Q(0.5)$ minimizes the mean absolute error ([45]),

$Q(0.5)=\underset{c\in\mathbb{R}}{\arg\min}\,\mathbb{E}|X-c|=\underset{c\in\mathbb{R}}{\arg\min}\,\mathbb{E}(|X-c|-|X|),$

where the last expression is to be used in cases when the expectation does not
exist. By extension, the minimizer of the multivariate mean absolute error is
typically defined as the multivariate median:

$\text{Med}_{p}=\arg\min_{s\in\mathbb{R}^{d}}\mathbb{E}(\|X-s\|_{p}-\|X\|_{p}),$
(3)

where $\|\cdot\|_{p}$ denotes the $L^{p}$-norm in $\mathbb{R}^{d}$. In the
literature, one finds both the $L^{2}$-norm (Euclidean distance) and
$L^{1}$-norm (Manhattan distance) used in (3); [37] calls the former the
_spatial median_ , and the latter the _marginal median_ (vector of marginal
medians). The marginal median, first studied by [38] and later investigated
for outlier detection in [7], suffers from lack of invariance with respect to
coordinate rotations, or more precisely, actions of the special orthogonal
group $SO(d)$. Other names in common usage when the 2-norm is employed in (3),
include _geometric median_ and _$L^{1}$ median_ ([44, 30, 39]). It is easily
shown that geometric median is invariant with respect to coordinate rotations
and translations of the origin. Henceforth in this article, we will consider
the $L^{2}$-norm in Equation (3) as we attempt to generalize and clarify
various notions.

Chaudhuri [12] extended the notion of the geometric median and defined a
quantity that he referred to as the _geometric quantile_. For $v\in B_{1}(0)$,
the open unit ball in $\mathbb{R}^{d}$, this is given by

$Q_{G}(v)=\underset{s\in\mathbb{R}^{d}}{\arg\min}\;\mathbb{E}\\{\|X-s\|_{2}+\langle
v,X-s\rangle-\|X\|_{2}-\langle v,X\rangle\\},$ (4)

where $\langle\cdot,\cdot\rangle$ denotes the usual Euclidean inner product.
The justification given for referring to this as a quantile is that $Q_{G}(0)$
corresponds to the geometric median, i.e., (3) with $p=2$. Contemporaneously,
Koltchinskii [25] developed a function (henceforth _K-function_), which in
principle resembles the integral of a multivariate cdf,

$f(s)=\mathbb{E}[\|s-X\|_{2}-\|X\|_{2}].$ (5)

If we now consider the subdifferential function $\partial f(s)$, which we call
the _K-transform_ (defined in (9) with explicit reference to the underlying
measure $\mu$), Koltchinskii’s _M-quantile_ is then the inverse K-transform,
which constitutes his attempt at defining a multivariate quantile. To relate
Chaudhuri’s notion of a quantile to Koltchinskii’s, note that:

$Q_{G}(v)=\underset{s\in\mathbb{R}^{d}}{\arg\min}\;\mathbb{E}\\{\|X-s\|_{2}+\langle
v,X-s\rangle\\}=\underset{s\in\mathbb{R}^{d}}{\arg\max}\;\mathbb{E}\\{\langle
v,s-X\rangle-\|s-X\|_{2}\\}.$ (6)

Now, since $f(s)$ is convex in $s$, we see that the function $Q_{G}(v)$ is
simply the Fenchel conjugate of $f(s)$.

The noteworthy difference between these two approaches is that, whereas [12]
directly defines $Q_{G}(v)$ as a quantile without relating it to the inverse
of a suitably defined cdf, [25] starts with $f(s)$, defines $\partial f(s)$ as
an _M-distribution_ (which coincides with the signed cdf $\overline{F}(s)$ in
the unidimensional case), and uses results from convex analysis to obtain
$Q_{G}(v)=\partial f^{-1}(v)$. An immediate corollary of this fact is that the
geometric quantile uniquely determines the underlying distribution ([25]). [7]
claims that this property may be used for outlier detection, in that values of
$\|Q_{G}(v)\|_{2}$ closer to one are indicative of outliers. However, [13]
notes that this claim may not hold in $\mathbb{R}^{2}$, by demonstrating that
outliers may have a smaller value of $\|Q_{G}(v)\|$ than inliers for non-
standard distributions such as mixtures of Gaussians. However, $d=2$ is a
special case for the K-transform as may be seen in [7, Theorem 3.2], and it is
unclear whether outlier detection may be performed for higher dimensional
data.

The concept of geometric quantiles has been a very active area of research.
Although it is not very useful for outlier detection beyond two dimensions as
we show in this article, it is aesthetically pleasing due to connections with
convex analysis. The mathematical elegance also permits generalizations to
infinite dimensional Banach spaces. With equation (4) as the population
geometric quantile, the empirical geometric quantile for a $d$-dimensional
sample of $n$ data points $\\{x_{1},\ldots,x_{n}\\}$ is

$\widehat{Q}_{G}(v)=\underset{s\in\mathbb{R}^{d}}{\arg\min}\;\frac{1}{n}\,\sum_{i=1}^{n}\\{\|x_{i}-s\|_{2}+\langle
v,x_{i}-s\rangle-\|x_{i}\|_{2}-\langle v,x_{i}\rangle\\}.$ (7)

[12] and [7] have derived a Bahadur-type representation for
$\widehat{Q}_{G}(v)$, and proved Glivenko-Cantelli and Donsker-type results.
[6] have extended these results for data on infinite dimensional Banach
spaces. Although the quantiles in (4) are equivariant under rotations of the
data cloud, they are not equivariant under affine transformations. [8] develop
the notion of an affine equivariant multivariate median for finite dimensional
data, and [7] extend this notion to $L^{p}$ spaces. Computation of $Q_{G}(v)$
has however been a stumbling block beyond the simplest of distributions. In
Appendix A we present an algorithm for this purpose.

Several extensions and applications of these M-quantiles have appeared in the
literature. Their main use has been to effect outlier detection in areas as
varied as manufacturing processes, image processing, economic models, risk
assessment, and biostatistics. With very few exceptions, these are all closely
related to the geometric quantiles of [12] or the inverse K-transform of [25]
(which are essentialy identical as explained earlier). They were proposed as a
means to detect outliers in two-dimensional data in [7] and [10]. However, the
underlying distributions were only chosen to be the Gaussian, and hence this
paradigm may not generalize well. It has also been proposed that regions of
low probability density should be associated with outliers in data ([9, 1]).

Various tools called projection quantiles ([35]), quantile contour plots ([7,
10, 20]), multivariate quantile-quantile (QQ) plots ([7, 31, 16]), and spatial
depth functions ([43]), have been developed based on these M-quantiles. Some
applications of these tools include [13] who use them for outlier detection in
images, [15] for image denoising through spatial depth, and [46] who devised a
calibration methodology to improve the quality of the electrocardiograms. The
concept of M-quantiles has also found uses in risk assessment ([24]), control
charts in manufacturing processes ([29]), climate modeling ([17]), and
development of economic models ([2]).

In this article, we will investigate the connection between the probability
density function (pdf) and the K-function in order to determine the utility of
the K-transform for outlier detection in multivariate data. For $d=1$, note
that the K-transform $\nabla f(s)=f^{\prime}(x)=\overline{F}(x)$, which
reduces to the signed cdf. (This was the original motivation of [25] in
defining $\partial f(s)$ as an M-distribution.) Thus
$f^{\prime\prime}(s)=2\rho(s)$ almost everywhere, where $\rho(\cdot)$ is the
underlying pdf of $X$. As we will show, this is the only case where
differentials of the K-transform lead directly to the pdf. Specifically, and
with regard to continuous measures only for which the subdifferential reduces
to the ordinary differential (or gradient), $\partial f(s)=\nabla f(s)$, this
article makes the following contributions.

  * •

For odd dimensions $d\geq 3$, we show it is in fact the poly-Laplacian of
order $(d+1)/2$ of $f(s)$ that is (almost everywhere) equal to a scalar
multiple of $\rho(s)$. If one is to believe that regions of low probability
density are indicative of outliers and anomalies in data (e.g., [9, 1]), then
this shows that $\nabla f(s)$ does not directly lead to the discovery of
regions of low or high probability density in high odd dimensional spaces.
Thus the M-distribution of [25] cannot be considered as a distribution in any
sense in higher dimensional spaces.

  * •

For dimension $d=2$ we elucidate the reason for the observation of [13]
mentioned above, that values of $\|Q_{G}(v)\|_{2}\approx 1$ are not
necessarily indicative of outliers. It was shown by [25] that the K-transform
is in a one-to-one correspondence between $\mathbb{R}^{d}$ and the open unit
ball in $\mathbb{R}^{d}$. However, as we will show, the K-transform has a
zoom-in effect whereby regions of high density occupy a higher volume in the
co-domain, and conversely, regions of low density occupy a smaller volume. Due
to this effect, the distance of a point from the origin in the co-domain does
not necessarily imply that the point lies in a low density region of the
distribution (as would an outlier).

  * •

Exploiting the connections with convexity theory, we propose an algorithm to
compute the geometric quantile $Q_{G}(v)$. The algorithm, based on the MM
principle [27], allows one to easily verify the numerical results obtained by
[7] and [13], and to compute the inverse of the K-transform, even for highly
nonstandard distributions, and in any finite dimensional vector space.

These results lead to the main conclusion of this article: that the
K-transform cannot theoretically be considered to be a distribution for
dimensions $d\geq 2$, and consequently, that its inverse cannot be viewed as a
quantile function. The remainder of the paper is structured as follows. The
K-transform and its inverse are defined and discussed in Section 2. In Section
3, we present formal results of the points made above concerning the
K-transform and empirical geometric quantile, and discuss their significance.
In Section 4, we present numerical examples to illustrate the zoom-in effect
of the K-transform, and comment on its potential use for outlier detection in
two dimensional data. The proofs of all formal results are relegated to the
appendix.

Before we proceed to the formal definitions, we note that a related approach
of using a convex function to define center-outward distribution and quantiles
has seen attention in the literature lately [14, 23, 19]. This approach uses
the celebrated Brenier-McCann theorem that asserts the existence of a convex
function whose gradient transports a reference distribution absolutely
continuous with respect to the Lebesgue measure on $\mathbb{R}^{d}$ to the
population distribution [32]. The convex function from this method does not
have a formula and therefore cannot be directly analyzed as we have done for
the K-function. However, we note that the center-outward quantile contours are
proved to be connected, closed, nested, with continuous boundaries [23]. We
foresee that the connectedness of the quantile contour would cause
difficulties in detecting outliers for the non-standard distributions from
[48] depicted in Figure 3 and 4 respectively, but we do not investigate this
method any further in this article.

##  2 The K-transform: Definition and Outlier Detection

Let $\left(\mathbb{R}^{d},\mathcal{F},\mu\right)$ be a probability space. For
any $s\in\mathbb{R}^{d}$, the Koltchinskii function (K-function) is defined as
the integral transform of $f(s,x)=\|s-x\|_{2}-\|x\|_{2}$ with respect to
$\mu$:

$f_{\mu}\left(s\right)=\mathbb{E}f(s,x)=\int_{\mathbb{R}^{d}}\\!\left(\|x-s\|_{2}-\|x\|_{2}\right)\,\mu(dx).$
(8)

It is easily verified that $f_{\mu}$ is a convex function. The K-transform
([25]), is defined in as the sub-differential of $f_{\mu}$. For
$s\in\mathbb{R}^{d}$,

$\partial
f_{\mu}\left(s\right)=\int_{\mathbb{R}^{d}}\frac{s-x}{\|s-x\|_{2}}\,\mu(dx).$
(9)

The K-transform is in one-to-one correspondence between $\mathbb{R}^{d}$ and
the open unit ball $B_{1}(0)$. Now consider the conjugate function of
$f_{\mu}$, given by,

$f^{\ast}_{\mu}\left(v\right)=\underset{s\in\mathbb{R}^{d}}{\sup}\left[\langle
s,v\rangle-f_{\mu}(s)\right],\quad v\in B_{1}(0).$ (10)

The subdifferential of $f^{\ast}_{\mu}$ exists ([40]), and is given by
$\partial f^{\ast}_{\mu}=\left(\partial f_{\mu}\right)^{-1}$. The map
$\partial f^{\ast}_{\mu}$ is the inverse K-transform, and is referred to as an
M-quantile ([25]). For dimension $d=1$, the K-transform is the signed cdf
$\partial
f_{\mu}(s)=\mu\left(\left(-\infty,s\right]\right)-\mu([s,\infty):=\overline{F}(s)$.
In this case, the inverse K-transform is interpreted as a quantile function.
Also for $d=1$, there are well-established techniques for outlier detection in
data based on quantile plots, due to the facts that there is a natural
ordering on the real line, and the derivative of the K-transform is the
density.

The inverse K-transform has been proposed as a means to detect outliers in
two-dimensional data (see e.g., [7] and [10]). However, the underlying
distributions were chosen to be the Gaussian, which is one of the simplest of
cases. Figure 1 illustrates the potentially misleading result obtained in [7].
In Figure 1, we have plotted 3000 data points from a bivariate normal
distribution with zero mean, unit standard deviation, and a correlation
coefficient of $0.75$. The contours here are the inverse K-transform of the
contours in Figure 1. [7] proposes that data points falling outside the
inverse K-transform pertaining to contours of the K-transform with a large
radius (say $r\geq 0.9$) may be classified as outliers. In Section 4 we
present numerical simulations with more complex densities that reveals a very
interesting feature of the K-transform that we call the zoom-in effect. Due to
this effect, we cannot interpret the radial distance from the origin in Figure
1 as representing “quantiles”, a view that could lead to potentially erroneous
outlier detection results.

Figure 1: K-transform and its inverse for a bivariate Gaussian with zero mean,
unit standard deviation, and correlation coefficient of $0.75$. The
K-transform is shown in Figure 1, and the inverse K-transform in Figure 1,
along with 3000 data points simulated from the distribution.

Our main result in the next section will show that for odd $d\geq 3$, a poly-
Laplacian of the K-function is proportional to the pdf. The proof of this
result requires the following non-trivial measure theoretic version of the
Leibniz Rule. We state this here as a major result in its own right.

######  Theorem 2.1 (Lebesgue Measure Leibniz Rule).

Suppose $f(s,x)$ is bounded for all $x\in[a,b]$, and is both bounded and
integrable with respect to Lebesgue measure $\lambda$ for all
$s\in[\alpha,\beta]$. For any Lipschitz continuous functions
$g,h:[\alpha,\beta]\rightarrow[a,b]$, let:

$\phi(s)=\int_{h(s)}^{g(s)}\\!f(s,x)\,dx.$

For all $s\in(\alpha,\beta)$ and for all $x\in[a,b]\setminus E_{0},$ where
$\lambda(E_{0})=0$, let $\partial f(s,x)/\partial s$ be defined on
$(\alpha,\beta)\times[a,b]$ and bounded in absolute value. Then, for almost
all $s\in(\alpha,\beta)$, the derivative $d\phi(s)/ds$ exists and can be split
as follows:

$\frac{d}{ds}\phi(s)=f(s,g(s))\,\frac{dg(s)}{ds}-f(s,h(s))\,\frac{dh(s)}{ds}+\int_{[h(s),g(s)]\setminus
E_{0}}\\!\frac{\partial}{\partial s}f(s,x)\;dx.$ (11)

##  3 The K-transform in Relation to the Distribution

In this section, we will present our main result concerning the relationship
of the K-transform to the distribution of the underlying random vector $X$ in
a probability space $\left(\mathbb{R}^{d},\mathcal{F},\mu\right)$. We will
show that for any odd values of the dimension $d$, it is in fact the poly-
Laplacian of the K-function that corresponds to a scalar multiple of the pdf
$\rho(x)$ of $X$.

Let $s\in\mathbb{R}^{d}$ be a given point. Let $\sigma_{k}$ denote the
hypersurface area measure on the $k$-sphere. For example, when $d=2$ we have

$d\sigma_{1}(\varphi_{1})=d\varphi_{1},d\sigma_{2}(\varphi_{1},\varphi_{2})=\sin(\varphi_{2})\,d\varphi_{1}\,d\varphi_{2},$

and in general,

$d\sigma_{k}(\varphi_{1},\varphi_{2},\cdots,\varphi_{k})=\prod_{j=1}^{k}\,\sin^{j-1}(\varphi_{j})\,d\varphi_{j},\qquad
k=1,2,3,\ldots.$

For simplicity, write $\Phi_{k}$ for the angular tuple
$(\varphi_{1},\varphi_{2},\cdots,\varphi_{k})$. The volume element of
$x\in\mathbb{R}^{d}$ represented by the polar coordinates
$\left(r,\varphi_{1},\varphi_{2},\cdots,\varphi_{d-1}\right)$, where $r$ is
the distance of $x$ from $s$, is given by

$dV(x)=r^{d-1}\,dr\,d\sigma_{d-1}(\Phi_{d-1}).$

An essentially bounded function
$\rho:\mathbb{R}^{d}\rightarrow\mathbb{R}_{+}$, satisfying

$\displaystyle\int_{\mathbb{R}^{d}}\\!\rho(x)\,dV(x)=1,$

is said to be a _bounded density_ (bounded pdf). On the sigma-algebra of Borel
subsets of $\mathbb{R}^{d}$, we may use $\rho$ to define a Radon measure as
arising from a positive linear form on $C_{c}\left({\mathbb{R}^{d}}\right)$,
the set of continuous functions with compact support on $\mathbb{R}^{d}$.
Alternatively, one starts with an absolutely continuous Radon measure with
respect to the Lebesgue measure, and takes its Radon-Nikodym derivative to be
the density $\rho$.

######  Remark 3.1.

Note that if $\mu(dx)=\rho(x)dV(x)$, then the K-transform is simply the
gradient of the K-function, so that $\partial f_{\mu}(s)=\nabla f_{\mu}(s)$ in
(9).

The lemma below is a building block on which the theory is based. It shows
that the expectation of certain singular kernels is in fact well defined. The
only mild requirement is that the probability density should be bounded.

######  Lemma 3.2.

Suppose $\rho(x)$ is a bounded pdf for the random vector $X$ on
$\mathbb{R}^{d}$, with $d\geq 3$. Then, for any $s\in\mathbb{R}^{d}$, the
function defined as

$G(s)=\int_{\mathbb{R}^{d}}\\!\|s-x\|_{2}^{-k}\,\rho(x)\,dV(x),$ (12)

exists for any $k=0,\ldots,d-1.$

The second lemma shows that the gradient of $G(s)$ defined above exists for
$0\leq k\leq d-2.$

######  Lemma 3.3.

Suppose $\rho(x)$ is a bounded pdf on $\mathbb{R}^{d}$, for $d\geq 3$. Then,
for any $s\in\mathbb{R}^{d}$ and $k=0,\ldots,d-2$, the Frechet derivative of
$G(s)$ defined in (12) with respect to $s$, exists and is given by

$\nabla_{s}\,G(s)=-k\,\int_{\mathbb{R}^{d}}\\!\frac{(s-x)}{\|s-x\|_{2}^{k+2}}\,\rho(x)\,dV(x).$
(13)

The third lemma shows that the Hessian of $G(s)$ exists for $0\leq k\leq d-3.$
The Hessian computation is important as the Laplacian may be computed as the
trace of the Hessian matrix.

######  Lemma 3.4.

Suppose $\rho(x)$ is a bounded pdf on $\mathbb{R}^{d}$ for $d\geq 3$. Then,
for any $s\in\mathbb{R}^{d}$ and $k=0,\ldots,d-3$, the Hessian of $G(s)$
defined in (12) with respect to $s$, is given by:

$D_{s}G(s)=-k\int_{\mathbb{R}^{d}}\left(\frac{I_{d}}{\|s-x\|_{2}^{k+2}}-\left(k+2\right)\frac{\left(s-x\right)\left(s-x\right)^{t}}{\|s-x\|_{2}^{k+4}}\right)\rho\left(x\right)dV\left(x\right).$
(14)

We are now ready to produce results for the K-transform. The lemma below
considers the univariate case ($d=1$); the special case of odd dimension that
is omitted in the statement of Theorem 3.6. The result is proved in [25] using
convex analysis. Here we use real analysis techniques, because for dimensions
greater than one, they allow one to obtain expressions for the poly-Laplacian
of the K-function.

######  Lemma 3.5.

Consider the probability space $\left(\mathbb{R},\mathcal{F},\mu\right)$,
where the measure $\mu$ of the random element $X$ has a bounded density,
$\rho(x)$, on $\mathbb{R}$. Consider the $d=1$ version of the K-function
defined in (8):

$f_{\mu}(s)=\int_{\mathbb{R}}(|s-x|-|x|)\rho(x)dx.$

Then, for almost every $s\in\mathbb{R}$, the second derivative of $f_{\mu}(s)$
exists, and is given by,

$f_{\mu}^{\prime\prime}(s)=2\rho(s).$ (15)

With the above lemmas in place, we can now state the main theoretical result
of this article. Taken together, we show that for odd dimensions $d\geq 3$,
the poly-Laplacian of order $(d+1)/2$ of the K-function is equal to a scalar
multiple of the density (almost everywhere), with the scalar depending only on
$d$. As regions of low probability density are commonly associated with
outliers and anomalies in data (see [9, 1]), this result shows that, contrary
to emerging belief, it is not the K-transform $\nabla f_{\mu}(s)$ that is
directly related to the density, but the poly-Laplacian
$\Delta^{(d+1)/2}f_{\mu}(s)$ instead. Hence, it is dubious to use the
K-transform to compute “quantiles” in odd dimensions $d\geq 3$; a situation
that is completely different from the $d=1$ case.

Using Lemma 3.4, we can obtain a different type of result for even dimensions.
However, we can only establish a lower bound for $\|\nabla\,f_{\mu}(s)\|_{2}$
for $d=2$, and an upper bound for $\|\nabla\,\Delta^{(d-2)/2}f_{\mu}(s)\|_{2}$
for $d\geq 4$, in terms of $\rho(s)$. Unfortunately these results do not
appear to have as clear an interpretation as those for odd dimensions.

######  Theorem 3.6.

Consider the probability space $\left(\mathbb{R}^{d},\mathcal{F},\mu\right)$,
where the probability measure $\mu$ has a bounded density $\rho(x)$ in
$\mathbb{R}^{d}$, and whose K-function $f_{\mu}(s)$ is as defined in (8). Let
$\Delta^{j}f_{\mu}(s)$ be the poly-Laplacian of order $j$ of $f_{\mu}(s)$,
defined as:

$\Delta^{j}f_{\mu}(s)=(-1)^{j+1}(d-1)(d-3)\dots(d-2j+1)\,\mathbb{E}\left(\frac{1}{\|s-X\|_{2}^{2j-1}}\right).$
(16)

Then, for any $s\in\mathbb{R}^{d}$ and odd values of the dimension $d\geq 3$,
$\Delta^{j}f_{\mu}(s)$ exists for all $j=1,\ldots,N$, where $N=(d-1)/2$.

######  corollary 1.

For any odd $d\geq 3$, and for almost every $s\in\mathbb{R}^{d}$, we have
that:

$\Delta^{N+1}f_{\mu}(s)=\alpha\,\rho(s),$ (17)

where
$\alpha=(-1)^{N-1}\,(d-1)\,\cdots\,(2)\,d\,(d-2)\pi^{d/2}/\Gamma\left(\frac{d}{2}+1\right)$,
and $\Gamma$ is the gamma function.

Note that the first few poly-Laplacians are (assuming odd $d>5$):

$\Delta_{s}f_{\mu}(s)=(d-1)\,\mathbb{E}\left(\frac{1}{\|s-X\|_{2}}\right),\qquad\Delta_{s}^{2}f_{\mu}(s)=(-1)\,(d-1)\,(d-3)\,\mathbb{E}\left(\frac{1}{\|s-X\|_{2}^{3}}\right),$

and

$\Delta_{s}^{3}f_{\mu}(s)=(-1)^{2}\,(d-1)\,(d-3)\,(d-5)\,\mathbb{E}\left(\frac{1}{\|s-X\|_{2}^{5}}\right).$

The final polyharmonic is:

$\Delta_{s}^{N}f_{\mu}(s)=(-1)^{N-1}\,(d-1)\,\cdots\,(2)\,\mathbb{E}\left(\frac{1}{\|s-X\|_{2}^{d-2}}\right).$
(18)

##  4 The zoom-in effect for $\mathbb{R}^{2}$

The inverse K-transform is useful for detecting outliers in simple
distributions on $\mathbb{R}^{2}$, for instance the Gaussian, as demonstrated
by [7]. However, [13] cautioned that this claim is incorrect for non-standard
distributions in $\mathbb{R}^{2}$ such as a mixture of Gaussians. In this
section, we demonstrate why this happens by highlighting a special feature of
the K-transform that we call the “zoom-in” effect.

Consider the unit disk in $\mathbb{R}^{2}$ centered at the origin with a
uniform density of $\pi^{-1}$. By symmetry of the disk, the K-transform at the
point $s=(s_{1},s_{2})$ is given by:

$\nabla
f_{\mu}(s_{1},s_{2})=\int_{B_{1}(0)}\\!\frac{(s_{1},s_{2})-(x_{1},x_{2})}{\|(s_{1},s_{2})-(x_{1},x_{2})\|_{2}}\,\mu(dx).$
(19)

As discusssed in Section 2, the K-transform will map all of $\mathbb{R}^{2}$
onto the unit disk in the codomain. Figures 2, 3, and 4 depict the zoom-in
effect on the banana-shaped [23], spiral-shaped [21, 48], and square-shaped
[21, 48] distributions. Note that the _zoom-in_ effect refers to the tendency
of the K-transform to magnify regions of the support, so that they occupy a
disproportionately larger area of the unit disk.

The banana-shaped distribution is a mixture of three Gaussian distributions
with $n=20,000$ sample points, given by,

$\frac{3}{8}N(-3\tau_{1},\Sigma_{1})+\frac{3}{8}N(3\tau_{1},\Sigma_{2})+\frac{1}{4}N(-\frac{5}{2}\tau_{2},\Sigma_{3}),$

where $\tau_{1}=\begin{pmatrix}1\\\ 0\end{pmatrix}$,
$\tau_{2}=\begin{pmatrix}0\\\ 1\end{pmatrix}$,
$\Sigma_{1}=\begin{pmatrix}5&-4\\\ -4&5\end{pmatrix}$,
$\Sigma_{2}=\begin{pmatrix}5&4\\\ 4&5\end{pmatrix}$, and,
$\Sigma_{3}=\begin{pmatrix}4&0\\\ 0&1\end{pmatrix}$. The contour lines in
Figure 2, 3, and 4 are the inverse K-transform of the contours with radial
distances of 0.5 (blue), 0.75 (red), and 0.90 (maroon) in the codomain. We
have used our developed MM-type algorithm to compute the inverse K-transform
of the contours that are in the codomain. [7, 20] described the inverse
K-contour lines as quantile contours and [7] classified points that are
outside of radial distance 0.9 as outliers. From these figures, we can observe
that because of the zoom-in effect, inlier points become close to the boundary
of the contour lines with a radial distance of 0.9, and outlier points get
close to the inliers.

Figure 2: Banana-shaped distribution with $n=20,000$ sample points [23].
Concentric circles of radius 0.50 (blue), 0.75 (red), and 0.90 (maroon) are
shown in Figure 2. Their corresponding K-transforms are in Figure 2. This
example illustrates the zoom-in effect of the K-transform. Points outside the
support of the distribution in the domain occupy a small region close to the
boundary of the unit circle in the co-domain. Hence, $1-\|\nabla f\|_{2}$
cannot be interpreted as a depth function.

Figure 3: Spiral-shaped distribution with $n=2259$ sample points [48]. Figure
3 shows the K-transform and 3 shows the inverse K-transform. This example
illustrates zoom-in effect and the low power of the K-transform method. Points
in the empty spaces are within the $r=0.5$ M-quantile while points in the
outer spiral of the data are outside the circle of radius $0.9$ in the co-
domain.

Figure 4: This example illustrates zoom-in effect and the low power of the
K-transform method for a non-standard square-shaped distribution with $n=1242$
sample points [48]. Figure 4 shows the K-transform and 4 shows the inverse
K-transform. Points in the empty spaces within the $r=0.5$ M-quantile should
be outliers but are not classified as such by the method.

##  5 Conclusion

Regions of low probability density are commonly associated with outliers and
anomalies in data. In this article, we have shown that for multivariate data
on spaces with odd dimension $d$ greater than or equal to three, it is not the
K-transform $\nabla f_{\mu}(s)$ that is directly proportional to the
probability density, but the poly-Laplacian $\Delta^{(d+1)/2}\,f_{\mu}(s)$
instead. This leads to the natural conclusion that the K-transform is not a
useful object from which to attempt to compute “quantiles” in higher
dimensions; a situation that is very different from the one-dimensional case.
For the lowest even dimensional case of $d=2$, we elucidate a zoom-in feature
of the K-transform whereby regions of higher probability density are magnified
when mapped to the co-domain. Due to this effect, it may not be appropriate to
think of contours of equal radial distance from the origin in the co-domain as
corresponding to quantiles. In summary, the established lore of interpreting
the K-transform as a distribution and its inverse as a quantile function,
applies only to the one-dimensional case.

##  Appendix A Algorithm for Computation of Geometric Quantiles

In the following steps, we describe an MM type algorithm ([27]) for the
computation of the geometric quantile in (6).

Step 0:

Randomly choose an initial value of $s$, say $s_{0}$.

Step 1:

Specify the surrogate objective function to be

$\underset{s\in\mathbb{R}^{d}}{\arg\min}\;\frac{1}{2}\mathbb{E}\left(\frac{\|s-X\|_{2}^{2}}{\|s_{0}-X\|_{2}}+\|s_{0}-X\|_{2}-2\|X\|_{2}\right)-\langle
s,v\rangle,v\in B_{1}(0).$ (20)

Step 2:

Compute the gradient of (20) as,

$v=\mathbb{E}\left(\frac{s-X}{\|s_{0}-X\|_{2}}\right).$ (21)

Step 3:

For $k\geq 1$, iteratively update the value of $s$ until convergence is
reached.

$s_{k}=\frac{v+\mathbb{E}\left(X\|s_{k-1}-X\|_{2}^{-1}\right)}{\mathbb{E}\left(\|s_{k-1}-X\|_{2}^{-1}\right)}.$
(22)

##  Appendix B Proof of Theorem 2.1

###### Proof.

Suppose $dg(t)/dt$ and $dh(t)/dt$ exist for all $t\in(\alpha,\beta),$ and let
$\\{t+t_{k}\\}_{k=1}^{\infty}$ be any sequence contained in $(\alpha,\beta)$
converging to $t.$ We then have:

$\displaystyle\frac{\phi(t+t_{k})-\phi(t)}{t_{k}}$
$\displaystyle=\frac{1}{t_{k}}\left[\int_{h(t+t_{k})}^{g(t+t_{k})}\\!\\!\\!f(x,t+t_{k})\,dx-\int_{h(t)}^{g(t)}\\!\\!\\!f(x,t)\,dx\right]$
$\displaystyle=\frac{1}{t_{k}}\left[\int_{h(t+t_{k})}^{g(t+t_{k})}\\!\\!\\!f(x,t+t_{k})\,dx-\int_{h(t+t_{k})}^{g(t)}\\!\\!\\!f(x,t+t_{k})\,dx\right]+$
$\displaystyle\frac{1}{t_{k}}\left[\int_{h(t+t_{k})}^{g(t)}\\!\\!\\!f(x,t+t_{k})\,dx-\int_{h(t)}^{g(t)}\\!\\!\\!f(x,t+t_{k})\,dx\right]+$
$\displaystyle\frac{1}{t_{k}}\left[\int_{h(t)}^{g(t)}\\!\\!\\!f(x,t+t_{k})\,dx-\int_{h(t)}^{g(t)}\\!\\!\\!f(x,t)\,dx\right]$
$\displaystyle=\frac{1}{t_{k}}\left[\int_{g(t)}^{g(t+t_{k})}\\!\\!\\!f(x,t+t_{k})\,dx-\int_{h(t)}^{h(t+t_{k})}\\!\\!\\!f(x,t+t_{k})\,dx\right]$
$\displaystyle+\frac{1}{t_{k}}\left[\int_{h(t)}^{g(t)}\\!\\!\\![f(x,t+t_{k})\,dx-f(x,t)]\,dx\right].$
(23)

Let $|f(t,x)|\leq M$ on $[\alpha,\beta]\times[a,b]$ and $|\partial
f(t,x)/\partial t|\leq K$ on $(\alpha,\beta)\times(a,b)$. We note that for
$\varepsilon>0$ , there exists $\delta>0$ such that for all $|t_{k}|<\delta$
and $x\in[a,b]$, $|f(x,t+t_{k})-f(x,t)|<K\,|t-t_{k}|$. In addition, since $g$
is Lipschitz continuous, we have that

$\displaystyle\Big{|}\frac{1}{t_{k}}\,\int_{g(t)}^{g(t+t_{k})}\\!\\!\\!\left(f(x,t+t_{k})-f(x,t)\right)\,dx\Big{|}$
$\displaystyle\leq\Big{|}\frac{1}{t_{k}}\,\int_{g(t)}^{g(t+t_{k})}\\!\\!\\!|f(x,t+t_{k})-f(x,t)|\,dx\Big{|}$
$\displaystyle\leq\Big{|}\frac{K}{t_{k}}\,\int_{g(t)}^{g(t+t_{k})}\\!\\!\\!|t+t_{k}-t|\,dx\Big{|}$
$\displaystyle\leq K\,|g(t+t_{k})-g(t)|$ $\displaystyle\leq
K\,\|g\|_{\infty}\,|t_{k}|.$

Therefore,

$\lim_{t_{k}\rightarrow
0}\,\Big{|}\frac{1}{t_{k}}\,\int_{g(t)}^{g(t+t_{k})}\\!\\!\\![f(x,t+t_{k})-f(x,t)]\,dx\Big{|}=0.$
(24)

Now, by the change of variable formula (see e.g., [33, Sec. 38.3]):

$\displaystyle\frac{1}{t_{k}}\int_{g(t)}^{g(t+t_{k})}\\!\\!\\!\\!f(x,t)\,dx$
$\displaystyle=\frac{1}{t_{k}}\int_{t}^{t+t_{k}}\\!\\!\\!\\!f(g(\tau),t)\,g^{\prime}(\tau)\,d\tau.$
(25)

As this integrand is integrable, we then have by the Lebesgue Differentiation
Theorem that for almost every $t\in[\alpha,\beta]$,

$\lim_{k}\,\frac{1}{t_{k}}\int_{t}^{t+t_{k}}\\!\\!\\!f(g(\tau),t)\,g^{\prime}(\tau)\,d\tau=f(g(t),t)\,g^{\prime}(t).$
(26)

By (24) and (26),

$\displaystyle\lim_{k}\,\frac{1}{t_{k}}\,\int_{g(t)}^{g(t+t_{k})}\\!\\!\\!f(x,t+t_{k})\,dx$
$\displaystyle=\lim_{k}\,\frac{1}{t_{k}}\,\int_{g(t)}^{g(t+t_{k})}\\!\\!\\![f(x,t+t_{k})-f(x,t)]\,dx$
$\displaystyle\quad-\lim_{k}\,\frac{1}{t_{k}}\,\int_{g(t)}^{g(t+t_{k})}\\!\\!\\!f(x,t)\,dx$
$\displaystyle=f(g(t),t)\,g^{\prime}(t).$ (27)

By the same approach, the second integral in (23) can be shown to be for
almost every $t\in[\alpha,\beta]$

$\lim_{k}\,\frac{1}{t_{k}}\int_{t}^{t+t_{k}}\\!\\!\\!f(h(\tau),t+t_{k})\,h^{\prime}(\tau)\,d\tau=f(h(t),t)\,h^{\prime}(t).$
(28)

Consider now the final term in (23). Let $x\notin E_{0}$. By the Mean Value
Theorem applied to the interval $[t,t+t_{k}]$ or $[t+t_{k},t]$, depending on
the sign of $t_{k}$, we then have that:

$\Big{|}\frac{f(x,t+t_{k})\,dx-f(x,t)}{t_{k}}\Big{|}=\Big{|}\frac{\partial
f}{\partial t}(\bar{t},x)\Big{|}\leq K,$

for some $\bar{t}$ in the segment joining $t$ and $t+t_{k}$. Now, invoking the
(Lebesgue) Dominated Convergence Theorem, we obtain:

$\displaystyle\lim_{k}\,\int_{h(t)}^{g(t)}\\!\\!\\!\frac{f(x,t+t_{k})\,dx-f(x,t)}{t_{k}}\,dx$
$\displaystyle=$
$\displaystyle\int_{h(t)}^{g(t)}\\!\\!\\!\lim_{k}\,\frac{f(x,t+t_{k})\,dx-f(x,t)}{t_{k}}\,dx$
(29) $\displaystyle=$ $\displaystyle\int_{[h(t),g(t)]\setminus
E_{0}}\\!\\!\\!\frac{\partial f}{\partial t}(\bar{t},x)\,dt.$

Collecting (27) – (29), we obtain finally that for almost every
$t\in(\alpha,\beta)$:

$\displaystyle\frac{d}{dt}\phi(t)$ $\displaystyle=$
$\displaystyle\lim_{k}\,\frac{1}{t_{k}}\left[\int_{h(t+t_{k})}^{g(t+t_{k})}\\!\\!\\!f(x,t+t_{k})\,dx-\int_{h(t)}^{g(t)}\\!\\!\\!f(x,t)\,dx\right]$
(30) $\displaystyle=$ $\displaystyle
f(t,g(t))\,\frac{dg(t)}{dt}-f(t,h(t))\,\frac{dh(t)}{dt}+\int_{[h(t),g(t)]\setminus
E_{0}}\\!\frac{\partial}{\partial t}f(t,x)\;dx.$

∎

##  Appendix C Proof of Lemma 3.2

###### Proof.

Set $z=x-s$. Denote $r=\|z\|$. We will use the volume form for
$\mathbb{R}^{d}$ in the coordinates $z$.

$\displaystyle\int_{\mathbb{R}^{d}}\\!\|s-x\|_{2}^{-k}\,\rho(x)\,dV(x)$
$\displaystyle=\int_{B_{1}(0)}\\!\|z\|_{2}^{-k}\,\rho(z+s)\,dV(z)+\int_{B_{1}^{c}(0)}\\!\|z\|_{2}^{-k}\,\rho(z+s)\,dV(z)$
$\displaystyle=\int_{B_{1}(0)}\\!\\!\\!\\!\\!\\!\\!\\!\\!\rho(z+s)\,r^{d-1-k}\,dr\,d\sigma_{d-1}(\Phi_{d-1})+\int_{\mathbb{R}^{d}\setminus
B_{1}(0)}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\rho(z+s)\,r^{d-1-k}\,dr\,d\sigma_{d-1}(\Phi_{d-1}).$

As $k\geq 0$,

$r^{d-1-k}\leq\left\\{\begin{array}[]{rl}1;&r\leq 1\\\
r^{d-1};&r>1.\end{array}\right.$

Hence,

$\displaystyle\int_{\mathbb{R}^{d}}\\!\|s-x\|_{2}^{-k}\,\rho(x)\,dV(x)$
$\displaystyle\leq\int_{B_{1}(0)}\\!\\!\\!\\!\rho(z+s)\,dr\,d\sigma_{d-1}(\Phi_{d-1})+\int_{\mathbb{R}^{d}\setminus
B_{1}(0)}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\rho(z+s)\,r^{d-1}\,dr\,d\sigma_{d-1}(\Phi_{d-1})$
$\displaystyle\leq\|\rho\|_{\infty}\,\int_{B_{1}(0)}\\!\\!\\!\,dr\,d\sigma_{d-1}(\Phi_{d-1})+\int_{\mathbb{R}^{d}}\\!\\!\rho(z+s)\,r^{d-1}\,dr\,d\sigma_{d-1}(\Phi_{d-1})$
$\displaystyle<C+1$ $\displaystyle<\infty.$

∎

##  Appendix D Proof of Lemma 3.3

###### Proof.

Write $\|\cdot\|$ to mean the usual Euclidean norm $\|\cdot\|_{2}$ in order to
simplify the notation in the rest of the proof. There is a need for subtlety
because the integrands of $G$ and its proposed gradient are undefined for
$x=s$. We need to excise a small region around the point $s$ in the variable
$x$ while dealing with the integrals. To this end, let $R>0$. Since the
integrand of $G(s)$ given by $g(s,x)=\|s-x\|^{-k}$ is a differentiable and
bounded function on $\mathbb{R}^{d}\setminus B_{R}(s)$, then for
$h\in\mathbb{R}^{d}$ satisfying $0<\|h\|<R/2$, we have by the Mean Value
Theorem that for all $x\in\mathbb{R}^{d}\setminus B_{R}(s)$,

$g(s+h,x)-g(s,x)\;=\;\int_{0}^{1}\\!\nabla_{s}g(s+c\,h,x)\,dc\,\cdot\,h.$ (31)

where for $a,b\in\mathbb{R}^{d}$, $a\cdot b=\langle a,b\rangle$ denotes the
usual vector dot product. Hence,

$\displaystyle\frac{|g(s+h,x)-g(s,x)-\nabla_{s}g(s,x)\cdot h|}{\|h\|}$
$\displaystyle\leq\left\|\int_{0}^{1}\\!\nabla_{s}g(s+c\,h,x)\,dc\right\|+\|\nabla_{s}g(s,x)\|$
$\displaystyle=\left\|-k\,\int_{0}^{1}\\!\frac{s+c\,h-x}{\|s+c\,h-x\|^{k+2}}\,dc\right\|+k\,\left\|\frac{s-x}{\|s-x\|^{k+2}}\right\|$
$\displaystyle\leq
k\,\int_{0}^{1}\\!\|s+c\,h-x\|^{-k-1}\,dc+k\,\|s-x\|^{-k-1}.$

Integrating over $\mathbb{R}^{d}\setminus B_{R}(s)$, we obtain,

$\displaystyle\int_{\mathbb{R}^{d}\setminus
B_{R}(s)}\\!\frac{|g(s+h,x)-g(s,x)-\nabla_{s}g(s,x)\cdot
h|}{\|h\|}\,\rho(x)\,dV(x)$ (32)
$\displaystyle\leq\int_{\mathbb{R}^{d}\setminus
B_{R}(s)}\\!k\,\left(\int_{0}^{1}\\!\|s+c\,h-x\|^{-k-1}\,dc+\,\|s-x\|^{-k-1}\right)\;\rho(x)\,dV(x)$
(33)
$\displaystyle\leq\int_{\mathbb{R}^{d}}\\!k\,\left(\int_{0}^{1}\\!\|s+c\,h-x\|^{-k-1}\,dc+\,\|s-x\|^{-k-1}\right)\;\rho(x)\,dV(x).$
(34)

By Lemma 3.2, we have that for all $c\in[0,1]$,

$\int_{\mathbb{R}^{d}}\\!k\,\|s+c\,h-x\|^{-k-1}\,\rho(x)\,dV(x)<\infty.$

Defining

$\theta=\max_{c\in[0,1]}\,\int_{\mathbb{R}^{d}}\\!k\,\|s+c\,h-x\|^{-k-1}\,\rho(x)\,dV(x),$

we see that $\theta$ is a real value (the maximum of a continuous function on
a compact interval is achieved at a point in the interval). Hence we obtain
that:

$\int_{0}^{1}\int_{\mathbb{R}^{d}}\\!k\,\left(\|s+c\,h-x\|^{-k-1}+\,\|s-x\|^{-k-1}\right)\,\rho(x)\,dV(x)\,dc\leq
2\,\theta.$

Now, by (34) and Tonelli, note that

$\displaystyle\int_{\mathbb{R}^{d}\setminus
B_{R}(s)}\\!\frac{|g(s+h,x)-g(s,x)-\nabla_{s}g(s,x)\cdot
h|}{\|h\|}\,\rho(x)\,dV(x)$ $\displaystyle\leq 2\,\theta.$

By the Mean Value Theorem applied to (31), and for each $s$, we have the
following pointwise convergence as $h\rightarrow 0$:

$\displaystyle\lim_{\|h\|\rightarrow
0}\frac{|g(s+h,x)-g(s,x)-\nabla_{s}g(s,x)\cdot h|}{\|h\|}$
$\displaystyle\leq\lim_{\|h\|\rightarrow
0}\|\nabla_{s}g(s+c\,h,x)-\nabla_{s}g(s,x)\|=0.$

Therefore, by the Dominated Convergence Theorem,

$\displaystyle\lim_{h\rightarrow
0}\,\frac{1}{\|h\|}\,\int_{\mathbb{R}^{d}\setminus
B_{R}(s)}\\!\left(g(s+h,x)-g(s,x)-\nabla_{s}g(s,x)\cdot
h\right)\,\rho(x)\,dV(x)$ $\displaystyle=0.$

By the definition of the Frechet derivative, we have:

$\displaystyle\nabla_{s}\int_{\mathbb{R}^{d}\setminus
B_{R}(s)}\\!g(s,x)\,\,\rho(x)\,dV(x)$
$\displaystyle=\int_{\mathbb{R}^{d}\setminus
B_{R}(s)}\\!\nabla_{s}g(s,x)\,\rho(x)\,dV(x)$
$\displaystyle=-k\,\int_{\mathbb{R}^{d}\setminus
B_{R}(s)}\\!\frac{(s-x)}{\|s-x\|^{k+2}}\,\rho(x)\,dV(x).$

Noting that this result does not depend on $R>0$, let $R\rightarrow 0$ to
obtain:

$\lim_{R\rightarrow 0}\,\nabla_{s}\int_{\mathbb{R}^{d}\setminus
B_{R}(s)}\\!g(s,x)\,\rho(x)\,dV(x)=-k\,\lim_{R\rightarrow
0}\,\int_{\mathbb{R}^{d}\setminus
B_{R}(s)}\\!\frac{(s-x)}{\|s-x\|^{k+2}}\,\rho(x)\,dV(x).$

Due to Lemma 3.2, the integral on the right hand side converges absolutely and
uniformly as $R\rightarrow 0$, and for any $0\leq k\leq d-2$. Due to the same
lemma, the integral on the left hand side also converges for each $s$ to
$\int_{\mathbb{R}^{d}}\\!g(s,x)\,\rho(x)\,dV(x)$. Therefore, by a standard
result in analysis (e.g., Theorem 7.17 in [41]), the above integrals are equal
over all of $\mathbb{R}^{d}$:

$\nabla_{s}\int_{\mathbb{R}^{d}}\\!g(s,x)\,\rho(x)\,dV(x)=-k\,\int_{\mathbb{R}^{d}}\\!\frac{(s-x)}{\|s-x\|^{k+2}}\,\rho(x)\,dV(x).$

∎

(A simple modification of the proof also works for integrands the type
$\phi(g(s,x))$ where $\phi$ is a differentiable function.)

##  Appendix E Proof of Lemma 3.4

###### Proof.

Once again write $\|\cdot\|$ to mean the usual Euclidean norm $\|\cdot\|_{2}$.
In Lemma 3.3, we showed that under the given conditions on $\rho$, the Frechet
derivative of $G$ with respect to $s$ is given by,

$\nabla_{s}G(s)=-k\int_{\mathbb{R}^{d}}\frac{s-x}{\|s-x\|^{k+2}}\rho(x)dV(x).$
(35)

Below we repeat the steps of the proof of Lemma 3.4, which hold even in the
presence of the additional complication that the integrand in (35) is a vector
function of $s$. Let $R>0$ and $g_{1}(s,x)=-k(s-x)\|s-x\|^{-(k+2)}$, which is
bounded and differentiable on $\mathbb{R}^{d}\backslash B_{R}(s)$. Now, for
$h\in\mathbb{R}^{d}$ satisfying $0<\|h\|<R/2$, an application of the Mean
Value Theorem yields that for all $x\in\mathbb{R}^{d}\backslash B_{R}(s)$,

$g_{1}(s+h,x)-g_{1}(s,x)=\left[\int_{0}^{1}D_{s}g_{1}\left(s+ch,x\right)dc\right]\cdot
h$ (36)

Noting that,

$\begin{split}&\Bigg{\|}\left[-k\int_{0}^{1}\left(\frac{I_{d}}{\|s+ch-x\|^{k+2}}-\frac{(k+2)\left(s+ch-x\right)\left(s+ch-x\right)^{t}}{\|s+ch-x\|^{k+4}}\right)dc\right]\cdot
h\Bigg{\|}\\\
&\leq\int_{0}^{1}k\left(\frac{\sqrt{d}}{\|s+ch-x\|^{k+2}}+(k+2)\frac{\|s+ch-x\|^{2}}{\|s+ch-x\|^{k+4}}\right)dc\\\
&\leq k(\sqrt{d}+k+2)\int_{0}^{1}\|s+ch-x\|^{-k-2}dc.\end{split}$

we obtain,

$\displaystyle\frac{|g_{1}(s+h,x)-g_{1}(s,x)-D_{s}g_{1}(s,x).h|}{\|h\|}\leq$
$\displaystyle\qquad\qquad
k(\sqrt{d}+k+2)\left[\int_{0}^{1}\|s+ch-x\|^{-k-2}dc+\|s-x\|^{-k-2}\right].$
(37)

Integrating both sides of (37) over $\mathbb{R}^{d}\backslash B_{R}(s)$
yields,

$\begin{split}&\int_{\mathbb{R}^{d}\backslash
B_{R}(s)}\frac{|g_{1}(s+h,x)-g_{1}(s,x)-D_{s}g_{1}(s,x).h|}{\|h\|}\rho(x)dV(x)\\\
&\leq k(\sqrt{d}+k+2)\int_{\mathbb{R}^{d}\backslash
B_{R}(s)}\left[\int_{0}^{1}\|s+ch-x\|^{-k-2}dc+\|s-x\|^{-k-2}\right]\rho(x)dV(x)\\\
&\leq
k(\sqrt{d}+k+2)\int_{\mathbb{R}^{d}}\left[\int_{0}^{1}\|s+ch-x\|^{-k-2}dc+\|s-x\|^{-k-2}\right]\rho(x)dV(x).\end{split}$
(38)

By Lemma 3.3, we have that for all $c\in\left[0,1\right]$,

$\int_{\mathbb{R}^{d}}k(\sqrt{d}+k+2)\|s+ch-x\|^{-k-2}\rho(x)dV(x)<\infty.$
(39)

Denoting

$\nu=\underset{c\in[0,1]}{\max}\int_{\mathbb{R}^{d}}k(\sqrt{d}+k+2)\|s+ch-x\|^{-k-2}\rho(x)dV(x),$
(40)

we have immediately that,

$\int_{0}^{1}\int_{\mathbb{R}^{d}}k(\sqrt{d}+k+2)\left(\|s+ch-x\|^{-k-2}+\|s-x\|^{-k-2}\right)\rho(x)dV(x)dc\leq
2\nu.$ (41)

Now, from (38) we have from Fubini’s Theorem that,

$\int_{\mathbb{R}^{d}\backslash
B_{R}(s)}\,\frac{|g_{1}(s+h,x)-g_{1}(s,x)-D_{s}g_{1}(s,x)\cdot
h|}{\|h\|}\,\rho(x)\,dV(x)\leq 2\nu.$ (42)

By (36), for each $s$, the function
$(g_{1}(s+h,x)-g_{1}(s,x)-D_{s}g_{1}(s,x)\cdot h)/\|h\|$ converges pointwise
to $0$ as $h\rightarrow 0$, since

$\displaystyle\lim_{\|h\|\rightarrow
0}\frac{|g_{1}(s+h,x)-g_{1}(s,x)-D_{s}g_{1}(s,x)\cdot h|}{\|h\|}$
$\displaystyle\leq\lim_{\|h\|\rightarrow
0}\Bigg{\|}\int_{0}^{1}\\!D_{s}g_{1}(s+c\,h,x)\,dc-D_{s}g_{1}(s,x)\Bigg{\|}$
$\displaystyle=\Bigg{\|}\int_{0}^{1}\\!D_{s}g_{1}(s,x)\,dc-
D_{s}g_{1}(s,x)\Bigg{\|}$ $\displaystyle=0,$

whence, by Lebesgue’s Dominated Convergence Theorem,

$\underset{h\rightarrow 0}{\lim}\int_{\mathbb{R}^{d}\backslash
B_{R}(s)}\frac{|g_{1}(s+h,x)-g_{1}(s,x)-D_{s}g_{1}(s,x)\cdot
h|}{\|h\|}\rho(x)dV(x)=0.$ (43)

By the definition of the Frechet derivative,

$\displaystyle D_{s}\int_{\mathbb{R}^{d}\backslash
B_{R}(s)}g_{1}(s,x)\,\rho(x)\,dV(x)=\int_{\mathbb{R}^{d}\backslash
B_{R}(s)}\\!D_{s}g_{1}(s,x)\,\rho(x)\,dV(x)$
$\displaystyle=-k\int_{\mathbb{R}^{d}\backslash
B_{R}(s)}\\!\\!\\!\left(\frac{I_{d}}{\|s-x\|^{k+2}}-(k+2)\frac{(s-x)(s-x)^{t}}{\|s-x\|^{k+4}}\right)\rho(x)dV(x).$

Finally, letting $R\rightarrow 0$ yields the desired result by the same
argument used at the end of Lemma 3.3. ∎

##  Appendix F Proof of Lemma 3.5

###### Proof.

For $0<\varepsilon<1$, there exists $[t_{1},t_{2}]$ such that
$\mu([t_{1},t_{2}])>1-\varepsilon.$ Let $s\in[t_{1},t_{2}]$ and
$\\{s_{k}\\}_{k=1}^{\infty}$ be a sequence converging to $0$. Define,

$\displaystyle f_{\mu}(s)$
$\displaystyle=\int_{\mathbb{R}}(|s-x|-|x|)\rho(x)dx,$ (44)
$\displaystyle\phi(s;t_{1},t_{2})$
$\displaystyle=\int_{[t_{1},t_{2}]}(|s-x|-|x|)\rho(x)dx.$ (45)

As $\rho$ is a pdf, for each $s$ the function
$\theta(x):=(|s-x|-|x|)\,\rho(x)$ $\in L^{1}([t_{1},t_{2}])$. Therefore,

$\frac{\phi(s+s_{k};t_{1},t_{2})-\phi(s;t_{1},t_{2})}{s_{k}}=\frac{1}{s_{k}}\,\int_{[t_{1},t_{2}]}(|s+s_{k}-x|-|s-x|)\rho(x)dx:=A+B,$

where

$A=\frac{1}{s_{k}}\,\left[\int_{t_{1}}^{s+s_{k}}\\!(s+s_{k}-x)\,\rho(x)\,dx-\int_{t_{1}}^{s}\\!(s-x)\,\rho(x)\,dx\right],$

and

$B=\frac{1}{s_{k}}\,\left[\int_{s+s_{k}}^{t_{2}}\\!(s+s_{k}-x)\,\rho(x)\,dx-\int_{s}^{t_{2}}\\!(s-x)\,\rho(x)\,dx\right].$

Now note that

$\displaystyle A$
$\displaystyle=\frac{1}{s_{k}}\,\int_{t_{1}}^{s}\,s_{k}\rho(x)dx+\frac{1}{s_{k}}\int_{s}^{s+s_{k}}\\!\\!\\!\\!\\!(s+s_{k}-x)\,\rho(x)dx$
$\displaystyle=\mu([t_{1},s])+\frac{1}{s_{k}}\int_{s}^{s+s_{k}}\\!\\!\\!\\!\\!(s+s_{k}-x)\,\rho(x)\,dx,$

and similarly,

$B=-\mu([s,t_{2}])+\frac{1}{s_{k}}\int_{s}^{s+s_{k}}\\!\\!\\!\\!\\!(s+s_{k}-x)\,\rho(x)\,dx.$

Putting things together,

$\frac{\phi(s+s_{k};t_{1},t_{2})-\phi(s;t_{1},t_{2})}{s_{k}}=\mu([t_{1},s])-\mu([s,t_{2}])+\frac{2}{s_{k}}\int_{s}^{s+s_{k}}\\!\\!\\!\\!\\!(s+s_{k}-x)\,\rho(x)\,dx.$

By Theorem 2.1 applied to $(s,y)\in[t_{1},t_{2}]\times[t_{1},t_{2}]$, yields
that, for almost all $s\in[t_{1},t_{2}]$,

$\frac{d}{dy}\Big{|}_{y=s}\,\int_{s}^{y}\\!(y-x)\,\rho(x)\,dx=0\;\;\;\Longrightarrow\;\;\lim_{s_{k}\rightarrow
0}\,\frac{1}{s_{k}}\int_{s}^{s+s_{k}}\\!(s+s_{k}-x)\,\rho(x)\,dx=0.$ (46)

It follows that for almost all $s\in[t_{1},t_{2}]$,

$\phi^{\prime}(s)=\mu([t_{1},s])-\mu([s,t_{2}])=\mu([-\infty,s])-\mu([s,\infty])+2\,\varepsilon.$

Letting $\varepsilon\rightarrow 0$, we see that for almost all
$s\in\mathbb{R}$, $f^{\prime}_{\mu}(s)=\mu([-\infty,s])-\mu([s,\infty])$. Now,
applying Leibnitz’s rule (2nd fundamental theorem of calculus), yields that,
for almost all $s\in\mathbb{R}$,

$f^{\prime\prime}_{\mu}(s)=\rho(s)-(-\rho(s))=2\,\rho(s),$ (47)

thus proving the desired result. If $\rho$ is a continuous function, then (the
same Riemann integral version of) Leibnitz’s rule may be employed in (46)
instead of Theorem 2.1 to conclude (47) for all $s\in\mathbb{R}$. ∎

##  Appendix G Proof of Theorem 3.6

By Lemma 3.3, we can obtain the gradient of $f_{\mu}(s)$ as,

$\nabla_{s}f_{\mu}(s)=\mathbb{E}\left(\frac{s-X}{\|s-X\|_{2}}\right),$ (48)

and by Lemma 3.4 we can compute its Hessian as

$D_{s}\,\nabla_{s}f_{\mu}(s)=\mathbb{E}\left(\frac{I_{d}}{\|s-X\|_{2}}-\frac{(s-X)(s-X)^{t}}{\|s-X\|_{2}^{3}}\right),$
(49)

which we note is symmetric and positive definite. Hence, we can compute the
Laplacian of $f_{\mu}$ as

$\Delta_{s}f_{\mu}(s)=\text{trace}\left(D_{s}\,\nabla_{s}f_{\mu}(s)\right)=(d-1)\mathbb{E}\left(\frac{1}{\|s-X\|_{2}}\right).$
(50)

Now, by Lemma 3.2, $\Delta_{s}f_{\mu}(s)$ exists for $d\geq 3$, so that (16)
holds for $j=1$. Proceeding by induction, we suppose that (16) holds for all
$j\,\in\\{1,\ldots,N-1\\}$, and aim to show that (16) holds also for $j+1$.

To this effect, set $k=2\,j-1$, and using Lemmas 3.2 and 3.4, we compute,
respectively, the Frechet gradient and Hessian of $\mathbb{E}\|s-X\|_{2}^{-k}$
as:

$\nabla_{s}\mathbb{E}\left(\frac{1}{\|s-X\|_{2}^{k}}\right)=-k\,\mathbb{E}\left(\frac{s-X}{\|s-X\|_{2}^{k+2}}\right),$
(51)

and

$D_{s}\nabla_{s}\mathbb{E}\left(\frac{1}{\|s-X\|_{2}^{k}}\right)=-k\,\mathbb{E}\left(\frac{I_{d}}{\|s-X\|_{2}^{k+2}}-(k+2)\frac{(s-X)(s-X)^{t}}{\|s-X\|_{2}^{k+4}}\right).$
(52)

Hence, the Laplacian of $\mathbb{E}\|s-X\|_{2}^{-k}$ is given by:

$\Delta_{s}\mathbb{E}\left(\frac{1}{\|s-X\|_{2}^{k}}\right)=-k\,(d-k-2)\,\mathbb{E}\left(\frac{1}{\|s-X\|_{2}^{k+2}}\right).$
(53)

Now, the expectation on the right hand side of (53) exists by Lemma 3.3, as
long as $k\,\leq\,d-2$, or equivalently, $j\leq N$. Converting back to the $j$
variable, as
$\Delta^{j}_{s}f_{\mu}(s)=\Delta_{s}\,(\Delta^{j-1}_{s}f_{\mu}(s))$, yields
the final result.

##  Appendix H Proof of Corollary 1

###### Proof.

Consider the Poisson equation,

$\Delta\phi(s)=\rho(s),\quad\text{a.e.},$ (54)

whose fundamental solution is

$\phi(s)=\frac{\Gamma\left(\frac{d}{2}+1\right)}{d\,(d-2)\pi^{d/2}}\mathbb{E}\left(\frac{1}{\|s-X\|_{2}^{d-2}}\right).$

Putting (18) and (54) together, we obtain:

$\displaystyle\Delta^{N}f_{\mu}(s)$
$\displaystyle=(-1)^{N-1}\,(d-1)\,\cdots\,(2)\mathbb{E}\left(\frac{1}{\|s-X\|_{2}^{d-2}}\right)$
$\displaystyle=(-1)^{N-1}\,(d-1)\,\cdots\,(2)\,d\,(d-2)\,\frac{\pi^{d/2}}{\Gamma\left(\frac{d}{2}+1\right)}\,\phi(s),$

so that

$\displaystyle\Delta\,\Delta^{N}f_{\mu}(s)$
$\displaystyle=(-1)^{N-1}\,(d-1)\,\cdots\,(2)\mathbb{E}\left(\frac{1}{\|s-X\|_{2}^{d-2}}\right)$
$\displaystyle=(-1)^{N-1}\,(d-1)\,\cdots\,(2)\,d\,(d-2)\,\frac{\pi^{d/2}}{\Gamma\left(\frac{d}{2}+1\right)}\,\rho(s),\quad\text{a.e},$

thus proving the result. ∎

## References

  * [1] [author] Aggarwal, CharuC. (2015). Data mining. Springer.
  * [2] [author] Aryal, GaurabG. (2016). Identifying a model of screening with multidimensional consumer heterogeneity. Available at SSRN 2531188.
  * [3] [author] Ben-Gal, IradI. (2005). Outlier detection. In Data mining and knowledge discovery handbook 131–146. Springer.
  * [4] [author] Breckling, JensJ. and Chambers, RayR. (1988). M-quantiles. Biometrika 75 761–771.
  * [5] [author] Brown, BMB. (1983). Statistical uses of the spatial median. Journal of the Royal Statistical Society: Series B (Methodological) 45 25–30.
  * [6] [author] Chakraborty, AnirvanA. and Chaudhuri, ProbalP. (2014). The spatial distribution in infinite dimensional spaces and related quantiles and depths. The Annals of Statistics 42 1203–1231.
  * [7] [author] Chakraborty, BimanB. (2001). On affine equivariant multivariate quantiles. Annals of the Institute of Statistical Mathematics 53 380–403.
  * [8] [author] Chakraborty, BimanB. and Chaudhuri, ProbalP. (1996). On a transformation and re-transformation technique for constructing an affine equivariant multivariate median. Proceedings of the American mathematical society 124 2539–2547.
  * [9] [author] Chandola, VarunV., Banerjee, ArindamA. and Kumar, VipinV. (2009). Anomaly detection: a survey. ACM computing surveys 41 15\.
  * [10] [author] Chaouch, MohamedM. and Goga, CameliaC. (2010). Design-based estimation for geometric quantiles with application to outlier detection. Computational Statistics & Data Analysis 54 2214–2229.
  * [11] [author] Chaudhuri, ProbalP. (1992). Multivariate location estimation using extension of R-estimates through U-statistics type approach. The Annals of Statistics 20 897–916.
  * [12] [author] Chaudhuri, ProbalP. (1996). On a geometric notion of quantiles for multivariate data. Journal of the American statistical association 91 862–872.
  * [13] [author] Chen, YixinY., Dang, XinX., Peng, HanxiangH. and Bart Jr. , HenryH. (2008). Outlier detection with the kernelized spatial depth function. IEEE Transactions on Pattern Analysis and Machine Intelligence 31 288-305.
  * [14] [author] Chernozhukov, VictorV., Galichon, AlfredA., Hallin, MarcM. and Henry, MarcM. (2017). Monge–Kantorovich depth, quantiles, ranks and signs. The Annals of Statistics 45 223–256.
  * [15] [author] Church, James CJ. C., Chen, YixinY. and Rice, Stephen VS. V. (2008). A spatial median filter for noise removal in digital images. In IEEE SoutheastCon 2008 618–623. IEEE.
  * [16] [author] Dhar, Subhra SankarS. S., Chakraborty, BimanB. and Chaudhuri, ProbalP. (2014). Comparison of multivariate distributions using quantile–quantile plots and related tests. Bernoulli 20 1484 – 1506.
  * [17] [author] Field, ChristopherC. and Genton, Marc GM. G. (2006). The multivariate g-and-h distribution. Technometrics 48 104–111.
  * [18] [author] Fraiman, RicardoR. and Pateiro-López, BeatrizB. (2012). Quantiles for finite and infinite dimensional data. Journal of Multivariate Analysis 108 1–14.
  * [19] [author] Ghosal, PromitP. and Sen, BodhisattvaB. (2022). Multivariate ranks and quantiles using optimal transport: Consistency, rates and nonparametric testing. The Annals of Statistics 50 1012–1037.
  * [20] [author] Girard, StéphaneS. and Stupfler, GillesG. (2017). Intriguing properties of extreme geometric quantiles. REVSTAT-Statistical Journal 15 107–139.
  * [21] [author] Ha, JihyunJ., Seok, SeulgiS. and Lee, Jong-SeokJ.-S. (2015). A precise ranking method for outlier detection. Information Sciences 324 88–107.
  * [22] [author] Haldane, JBSJ. (1948). Note on the median of a multivariate distribution. Biometrika 35 414–417.
  * [23] [author] Hallin, MarcM., Del Barrio, EustasioE., Cuesta-Albertos, JuanJ. and Matrán, CarlosC. (2021). Distribution and quantile functions, ranks and signs in dimension d: A measure transportation approach. The Annals of Statistics 49 1139–1165.
  * [24] [author] Herrmann, KlausK., Hofert, MariusM. and Mailhot, MélinaM. (2018). Multivariate geometric expectiles. Scandinavian Actuarial Journal 2018 629–659.
  * [25] [author] Koltchinskii, Vladimir IV. I. (1997). M-estimation, convexity and quantiles. The Annals of Statistics 25 435–477.
  * [26] [author] Kong, LinglongL. and Mizera, IvanI. (2012). Quantile tomography: using quantiles with multivariate data. Statistica Sinica 1589–1610.
  * [27] [author] Lange, KennethK. (2016). MM optimization algorithms. SIAM, Philadelphia, PA.
  * [28] [author] Li, JunJ. (2018). EDF goodness-of-fit tests based on centre-outward ordering. Journal of Nonparametric Statistics 30 973–989.
  * [29] [author] Li, JunJ., Zhang, XinX. and Jeske, Daniel RD. R. (2013). Nonparametric multivariate CUSUM control charts for location and scale changes. Journal of Nonparametric Statistics 25 1–20.
  * [30] [author] Lopuhaa, H.H. and Rouseeuw, P.P. (1991). Breakdown points of affine equivariant estimators of multivariate location and covariance matrices. The Annals of Statistics 19 229–248.
  * [31] [author] Marden, John IJ. I. (2004). Positions and QQ plots. Statistical Science 606–614.
  * [32] [author] McCann, R. J.R. J. (2001). Existence and uniqueness of monotone measure-preserving maps. Duke Mathematical Journal 80 309–323.
  * [33] [author] McShane, E. J.E. J. and Botts, T. A.T. A. (2005). Real Analysis. Dover reprinted from 1959.
  * [34] [author] Mosler, KarlK. (2013). Depth Statistics In Robustness and Complex Data Structures: Festschrift in Honour of Ursula Gather 17–34. Springer, Berlin, Heidelberg.
  * [35] [author] Mukhopadhyay, Nitai DN. D. and Chatterjee, SnigdhansuS. (2011). High dimensional data analysis using multivariate generalized spatial quantiles. Journal of multivariate analysis 102 768–780.
  * [36] [author] Oja, HannuH. (1983). Descriptive statistics for multivariate distributions. Statistics & Probability Letters 1 327–332.
  * [37] [author] Oja, HannuH. (2013). Multivariate Median In Robustness and Complex Data Structures: Festschrift in Honour of Ursula Gather 3–15. Springer, Berlin, Heidelberg.
  * [38] [author] Puri, M. L.M. L. and Sen, P. K.P. K. (1971). Nonparametric methods in multivariate analysis. John Wiley & sons.
  * [39] [author] Rizzi, A.A., Vichi, M.M. and Bock, H. H.H. H. (1998). Multivariate $L^{1}$ median. In Proc. 6th conf. of the International Federation of Classification Societies 539–546.
  * [40] [author] Rockafellar, R TyrrellR. T. (1970). Convex analysis 18\. Princeton university press.
  * [41] [author] Rudin, WalterW. (1976). Principles of mathematical analysis, 3 ed. McGraw Hill International.
  * [42] [author] Serfling, RobertR. (2002). Quantile functions for multivariate analysis: approaches and applications. Statistica Neerlandica 56 214–232.
  * [43] [author] Serfling, RobertR. (2002). A depth function and a scale curve based on spatial quantiles. In Statistical data analysis based on the L1-norm and related methods 25–38. Springer.
  * [44] [author] Small, C. G.C. G. (1990). A survey of multidimensional medians. International Statistical Review 58 263-277.
  * [45] [author] Stroock, DanielD. (2011). Probability Theory: an analytic view. Cambridge University Press.
  * [46] [author] Tarabelloni, NicholasN., Schenone, ElisaE., Collin, AnnabelleA., Ieva, FrancescaF., Paganoni, Anna MariaA. M., Gerbeau, Jean-FredericJ.-F. et al. (2018). Statistical assessment and calibration of numerical ecg models. JP Journal of Biostatistics 15 151–173.
  * [47] [author] Tukey, John WJ. W. (1977). Exploratory data analysis. Addison-Wesley Publishing Company.
  * [48] [author] Xie, JiangJ., Xiong, ZhongyangZ., Dai, QizhuQ., Wang, XiaoxiaX. and Zhang, YufangY. (2020). A local-gravitation-based method for the detection of outliers and boundary points. Knowledge-Based Systems 192 105331\.
  * [49] [author] Zuo, YijunY. and Serfling, RobertR. (2000). General notions of statistical depth function. The Annals of statistics 461–482.