Marc Hallin, Université libre de Bruxelles, ECARES and Département de
Mathématique, Brussels, B-1050, Belgium<EMAIL_ADDRESS>Thomas Verdebout’s
research is supported by an ARC grant of the ULB and the Communauté Française
de Belgique, a Projet de Recherche (PDR) grant of the FNRS, and the Fonds
Thelam of the Fondation Roi Baudouin.

# Nonparametric Measure-Transportation-Based Methods for Directional Data

Hallin M  Liu H  International Institute of Finance, School of Management,
University of Science and Technology of China  Verdebout T

###### Abstract

This paper proposes various nonparametric tools based on measure
transportation for directional data. We use optimal transports to define new
notions of distribution and quantile functions on the hypersphere, with
meaningful quantile contours and regions yielding closed-form formulas under
the classical assumption of rotational symmetry. The empirical versions of our
distribution functions enjoy the expected Glivenko-Cantelli property of
traditional distribution functions. They provide fully distribution-free
concepts of ranks and signs and define data-driven systems of (curvilinear)
parallels and (hyper)meridians. Based on this, we also construct a universally
consistent test of uniformity and a class of fully distribution-free tests for
directional MANOVA which, in simulations, outperform all their existing
competitors. Two real-data examples involving the analysis of sunspots and the
structure of proteins conclude the paper.

###### keywords:

directional statistics, directional quantiles, ranks, signs, optimal
transport, directional MANOVA

##  1 Introduction

Directional data analysis is dealing with random directions—taking values on
circles and (hyper)spheres or, more generally, on manifolds. Directional data
can be found in a variety of fields, including astronomy [66, 65],
environmetrics [72, 4, 54], biology and medicine [29, 45, 28], to cite only a
few. Directions generally are represented as points on the unit hypersphere
$\mathcal{S}^{d-1}\coloneqq\\{\textbf{x}\in{\mathbb{R}}^{d}:\|\textbf{x}\|^{2}\coloneqq\mathbf{x}^{\top}\mathbf{x}=1\\}$
in ${\mathbb{R}}^{d}$. The special nature of their sample spaces gives
directional statistics a specific flavor. The models of interest are of the
standard statistical types (one-sample and multi-sample location, analysis of
variance, regression, etc.) and inference is carried out according to the
usual principles, but their implementation takes special forms and requires
distinctive techniques. A comprehensive exposition of the theoretical
background and inference methods for circular and spherical data can be found
in [64] or [74]. More recent advances are presented in [60, 61], while
statistical tools for data on more general spaces are studied in [22] and [15,
16].

The most popular parametric model in directional statistics, which can be
traced back to the first decades of the 20th century, is the von Mises-Fisher
model, characterized by the family of von Mises distributions, with densities
(with respect to the classical surface area measure on $\mathcal{S}^{d-1}$) of
the form

$\mathbf{z}\in\mathcal{S}^{d-1}\mapsto c_{\kappa}{\rm
exp}(\kappa\mathbf{z}^{\top}{\boldsymbol{\theta}}),\vspace{-1mm}$ (1.1)

where ${\boldsymbol{\theta}}\in\mathcal{S}^{d-1}$ plays the role of a location
parameter, $\kappa\in\mathbb{R}^{+}$, which drives the probability mass in the
vicinity of ${\boldsymbol{\theta}}$, is the so-called concentration parameter,
$c_{\kappa}$ is a normalizing constant, and $\mathbf{z}^{\top}$ denotes the
transpose of $\mathbf{z}$. The von Mises distribution is often seen as the
“directional Gaussian” distribution due to the fact that the maximum
likelihood estimator of ${\boldsymbol{\theta}}$ is the normalized sample
average. The von Mises distributions belong to the more general class of
rotationally symmetric distributions, which contains all distributions with
densities (still with respect to the surface area measure on
$\mathcal{S}^{d-1}$) of the form

$\mathbf{z}\in\mathcal{S}^{d-1}\mapsto
c_{f}f(\mathbf{z}^{\top}{\boldsymbol{\theta}}),\vspace{-1mm}$ (1.2)

where ${\boldsymbol{\theta}}\in\mathcal{S}^{d-1}$, $f$ is some positive
angular function, and $c_{f}$ is a norming constant. The projection $({\bf
I}_{p}-{\boldsymbol{\theta}}{\boldsymbol{\theta}}^{\top})\mathbf{Z}$ of a
rotationally symmetric random vector $\mathbf{Z}$ has a spherical distribution
in the tangent space (to $\mathcal{S}^{d-1}$) at ${\boldsymbol{\theta}}$.
Inference for the parameters of rotationally symmetric distributions has
recently been considered in [23], [47], and [70, 71].

Extensions of rotational symmetry yielding, after projection in the tangent
space (to $\mathcal{S}^{d-1}$) at ${\boldsymbol{\theta}}$, distributions with
elliptical contours have been proposed in [48], [78], and [33]; see also [53],
[54], and [49]. Rank-based methods (including MANOVA) for directional data
have been proposed in [58, 59] and [87]. While the concepts of ranks
considered there enjoy several attractive features, their distribution-
freeness, hence also their applicability is limited to the class of
rotationally symmetryic distributions Based on sample projections onto the
mean direction, the quantiles studied in [57] are canonical for rotationally
symmetric distributions on $\mathcal{S}^{d-1}$ but fail to be informative for
general distributions on $\mathcal{S}^{d-1}$. The assumptions of rotational
symmetry and their “elliptical extensions” are pervasive in the literature on
directional data; they can be compared, in this respect, to the assumptions of
spherical and elliptical symmetry in traditional multivariate analysis. Just
as the latter, they drastically simplify the mathematical structure of
inference problems but also very severely restrict their validity: in most
applications, rotational symmetry, indeed, is extremely unlikely to hold.

The situation in the directional context, thus, is quite comparable to the
situation in the general multivariate context: the absence of a canonical
ordering, at first sight, precludes the canonical definition of essential
order-related statistical tools as ranks, signs, and quantiles unless one is
ready to make a very strong assumption of rotational symmetry (or its
“elliptical extensions") in the directional case, of spherical or elliptical
symmetry in the general multivariate case. Based on measure transportation
ideas and novel concepts of multivariate distribution and quantile functions
in ${\mathbb{R}}^{d}$, new non- and semi-parametric methods have been proposed
very recently in [21] and [40], which have demonstrated their efficiency in a
variety of multivariate inference problems: see [26], [44], [35], [41, 42,
43], [82, 83], [8], [27]. [85] use optimal transport to perform causal
inference, while machine learning methods based on optimal transport have been
proposed, e.g., in [51]. See [39] for a nontechnical survey.

The objective of this paper is to develop a measure-transportation-based
approach for directional data similar to the one adopted in the above
references, all of which are dealing with observations in ${\mathbb{R}}^{d}$.
Measure transportation on hyperspheres and more general manifolds, however, is
not a cosmetic adaptation of the Euclidean case, as it involves non-Euclidean
topologies and metrics. As a consequence, the measure-transportation-based
concepts and methods developed here do not follow along the same lines as in
[21] and [40]. Below, we mainly build upon [68], where the celebrated results
of [67] on measure-preserving monotone mappings in ${\mathbb{R}}^{d}$ are
extended to Riemanniann manifolds, to propose directional concepts of
distribution and quantile functions inducing a distribution-specific system of
curvilinear parallels and hypermeridians. These directional distribution
functions can be seen as canonical transformations in the sense of [46]. The
empirical counterparts of these concepts yield directional ranks, signs, and
empirical quantiles enjoying the properties one is expecting from such notions
(in particular, distribution-freeness of ranks and signs and quantile regions
with prescribed probability contents irrespective of the underlying
distribution); empirical distribution functions, moreover, are shown to
satisfy a Glivenko-Cantelli consistency property.

The motivation for these definitions are the same as in ${\mathbb{R}}$ and
${\mathbb{R}}^{d}$: conducting distribution-free inference based on ranks and
signs, constructing quantile contours and quantile regions (closed, connected,
and nested, with probability contents that do not depend on the underlying
distribution), hence directional outliers and values at risk, performing
directional quantile regression, etc.

Earlier attemps have been made with similar objectives, which are only
partially successful: the notion of ranks studied in [58, 59] and [87],
indeed, is not enjoying distribution-freeness unless the underlying
distribution is rotationally symmetric; similarly, the quantile contours
proposed in [57] are satisfactory under rotational symmetry only.

In order to show that, contrary to these earlier proposals, our measure-
transportation-based concepts are achieving these objectives, we consider two
very classical problems in directional inference:

  1. (i)

(testing for directional uniformity) the classical problem of testing
uniformity on $\mathcal{S}^{d-1}$. This GOF problem is probably the oldest and
most fundamental one in directional statistics, which can be traced back to
the discussion by [12] on whether the closeness of the orbital planes of
various planets arose by chance or not. To cite only a few recent works, tests
of uniformity in the context of noisy directional data were studied in [55]
and [50] and, in a high-dimensional context, in [20], [19], and [25].
Projection-based tests were proposed by [24] and [32]. Tests of uniformity
were used in [33] to construct tests for rotational symmetry and inference for
location in the vicinity of uniformity was considered in [69]. We refer to
[34] for a recent review of this topic.

  2. (ii)

(directional MANOVA) MANOVA on the hypersphere—testing the hypothesis of no-
treatment effect, where the treatment can impact location, concentration, or
skewness, is another very classical problem in directional inference. Due to
the absence of distribution-free tools, the traditional approach typically
reduces to a pseudo-von Mises MANOVA procedure—see for instance the recent
contributions [30], [59], [81] and [52] and the references therein; a major
weakness of pseudo-von Mises methods, however, is that their validity, which
requires the classical Fisher consistency property, is guaranteed under
rotational symmetry only.

In both cases, our directional measure-transportation-based concepts provide
natural nonparametric and fully distribution-free solutions which (in
simulations) perform equally well as some of their competitors under
rotational symmetry, and outperform them all in the non-rotationally symmetric
case:

  1. (i)

(testing for directional uniformity) the directional distribution function for
the uniform over $\mathcal{S}^{d-1}$ is the identity function; a (fully
distribution-free) Cramér-von Mises-type test of uniformity on
$\mathcal{S}^{d-1}$ based on the empirical directional distribution function
is shown to be universally consistent while outperforming (in simulations) its
traditional competitors (the Rayleigh, Bingham, Ajne, Giné, and Bakshaev
tests), as well as the more recent projection-based tests by [32]. In sharp
contrast with the universal consistency of our measure-transportation-based
Cramér-von Mises test, these competitors, moreover, typically have
undetectable blind spots. Also note that our measure-transportation-based
Cramér-von Mises tests straightforwardly extend from testing uniformity to the
general GOF problem;

  2. (ii)

(directional MANOVA) the class of fully distribution-free MANOVA tests we are
proposing is based on directional ranks and signs; as in the traditional case
of univariate rank tests, it involves the choice of a score function, which
can target various types of alternatives. We provide the asymptotic
(noncentral chi-square) distributions of the test statistics under the null
and under local alternatives and show via simulations how our tests outperform
their classical pseudo-von Mises competitors (the validity of which, moreover,
is limited to the rotationally symmetric case).

The paper is organized as follows. In Section 2, we provide the relevant
mathematical notions and briefly discuss the concept of measure transportation
on Riemannian manifolds and hyperspheres. We define the population versions
$\bf F$ and $\bf Q$ of our directional distribution and quantile functions in
Section 3, then derive their explicit forms in the rotationally symmetric
case. In Section 4, we introduce the empirical version ${\bf F}^{(n)}$ of $\bf
F$ and establish its Glivenko-Cantelli asymptotic behavior. We also introduce
the resulting ranks and signs and discuss (see also Supplementary Material A)
their distributional and equivariance properties. Section 4.3 illustrates the
construction of empirical quantile contours in simulated data. In Section 5,
based on ${\bf F}^{(n)}$, we construct Cramér-von Mises-type tests for
goodness-of-fit problems, establish their universal consistence against fixed
alternatives, and compare their performance with that of its many competitors.
In Section 6, we propose a class of distribution-free MANOVA testing
procedures and study their asymptotic distributions under the null and under
local alternatives. Real-data applications are presented in Section 7 where we
implement our MANOVA procedure on sunspots data and illustrate our novel
concept of directional quantiles in an analysis of proteins structures.
Finally, Supplementary Material A and B are collecting the basic properties of
our ranks and signs and the proofs of the various results provided in the
paper, respectively.

##  2 Measure transportation on Riemannian manifolds

Throughout the paper, we denote by
$d(\mathbf{y},\mathbf{z})\coloneqq\left|\arccos(\mathbf{y}^{\top}\mathbf{z})\right|$
the geodesic distance and by $c(\mathbf{y},\mathbf{z})\coloneqq
d^{2}(\mathbf{y},\mathbf{z})/2$ the squared Riemannian distance between two
points $\mathbf{y}$ and $\mathbf{z}$ on the unit hypersphere
$\mathcal{S}^{d-1}$; both $c$ and $d$ are continuous and bounded. When
equipped with the geodesic distance, $\mathcal{S}^{d-1}$ is a separable
complete metric space, hence a Polish metric space, with Borel $\sigma$-field
${\cal B}^{d-1}$. Denote by $\lambda^{d-1}$ the surface area measure over
$\mathcal{S}^{d-1}$. Throughout, measurability tacitly is understood with
respect to ${\cal B}^{d-1}$ and ${\rm P}$ is assumed to be in the family
$\mathfrak{P}_{d}$ of $\lambda^{d-1}$-absolutely continuous distributions with
densities bounded away from 0 and $\infty$.

Measure transportation on Polish spaces and, more particularly, on Riemannian
manifolds, is a well-studied subject; below, we mainly follow [68], where the
celebrated results of [67] on measure-preserving monotone mappings in
${\mathbb{R}}^{d}$ are extended to Riemanniann manifolds, [3], and [79]. Let
${\rm P}$ and ${\rm Q}$ denote two probability measures on
$\mathcal{S}^{d-1}$. Monge’s optimal transport problem on ${\cal S}^{d-1}$
consists in minimizing, among the set ${\mathcal{S}}({\rm P},{\rm Q})$ of all
measurable transport maps ${\bf
S}\\!:\mathcal{S}^{d-1}\\!\to\\!\mathcal{S}^{d-1}$ suchthat
$(\mathbf{S}\\#{\rm P})(V)\coloneqq{\rm P}(\mathbf{S}^{-1}(V))={\rm Q}({V})$
for all $V\\!\in\\!{\cal B}^{d-1}$ (in the measure transportation notation and
termino-logy, $\mathbf{S}\\!\in\\!{\mathcal{S}}({\rm P},{\rm Q})$ is pushing
${\rm P}$ forward to ${\rm Q}$), the transportation cost

$C_{\rm M}({\bf S})\coloneqq\int_{{\cal S}^{d-1}}c\big{(}\mathbf{z},{\bf
S}(\mathbf{z})\big{)}{\mathrm{d}}{\rm P}(\mathbf{z})={\rm E}_{\rm
P}\big{[}c(\mathbf{Z},{\bf S}(\mathbf{Z}))\big{]}.$ (2.1)

Due to the fact that ${\mathcal{S}}({\rm P},{\rm Q})$ is not convex, Monge’s
problem is an uneasy one.

Closely related to Monge’s problem is the so-called Kantorovich problem.
Denote by $\Gamma({\rm P},{\rm Q})$ the set of all probability measures
$\gamma$ on $\mathcal{S}^{d-1}\\!\times~{}\\!\mathcal{S}^{d-1}\\!$ with
“marginals" ${\rm P}$ and ${\rm Q}$; an element $\gamma$ of $\Gamma({\rm
P},{\rm Q})$ is called a transport plan. The Kantorovich problem on ${\cal
S}^{d-1}$ consists in minimizing, among all transport plans
$\gamma\in\Gamma({\rm P},{\rm Q})$, the expectation

$C_{\rm
K}(\gamma)\coloneqq\int_{\mathcal{S}^{d-1}\times\mathcal{S}^{d-1}}c(\mathbf{y},\mathbf{z}){\mathrm{d}}\gamma(\mathbf{y},\mathbf{z})\vspace{-1mm}$
(2.2)

of the transportation cost $c(\mathbf{Y},\mathbf{Z})$ under
$(\mathbf{Y},\mathbf{Z})\sim\gamma$. Unlike ${\mathcal{S}}({\rm P},{\rm Q})$,
the set $\Gamma({\rm P},{\rm Q})$ is a convex subset of a Banach space; since
$c$ is bounded, the existence of a solution is guaranteed.

The Kantorovich problem is a relaxation of Monge’s. Letting $({\bf I}_{\rm
d}~{}\\!\times~{}\\!{\bf S}){\bf Z}\coloneqq~{}\\!({\bf Z},{\bf S}({\bf Z}))$,
a transport map ${\bf S}\in{\mathcal{S}}({\rm P},{\rm Q})$ indeed induces a
transport plan $\gamma=({\bf I}_{\rm d}\times{\bf S})\\#{\rm P}\in\Gamma({\rm
P},{\rm Q})$ which is concentrated on the graph of $\mathbf{S}$ (i.e., such
that $\gamma\big{(}({\bf Z},{\bf S}({\bf
Z}))\in\big{\\{}(\mathbf{z},\mathbf{S}(\mathbf{z})):\mathbf{z}\in\mathcal{S}^{d-1}\\}\big{)}~{}\\!\\!=~{}\\!\\!1$)
and satisfying $C_{\rm K}(\gamma)=C_{\rm M}(\mathbf{S})$. As we shall see
(Proposition 2.3 below), the solution of Kantorovich’s problem is precisely of
the form $({\bf I}_{\rm d}\times{\bf F})\\#{\rm P}$ where ${\bf F}$, thus, is
a solution of Monge’s problem, and

