# Change Point Detection with Optimal Transport and Geometric Discrepancy

Pronko Nikita Konstantinovich

###### Abstract

We present novel retrospective change point detection approach based on
optimal transport and geometric discrepancy. The method does not require any
parametric assumptions about distributions separated by change points. It can
be used both for single and multiple change point detection and estimation,
while the number of change points is either known or unknown. This result is
achieved by construction of a certain sliding window statistic from which
change points can be derived with elementary convex geometry in a specific
Hilbert space. The work is illustrated with computational examples, both
artificially constructed and based on actual data.

## Introduction

Change point problem was firstly mentioned by Stewhart [28] as a problem of
industrial quality control in the year 1928. This problem was solved by
Girshick and Rubin [10] in 1953. The optimal solutions for parametric online
formulation of the problem were provided by Shiryaev and Pollak [29],[25].
Asymptotically optimal posteriori change point detection method was developed
by Borovkov [2]. This method requires knowledge of parametric likelihood
functions, however it can be translated to nonparametric case with empirical
likelihood [15]. The fundamental results in nonparametric study of change-
point problem are due Brodsky and Darkhovsky [3], and Horvath et al. [6].
These methods are generally tailored for one-dimensional data. Nonparametric
change point inference with multidimensional data is still an open problem.
New methods are still being developed usually with applications of ideas from
other scopes to change point data. For example, [21], where divergence from
cluster network analysis is used in order to estimate both the number and the
locations of change points or [20], where multivariate version of Wilcoxon
rank statistic is computed.

This work is focused on nonparametric change point detection in multivariate
data. This problem has large scope of applications including finance
[31],genomics [23] and signal processing [16].

The problem of optimal transportation were introduced by Monge [22] in the
year 1781. The modern statement of the problem is due to Kantorovich as it was
stated in the seminal paper in 1941 [14]. We recommend [33] and [26] as modern
references to the subject

Vector rank and quantile functions were introduced by Serfiling [27] as part
Depth-Outlyingness-Quantile-Rank (DOQR) paradigm. Monge-kantorovich vector
ranks were introduced by Chernozhukov et al. [4].

Study of discrepancy were started by H. Weyl [34] in 1916 with number theory
as intended area of application, more computational aspects of discrepancy
were pioneered by Niederreiter [24] (quasi- and pseudo- random number
generation) and Hlawka [13] (numeric integration). We say that discrepancy is
geometric when we talk about continuous distributions with fixed well-
identified support, this is the opposition to the case of combinatorial
discrepancy which study nonuniformity of finite distributions. The main
contribution to the theory of quadratic generalized geometric discrepancy is
by Hickernel [11]. This is a work of a highly synthetic nature. This Synthetic
nature is coming from the combination of vector ranks and geometric
discrepancy, optimal transport and change point problems. All this concepts
are rarely to never seen together in scientific publications, however in this
paper they are intertwined in a singular computational methodology.

But why these results come to existence only now? The answer to this questions
in our case is rooted into the fact that computing discrepancy used to be
computationally challenging task for non-uniform probability distributions.
However, in the recent publication [4] by Chernozhukov,Carlie, Galichon and
Halin the Monge-Kantorovich vector rank functions were introduced. Vector rank
function can be understood in the context of this work as homeomorphic
extension of classical one-dimensional cumulative distribution functions to
multidimensional spaces. The idea of using low-discrepancy sequences as a tool
for the estimation of the vector rank functions was already present in the
[4]. As the name suggests discrepancy is the natural measure of consistency of
the sequence to be low-discrepancy sequence in the same way as the order
statistic is the natural measure of the sequence to be ordered. Thus, as
classical one-dimensional ranks and order statistic is used in the
construction of the classical Kolmogorov-Smirnov and Crámer-Von Mizes tests of
goodness-of-fit, the vector rank function and discrepancy can be used in order
to construct natural extension of this tests in the multidimensional context.

Thus, the main objective of the work is the exploration of vector rank
discrepancy as methodology for nonparametric high-dimensional statistic. We
select change point problem as the test subject for this inquiry.

The strict measure theoretic approach mentioned above is strengthened in this
work as we speak about probability measures instead of probability
distributions and densities whenever is possible which contradicts current
trends in the statistical science. This can be justified by convenience of
understanding measures as points in certain multidimensional spaces bearing
intrinsic geometric structure. On the other hand, we view all mathematical
constructions in this work as theoretical prototypes of certain concrete
algorithms ad data structures which can be implemented in order to solve
practical tasks. Our approach to the problem is highly inspired by [19] and we
use convex geometric intuition whenever possible.

### acknowledgement

Underlying research was conducted during author’s work for master’s degree in
”Mathematical methods of Optimization and Stochastics” program in NRU Higher
School of Economics, Moscow, Russia. Author wants to thank head of the program
Vladimir Spokoiny, his scientific advisor Andrei Sobolevski for suggesting the
research topic and support, Alexey Naumov and Maxim Panov for the discussion
in IITP RAS; Yuri Yanovich,Alexandra Suvorikova , Elena Chernousova and Alexey
Kroshnin for fruitful discussion at the IUM seminar.

## Notation

This is notation for the use in the sequel.

$\mathbb{R}^{d}$ is a $d$-dimensional euclidean space with $d<\infty$ .
$I=[0,1]$ is a unit interval and $I^{d}$ is a unit hypercube in
$\mathbb{R}^{d}$. In line with Kolmogorov’s axiomatics triple
$(\Omega,\mathcal{F},\mathbb{P})$ is a probability space, and every
observation as a random variable is assumed to be a Borel measurable map from
$\Omega$ to $\mathbb{R}^{d}$. For brevity $\Omega$ denotes the whole
probability measure structure $(\Omega,\mathcal{F},\mathbb{P})$.