$C_{\rm M}(\mathbf{F})=\min_{\mathbf{S}\in{\cal S}({\rm P},{\rm Q})}C_{\rm
M}(\mathbf{S})=\min_{\gamma\in{\Gamma}({\rm P},{\rm Q})}C_{\rm
K}(\gamma)=C_{\rm K}(({\bf I}_{\rm d}\times{\bf F})\\#{\rm P}),\vspace{-1mm}$

where $\min_{\gamma\in{\Gamma}({\rm P},{\rm Q})}C^{1/2}_{\rm
K}(\gamma)\eqqcolon{\cal W}_{2}({\rm P},{\rm Q})$ is the Wasserstein distance
between ${\rm P}$ and ${\rm Q}$.

Before providing more precise statements about the solutions of the Monge and
Kantorovich problems, let us introduce some notation from differential
geometry. Recall that the tangent space (tangent to $\mathcal{S}^{d-1}$) at
$\mathbf{x}\in\mathcal{S}^{d-1}$ is the $(d-1)$-dimensional linear subspace
${\cal
T}_{\mathbf{x}}\mathcal{S}^{d-1}\coloneqq\\{\mathbf{z}\in\mathbb{R}^{d},\mathbf{z}^{\top}\mathbf{x}=0\\}.$
Letting $\mathbf{x}\in\mathcal{S}^{d-1}$ and ${\bf v}\in{\cal
T}_{\mathbf{x}}\mathcal{S}^{d-1}$, the directional derivative at
$\mathbf{x}\in\mathcal{S}^{d-1}$ of a smooth function
$\varsigma:\mathcal{S}^{d-1}\to\mathbb{R}$ in direction ${\bf v}$ is defined
as ${v}({\varsigma})\coloneqq(\varsigma\circ\nu)^{\prime}(0),$ where
$\nu~{}\\!:~{}\\![0,1]~{}\\!\to~{}\\!\mathcal{S}^{d-1}$ is a differentiable
path such that $\nu(0)=\mathbf{x}$ and $\nu^{\prime}(0)={\bf v}$. The gradient
$\nabla\varsigma(\mathbf{x})$ of $\varsigma$ at
$\mathbf{x}\in\mathcal{S}^{d-1}$ then is defined as the vector in ${\cal
T}_{\mathbf{x}}\mathcal{S}^{d-1}$ such that ${\bf
v}^{\top}\nabla\varsigma(\mathbf{x})={v}({\varsigma})$ for all ${\bf
v}\in{\cal T}_{\mathbf{x}}\mathcal{S}^{d-1}$ and
$\mathbf{x}\in\mathcal{S}^{d-1}$. The exponential map at a point
$\mathbf{x}\in~{}\\!\mathcal{S}^{d-1}$ is a map from the tangent space ${\cal
T}_{\mathbf{x}}\mathcal{S}^{d-1}$ to $\mathcal{S}^{d-1}$. More precisely,
denoting by $\nu_{\bf v}(t)$ the unique geodesic with tangent vector ${\bf
v}\in{\cal T}_{\mathbf{x}}\mathcal{S}^{d-1}$ running through
$\mathbf{x}\in\mathcal{S}^{d-1}$ ($\nu_{\bf v}(0)=\mathbf{x}$), then
$\exp_{\mathbf{x}}({\bf v})\coloneqq\nu_{\bf v}(1)$. Finally, consider the
following concepts of $c$-cyclical monotonicity and $c$-convexity.

######  Definition 2.1.

A subset $\Omega\in{\cal S}^{d-1}\times{\cal S}^{d-1}$ is called
$c$-cyclically monotone ($c$ the squared Riemannian distance on ${\cal
S}^{d-1}$) if, denoting by $\Sigma(k)$ the set of permutations of
$\\{1,\ldots,k\\}$, for all $k\in{\mathbb{N}}$, all $\sigma\in\Sigma(k)$, and
all
$(\mathbf{x}_{1},\mathbf{y}_{1}),\ldots,(\mathbf{x}_{k},\mathbf{y}_{k})\in\Omega$,

$\sum_{i=1}^{k}c(\mathbf{x}_{i},\mathbf{y}_{i})\leq\sum_{i=1}^{k}c(\mathbf{x}_{\sigma(i)},\mathbf{y}_{i}).$

######  Definition 2.2.

Two functions $\psi$ and $\phi$ from $\mathcal{S}^{d-1}$ to $\mathbb{R}$ such
that

$\psi(\mathbf{x})=\\!\\!\sup_{\mathbf{y}\in\mathcal{S}^{d-1}}\\{\phi(\mathbf{y})-c(\mathbf{x},\mathbf{y})\\}\text{
and }\
\phi(\mathbf{y})=\\!\\!\inf_{\mathbf{x}\in\mathcal{S}^{d-1}}\\{\psi(\mathbf{x})+c(\mathbf{x},\mathbf{y})\\}$
(2.3)

are called $c$-convex and $c$-concave, respectively; call $\phi$ (resp.
$\psi$) the c-transform of $\psi$ (resp. $\phi$).

Since the squared Riemannian distance $c$ satisfies the conditions of Theorem
10.26(i) of [88], the functions $\psi$ and $\phi$ in (2.3) are a.e.
differentiable. The following theorem summarizes, for the Monge and
Kantorovich problems with transportation cost the squared Riemannian distance
on ${\cal S}^{d-1}$, some of the results contained in Chapter 5 of [88] (for
the necessity part of (i), first established by [77]), Theorem 1 of [79] (for
the sufficiency part of (i); see also [73]), Theorems 8 and 9 and Corollary 10
of [68] (for (ii), (iii), and (iv)). The assumptions made in these references
are automatically satisfied here, since ${\cal S}^{d-1}$ equipped with the
geodesic distance is a Polish metric space and $c$ is a continuous bounded
cost function.

######  Proposition 2.3.

Let $\rm P\in\mathfrak{P}_{d}$ and $\rm Q$ denote two probability measures on
$\mathcal{S}^{d-1}$. Then,

  1. (i)

a transport plan $\gamma\in\Gamma({\rm P},{\rm Q})$ is the solution
(minimizing $C_{\rm K}$ in (2.2) over $\Gamma({\rm P},{\rm Q})$) of the
Kantorovich problem if and only if it is supported on a $c$-cyclically
monotone subset of ${\cal S}^{d-1}\times{\cal S}^{d-1}$;

  2. (ii)

this solution of the Kantorovich problem exists, is unique, and is of the form
$({\bf I}_{\rm d}\times{\bf F})\\#{\rm P}$ where ${\bf F}\in{\mathcal{S}}({\rm
P},{\rm Q})$ is the $\rm P$-a.s. unique solution of the corresponding Monge
problem (minimizing $C_{\rm M}$ in (2.1) over ${\mathcal{S}}({\rm P},{\rm
Q})$);

  3. (iii)

there exist $c$-concave differentiable mappings $\psi$ from
$\mathcal{S}^{d-1}$ to $\mathbb{R}$ such that ${\bf
F}(\mathbf{x})=\exp_{\mathbf{x}}(-\nabla\psi(\mathbf{x}))$,
$\lambda^{d-1}\text{-a.e.}$

If, moreover, $\rm Q\in\mathfrak{P}_{d}$, then

  1. (iv)

${\bf Q}(\mathbf{x})\coloneqq\exp_{\mathbf{x}}(-\nabla\phi(\mathbf{x}))$, with
$\phi$ the $c$-transform of $\psi$, belongs to ${\mathcal{S}}({\rm Q},{\rm
P})$, is the $\rm Q$-a.s. unique minimizer of $C_{\rm M}$ in (2.1) over
${\mathcal{S}}({\rm Q},{\rm P})$ (i.e., the $\rm Q$-a.s. unique solution of
the corresponding Monge problem), and satisfies

${\bf Q}({\bf F}(\mathbf{x}))=\mathbf{x}\quad\rm P\text{-a.s.\quad and
\quad}{\bf F}({\bf Q}(\mathbf{x}))=\mathbf{x}\quad\rm Q\text{-a.s.}$

Proposition 2.3 entails the existence of a mapping ${\bf F}$ pushing any
probability measure ${\rm P}\in\mathfrak{P}_{d}$ forward to any other
probability measure ${\rm Q}\in\mathfrak{P}_{d}$ and minimizing, among all
such mappings, the transportation cost (2.1). In the next Section, we will use
this optimality result to define directional distribution and quantile
functions.

Our results all require densities that are bounded away from zero and
$\infty$. It would be interesting to see if our results and our methods could
be extended (i) to densities supported on a subset of the hypersphere and (ii)
to other Riemanian manifolds. Encouraging results on the uniform consistency
of empirical transports under (i) have been obtained in [7], [6], and [80].
All these results are limited to the Euclidean case, though, and extending
them to the hypersphere is beyond the scope of this paper. As for (ii),
whether our methodology adapts to more general Riemannian manifolds, and how,
is a challenging question which, definitely, should be left for future
research.

##  3 Directional distribution and quantile functions

###  3.1 Definitions and some basic properties

The concepts of directional distribution and quantile functions we are
proposing are based on the optimal transport ${\bf F}$ from a generic
probability measure ${\rm P}\in\mathfrak{P}_{d}$ to the uniform probability
measure ${\rm P}^{\bf U}$ on $\mathcal{S}^{d-1}$. The idea is inspired by the
so-called center-outward distribution and quantile functions
${\mathbf{F}_{\pm}}$ and ${\mathbf{Q}_{\pm}}\coloneqq\mathbf{F}_{\pm}^{-1}$
introduced in [40] for distributions on $\mathbb{R}^{d}$, which involve
optimal transports to the spherical uniform over the unit ball. On the
hypersphere, the definition of the directional distribution function of a
random vector $\mathbf{Z}\sim{\rm P}^{\mathbf{Z}}\in\mathfrak{P}_{d}$ with
values on ${\cal S}^{d-1}$ is as follows (we indifferently use the expressions
directional distribution and quantile functions of $\mathbf{Z}$ and
directional distribution and quantile functions of ${\rm P}^{\mathbf{Z}}$).

######  Definition 3.1.

Call (directional) distribution function of ${\rm
P}^{\mathbf{Z}}\in\mathfrak{P}_{d}$ the a.s. unique optimal transport map
$\mathbf{F}$ from $\mathcal{S}^{d-1}$ to $\mathcal{S}^{d-1}$ such that
$\mathbf{F}\\#{\rm P}^{\mathbf{Z}}={\rm P}^{\bf U}$.

This definition is very natural, since $\mathbf{F}$ of $\mathbf{Z}\sim{\rm
P}^{\mathbf{Z}}$ is such that $\mathbf{F}(\mathbf{Z})\sim{\rm P}^{\bf U}$—so
that $\bf F$ has the same probability integral transformation flavor as
classical univariate distribution functions. We then have the following
important property (see the Appendix for a proof).

######  Proposition 3.2.

The distribution function $\mathbf{F}$ of $\mathbf{Z}\sim{\rm
P}^{\mathbf{Z}}\in\mathfrak{P}_{d}$ is a homeomorphism from
$\mathcal{S}^{d-1}$ to $\mathcal{S}^{d-1}$.

It directly follows from Proposition 3.2 that the inverse
$\mathbf{Q}\coloneqq\mathbf{F}^{-1}$ of $\mathbf{F}$ is well defined, which
yields the following definition of a directional quantile function.

######  Definition 3.3.

Call $\mathbf{Q}\coloneqq\mathbf{F}^{-1}$ the directional quantile function of
${\rm P}^{\mathbf{Z}}\in\mathfrak{P}_{d}$.

In $\mathbb{R}^{d}$, the center-outward distribution function
$\mathbf{F}_{\pm}$ was defined as an optimal transport to the spherical
uniform ${\rm U}_{d}$ over the unit ball of $\mathbb{R}^{d}$. For ${\rm
U}_{d}$, nested balls of the form $\tau{\mathcal{S}}^{d-1}$, centered at the
origin and with ${\rm U}_{d}$-probability contents $\tau$, naturally play the
role of quantile regions of order $\tau\in[0,1)$; the origin thus naturally
qualifies as the median of ${\rm U}_{d}$. The center-outward quantile regions
then are obtained as the images by
$\mathbf{Q}_{\pm}\coloneqq\mathbf{F}_{\pm}^{-1}$ of these balls.

The basic idea is quite similar in $\mathcal{S}^{d-1}$. If, however, nested
regions with ${\rm P}^{\bf U}$-probability contents $\tau\in[0,1]$ are to be
defined on $\mathcal{S}^{d-1}$, a central point or pole
${\boldsymbol{\theta}}\in\mathcal{S}^{d-1}$ has to be chosen, playing the role
of a directional median for ${\rm P}^{\bf Z}$. Since the uniform ${\rm P}^{\bf
U}$ automatically enjoys rotational symmetry with respect to the image ${\bf
F}({\boldsymbol{\theta}})$ of that pole, the nested regions playing the role
of ${\rm P}^{\bf U}$-quantile regions should be invariant with respect to
rotations with axis ${\bf F}({\boldsymbol{\theta}})$: hence, the collection of
spherical caps with ${\rm P}^{\bf U}$-probability contents $\tau$ centered at
${\bf F}({\boldsymbol{\theta}})$ naturally qualifies as the family of (nested)
quantile regions of orders $\tau\in~{}\\!\\![0,1]$ of ${\rm P}^{\bf U}$.

The choice of a directional median or pole ${\boldsymbol{\theta}}$, of course,
should be ${\rm P}^{\mathbf{Z}}$-specific. The literature on directional data
provides several concepts of directional medians, the most convenient of which
is the so-called Fréchet mean of $\mathbf{Z}\sim{\rm P}^{\mathbf{Z}}$, defined
as

${\boldsymbol{\theta}}_{\rm Fr}\coloneqq{\rm
argmin}_{\mathbf{z}\in\mathcal{S}^{d-1}}{\rm E}_{{\rm
P}^{\mathbf{Z}}}[c(\mathbf{Z},\mathbf{z})]$ (3.1)

(recall that $c(\mathbf{z}_{1},\mathbf{z}_{2})$ denotes the squared Riemannian
distance between $\mathbf{z}_{1}$ and $\mathbf{z}_{2}$). The Fréchet mean
always exists but is not necessarily unique. In case (3.1) has multiple
solutions, one can randomly select ${\boldsymbol{\theta}}$ as one of them.

Another possible choice can be based on the transport map $\mathbf{F}$
(equivalently, on the quantile function $\mathbf{Q}$) itself. For each
$\mathbf{z}\in{\cal S}^{d-1}$, consider the hyper-hemi-sphere ${\cal
H}_{\mathbf{z}}\coloneqq\\{{\bf u}\in{\cal S}^{d-1}:\mathbf{u}^{\top}{\bf
F}(\mathbf{z})\geq 0\\}$ centered at ${\bf F}(\mathbf{z})$. Since ${\bf F}$
pushes ${\rm P}^{\mathbf{Z}}$ forward to ${\rm P}^{\mathbf{U}}\\!$, we have
that ${\rm P}^{\mathbf{U}}({\cal H}_{\mathbf{z}})={\rm P}^{\mathbf{Z}}({\bf
Q}({\cal H}_{\mathbf{z}}))=1/2$. In general, however, ${\bf Q}({\cal
H}_{\mathbf{z}})$ is not a hyper-hemisphere and $1/2={\rm
P}^{\mathbf{U}}({\cal H}_{\mathbf{z}})\neq{\rm P}^{\mathbf{U}}({\bf Q}({\cal
H}_{\mathbf{z}}))$. A transport-based concept of median for ${\rm
P}^{\mathbf{Z}}$ is thus

${\boldsymbol{\theta}}_{\rm Tr}\coloneqq{\rm argmin}_{\mathbf{z}\in{\cal
S}^{d-1}}{\rm P}^{\mathbf{U}}({\bf Q}({\cal H}_{\mathbf{z}})).$ (3.2)

This directional median or pole ${\boldsymbol{\theta}}_{\rm Tr}$ indeed is
such that the area of the image by ${\bf Q}$ of the hyper-hemisphere centered
at ${\bf F}({\boldsymbol{\theta}}_{\rm Tr})$ is minimal among the areas of all
images ${\bf Q}({\cal H}_{\mathbf{z}})$ of hyper-hemispheres, indicating a
concentration of the ${\rm P}^{\mathbf{Z}}$ probability mass around
${\boldsymbol{\theta}}_{\rm Tr}$. This pole always exists. Again, it may not
be unique and one can randomly select ${\boldsymbol{\theta}}$ within the set
of points satisfying (3.2).

Assume that a pole ${\boldsymbol{\theta}}_{\rm M}$, say, has been selected for
${\rm P}^{\mathbf{Z}}$ (${\boldsymbol{\theta}}_{\rm Fr}$ or
${\boldsymbol{\theta}}_{\rm Tr}$ being possible choices). The spherical cap
with ${\rm P}^{\bf U}$-probability $\tau$ centered at ${\bf
F}({\boldsymbol{\theta}}_{\rm M})$ is

${\mathbb{C}}_{\tau}^{\mathbf{U}}\coloneqq\left\\{{\bf u}\in{\cal
S}^{d-1}:F_{*}\left({\bf u}^{\top}{\bf F}({\boldsymbol{\theta}}_{\rm
M})\right)\geq 1-\tau\right\\}\quad 0\leq\tau\leq 1,$

with boundary

${\mathcal{C}}_{\tau}^{\mathbf{U}}\coloneqq\left\\{{\bf u}\in{\cal
S}^{d-1}:F_{*}\left({\bf u}^{\top}{\bf F}({\boldsymbol{\theta}}_{\rm
M})\right)=1-\tau\right\\}\quad 0\leq\tau\leq 1,$ (3.3)

(a $(d-2)$-dimensional hypersphere) where

$F_{*}(u)\coloneqq{\int_{-1}^{u}(1-s^{2})^{(d-3)/2}\;{\mathrm{d}}s}{\Big{/}}{\int_{-1}^{1}(1-s^{2})^{(d-3)/2}\;{\mathrm{d}}s},\quad-1\leq
u\leq 1.$ (3.4)

Indeed, it is easy to see that $u\mapsto F_{*}(u)$ is the distribution
function of ${\bf U}^{\top}{\bf F}({\boldsymbol{\theta}}_{\rm M})$ where ${\bf
U}\sim{{\rm P}^{\mathbf{U}}}$. Hence, $F_{*}({\bf u}^{\top}{\bf
F}({\boldsymbol{\theta}}_{\rm M}))$ is the ${\rm P}^{\mathbf{U}}$-probability
of the ${\bf F}({\boldsymbol{\theta}}_{\rm M})$-centered spherical cap running
through $\mathbf{u}$ and a measure of $\mathbf{u}$’s latitude, ranging from 0
(for $\mathbf{u}=-{\bf F}({\boldsymbol{\theta}}_{\rm M})$) to 1 (for
$\mathbf{u}={\bf F}({\boldsymbol{\theta}}_{\rm M})$) and scaled in such a way
that the latitude of ${\bf U}\sim{{\rm P}^{\mathbf{U}}}$ is uniform over
$[0,1]$. Therefore, for the hypersphere ${\cal S}^{d-1}$ equipped with the
uniform distribution ${{\rm P}^{\mathbf{U}}}$ and the pole ${\bf
F}({\boldsymbol{\theta}}_{\rm M})$, the contour
${\mathcal{C}}_{\tau}^{\mathbf{U}}$ plays the role of a parallel of order
$\tau$ (${\mathcal{C}}_{1/2}^{\mathbf{U}}$, thus, is an equatorial
hypersphere).

Accordingly, define the quantile contour and quantile region of order $\tau$
of $\mathbf{Z}\sim{\rm P}^{\mathbf{Z}}$ as the images

${\mathcal{C}}_{\tau}\coloneqq{\mathbf{Q}}\left({\mathcal{C}}_{\tau}^{\mathbf{U}}\right)=\left\\{\mathbf{z}\in\mathcal{S}^{d-1}:F_{*}\left((\mathbf{F}(\mathbf{z})^{\top}\mathbf{F}({\boldsymbol{\theta}}_{\rm
M})\right)=1-\tau\right\\}\vspace{-2mm}$ (3.5)

and

${\mathbb{C}}_{\tau}\coloneqq{\mathbf{Q}}\left({\mathbb{C}}_{\tau}^{\mathbf{U}}\right)=\left\\{\mathbf{z}\in\mathcal{S}^{d-1}:F_{*}\left((\mathbf{F}(\mathbf{z}))^{\top}\mathbf{F}({\boldsymbol{\theta}}_{\rm
M})\right)\geq 1-\tau\right\\},\vspace{-2mm}$ (3.6)

by $\bf Q$ of ${\mathcal{C}}_{\tau}^{\mathbf{U}}$ and
${\mathbb{C}}_{\tau}^{\mathbf{U}}$, respectively. Since $\mathbf{Q}$ is a
measure-preserving transformation pushing ${\rm P}^{\mathbf{U}}$ forward to
${\rm P}^{\mathbf{Z}}$, the ${\rm P}^{\mathbf{Z}}$ probability content of
${\mathbb{C}}_{\tau}$ is $\tau$. The following proposition summarizes the main
properties of $\mathbf{F}$, $\mathbf{Q}$, the quantile contours
$\mathcal{C}_{\tau}$ and quantile regions $\mathbb{C}_{\tau}$. They all follow
from the definitions and the continuity of $\mathbf{F}$ and $\mathbf{Q}$.

######  Proposition 3.4.

Let $\mathbf{Z}\sim{\rm P}^{\mathbf{Z}}\in\mathfrak{P}_{d}$ have distribution
and quantile functions $\mathbf{F}$ and $\mathbf{Q}$, respectively. Then,

  1. (i)

$\mathbf{F}$ entirely characterizes ${\rm P}^{\mathbf{Z}}$,
$\mathbf{F}(\mathbf{Z})\sim\mathbf{F}\\#{\rm P}^{\mathbf{Z}}=~{}\\!{\rm
P}^{\mathbf{U}}$, and
$F_{*}\left((\mathbf{F}(\mathbf{Z}))^{\top}\mathbf{F}({\boldsymbol{\theta}}_{\rm
M})\right)\sim~{}\\!{\rm U}_{[0,1]}$;

  2. (ii)

$\mathbf{Q}$ entirely characterizes ${\rm P}^{\mathbf{Z}}$ and
$\mathbf{Q}(\mathbf{U})\sim\mathbf{Q}\\#{\rm P}^{\mathbf{U}}=~{}\\!{\rm
P}^{\mathbf{Z}}$ ;

  3. (iii)

the quantile contours $\mathcal{C}_{\tau}$, $\tau\in[0,1]$ are continuous; the
quantile regions ${\mathbb{C}}_{\tau}$ are closed, connected, and nested;
their intersection $\bigcap_{\tau\in[0,1]}\mathbb{C}_{\tau}$ is the
directional median ${\boldsymbol{\theta}}_{\rm M}$;

  4. (iv)

the probability content ${\rm P}^{\mathbf{Z}}\left(\mathbb{C}_{\tau}\right)$
of ${\mathbb{C}}_{\tau}$, $\tau\in[0,1]$, is $\tau$, irrespective of $\rm
P^{\mathbf{Z}}\in\mathfrak{P}_{d}$.

###  3.2 A natural distribution-specific coordinate system

The quantile function $\mathbf{Q}$ of ${\rm P}^{\mathbf{Z}}$ actually creates
over $\mathcal{S}^{d-1}$ a coordinate system with latitudes and
(hyper)longitudes adapted to the distribution of $\mathbf{Z}$.

The usual latitude-(hyper)longitude coordinate system is based on the
classical tangent-normal decomposition, with respect to a pole ${\bf
F}({\boldsymbol{\theta}}_{\rm M})$, say, of a point
$\mathbf{u}\in\mathcal{S}^{d-1}$ into the sum of two mutually orthogonal terms

${\bf u}=({\bf u}^{\top}{\bf F}({\boldsymbol{\theta}}_{\rm M})){\bf
F}({\boldsymbol{\theta}}_{\rm M})+({\bf I}_{p}-{\bf
F}({\boldsymbol{\theta}}_{\rm M})({\bf F}({\boldsymbol{\theta}}_{\rm
M}))^{\top}){\bf u}=({\bf u}^{\top}{\bf F}({\boldsymbol{\theta}}_{\rm M})){\bf
F}({\boldsymbol{\theta}}_{\rm M})+\sqrt{1-({\bf u}^{\top}{\bf
F}({\boldsymbol{\theta}}_{\rm M}))^{2}}{\bf S}_{{\bf
F}({\boldsymbol{\theta}}_{\rm M})\\!}({\bf u})$ (3.7)

(see the left panel of Figure 1 for $d=3$) where ${\bf u}^{\top}{\bf
F}({\boldsymbol{\theta}}_{\rm M})$, being constant over the parallel
$\mathcal{C}^{\mathbf{U}}_{\tau}$ ($\tau=1-F_{*}({\bf u}^{\top}{\bf
F}({\boldsymbol{\theta}}_{\rm M}))$) defined in Section 3.1, is a latitude
while the unit vector (a directional sign)

$\mathbf{S}_{{\bf F}({\boldsymbol{\theta}}_{\rm
M})}({\mathbf{u}})\coloneqq{(\mathbf{u}-(\mathbf{u}^{\top}{\bf
F}({\boldsymbol{\theta}}_{\rm M})){\bf F}({\boldsymbol{\theta}}_{\rm
M}))}/{\|\mathbf{u}-(\mathbf{u}^{\top}{\bf F}({\boldsymbol{\theta}}_{\rm
M})){\bf F}({\boldsymbol{\theta}}_{\rm M})\|}$ (3.8)

(with the convention ${\bm{0}}/0={\bm{0}}$ for $\mathbf{u}=\pm{\bf
F}({\boldsymbol{\theta}}_{\rm M})$) with values on the equatorial hypersphere
${\mathcal{C}}_{1/2}^{\mathbf{U}}$—a hyperlongitude, thus—characterizes
hypermeridians ${\mathcal{M}}_{\bf
s}^{\mathbf{U}}\coloneqq\left\\{\mathbf{u}\in\mathcal{S}^{d-1}:\mathbf{S}_{{\bf
F}({\boldsymbol{\theta}}_{\rm M})}({\mathbf{u}})={\bf s}\right\\}$, ${\bf
s}\in{\mathcal{C}}_{1/2}^{\mathbf{U}}$.

This parallel/hypermeridian system is well adapted to ${\mathcal{S}}^{d-1}$
when equipped with the surface area measure or the uniform measure ${\rm
P}^{\mathbf{U}}$: parallels, for instance, then coincide with ${\rm
P}^{\mathbf{U}}$’s quantile contours $\mathcal{C}^{\mathbf{U}}_{\tau}$, while
the hyperlongitudes $\mathbf{S}_{{\bf F}({\boldsymbol{\theta}}_{\rm
M})}({\mathbf{U}})$, for $\mathbf{U}\sim{\rm P}^{\mathbf{U}}$, are uniform
over the equatorial hypersphere. Its image by $\mathbf{Q}$ is more natural
under the probability measure ${\rm P}^{\mathbf{Z}}$: the same properties,
indeed, then hold for the “curvilinear” parallels
$\mathbf{Q}(\mathcal{C}^{\mathbf{U}}_{\tau})={\mathcal{C}}_{\tau}$ and the
“curvilinear” hypermeridians

${\mathcal{M}}_{\bf s}\coloneqq\mathbf{Q}({\mathcal{M}}_{\bf
s}^{\mathbf{U}})\coloneqq\left\\{\mathbf{z}\in\mathcal{S}^{d-1}:\mathbf{S}_{{\bf
F}({\boldsymbol{\theta}}_{\rm M})}(\mathbf{F}(\mathbf{z}))={\bf
s}\right\\},\quad{\bf s}\in{\mathcal{C}}_{1/2}^{\mathbf{U}}.\vspace{-4mm}$
(3.9)

${\bf F}({\boldsymbol{\theta}}_{\rm M})$${\bf F}({\bf z})$$(({\bf
F}({\boldsymbol{\theta}}_{\rm M}))^{\top}{\bf F}(\mathbf{z})){\bf
F}({\boldsymbol{\theta}}_{\rm M})$${\bf S}_{{\bf F}({\boldsymbol{\theta}}_{\rm
M})}{\bf F}(\mathbf{z})$

${{\boldsymbol{\theta}}_{{\rm M}}}$${\bf z}$${\bf F}({\bf
z})$$(\mathbf{z}^{\top}{\boldsymbol{\theta}}_{\rm
M}){\boldsymbol{\theta}}_{\rm
M}$$F_{f}^{*}(\mathbf{z}^{\top}{\boldsymbol{\theta}}_{\rm
M}){\boldsymbol{\theta}}_{\rm M}$${\bf S}_{{\boldsymbol{\theta}}_{\rm
M}}(\mathbf{z})$

Figure 1: Left: the tangent-normal decomposition of ${\bf F}$ with respect to
${\bf F}({\boldsymbol{\theta}}_{\rm M})$. Right: in the rotationally symmetric
case, the optimal transport $\mathbf{F}$ reduces to an optimal univariate
transport acting on the projection ${\bf
z}^{\top}{\boldsymbol{\theta}}_{\text{\rm M}}$ of $\mathbf{z}$ along the
rotation axis ${\boldsymbol{\theta}}_{\text{\rm M}}$.

###  3.3 Rotational symmetry

Optimal transports seldom admit closed-form expressions. In the particular
case of rotational symmetry, however, explicit expressions for ${\bf F}$ and
$\mathbf{Q}$ are possible.

Consider the rotationally symmetric variable $\mathbf{Z}\sim{\rm
P}^{\mathbf{Z}}$ with uniquely defined axis $\pm\,{\boldsymbol{\theta}}$; all
sensible directionalmedians ${\boldsymbol{\theta}}_{\rm M}$ (including
${\boldsymbol{\theta}}_{\rm Fr}$ and ${\boldsymbol{\theta}}_{\rm Tr}$) then
are lying along that axis, and ${\rm P}^{\mathbf{Z}}$ has density (1.2) with
${\boldsymbol{\theta}}={\boldsymbol{\theta}}_{\rm M}$: write ${\rm
P}^{\mathbf{Z}}\eqqcolon{\rm P}_{{\boldsymbol{\theta}}_{\rm M},f}$. Denote by

$F_{f}(u)\coloneqq{\int_{-1}^{u}f(s)(1-s^{2})^{(p-3)/2}\;{\mathrm{d}}s}{\big{/}}{\int_{-1}^{1}f(s)(1-s^{2})^{(p-3)/2}\;{\mathrm{d}}s},\quad-1\leq
u\leq 1,$

the distribution function of $\mathbf{Z}^{\top}{\boldsymbol{\theta}}_{\rm M}$,
by $Q_{f}\coloneqq F_{f}^{-1}$ its quantile function. Under rotational
symmetry, we have that

$\mathbf{S}_{{\boldsymbol{\theta}}_{\rm
M}}(\mathbf{\mathbf{Z}})={(\mathbf{Z}-(\mathbf{Z}^{\top}{\boldsymbol{\theta}}_{\rm
M}){\boldsymbol{\theta}}_{\rm
M})}/{\|\mathbf{Z}-(\mathbf{Z}^{\top}{\boldsymbol{\theta}}_{\rm
M}){\boldsymbol{\theta}}_{\rm M}\|}$ (3.10)

(with the convention that ${\bm{0}}/0={\bm{0}}$) is uniformly distributed over
the equatorial hypersphere

$\mathcal{S}^{d-2}_{{\boldsymbol{\theta}}^{\perp}_{\rm
M}}\coloneqq\left\\{\mathbf{u}\in\mathbb{R}^{d}:\mathbf{u}^{\top}{\boldsymbol{\theta}}_{\rm
M}=0,\mathbf{u}^{\top}\mathbf{u}=1\right\\}$ (3.11)

(see, e.g., [69]). Actually, $\mathbf{S}_{{\boldsymbol{\theta}}_{\rm
M}}(\mathbf{\mathbf{z}})$ is the point in the equatorial sphere with the same
hyperlongitudes as $\mathbf{z}$ and can be interpreted as the directional sign
of $\mathbf{z}$. Since $\mathbf{F}(\mathbf{Z})$, by definition, is uniformly
distributed over $\mathcal{S}^{d-1}$,
$(\mathbf{F}(\mathbf{Z}))^{\top}{\boldsymbol{\theta}}_{\text{\rm M}}$ has
distribution function $F_{*}$ as defined in (3.4) and quantile function
$Q_{*}=F_{*}^{-1}$. We then have an explicit form for the directional
distribution function ${\bf F}$ of $\mathbf{Z}\sim{\rm
P}_{{\boldsymbol{\theta}}_{\rm M},f}$ (see Appendix B for a proof).

######  Proposition 3.5.

Let $\mathbf{Z}$ have rotationally symmetric distribution ${\rm
P}_{{\boldsymbol{\theta}}_{\rm M},f}$. Then, letting $F_{{f}}^{*}(u)\coloneqq
Q_{*}(F_{f}(u))$,

$\mathbf{F}(\mathbf{z})=F_{f}^{*}({\bf z}^{\top}{\boldsymbol{\theta}}_{\rm
M}){\boldsymbol{\theta}}_{\rm M}+\sqrt{1-(F^{*}_{{f}}({\bf
z}^{\top}{\boldsymbol{\theta}}_{\rm
M}))^{2}}\,\mathbf{S}_{{\boldsymbol{\theta}}_{\rm M}}(\mathbf{z}).$ (3.12)

Proposition 3.5 tells us that, in the rotationally symmetric case, the optimal
transport $\mathbf{F}$ essentially reduces to the optimal univariate
transport, the function $F_{{f}}^{*}$ acting on the projection ${\bf
z}^{\top}{\boldsymbol{\theta}}_{\text{\rm M}}$ of $\mathbf{z}$ onto the axis
${\boldsymbol{\theta}}_{\text{\rm M}}$. This is illustrated in Figure 1 (right
panel). Again, this explains—but does not justify—the ubiquity, in directional
statistics, of the rotational symmetry assumption: just as the assumption of
ellipticity in $\mathbb{R}^{d}$, rotational symmetry, when it holds true,
actually reduces complex multivariate problems to much simpler univariate
ones.

##  4 Empirical distribution and quantile functions, ranks,
and signs

So far, we have been dealing with population concepts of distribution and
quantile functions. If a statistical analysis is to be performed, we need
empirical versions of the same; these involve concepts of ranks and signs,
which are completely new in the context of directional observations.

###  4.1 Empirical directional distribution function

Denoting by $\mathbf{Z}_{1}^{(n)},\ldots,\mathbf{Z}_{n}^{(n)}$ a sample of $n$
i.i.d. observations with distribution ${\rm P}^{\mathbf{Z}}$, consider a
“regular” $n$-point grid
$\mathfrak{G}^{(n)}=~{}\\!\\{{\scriptstyle{\mathfrak{G}}}^{(n)}_{1},\ldots,{\scriptstyle{\mathfrak{G}}}^{(n)}_{n}\\}$
over the unit hypersphere ${\mathcal{S}}^{d-1}$. In this section, the only
property required from $\mathfrak{G}^{(n)}$ is that the sequence of uniform
distributions over the $n$-tuple
${\scriptstyle{\mathfrak{G}}}^{(n)}_{1},\ldots,{\scriptstyle{\mathfrak{G}}}^{(n)}_{n}$
of gridpoints converges weakly to the uniform distribution ${\rm P}^{\bf U}$
over ${\mathcal{S}}^{d-1}$ as $n\to\infty$.

Since $\mathbf{F}$ minimizes the transportation cost (2.1), a plugin estimator
is obtained as the solution of an optimal coupling problem between the
observations and the grid $\mathfrak{G}^{(n)}$. More precisely, let
$\mathcal{T}^{(n)}$ denote the set of all permutations of the integers
$\\{1,\ldots,n\\}$: each permutation $T\in~{}\\!\mathcal{T}^{(n)}$ defines a
bijection
$\mathbf{Z}^{(n)}_{i}\mapsto{\scriptstyle{\mathfrak{G}}}^{(n)}_{T(i)}$ between
$\\{\mathbf{Z}_{1}^{(n)},\ldots,\mathbf{Z}_{n}^{(n)}\\}$ and
$\mathfrak{G}^{(n)}$. The empirical directional distribution function
$\mathbf{F}^{(n)}$ is then defined as the mapping

$\mathbf{F}^{(n)}:\mathbf{Z}^{(n)}\coloneqq(\mathbf{Z}_{1}^{(n)},\ldots,\mathbf{Z}_{n}^{(n)})\mapsto(\mathbf{F}^{(n)}(\mathbf{Z}_{1}^{(n)}),\ldots,\mathbf{F}^{(n)}(\mathbf{Z}_{n}^{(n)}))$
(4.1)

satisfying (with $c(\mathbf{z}_{1},\mathbf{z}_{2})$ the squared Riemannian
distance)

$\sum_{i=1}^{n}c(\mathbf{Z}_{i}^{(n)},\mathbf{F}^{(n)}(\mathbf{Z}_{i}^{(n)}))=\underset{T\in\mathcal{T}^{(n)}}{\min}\sum_{i=1}^{n}c(\mathbf{Z}_{i}^{(n)},{\scriptstyle{\mathfrak{G}}}^{(n)}_{T(i)}).$
(4.2)

The terminology empirical directional distribution function is justified by
the following Glivenko-Cantelli result (see Appendix B for a proof).

######  Proposition 4.1.

(Glivenko-Cantelli). Let $\mathbf{Z}_{1}^{(n)},\ldots,\mathbf{Z}_{n}^{(n)}$ be
i.i.d. with distribution ${\rm P}^{\mathbf{Z}}$ over $\mathcal{S}^{d-1}$.
Then, provided that the sequence of uniform discrete distributions over the
$n$ gridpoints
$\\{{\scriptstyle{\mathfrak{G}}}^{(n)}_{1},\ldots,{\scriptstyle{\mathfrak{G}}}^{(n)}_{n}\\}$
of ${\mathfrak{G}}^{(n)}$ converges weakly, as $n\to\infty$, to the uniform
distribution ${\rm P}^{\mathbf{U}}$ over $\mathcal{S}^{d-1}$,

${\rm max}_{1\leq i\leq
n}\big{\|}\mathbf{F}^{(n)}(\mathbf{Z}_{i}^{(n)})-\mathbf{F}(\mathbf{Z}_{i}^{(n)})\big{\|}\longrightarrow
0\qquad\text{almost surely as $n\rightarrow\infty$. }$ (4.3)

If a consistent estimation of $\mathbf{F}$ is the objective, the empirical
distribution $\mathbf{F}^{(n)}$ defined in (4.1), in view of (4.3), offers a
perfect solution. If empirical counterparts of the quantile contours and
regions $\mathcal{C}_{\tau}$ and $\mathbb{C}_{\tau}$, hence of the parallels
${\mathcal{C}}_{\tau}$ and hypermeridians ${\mathcal{M}}_{\bf s}$ (see (3.3)
and (3.9)) of Section 3.2 are to be constructed, or if ranks and signs
generating a maximal ancillary sigma-field are to be defined, however, more
structure is required from the grid $\mathfrak{G}^{(n)}$.

###  4.2 Directional ranks, signs, and empirical quantiles

In this subsection, we show how imposing some additional finite-$n$ structure
on the grid $\mathfrak{G}^{(n)}$ yields empirical versions of $\mathbf{F}$
with natural data-driven concepts of

  1. (i)

empirical quantile contours and regions consistently estimating the actual
ones,

  2. (ii)

a data-driven coordinate system of empirical parallels and hypermeridians
adapting to the underlying distribution ${\rm P}^{\mathbf{Z}}$ of the
observations, and

  3. (iii)

distribution-free signs and ranks, paving the way to a theory of rank-based
inference for directional data with unspecified density.

All these concepts involve a pole ${\boldsymbol{\theta}}_{\rm M}$ (the
population Fréchet mean, for instance) or a consistent estimator
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}_{\rm
M}{\phantom{\theta}}^{{(n)}}$ thereof (e.g., the empirical Fréchet mean). That
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}_{\rm
M}{\phantom{\theta}}^{{(n)}}$ in turn will define the pole
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}\\!$
of a structured grid
$\mathfrak{G}^{(n)}({{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}}\\!)$
to be used in the estimation of $\mathbf{F}$. Accordingly, we are proceeding
in two steps: a first step yielding an empirical pole
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}\\!$,
then a second step estimating $\mathbf{F}$ on the basis of the grid
$\mathfrak{G}^{(n)}({\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}})$
exploiting the role of
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$ as
an empirical pole.

Step 1. Construct an empirical version
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$ of
the population quantity ${\bf F}({\boldsymbol{\theta}}_{\rm M})$. Let
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}_{\rm
M}{\phantom{\theta}}^{{(n)}}$ be a consistent estimator of
${\boldsymbol{\theta}}_{\rm M}$ and consider an $n$-point grid
$\mathfrak{G}^{(n)}_{0}$ satisfying the assumptions of Proposition 4.1;
denoting by ${\bf F}^{(n)}_{0}(\mathbf{Z}_{i}^{(n)})$ the resulting estimator
of ${\bf F}(\mathbf{Z}_{i}^{(n)})$, choose
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}\coloneqq{\bf
F}^{(n)}_{0}(\mathbf{Z}_{i}^{(n)})$ where $\mathbf{Z}_{i}^{(n)}$ is the sample
point closest to ${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}_{\rm
M}{\phantom{\theta}}^{{(n)}}$ in the $d(\cdot,\cdot)$ metric: Proposition 4.1
and the continuity of $\mathbf{F}$ imply that
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$ a.s.
converges to ${\bf F}({\boldsymbol{\theta}}_{\rm M})$.

Step 2. (2a) Construct a further regular grid
$\mathfrak{G}^{(n)}({{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}})$
over ${\cal S}^{d-1}$. Factorizing $n$ into $n=n_{R}n_{S}+n_{0}$ where
$n_{R},n_{S},n_{0}\in\mathbb{N}$ and $0\leq n_{0}<\min\\{n_{R},n_{S}\\},$
define that grid
$\mathfrak{G}^{(n)}({{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}})$
as the product of two independent grids:

  1. (i)