Every random variable $\xi:\Omega\to\mathbb{R}^{d}$ defines the probability
measure $P=\xi_{\\#}\mathbb{P}$ that is referred to as the probability law of
$\xi$, where $\xi_{\\#}\mathbb{P}$ denotes a pushforward of the measure
$\mathbb{P}$ by $\xi$ . That is, for every Borel set $B\subset\mathbb{R}^{d}$
it holds that $P(B)=\mathbb{P}(\xi^{-1}(B))$. In this case the notation
$\xi\sim P$ is used. In case $\xi$ is a random process over domain
$\mathcal{T}$ notation $\xi\sim P$ means that $P$ defines the hierarchical
probability distribution system of $\xi$ as a whole. In this case for all
$t,s\in\mathcal{T}$ the measure $P_{t}$ is the distribution of a single
observation $\xi_{t}$ and $P_{t,s}$ is the joint distribution of $\xi_{t}$ and
$\xi_{s}$ i.e. $\xi_{t}\sim P_{t},(\xi_{t},\xi_{s})\sim P_{t,s}$.
$\xi\underset{\mathrm{i.i.d}}{\sim}P$ means that $\xi$ is independently
identically distributed sample (i.i.d) with probability law $P$.

A set of all probability laws with finite first and second moments is denoted
by $\mathcal{P}\left(\mathbb{R}^{d}\right)$. A subset of
$\mathcal{P}\left(\mathbb{R}^{d}\right)$ of probability measures absolutely
continuous with respect to Lebesgue measure $\lambda$ is denoted by
$\mathcal{P}_{\ll\lambda}\left((\right)\mathbb{R})$. $\mathcal{U}^{d}$ is the
uniform distribution over $I^{d}$. It is obvious that
$\mathcal{U}^{d}\in\mathcal{P}_{\ll\lambda}\left(\mathbb{R}^{d}\right)$.In the
sequel $I^{n}$ is assumed to be implicitly equipped with $\mathcal{U}^{n}$ for
each $n\in\mathbb{N}$, which makes it into a probability space, and
$I^{0}=\\{\emptyset\\}$ is equipped with the counting measure.

If $\xi_{n}$ is an sequence of random elements in some metric space with
distance metric $\rho$, when it is said that $\xi$ converges to $x$ in
probability if for each $\varepsilon\in\mathbb{R}_{++}$ it holds
$\lim_{n\to\infty}\mathbb{P}(\rho(\xi_{n},x)>\varepsilon)=0$, which is denoted
as $\xi_{n}\to_{\mathbb{P}}x$. In case $\rho(\xi_{n},x)$ is not a measurable
function for all $n$ the convergence in outer probability may be introduced
and denoted by $\xi\to_{\mathbb{P}\bullet}x$. Having
$A^{\varepsilon}_{n}=\\{\omega\in\Omega:\rho(\xi_{n}(\omega),x)>\varepsilon\\}$
the convergence in outer probability is equivalent to
$\lim_{n\to\infty}\inf\\{\mathbb{P}(B)|B\in\mathcal{F}:A_{n}^{\varepsilon}\subset
B\\}=0$ for all $\varepsilon>0$.

If $\xi$ is a discrete time ergodic process with convergence in probability,

$\frac{1}{n}\sum^{n}_{i=1}f(\xi_{i})\to_{\mathbb{P}}\int_{\mathbb{R}^{d}}f\,\mathrm{d}\bar{P}$
(1)

for every bounded Lipschitz function $f$, we say that $\bar{P}$ is the weak
limit distribution of $\xi$.

## Problem Formulation

In this section a brief review for common variations of change point problem
is provided.

Ordered Set $\mathcal{T}$ is associated with time. It is possible to discuss
change point problem with continuous time [9]. However, this work deals only
with case of discrete time. It is also possible to consider a change point
detection for random fields [3].

Two main types of change point problem are offline and online detection. In
offline detection $\mathcal{T}$ is assumed to be a finite set of size
$|\mathcal{T}|=T$. Without loss of generality, it is assumed that
$\mathcal{T}=\\{1,\ldots,T\\}$. Initial observations $X$ are treated as a time
series indexed by $\mathcal{T}$. Firstly, we discuss the case of a single
change point.

######  Definition 1.

The change point of $X$ is an unknown moment of time $\theta\in\mathcal{T}$
such that $(X_{t})_{t=1}^{\theta}\sim P^{1}$ and $(X_{t})^{T}_{t=\theta+1}\sim
P^{2}$. There are two possibilities. Firstly, both subsamples may be i.i.d
with $(X_{t})^{\theta}_{t=1}\underset{\mathrm{i.i.d}}{\sim}\bar{P}^{1}$ and
$(X_{t})^{T}_{t=\theta+1}\underset{\mathrm{i.i.d}}{\sim}\bar{P}^{2}$. It is
also possible to consider unknowns $\bar{P}^{1}$ and $\bar{P}^{2}$ to be
ergodic weak limit distributions of $(X_{t})^{\theta}_{t=1}$ and
$(X_{t})^{T}_{t=\theta+1}$ respectively.

Then offline single change point detection problem is a hypothesis test for
$H_{0}:\bar{P}_{1}=\bar{P}_{2}$ versus $H_{1}:\bar{P}_{1}\neq\bar{P}_{2}$ for
every possible estimate $\hat{\theta}$ of $\theta$ having $1<\hat{\theta}<T$.
Retrieving valid value of $\hat{\theta}$ is a related problem which will be
referred as a change point estimation.

######  Definition 2.

In case of multiple change points existence of up to $N$ change moments
$\theta_{1}<\ldots<\theta_{n}<\ldots\theta_{N}\in\mathcal{T}$ is assumed. In
simpler formulation of problem value of $N$ is assumed to be known, while in
more complicated one $N$ is also an object of estimation. As the previous
variation the problem is a hypothesis test about distribution of samples
$(X_{t})_{t=\theta_{n-1}+1}^{\theta_{n}}\sim\bar{P}_{n}$ where $\theta_{0}=1$
and $\theta_{N+1}=T$ for simplicity. The problem splits into multiple
hypothesis tests for $H^{n}_{0}:\bar{P}_{n}=\bar{P}_{n+1}$ versus
$H^{n}_{1}:\bar{P}_{n}\neq\bar{P}_{n+1}$ for all possible estimates
$\hat{\theta}$ of $\theta$ in case of known $N$. For unknown $N$ the problem
is structured around testing against
$H_{0}:\bar{P}_{1}\neq\bar{P}_{2}\neq\ldots\neq\bar{P}_{\hat{N}}$ for all
estimates $(\hat{N},\hat{\theta})$ of $(N,\theta)$ where possible values of
$\hat{N}$ are constricted to some meaningful finite set.

For online change point problem $\mathcal{T}\cong\mathbb{Z}$ as an ordered
set, and $X$ can be thought as infinite times series unfolding in real time.
In this case the goal is to identify change point $\theta$ as soon as possible
which means using minimal amount of observations $X_{t}$ with $t\geq\theta$.
Branding every moment of time as a change point will achieve zero delay.
However, this detection procedure is unacceptable as it achieves maximal false
alarm rate. This means that change-point detection algorithm needs to be tuned
in for minimal delay with fixed false alarm rate.

In a more theoretical framework a single change point $\theta_{0}$ can be
considered, which implies that infinite set $\\{X_{t}:t<0\\}$ of observations
with prior distribution can be used for detection. However, in more practical
situation the stream of chanege ponts $\theta:\mathbb{N}\to\mathcal{T}$ is
considered, with only finite sample
$\\{X_{t}:\hat{\theta}_{n-1}<t\leq\theta_{n}\\}$ available for detection of
$\theta_{n}$.

Generally, methods of change point detection are focused on detection of
change in such properties of statistical distributions as mean, variance and
spectral density. The method presented in this work is focused on change in
ergodic weak limit distribution of data. Which means that change in mean and
variance must be detectable, although, change in spectral density is not.
Moreover, changes restricted to such properties of distribution as median or
higher moments also must be traceable. Another restriction of proposed
solution is a demand for all measures $\bar{P}_{n}$ to belong to
$\mathcal{P}_{\ll\lambda}\left(\mathbb{R}^{d}\right)$.

Construction of proposed change point detection starts with offline single
detection problem. We restrict our attention to models with i.i.d distributed
observations along one side of change point. In order to give meaningful
characteristic to point estimates of change point this problem can be
reformulated as certain probabilistic game against Nature. Two sample spaces
$(\Omega^{1},\mathcal{F}^{1},\mathbb{P}^{1})$ and
$(\Omega^{2},\mathcal{F}^{2},\mathbb{P}^{2})$ with corresponding families of
random variables $\xi^{1}$ and $\xi^{2}$ taking values in $\mathbb{R}^{d}$ are
assumed to exist. We assume that both families are intrinsically i.i.d which
means that $\xi^{1}\underset{\mathrm{i.i.d}}{\sim}P_{1}$ and
$\xi^{2}\underset{\mathrm{i.i.d}}{\sim}P_{2}$ and the same holds for
$\xi^{2}$. This means that both sample probability spaces can be decomposed as
$\Omega^{1}=\prod^{\infty}_{t=1}\Omega^{1}_{t}$ and
$\Omega^{2}=\prod^{\infty}_{t=1}\Omega^{2}_{t}$ such that for each
$\omega\in\Omega^{1}$ the random variable $\xi_{t}^{1}(\omega)$ depends only
on $\omega^{t}$. The common sample space is constructed by the rule

$\Omega_{t}=\Omega^{1}_{t}\sqcup\Omega^{2}_{t}=\big{(}\\{1\\}\times\Omega^{1}_{t}\big{)}\cup\big{(}\\{2\\}\times\Omega^{2}_{t}\big{)},\quad\Omega=\prod^{\infty}_{t=1}\Omega_{t}$
(2)

with $\sigma$-algebras generated by a set $\\{\\{1\\}\times
A:A\in\mathcal{F}_{1}\\}\cup\\{\\{2\\}\times B:B\in\mathcal{F}_{2}\\}$ . Note,
that the only possible probability distribution on the set $\\{1,2\\}$ is a
Bernoulli distribution $\mathrm{Bern}(r)$ with a parameter $r\in[0,1]$. So
probability measure on $\Omega$ can be defined by
$\mathbb{P}(S)=r\mathbb{P}_{1}(\iota^{-1}_{1}S)+(1-r)\mathbb{P}_{2}(\iota^{-1}_{2}S)$,
where $\iota_{j}$ is a natural injection defined by $\omega\mapsto(j,\omega)$.
Nature defines value of $r$ and distribution laws $P_{1}$ and $P_{2}$. Random
variables $(b_{t},X_{t})$ are constructed by letting
$(b_{t},X_{t})(j,\omega)=(j,\xi^{j}_{t}(\omega))$. Then, nature generates
sample of $T$ observations of $(b,X)$, sorts it by $b$, then erases values of
$b$ which retains observable sample of $X$. Observer either claims that there
were no change point which leads to a victory in case $P_{1}=P_{2}$ or
equivalently $r$ is too close to either $0$ or $1$. Otherwise, the observer
gives point-estimate $\hat{r}$ of $r$ which can be converted to a classical
change-point estimate $\hat{\theta}=\hat{r}T$.

This construction also can be understood as a Bayesian assumption
$\theta\sim\mathrm{B}(T,r)$, where $\mathrm{B}(T,r)$ stands for binomial
distribution. However, the change point problem is a single element sample
problem, which means that only minimal properties of the parameter can be
inferred from the data. This does not change the problem except for treating
$\hat{\theta}$ as an estimate of the mean.

Obviously,the observer chooses strategy which maximizes chances of victory.
This can be interpreted in terms of loses $L$, which are equal to
$L(0)=\min\\{1-r,r\\}=L(1)$ in case observer claims no change point and
$L(\hat{r})=|\hat{r}-r|$ otherwise. Let $L(\hat{r})=\min{1-\hat{r},\hat{r}}$
in case the observer gives estimate when no real change point exists.

If one preserves $r$ while increasing $T$ the empirical distribution of $X$
will converge to a measure $\mu_{r}=r\bar{P}_{1}+(1-r)\bar{P}_{2}$, the convex
combination of measures. It possible to think of $\mu_{r}$ as contained inside
one-dimensional interval $[\bar{P}_{1},\bar{P}_{2}]$ embedded into the
infinite-dimensional convex space of probability measures. In order to
construct such measures from random variables we will use random mixing map
$I^{s}_{Z}:Y\mapsto b_{s}Y+(1-b_{s})Z$ in the sequel. Here
$b_{s}\sim\mathrm{Bern}(s)$ and $s\in[0,1]$.

In order to propose the solution strategy for the change point problem we
introduce vector rank function $R:\mathbb{R}^{d}\to\mathbb{R}^{d}$ such that
$R_{\\#}\mu_{r}=\mathcal{U}^{d}$, a known reference distribution of simple
structure. For our application $\mathcal{U}^{d}$ is a uniform distribution
over a unit hypercube $I^{d}$. It was shown in [4] and [8] that $\hat{R}$, an
estimate of $R$ , can be recovered from data without any parametric
assumptions.

If data is i.i.d distributed with laws $P_{1}$ and $P_{2}$ on each side of
change point, when figure 1 will exhibit relations between probability
distribution laws discussed so far.

$\begindc{\commdiag}[30]\obj(0,40)[PP1]{\mathbb{P}^{1}}\obj(40,40)[PP2]{\mathbb{P}^{2}}\obj(20,40)[PP]{\mathbb{P}}\obj(0,20)[P1]{P_{1}}\obj(20,20)[M]{\mu_{r}}\obj(40,20)[P2]{P_{2}}\obj(0,0)[RP1]{R_{\\#}P_{1}}\obj(40,0)[RP2]{R_{\\#}P_{2}}\obj(20,0)[U]{\mathcal{U}^{d}}\mor{P1}{M}{I^{r}_{\xi^{2}}}\mor{P2}{M}{I^{1-r}_{\xi^{1}}}[\atright,\solidarrow]\mor{PP1}{PP}{I^{r}_{\bullet}}\mor{PP2}{PP}{I^{1-r}_{\bullet}}[\atright,\solidarrow]\mor{P1}{RP1}{R}\mor{P2}{RP2}{R}\mor{RP1}{U}{I^{r}_{R(\xi^{2})}}\mor{RP2}{U}{I^{1-r}_{R(\xi^{1})}}[\atright,\solidarrow]\mor{M}{U}{R}\mor{PP1}{P1}{\xi^{1}}\mor{PP2}{P2}{\xi^{2}}\mor{PP}{M}{X}\enddc$
Figure 1:  Diagram with measures as objects and pushforwards as arrows . The
top row represents abstract measure spaces which are pushed forward to unknown
measures by observed Random variables. The bottom row exhibits state of the
probability laws after application of the vector rank function of $\mu_{r}$.

The diagram shows that application of vector rank function can be thought as a
probabilist’s ”change of basis”. The pragmatic value for change point problem
is in elimination of unknown distribution $\mu_{r}$ from the model.

It is obvious that the diagram depicted at fig. 1 commute. As the pushforward
acts as a linear map of signed measures defined over two different measurable
spaces. If $\alpha$ and $\beta$ are measures, $x,y\in\mathbb{R}$, $B$ is a
measurable set and $f$ is a measurable function, then

$f_{\\#}(x\alpha+y\beta)(B)=x\alpha(f^{-1}(B))+y\beta(f^{-1}(B))=xf_{\\#}\alpha(B)+yf_{\\#}\beta(B).$

Therefore, by linearity of the pushforward

$\mathcal{U}^{d}=R_{\\#}\mu_{r}=R_{\\#}(rP_{1}+(1-r)P_{2})=rR_{\\#}P_{1}+(1-r)R_{\\#}P_{2}.$

Thus, the diagram is correct.

The real result of change-point estimation will depend on measurable
difference between $P_{1}$ and $P_{2}$. We introduce notion of Kullback-
Leibler divergence in order to quantify this difference.

######  Definition 3.

Kullback-Leibler divergence between measures $P_{1}$ and $P_{2}$ admitting
densities $p_{1}$ and and $p_{2}$ respectively is

$\mathrm{KL}(P_{1},P_{2})=\int_{\Omega_{1}}\log\frac{p_{1}}{p_{2}}\,\mathrm{d}P_{1}$

Note that while
$P_{1},P_{2}\in\mathcal{P}_{\ll\lambda}\left(\mathbb{R}^{d}\right)$ they admit
densities. Now we formulate a sufficient condition on $R$ that ensures that
our transformations preserve divergence between distributions.

If $R$ is invertible almost everywhere $P_{1}$, when
$\mathrm{KL}(R_{\\#}P_{1},R_{\\#}P_{2})=\mathrm{KL}(P_{1},P_{2})$ Let $X$ be a
random variable such that $X\sim P_{1}$. Then, $R(X)\sim R_{\\#}P_{1}$ and by
invertibility a.e. of $R$ distribution $R_{\\#}P_{1}$ has density $p_{1}\circ
R^{-1}$ a.e. ; similar fact is also true for $P_{2}$ . Then,

$\mathrm{KL}(R_{\\#}P_{1},R_{\\#}P_{2})=\operatorname*{\mathbb{E}}\log\frac{p_{1}(R^{-1}RX)}{p_{2}(R^{-1}RX)}=\operatorname*{\mathbb{E}}\log\frac{p_{1}(X)}{p_{2}(X)}=\mathrm{KL}(P_{1},P_{2})$

## Formulation of the method

In this chapter our approach to construction of the vector rank function’s
estimates $\hat{R}_{n}$ is presented. The main goal of This approach is the
avoidance of approximation of the entire function $R$. Note, that even the
estimation that the aim of this task is actually not a precise estimation of
values $R(X_{i})$ for each observation $X_{i}$ but a construction of an
estimate that preserves the relative geometry of a sample $X$ along its change
points.

The work [4] presents continuous, semidiscrete, and discrete methods of
estimation of $R$. The Implementation of change point detection discussed in
this section is based solely on discrete approach. The reason behind this
choice is computational simplicity.

Vector rank function in general can be understood as $(P,\,mathbb{U})$-almost
surely homeomorphisms between supports of sample probability distribution $P$
and reference probability law of choice $\mathbb{U}$ ($R$ is defined almost
everywhere $P$ and continuous with respect to subset topology, and admits
similar continuous inverse existing almost everywhere $\mathbb{U}$) , such
that $R_{\\#}P=\mathbb{U}$ and both depth median and symmetries are preserved.
This symmetries and the concept of depth median need to be additionally
defined, which goes beyond of the scope of this paper. In our application
uniform measure over unit hypercube $I^{d}$ is selected as the reference $U$.
In case $d=1$ the cumulutive distribution function (cdf) of $P$ is also the
’vector’ rank function of the probability distribution of $P$.

The vector rank function is not unique in general. We use Monge-Kantarovich
vector rank developed in [4], which is defined as optimal transport between
$P$ and $U$

$R=\operatorname*{arg\,min}_{R:R_{\\#}P=\mathcal{U}^{d}}\int_{\mathbb{R}^{d}}d^{2}(R(x),x)\,\mathrm{d}P(x).$

Vector rank role of the optimal transport map $R$ can be intuitively justified
by the equivalence of the above optimization problem to the maximization of
the integral

$\int_{\mathbb{R}^{d}}\langle
R(x)-m(\mathcal{U}^{d}),x-m(P)\rangle\,\mathrm{d}P(x),$

where $m(\mathcal{U}^{d})$ and $m(P)$ stands for depth medians of
distributions $\mathcal{U}^{d}$ and $P$ respectively. Thus, the optimal
transport preserves geometric properties of distribution $P$ which can be
expressed by inner products of points relative to its center of symmetry. This
is what is understood as the relative geometry of the data ,and what we try to
preserve during vector rank estimation.

In this work we are focused on the discrete estimation of vector rank
function. Let $u$ be an equidistributed sequence of points in $I^{d}$. That
is, for every Lipschitz continuous function $f$ defined on $I^{d}$ it holds

$\lim_{n\to\infty}\frac{1}{n}\sum^{n}_{i=1}f(u_{i})=\int_{I^{d}}f(x)\,\mathrm{d}x.$
(3)

Then, for $T$ observations define
$Y_{i}=\hat{R}_{T}(X_{i})=u_{\sigma^{*}(i)}$, where

$\sigma^{*}=\operatorname*{arg\,min}_{\sigma\in
S^{T}}\sum^{T}_{i=1}\|X_{i}-u_{\sigma(i)}\|_{2}^{2}.$ (4)

In case convergence in (3) is understood as convergence a.s or even as
convergence in probability the sequence $u$ can be taken as an i.i.d sequence
of random variables or as an ergodic process with a uniform weak limit
distribution. However, this implementation is more suited for deterministic
form of $u$. Discrepancy can be understood as a natural measure of slackness
of condition (3) for a fixed value of $T$.

Classical (one-sided) geometric discrepancy of the sequence
$(u_{i})^{n}_{i=1}$ is described by the formula

$D(u)=\max_{x\in I^{d}}\left|\frac{\big{|}\\{u_{i}|1\leq i\leq
n\\}\cap[0,x]\big{|}}{T}-\prod^{d}_{j=1}x^{j}\right|.$ (5)

Note that (5) admits representation $D(u)=\|\varphi_{u}\|_{\infty}$ for a
certain function $\varphi_{u}$. This idea were used in [11] to introduce
generalized quadratic discrepancy with supremum-norm replaced by a certain
quadratic norm in a certain Sobolev Space. The result can be expressed as

$\Big{(}D_{2}^{\kappa}(u)\Big{)}^{2}=\\\
=\iint\eta(x,y)\,\mathrm{d}x\mathrm{d}y-\frac{2}{n}\sum^{n}_{i=1}\int\eta(x,u_{i})\,\mathrm{+}\frac{1}{n^{2}}\sum^{n}_{i,j=1}\eta(u_{j},u_{i}),$
(6)

where $\eta$ is the reproducing Hilbert kernel of the Sobolev space

$\eta(x,y)=\prod^{d}_{i=1}\Big{(}M+\beta^{2}\big{(}\kappa(x^{i})+\kappa(y^{i})+\frac{1}{2}B_{2}\big{(}x^{i}-y^{i}\mod
1\ \big{)}+B_{1}(x^{i})B_{1}(y^{i})\Big{)}$ (7)

with $\beta\in\mathbb{R}$ standing for a scale parameter and $\kappa$ standing
for a functional parameter with square-integrable derivative with
$\int^{1}_{0}\kappa=0$ and the value $M$ defined by

$M=1-\beta^{2}\int^{1}_{0}(\kappa^{\prime})^{2},$

and $B_{i}$ is the ith Bernoulli polynomial. Note, that in case $d=1$
statistic (5) turns into Kolmogorov-Smirnov statistic and (6) turns into
Cr/‘amer-Von Mizes statistic for uniform distribution test.

The function $\kappa$ is a functional parameter defining exact form of the
discrepancy. In this text we use star discrepancy produced by

$\kappa^{*}=\frac{1}{6}-\frac{1}{x^{2}}$

and the centred discrepancy produced by selecting

$\kappa^{c}=-\frac{1}{2}B_{2}\left(x-\frac{1}{2}\mod 1\right)$

The case of $u$ being a low-discrepancy sequence is a particularly well suited
for our application.

######  Definition 4.

low-discrepancy sequence is a deterministic sequence $u$ in $I^{d}$ designed
with a goal of minimizing value the $D^{\kappa}_{p}(u_{i})^{n}_{i=1}$ for each
natural number $n$.

Definition above is rather informal. However, it is postulated firmly that for
any low-discrepancy sequence $u$ property (3) holds and
$\lim_{n\to\infty}D^{\kappa}_{p}(u_{i})^{n}_{i=1}=0$ for any choice of $p$ and
$\kappa$. Thus, the rate of convergence of $D^{\kappa}_{p}(u_{i})^{n}_{i=1}$
can be understood as a measure of efficiency of a low-discrepancy sequence
$u_{i}$.

In our application we use Sobol sequence with grey code implementation. Sobol
sequence were introduced by Sobol [30] and the grey code implementation is due
[1]. Detailed investigation of the nature of this sequence goes beyond the
scope of this work. However, any other low-discrepancy sequence can be used.

For a Sobol sequence value $D^{\kappa}_{2}(u_{i})^{n}_{i=1}$ converges to zero
as $O(n^{-1+\varepsilon})$, where $\varepsilon$ depends on $d$ implicitly.
While convergence rate of discrepancy for a random uniform sequence is
$O(1/\sqrt{n})$, for $d\ll 200$ value of $\varepsilon<1/2$, which makes low-
discrepancy sequence rather efficient. However, scrambling and effective
dimension techniques can increase rate of convergence [12].

It can be established that the sequences $u$ and $X$ have no repeating values
almost surely, As our data is assumed to come from atomless distributions.
This means that the problem of finding Permutation $\sigma^{*}$ in (4) is the
optimal assignment problem. The Optimal assignment problem is the special case
the linear programming, which can be solved in $O(n^{3})$ time by the
Hungarian algorithm [17][32]. Amortizations based on applications of parallel
programming can improve computation complexity to $O(n^{2}\log n)$.
Approximate algorithms can be used for faster computations, for example [7].

Note, that even if sample $X$ has an i.i.d distribution, then the resulting
transport $Y$ is not independent itself. However, the covariance of the
elements is converging to zero as $T$ goes to infinity. Let $n\neq m$ be two
indices less or equal to $T$ and $i,j$ be coordinate indices. As
$(X_{k})^{T}_{k=1}$ is assumed to be i.i.d, it follows that $Y^{i}_{n}$ has a
discrete uniform distribution over $\\{u^{i}_{k}\\}^{T}_{k=1}$. Then,

$\operatorname*{Cov}(Y^{i}_{n},Y^{j}_{m})=\operatorname*{\mathbb{E}}Y^{i}_{n}Y^{j}_{m}-\operatorname*{\mathbb{E}}Y^{i}_{n}\operatorname*{\mathbb{E}}Y^{j}_{m}=\sum^{T}_{k=1}\sum^{T}_{l\neq
k}\frac{u^{i}_{k}u^{j}_{l}}{T(T-1)}-\sum^{T}_{k=1}\sum^{T}_{l=1}\frac{u_{k}^{i}u_{l}^{j}}{T^{2}}=\\\
=\sum^{T}_{k=1}\sum^{T}_{l\neq
k}\frac{u^{i}_{k}u^{j}_{l}}{T^{2}(T-1)}-\sum^{T}_{k=1}\frac{u^{i}_{k}u^{j}_{k}}{T^{2}},$

so by bounding from above and below with equidistribution property of $u$ in
the limit case

$\operatorname*{Cov}(Y^{i}_{n},Y^{j}_{m})\leq\sum^{T}_{k=1}\sum^{T}_{l=1}\frac{u^{i}_{k}u^{i}_{l}}{T^{2}(T-1)}\xrightarrow[T\to\infty]{}\lim_{T\to\infty}\frac{\operatorname*{\mathbb{E}}U\operatorname*{\mathbb{E}}U}{T-1}=\lim_{T\to\infty}\frac{1}{4(T-1)}=0,$
$\operatorname*{Cov}(Y^{i}_{n},Y^{j}_{m})\geq-\sum^{T}_{k=1}\frac{u^{i}_{k}u^{j}_{k}}{T^{2}}\geq-\sum^{T}_{k=1}\frac{u^{i}_{k}}{T^{2}}\xrightarrow[T\to\infty]{}\lim_{T\to\infty}-\frac{\operatorname*{\mathbb{E}}U}{T}=\lim_{T\to\infty}-\frac{1}{2T}=0;$

where $U$ is a uniform random variable on $[0,1]$. Thus, the value
$\operatorname*{Cov}(Y^{i}_{n},Y^{j}_{m})$ converges to zero. As variance of
$Y^{i}_{n}$ converges to the variance of a standard uniform distribution on
$[0,1]$ we will treat sample $Y$ as uncorrelated in asymptotic context under
assumption of i.i.d. distribution of $X$.

It can be easily seen that that quadratic discrepancy admits a degenerate
V-statistic representation with the kernel $\mathcal{K}$:

$\Big{(}D^{\kappa}_{2}(y)\Big{)}^{2}=\frac{1}{n^{2}}\sum^{n}_{i,j=1}\mathcal{K}(y_{i},y_{j}),$
(8)

assuming sample $y$ of $n$ elements. Properties of V-statistic produces the
asymptotic result

$n\Big{(}D^{\kappa}_{2}(y)\Big{)}^{2}\xrightarrow[n\to\infty]{d}V=\sum^{\infty}_{i=1}\lambda_{i}Z_{i}^{2}\sim\mathcal{V}(\lambda),$

where $Z\underset{\mathrm{i.i.d}}{\sim}\mathcal{N}(0,1)$ and $\lambda$ are
non-zero eigenvalues of the integral operator $\mathcal{A}$ defined by the
relation

$\mathcal{A}(f)(x)=\int^{1}_{0}f(y)\mathcal{K}(x,y)\mathrm{d}y$

It can be shown that $\mathcal{A}$ is in fact positive-semidefinite and trace-
class, which means that all $\lambda_{i}>0$ and that

$\sum^{\infty}_{n=1}\lambda_{i}<\infty$

Eigenvalues of $\mathcal{A}$ can be approximated by a finite collection of $N$
numbers $\hat{\lambda}$ with Nyström method. As it was shown in [5] the
twofold approximation of cdf of
$V_{|N}=\sum^{N}_{i=1}\widehat{\lambda}_{i}^{N}Z_{i}^{2}$ is possible for
fixed natural numbers $N$, $K$ and parameter $\alpha\in(0,1]$:

$\mathbb{P}(V_{|N}<x)\approx\\\
\approx\frac{1}{2}-\sum^{K}_{k=0}\frac{\sin\left(\frac{1}{2}\sum^{N}_{1=1}\arctan\left(2\left(k+\dfrac{1}{2}\right)\alpha\widehat{\lambda}^{N}_{i}\right)-\left(k+\frac{1}{2}\right)\alpha
x\right)}{\uppi\left(k+\frac{1}{2}\right)\prod^{N}_{i=1}\sqrt[4]{1+\left(2\left(k+\frac{1}{2}\right)\alpha\widehat{\lambda}^{N}_{i}\right)^{2}}}.$
(9)

This formula can be used for computing quantiles and critical values of
$\mathcal{V}(\lambda)$ with simple one-dimensional zero-seeking algorithm.

### Case of A Single Change Point

In this chapter an approach for detecting a single change point is presented.
We impose a model assumptions that for a fixed $\tau$ change point $\theta$
either satisfy $\tau<\theta<T-\tau$ or it does not exist . This value $\tau$
can be selected in a such way that $1\ll\tau\ll T$ and be used as a sliding
window bandwidth defining a ’control chart’-like object which we refer to as
diphoragram.

######  Definition 5.

empirical sliding diphoragram for change point data $(X_{t})^{T}_{i=1}$ is
defined by

$\widehat{\Delta}_{t}^{T,\tau}=\Big{(}D^{\kappa}_{2}\big{(}\hat{R}_{T}(X_{i})\big{)}^{t+\tau}_{i=t}\Big{)}^{2}=\Big{(}D^{\kappa}_{2}(Y_{i})^{t+\tau}_{i=t}\Big{)}^{2}$

and ideal sliding diphoragram by

$\Delta_{t}^{T,\tau}=\Big{(}D^{\kappa}_{2}\big{(}R_{\mu_{r}}(X_{i})\big{)}^{t+\tau}_{i=t}\Big{)}^{2}$

for $t$ in range from $1$ to $T^{\prime}_{\tau}=T-\tau$.

If $\kappa$ is a continuous function, then
$\widehat{\Delta}_{t}^{T,\tau}\xrightarrow[]{P}\Delta_{t}^{T,\tau}$ as
$T\to\infty$. With this condition discrepancy is a continuous function of
data. So, convergence of vector ranks proved in [4] implies convergence of
diphoragrams. Note, that with kernel representation charts admit a difference
representation

$\widehat{\Delta}_{t+1}^{T,\tau}-\widehat{\Delta}_{t}^{T,\tau}=\frac{1}{\tau^{2}}\left(\sum^{t+1+\tau}_{i=t+1}\mathcal{K}(Y_{i},Y_{t+\tau+1})-\sum^{t+\tau}_{i=t}\mathcal{K}(Y_{i},Y_{t})\right).$

From the computational standpoint this means that computation of the whole
time-series $\widehat{\Delta}^{T}$ and $\widehat{\Delta}^{T,\tau}$ takes only
quadratic time $O(dT^{2})$ in number in observation. Moreover, if crisp
optimal transport with low-discrepancy sequence $u$ is used in estimation of
vector rank function, then $A_{Y}=A_{u}$, so all values can be precomputed.

For the application to the change point problem consider two increasing
sequences of integers $T_{n}$ and $\tau_{n}$ such that
$\frac{T_{n}}{\tau_{n}}=a>1$ for all $n\in\mathbb{N}$ are constructed. Then as
$Y_{i}$ are assumed to be independent $\widehat{\Delta}_{t}^{T_{n},\tau_{n}}$
and $\widehat{\Delta}_{t+\tau_{n}}^{T_{n},\tau_{n}}$ are also independent
random variables.

######  Definition 6.

mean sliding discrepancy for set sample $(Y_{i})_{i=1}^{T_{n}}$ is computed by

$\overline{\Delta}_{n}=\frac{\tau_{n}^{2}}{T_{n}}\sum^{a-1}_{j=0}\widehat{\Delta}_{1+j\tau_{n}}^{T_{n},\tau_{n}}$

Note, that the mean sliding discrepancy is not the same as the discrepancy of
the whole sample. By independence, in case $H_{0}$ holds, as $n\to\infty$

$\overline{\Delta}_{n}\xrightarrow[T\to\infty]{d}\sum^{a}_{j=1}\sum^{\infty}_{i=1}\frac{\lambda_{i}}{a}Z_{i,j}^{2}\sim\mathcal{V}\left(\left(\frac{\lambda_{i}}{a}\right)^{a}_{j=1}\right)^{\infty}_{i=1},$

where $Z\underset{\mathrm{i.i.d}}{\sim}\mathcal{N}(0,1)$. Otherwise, there
will be a sliding discrepancies $\widehat{\Delta}_{t}^{T_{n},\tau_{n}}$
sampled form the non-uniform data which goes to $\infty$ in probability as
$n\to\infty$. So, the whole sum $\overline{\Delta}_{n}\xrightarrow[n\to
infty]{\mathbb{P}}\infty$.

In case the single change point exists the statistic
$\widehat{\Delta}^{T_{n},\tau_{n}}$ is expected to attain the minimal value at
such moment of time $t^{*}$ when the empirical distribution of
$(Y_{t})^{t^{*}+\tau_{n}}_{t=t^{*}}$ approaches the empirical distribution of
the whole sample, as the empirical distribution of the whole sample is
converging to the $\mathcal{U}^{d}$. Let

$t^{*}_{n}=\operatorname*{arg\,min}_{1\leq t\leq
T_{n}-\tau_{n}}\widehat{\Delta}_{t}^{T_{n},\tau_{n}}.$

It is expected that the ratio of numbers of elements from both sides of change
point in the subsample used in the computation of
$\widehat{\Delta}_{t^{*}_{n}}^{T_{n},\tau_{n}}$ and in the whole sample are
equal. This can be represented in the algebraic relation

$\frac{\tilde{\theta}_{n}-t^{*}_{n}}{\tau_{n}}=\frac{\tilde{\theta}_{n}}{T_{n}}=\hat{r}_{n},$

where $\tilde{\theta}_{n}$ is the produced estimate of the change point. This
provides the expression for the estimate

$\tilde{\theta}_{n}=\frac{t^{*}_{n}}{1-a^{-1}}.$

We accept $H_{0}$ for a fixed $\gamma$-level in case

$p_{n}=\mathbb{P}\Big{(}V_{|N}\leq\overline{\Delta}_{n}\Big{)}<1-\gamma,$ (10)

where the cdf is numerically estimated as in (9) for some fixed parameters $N$
and $K$. Otherwise reject $H_{0}$ and state that there was a change most
probable at $\hat{\theta}_{n}$.

In order to reason about properties of the estimate $\hat{\theta}_{n}$. Let
$\mathcal{M}(I^{d})$ denote space of finite signed measures other the
hypercube $I^{d}$.

######  Definition 7.

Space of nonuniformities $\mathcal{D}$ is defined as a quotient of the real
vector spaces

$\mathcal{D}=\frac{\mathcal{M}(I^{d})}{\mathbb{R}\mathcal{U}^{d}}$

endowed with a Hilbert space structure by inner product defined for
$[\nu],[\mu]\in\mathcal{D}$ by

$\big{\langle}[\nu],[\mu]\big{\rangle}=\int\int\mathcal{K}(x,y)\,\mathrm{d}\nu(x)\,\mathrm{d}\mu(y)$
(11)

Note, that the relation (11) is well defined as for all $a,b\in\mathbb{R}$

$\Big{\langle}\big{[}\nu+a\mathcal{U}^{d}\big{]},\big{[}\mu+b\mathcal{U}^{d}\big{]}\Big{\rangle}=\int\int\mathcal{K}(x,y)\,\mathrm{d}\nu(x)\,\mathrm{d}\mu(y)+b\int\int\mathcal{K}(x,y)\,\mathrm{d}\nu(x)\,\mathrm{d}y+\\\
+a\int\int\mathcal{K}(x,y)\,\mathrm{d}x\,\mathrm{d}\mu(y)+ab\int\int\mathcal{K}(x,y)\,\mathrm{d}x\,\mathrm{d}y=\int\int\mathcal{K}(x,y)\,\mathrm{d}\nu(x)\,\mathrm{d}\mu(y),$

as $\int\mathcal{K}(x,y)\,\mathrm{d}x=0$ for any value of $y$ and
$\mathcal{K}$ is symmetric and is indeed an inner product as $\mathcal{K}$ is
a positive-definite kernel. As the line $\mathbb{R}\mathcal{U}^{d}$ intersects
simplex of probability measures $\mathcal{P}(I^{d})$ only in one point
($\mathcal{U}^{d}$ itself) the natural projection $\mu\mapsto[\mu]$ is
injective on $\mathcal{P}(I^{d})$, so we denote a nonuniformity arising from
each $\mu\in\mathcal{P}(I^{d})$ just as $\mu$. If every subsample
$(Y_{i})_{i=t}^{t+\tau}$ is associated with an empirical measures
$\hat{\mu}_{t}=\sum_{i=t}^{t+\tau_{n}}\delta_{Y_{i}}$, then

$\widehat{\Delta}_{t}^{T_{n},\tau_{n}}=\|\hat{\mu}_{t}\|^{2}_{\mathcal{D}}=\langle\hat{\mu}_{t},\hat{\mu}_{t}\rangle.$

By construction of the vector rank function
$\|rR_{\\#}P_{1}+(1-r)R_{\\#}P_{2}\|_{\mathcal{D}}=0$, hence
$rR_{\\#}P_{1}+(1-r)R_{\\#}P_{2}=0$ in $\mathcal{D}$. This produces result

$R_{\\#}P_{2}=-\frac{r}{1-r}R_{\\#}P_{1}\Longrightarrow\Big{\|}R_{\\#}P_{2}\Big{\|}_{\mathcal{D}}^{2}=\frac{r^{2}}{(1-r)^{2}}\Big{\|}R_{\\#}P_{1}\Big{\|}^{2}_{\mathcal{D}}\Longrightarrow\\\
\Longrightarrow
r=\frac{\|R_{\\#}P_{2}\|_{\mathcal{D}}}{\|R_{\\#}P_{1}\|_{\mathcal{D}}+\|R_{\\#}P_{2}\|_{\mathcal{D}}}$
(12)

Considering that $r=1/2$ it follows that $[R_{\\#}P_{1}]=-[R_{\\#}P_{2}]$, so
magnitude of the discrepancy will have the same distribution for sample with
equal proportions of elements with distributions of $R_{\\#}P_{1}$ and
$R_{\\#}P_{2}$. If $t^{*}=(1-a^{-1})\theta$, then it can be shown that
estimate $\hat{\theta}$ is unbiased:

$\operatorname*{\mathbb{E}}\hat{\theta}_{n}=\dfrac{\sum^{T_{n}-\tau_{n}}_{t=1}t\mathbb{P}\big{(}t=\operatorname*{arg\,min}_{t^{\prime}}\widehat{\Delta}_{t^{\prime}}^{T_{n},\tau_{n}}\big{)}}{1-a^{-1}_{n}}=\frac{t^{*}+B^{+}_{n}-B^{-}_{n}}{1-a^{-1}_{n}}=\frac{t^{*}}{1-a^{-1}_{n}}=\theta=\frac{T_{n}}{2},$

where

$B^{-}_{n}=\sum^{\lfloor
r(T_{n}-\tau_{n})\rfloor}_{t=1}t\mathbb{P}\Big{(}t^{*}-t=\operatorname*{arg\,min}_{t^{\prime}}\widehat{\Delta}_{t^{\prime}}^{T_{n},\tau_{n}}\Big{)},\quad\\\
B^{+}_{n}=\sum^{\lfloor(1-r)(T_{n}-\tau_{n})\rfloor}_{t=1}t\mathbb{P}\Big{(}t^{*}+t=\operatorname*{arg\,min}_{t^{\prime}}\widehat{\Delta}_{t^{\prime}}^{T_{n},\tau_{n}}\Big{)}.$

However, then $r\neq 1/2$ the estimate $\hat{r}$ is projected to be biased
towards $1/2$ as the value $B^{+}$ will only increase and the value $B^{-}$
will only decrease as $r$ decreases. Otherwise, increase of $r$ increases the
value $B^{-}$ too, and decreases the value of $B^{+}$. In order to proof
consistency of estimator $\hat{r}_{n}$ we introduce a family of sequences for
each $s\in((1-a^{-1})^{-1},1)$

$\tilde{\Delta}^{s}_{n}=\widehat{\Delta}_{t^{\prime}}^{T_{n},\tau_{n}}\quad\text{having}\quad
t^{\prime}=\operatorname*{arg\,min}_{1\leq t\leq
T_{n}-\tau_{n}}\left|s-\frac{t}{T_{n}(1-a^{-1})}\right|$

Then by the weak convergence of vector ranks for every $s$ the sequence
converges to a value:

$\tilde{\Delta}^{s}_{n}\xrightarrow[n\to\infty]{\mathbb{P}}\Delta^{s}=\Bigg{\|}1_{s<\frac{r-a^{-1}}{1-a^{-1}}}(s)R_{\\#}P_{1}+1_{s>\frac{r}{1-a^{-1}}}(s)R_{\\#}P_{2}+\\\
+1_{\frac{r-a^{-1}}{1-a^{-1}}\leq
s\leq\frac{r}{1-a^{-1}}}(s)\left(\frac{r-s(1-a^{-1})}{a^{-1}}R_{\\#}P_{1}+\frac{a^{-1}-r+s(1-a^{-1})}{a^{-1}}R_{\\#}P_{2}\right)\Bigg{\|}_{\mathcal{D}}.$

Otherwise $\tilde{\Delta}^{r}_{n}\xrightarrow{\mathbb{P}}0$. We can assert by
structural uniformity of $\tilde{\Delta}_{n}$ that
$\tilde{\Delta}\xrightarrow{\mathbb{P}}\Delta$ in Skorohod’s topology. Thus,
$\hat{r}_{n}=\operatorname*{arg\,min}_{s}\Delta_{n}^{s}\xrightarrow{\mathbb{P}}\operatorname*{arg\,min}_{s}\Delta^{s}=r$
providing convergence in probability. This suggests that the bias of
$\hat{r}_{n}$ can be bounded by some sequence $\beta$ with convergence
$\beta_{n}\xrightarrow[n\to\infty]{}0$:

$\frac{|B^{+}_{n}-B^{-}_{n}|}{T_{n}}\leq\beta_{n}$

Without loss of generality assume that $r<1/2$. Then, we separate positive
bias into mixing and non-mixing parts:

$\frac{|B^{+}_{n}-B^{-}_{n}|}{T_{n}}<\frac{B^{+}_{n}}{T_{n}}<\\\
<\sum_{t=1}^{\lfloor(1-r)\tau_{n}\rfloor}\frac{t\mathbb{P}\Big{(}\widehat{\Delta}_{t^{\prime}}^{T_{n},\tau_{n}}<\tilde{\Delta}^{r}_{n}\Big{)}}{T_{n}}+\sum_{t=\lfloor(1-r)\tau_{n}\rfloor+1}^{T_{n}-\tau_{n}-t^{*}}\frac{t\mathbb{P}\Big{(}\widehat{\Delta}_{t^{\prime}}^{T_{n},\tau_{n}}<\tilde{\Delta}^{r}_{n}\Big{)}}{T_{n}}.$

For non-mixing part apply Markov inequality for each fixed value of
$\tilde{\Delta}^{r}_{n}$ and some value $\zeta>1$

$\mathbb{P}\Big{(}\widehat{\Delta}_{t^{\prime}}^{T_{n},\tau_{n}}<\tilde{\Delta}^{r}_{n}\Big{)}=\mathbb{P}\Big{(}\big{(}\widehat{\Delta}_{t^{\prime}}^{T_{n},\tau_{n}}\big{)}^{-\zeta}>\big{(}\tilde{\Delta}^{r}_{n}\big{)}^{-\zeta}\bigg{)}\leq\big{(}\tilde{\Delta}^{r}_{n}\big{)}^{\zeta}\operatorname*{\mathbb{E}}\bigg{(}\big{(}\widehat{\Delta}_{t^{\prime}}^{T_{n},\tau_{n}}\big{)}^{-\zeta}\bigg{)}.$

It is possible to use inverse as $\tilde{\Delta}_{n}>0$ for each $n$. With
this inequality the non-mixing part can be bounded

$\sum_{t=\lfloor(1-r)\tau_{n}\rfloor+1}^{T_{n}-\tau_{n}-t^{*}}\frac{t\mathbb{P}\Big{(}\widehat{\Delta}_{t^{\prime}}^{T_{n},\tau_{n}}<\tilde{\Delta}^{r}_{n}\Big{)}}{T_{n}}\leq
T_{n}\operatorname*{\mathbb{E}}\bigg{(}\big{(}\tilde{\Delta}^{r}_{n}\big{)}^{\zeta}\bigg{)}\operatorname*{\mathbb{E}}\bigg{(}\big{(}\tilde{\Delta}^{1}_{n}\big{)}^{-\zeta}\bigg{)}\leq\\\
\leq
C_{r}T_{n}\operatorname*{\mathbb{E}}\bigg{(}\big{(}\tilde{\Delta}^{r}_{n}\big{)}^{\zeta}\bigg{)}.$

As $\tilde{\Delta}^{r}_{n}\to 0$ with the rate $O(\tau_{n}^{-1+\varepsilon})$
it is possible to specify a constant $\zeta$ in such a way that
$T_{n}\operatorname*{\mathbb{E}}((\tilde{\Delta}^{r}_{n})^{\zeta})\downarrow
0$ and the expectation
$\operatorname*{\mathbb{E}}((\tilde{\Delta}^{1}_{n})^{-\zeta})$ approaches
some constant value by weak convergence and therefore can be bounded. Thus,
the non-mixing bias approaches zero as $n$ goes to infinity. Note, that all
moments of discrepancy exists as it is bounded on a positive interval for each
$n$. This method can be extended in order to show that the estimate is
asymptotically unbiased under an assumption of change point slackness,

######  Definition 8.

Sequence of diphoragrams $\widehat{\Delta}_{t}^{T_{n},\tau_{n}}$ has slackness
property iff for some constant $\gamma\in(0,a^{-1})$ and for all $n$ large
enough there are time points $t^{**}_{n}>t^{*}_{n}$ such that :

$\frac{t^{**}_{n}-t^{*}_{n}}{T_{n}}\geq\gamma,\quad\sum_{t=1}^{t^{**}_{n}}t\mathbb{P}\Big{(}t=\operatorname*{arg\,min}_{t^{\prime}}\widehat{\Delta}_{t^{\prime}}^{T_{n},\tau_{n}}\Big{)}-B^{-}_{n}\leq
0.$

Then, it is possible to construct a similar bound

$\frac{|B^{+}_{n}-B^{-}_{n}|}{T_{n}}\leq
T_{n}\operatorname*{\mathbb{E}}\Bigg{(}\big{(}\tilde{\Delta}^{r}_{n}\big{)}^{\lambda}\operatorname*{\mathbb{E}}\bigg{(}\big{(}\tilde{\Delta}^{\gamma+r}_{n}\big{)}^{-\lambda}\Big{|}\tilde{\Delta}^{r}_{n}\bigg{)}\Bigg{)},$

which converges to zero at infinity.

### alternative methods

One of the negative properties of the method described in the previous chapter
is the requirement of specification of bandwidth $\tau$, which prevents change
point detection in the proximity of the limit points $1$ and $T$. However, as
space of nonuniformities $\mathcal{D}$ is a metric space, it is possible to
determine change point by maximizing distance between two empirical measures
$\mathrm{dist}(\hat{\theta})=\|\hat{\mu}^{+}_{\hat{\theta}}-\hat{\mu}^{-}_{\hat{\theta}}\|^{2}_{\mathcal{D}}$,
where

$\hat{\mu}^{-}_{\hat{\theta}}=\frac{1}{\widehat{\theta}}\sum^{\hat{\theta}}_{t=1}\delta_{Y_{t}},\quad\hat{\mu}^{+}_{\hat{\theta}}=\frac{1}{T-\widehat{\theta}}\sum^{T}_{t=\hat{\theta}+1}\delta_{Y_{t}}.$
(13)

Then, by definition of the norm the distance is computed as

$\mathrm{dist}\Big{(}\hat{\theta}\Big{)}=\frac{1}{(\hat{\theta})^{2}}\sum^{\hat{\theta}}_{n,m=1}\mathcal{K}(Y_{n},Y_{m})-\\\
-\frac{2}{\hat{\theta}(T-\hat{\theta})}\sum^{\hat{\theta}}_{n=1}\sum^{T\theta}_{m=\hat{\theta}+1}\mathcal{K}(Y_{n},Y_{m})+\frac{1}{(T-\hat{\theta})^{2}}\sum^{T}_{n,m=\hat{\theta}+1}\mathcal{K}(Y_{n},Y_{m}),$

leading to a change point estimate
$\hat{\theta}=\operatorname*{arg\,max}_{\tilde{\theta}}\mathrm{dist}(\tilde{\theta})$.
As empirical measures of $\eqref{measures}$ will converge to some points of
interval $[R_{\\#}P_{1},R_{\\#}P_{2}]$ as $T_{n}$ goes to infinity, the
estimate $\hat{\theta}_{n}$ converges to true value $\theta$ in probability.

Change points detection methods of this forms were explored in the work [21].
Thus, we will not explore it in further depth. The important properties of
this statistic is that it still can be computed in $O(T^{2})$ time and that
$\mathrm{dist}(\hat{\theta})$ is a U-statistic.

Geometric idea of (12) suggests that the value

$\varsigma\Big{(}\hat{\theta}\Big{)}=\frac{\hat{\theta}}{T}-\frac{\|\hat{\mu}^{+}_{\hat{\theta}}\|_{\mathcal{D}}}{\|\hat{\mu}^{-}_{\hat{\theta}}\|_{\mathcal{D}}+\|\hat{\mu}^{+}_{\hat{\theta}}\|_{\mathcal{D}}}$

approaches zero as data size grows to infinity iff $\hat{r}$ approaches the
true ratio $r$. Thus, if change point exists, then change point can be
estimated as
$\hat{\theta}=\operatorname*{arg\,min}_{\tilde{\theta}}|\varsigma(\tilde{\theta})|$.
Or alternatively compute $\hat{r}^{m}$ by applying iteration

$\hat{r}^{m}=\frac{\|\hat{\mu}^{+}_{r^{m-1}}\|_{\mathcal{D}}}{\|\hat{\mu}^{-}_{r^{m-1}}\|_{\mathcal{D}}+\|\hat{\mu}^{+}_{r^{m-1}}\|_{\mathcal{D}}},$

where $\mu^{+}_{r^{m-1}}$ and $\mu^{-}_{r^{m-1}}$ are empirical measures
corresponding to the estimates obtained at previous iterations. The initial
value $\hat{r}^{0}$ can be selected to be equal to $1/2$.

## Multiple Change Points

In this chapter the situation of $K$ possible consecutive change points
$\theta=(\theta^{1},\ldots,\theta^{K})$ in the data is considered. Now, $r$
denotes a list of positive values

$r=(r^{1},\ldots,r^{K})=\left(\frac{\theta^{1}}{T},\frac{\theta^{2}-\theta^{1}}{T},\ldots,\frac{\theta^{K}-\sum^{K-1}_{k=1}\theta^{k}}{T}\right),$

which can be understood as the first $K$ coefficients in the convex
combination of $K+1$ probability distributions $P_{1},\ldots,P_{K+1}$. That
is, define

$\mu_{r}=\sum^{K}_{k=1}r^{k}P_{k}+\left(1-\sum^{K}_{k=1}r^{k}\right)P_{K+1}.$

Furthermore, estimates $\hat{\theta}$ and $\hat{r}$ are treated as lists of
corresponding structure.

By definition of vector rank function $R_{\\#}\mu_{r}=\mathcal{U}^{d}$.
However, in order to apply methods similar to ones developed in the previous
chapter we need one more property of the model.

######  Definition 9.

$K$ probability measures $P_{1},\ldots P_{K}$ are said to be convexly
independent if for all lists of $K$ coefficient
$\lambda\in\mathbb{R}^{K}_{+}$, such that $\sum^{K}_{k=1}\lambda_{k}=1$, for
each $i$ equality $P_{i}=\sum^{K}_{k=1}\lambda_{k}P_{k},$ holds only if
$\lambda_{j}=\delta_{i,j}$ for each $j$, where $\delta_{i,j}$ is the Kronecker
delta.

Assume that the true value of $K\in\mathfrak{K}$. Then, by construction
$\widehat{K}^{\mathrm{SMA}}\leq K$ with probability $\alpha$.

It can be postulated, that $0$ of $\mathcal{D}$ lies in the convex hull of
$R_{\\#}P_{i}$. That is $0\in\mathrm{conv}\\{R_{\\#}P_{k}\\}^{K+1}_{k=1}.$ It
suggests that there are projections of $0$ to the edges of the convex polytope
$\pi_{k}0\in[R_{\\#}P_{k},R_{\\#}P_{k+1}]$ which minimizes
$\|\pi_{k}0\|_{\mathcal{D}}$. Hence, if $\tau$ is taken to be small enough
then local minimas of $\widehat{\Delta}^{T,\tau}$ will happen in the proximity
of the change points with high probability.

The problem with this method is that proportions of points from different
sides of a change point at the local minimise are given by the relation

$\pi_{k}0=\frac{\|R_{\\#}P_{k+1}\|^{2}_{\mathcal{D}}-\langle
R_{\\#}P_{k},R_{\\#}P_{k+1}\rangle_{\mathcal{D}}}{d^{2}_{\mathcal{D}}(R_{\\#}P_{k},R_{\\#}P_{k+1})}R_{\\#}P_{k}+\\\
+\frac{\|R_{\\#}P_{k}\|^{2}_{\mathcal{D}}-\langle
R_{\\#}P_{k},R_{\\#}P_{k+1}\rangle_{\mathcal{D}}}{d^{2}_{\mathcal{D}}(R_{\\#}P_{k},R_{\\#}P_{k+1})}R_{\\#}P_{k+1},$

which can not be recovered from the diphoragram. For this reason we propose an
iterative procedure. Let $\hat{\mu}_{\hat{\theta}^{k},\hat{\theta}^{k+1}}^{n}$
denote empirical measures of points from the interval bounded by change point
estimates $[\hat{\theta}^{k}]^{n}$ and $[\hat{\theta}^{k+1}]^{n}$ .

  1. 1.

start with $\mathcal{T}_{1}=\mathcal{T}$.

  2. 2.

for each $k\leq K$ select
$t^{*}_{k}=\operatorname*{arg\,min}_{t\in\mathcal{T}_{k}}\tilde{\Delta}^{T,\tau}_{t}$
and update $T_{k+1}=\\{t\in\mathcal{T}_{k}:|t-t^{*}_{k}|>\tau\\}$.

  3. 3.

make initial change point estimation with blind adjustment
$\Big{[}\hat{\theta}^{k}\Big{]}_{1}=t^{*}_{k}+\tau/2$.

  4. 4.

readjust change points for a fixed number of iterations $N$ with nth plus one
readjustment being

$\Big{[}\hat{\theta}^{k}\Big{]}_{n+1}=t^{*}_{k}+\Big{[}\hat{\lambda}^{k}_{1}\Big{]}_{n}\tau,$

where

$\Big{[}\hat{\lambda}^{k}_{1}\Big{]}_{n}=\frac{\Big{\|}\mu^{n}_{\hat{\theta}_{k},\hat{\theta}_{k+1}}\Big{\|}^{2}_{\mathcal{D}}-\Big{\langle}\mu^{n}_{\hat{\theta}_{k-1},\hat{\theta}_{k}},\mu^{n}_{\hat{\theta}_{k},\hat{\theta}_{k+1}}\Big{\rangle}_{\mathcal{D}}}{d^{2}_{\mathcal{D}}\Big{(}\mu^{n}_{\hat{\theta}_{k-1},\hat{\theta}_{k}},\mu^{n}_{\hat{\theta}_{k},\hat{\theta}_{k+1}}\Big{)}}$

with surrogate change points being $\Big{[}\hat{\theta}^{0}\Big{]}_{n-1}=0$
and $\Big{[}\hat{\theta}^{K+1}\Big{]}_{n-1}=T$ .

In order to approach problem of model misspecification we apply smallest
accepted model (SMA) methodology. For a collection of $\widehat{K}$ we acquire
minimizing time points $t^{*}$. Then, test for change points in the intervals
bounded by $\hat{\theta}_{i}$ and $\hat{\theta}_{i+1}$ with surrogate values
as above by computing a p-value approximations.

Thus, the model can be estimated while new local minimizers are being
recovered and the process can be stopped as minimal accepted value of
$\hat{K}$ has been achieved.

## Computational Results

In order to conduct computational experiments the methods discussed in the
previous chapter were implemented in Python programming language with use of
numpy and scipy libraries. We use sobol_seq package in order to generate Sobol
sequences.

### Simulations with Zero Change Points

Figure 2: Example of a control chart for 1000 observations of 5-dimensional
data with standard normal distribution without change points computed by star
kernel with $\tau=100$. Abscissa represents values of shifted time $t$ in a
range from 1 to 900 and ordinate stands for value of $\hat{\Delta}_{t}$.
Successfully, no change points were detected. Note, that whole process
$\hat{\Delta}_{t}$ behaves as a martingale with ’small’ expected value.

For simulations with zero change points we are interested in measuring
statistical significance or confidence of the test statistic $T$ which can be
understood as

$\mathrm{significance}(T)=\mathbb{P}(T\,\textrm{rejects}\,H_{0}|H_{0}\,\textrm{is
true}),\quad\mathrm{confidence}(T)=\mathbb{P}(T\,\textrm{accepts}\,H_{0}|H_{0}\,\textrm{is
true}).$

In simulations with zero change points $H_{0}$ obviously is true. So, for a
run of $n$ simulations we estimate confidence as

$\widehat{\mathrm{conf}}(T,n)=\frac{\\#\\{\textrm{simulations with no change
points detected}\\}}{n}.$

We run simulations without change points and vary certain fixed parameter
while measuring confidence and inverse $p$-value for each value of parameter.
The only nontrivial results were obtained for change in data dimension $d$.

Figure 3: Change of mean inverse $p$-value with the growth of dimensionality.
We ran 20 simulations for each value $d$ of dimension in range from 1 to 12.
Each simulation produced sample 200 standard normal variables and no change
points. Then hypothesis were tested with significance parameter $\alpha=0.1$
and diphoragrams built with $\tau=30$. Observe that the inverse $p$-value
decreases from $1$ to $6$ which can be understood as increase in statistical
confidence in the absence of change points. Then it goes to the value of $0.5$
which can be interpreted as complete statistical indetermination due to
complete lack of information caused by high data dimension relative to the
sliding window bandwidth.

Additional experiments were conducted with growing variance and changes in
covariance structure of the observations, however no dependencies were
identified. This could be due to stabilizing effect of vector rank functions.

## Simulations with One Change Point

Figure 4: Example of a diphoragram for 2000 observations of 5-dimensional data
with standard normal distribution from one side of change point at
$\theta=1000$ and normal distribution with mean of value 5 computed by star
kernel with $\tau=100$. Abscissa represents values of shifted time $t$ in a
range from 1 to 1900 and ordinate stands for value of $\hat{\Delta}_{t}$.
Successfully, change point was estimated at $\hat{\theta}=1000$ with
$t^{*}=950$. Note, that the whole process $\hat{\Delta}_{t}$ behaves as a
martingale with ’big’ expected value in the distance of change point and goes
for a ’swoop’ in proximity of $\theta$.

While running a simulations with one change point we are naturally interested
in the estimation of statistical power of the test $T$ which can be understood
as

$\mathrm{power}(T)=\mathbb{P}(T\,\textrm{rejects}\,H_{0}|H_{0}\,\textrm{is
false}),$

and can be estimated for $n$ simulations as

$\widehat{\mathrm{pow}}(T,n)=\frac{\\#\\{\textrm{simulations with change
points detected}\\}}{n}.$

We investigate dependence of $\widehat{\mathrm{pow}}(T,n)$ on differences
between distributions from opposite sides of the change point.

Figure 5: Change of inverse p-value and statistical power with the growth of
difference in the mean. We ran 20 simulations for each values of the
difference in the range from 0.2 to 1 with the step equal to 0.1. Each
simulation contained a sample of 200 normally distributed variables with $d=3$
and $\theta=100$ with corresponding difference in the mean of distributions
from both sides of the change point. The hypothesis was tested with
$\alpha=0.1$ and bandwidth $\tau=30.$ The graph shows that the star kernel
outperforms symmetric kernel at all values of the difference.  Figure 6:
Change of inverse p-value ans statistical power with the growth of difference
in the variance. We ran 20 simulations for each values of the difference in
the range from 2 to 10 with the step equal to 0.5. Each simulation contained a
sample of 200 normally distributed variables with $d=5$ and $\theta=100$ with
corresponding difference in the variance of distributions from both sides of
the change point. The hypothesis was tested with $\alpha=0.1$ and bandwidth
$\tau=30.$ The graph shows that the star kernel outperforms symmetric kernel
at all values of the difference.

The figures shows that star discrepancy outperforms symmetric discrepancy in
detecting change both in expectation and in variance. However empirical
results in [18] indicates that in some situations symmetric discrepancy may
turn out tob be a better tool.

In the situation of existing change point not only power of the test is of
interest but also a precision of change point estimations. As it was shown in
previous chapter bias of a change point estimate increases as true ratio $r$
of distributions in the sample diverges from the value of $1/2$. We provide a
computational illustrations:

Figure 7:  Discrepancy diphoragram for 300 normal i.i.d observations with
$d=3$ with change in expectation at $\theta=200$. Abscissa represents values
of shifted time $t$ in a range from 1 to 250 and ordinate stands for value of
$\hat{\Delta}_{t}$ constructed with $\tau=50$. Change point was detected at
location $\hat{\theta}=197$, which can be interpreted as a weak drift towards
midpoint. Note, that the relation of mean values in martingale parts of
$\hat{\Delta}_{t}$ at different sides of the change point corresponds to the
nonuniformity space theory presented in the previous chapter.  Figure 8:
Change of error $\theta-\hat{\theta}$ and absolute error
$|\theta-\hat{\theta}|$ with shift of value of $r$. We ran 20 simulations for
each values of the ration in the range from $1/2$ to $13/60$ with the step
equal to $1/60$. Each simulation contained a sample of $r*300$ normally
distributed variables with zero expectation at one side of change point and
$(1-r)*300$ variables uniformly distributed in $[-1,1]$ square with $d=2$ and
$\theta=r*300$ . Change points were estimated with $\tau=50$. The graph shows
that error indeed grows in the distance of the midpoint and that estimate has
strong drift towards $1/2$ for such ratios.

### Simulations with Multiple Change Points

For case of multiple change point we provide only example with diphoragrams:

Figure 9:  Discrepancy diphoragrams for 300 observations with $d=3$. Sample
has distribution model
$P_{1}=\mathcal{N}((3,3,3),I),P_{2}=\mathcal{U}[10,20]^{3}$ and
$P_{3}=\mathcal{N}(-(3,3,3),I)$ with $\theta_{1}=120$ and $\theta_{2}=240$.
Abscissa represents values of shifted time $t$ in a range from 1 to 250 and
ordinate stands for value of $\hat{\Delta}_{t}$ constructed with $\tau=50$.
Change point was detected at locations $\hat{\theta}_{1}=111$ and
$\hat{\theta}_{2}=237$. However, after application of iterative correction
improved estimate $\hat{\theta}_{1}^{\prime}=119$ and
$\hat{\theta}_{2}^{\prime}=241$ were acquired.

## Examples with financial data:

In order to provide examples, which are not artificially constructed,
financial data similar to one in [31]. This data was acquired from Data
library of Keneth R. French. It contains mean monthly returns from five
portfolios each composed of one of five industrial sectors in the USA, which
are:

  1. (A)

Finance,

  2. (B)

Manufacturing,

  3. (C)

Retail, wholesale and some services,

  4. (D)

Utilities,

  5. (E)

Other.

This provides data dimension of $d=5$ and total number of observations of
$T=1091$ as they were recorded monthly from July of the year 1925 to the may
of the year 2017.

(A)  (B)  (C)

(D)  (E)

Figure 10:  We run SMA-based change point estimation with star
kernel,bandwidth equal to 50 and the false alarm rate equal to 0.2. The
procedure provides 5 change points at moment March of 1934, April of 1955,
April of 1964, December of 1984, September of 1988.

In the original paper [31] the parametric Bayesian inference was used too
estimate nine and seven change points, and the total number of change point
was only guessed and not inferred. Our results differ significantly from this
original results, however different time range was used.

(A)

(B)

(C)

Figure 11:  By running SMA-procedure with symmetric kernel, bandwidth equal to
40 and the false alarm rate of 0.2, we get less conservative estimate. Change
points are February 1927, October 1940, July 1954, April 1964, January 1969,
May 1976, May 1990, July 2012.

(D)  (E)

It can be projected that all above change points can be attributed to the
important events in the economic history. For example, change point at July
1954 can be related to the end of recession of 1953, which itself can be
explained by the change of interrelations of the industrial sectors leading to
the change in the structure of the observed distribution. Hence, we can
describe performance of the SMA-procedure as satisfactory.

## Conclusion and Discussion

As result of the work a collection of methods of change point detection
methods was developed. All this methods are based on interaction of vector
ranks and geometric discrepancy which is a novel result. Certain basic
consistency result were proved for the method in its basic form designed for
detecting a single change point. However, they also can be applied for
detecting and estimating multiple change points.

Computational results shows applicability of the method both for simple
artificial change point problems an problems concerning actual data from
applications. It is strictly indicated by resulting experience that the method
is much more potent in situations then the distributions are not concentric.
Empirical evidences of our theoretical findings were also observed. These
theoretical results include expression of relation on different sides of
change point through inner product in the Hilbert space $\mathcal{D}$. Another
Theoretical result concerns dependence of estimate’s bias on ratio in which
change point separates sample.

Another positive results is the discovery of $\mathcal{D}$, which can be used
for establishing alternative iterative procedures defined by relations in this
Hilbert space. Furthermore, this Hilbert space structure on empirical measures
can be used for proving more theoretical results in the future.

Better proofs which cover ergodic processes and provide exact rates of
convergence a still need to be worked out for the methods. Furthermore
concentration results for change point estimates might render this algorithms
interesting for practical application. However, they are absent in the current
work.

Finally, online version of method can be implemented. The only requirement for
this algorithm is fast online computation of the optimal assignment problem.
It can be projected that such algorithm can be derived from the Sinkhorn
distance regularised algorithm which were designed by Cuturi [7].

## References

  * [1] Ilya A Antonov and VM Saleev  “An economic method of computing LP$\tau$-sequences”  In _USSR Computational Mathematics and Mathematical Physics_ 19.1 Elsevier, 1979, pp. 252–256
  * [2] Aleksandr Alekseevich Borovkov  “Asymptotically optimal solutions in the change-point problem”  In _Theory of Probability & Its Applications_ 43.4 SIAM, 1999, pp. 539–561
  * [3] E Brodsky and BS Darkhovsky  “Nonparametric Methods in Change Point Problems”  Springer Science & Business Media, 1993
  * [4] Victor Chernozhukov, Alfred Galichon, Marc Hallin and Marc Henry  “Monge-Kantorovich depth, quantiles, ranks, and signs”, 2014  arXiv:1412.8434 [stat.ST]
  * [5] Christine Choirat and Raffaello Seri  “The asymptotic distribution of quadratic discrepancies”  In _Monte Carlo and quasi-Monte Carlo methods 2004_ Springer, 2006, pp. 61–76
  * [6] Miklós CsörgHo and Lajos Horvath  “20 Nonparametric methods for changepoint problems”  In _Handbook of statistics_ 7 Elsevier, 1988, pp. 403–425
  * [7] Marco Cuturi  “Sinkhorn distances: Lightspeed computation of optimal transport”  In _Advances in Neural Information Processing Systems_ , 2013, pp. 2292–2300  arXiv:1306.0895 [stat.ML]
  * [8] Alexis Decurninge  “Multivariate quantiles and multivariate L-moments”, 2014  arXiv:1409.6013 [math.ST]
  * [9] A Dvoretzky, Jack Kiefer and Jacob Wolfowitz  “Sequential decision problems for processes with continuous time parameter. Testing hypotheses”  In _The Annals of Mathematical Statistics_ 24.2 Institute of Mathematical Statistics, 1953, pp. 254–264
  * [10] Meyer A Girshick and Herman Rubin  “A Bayes approach to a quality control model”  In _The Annals of mathematical statistics_ JSTOR, 1952, pp. 114–125
  * [11] Fred Hickernell  “A generalized discrepancy and quadrature error bound”  In _Mathematics of Computation of the American Mathematical Society_ 67.221, 1998, pp. 299–322
  * [12] Fred J Hickernell  “Error Analysis for Quasi-Monte Carlo Methods”, 2017  arXiv:1702.01487 [math.NA]
  * [13] Edmund Hlawka  “Discrepancy and Riemann integration”  In _Studies in Pure Mathematics (New York)_ 3, 1971
  * [14] L M Kantorovich  “On the transfer of masses”  In _Dokl. Akad. Nauk. SSSR_ 37, 1942, pp. 227–229
  * [15] Yoshinobu Kawahara and Masashi Sugiyama  “Sequential change-point detection based on direct density-ratio estimation”  In _Statistical Analysis and Data Mining_ 5.2 Wiley Online Library, 2012, pp. 114–127
  * [16] Albert Y Kim, Caren Marzban, Donald B Percival and Werner Stuetzle  “Using labeled data to evaluate change detectors in a multivariate streaming environment”  In _Signal Processing_ 89.12 Elsevier, 2009, pp. 2529–2536
  * [17] Harold W Kuhn  “The Hungarian method for the assignment problem”  In _Naval research logistics quarterly_ 2.1-2 Wiley Online Library, 1955, pp. 83–97
  * [18] Jia-Juan Liang, Kai-Tai Fang, Fred Hickernell and Runze Li  “Testing multivariate uniformity and its applications”  In _Mathematics of Computation_ 70.233, 2001, pp. 337–355
  * [19] Bruce G Lindsay  “Mixture models: theory, geometry and applications”  In _NSF-CBMS regional conference series in probability and statistics_ , 1995, pp. i–163  JSTOR
  * [20] Alexandre Lung-Yut-Fong, Céline Lévy-Leduc and Olivier Cappé  “Homogeneity and change-point detection tests for multivariate data using rank statistics”, 2011  arXiv:1107.1971 [math.ST]
  * [21] David S Matteson and Nicholas A James  “A nonparametric approach for multiple change point analysis of multivariate data”  In _Journal of the American Statistical Association_ 109.505 Taylor & Francis, 2014, pp. 334–345  arXiv:1306.4933 [stat.ME]
  * [22] Gaspard Monge  “Mémoire sur la théorie des déblais et des remblais”  In _Histoire de l‘Académie Royale des Sciences de Paris, avec les Mémoires de Mathématique et de Physique pour la méme année_ De l’Imprimerie Royale, 1781, pp. 666–704
  * [23] Vito MR Muggeo and Giada Adelfio  “Efficient change point detection for genomic sequences of continuous measurements”  In _Bioinformatics_ 27.2 Oxford University Press, 2010, pp. 161–166
  * [24] Harald Niederreiter  “Quasi-Monte Carlo methods and pseudo-random numbers”  In _Bulletin of the American Mathematical Society_ 84.6, 1978, pp. 957–1041
  * [25] Moshe Pollak  “Optimal detection of a change in distribution”  In _The Annals of Statistics_ JSTOR, 1985, pp. 206–227
  * [26] Filippo Santambrogio  “Optimal transport for applied mathematicians”  In _Birkäuser, NY (due in September 2015)_ Springer, 2015
  * [27] Robert Serfling  “Quantile functions for multivariate analysis: approaches and applications”  In _Statistica Neerlandica_ 56.2 Wiley Online Library, 2002, pp. 214–232
  * [28] WA Shewhart  “The application of statistics as an aid in maintaining quality of a manufactured product”  In _Journal of the American Statistical Association_ 20.152 Taylor & Francis Group, 1925, pp. 546–548
  * [29] Albert N Shiryaev  “On optimum methods in quickest detection problems”  In _Theory of Probability & Its Applications_ 8.1 SIAM, 1963, pp. 22–46
  * [30] I Sobol and Yu L Levitan  “The production of points uniformly distributed in a multidimensional cube”  In _Preprint IPM Akad. Nauk SSSR_ 40.3, 1976
  * [31] Makram Talih and Nicolas Hengartner  “Structural learning with time-varying components: tracking the cross-section of financial time series”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 67.3 Wiley Online Library, 2005, pp. 321–341
  * [32] N Tomizawa  “On some techniques useful for solution of transportation network problems”  In _Networks_ 1.2 Wiley Online Library, 1971, pp. 173–194
  * [33] Cédric Villani  “Optimal transport: old and new”  Springer Science & Business Media, 2008
  * [34] Hermann Weyl  “Über die gleichverteilung von zahlen mod. eins”  In _Mathematische Annalen_ 77.3 Springer, 1916, pp. 313–352