a reference grid
$\mathfrak{S}^{(n_{S})}\coloneqq\\{{\mathbf{s}_{1},\ldots,\mathbf{s}_{n_{S}}}\\}$
over ${\mathcal{S}}^{d-2}$; again, this grid should be as uniform as possible
but the only requirement is the weak convergence, as $n_{S}\to~{}\\!\infty$,
of the uniform distribution over $\mathfrak{S}^{(n_{S})}$ to the uniform over
${\mathcal{S}}^{d-2}$. Note that for $d=~{}\\!3$, a fully regular grid
$\mathbf{s}_{1},\ldots,\mathbf{s}_{n_{S}}$ is obtained by dividing the unit
circle ${\cal S}^{1}$ into $n_{S}$ equal parts (see Figure 2);

  2. (ii)

a grid of $n_{R}$ points over the unit interval, of the form
${i}/{(n_{R}+1)}$, $i=1,\ldots,n_{R}$.

${\bf s}_{1}=(1,0)^{\top}$${\bf s}_{2}$${\bf
s}_{n_{S}}$$\bullet$$\bullet$$\bullet$$1/(n_{R}+1)$$2/(n_{R}+1)$$\vdots$$n_{R}/(n_{R}+1)$01

${\bf F}^{(n)}(\widehat{{\boldsymbol{\theta}}}_{\rm M})$grid of $n_{S}$ points
determiningthe meridians (the longitudes)grid of $n_{R}$ points determiningthe
parallels (the latitudes)

Figure 2: The grid used to define signs and ranks on ${\cal S}^{d-1}$ for
$d=3$. The final grid is obtained as the product of a reference grid
$\mathfrak{S}^{(n_{S})}\coloneqq\\{{\mathbf{s}_{1},\ldots,\mathbf{s}_{n_{S}}}\\}$
over ${\mathcal{S}}^{d-2}$ (here, $n_{S}$ equispaced points on the circle
${\mathcal{S}}^{1}$) and $n_{R}$ equispaced points on the unit interval.

(2b) Import the grid $\mathfrak{S}^{(n_{S})}$ of step (2a) to the equatorial
space defined by the pole
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$
computed in step 1. More precisely, construct the grid
$\mathfrak{G}^{(n)}({{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}})$
that consists in $n_{0}$ copies of
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$ (if
$n_{0}\neq 0$) and the $n_{R}n_{S}$ points
${\scriptstyle{\mathfrak{G}}}^{(n)}_{ij}({{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}})$
such that

$1-F_{*}(({\scriptstyle{\mathfrak{G}}}^{(n)}_{ij}({{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}}))^{\top}{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}})=\frac{i}{n_{R}+1}\quad\text{
and
}\quad\mathbf{S}_{\,{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}}({\scriptstyle{\mathfrak{G}}}^{(n)}_{ij}({{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}}))={\boldsymbol{\Gamma}}_{{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}}{{\bf
s}}_{j},\quad i=1,\ldots,n_{R},\,j=1,\ldots,n_{S}\vspace{-1mm}$

where ${\boldsymbol{\Gamma}}_{\boldsymbol{\theta}}$ denotes a $d\times(d-1)$
semi-orthogonal matrix such that
${\boldsymbol{\Gamma}}_{\boldsymbol{\theta}}{\boldsymbol{\Gamma}}_{\boldsymbol{\theta}}^{\top}={\bf
I}_{d}-{\boldsymbol{\theta}}{\boldsymbol{\theta}}^{\top}$ and
${\boldsymbol{\Gamma}}_{\boldsymbol{\theta}}^{\top}{\boldsymbol{\Gamma}}_{\boldsymbol{\theta}}={\bf
I}_{d-1}$; the columns of
$({\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}},{\boldsymbol{\Gamma}}_{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}})$,
thus, constitute an orthonormal coordinate system of ${\mathbb{R}}^{d}$, the
columns of
${\boldsymbol{\Gamma}}_{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$
an arbitrary orthonormal coordinate system of the equatorial hyperplane
determined by
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}\\!$.
Note that, due to its dependence on
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}\\!$,
this grid
$\mathfrak{G}^{(n)}({{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}})$
is random.

(2c) Denote by
$\mathbf{F}^{(n)}(\mathbf{Z}_{1}^{(n)}),\ldots,\mathbf{F}^{(n)}(\mathbf{Z}_{n}^{(n)})$
the solutions of the optimal coupling problem (4.2) based on this second grid
$\mathfrak{G}^{(n)}({{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}})$.

The (directional) empirical sign ${\bf S}_{i}^{(n)}$ and the (directional)
rank $R_{i}^{(n)}$ of $\mathbf{Z}_{i}^{(n)}$ then are naturally defined as

${\bf S}_{i}^{(n)}\coloneqq{\bf
S}_{\,{\,\widehat{\boldsymbol{\theta}}}^{\;\phantom{n}}{\phantom{\theta}}^{{(n)}}}({\bf
F}^{(n)}(\mathbf{Z}_{i}^{(n)}))\vspace{-4.5mm}$ (4.4)

and

$R^{(n)}_{i}=R^{(n)}(\mathbf{Z}_{i}^{(n)})\coloneqq(n_{R}+1)\left[1-F_{*}(({\bf
F}^{(n)}(\mathbf{Z}_{i}^{(n)}))^{\top}{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}))\right],\quad
i=1,\ldots,n\vspace{-2.5mm}$ (4.5)

(with values in $\\{1,\ldots,n_{R}\\}$), respectively, provided that ${\bf
F}^{(n)}(\mathbf{Z}_{i}^{(n)})\neq{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$;
for $\mathbf{Z}_{i}^{(n)}$ such that ${\bf
F}^{(n)}(\mathbf{Z}_{i}^{(n)})={\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$,
let ${\bf S}_{i}^{(n)}\coloneqq{\bm{0}}$ (the pole has no specific longitude)
and $R^{(n)}_{i}\coloneqq 0\vspace{0mm}$.

The empirical versions of the quantile contours and regions (3.5) and (3.6)
are the collections of observations

$\mathcal{C}^{(n)}_{j/(n_{R}+1)}\\!\coloneqq\\!\left\\{\mathbf{Z}_{i}^{(n)}\\!\\!:R_{i}^{(n)}\\!=j\right\\}\quad\\!\\!\\!\text{and}\\!\quad\\!\\!\mathbb{C}^{(n)}_{j/(n_{R}+1)}\\!\coloneqq\\!\left\\{\mathbf{Z}_{i}^{(n)}\\!\\!:R_{i}^{(n)}\leq
j\right\\},\quad\\!j=1,\ldots,n_{R},$ (4.6)

respectively. An empirical quantile contour of order $j/(n_{R}+1)$ (the
empirical parallel of order $j/(n_{R}+1)$) thus consists of the $n_{S}$ sample
points with given rank $j$. Similarly, the empirical hypermeridian with
longitude ${\bf s}$ consists of the $n_{R}$ observations with given sign
(hyperlongitude) ${\bf s}$:

${\mathcal{M}}^{(n)}_{{\bf s}}\\!\coloneqq\\!\left\\{\mathbf{Z}_{i}^{(n)}\\!:\
\mathbf{F}^{(n)}(\mathbf{Z}_{i}^{(n)})\neq{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}\text{
and }{\bf S}_{i}^{(n)}={\bf s}\right\\},\quad{{\bf
s}}\in\\{{\boldsymbol{\Gamma}}_{{\,\widehat{\boldsymbol{\theta}}}^{\;\phantom{n}}{\phantom{\theta}}^{{(n)}}}{{\bf
s}}_{1},\ldots,{\boldsymbol{\Gamma}}_{{\,\widehat{\boldsymbol{\theta}}}^{\;\phantom{n}}{\phantom{\theta}}^{{(n)}}}{{\bf
s}}_{n_{S}}\\}.$ (4.7)

Because of the data-driven choice of
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$, the
grid
$\mathfrak{G}^{(n)}({{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}})$
is random. If $n_{R}\to\infty$ and $n_{S}\to\infty$, the assumptions in
Proposition 4.1 are satisfied a.s., and the Glivenko-Cantelli result (4.3)
holds a.s. as well. Similar consistency properties then follow for quantile
contours and regions, parallels, and hypermeridians. In Supplementary Material
A, we establish and discuss several properties—distribution-freeness,
ancillarity, and equivariance—of directional signs and ranks.

###  4.3 Numerical illustration: simulations

To illustrate the concept of quantile contours in (4.6), we generated
sequences of $n=2001$ i.i.d. unit random vectors from three types of
distributions on $\mathcal{S}^{2}$:

  1. (i)

the von Mises-Fisher (vMF) distribution (see (1.1)) with concentration
$\kappa=10$ and location ${\boldsymbol{\theta}}=(0,0,1)^{\top}$. Below,
$\mathcal{M}_{d}({\boldsymbol{\theta}},\kappa)$ denotes the von Mises
distribution on ${\cal S}^{d-1}$ with location ${\boldsymbol{\theta}}$ and
concentration $\kappa$;

  2. (ii)

the tangent vMF distribution as defined in [33]. The tangent vMF distribution
with location ${\boldsymbol{\theta}}$, angular function $G$, skewness
direction $\mu$, and skewness intensity $\kappa$ is the distribution of

$\mathbf{Z}\coloneqq
V{\boldsymbol{\theta}}+\sqrt{1-V^{2}}{\bm{\Gamma}}_{{\boldsymbol{\theta}}}{\bf
U},\vspace{-2mm}$ (4.8)

where ${\boldsymbol{\Gamma}}_{\boldsymbol{\theta}}$ is the $d\times(d-1)$
semi-orthogonal matrix described in Part (2b) of Step 2 of the construction
(Section 4.2) of the grid, $V$ (an absolutely continuous scalar random
variable with values in [-1,1]) and ${\bf
U}\sim\mathcal{M}_{2}({\mbox{\boldmath$\mu$\unboldmath}},\kappa)$ are mutually
independent; in the simulation, we set ${\boldsymbol{\theta}}=(0,0,1)^{\top}$,
${\mbox{\boldmath$\mu$\unboldmath}}=(0.7,\sqrt{0.51})^{\top}$, $\kappa=10$,
and $V=2\widetilde{V}-1$ with $\widetilde{V}\sim\text{\rm Beta}(2,8)$;

  3. (iii)

a mixture of two vMF distributions—the distribution of $I[U\leq
0.3]\mathbf{Z}_{1}+I[U>0.3]\mathbf{Z}_{2},$ where $U\sim{\rm U}[0,1]$, $I[\>]$
denotes the indicator function, and
$\mathbf{Z}_{1}\sim\mathcal{M}_{3}({\boldsymbol{\theta}}_{1},\kappa_{1})$ and
$\mathbf{Z}_{2}\sim\mathcal{M}_{3}({\boldsymbol{\theta}}_{2},\kappa_{2})$ are
mutually independent. In the simulation, we set $\kappa_{1}=20$,
${\boldsymbol{\theta}}_{1}=(0,-0.5,\sqrt{0.75})^{\top}$, $\kappa_{2}=20$, and
${\boldsymbol{\theta}}_{2}=(0,0.5,\sqrt{0.75})^{\top}$.

For each simulation scheme, we computed the Fréchet mean via the rbase.robust
function from the R package RiemBase. The optimal coupling between the sample
points and the regular grid has been obtained by using the fast network
simplex algorithm (FNSA) of [18] as implemented in the transport function of
the R package transport. Many efficient algorithms have been proposed for
solving the optimal transport problem, e.g., auction algorithm and its refined
version [13, 14]. We used the FNSA which is reported in [18] to reduce the
complexity of the optimal transport problem from $O(n^{3})$ to $O(n^{2})$. We
factorized $n$ into $n_{R}n_{S}+n_{0}$ with $n_{R}=40$, $n_{S}=50$ and
$n_{0}=1$: each empirical quantile contour then consists in 50 points. Plots
of the median
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$, the
empirical quantile contours for probability contents 12.2% (rank = 5), 48.8%
(rank = 20), and 70.7% (rank = 29) are shown in Figures 3 (left, middle, and
right panels for the simulation schemes (i), (ii), and (iii), respectively).
We also provide in the same figures the corresponding empirical meridians
(points with the same signs, same color in the figures).

Figure 3: Upper panel: empirical quantile contours (probability contents
12.2%, 48.8%, and 70.7%, respectively) computed from $n=2001$ ($n_{R}=40$,
$n_{S}=50$ and $n_{0}=1$) points drawn from a von Mises-Fisher distribution
(left), a mixture of two von Mises-Fisher distributions (middle) and a tangent
von Mises-Fisher distribution (right). Lower panel: the corresponding
empirical meridians; points with the same color have the same signs.

Under the vMF distribution (Figure 3, left panel), the empirical quantile
contours are symmetrically distributed around the median, as expected from the
rotational symmetry of the vMF distribution. Under the mixture distribution
(Figure 3, middle panel), quantile contours adapt to the underlying
multimodality, and the median is located in the mixture component with larger
probability weight. Under the tangent vMF distribution, where skewness is
involved (Figure 3, right panel), the empirical quantile contours exhibit a
distinctive skewed shape. Our empirical quantile contours thus nicely pick up
the shapes of the underlying distributions while controlling the probability
contents of the corresponding quantile regions.

##  5 Directional Goodness-of-fit

Classical goodness-of-fit tests (Kolmogorov-Smirnov, Cramér-von Mises, etc.)
are based on distances between distribution functions. Our concepts of
population and empirical directional distribution functions quite naturally
lead to the construction of directional versions based on distances between
${\bf F}^{(n)}$ and ${\bf F}$. Assuming that a sample of i.i.d. random unit
vectors $\mathbf{Z}_{1}^{(n)},\ldots,\mathbf{Z}_{n}^{(n)}$ have common
distribution ${\rm P}^{\mathbf{Z}}$ the goodness-of-fit problem consists in
testing the null hypothesis ${\cal H}_{0}:{\rm P}^{\mathbf{Z}}={\rm P}_{0}$
against ${\cal H}_{1}:{\rm P}^{\mathbf{Z}}\neq{\rm P}_{0}$, where ${\rm
P}_{0}$ is some specified distribution. This has been studied extensively in
directional statistics, with special attention to the problem of testing
uniformity—in which case ${\rm P}_{0}$ is the uniform distribution ${\rm
P}^{\bf U}$ over $\mathcal{S}^{d-1}$. In view of Proposition 3.4, testing
${\cal H}_{0}\\!:\\!{\rm P}^{\mathbf{Z}}\\!=~{}\\!\\!{\rm P}_{0}$ against
${\cal H}_{1}\\!:\\!{\rm P}^{\mathbf{Z}}\\!\neq~{}\\!\\!{\rm P}_{0}$ is
equivalent to testing ${\cal H}_{0}\\!:\\!{\bf F}\\!=\\!{\bf F}_{0}$ against
${\cal H}_{1}\\!:~{}\\!\\!{\bf F}\\!\neq~{}\\!\\!{\bf F}_{0}$ where ${\bf
F}_{0}$ denotes the distribution function of ${\rm P}_{0}$—i.e., ${\bf
F}_{0}(\mathbf{z})=\mathbf{z}$ in case ${\rm P}_{0}={\rm P}^{\bf U}$.

The test we are proposing in this Section is a Cramér-von Mises-type test that
rejects the null hypothesis for large values of the test statistic

$T_{n}\coloneqq n^{-1}\sum_{i=1}^{n}\left\|{\bf
F}^{(n)}(\mathbf{Z}_{i}^{(n)})-{\bf
F}_{0}(\mathbf{Z}_{i}^{(n)})\right\|^{2},\vspace{-2mm}$ (5.1)

where ${\bf F}^{(n)}$ is obtained as in Section 4 (the construction considered
in Section 4.1 here is sufficient). That test statistic is the empirical
counterpart of the (squared) $L_{2}$ distance ${\rm E}\left[\|{\bf
F}(\mathbf{Z}_{i}^{(n)})-{\bf F}_{0}(\mathbf{Z}_{i}^{(n)})\|^{2}\right].$ As
always in goodness-of-fit tests, the critical value $c_{\alpha}$ such that
${\rm P}[T_{n}>c_{\alpha}]=\alpha$ under ${\cal H}_{0}$ is easily approximated
via Monte Carlo simulations since the null hypothesis is simple; the resulting
test thus has exact size $\alpha$. The following result establishes its
consistency.

######  Proposition 5.1.

Assume that $\mathbf{Z}_{1}^{(n)},\ldots,\mathbf{Z}_{n}^{(n)}$ are i.i.d. with
common optimal transport map ${\bf F}$. Then, provided that the uniform
discrete distribution over the $n$-points grid ${\mathfrak{G}}^{(n)}$
converges weakly, as $n\to\infty$, to the uniform distribution ${\rm
P}^{\mathbf{U}}$ over $\mathcal{S}^{d-1}$,

  1. (i)

$T_{n}=o_{\rm P}(1)$ as $n\rightarrow\infty$ under ${\cal H}_{0}$, while

  2. (ii)

$T_{n}$ converges in probability to a strictly positive constant as
$n\rightarrow\infty$ if ${\bf F}\neq{\bf F}_{0}$.

It follows from Proposition 5.1 that the test $\phi_{n}$ is asymptotically
consistent against any fixed alternative and therefore qualifies as a
universally consistent omnibus test.

Simulation for $d=3$. In order to investigate the finite-sample performances
of our test of uniformity based on (5.1), we performed Monte Carlo size and
power comparisons of our test and the projected Cramér-von Mises (PCvM),
projected Anderson-Darling(PAD), and projected Rothman (PRt) tests recently
proposed in [32] and the Rayleigh, Bingham, Ajne, Giné and Bakshaev tests of
uniformity [75, 17, 2, 36, 5]; these tests were implemented in the R package
sphunif. To do so, we generated $N=1000$ independent samples of $n=400$ i.i.d.
unit random vectors with the following distributions:

  1. (i)

the uniform distribution;

  2. (ii)

vMF distributions, all with location parameter
${\boldsymbol{\theta}}=(0,0,1)^{\top}$, with concentration parameters
$\kappa=0.05$, $\kappa=0.1$, and $\kappa=0.5$;

  3. (iii)

tangent vMF distributions as in (4.8), all with angular distribution $F_{*}$,
location ${\boldsymbol{\theta}}=(0,0,1)^{\top}$, skewness direction
${\mbox{\boldmath$\mu$\unboldmath}}=(0,1)^{\top}$, and skewness intensities
$\kappa=0.05$, $\kappa=0.1$, and $\kappa=0.2$;

  4. (iv)

mixtures of two vMFs, of the form $I[U\leq
0.5]\mathbf{Z}_{1}+I[U>0.5]\mathbf{Z}_{2}$, where $U\sim{\rm U}[0,1]$,
$\mathbf{Z}_{1}\sim\mathcal{M}_{3}({\boldsymbol{\theta}}_{1},\kappa_{1})$, and
$\mathbf{Z}_{2}\sim\mathcal{M}_{3}({\boldsymbol{\theta}}_{2},\kappa_{2})$ are
mutually independent, with
${\boldsymbol{\theta}}_{1}=(0,-0.3,\sqrt{0.91})^{\top}$ and
${\boldsymbol{\theta}}_{2}=(0.3,\sqrt{0.66},0.5)^{\top}$, and concentration
parameters $\kappa_{1}=\kappa_{2}=0.1$, $\kappa_{1}=\kappa_{2}=0.2$, and
$\kappa_{1}=\kappa_{2}=0.3$;

  5. (v)

mixtures of two vMFs and a tangent vMF, of the form $I[U\leq
0.5]\mathbf{Z}_{1}+I[0.5<U<0.75]\mathbf{Z}_{2}+I[U\geq 0.75]\mathbf{Z}_{3}$,
where $U\sim{\rm U}[0,1]$, $\mathbf{Z}_{1}$, $\mathbf{Z}_{2}$ and
$\mathbf{Z}_{3}$ are mutually independent;, with ${\mathbf{Z}}_{1}$ following
the tangent vMF distribution as in (iii) with skewness intensities
$\kappa_{1}=0.07$, $\kappa_{1}=0.1$, and $\kappa_{1}=0.2$, and
$\mathbf{Z}_{2}\sim\mathcal{M}_{p}({\boldsymbol{\theta}}_{2},\kappa_{2})$ and
$\mathbf{Z}_{3}\sim\mathcal{M}_{p}({\boldsymbol{\theta}}_{3},\kappa_{3})$ as
in scheme (iv) with concentration parameters $\kappa_{2}=\kappa_{3}=0.07$,
$\kappa_{2}=\kappa_{3}=0.1$; and $\kappa_{2}=\kappa_{3}=0.2$.

In Table 1, we report the rejection frequencies (at nominal level
$\alpha=0.05$), out of $N=~{}\\!1000$ replications, of the optimal transport-
based Cramér-von Mises test described in (5.1) (OT), the PCvM, PAD, PRt,
Rayleigh, Bingham, Ajne, Giné, and Bakshaev tests of uniformity. The critical
values are obtained through $2000$ Monte Carlo replications. Inspection of
Table 1 reveals that all rejection frequencies under the null hypothesis of a
uniform distribution are close to the nominal level $0.05$. Under the
rotationally symmetric vMF alternatives, all tests show similar performances,
except for the Bingham and Giné tests, which have significantly lower powers.
Under non-rotationally symmetric alternatives such as the tangent vMF or
mixture distributions, our test uniformly outperforms all its competitors.

Table 1: Rejection frequencies of the OT, PCvM, PAD, PRt, Rayleigh, Bingham, Ajne, Giné and Bakshaev tests of uniformity over $\mathcal{S}^{2}$ under the uniform, vMF, tangent vMF, mixtures of two vMF distributions and mixtures of two vMFs and a tangent vMF, with various values of concentration or intensity parameter (simulation settings: $N=1000$ replications, sample size $n=400$) | OT | PCvM | PAD | PRt | Rayleigh | Bingham | Ajne | Giné | Bakshaev
---|---|---|---|---|---|---|---|---|---
Uniform | 0.045 | 0.054 | 0.054 | 0.054 | 0.053 | 0.060 | 0.051 | 0.059 | 0.054
vMF ($\kappa=0.05$) | 0.078 | 0.076 | 0.073 | 0.072 | 0.070 | 0.052 | 0.072 | 0.053 | 0.076
vMF ($\kappa=0.1$) | 0.134 | 0.137 | 0.139 | 0.140 | 0.139 | 0.053 | 0.137 | 0.055 | 0.137
vMF ($\kappa=0.5$) | 0.990 | 0.998 | 0.998 | 0.998 | 0.998 | 0.082 | 0.998 | 0.082 | 0.998
tangent vMF ($\kappa=0.05$) | 0.073 | 0.042 | 0.042 | 0.044 | 0.050 | 0.031 | 0.051 | 0.03 | 0.042
tangent vMF ($\kappa=0.1$) | 0.202 | 0.138 | 0.133 | 0.138 | 0.137 | 0.027 | 0.143 | 0.029 | 0.138
tangent vMF ($\kappa=0.2$) | 0.669 | 0.578 | 0.575 | 0.573 | 0.580 | 0.030 | 0.586 | 0.028 | 0.578
| Mixture of two vMFs
---
($\kappa_{1}=\kappa_{2}=0.1$)
0.122 | 0.099 | 0.100 | 0.098 | 0.105 | 0.055 | 0.100 | 0.059 | 0.099
| Mixture of two vMFs
---
($\kappa_{1}=\kappa_{2}=0.2$)
0.325 | 0.303 | 0.302 | 0.297 | 0.295 | 0.058 | 0.297 | 0.061 | 0.303
| Mixture of two vMFs
---
($\kappa_{1}=\kappa_{2}=0.3$)
0.628 | 0.615 | 0.611 | 0.615 | 0.621 | 0.042 | 0.624 | 0.046 | 0.615
| Mixture of two vMFs
---
and a tangent vMF
($\kappa_{1}=\kappa_{2}=\kappa_{3}=0.07$)
0.117 | 0.108 | 0.112 | 0.108 | 0.102 | 0.074 | 0.104 | 0.077 | 0.108
| Mixture of two vMFs
---
and a tangent vMF
($\kappa_{1}=\kappa_{2}=\kappa_{3}=0.1$)
0.272 | 0.203 | 0.205 | 0.206 | 0.201 | 0.064 | 0.202 | 0.061 | 0.203
| Mixture of two vMFs
---
and a tangent vMF
($\kappa_{1}=\kappa_{2}=\kappa_{3}=0.2$)
0.707 | 0.609 | 0.612 | 0.606 | 0.598 | 0.083 | 0.604 | 0.082 | 0.609

Simulation for $d=2$. We investigate, for $d=2$, the finite-sample performance
of the same tests by generating $N=1000$ replications of $n=100$ i.i.d. random
vectors from

  1. (i)

the uniform distribution;

  2. (ii)

vMF distributions with location parameter
${\boldsymbol{\theta}}=(0,1)^{\top}$, and concentration parameters
$\kappa=0.05$, $\kappa=0.1$, and $\kappa=0.5$;

  3. (iii)

mixtures of two vMFs, of the form $I[U\leq
0.7]\mathbf{Z}_{1}+I[U>0.7]\mathbf{Z}_{2}$, where $U\sim{\rm U}[0,1]$,
$\mathbf{Z}_{1}\sim\mathcal{M}_{2}({\boldsymbol{\theta}}_{1},\kappa_{1})$, and
$\mathbf{Z}_{2}\sim\mathcal{M}_{2}({\boldsymbol{\theta}}_{2},\kappa_{2})$ are
mutually independent, with
${\boldsymbol{\theta}}_{1}=(-0.3,\sqrt{0.91})^{\top}$ and
${\boldsymbol{\theta}}_{2}=(0.6,0.8)^{\top}$, and concentration parameters
$\kappa_{1}=\kappa_{2}=0.1$, $\kappa_{1}=\kappa_{2}=0.25$, and
$\kappa_{1}=\kappa_{2}=0.5$;

  4. (iv)

the sine-skew distribution proposed by [86] and [1], the density of which
takes form $\phi\mapsto
f(\phi-\mu)(1+\lambda\sin(\phi-\mu),\quad\phi\in[-\pi,\pi),$ where $f$ is a
density (over $[-\pi,\pi)$) symmetric about zero, $\mu\in[-\pi,\pi)$ is an
angular location, and $\lambda\in(-1,1)$ a skewness parameter. See Ley and
Verdebout [60, Chapter 2] for the data-generating process. In the simulations,
we chose $f$ to be a vMF density symmetric about $0$ with concentration
parameters $\kappa=0.1$, set $\mu=0$ and $\lambda=0.1$, 0.3, and 0.35.

In Table 2, we report the rejection frequencies (at nominal level
$\alpha=0.05$) of the OT, PCvM, PAD, PRt, Rayleigh, Bingham, Ajne, Giné, and
Bakshaev tests of uniformity over $\mathcal{S}^{1}$. The results are
consistent with those for $d=3$ in Table 1. Specifically, all tests have
rejection frequencies close to the nominal level under the null hypothesis of
uniformity. Under the vMF alternatives, the Bingham and Giné tests have less
power than all other ones. Our test outperforms all its competitors whenever
skewness or multi-modality are present.

Table 2: Rejection frequencies of the OT, PCvM, PAD, PRt, Rayleigh, Bingham, Ajne, Giné, and Bakshaev tests of uniformity over $\mathcal{S}^{1}$ under the uniform, vMF, mixtures of two vMF distributions and sine-skew distributions (simulation settings: $N=1000$ replications, sample size $n=100$) | OT | PCvM | PAD | PRt | Rayleigh | Bingham | Ajne | Giné | Bakshaev
---|---|---|---|---|---|---|---|---|---
Uniform | 0.058 | 0.053 | 0.055 | 0.052 | 0.054 | 0.039 | 0.054 | 0.042 | 0.054
vMF ($\kappa=0.05$) | 0.080 | 0.069 | 0.068 | 0.068 | 0.067 | 0.054 | 0.067 | 0.051 | 0.070
vMF ($\kappa=0.1$) | 0.105 | 0.098 | 0.097 | 0.096 | 0.096 | 0.045 | 0.098 | 0.048 | 0.097
vMF ($\kappa=0.5$) | 0.877 | 0.873 | 0.867 | 0.876 | 0.875 | 0.064 | 0.873 | 0.069 | 0.875
| Mixture of two vMFs
---
($\kappa_{1}=\kappa_{2}=0.1$)
0.092 | 0.080 | 0.079 | 0.079 | 0.083 | 0.045 | 0.084 | 0.043 | 0.078
| Mixture of two vMFs
---
($\kappa_{1}=\kappa_{2}=0.25$)
0.288 | 0.256 | 0.249 | 0.256 | 0.264 | 0.055 | 0.262 | 0.053 | 0.260
| Mixture of two vMFs
---
($\kappa_{1}=\kappa_{2}=0.5$)
0.847 | 0.829 | 0.822 | 0.831 | 0.831 | 0.057 | 0.830 | 0.058 | 0.827
Sine-skew ($\lambda=0.1$) | 0.093 | 0.080 | 0.079 | 0.084 | 0.081 | 0.041 | 0.081 | 0.041 | 0.080
Sine-skew ($\lambda=0.3$) | 0.498 | 0.465 | 0.457 | 0.465 | 0.473 | 0.035 | 0.470 | 0.038 | 0.471
Sine-skew ($\lambda=0.35$) | 0.637 | 0.618 | 0.610 | 0.624 | 0.627 | 0.046 | 0.627 | 0.049 | 0.622

##  6 Directional MANOVA

Let
$\mathbf{X}_{i1}^{(n_{i})},\ldots,\mathbf{X}_{in_{i}}^{(n_{i})},i=1,\ldots,m,$
denote $m(\geq 2)$ independent samples on the unit hypersphere
$\mathcal{S}^{d-1}$. For each $i$, we assume that
$\mathbf{X}_{i1}^{(n_{i})},\ldots,\mathbf{X}_{in_{i}}^{(n_{i})}$ are i.i.d.
with common absolutely continuous distribution ${\rm P}_{i}$, density $f_{i}$,
and directional distribution function $\mathbf{F}_{i}$, $i=1,\ldots,m$. The
null hypothesis of interest is the hypothesis
$\mathcal{H}_{0}:\mathbf{F}_{1}=\ldots=\mathbf{F}_{m}=:\mathbf{F}$
(unspecified $\mathbf{F}$) of no treatment effect.

Let $\mathbf{F}^{(n)}$ denote the empirical directional distribution function
computed from the pooled sample

$\\{\mathbf{Y}^{(n)}_{1},\ldots,\mathbf{Y}^{(n)}_{n}\\}\coloneqq\\{\mathbf{X}_{11}^{(n_{1})},\ldots,\mathbf{X}_{1n_{1}}^{(n_{1})},\ldots,\mathbf{X}_{m1}^{(n_{m})},\ldots,\mathbf{X}_{mn_{m}}^{(n_{m})}\\}.$

The tests we propose here are based on a $(md_{\bf J})$-dimensional statistic
${\mbox{\boldmath$\Delta$}}_{\mathbf{J}}^{(n)}\coloneqq(({\mbox{\boldmath$\Delta$}}_{1;\mathbf{J}}^{(n_{1})})^{\top},\ldots,({\mbox{\boldmath$\Delta$}}_{m;\mathbf{J}}^{(n_{m})})^{\top})^{\top}$
where, for some score function
$\mathbf{J}:\mathcal{S}^{d-1}\to\mathbb{R}^{d_{\bf J}}$,

${\mbox{\boldmath$\Delta$}}_{i;\mathbf{J}}^{(n_{i})}\coloneqq
n_{i}^{-1/2}\sum_{j=1}^{n_{i}}\mathbf{J}(\mathbf{F}^{(n)}(\mathbf{X}_{ij}^{(n_{i})}))-n^{-1/2}\sum_{\ell=1}^{n}\frac{n_{i}^{1/2}}{n^{1/2}}\mathbf{J}(\mathbf{F}^{(n)}(\mathbf{Y}^{(n)}_{\ell}))=n^{-1/2}\sum_{\ell=1}^{n}(a_{i\ell}^{(n)}-\bar{a}_{i}^{(n)})\mathbf{J}(\mathbf{F}^{(n)}(\mathbf{Y}^{(n)}_{\ell})),\
\ i=1,\ldots,m$

with
$a_{i\ell}^{(n)}\\!\coloneqq(n^{1/2}/n_{i}^{1/2})I[{n_{1}+\ldots+n_{i-1}+1\leq\ell\leq
n_{1}+\ldots+n_{i-1}+n_{i}}]$ and
$\bar{a}_{i}^{(n)}\\!=n^{-1}\sum_{\ell=1}^{n}a_{i\ell}^{(n)}=n_{i}^{1/2}/n^{1/2}$
(as usual, $I[{.}]$ denotes the indicator function). Below, we assume that
$r_{i}^{(n)}\coloneqq n_{i}/n$ converges to a constant $r_{i}$ as
$n\rightarrow\infty$. Letting $\mathbf{D}_{\mathbf{J}}\coloneqq{\rm
Var}(\mathbf{J}(\mathbf{U}))$ where $\mathbf{U}$ is uniform over
$\mathcal{S}^{d-1}\\!$, our tests reject ${\cal H}_{0}$ for large values of
${Q}_{\mathbf{J}}^{(n)}\coloneqq({\mbox{\boldmath$\Delta$}}_{\mathbf{J}}^{(n)})^{\top}({\bf
I}_{m}\otimes{\bf
D}_{\mathbf{J}}^{-}){\mbox{\boldmath$\Delta$}}_{\mathbf{J}}^{(n)}$, where
${\bf D}_{\mathbf{J}}^{-}$ denotes the Moore-Penrose inverse of ${\bf
D}_{\mathbf{J}}$. The next result shows that ${Q}_{\mathbf{J}}^{(n)}$ is
asymptotically chi-square under ${\mathcal{H}}_{0}$ and some mild regularity
assumptions on the score function ${\bf J}$.

######  Proposition 6.1.

Assume that (i) $\mathbf{J}$ is continuous over $\mathcal{S}^{d-1}$ and (ii)
$\mathbf{J}$ is square-integrable, that is, $\int_{\mathcal{S}^{d-1}}\|{\bf
J}(\mathbf{u})\|^{2}{\rm dP}^{\mathbf{U}}(\mathbf{u})<\infty$, and (iii) for
any sequence
$\mathfrak{G}^{(n)}\coloneqq\\{\mathfrak{G}_{1}^{(n)},\ldots,\mathfrak{G}_{n}^{(n)}\\}$
of $n$-tuples in $\mathcal{S}^{d-1}$ such that the uniform discrete
distribution over $\mathfrak{G}^{(n)}$ converges weakly to ${\rm
P}^{\mathbf{U}}$ as $n\rightarrow\infty$,
$\underset{n\rightarrow\infty}{\lim}n^{-1}\sum_{{\ell}=1}^{n}\|{\bf
J}(\mathfrak{G}_{\ell}^{(n)})\|^{2}=\int_{\mathcal{S}^{d-1}}\|{\bf
J}(\mathbf{u})\|^{2}{\rm dP}^{\mathbf{U}}(\mathbf{u}).$ Then, under ${\cal
H}_{0}$ as $n\rightarrow\infty$, ${Q}_{\mathbf{J}}^{(n)}$ is asymptotically
chi-square with $(m-1)d^{*}$ degrees of freedom, where $d^{*}$ is the rank of
the $d_{\bf J}\times d_{\bf J}$ matrix ${\bf D}_{\mathbf{J}}$.

The MANOVA tests $\phi_{\mathbf{J}}^{(n)}$ we propose reject the hypothesis of
no treatment effect at asymptotic level $\alpha$ whenever

${Q}_{\mathbf{J}}^{(n)}>\chi^{2}_{(m-1)d^{*};1-\alpha},$

where $\chi^{2}_{p;\tau}$ stands for the quantile of order $\tau$ of a chi-
square distribution with $p$ degrees of freedom. Being based on the
$\mathbf{F}^{(n)}$-measurable test statistic ${Q}_{\mathbf{J}}^{(n)}$, these
tests are fully distribution-free. Below, we investigate their asymptotic
properties under local alternatives. To do so, let us consider a parametric
framework where the underlying distributions are indexed by some finite-
dimensional parameter
$\bm{\omega}\in\mbox{\boldmath$\Omega$}\subset\mathbb{R}^{k}$. More precisely,
denote by ${\rm P}^{(n)}_{\underline{{\scriptstyle{\bm{\omega}}}}}$,where
$\underline{\bm{\omega}}\coloneqq(\bm{\omega}_{1},\ldots,\bm{\omega}_{m})\in\mbox{\boldmath$\Omega$}^{m}\subset\mathbb{R}^{mk}$,
the joint distribution of the pooled sample
$(\mathbf{Y}^{(n)}_{1},\ldots,\mathbf{Y}^{(n)}_{n})$ when the $i$th sample
$(\mathbf{X}_{i1}^{(n_{i})},\ldots,\mathbf{X}_{in_{i}}^{(n_{i})})$ has
distribution ${\rm P}^{(n_{i})}_{{\scriptstyle{\bm{\omega}}}_{i}}$,
$i=1,\ldots,m$. Clearly, under the null hypothesis, ${\rm
P}^{(n)}_{\underline{{\scriptstyle{\bm{\omega}}}}}$ is of the form ${\rm
P}^{(n)}_{({\scriptstyle{\bm{\omega}}}_{0},\ldots,{\scriptstyle{\bm{\omega}}}_{0})}$
for some $\bm{\omega}_{0}\in\mbox{\boldmath$\Omega$}$: write
$\underline{\bm{\omega}}_{0}$ for $(\bm{\omega}_{0},\ldots,\bm{\omega}_{0})$.

In order to study the local power of our tests $\phi_{\mathbf{J}}^{(n)}$,
consider local alternatives of the form
$\underline{\bm{\omega}}_{0}+n^{-1/2}\mbox{\boldmath$\nu$}^{(n)}\mbox{\boldmath$\tau$}^{(n)}\in~{}\\!\mbox{\boldmath$\Omega$}^{m}$,
where $\mbox{\boldmath$\nu$}^{(n)}\coloneqq{\rm diag}((r_{1}^{(n)})^{-1/2}{\bf
I}_{k},\ldots,(r_{m}^{(n)})^{-1/2}{\bf I}_{k})$ and
$\mbox{\boldmath$\tau$}^{(n)}\coloneqq((\mbox{\boldmath$\tau$}_{1}^{(n_{1})})^{\top},\ldots,(\mbox{\boldmath$\tau$}_{m}^{(n_{m})})^{\top})^{\top}$
with $\mbox{\boldmath$\tau$}_{i}^{(n_{i})}$ a bounded sequence of
$\mathbb{R}^{k}$. We will assume that the underlying sequence of experiments
is sufficiently regular in the sense that, under ${\rm
P}^{(n)}_{{\underline{{\scriptstyle{\bm{\omega}}}}}_{0}}$, as
$n\rightarrow\infty$,

$\Lambda^{(n)}\coloneqq\log\frac{{\rm d}{\rm
P}^{(n)}_{{\underline{{\scriptstyle{\bm{\omega}}}}}_{0}+n^{-1/2}{\scriptstyle{\bm{\nu}}}^{(n)}{\scriptstyle{\bm{\tau}}}^{(n)}}}{{\rm
d}{\rm
P}^{(n)}_{{\underline{{\scriptstyle{\bm{\omega}}}}}_{0}}}=(\mbox{\boldmath$\tau$}^{(n)})^{\top}\mbox{\boldmath$\Delta$}^{(n)}_{\underline{{\scriptstyle{\bm{\omega}}}}_{0}}-\frac{1}{2}(\mbox{\boldmath$\tau$}^{(n)})^{\top}{\bm{\mathcal{I}}}_{\underline{{\scriptstyle{\bm{\omega}}}}_{0}}\mbox{\boldmath$\tau$}^{(n)}+o_{\rm
P}(1),$ (6.1)

where
$\mbox{\boldmath$\Delta$}^{(n)}_{{\underline{{\scriptstyle{\bm{\omega}}}}}_{0}}\coloneqq((\mbox{\boldmath$\Delta$}^{(n_{1})}_{1;{\scriptstyle{\bm{\omega}}}_{0}})^{\top},\ldots,(\mbox{\boldmath$\Delta$}^{(n_{m})}_{m;{\scriptstyle{\bm{\omega}}}_{0}})^{\top})^{\top}$,
with $\mbox{\boldmath$\Delta$}^{(n_{i})}_{i;{\scriptstyle{\bm{\omega}}}_{0}}$
of the form

$\mbox{\boldmath$\Delta$}^{(n_{i})}_{i;{\scriptstyle{\bm{\omega}}}_{0}}\coloneqq
n_{i}^{-1/2}\sum_{j=1}^{n_{i}}\mbox{\boldmath$\varphi$}_{{\scriptstyle{\bm{\omega}}}_{0}}(\mathbf{X}_{ij}^{(n_{i})}),\quad
i=1,\ldots,m$ (6.2)

for some square-integrable function
$\mbox{\boldmath$\varphi$}_{{\scriptstyle{\bm{\omega}}}_{0}}$, and
${\bm{\mathcal{I}}}_{{\underline{{\scriptstyle{\bm{\omega}}}}}_{0}}\coloneqq{\rm
diag}({\bm{\mathcal{I}}}_{1;{\scriptstyle{\bm{\omega}}}_{0}},\ldots,{\bm{\mathcal{I}}}_{m;{\scriptstyle{\bm{\omega}}}_{0}})$
are such that, still under ${\rm
P}^{(n)}_{{\underline{{\scriptstyle{\bm{\omega}}}}}_{0}}$,
$\mbox{\boldmath$\Delta$}^{(n)}_{{\underline{{\scriptstyle{\bm{\omega}}}}}_{0}}$
is asymptotically normal with mean ${\bm{0}}$ and covariance matrix
${\bm{\mathcal{I}}}_{{\underline{{\scriptstyle{\bm{\omega}}}}}_{0}}$. That is,
we assume that the underlying model is LAN with central sequence
$\mbox{\boldmath$\Delta$}^{(n)}_{{\underline{{\scriptstyle{\bm{\omega}}}}}_{0}}$
and information matrix
${\bm{\mathcal{I}}}_{\underline{{\scriptstyle{\bm{\omega}}}}_{0}}$. Letting
${\bf K}_{{\bf
J},{\scriptstyle{\bm{\omega}}}_{0}}\coloneqq\int_{\mathcal{S}^{d-1}}{\bf
J}(\mathbf{u})\mbox{\boldmath$\varphi$}_{{\scriptstyle{\bm{\omega}}}_{0}}({\bf
F}^{-1}(\mathbf{u})){\rm dP}^{\mathbf{U}}(\mathbf{u}),$ we have the following
result.

######  Proposition 6.2.

Under ${\rm
P}^{(n)}_{{\underline{{\scriptstyle{\bm{\omega}}}}}_{0}+n^{-1/2}{\scriptstyle{\bm{\nu}}}^{(n)}{\scriptstyle{\bm{\tau}}}^{(n)}}$,
${Q}_{\mathbf{J}}^{(n)}$ is asymptotically non-central chi-square with
$(m-1){d}^{*}$ degrees of freedom and, letting
$\mbox{\boldmath$\tau$}_{j}\coloneqq\lim_{n\rightarrow\infty}\mbox{\boldmath$\tau$}_{j}^{(n)}$,
non-centrality parameter
$\displaystyle{\sum_{j=1}^{m}(1-r_{j})\mbox{\boldmath$\tau$}_{j}^{\top}{\bf
K}_{{\bf J},{\scriptstyle{\bm{\omega}}}_{0}}^{\top}{\bf D}_{{\bf J}}^{-}{\bf
K}_{{\bf
J},{\scriptstyle{\bm{\omega}}}_{0}}\sum_{j=1}^{m}(1-r_{j})\mbox{\boldmath$\tau$}_{j}.}$

Proposition 6.2 and Le Cam’s third lemma readily yield integral expressions
(involving the scores $\bf J$, the directional distribution function ${\bf
F}^{(n)}$, and $\mbox{\boldmath$\varphi$}_{{\scriptstyle{\bm{\omega}}}_{0}}$)
for local and asymptotic powers in this class of parametric models. They also
help select a score function $\bf J$: the test based on ${\bf
J}=\mbox{\boldmath$\varphi$}_{{\scriptstyle{\bm{\omega}}}_{0}}\circ\left({\bf
F}^{(n)}\right)^{-1}$, for instance, achieves optimality in the parametric LAN
model just described. Rather than listing such theoretical expressions for
specific LAN models and specific alternatives (location, concentration,
multimodality, skewness, … ), we conclude this Section by investigating the
corresponding finite-sample performances for various scores:

  1. (i)

the uniform score ${\mathbf{J}}({\mathbf{u}})\coloneqq\mathbf{u}$;

  2. (ii)

the estimated vMF-location score
$\mathbf{J}(\mathbf{F}^{(n)}(\mathbf{Y}^{(n)}_{\ell}))\coloneqq\hat{\kappa}\sqrt{1-\Big{(}G_{\hat{\kappa}}^{-1}\Big{(}1-\frac{R^{(n)}_{\ell}}{n_{R}+1}\Big{)}\Big{)}^{2}}\,\mathbf{S}_{\ell}^{(n)}$,
where $\hat{\kappa}$ is the vMF maximum-likelihood estimator of the
concentration parameter, $G_{\kappa}$ the distribution function of
$\mathbf{Z}^{\top}{\boldsymbol{\theta}}$ with $\mathbf{Z}\sim{\cal
M}_{d}(\mbox{\boldmath$\theta$},\kappa)$, $R^{(n)}_{\ell}\coloneqq
R^{(n)}(\mathbf{Y}_{\ell}^{(n)})=(n_{R}+1)\left[1-F_{*}(({\bf
F}^{(n)}(\mathbf{Y}_{\ell}^{(n)}))^{\top}{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}))\right]$,
and
$\mathbf{S}_{\ell}^{(n)}\coloneqq\mathbf{S}_{{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}}(\mathbf{F}^{(n)}(\mathbf{Y}^{(n)}_{\ell}))$;
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$ is
obtained as in Section 4.2, step 1 of the construction of the grid;

  3. (iii)

the estimated vMF-concentration score
$\mathbf{J}(\mathbf{F}^{(n)}(\mathbf{Y}^{(n)}_{\ell}))\coloneqq
G_{\hat{\kappa}}^{-1}\Big{(}1-\frac{R^{(n)}_{\ell}}{n_{R}+1}\Big{)},$ and

  4. (iv)

the estimated vMF-location-concentration score (a linear combination of the
scores in (ii) and (iii))

$\mathbf{J}(\mathbf{F}^{(n)}(\mathbf{Y}^{(n)}_{\ell}))\coloneqq\hat{\kappa}\Bigg{\\{}G_{\hat{\kappa}}^{-1}\Big{(}1-\frac{R^{(n)}_{\ell}}{n_{R}+1}\Big{)}{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}+\sqrt{1-\Big{(}G_{\hat{\kappa}}^{-1}\Big{(}1-\frac{R^{(n)}_{\ell}}{n_{R}+1}\Big{)}\Big{)}^{2}}\,\mathbf{S}_{\ell}^{(n)}\Bigg{\\}}.\vspace{-0mm}$

In the simulation exercise below, we compare the resulting tests with their
pseudo-vMF (pvMF) counterpart, which rejects the null hypothesis at asymptotic
level $\alpha$ whenever $Q^{(n)}>\chi^{2}_{(m-1)(d-1);1-\alpha}$, where

$\displaystyle Q^{(n)}$
$\displaystyle\coloneqq(d-1)\Bigg{(}\sum_{i=1}^{m}\frac{n_{i}D_{i}}{E_{i}}(\bar{\mathbf{X}}_{i}^{(n_{i})})^{\top}({\bf
I}_{d}-\hat{\mbox{\boldmath$\theta$}}\hat{\mbox{\boldmath$\theta$}}^{\top})\bar{\mathbf{X}}_{i}^{(n_{i})}-\sum_{i,j}^{m}\frac{n_{i}n_{j}}{n}\frac{D_{i}D_{j}}{H}(\bar{\mathbf{X}}_{i}^{(n_{i})})^{\top}({\bf
I}_{d}-\hat{\mbox{\boldmath$\theta$}}\hat{\mbox{\boldmath$\theta$}}^{\top})\bar{\mathbf{X}}_{j}^{(n_{j})}\Bigg{)},$

with $\hat{\mbox{\boldmath$\theta$}}$ the sample Fréchet mean,
$\bar{\mathbf{X}}_{i}^{(n_{i})}\coloneqq
n_{i}^{-1}\sum_{j=1}^{n_{i}}\mathbf{X}_{ij}^{(n_{i})}$, $E_{i}\coloneqq
n_{i}^{-1}\sum_{j=1}^{n_{i}}(\mathbf{X}_{ij}^{(n_{i})})^{\top}\hat{\mbox{\boldmath$\theta$}}$,
$H\coloneqq\sum_{i=1}^{m}r_{i}^{(n)}D_{i}^{2}B_{i}$, $D_{i}\coloneqq
E_{i}/B_{i}$, and $B_{i}\coloneqq
1-n_{i}^{-1}\sum_{j=1}^{n_{i}}((\mathbf{X}_{ij}^{(n_{i})})^{\top}\hat{\mbox{\boldmath$\theta$}})^{2}$;
see [59] for details.

Figure 4: The rejection rates (at $0.05$ nominal level) of the pseudo-vMF and
rank-based tests with uniform, vMF-location, vMF-concentration, and vMF-
location-concentration scores, for Cases (1) (top-left panel), (2) (top-right
panel), (3) (bottom-left panel), and (4) (bottom-right panel).

In the Monte-Carlo simulation, we set $d=3$, $m=2$, $n_{1}=500$, $n_{2}=600$,
hence $n=1100$. The structured regular grid for the computation of the vMF-
location, vMF-concentration, and vMF-location-concentration scores was based
on $n_{R}=44$, $n_{S}=25$, and $n_{0}=0$. The following data-generating
processes were considered.

  1. Case (1)

(location alternatives): $\mathbf{X}_{1j}^{(n_{1})}\\!\\!\sim\\!{\cal
M}_{d}(\mbox{\boldmath$\theta$},\kappa)$, $j=1,\ldots,n_{1}$ and
$\mathbf{X}_{2\ell}^{(n_{2})}\\!\\!\sim\\!{\cal M}_{d}({\bf
O}_{\xi}\mbox{\boldmath$\theta$},\kappa)$, $\ell=1,\ldots,n_{2}$ where
$\kappa=3$,

$\mbox{\boldmath$\theta$}=\left(\begin{array}[]{c}1\\\ 0\\\
0\end{array}\right),\quad\text{ and}\quad{\bf
O}_{\xi}=\begin{pmatrix}\cos(\pi\xi/15)&-\sin(\pi\xi/15)&0\\\
\sin(\pi\xi/15)&\cos(\pi\xi/15)&0\\\ 0&0&1\end{pmatrix}\quad\text{ with
}\xi=0,\ 0.2,\ 0.4,\ 0.6,\ 0.8;\vspace{-1mm}$ (6.3)

  2. Case (2)

(concentration alternatives): $\mathbf{X}_{1j}^{(n_{1})}\\!\\!\sim\\!{\cal
M}_{d}(\mbox{\boldmath$\theta$},\kappa)$, $j=1,\ldots,n_{1}$ and
$\mathbf{X}_{2\ell}^{(n_{2})}\\!\\!\sim\\!{\cal
M}_{d}(\mbox{\boldmath$\theta$},\xi+\kappa)$, $\ell=1,\ldots,n_{2}$ where
$\kappa=3$ and $\mbox{\boldmath$\theta$}=(1,0,0)^{\top}$, with $\xi=0,\ 0.5,\
1,\ 1.5,\ 2$;

  3. Case (3)

(multimodal alternatives): $\mathbf{X}_{1j}^{(n_{1})}$, $j=1,\ldots,n_{1}$, is
a mixture, with mixing probabilities $\frac{3}{8}$, $\frac{3}{8}$, and
$\frac{1}{4}$, of ${\cal M}_{d}(\mbox{\boldmath$\theta$}_{1},\kappa_{1})$,
${\cal M}_{d}(\mbox{\boldmath$\theta$}_{2},\kappa_{2})$, and ${\cal
M}_{d}(\mbox{\boldmath$\theta$}_{3},\kappa_{3})$, where
$\mbox{\boldmath$\theta$}_{1}=(1,0,0)^{\top}$, $\kappa_{1}=3$,
$\mbox{\boldmath$\theta$}_{2}=(-0.8,0.3,\sqrt{0.27})^{\top}$, $\kappa_{2}=2$,
$\mbox{\boldmath$\theta$}_{3}=(0,-0.7,\sqrt{0.51})^{\top}$, $\kappa_{3}=3$,
and $\mathbf{X}_{2\ell}^{(n_{2})}={\bf O}_{\xi}\mathbf{Z}_{\ell}^{(n_{2})}$,
$\ell=1,\ldots,n_{2}$, where $\mathbf{Z}_{\ell}^{(n_{2})}$ has the same
distribution as $\mathbf{X}_{1j}^{(n_{1})}$ and ${\bf O}_{\xi}$ is as in (6.3)
with $\xi=0,1,\ldots,5$;

  4. Case (4)

(skewed alternatives): $\mathbf{X}_{1j}^{(n_{1})}$, $j=1,\ldots,n_{1}$ and
$\mathbf{X}_{1j\ell}^{(n_{2})}$, $\ell=1,\ldots,n_{1}$ are generated from a
tangent vMF distribution of the form (4.8) with
${\boldsymbol{\theta}}=(0,0,1)^{\top}$,
${\mbox{\boldmath$\mu$\unboldmath}}=(0.7,\sqrt{0.51})^{\top}$, $\kappa=1$,
with $\widetilde{V}\sim\text{\rm Beta}(2,5)$ for $\mathbf{X}_{1j}^{(n_{1})}$
and $\widetilde{V}\sim\text{\rm Beta}(2,5+\xi)$, $\xi=0,0.2,\ldots,1$ for
$\mathbf{X}_{2\ell}^{(n_{2})}$.

For each case, the experiments were repeated $500$ times for the directional
rank-based tests and the pvMF test at $0.05$ nominal level, for the various
values of $\xi$ (where $\xi=0$ yields ${\mathcal{H}}_{0}$). Figure 4 shows a
plot of the rejection frequencies of these tests against $\xi$, for Cases
(1)–(4). Under the null ($\xi=0$), all directional rank-based tests yield
rejection rates close to the nominal level while the pvMF test exhibits a
severe over-rejection frequency of $0.122$ under Case (3), revealing a poor
behavior under multi-modality. In terms of power, for the vMF location problem
(Case (1)), the vMF-location score, uniform score and vMF-location-
concentration score rank-based tests all have rejection rates very close to
the pvMF test which is the optimal test in this case. For the vMF-
concentration parameter problem (Case (2)), the pvMF test has no power and is
outperformed by all the directional rank-based tests, the power of which
increases with $\xi$. For the mixture of three vMF distributions (Case (3)),
the rejection frequencies of the pvMF test under the alternative are
meaningless in view of its size problem; the uniform, vMF-location, and vMF-
location-concentration score rank-based tests have greater power than the vMF-
concentration score one. When skewness is present (Case (4)), all the rank-
based tests outperform the classical pvMF test. In general, the superiority of
the vMF-location-concentration score and uniform score rank-based tests is
clear, considering their consistency, sizeable power against multi-modal
alternatives, skewness and vMF distributions (both for location and
concentration treatment effects). Due to its simplicity, the uniform score is
a good choice in terms of computational efficiency.

##  7 Applications: two real-data examples

###  7.1 MANOVA for Sunspots data

Sunspots are regions on the Sun’s photosphere that are darker and cooler than
the surrounding areas. They are caused by concentrations of magnetic flux that
inhibit convection and are temporary phenomena that experience continuous
changes, lasting for hours to days. Their population first increases rapidly
and then decreases slowly over a period of approximately 11 years, which is
referred to as solar cycle. Early in a solar cycle, sunspots appear at higher
latitudes and then move towards the equator as the cycle approaches maximum, a
phenomenon known as Spörer’s law. Sunspots are widely used to study and
measure solar activity, whose effects may affect earth’s long-term climate
[38].

The data we analyze is based on the Debrecen Photoheliographic Data (DPD)
sunspot catalogue, which contains locations of sunspots since 1974 and is a
continuation of the Greenwich Photoheliographic Results (GPR) catalogue
(spanning 1874-1976). The data is available from the R package rotasym.

Figure 5: Sunspots data: plot of the sunspots of the 22nd solar cycle
($n_{1}=4551$ red points) and the 23rd ($n_{2}=5373$ green points) solar
cycle.

The sunspots of the 22nd (September 1986 to July 1996; $n_{1}=4551$ points, in
red) and 23rd (August 1996 to November 2008; $n_{2}=5373$ points, in green)
solar cycles are shown in in Figure 5. Visual inspection hardly help decide
whether these two samples are from the same distribution or not. According to
[33], various tests suggest rotational symmetry around the north pole for the
23rd cycle while rejecting that hypothesis ($p$-values smaller than $0.02$)
for the 22nd cycle.

We performed the pvMF and the directional rank-based tests based on the
uniform, vMF-location and vMF-concentration scores described in Section 6 for
the null hypothesis of equal distributions in the two samples. For the vMF-
location and vMF-concentration scores, we factorize $n=9924$ into $n_{S}=121$,
$n_{R}=82$, and $n_{0}=2$. The $p$-values of these tests are shown in Table 3.
The pvMF test, with $p$-value 0.14, does not reject the null of the equality
of distributions, even at significance level $0.10$. The rank-based test with
vMF-concentration scores rejects the same hypothesis at level $0.10$, but not
at level $0.05$. The rank-based tests with uniform, vMF-location and vMF-
location-concentration scores, with $p$-values at most 0.036, all suggest that
the two samples have different distributions. This result is in line with the
Monte Carlo experiments (see, e.g., Case (4)), where the directional rank-
based tests significantly outperform the pvMF test when the underlying
distribution are not rotationally symmetric.

scores | pvMF | Unif | vMF-location | vMF-concentration | vMF-location-concentration
---|---|---|---|---|---
$p$-values | 0.140 | 0.036 | 0.005 | 0.052 | 0.005
Table 3: Sunspots data: $p$-values of the pvMF and directional rank-based
tests of the equality of distributions of the sunspots of the 22nd and 23rd
solar cycles.

###  7.2 Quantile contours for Protein data

Proteins are polypeptide chains built up by amino acids, each of them having a
central carbon atom, denoted Cα. Motivated by the key role of Cα atoms in the
protein’s backbone, [56] proposed a representation of the backbone that
encodes, from the position of a Cα atom, the location of the next Cα using the
pseudo-bond joining them (the term pseudo emphasizes that the atoms are not
linked by a single chemical bond but rather by several of them). Since the
length of pseudo-bonds can be considered constant, the chain of Cα atoms can
be represented as a sequence of vectors on the sphere with the parametrization
$\left(\cos(\theta),\sin(\theta)\cos(\tau),\sin(\theta)\sin(\tau)\right)^{\top}$,
with $\theta\in[0,\pi)$and $\tau\in[-\pi,\pi)$, where the origin is set as the
previous Cα atom. This codification is employed in the hidden Markov model of
[45]. Directional models for evolution of proteins structures have been
studied, in particular in [37]. Due to clash-avoiding atom constraints,
$\theta$ usually lies in $[80^{\circ},150^{\circ}]$ whereas $\tau$ can take
any values in $[-180^{\circ},180^{\circ})$. In Figure 6, we provide the
quantile contours associated with a sample of $n=2001$ (factorized into
$n_{R}=40$, $n_{S}=50$, and $n_{0}=1$) pseudo-bond directions of the Cα atom.
These $n=2001$ directions have been extracted from the top 500 dataset in [89]
by taking $n=2001$ points randomly within the points for which both $\theta$
and $\tau$ are available (no missing data). It is evident from the left panel
of Figure 6 that the overall distribution of the spherical vectors is highly
non-symmetric. Our quantile contours adapt impressively well to this
asymmetry.

Figure 6: Protein data: the empirical quantile contours (probability contents
12.2%, 24.4%, 36.6% respectively) computed from the protein data along with
the corresponding empirical meridians ($n_{R}=40$, $n_{S}=50$ and $n_{0}=1$).

##  8 Conclusions

In the present paper, we propose various nonparametric tools for directional
data based on optimal transport maps from the underlying distribution to the
uniform distribution over the hypersphere. More precisely, (i) we propose
concepts of directional distribution and quantile functions, distribution-free
directional signs and ranks, (ii) based empirical distribution functions, we
construct a distribution-free and universally consistent Cramér-von Mises-type
test of uniformity on $\mathcal{S}^{d-1}$, and (iii) based on directional
ranks and signs, we develop a class of fully distribution-free MANOVA
procedures. Contrary to the pseudo-von Mises methods proposed in the
literature—the (asymptotic) validity of which is restricted to rotationally
symmetric distributions—our procedures are unrestrictedly finite-sample valid,
while outperforming their natural competitors (particularly so in the absence
of rotational symmetry.)

## Acknowledgements

The authors express their warm thanks to Alberto González-Sanz for insightful
comments and suggestions on the convergence of empirical transport maps. They
also warmly thank the reviewers for their their careful reading of the
manuscript and for their valuable comments that improved the manuscript.

## Funding

Marc Hallin acknowledges the support of the COST (European Cooperation in
Science and Technology) Action HiTEc CA21163 and the Czech Science Foundation
grant GAČR22036365. Hang Liu’s research is supported by the USTC grants
WK2040000055 and YD2040002016. Thomas Verdebout’s research is supported by an
ARC grant of the ULB and Communauté Française de Belgique, a Projet de
Recherche (PDR) grant of the FNRS, and the Fonds Thelam of the Fondation Roi
Baudouin.

## References

  * Abe et al. [2011] Abe, T., Pewsey, A. et al. (2011) Sine-skewed circular distributions.  Statistical Papers, 52, 683–707.
  * Ajne [1968] Ajne, B. (1968) A simple test for uniformity of a circular distribution.  Biometrika, 55, 343–354.
  * Ambrosio and Pratelli [2003] Ambrosio, L. and Pratelli, A. (2003) Existence and Stability Results in the $L_{1}$ Theory of Optimal Transportation, vol. 1813 of CIME Course, Lecture Notes in Mathematics.  Springer-Verlag Berlin Heidelberg.
  * Ameijeiras-Alonso et al. [2018] Ameijeiras-Alonso, J., Crujeiras, R. M. and Rodríguez-Casal, A. (2018) Directional statistics for wildfires.  In Applied Directional Statistics (eds. C. Ley and T. Verdebout), Chapman & Hall/CRC Interdisciplinary Statistics Series, 187–210. Boca Raton: CRC Press.
  * Bakshaev [2010] Bakshaev, A. (2010) $N$-distance tests of uniformity on the hypersphere.  Nonlinear Anal. Model. Control, 15, 15–28.
  * del Barrio and González-Sanz [2023] del Barrio, E. and González-Sanz, A. (2023) Regularity of center-outward distribution functions in non-convex domains.  arXiv2303.16862.
  * del Barrio et al. [2020] del Barrio, E., González-Sanz, A. and Hallin, M. (2020) A note on the regularity of optimal-transport-based center-outward distribution and quantile functions.  Journal of Multivariate Analysis, 180.
  * del Barrio et al. [2022] — (2022) Nonparametric multiple-output center-outward quantile regression.  arXiv:2204.11756.
  * del Barrio et al. [2021] del Barrio, E., González-Sanz, A. and Loubes, J.-M. (2021) Central limit theorems for general transportation costs.  arXiv:2102.06379.
  * Basu [1959] Basu, D. (1959) The family of ancillary statistics.  Sankhyā, 21, 247–256.
  * Beiglböck [2015] Beiglböck, M. (2015) Cyclical monotonicity and the ergodic theorem.  Ergodic Theory Dynam. Systems, 35, 710–713.
  * Bernoulli [1735] Bernoulli, D. (1735) Quelle est la cause physique de l’inclinaison des plans des orbites des planètes par rapport au plan de l’équateur de la révolution du soleil autour de son axe ; et d’où vient que les inclinaisons de ces orbites sont différentes en elles.  In Recueil des pièces qui ont remporté le prix de l’Académie Royale des Sciences (ed. A. R. des Sciences), vol. 3, 93–122. Paris: Académie Royale des Sciences.
  * Bertsekas [1988] Bertsekas, D. P. (1988) The auction algorithm: A distributed relaxation method for the assignment problem.  Annals of operations research, 14, 105–123.
  * Bertsekas [1992] — (1992) Auction algorithms for network flow problems: A tutorial introduction.  Computational optimization and applications, 1, 7–66.
  * Bhattacharya and Bhattacharya [2008] Bhattacharya, A. and Bhattacharya, R. (2008) Statistics on Riemannian manifolds: asymptotic distributions and curvature.  Proceedings of the American Mathematical Society, 136, 2959–2967.
  * Bhattacharya and Bhattacharya [2012] — (2012) Nonparametric Inference on Manifolds, with Applications to Shape Spaces.  Cambridge University Press.
  * Bingham [1974] Bingham, C. (1974) An antipodally symmetric distribution on the sphere.  Ann. Statist., 2, 1201–1225.
  * Bonneel et al. [2011] Bonneel, N., Van De Panne, M., Paris, S. and Heidrich, W. (2011) Displacement interpolation using lagrangian mass transport.  In Proceedings of the 2011 SIGGRAPH Asia conference, 1–12.
  * Cai et al. [2013] Cai, T., Fan, J. and Jiang, T. (2013) Distributions of angles in random packing on spheres.  J. Mach. Learn. Res., 14, 1837–1864.
  * Cai and Jiang [2012] Cai, T. and Jiang, T. (2012) Phase transition in limiting distributions of coherence of high-dimensional random matrices.  J. Multivar. Anal., 107, 24–39.
  * Chernozhukov et al. [2017] Chernozhukov, V., Galichon, A., Hallin, M. and Henry, M. (2017) Monge-Kantorovich depth, quantiles, ranks and signs.  Ann. Statist., 45, 223–256.
  * Chikuse [2003] Chikuse, Y. (2003) Statistics on Special Manifolds, vol. 174 of Lecture Notes in Statistics.  Heidelberg: Springer.
  * Christie [2015] Christie, D. (2015) Efficient von Mises-Fisher concentration parameter estimation using Taylor series.  J. Stat. Comput. Simul., 85, 3259–3265.
  * Cuesta-Albertos et al. [2009] Cuesta-Albertos, J. A., Cuevas, A. and Fraiman, R. (2009) On projection-based tests for directional and compositional data.  Stat. Comput., 19, 367–380.
  * Cutting et al. [2017] Cutting, C., Paindaveine, D. and Verdebout, T. (2017) Testing uniformity on high-dimensional spheres against monotone rotationally symmetric alternatives.  Ann. Statist., 45, 1024–1058.
  * Deb et al. [2021] Deb, N., Bhattacharya, B. B. and Sen, B. (2021) Efficiency lower bounds for distribution-free Hotelling-type two-sample tests based on optimal transport.  arXiv:2104.01986v2.
  * Deb and Sen [2022] Deb, N. and Sen, B. (2022) Multivariate rank-based distribution-free nonparametric testing using measure transportation.  J. Amer. Statist. Assoc., to appear.
  * Dortet-Bernadet and Wicker [2008] Dortet-Bernadet, J.-L. and Wicker, N. (2008) Model-based clustering on the unit sphere with an illustration using gene expression profiles.  Biostatistics, 9, 66–80.
  * Dryden [2005] Dryden, I. L. (2005) Statistical analysis on high-dimensional spheres and shape spaces.  Ann. Stat., 33, 1643–1665.
  * Figueiredo [2006] Figueiredo, A. (2006) Two-way analysis of variance for data from a concentrated bipolar watson distribution.  Journal of Applied Statistics, 33, 575–581.
  * Gangbo and McCann [1996] Gangbo, W. and McCann, R. J. (1996) The geometry of optimal transportation.  Acta Math., 177, 113 – 161.
  * García-Portugués et al. [2022] García-Portugués, E., Navarro-Esteban, P. and Cuesta-Albertos, J. A. (2022) On a projection-based class of uniformity tests on the hypersphere.  Bernoulli, to appear.
  * García-Portugués et al. [2020] García-Portugués, E., Paindaveine, D. and Verdebout, T. (2020) On optimal tests for rotational symmetry against new classes of hyperspherical distributions.  J. Amer. Statist. Assoc., 115, 1873–1887.
  * García-Portugués and Verdebout [2018] García-Portugués, E. and Verdebout, T. (2018) A review of uniformity tests on the hypersphere.  arXiv:1804.00286.
  * Ghosal and Sen [2022] Ghosal, P. and Sen, B. (2022) Multivariate ranks and quantiles using optimal transport: Consistency, rates and nonparametric testing.  Ann. Statist., 50, 1012 – 1037.  URL: https://doi.org/10.1214/21-AOS2136.
  * Giné [1975] Giné, E. (1975) Invariant tests for uniformity on compact Riemannian manifolds based on Sobolev norms.  Ann. Statist., 3, 1243–1266.
  * Golden et al. [2017] Golden, M., García-Portugués, E., Sørensen, M., Mardia, K. V., Hamelryck, T. and Hein, J. (2017) A generative angular model of protein structure evolution.  Molecular biology and evolution, 34, 2085–2100.
  * Haigh [2007] Haigh, J. D. (2007) The Sun and the Earth’s climate.  Living Rev. Sol. Phys., 4, 2.
  * Hallin [2022] Hallin, M. (2022) Measure transportation and statistical decision theory.  Annual Review of Statistics and Its Application, 9, 401–424.  URL: https://doi.org/10.1146/annurev-statistics-040220-105948.
  * Hallin et al. [2021a] Hallin, M., del Barrio, E., Cuesta-Albertos, J. A. and Matrán, C. (2021a) Center-outward distribution functions, quantiles, ranks, and signs in $\mathbb{R}^{d}$.  Ann. Statist., 48, 1139–1165.
  * Hallin et al. [2022a] Hallin, M., Hlubinka, D. and Hudecová, Š. (2022a) Fully distribution-free center-outward rank tests for multiple-output regression and MANOVA.  J. Amer. Statist. Assoc., to appear.
  * Hallin et al. [2022b] Hallin, M., La Vecchia, D. and Liu, H. (2022b) Center-outward R-estimation for semiparametric VARMA models.  J. Amer. Statist. Assoc., 117.
  * Hallin et al. [2022c] — (2022c) Rank-based testing for semiparametric var models: a measure transportation approach.  Bernoulli, 925–938.
  * Hallin et al. [2021b] Hallin, M., Mordant, G. and Segers, J. (2021b) Multivariate goodness-of-fit tests based on Wasserstein distance.  Electron. J. Statist., 15, 1328–1371.  URL: https://doi.org/10.1214/21-EJS1816.
  * Hamelryck et al. [2006] Hamelryck, T., Kent, J. T. and Krogh, A. (2006) Sampling realistic protein conformations using local structural bias.  PLoS Comput. Biol., 2, e131.
  * Jupp and Kume [2020] Jupp, P. E. and Kume, A. (2020) Measures of goodness of fit obtained by almost-canonical transformations on Riemannian manifolds.  J. Multivar. Anal., 176, 104579.
  * Kanika et al. [2015] Kanika, Kumar, S. and SenGupta, A. (2015) A unified approach to decision-theoretic properties of the MLEs for the mean directions of several Langevin distributions.  J. Multivar. Anal., 133, 160–172.
  * Kent [1982] Kent, J. T. (1982) The Fisher-Bingham distribution on the sphere.  J. R. Stat. Soc. Ser. B Methodol., 44, 71–80.
  * Kent et al. [2018] Kent, J. T., Bhattacharjee, S., Hussein, I. I., Faber, W. R. and Jah, M. (2018) Fisher-Bingham-Kent mixture models for angles-only observation processing.  In 2018 Space Flight Mechanics Meeting. Reston: American Institute of Aeronautics and Astronautics.
  * Kim et al. [2016] Kim, P. T., Koo, J.-Y. and Pham Ngoc, T. M. (2016) Supersmooth testing on the sphere over analytic classes.  J. Nonparametr. Stat., 28, 84–115.
  * Kolouri et al. [2017] Kolouri, S., Park, S. R., Thorpe, M., Slepcev, D. and Rohde, G. K. (2017) Optimal mass transport: Signal processing and machine-learning applications.  IIEEE Signal Process. Mag., 34, 43–59.
  * Kulkarni and SenGupta [2022] Kulkarni, H. V. and SenGupta, A. (2022) An efficient test for homogeneity of mean directions on the hyper-sphere.  International Statistical Review, 90, 41–61.
  * Kume et al. [2013] Kume, A., Preston, S. P. and Wood, A. T. A. (2013) Saddlepoint approximations for the normalizing constant of Fisher-Bingham distributions on products of spheres and Stiefel manifolds.  Biometrika, 100, 971–984.
  * Kume and Sei [2018] Kume, A. and Sei, T. (2018) On the exact maximum likelihood inference of Fisher-Bingham distributions using an adjusted holonomic gradient method.  Stat. Comput., 28, 835–847.
  * Lacour and Pham Ngoc [2014] Lacour, C. and Pham Ngoc, T. M. (2014) Goodness-of-fit test for noisy directional data.  Bernoulli, 20, 2131–2168.
  * Levitt [1976] Levitt, M. (1976) A simplified representation of protein conformations for rapid simulation of protein folding.  J. Mol. Biol., 104, 59–107.
  * Ley et al. [2014] Ley, C., Sabbah, C. and Verdebout, T. (2014) A new concept of quantiles for directional data and the angular Mahalanobis depth.  Electron. J. Stat., 8, 795–816.
  * Ley et al. [2013] Ley, C., Swan, Y., Thiam, B. and Verdebout, T. (2013) Optimal R-estimation of a spherical location.  Statist. Sinica, 23, 305–332.
  * Ley et al. [2017] Ley, C., Swan, Y. and Verdebout, T. (2017) Efficient ANOVA for directional data.  Ann. Inst. Statist. Math., 69, 39–62.
  * Ley and Verdebout [2017] Ley, C. and Verdebout, T. (2017) Modern Directional Statistics.  Chapman & Hall/CRC. Boca Raton: CRC Press.
  * Ley and Verdebout [2018] Ley, C. and Verdebout, T. (eds.) (2018) Applied Directional Statistics.  Chapman & Hall/CRC. Boca Raton: CRC Press.
  * Loeper [2009] Loeper, G. (2009) On the regularity of solutions of optimal transportation problems.  Acta Math., 202, 241–283.
  * Loeper [2011] Loeper, G. e. (2011) Regularity of optimal maps on the sphere: the quadratic cost and the reflector antenna.  Arch. Ration. Mech. Anal., 199, 269–289.
  * Mardia and Jupp [1999] Mardia, K. V. and Jupp, P. E. (1999) Directional Statistics.  Wiley Series in Probability and Statistics. Chichester: Wiley.
  * Marinucci and Peccati [2011] Marinucci, D. and Peccati, G. (2011) Random Fields on the Sphere.  London Mathematical Society Lecture Note Series. Cambridge: Cambridge University Press.
  * Marinucci et al. [2008] Marinucci, D., Pietrobon, D., Balbi, A., Baldi, P., Cabella, P., Kerkyacharian, G., Natoli, P., Picard, D. and Vittorio, N. (2008) Spherical needlets for cosmic microwave background data analysis.  Mon. Not. R. Astron. Soc., 383, 539–545.
  * McCann [1995] McCann, R. J. (1995) Existence and uniqueness of monotone measure-preserving maps.  Duke Mathematical Journal, 80, 309–323.
  * McCann [2001] — (2001) Polar factorization of maps on Riemannian manifolds.  Geom. Funct. Anal., 11, 589–608.
  * Paindaveine and Verdebout [2017] Paindaveine, D. and Verdebout, T. (2017) Inference on the mode of weak directional signals: a Le Cam perspective on hypothesis testing near singularities.  Ann. Statist., 45, 800–832.
  * Paindaveine and Verdebout [2020a] — (2020a) Detecting the direction of a signal on high-dimensional spheres: Non-null and le cam optimality results.  Probab. Theory Relat. Fields, 176, 1165–1216.
  * Paindaveine and Verdebout [2020b] — (2020b) Inference for spherical location under high concentration.  Ann. Statist., 48, 2982–2998.
  * Garcí a Portugués et al. [2014] Garcí a Portugués, E., Barros, A. M. G., Crujeiras, R. M., González-Manteiga, W. and Pereira, J. (2014) A test for directional-linear independence, with applications to wildfire orientation and size.  Stoch. Environ. Res. Risk Assess., 28, 1261–1275.
  * Pratelli [2008] Pratelli, A. (2008) On the sufficiency of $c$-cyclical monotonicity for optimality of transport plans.  Mathematische Zeitschrift, 258, 677–690.
  * Rao and SenGupta [2001] Rao, J. and SenGupta, A. (2001) Topics in Circular Statistics, vol. 5 of Series on Multivariate Analysis.  Singapore: World Scientific.
  * Rayleigh [1919] Rayleigh, L. (1919) On the problem of random vibrations, and of random flights in one, two, or three dimensions.  Philosophical Magazine, 37, 321–347.
  * Rudin [1976] Rudin, W. (1976) Principles of mathematical analysis, vol. 3.  McGraw-hill New York.
  * Rüschendorf [1996] Rüschendorf, L. (1996) On $c$-optimal random variables.  Statist. Probab. Lett., 27, 267–270.
  * Scealy and Wood [2019] Scealy, J. L. and Wood, A. T. A. (2019) Scaled von Mises–Fisher distributions and regression models for paleomagnetic directional data.  J. Amer. Statist. Assoc., 114, 1547–1560.
  * Schachermayer and Teichmann [2008] Schachermayer, W. and Teichmann, J. (2008) Characterization of optimal transport plans for the Monge-Kantotovich problem.  Proc. Amer. Math. Soc., 137, 519–529.
  * Segers [2023] Segers, J. (2023) Graphical and uniform consistency of estimated optimal transport plans.  arXiv preprint arXiv:2208.02508.
  * SenGupta and Kulkarni [2020] SenGupta, A. and Kulkarni, H. V. (2020) Universal and efficient tests for homogeneity of mean directions of circular populations.  Statistica Sinica, 30, 1995–2021.
  * Shi et al. [2022a] Shi, H., Drton, M. and Han, F. (2022a) Distribution-free consistent independence tests via center-outward ranks and signs.  J. Amer. Statist. Assoc., to appear.  URL: https://doi.org/10.1080/01621459.2020.1782223.
  * Shi et al. [2022b] Shi, H., Hallin, M., Drton, M. and Han, F. (2022b) On universally consistent and fully distribution-free rank tests of vector independence.  Ann. Statist.,, 50, 1933–1959.
  * Shorack [2000] Shorack, G. (2000) Probability for Statisticians.  Springer.
  * Torous et al. [2022] Torous, W., Gunsilius, F. and Rigollet, P. (2022) An optimal transport approach to causal inference.  arXiv:2108.05858.
  * Umbach and Jammalamadaka [2009] Umbach, D. and Jammalamadaka, S. R. (2009) Building asymmetry into circular distributions.  Statistics & probability letters, 79, 659–663.
  * Verdebout [2017] Verdebout, T. (2017) On the efficiency of some rank-based test for the homogeneity of concentrations.  J. Stat. Plann. Inference, 191, 101–109.
  * Villani [2009] Villani, C. (2009) Optimal Transport: Old and New.  Berlin: Springer-Verlag.
  * Word et al. [1999] Word, J. M., Lovell, S. C., LaBean, T. H., Taylor, H. C., Zalis, M. E.$,\ldots,$ and Richardson, D. C. (1999) Visualizing and quantifying molecular goodness-of-fit: small-probe contact dots with explicit hydrogen atoms.  J. Mol. Biol., 285, 1711–1733.

##  Appendix A Supplementary material – Distribution-freeness, ancillarity,
and equivariance

Denote by ${\rm P}_{\mathbf{Z}}^{(n)}$ the joint distribution of the $n$-tuple
$\mathbf{Z}^{(n)}\coloneqq\left(\mathbf{Z}_{1}^{(n)},...,\mathbf{Z}_{n}^{(n)}\right)$
of i.i.d. observations with common distribution ${\rm
P}^{\mathbf{Z}}\in~{}\\!\mathfrak{P}_{d}$ and let
$\mathfrak{P}_{d}^{(n)}\coloneqq\\{{\rm P}_{\mathbf{Z}}^{(n)}|\,{\rm
P}^{\mathbf{Z}}\in\mathfrak{P}_{d}\\}$. Denote by
$\mathbf{Z}^{(n)}_{(\cdot)}=\left(\mathbf{Z}_{(1)}^{(n)},...,\mathbf{Z}_{(n)}^{(n)}\right)$
the order statistic of $\mathbf{Z}^{(n)}$, obtained by ordering the
observations via some arbitrary criterion, e.g., by increasing values of first
components. As in ${\mathbb{R}}^{d}$, the order statistic generates the sigma-
field of permutationally invariant measurable functions of $\mathbf{Z}^{(n)}$,
and it readily follows from the factorization criterion that it is sufficient
for $\mathfrak{P}_{d}^{(n)}$ (see Proposition A.1 below).

It is quite natural, thus, to assume that the estimator
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$ is a
symmetric function of the $\mathbf{Z}^{(n)}$’s, hence is measurable with
respect to $\mathbf{Z}^{(n)}_{(\cdot)}$: were it not, the Rao-Blackwell
theorem, indeed, would recommend to replace it with ${\rm
E}[{{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}}|\
\mathbf{Z}^{(n)}_{(\cdot)}]$ which, by construction is
$\mathbf{Z}^{(n)}_{(\cdot)}$-measurable. The random grid
$\mathfrak{G}^{(n)}({{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}})$
then is also $\mathbf{Z}^{(n)}_{(\cdot)}$-measurable, and so is the $n$-tuple
$\left(\mathbf{F}^{(n)}(\mathbf{Z}_{(1)}^{(n)}),...,\mathbf{F}^{(n)}(\mathbf{Z}_{(n)}^{(n)})\right)$.
The following results extend to the spherical context the essential properties
of univariate and multivariate order statistics, ranks, and signs.

######  Proposition A.1.

Let $\\{\mathbf{Z}_{1}^{(n)},...,\mathbf{Z}_{n}^{(n)}\\}$ be i.i.d. with
distribution ${\rm P}^{\mathbf{Z}}\in\mathfrak{P}_{d}$. Assume that
$\widehat{{\boldsymbol{\theta}}}_{\rm M}$ is
$\mathbf{Z}^{(n)}_{(\cdot)}$-measurable. Denote by $\mathbf{F}^{(n)}$ its
empirical distribution function, computed as described in Section 4.2. Then,
in the notation of Section 4.2,

  1. (i)

the order statistic $\mathbf{Z}_{(\cdot)}^{(n)}$ is sufficient and complete,
hence minimal sufficient, for ${\rm P}^{\mathbf{Z}}\in~{}\\!\mathfrak{P}_{d}$;

  2. (ii)

for all ${\rm P}^{\mathbf{Z}}\in\mathfrak{P}_{d}$, conditional on
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$,
$\mathbf{F}^{(n)}(\mathbf{Z}^{(n)})\coloneqq(\mathbf{F}^{(n)}(\mathbf{Z}_{1}^{(n)}),\ldots,\mathbf{F}^{(n)}(\mathbf{Z}_{n}^{(n)}))$
is uniformly distributed over the $n!/n_{0}!$ permutations (with repetitions,
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$being
counted as $n_{0}$ indistinguishable points in case $n_{0}>~{}\\!1$) of the
grid
$\mathfrak{G}^{(n)}({\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}})$;

  3. (iii)

the vector of ranks ${\bf
R}^{(n)}\coloneqq(R^{(n)}_{1},\ldots,R^{(n)}_{n_{S}})$ is distribution-free
and uniformly distributed over the $n!/n_{0}!(n_{S}!)^{n_{R}}$ permutations
with repetitions of $\\{0\text{ ($n_{0}$ copies), }1\text{ ($n_{S}$ copies),
}\ldots,n_{R}\text{ ($n_{S}$ copies)}\\}$; the rank $R^{(n)}_{i}$,
$i=1,\ldots,n$ takes value 0 with probability $n_{0}/n$ and value $j$ with
probability $n_{S}/n$, $j=1,\ldots,n_{R}$;

  4. (iv)

conditionally on
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$, the
vector of signs $({\bf S}^{(n)}_{1},\ldots,{\bf S}^{(n)}_{n})$ is uniformly
distributed over the $n!/n_{0}!(n_{R}!)^{n_{S}}$ permutations with repetitions
of $\\{{\bm{0}}\text{ ($n_{0}$ copies),
}{\boldsymbol{\Gamma}}_{{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}}{{\bf
s}}_{1}\text{ ($n_{R}$ copies),
}\ldots,{\boldsymbol{\Gamma}}_{{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}}{{\bf
s}}_{n_{S}}\text{ ($n_{R}$ copies)}\\}$;

  5. (v)

conditional on
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$ and
for $n_{0}=0$, the vector of ranks $(R_{1}^{(n)},\ldots,R_{n}^{(n)})$ and the
vector of signs $({\bf S}_{1}^{(n)},\ldots,{\bf S}_{n}^{(n)})$ are mutually
independent for all ${\rm P}^{\mathbf{Z}}\in\mathfrak{P}_{d}$.

Note that part (v) of this proposition only holds for $n_{0}=0$. For
$n_{0}>0$, indeed, $R^{(n)}_{i}=~{}\\!0$ automatically implies
${\mathbf{S}^{(n)}_{i}={\bm{0}}}$. However, since $n_{0}/n\rightarrow 0$ as
$n\to\infty$, this dependence between ranks and signs induced by the $n_{0}$
observations mapped to
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$ is
asymptotically negligible. A tie-breaking procedure similar to that proposed
in [40] easily allows for removing this restriction; details are skipped.

An essential property of traditional (univariate) ranks and the center-outward
ranks and signs introduced in [40] is their maximal ancillarity, which is the
source of their success in inferential problems. It follows from Proposition
A.1 that, due to the dependence on
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$ of
the signs $\mathbf{S}^{(n)}_{i}$, this property cannot hold for the
directional ranks and signs defined in (4.5) and (4.4). However, defining
$\overline{\mathbf{S}}^{(n)}_{i}\coloneqq{\boldsymbol{\Gamma}}^{\top}_{{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}}\mathbf{S}^{(n)}_{i}\quad
i=1,\ldots,n,$ with values in ${\mathfrak{S}}^{(n_{S})}$ (in
$\\{\bm{0}\\}\cup{\mathfrak{S}}^{(n_{S})}$ if $n_{0}>0$), we recover this
highly desirable property: call $\overline{\mathbf{S}}^{(n)}_{i}$ the absolute
sign of $\mathbf{Z}^{(n)}_{i}$. Absolute signs are slightly less informative
than the signs; actually, they characterize the signs up to an orthogonal
transformation of the equatorial hypersphere—i.e., up to the (arbitrary)
choice of an orthonormal coordinate system (the columns of
${\boldsymbol{\Gamma}}_{{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}}$)
for the equatorial hyperplane defined by
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}\\!$.
In particular, all cosines of the form
$(\mathbf{S}^{(n)}_{i})^{\top}\mathbf{S}^{(n)}_{j}$ are preserved in the
transformation from signs to absolute signs:

$(\mathbf{S}^{(n)}_{i})^{\top}\mathbf{S}^{(n)}_{j}=(\overline{\mathbf{S}}^{(n)}_{i})^{\top}\overline{\mathbf{S}}^{(n)}_{j}\quad
i,j=1,\ldots,n.$ (A.1)

######  Proposition A.2.

Under the assumptions of Proposition A.1,

  1. (i)

the vector of absolute signs
$(\overline{\mathbf{S}}^{(n)}_{1},\ldots,\overline{\mathbf{S}}^{(n)}_{n})$ is
distribution-free under $\mathfrak{P}_{d}^{(n)}$ and uniformly distributed
over the $n!/n_{0}!(n_{R}!)^{n_{S}}$ permutations with repetitions of
$\\{{\bm{0}}\text{ ($n_{0}$ copies), }{\mathbf{s}}_{1}\text{ ($n_{R}$
copies),}\ldots,\mathbf{s}_{n_{S}}\text{ ($n_{R}$ copies)}\\}$;

  2. (ii)

for $n_{0}=0$, the vector of ranks $(R_{1}^{(n)},\ldots,R_{n}^{(n)})$ and the
vector
$(\overline{\mathbf{S}}^{(n)}_{1},\ldots,\overline{\mathbf{S}}^{(n)}_{n})$ of
absolute signs are mutually independent for all ${\rm
P}^{\mathbf{Z}}\in\mathfrak{P}_{d}$;

  3. (iii)

the $\sigma$-field generated by the ranks $R^{(n)}_{1},\ldots,R^{(n)}_{n}$ and
the absolute signs
$\overline{\mathbf{S}}^{(n)}_{1},\ldots,\overline{\mathbf{S}}^{(n)}_{n}$ is
maximal ancillary for the family of distributions $\mathfrak{P}_{d}^{(n)}$.

Here again, the restriction of point (ii) to $n_{0}=0$ can be waived via the
tie-breaking procedure described in [40].

The importance of absolute signs is of a theoretical nature and stems from the
maximal ancillarity property (Proposition A.2 (iii)). They seldom are needed
in practice, where either a coordinate system of the parallel/meridian type is
to be constructed—in which case the signs as defined in (4.4) are fine—or
distribution-free tests are to be performed, which typically are based on
quadratic test statistics involving signs via scalar products of the form
(A.1) only, hence do not depend on the type of signs adopted.

Equivariance and invariance properties with respect to orthogonal
transformations play a crucial role in directional statistics. The following
proposition provides such invariance/equivariance properties of the empirical
distribution functions, signs, and ranks.

######  Proposition A.3.

Denote by $\mathfrak{G}^{(n)}$ an $n$-point grid and by ${\bf
z}\coloneqq\left({\bf z}_{1},\ldots,{\bf z}_{n}\right)$ an arbitrary $n$-tuple
of points on ${\mathcal{S}}^{d-1}$, by ${\bf O}$ a $d\times d$ orthogonal
matrix, and by ${\bf F}^{(n)}_{{\bf Oz}}$ the empirical distribution function
constructed from ${\bf Oz}\coloneqq\left({\bf Oz}_{1},\ldots,{\bf
Oz}_{n}\right)$ and the grid ${\bf O}\mathfrak{G}^{(n)}$ as described in
Section 4.2. Assume that
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}$ is
such that
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}({\bf
Oz})={\bf
O}\,{\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}({\bf
z})$. Then,

  1. (i)

${\bf F}^{(n)}_{{\bf Oz}}({\bf Oz}_{i})={\bf O}{\bf F}^{(n)}_{\bf z}({\bf
z}_{i})$, $i=1,\ldots,n$;

  2. (ii)

denoting by ${R}_{\bf Oz}^{(n)}({\bf Oz}_{i})$ and $\overline{\bf S}_{\bf
Oz}^{(n)}({\bf Oz}_{i})$, respectively, the rank and the absolute sign of
${\bf Oz}_{i}$ computed from ${\bf F}^{(n)}_{{\bf Oz}}({\bf Oz}_{i})$,

${R}_{\bf Oz}^{(n)}({\bf Oz}_{i})={R}_{\bf z}^{(n)}({\bf
z}_{i})\quad\text{and}\quad\overline{\bf S}_{\bf Oz}^{(n)}({\bf
Oz}_{i})=\overline{\bf S}_{\bf z}^{(n)}({\bf z}_{i}),\quad i=1,\ldots,n.$

The empirical directional distribution function thus is equivariant under
orthogonal transformations; note that the assumption of orthogonal
equivariance property made on
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}\\!({\bf
z})$ is quite natural, and satisfied by all sensible estimators, among which
the Fréchet mean and the corresponding
${\widehat{\boldsymbol{\theta}}}^{\phantom{n}}{\phantom{\theta}}^{{(n)}}\\!({\bf
z})$.

##  Appendix B Supplementary material – Proofs

Proof of Proposition 2.3. First note that, for the squared Riemannian distance
$c$ used here, the a.e. uniqueness of the optimal transport plan is guaranteed
by [68]. Moreover, it follows from Section 2 in [31] that, since the cost
function $c$ belongs to the class of continuous mappings with domain ${\cal
S}^{d-1}\\!\times{\cal S}^{d-1}$, the support of the (unique) optimal
transport plan $({\bf I}_{\rm d}\times{\bf F})\\#{\rm P}^{\mathbf{Z}}$ is
$c$-cyclically monotone. Conversely, it directly follows from Theorem 1 in
[11] that a transport plan with $c$-cyclically monotone support is optimal.
$\square$

Proof of Proposition 3.2. In view of Corollary 10 of [68],
$\mathbf{F}=\exp_{\mathbf{x}}(\nabla\psi(\mathbf{x}))$ is a bijection from
$\mathcal{S}^{d-1}$ to $\mathcal{S}^{d-1}$. Moreover, $\mathbf{F}$ is
continuous over $\mathcal{S}^{d-1}$ according to Theorem 2.4 of [63]. Since
$\mathcal{S}^{d-1}$ equipped with $d$ is a compact metric space, it follows
(see e.g. Theorem 4.17 of [76]) that $\mathbf{F}^{-1}$ is continuous. The
result follows. $\square$

Proof of Proposition 3.5. We first prove that, if the distribution of
$\mathbf{Z}$ is rotationally symmetric with location
${\boldsymbol{\theta}}_{\rm M}$ and angular function $f$, then

$F_{{f}}^{*}({\bf Z}^{\top}{\boldsymbol{\theta}}_{\rm
M}){\boldsymbol{\theta}}_{\rm M}+\sqrt{1-(F^{*}_{{f}}({\bf
Z}^{\top}{\boldsymbol{\theta}}_{\rm
M}))^{2}}\mathbf{S}_{{\boldsymbol{\theta}}_{\rm M}}(\mathbf{Z})$ (B.1)

is uniform over $\mathcal{S}^{d-1}$. Indeed, from the classical tangent-normal
decomposition (see e.g. Paindaveine and Verdebout, 2017), we have that

$\mathbf{Z}=(\mathbf{Z}^{\top}{\boldsymbol{\theta}}_{\rm
M}){\boldsymbol{\theta}}_{\rm
M}+(1-(\mathbf{Z}^{\top}{\boldsymbol{\theta}}_{\rm
M})^{2})^{1/2}{\mathbf{S}}_{{\boldsymbol{\theta}}_{\rm M}}(\mathbf{Z})$

is such that ${\mathbf{S}}_{{\boldsymbol{\theta}}_{\rm M}}(\mathbf{Z})$ is
uniformly distributed over
$\mathcal{S}^{p-2}_{{\boldsymbol{\theta}}^{\perp}_{\rm M}}$ and is independent
of $\mathbf{Z}^{\top}{\boldsymbol{\theta}}_{\rm M}$, which has distribution
function $F_{f}$. It directly follows from the definition of $F_{{f}}^{*}$
that $F_{{f}}^{*}({\bf Z}^{\top}{\boldsymbol{\theta}}_{\rm M})$ has
distribution function $F_{*}$ so that (B.1) is uniform over
$\mathcal{S}^{d-1}$. Now, to prove the $c$-cyclical monotonicity of the
corresponding transport plan, it is enough show that, for any
$\mathbf{z}_{1},\mathbf{z}_{2}\in~{}\\!\mathcal{S}^{d-1}$,

$c(\mathbf{z}_{1},\mathbf{F}(\mathbf{z}_{1}))+c(\mathbf{z}_{2},\mathbf{F}(\mathbf{z}_{2}))\leq
c(\mathbf{z}_{2},\mathbf{F}(\mathbf{z}_{1}))+c(\mathbf{z}_{1},\mathbf{F}(\mathbf{z}_{2})).$
(B.2)

Clearly, the optimal transport $\mathbf{T}$ must satisfy
$\mathbf{S}_{{\boldsymbol{\theta}}_{\rm
M}}(\mathbf{T}(\mathbf{z}))=\mathbf{S}_{{\boldsymbol{\theta}}_{\rm
M}}(\mathbf{z})$. Define the class

$\mathcal{F}\coloneqq\\{\mathbf{T}\in S({\rm P}_{{\boldsymbol{\theta}}_{\rm
M},f},{\rm P}^{\bf U}),\mathbf{S}_{{\boldsymbol{\theta}}_{\rm
M}}(\mathbf{T}(\mathbf{z}))=\mathbf{S}_{{\boldsymbol{\theta}}_{\rm
M}}(\mathbf{z})\\}$

of transport maps such that $\mathbf{S}_{{\boldsymbol{\theta}}_{\rm
M}}(\mathbf{T}(\mathbf{z}))=\mathbf{S}_{{\boldsymbol{\theta}}_{\rm
M}}(\mathbf{z})$. It readily follows from (3.12) that
$\mathbf{F}\in\mathcal{F}$. Now, it remains to prove that $\mathbf{F}$ is
optimal within $\mathcal{F}$. Using again the equivalence between optimality
and $c$-cyclical monotonicity, we only need to prove that (B.2) holds for
$\mathbf{z}_{1},\mathbf{z}_{2}$ satisfying
$\mathbf{S}_{{\boldsymbol{\theta}}_{\rm
M}}(\mathbf{z}_{1})=\mathbf{S}_{{\boldsymbol{\theta}}_{\rm
M}}(\mathbf{z}_{2})$. The problem is essentially equivalent to showing the
$c$-cyclical monotonicity property of the transport on a unit circle with
${\boldsymbol{\theta}}_{M}$ as the origin. Without loss of generality, let us
assume that $d(\mathbf{z}_{1},{\boldsymbol{\theta}}_{\rm
M})<d(\mathbf{z}_{2},{\boldsymbol{\theta}}_{\rm M})$ and
$d(\mathbf{z}_{2},{\boldsymbol{\theta}}_{\rm
M})<d(\mathbf{F}(\mathbf{z}_{1}),{\boldsymbol{\theta}}_{\rm M})$. Let

$\alpha_{1}\coloneqq d(\mathbf{z}_{1},\mathbf{z}_{2}),\quad\alpha_{2}\coloneqq
d(\mathbf{z}_{2},\mathbf{F}(\mathbf{z}_{1})),\ \ \text{and}\ \
\alpha_{3}\coloneqq d(\mathbf{F}(\mathbf{z}_{1}),\mathbf{F}(\mathbf{z}_{2})).$

Since, in view of (3.12),
$\mathbf{z}_{1}^{\top}{\boldsymbol{\theta}}>\mathbf{z}_{2}^{\top}{\boldsymbol{\theta}}$
implies $(\mathbf{F}(\mathbf{z}_{1}))^{\top}{\boldsymbol{\theta}}_{\rm
M}>(\mathbf{F}(\mathbf{z}_{2}))^{\top}{\boldsymbol{\theta}}_{\rm M}$, we have
that

$d(\mathbf{z}_{1},\mathbf{F}(\mathbf{z}_{2}))=\alpha_{1}+\alpha_{2}+\alpha_{3}.$

The problem thus reduces to proving that

$(\alpha_{1}+\alpha_{2})^{2}+(\alpha_{2}+\alpha_{3})^{2}\leq(\alpha_{1}+\alpha_{2}+\alpha_{3})^{2}+\alpha_{2}^{2},$

which obviously holds true. $\square$

Proof of Proposition 4.1. In the population version, the transport map ${\bf
F}$ induces the transport plan

$\gamma\coloneqq({\bf I}_{\rm d}\times{\bf F})\\#{\rm
P}^{\mathbf{Z}}\in\Gamma({\rm P}^{\mathbf{Z}},{\rm P}^{\mathbf{U}}))$ (B.3)