Nonparametric testing via partial sorting Nonparametric testing via partial
sorting K. Bisewski111University of Lausanne, Lausanne, Switzerland. and H.M.
Jansen222University College Roosevelt, Middelburg, the Netherlands. and Y.
Nazarathy333The University of Queensland, St. Lucia, Queensland, Australia.
Nonparametric statistics; goodness-of-fit; Brownian bridge; partial sorting
62G10; 62G30 62G20; 60F99 X Y 0 0 0 10.1214/YY-TN In this paper we introduce
the idea of partially sorting data to design nonparametric tests. This
approach gives rise to tests that are sensitive to both the order and the
underlying distribution of the data. We focus in particular on a test that
uses the bubble sort algorithm to partially sort the data. We show that a
function of the data, referred to as the empirical bubble sort curve,
converges uniformly to a limiting curve. We define a goodness-of-fit test
based on the distance between the empirical curve and its limit. The
asymptotic distribution of the test statistic is a generalization of the
Kolmogorov distribution. We apply the test in several examples and observe
that it outperforms classical nonparametric tests for appropriately chosen
sorting levels.

##  1 Introduction

Determining whether a data set constitutes a random sample from a given
distribution is one of the central problems in nonparametric statistics. A
common procedure is the standard Kolmogorov–Smirnov (KS) test for real-valued
observations $X_{1},\ldots,X_{n}$. The KS test considers the hypotheses

$H_{0}\colon X_{i}\underset{\text{iid}}{\sim}F_{0},\qquad\text{vs.}\qquad
H_{1}\colon\text{otherwise}.$ (1)

Here $F_{0}$ is some distribution specified as part of $H_{0}$. The procedure
computes the KS test statistic

$\widehat{D}_{n}=\sqrt{n}\sup_{x\in{\mathbb{R}}}\big{|}\widehat{F}_{n}(x)-F_{0}(x)\big{|}\qquad\text{with}\qquad\widehat{F}_{n}(x)=\frac{1}{n}\sum_{i=1}^{n}\mathds{1}{\mathopen{}\mathclose{{}\left\\{X_{i}\leq
x}\right\\}},$

where $F_{0}$ is the cdf and $\widehat{F}_{n}$ is the empirical cumulative
distribution function (ecdf).

The KS test does not prescribe or require any sorting of the observations.
Even if $X_{1},\ldots,X_{n}$ are fully sorted or partially sorted, this does
not affect the test statistic $\widehat{D}_{n}$. In other words, the KS test
is insensitive to the order of the data. In this paper we generalize KS by
introducing and studying a modified procedure that is based on partially
sorting $X_{1},\ldots,X_{n}$. The advantage of using partial sorting is that
it makes the modified procedure sensitive to the order as well as the
underlying distribution of the data. We study the properties of the modified
procedure and provide examples in which KS is outperformed by our
generalization. Before describing the main idea in more detail, let us review
some aspects of KS.

For any continuous $F_{0}$ under $H_{0}$, as $n$ grows, the sequence
$\widehat{D}_{n}$ converges in distribution to a random variable $D$ that
follows the Kolmogorov distribution with cdf

${\mathbb{P}}(D\leq x)=1-2\sum_{k=1}^{\infty}(-1)^{k-1}e^{-2k^{2}x^{2}}.$ (2)

The quantiles of the Kolmogorov distribution are easily computed numerically,
so the $p$-value or the test decision for a pre-specified confidence level
$\alpha$ can be easily determined. The distribution of $D$ is derived based on
a Brownian bridge approximation for the process
$\sqrt{n}\big{(}\widehat{F}_{n}(x)-F_{0}(x)\big{)}$. See for example [27, Sec.
19.3] for details or the initial publications [20, 24, 25]. An overview of
additional recent work related to KS tests is provided below.

Figure 1: iid normal data ($n=1,000$): Increasing levels of partial sorting
with no sorted data (top left) to fully sorted data (bottom right). The scaled
normal cdf is plotted in red.

Importantly, KS tests are not only used for testing the distribution $F_{0}$
but also for the independence stated under $H_{0}$ or lack thereof. As an
extreme example where this independence certainly does not hold, assume that
$X_{1}=X_{2}=\ldots=X_{n}$, with $X_{1}$ distributed according to $F_{0}$. In
this case

$\widehat{D}_{n}=\sqrt{n}\max\big{\\{}F_{0}(X_{1}),1-F_{0}(X_{1})\big{\\}},$

which is distributed uniformly in the range $[\sqrt{n}/2,\sqrt{n}]$. Now with
a specified type I error of $\alpha$, the null hypothesis $H_{0}$ is rejected
almost with certainty as $n$ grows. For example, with $\alpha=0.01$, the
$0.99$th quantile of the Kolmogorov distribution is approximately $1.68$. Then
$H_{0}$ is already rejected with certainty for $n=12$, since
$\sqrt{12}/2\approx 1.73$. This extreme example indicates that in general KS
statistics may be used to test for independence. Similarly, even if the lack
of independence is not as extreme as having all observations identical, a KS
test may often reject $H_{0}$ if the observations exhibit a form of
dependence. It is thus common to use KS tests for testing lack of
independence.

Another important statistical test, focusing specifically on independence, is
the Wald–Wolfowitz (WW) test [28]. It is based on the runs statistic, which
counts the number of times a data point above the median is followed by a data
point below the median or vice versa. While both KS and WW may be used to test
for lack of independence, these two tests focus on two different aspects of
the data. The KS test takes the distribution of the data into account but is
oblivious to its order. This means that it is sensitive to changes in the
values of the data but not to changes in the order of the data. In contrast,
the WW test relies on a statistic that takes the order of the data into
account but is oblivious to the distribution. This means that it is sensitive
to changes in the order of the data but not to changes in the values of the
data, especially if data points remain on the same side of the median.

Ideally, we would like to use a statistic that is sensitive to both aspects.
For this we present a new nonparametric paradigm based on partial sorting. We
believe that the ideas that we present here may lead to a new suite of tests
that may detect in certain cases a lack of independence better than current
methods. We develop and analyze one such specific test in this paper. We also
present examples of cases where partial sorting based tests are valuable.

The overarching idea of partial sorting based tests is to apply a sorting
algorithm to the data without running the algorithm to completeness. This
leads to a partially sorted sample, say $X_{1}^{\beta},\ldots,X_{n}^{\beta}$,
where $\beta\in(0,1]$ indicates the level of partial sorting and is referred
to as the sorting level. In the extreme case of $\beta=1$ the sample is fully
sorted, while for lower values of $\beta$ the sample is less sorted. A test
statistic is then computed based on the partially sorted sample
$X_{1}^{\beta},\ldots,X_{n}^{\beta}$. For the test statistic that we present
here, there is a known asymptotic distribution under $H_{0}$ that is not
influenced by $F_{0}$. Our results then yield a statistical test similar in
nature to KS, but based on the partially sorted sample and often more
sensitive to a lack of independence.

Figure 2: iid uniform data ($\beta=0.25$): The empirical bubble sort curve is
the right frontier of the points. It converges to a deterministic curve as
$n\to\infty$.

With $\beta<1$ the partial sorting approach halts the sorting algorithm in mid
flight. Thus the specific sorting algorithm at hand plays a key role in the
analysis: different sorting algorithms generate different probability laws for
$X_{1}^{\beta},\ldots,X_{n}^{\beta}$. In this paper our sole focus is on
partial sorting with the classic and simple bubble sort algorithm, which is a
strong contender for being the most intuitive sorting algorithm. See for
example [19, Section 5.2.2] or [2] for an historical overview. As its name may
suggest, bubble sort works by ‘bubbling up’ observations based on pairwise
comparisons of adjacent values. In each case where adjacent observations are
not in order they are swapped. Each such pass on the data is called an
iteration and the algorithm applies $n$ iterations to yield a sorted sample.
With bubble sort, a natural interpretation of $\beta$ is that the number of
iterations performed is the nearest integer to $\beta n$.

Our bubble sort based procedure generalizes the KS test, since it exactly
agrees with KS if $\beta=1$. Importantly, we show examples where using
$\beta<1$ outperforms both KS and WW. Our purpose with this paper is to show
viability of this method together with a complete derivation of the properties
of associated stochastic processes needed to obtain the asymptotic
distribution of the test statistic. In doing so, we derive a generalization of
the Kolmogorov distribution which we call the generalized Kolmogorov
distribution with random variable denoted $D^{\beta}$. Tabulation of the cdf
and/or quantiles of $D^{\beta}$ is simply a matter of computation similar to
the Kolmogorov distribution, and like the Kolmogorov distribution our
generalized Kolmogorov distribution is not influenced by $F_{0}$. Further, the
underlying stochastic process that yields the derivation of $D^{\beta}$ is
based on a generalization of the Brownian bridge process used for KS. Our
analysis includes a Glivenko–Cantelli like (strong law of large numbers)
theorem for this process, which is followed by a second order analysis
establishing convergence of probability measures.

To get a feel for potential applicability of partial sorting with bubble sort
consider Figure 1 which is based on an iid sequence of $n=1,000$ observations
from a standard normal distribution. The figure presents the original data,
followed by a display of the partially sorted data for $\beta=0.25$,
$\beta=0.5$, and finally $\beta=1$ (fully sorted). The nature of bubble sort
is that with sorting level $\beta$ the highest valued $\beta n$ (rounded to an
integer) observations are perfectly sorted in the respective top positions of
the array but the remaining $(1-\beta)n$ observations are only partially
sorted. Nevertheless as is evident from the figure a pattern emerges.
Specifically observe the right frontier in the case of $\beta=0.25$ and
$\beta=0.5$ and note that it follows a regular shape different from the cdf of
the observations. We use this regular pattern to define the bubble sort curve
in the next section and make use of it throughout the paper. In fact, we are
able to characterize the bubble sort curve in terms of $F_{0}$ and $\beta$.

Our statistical method compares the bubble sort curve expected under $H_{0}$
with an empirical bubble sort curve obtained from the right frontier of the
data. This operates in a manner that resembles the KS test. Remarkably, like
in the case of the KS test we are able to characterize the deviations between
the empirical bubble sort curve and the theoretical curve under $H_{0}$.
Specifically, as $n\to\infty$, the scaled deviations between the curves
converge to a well defined Gaussian process. We use this process to develop a
test statistic with a distribution under $H_{0}$ which is agnostic to $F_{0}$
and thus directly generalizes the Kolmogorov distribution associated with the
KS test.

To get a feel for the large sample asymptotics consider Figure 2 which is
based on standard iid uniform observations, each partially sorted with
$\beta=0.25$. The figure hints that as $n\to\infty$ the empirical bubble sort
curve (the right frontier of the data) converges. In Theorem 3.1 we establish
the associated first order convergence and further in Theorem 3.2 we study the
fluctuations around this frontier and characterize a limiting stochastic
process of these fluctuations. Finally in Theorem 3.4 we compute the
distribution of a statistic arising from the fluctuations, similar to the KS
statistic.

Similarly to the analysis that surrounds the KS test, our analysis is based
both on a law of large numbers type theorem and a second order theorem dealing
with convergence of probability measures. The former is an analogue of the
celebrated Glivenko–Cantelli theorem, see [27, Ch. 19] and [12]. The latter is
analogous to Donsker’s theorem for cdfs; see [7]. Hence our work is an
extension of the basic development of the asymptotics of the KS test as in
[7], [8], and [9]. Specifically we develop process level convergence that is
similar in spirit to the process level convergence which appears in the
analysis of the KS test.

A good historical account of Kolmogorov’s original paper [20] can be found in
[26]. Since that initial work, there have been dozens of extensions and
generalizations, different in nature to our work here. We now mention a few
notable publications in this space. In [21] the authors consider Rényi-type
statistics (see [23]) and a adapt them to create modified KS tests that are
sensitive to tail alternatives. In [22] the authors continue investigation of
fitting tails of the distribution and study exact properties of the so-called
$M_{n}$ statistics of Berk and Jones [4]. Related is the recent work [10]
dealing with two-sample KS type tests, using local levels. In [16] the authors
consider goodness of fit tests via $\phi$-divergences. Beyond these works,
many other papers have explored variations and adaptations of KS tests; see
the notable papers [6, 11, 13, 14, 15, 17, 18]. Nevertheless, to the best of
our knowledge, the use of partial sorting to enhance goodness of fit
procedures, as we present here, is new.

The remainder of this paper is structured as follows. In Section 2 we define
the bubble sort algorithm and the goodness of fit test that uses the partially
sorted sample. This also includes the definition of the generalized Kolmogorov
distribution. In Section 3 the main theoretical results of this paper our
presented. Then sections 4, 5, and 6 contain the derivations and proofs of the
results. In Section 4 we establish the law of large numbers type convergence
results. In Section 5 we establish the weak convergence results. In Section 6
these results are assembled to obtain the distribution of $D^{\beta}$. In
Section 7 we present numerical examples of our proposed procedure and we
conclude in Section 8.

##  2 The Test Procedure

Figure 3:  Left: Densities of the generalized Kolmogorov distribution for
several sorting levels (the case $\beta=1$ is the Kolmogorov distribution).
Right: Criticial value quantiles as a function of sorting levels.

The bubble sort test considers the null hypothesis that the data consist of
independent observations from a given continuous distribution $F_{0}$. The
procedure starts by partially sorting the data and computing the ecdf of the
running maximum of the partially sorted data. We call this ecdf the empirical
bubble sort curve. Then, similarly to the KS test, insight into the stochastic
behavior of the empirical bubble sort curve under $H_{0}$ allows us to compare
it to the bubble sort curve associated with $F_{0}$ and $\beta$. The distance
between these two curves yields the bubble sort statistic. Finally, we compare
the value of this statistic to the quantiles of its theoretical distribution
under $H_{0}$, which is a generalized Kolmogorov distribution that only
depends on $\beta$ and not on the underlying distribution $F_{0}$. See Figure
3 where we present densities for a few selected $\beta$ values as well as
common quantiles as a function of $\beta$. Note that at $\beta=1$ the
distribution is the Kolmogorov distribution (2) and as is evident from the
quantile plots, for $\beta$ in the approximate range of $0.7$–$1.0$ the
generalized Kolmogorov distributions are very close to the Kolmogorov
distribution. This is due to the fact that with bubble sort and such high
$\beta$ values while the sample isn’t guaranteed to be fully sorted it almost
always is.

Steps 1–5 below describe the bubble sort test procedure for a data sample
$x_{1},\ldots,x_{n}$. Key objects in this description include the partially
sorted sample (3), the empirical bubble sort curve (4), the limiting bubble
sort curve (5), the bubble sort statistic (6), and the generalized Kolmogorov
distribution (7). This distribution is defined via probabilities of the form
(8) associated with Brownian bridges.

  * 1

Partial sorting: Take $k$ as the nearest integer to $\beta n$ and apply $k$
bubble sort iterations to the data. The result is the partially sorted sample
$x_{1}^{\beta},\dotsc,x_{n}^{\beta}$. A bubble sort iteration on a data array
$v$ with elements $v[1],\ldots,v[n]$ is defined via the operator $\mathsf{T}$,
which is the composition $\mathsf{S}_{n-1}\circ\dotsb\circ\mathsf{S}_{1}$ with

$\displaystyle\mathsf{S}_{i}(v)=\begin{cases}\big{(}v[1],\dotsc,v[i-1],v[i+1],v[i],v[i+2],\dotsc,v[n]\big{)}&\text{
if }v[i]>v[i+1],\\\ \big{(}v[1],\dotsc,v[n]\big{)}&\text{ else. }\end{cases}$

The partially sorted sample $x_{1}^{\beta},\dotsc,x_{n}^{\beta}$ is then given
by

$\displaystyle(x_{1}^{\beta},\ldots,x_{n}^{\beta})=\mathsf{T}^{k}(x_{1},\ldots,x_{n}).$
(3)

  * 2

Empirical bubble sort curve: Obtain the empirical bubble sort curve
$\widehat{B}_{n}^{\beta}$ via

$\displaystyle\widehat{B}_{n}^{\beta}(x)=\frac{1}{n}\sum_{i=1}^{n}\mathds{1}{\mathopen{}\mathclose{{}\left\\{\max\\{x_{1}^{\beta},\dotsc,x_{i}^{\beta}\\}\leq
x}\right\\}}.$ (4)

This is the ecdf of the running maximum of the partially sorted sample.

  * 3

Bubble sort curve: Determine the bubble sort curve $B_{0}^{\beta}$ associated
with the distribution $F_{0}$ and sorting level $\beta$ via

$\displaystyle
B_{0}^{\beta}(x)=\min\Big{\\{}\frac{\beta}{1-F_{0}(x)},1\Big{\\}}-\min\big{\\{}\beta,1-F_{0}(x)\big{\\}}.$
(5)

  * 4

Bubble sort statistic: Compute the bubble sort statistic
$\widehat{D}_{n}^{\beta}$ via

$\displaystyle\widehat{D}_{n}^{\beta}=\sup_{x\in\mathbb{R}}\sqrt{n}\big{|}\widehat{B}_{n}^{\beta}(x)-B_{0}^{\beta}(x)\big{|}.$
(6)

  * 5

Test conclusion: Compare the value of the bubble sort statistic
$\widehat{D}_{n}^{\beta}$ to the quantiles of its asymptotic distribution
$D^{\beta}$ under $H_{0}$. This is a generalized Kolmogorov distribution with
cdf

$\displaystyle\mathbb{P}(D^{\beta}\leq
x)=2\int_{0}^{\sqrt{\tfrac{\beta}{1-\beta}}x}\Psi\Big{(}x;\tfrac{1-\beta}{\beta},\sqrt{\tfrac{1-\beta}{\beta}}z\Big{)}\Psi\Big{(}x;\beta,\sqrt{\beta(1-\beta)}z\Big{)}\phi(z)\,\mathrm{d}z,$
(7)

where

$\displaystyle\Psi(x;T,a)=\begin{cases}\sum_{k\in\mathbb{Z}}(-1)^{k}\exp\Big{(}-\frac{2kx(kx-a)}{T}\Big{)}&\text{
if }x>|a|,\\\ 0&\text{ if }x\leq|a|,\end{cases}$ (8)

and $\phi$ is the standard normal density function. The resulting asymptotic
$p$-value is $\mathbb{P}(D^{\beta}>\widehat{D}_{n}^{\beta})$ provided $F_{0}$
is continuous.

In the boundary case $\beta=1$, this procedure reduces to the KS test.
Although the expression in (7) is not defined for $\beta=1$, it is readily
verified that, as $\beta\nearrow 1$, the test statistic
$\widehat{D}_{n}^{\beta}$ converges weakly to the Kolmogorov distribution with
cdf (2). Tabulation of the cdf and quantiles for generalized Kolmogorov
distributions of the form (7) is computationally straightforward. See also the
accompanying GitHub repository for this paper [5].

Figure 4 illustrates the procedure via $n=100$ iid standard uniform random
variables with sorting level $\beta=0.25$. The associated empirical bubble
sort curve and the limiting bubble sort curve are plotted in the left display.
The scaled difference between the curves is plotted in the right display and
is used to compute the test statistic. In this example the value of the test
statistic is $1.598$. Based on the generalized Kolomogorov distribution for
$\beta=0.25$ (see Figure 3), the corresponding asymptotic $p$-value is
$0.701$. Thus $H_{0}$ is not rejected. Further examples are in Section 7.

Figure 4: An illustration of the test procedure with $\beta=0.25$. Left: The
empirical bubble sort curve (green) together with the bubble sort curve
(blue). Right: The scaled difference between the curves with the maximal
absolute value marking the test statistic $\widehat{D}_{n}^{\beta}$.

##  3 Main Results

We now present the main results supporting the validity of the bubble sort
test procedure. Our context in this section and the three sections that follow
is probabilistic. We therefore avoid any specific references to the test
procedure and instead focus on the probabilistic results valid under $H_{0}$.
In the remainder we fix a sorting level $\beta$ in $(0,1]$ and a distribution
$F_{0}$ for the iid random variables $X_{1},\dotsc,X_{n}$. Some of the results
below hinge on assuming $F_{0}$ is continuous, and we state this explicitly
when needed.

The key object we investigate is the bubble sort curve (5), which can also be
more explicitly represented via

$\displaystyle
B_{0}^{\beta}(x)=\begin{cases}\beta\frac{F_{0}(x)}{1-F_{0}(x)}&\text{ if
}x<x^{*}_{\beta},\\\ F_{0}(x)&\text{ else,}\end{cases}\qquad\text{where}\qquad
x^{*}_{\beta}=F_{0}^{-1}(1-\beta).$ (9)

Observe that since for $x<x^{*}_{\beta}$ the bubble sort curve is the scaled
odds of the probability $F(x)$, we can recover the cdf $F_{0}(\cdot)$ from
$B_{0}(\cdot)$ using $F_{0}(x)=B_{0}^{\beta}(x)/(\beta+B_{0}^{\beta}(x))$ for
such $x$. Thus, analogous to the cdf, the bubble sort curve is a
nondecreasing, right-continuous function that fully characterizes its
underlying distribution.

The empirical bubble sort curve, $\widehat{B}_{n}^{\beta}(\cdot)$ defined in
(4), plays the same role with respect to $B_{0}^{\beta}(\cdot)$ as the ecdf
plays with respect to the cdf. Here it is constructed as a random process with
index $x\in{\mathbb{R}}$ based on the random sample $X_{1},\ldots,X_{n}$,

$\widehat{B}_{n}^{\beta}(x)=\frac{1}{n}\sum_{i=1}^{n}\mathds{1}{\mathopen{}\mathclose{{}\left\\{\max\\{X_{1}^{\beta},\dotsc,X_{i}^{\beta}\\}\leq
x}\right\\}},\quad\text{where}\quad
X_{1}^{\beta},\ldots,X_{n}^{\beta}=\mathsf{T}^{k}(X_{1},\ldots,X_{n}),$

and with $k$ the nearest integer to $\beta n$ and $\mathsf{T}$ defined in (3).

The following theorem is analogous to the Glivenko–Cantelli theorem. Note that
we recover the Glivenko–Cantelli theorem by taking $\beta=1$. Section 4 is
devoted to the proof.

######  Theorem 3.1.

The uniform distance
$\sup_{-\infty<x<\infty}|\widehat{B}^{\beta}_{n}(x)-B^{\beta}_{0}(x)|\to 0$
almost surely as $n\to\infty$.

Now, analogous to the Donsker theorem underpinning the Kolmogorov–Smirnov
theorem, we wish to characterize the fluctuations of $\widehat{B}_{n}^{\beta}$
around its limiting curve $B^{\beta}_{0}$. We do this in the following theorem
which coincides with the Donsker theorem for empirical cumulative distribution
functions and Brownian bridges if $\beta=1$, yet generalizes it with
$\beta<1$. Section 5 is devoted to the proof.

######  Theorem 3.2.

If $F_{0}$ is continuous then the process
$\sqrt{n}(\widehat{B}_{n}^{\beta}-B^{\beta}_{0})$ converges weakly to
$Y^{\beta}_{0}$ in the M1 topology where $Y^{\beta}_{0}$ is right-continuous,
zero mean, Gaussian process with covariance function

$\displaystyle
c(x,y)=\operatorname*{\mathbb{C}ov}\big{(}Y_{0}^{\beta}(x),Y_{0}^{\beta}(y)\big{)}=\begin{cases}\beta\frac{F_{0}(x)}{\big{(}1-F_{0}(x)\big{)}^{2}}&\text{
if }x\leq y<x^{*}_{\beta},\\\\[5.0pt]
\beta\frac{F_{0}(x)\big{(}1-F_{0}(y)\big{)}}{\big{(}1-F_{0}(x)\big{)}^{2}}&\text{
if }x<x^{*}_{\beta}\leq y,\\\\[5.0pt] F_{0}(x)\big{(}1-F_{0}(y)\big{)}&\text{
if }x^{*}_{\beta}\leq x\leq y.\end{cases}$ (10)

Although the process $Y^{\beta}_{0}$ might seem unfamiliar at first, it is
closely related to generalized Brownian bridges, which we discuss below. For
this also consider Figure 5. The covariance function of $Y^{\beta}_{0}$ and
the realizations in Figure 5 reveal the presence of a jump at $x^{*}_{\beta}$.
To the left of $x^{*}_{\beta}$, the process $Y^{\beta}_{0}$ behaves like a
time-scaled Brownian motion. To the right of $x^{*}_{\beta}$, the process
$Y^{\beta}_{0}$ behaves like a time-scaled Brownian bridge. The crucial
property of $Y^{\beta}_{0}$ is that, conditional on the jump at
$x^{*}_{\beta}$, on each side of $x^{*}_{\beta}$ it behaves like an
independent generalized Brownian bridge under a time scaling.

We say that a stochastic process $\\{W(t):t\in[0,T]\\}$ follows the law of a
generalized Brownian bridge from $(0,0)$ to $(T,a)$ and write $W\sim{\rm
BB}(T,a)$ if $W$ is a continuous Gaussian process with mean function
$\mathbb{E}\big{(}W(t)\big{)}=\frac{ta}{T}$ and covariance function
$\operatorname*{\mathbb{C}ov}\big{(}W(s),W(t)\big{)}=\frac{s(T-t)}{T}$ for
$0\leq s\leq t\leq T$. Note that a standard Brownian motion on $[0,T]$
conditioned on reaching $a$ at time $T$ follows the law of ${\rm BB}(T,a)$. In
particular, ${\rm BB}(1,0)$ is the law of a standard Brownian bridge on
$[0,1]$. An important probability associated with general Brownian bridges is

$\Psi(x;T,a)=\mathbb{P}\big{(}\sup_{t\in[0,T]}|W(t)|\leq x\big{)},$ (11)

which can be represented as the the double infinite series expression (8). See
for example [3, Eq. (4.12)] for a derivation. Note that $\Psi(x;1,0)$ recovers
the Kolmogorov distribution (2).

Figure 5: Three random outcomes of
$\sqrt{n}(\widehat{B}_{n}^{\beta}(x)-B^{\beta}_{0}(x))$ with $\beta=0.25$ and
$n=10,000$, where all elements of the initial array are iid, uniformly
distributed over $[-1,1]$. The vertical dashed line is at $x^{*}_{\beta}=0.5$.

With generalized Brownian bridges at hand we now construct a generalized
Kolmogorov distribution. For this let $Z$ be a random variable with a standard
normal distribution. Conditionally on $Z$, we let the processes $W_{1}$ and
$W_{2}$ be independent with laws

$W_{1}\sim{\rm
BB}\mathopen{}\mathclose{{}\left(\tfrac{1-\beta}{\beta},\sqrt{\tfrac{1-\beta}{\beta}}Z}\right)\quad\text{and}\quad
W_{2}\sim{\rm
BB}\mathopen{}\mathclose{{}\left(\beta,\sqrt{\beta(1-\beta)}Z}\right).$ (12)

We define the distribution of

$D^{\beta}=\sup_{t\in[0,(1-\beta)/\beta]}|W_{1}(t)|~{}~{}\vee~{}~{}\sup_{t\in[0,\beta]}|W_{2}(t)|$
(13)

as a generalized Kolmogorov distribution depending on the parameter $\beta$.
Using this construction we have an expression for the distribution of
$D^{\beta}$ in (7). The proof is in Section 6.

######  Lemma 3.3.

The cdf of $D^{\beta}$ as defined in (13) is given in (7).

Our final result deals with the asymptotic distribution of the bubble sort
statistic (6) and hinges on the process level convergence of Theorem 3.2.
Section 6 is devoted to the proof.

######  Theorem 3.4.

If $F_{0}$ is continuous then the bubble sort statistic
$\widehat{D}^{\beta}_{n}$ of (6) converges in distribution to $D^{\beta}$ of
(13).

We emphasize that the generalized Kolmogorov distribution depends only on the
sorting level $\beta$ and is insensitive to the underlying distribution
$F_{0}$. Similarly, the distribution of $\widehat{D}^{\beta}_{n}$ for any
finite $n$ also does not depend on the underlying distribution of $X_{i}$. As
an illustration of the convergence in distribution of Theorem 3.4 consider
Figure 6 focusing on the bubble sort statistic with three sorting levels
$\beta$. The black curves in all three plots are identical and are the cdf of
$D^{\beta}$ computed via (7). In all three cases we consider sample sizes of
$n=20$, $n=200$, and $n=2,000$. In each of these cases we simulate $30,000$
repetitions of the samples to obtain the bubble sort statistic (6) and plot
its estimated (finite $n$) distribution. It is visually apparent that as $n$
grows, the distributions converge to the generalized Kolmogorov distribution.
As is further apparent, for lower sorting levels $\beta$ the convergence is
slower.

Figure 6: Convergence in distribution of the bubble sort statistic
$\widehat{D}^{\beta}_{n}$ to the generalized Kolmogorov distribution
$D^{\beta}$ for various sorting levels. Note that the distributions of
$\widehat{D}^{\beta}_{n}$ are not influenced by the distribution of $X_{i}$.

##  4 The Law of Large Numbers

This section is devoted to the proof of Theorem 3.1. More precisely, we
consider the random array $V_{n}$ and study the frontier
$\widehat{B}_{n}^{\beta}$ associated with the partially sorted array
$V_{n}^{\beta n}$. Our goal is to show that the frontier
$\widehat{B}_{n}^{\beta}$ converges uniformly to its limiting bubble sort
curve $B_{0}^{\beta}$ defined in (9). In the remainder of this section we omit
the superscript $\beta$ and the subscript $0$, and with this notation our goal
is to show that $\widehat{B}_{n}$ converges uniformly to $B$.

Before we prove this result, we explore how bubble sort operates on binary
arrays, which contain only 0’s and 1’s. The next lemma explains how bubble
sort changes the positions of the 1’s. We demonstrate later that convergence
properties of general random arrays are connected to convergence properties of
related binary arrays.

######  Lemma 4.1.

Consider a binary array $U_{n}$. Assume that exactly $m$ of its entries equal
$1$ and let $I_{n,1}<\dotsc<I_{n,m}$ be the positions of the 1’s in $U_{n}$.
Denote by $I_{n,1}^{k}<\dotsc<I_{n,m}^{k}$ the positions of the 1’s in the
partially sorted array $U_{n}^{k}$. Then

$\displaystyle I_{n,i}^{k}$ $\displaystyle=\begin{cases}I_{n,i+k}-k&\text{ if
}\quad k\leq m-i,\\\ n-(m-i)&\text{ if }\quad
k>m-i,\end{cases}\qquad\text{for}\qquad i=1,\ldots,m.$

######  Proof 4.2.

We start by focusing on the 1 at position $I_{i}$ with $i<m$ in the array
$U_{n}$ during the first bubble sort iteration. Clearly, the 1 at $I_{i}$ does
not get swapped with any values to the left of position $I_{i}$, because all
values in $U_{n}$ are at most 1. We now distinguish two cases. If the value to
the right of $I_{i}$ equals 1, then $I_{i+1}=I_{i}+1$. In this case, the 1 at
position $I_{i}$ does not get swapped with any value, so
$I_{i}^{1}=I_{i}=I_{i+1}-1$. If the value to the right of $I_{i}$ does not
equal 1, then $I_{i+1}>I_{i}+1$. In this case, the 1 at position $I_{i}$ gets
swapped with all 0’s between position $I_{i}$ and $I_{i+1}$, so it reaches the
position on the left of $I_{i+1}$. However, it does not swap with the 1 at
position $I_{i+1}$, so $I_{i}^{1}=I_{i+1}-1$ in this case too. Therefore, the
statement of the lemma holds for $k=1$. This may be extended to all $k\leq
m-i$ by induction.

We now turn to the case $k>m-i$. The previous arguments also show that one
bubble sort iteration moves every 1 to the position of the last 0 in the
sequence of 0’s immediately right to it. That means that $I_{n,m}^{k}=n$ for
$k>0$ and $I_{n,i}^{k}=n-(m-i)$ for $k>m-i$ by induction.

Lemma 4.1 enables us to relate the position $I_{1}$ of the leftmost 1 in
$U_{n}$ to the position $I_{1}^{k}$ in $U_{n}^{k}$, which we obtain after
applying $k$ bubble sort iterations to $U_{n}$. This is important, because
understanding the dynamics of the leftmost 1 in $U_{n}$ is crucial for
understanding the behavior of its frontier.

As mentioned before, our goal is to connect the convergence of a general
random array to convergence of a related binary random array. More concretely,
given the random array $V_{n}$, we fix $x\in\mathbb{R}$ such that $0<F(x)<1$
and define the new binary array $U_{n}$ via

$\displaystyle U_{n}[i]=\begin{cases}0&\text{ if }\quad V_{n}[i]\leq x,\\\
1&\text{ if }\quad V_{n}[i]>x.\end{cases}$

The idea is that the position of the leftmost value in $V_{n}^{k}$ exceeding
$x$ has the same position as the leftmost $1$ in $U_{n}^{k}$, regardless of
the number $k$ of bubble sort iterations that we apply to both arrays. We
formalize this in the next lemma and use it in combination with Lemma 4.1 to
establish the convergence of $\widehat{B}_{n}$ to $B$.

######  Lemma 4.3.

For each $k=0,1,2,\ldots,n$, it holds that

$\displaystyle U_{n}^{k}[i]=\begin{cases}0&\text{ if }\quad V_{n}^{k}[i]\leq
x,\\\ 1&\text{ if }\quad V_{n}^{k}[i]>x,\end{cases}$

so in particular

$\displaystyle\inf\\{i=1,\dotsc,n~{}|~{}V_{n}^{k}[i]>x\\}$
$\displaystyle=\inf\\{i=1,\dotsc,n~{}|~{}U_{n}^{k}[i]=1\\}.$

######  Proof 4.4.

The statement of the lemma is trivially true if $V_{n}$ does not contain any
value exceeding $x$ or if $k=0$, so we assume that this is not the case. We
only show that the statement holds for $k=1$, since the case $k>1$ follows by
induction.

To establish the claim for $k=1$, it suffices to show that applying an
arbitrary swap operator $\mathsf{S}_{j}$ to both $V_{n}$ and $U_{n}$ yields
$\mathsf{S}_{j}(V_{n})[i]>x$ if and only if $\mathsf{S}_{j}(U_{n})[i]=1$.
Since by construction $V_{n}[i]>x$ if and only if $U_{n}[i]=1$, this can only
be possibly violated if $U_{n}[i]=1$, $U_{n}[i+1]=0$ and $j=i$. However, this
means that $V_{n}[i]>x$ and $V_{n}[i+1]\leq x$, so
$\mathsf{S}_{i}(V_{n})[i]=V_{n}[i+1]\leq x$ and
$\mathsf{S}_{i}(V_{n})[i+1]=V_{n}[i]>x$ with
$\mathsf{S}_{i}(U_{n})[i]=U_{n}[i+1]=0$ and
$\mathsf{S}_{i}(U_{n})[i+1]=U_{n}[i]=1$ in this case. This proves the claim.

We can use the previous two results to establish the convergence of
$\widehat{B}_{n}(x)$ to $B(x)$ as follows. We recall that $\beta$ is the
sorting level and that the entries of $V_{n}$ consist of the first $n$
elements of the sequence of iid random variables $X_{1},X_{2},\dotsc$. Thus
the entries of $U_{n}$ consist of the first $n$ elements of the sequence of
iid random variables
$\mathds{1}{\mathopen{}\mathclose{{}\left\\{X_{1}>x}\right\\}},\mathds{1}{\mathopen{}\mathclose{{}\left\\{X_{2}>x}\right\\}},\dotsc$.
The latter each follows a Bernoulli distribution with success probability
$p=1-F(x)\in(0,1)$.

Since
$\mathds{1}{\mathopen{}\mathclose{{}\left\\{X_{1}>x}\right\\}},\mathds{1}{\mathopen{}\mathclose{{}\left\\{X_{2}>x}\right\\}},\dotsc$
are iid, the index of the first 1 in that sequence has a geometric
distribution with success probability $p$. The same is true for the distance
between two consecutive 1’s in that sequence. Thus the positions of the 1’s in
the sequence
$\mathds{1}{\mathopen{}\mathclose{{}\left\\{X_{1}>x}\right\\}},\mathds{1}{\mathopen{}\mathclose{{}\left\\{X_{2}>x}\right\\}},\dotsc$
give rise to a sequence of iid random variables $G_{1},G_{2},\dotsc$ that have
a geometric distribution with parameter $p$.

The total number of 1’s in the binary array $U_{n}$ is given by

$\displaystyle
H_{n}=\sum_{i=1}^{n}\mathds{1}{\mathopen{}\mathclose{{}\left\\{X_{i}>x}\right\\}}$

and has a binomial distribution with parameter $p$. In line with the notation
used in Lemma 4.1, we also define

$\displaystyle I_{k}$ $\displaystyle=\sum_{i=1}^{k}G_{i}$

for $k=1,2,\dotsc$. The random variable $I_{k}$ simply denotes the position of
the $k$-th 1 in the sequence
$\mathds{1}{\mathopen{}\mathclose{{}\left\\{X_{1}>x}\right\\}},\mathds{1}{\mathopen{}\mathclose{{}\left\\{X_{2}>x}\right\\}},\dotsc$.
This should be distinguished from the random variable $I_{n,k}$, which denotes
the position of the $k$-th 1 in $U_{n}$.

Applying the results of Lemma 4.1 and Lemma 4.3 combined with our definition
of $I_{k}$, we see that the position $I_{n,1}^{\beta n}$ of the first 1 in the
partially sorted array $U_{n}^{\beta n}$ satisfies

$\displaystyle I_{n,1}^{\beta n}=\begin{cases}I_{1+\beta n}-\beta n&\text{ if
}\quad\beta n\leq H_{n}-1,\\\ n-H_{n}+1&\text{ if }\quad\beta n>H_{n}-1,\\\
\end{cases}$

so

$\displaystyle I_{n,1}^{\beta n}=\min\\{I_{1+\beta n},n\\}-\min\\{\beta
n,H_{n}-1\\}.$

Indeed, if $I_{1+\beta n}\leq n$, then there are at least $1+\beta n$ entries
of $U_{n}$ that contain a 1. Then the number of bubble sort iterations $\beta
n$ that we apply to $U_{n}$ is smaller than the total number of 1’s in
$U_{n}$, so Lemma 4.1 implies that $I_{n,1}^{\beta n}=I_{1+\beta n}-\beta n$
in this case. If $I_{1+\beta n}>n$, then $U_{n}$ has at most $\beta n$ entries
that contain a 1, so $H_{n}\leq\beta n$ and $I_{n,1}^{\beta n}=n-(H_{n}-1)$ by
Lemma 4.1 in this case.

We now combine these results to establish the convergence of
$\widehat{B}_{n}(x)$ to $B(x)$. As defined in (4), the frontier
$\widehat{B}_{n}$ is given as the empirical cdf of the running maximum of
$V_{n}^{\beta n}$. Thus $\widehat{B}_{n}(x)$ is equal to the relative number
of entries in the running maximum that do not exceed $x$. Clearly, the first
value that does exceed $x$ in the running maximum has the same position as the
first value exceeding $x$ in $V_{n}^{\beta n}$. Additionally, due to Lemma
4.3, the first value exceeding $x$ in $V_{n}^{\beta n}$ has the same position
as the first 1 in $U_{n}^{\beta n}$. Then $I_{n,1}^{\beta n}$ is also the
position of first value exceeding $x$ in the running maximum. Consequently,
$\widehat{B}_{n}(x)$ satisfies

$\displaystyle\widehat{B}_{n}(x)=\frac{1}{n}(I_{1,n}^{\beta
n}-1)=\min\Big{\\{}\frac{1}{n}I_{1+\beta
n},1\Big{\\}}-\min\Big{\\{}\frac{1+\beta n}{n},\frac{1}{n}H_{n}\Big{\\}}.$
(14)

Since $\frac{1}{n}I_{1+\beta n}\to\beta/p$ and $\frac{1}{n}H_{n}\to p$ almost
surely, the continuous-mapping theorem implies that

$\displaystyle\widehat{B}_{n}(x)\to\min\\{\beta/p,1\\}-\min\\{\beta,p\\}$

almost surely. Taking into account that $\beta<p$ if and only if $\beta/p<1$
and that $p=1-F(x)$, we conclude that $\widehat{B}_{n}(x)$ converges to $B(x)$
almost surely. This can be extended to almost sure uniform convergence of
$\widehat{B}_{n}$ to $B$ by the standard Glivenko–Cantelli arguments (cf. [27,
Th. 19.1]).

##  5 Weak Convergence

This section is devoted to the proof of Theorem 3.2. The proof has a classical
structure: we first establish weak convergence of the finite-dimensional
distributions (fdds) of the scaled and centered frontier
$\widetilde{B}_{n}=\sqrt{n}(\widehat{B}_{n}-B)$ on a dense subset of the real
line. It turns out that the fdds of $\widetilde{B}_{n}$ converge to the fdds
of a process that behaves like a time-scaled Brownian motion on
$(-\infty,x^{*}_{\beta})$ and like a time-scaled Brownian bridge on
$(x^{*}_{\beta},\infty)$, with a jump occuring at $x^{*}_{\beta}$. We then
proceed to show that $\widetilde{B}_{n}$ is tight in the M1-topology. These
two steps combined imply the weak convergence result of Theorem 3.2.

###  5.1 The finite-dimensional distributions

We start by establishing weak convergence of the fdds of $\widetilde{B}_{n}$.
For ease of exposition, we restrict ourselves to weak convergence of two-
dimensional distributions of the form
$(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$, where $x_{1}<x_{2}$ and
$x_{1},x_{2}\not=x^{*}_{\beta}$. We separately study three different cases:
Case I with $x^{*}_{\beta}<x_{1}<x_{2}$, Case II with
$x_{1}<x^{*}_{\beta}<x_{2}$, and Case III with $x_{1}<x_{2}<x^{*}_{\beta}$.

####  5.1.1 Notation

We introduce some new notation to prove the desired convergence results. In
line with the notation of the previous section, we define $I_{k}(x)$ as the
position of the $k$-th value exceeding $x$ in the underlying sequence of
random variables $X_{1},X_{2},\dotsc$. We also define
$H_{k}(x)=\sum_{i=1}^{k}\mathds{1}{\mathopen{}\mathclose{{}\left\\{X_{i}>x}\right\\}}$
as the total number of values that exceed $x$ in $X_{1},\dotsc,X_{k}$.

We let $G_{1}(x_{1}),G_{2}(x_{1}),\dotsc$ denote a collection of independent
random variables that follow a geometric distribution with parameter
$p_{1}=1-F(x_{1})$. The variable $G_{1}(x_{1})$ is the position of the first
value exceeding $x_{1}$ in the sequence $X_{1},X_{2},\dotsc$, while
$G_{1}(x_{1})+G_{2}(x_{1})$ is the position of the second value exceeding
$x_{1}$ in that sequence and so on.

Each value in $X_{1},X_{2},\dotsc$ that exceeds $x_{1}$ also exceeds $x_{2}$
with probability $p_{2}=(1-F(x_{2}))/(1-F(x_{1}))$. This gives rise to another
collection $G_{1}(x_{2}),G_{2}(x_{2}),\dotsc$ of independent random variables
having a geometric distribution with parameter $p_{2}$. Thus $G_{1}(x_{2})$ is
the position of the first value exceeding $x_{2}$ in the sequence of values
exceeding $x_{1}$, while $G_{1}(x_{2})+G_{2}(x_{2})$ is the position of the
second value exceeding $x_{2}$ in the sequence of values exceeding $x_{1}$ and
so on.

Given these definitions, we let $S_{k}(x)=\sum_{i=1}^{k}G_{i}(x)$. With this
notation, the position of the $k$-th value exceeding $x_{1}$ in
$X_{1},X_{2},\dotsc$ is given by

$\displaystyle I_{k}(x_{1})=\sum_{i=1}^{k}G_{i}(x_{1})=S_{k}(x_{1})$ (15)

and the position of the $k$-th value exceeding $x_{2}$ in $X_{1},X_{2},\dotsc$
is given by

$\displaystyle
I_{k}(x_{2})=\sum_{i=1}^{\sum_{j=1}^{k}G_{j}(x_{2})}G_{i}(x_{1})=S_{S_{k}(x_{2})}(x_{1}).$
(16)

The partial sum $S_{k}(x)$ gives rise to an associated counting process
$C_{t}(x)$ defined via $C_{t}(x)=\max\\{k\geq 0\,|\,S_{k}(x)\leq t\\}$, where
we take $t\in[0,\infty)$ and interpret $S_{0}(x)=0$. The random variable
$C_{k}(x_{1})$ simply counts the number of values exceeding $x_{1}$ in
$X_{1},\dotsc,X_{k}$, while $C_{k}(x_{2})$ counts how many of the first $k$
values that exceed $x_{1}$ also exceed $x_{2}$. Therefore we may represent

$\displaystyle H_{k}(x_{1})=C_{k}(x_{1})$ (17)

and

$\displaystyle H_{k}(x_{2})=C_{C_{k}(x_{1})}(x_{2}).$ (18)

Thus we can represent each of the variables $I_{n}(x_{1})$, $I_{n}(x_{2})$,
$H_{n}(x_{1})$, and $H_{n}(x_{2})$ in terms of nested partial sums or counting
processes. The convergence properties of these nested versions of partial sums
and their associated counting processes form the crux of the proof.

####  5.1.2 Proof of weak convergence

We are now in a position to establish the desired weak convergence results. We
start with convergence of the counting processes and their associated partial
sums, which we exploit to prove weak convergence of
$(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$.

With $t\in[0,\infty)$ and $[t]$ being equal to $t$ rounded down to the nearest
integer, we may regard $C_{t}(x_{i})$ and $S_{[t]}(x_{i})$ as right-continuous
stochastic processes on $[0,\infty)$. In view of [29, Th. 7.3.2], their scaled
and centered versions

$\displaystyle\frac{1}{\sqrt{n}}\big{(}C_{nt}(x_{i})-p_{i}nt\big{)}\quad\text{
and }\quad\frac{1}{\sqrt{n}}\big{(}S_{[nt]}(x_{i})-\frac{1}{p_{i}}nt\big{)}$

converge jointly to the processes

$\displaystyle-\sqrt{1-p_{i}}W_{p_{i}t}^{(i)}\quad\text{ and
}\quad\frac{\sqrt{1-p_{i}}}{p_{i}}W_{t}^{(i)},$

where $W_{t}^{(i)}$ is a standard Brownian motion with $W_{t}^{(1)}$ and
$W_{t}^{(2)}$ independent. Combined with the new representations of
$I_{n}(x_{1})$, $I_{n}(x_{2})$, $H_{n}(x_{1})$, and $H_{n}(x_{2})$ above, we
can use this result to derive weak convergence of
$(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$ for all three cases.

_Case I._ If $x^{*}_{\beta}<x_{1}<x_{2}$, then (14) implies that
$(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$ is asymptotically
equivalent to

$\displaystyle-\bigg{(}\sqrt{n}\bigg{(}\frac{1}{n}H_{n}(x_{1})-(1-F(x_{1})\bigg{)},\sqrt{n}\bigg{(}\frac{1}{n}H_{n}(x_{2})-(1-F(x_{2}))\bigg{)}\bigg{)}.$

The first element of this random vector equals
$\frac{1}{\sqrt{n}}(C_{n}(x_{1})-p_{1}n)$, while the second equals
$\frac{1}{\sqrt{n}}\big{(}C_{C_{n}(x_{1})}(x_{2})-p_{1}p_{2}\big{)}$. We can
rewrite the last expression as

$\displaystyle\frac{1}{\sqrt{n}}\big{(}C_{n\bar{C}_{n}(x_{1})}(x_{2})-p_{2}n\bar{C}_{n}(x_{1})\big{)}+p_{2}\frac{1}{\sqrt{n}}\big{(}C_{n}(x_{1})-p_{1}n\big{)},$
(19)

where we define $\bar{C}_{n}(x_{1})=\frac{1}{n}C_{n}(x_{1})$. Based on these
considerations, we conclude that
$(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$ is asymptotically
equivalent to

$\displaystyle-\bigg{(}\frac{1}{\sqrt{n}}(C_{n}(x_{1})-p_{1}n),p_{2}\frac{1}{\sqrt{n}}\big{(}C_{n}(x_{1})-p_{1}n\big{)}+\frac{1}{\sqrt{n}}\big{(}C_{n\bar{C}_{n}(x_{1})}(x_{2})-p_{2}n\bar{C}_{n}(x_{1})\big{)}\bigg{)}.$

Given the convergence of counting processes and their partial sums as well as
the convergence in probability of $\bar{C}_{n}(x_{1})$ to $p_{1}$, we see that
this vector converges weakly to

$\displaystyle\begin{bmatrix}1&0\\\
p_{2}&1\end{bmatrix}\begin{bmatrix}-\sqrt{1-p_{1}}W_{p_{1}}^{(1)}\\\
-\sqrt{1-p_{2}}W_{p_{2}p_{1}}^{(2)}\end{bmatrix}=-\begin{bmatrix}\sqrt{1-p_{1}}&0\\\
\sqrt{1-p_{1}}p_{2}&\sqrt{1-p_{2}}\end{bmatrix}\begin{bmatrix}W_{p_{1}}^{(1)}\\\
W_{p_{1}p_{2}}^{(2)}\end{bmatrix}.$

Since $W^{(1)}$ and $W^{(2)}$ are independent standard Brownian motions, this
implies that $(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$ converges
weakly to a zero-mean normal distribution with covariance matrix

$\displaystyle\begin{bmatrix}\sqrt{1-p_{1}}&0\\\
\sqrt{1-p_{1}}p_{2}&\sqrt{1-p_{2}}\end{bmatrix}\begin{bmatrix}p_{1}&0\\\
0&p_{1}p_{2}\end{bmatrix}\begin{bmatrix}\sqrt{1-p_{1}}&\sqrt{1-p_{1}}p_{2}\\\
0&\sqrt{1-p_{2}}\end{bmatrix}.$

Bearing in mind that $p_{1}=1-F(x_{1})$ and $p_{1}p_{2}=1-F(x_{2})$, some
straightforward computations reveal that the last matrix equals

$\displaystyle\begin{bmatrix}p_{1}(1-p_{1})&p_{1}(1-p_{1})p_{2}\\\
p_{1}(1-p_{1})p_{2}&p_{1}p_{2}(1-p_{1}p_{2})\end{bmatrix}=\begin{bmatrix}F(x_{1})(1-F(x_{1}))&F(x_{1})(1-F(x_{2}))\\\
F(x_{1})(1-F(x_{2}))&F(x_{2})(1-F(x_{2}))\end{bmatrix}.$

The entries of the last matrix are indeed characterized by the function $c$
specified in (10). Thus $(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$
converges weakly to a zero-mean normal distribution whose covariance matrix is
specified by (10). Not surprisingly, this is the limiting distribution that
one obtains in a Kolmogorov–Smirnov setting with full sorting.

_Case II._ If $x_{1}<x^{*}_{\beta}<x_{2}$, then (14) implies that
$(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$ is asymptotically
equivalent to

$\displaystyle\bigg{(}\sqrt{n}\bigg{(}\frac{1}{n}I_{\beta
n}(x_{1})-\beta-\beta\frac{F(x_{1})}{1-F(x_{1})}\bigg{)},-\sqrt{n}\bigg{(}\frac{1}{n}H_{n}(x_{2})-(1-F(x_{2}))\bigg{)}\bigg{)}.$

The first element of this random vector equals
$\frac{1}{\sqrt{n}}(S_{[n\beta]}(x_{1})-\frac{1}{p_{1}}n\beta)$, while the
second element can be rewritten as (19). Thus
$(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$ is asymptotically
equivalent to

$\displaystyle\bigg{(}\frac{1}{\sqrt{n}}(S_{[n\beta]}(x_{1})-\frac{1}{p_{1}}n\beta),-p_{2}\frac{1}{\sqrt{n}}\big{(}C_{n}(x_{1})-p_{1}n\big{)}+\frac{1}{\sqrt{n}}\big{(}C_{n\bar{C}_{n}(x_{1})}(x_{2})-p_{2}n\bar{C}_{n}(x_{1})\big{)}\bigg{)}.$

Again relying on the convergence of counting processes and their partial sums
as well as the convergence in probability of $\bar{C}_{n}(x_{1})$ to $p_{1}$,
we see that this vector converges weakly to

$\displaystyle\begin{bmatrix}1&0&0\\\
0&-p_{2}&1\end{bmatrix}\begin{bmatrix}\frac{\sqrt{1-p_{1}}}{p_{1}}W_{\beta}^{(1)}\\\
-\sqrt{1-p_{1}}W_{p_{1}}^{(1)}\\\
-\sqrt{1-p_{2}}W_{p_{2}p_{1}}^{(2)}\end{bmatrix}=\begin{bmatrix}\sqrt{1-p_{1}}{p_{1}}&0&0\\\
0&\sqrt{1-p_{1}}p_{2}&-\sqrt{1-p_{2}}\end{bmatrix}\begin{bmatrix}W_{\beta}^{(1)}\\\
W_{p_{1}}^{(1)}\\\ W_{p_{2}p_{1}}^{(2)}\end{bmatrix}.$

Since $W^{(1)}$ and $W^{(2)}$ are independent standard Brownian motions and
$\beta<p_{1}$, this implies that
$(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$ converges weakly to a
zero-mean normal distribution with covariance matrix

$\displaystyle\begin{bmatrix}\sqrt{1-p_{1}}{p_{1}}&0&0\\\
0&\sqrt{1-p_{1}}p_{2}&-\sqrt{1-p_{2}}\end{bmatrix}\begin{bmatrix}\beta&\beta&0\\\
\beta&p_{1}&0\\\
0&0&p_{1}p_{2}\end{bmatrix}\begin{bmatrix}\sqrt{1-p_{1}}{p_{1}}&0\\\
0&\sqrt{1-p_{1}}p_{2}\\\ 0&-\sqrt{1-p_{2}}\end{bmatrix},$

which equals

$\displaystyle\begin{bmatrix}\beta\frac{1-p_{1}}{p_{1}^{2}}&\beta\frac{1-p_{1}}{p_{1}}p_{2}\\\
\beta\frac{1-p_{1}}{p_{1}}p_{2}&p_{1}p_{2}(1-p_{1}p_{2})\end{bmatrix}=\begin{bmatrix}\beta\frac{F(x_{1})}{(1-F(x_{1}))^{2}}&\beta\frac{F(x_{1})}{(1-F(x_{1}))^{2}}(1-F(x_{2}))\\\
\beta\frac{F(x_{1})}{(1-F(x_{1}))^{2}}(1-F(x_{2}))&F(x_{2})(1-F(x_{2}))\end{bmatrix}.$

We conclude that $(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$
converges weakly to a zero-mean normal distribution whose covariance matrix is
specified by (10).

_Case III._ If $x_{1}<x_{2}<x^{*}_{\beta}$, then (14) implies that
$(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$ is asymptotically
equivalent to

$\displaystyle\bigg{(}\sqrt{n}\bigg{(}\frac{1}{n}I_{\beta
n}(x_{1})-\beta-\beta\frac{F(x_{1})}{1-F(x_{1})}\bigg{)},\sqrt{n}\bigg{(}\frac{1}{n}I_{\beta
n}(x_{2})-\beta-\beta\frac{F(x_{2})}{1-F(x_{2})}\bigg{)}\bigg{)}.$

The first element of this random vector equals
$\frac{1}{\sqrt{n}}(S_{[n\beta]}(x_{1})-\frac{1}{p_{1}}n\beta)$, while the
second element equals
$\frac{1}{\sqrt{n}}(S_{S_{[n\beta]}(x_{2})}(x_{1})-\frac{1}{p_{1}p_{2}}n\beta)$.
We can rewrite the last expression as

$\displaystyle\frac{1}{\sqrt{n}}\big{(}S_{n\bar{S}_{[n\beta]}(x_{2})}(x_{1})-\frac{1}{p_{1}}n\bar{S}_{[n\beta]}(x_{2})\big{)}+\frac{1}{p_{1}}\frac{1}{\sqrt{n}}(S_{[n\beta]}(x_{2})-\frac{1}{p_{2}}n\beta),$

where we define $\bar{S}_{[n\beta]}(x_{2})=\frac{1}{n}S_{[n\beta]}(x_{2})$.
Thus $(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$ is asymptotically
equivalent to

$\displaystyle\bigg{(}\frac{1}{\sqrt{n}}(S_{[n\beta]}(x_{1})-\frac{1}{p_{1}}n\beta),\frac{1}{\sqrt{n}}\big{(}S_{n\bar{S}_{[n\beta]}(x_{2})}(x_{1})-\frac{1}{p_{1}}n\bar{S}_{[n\beta]}(x_{2})\big{)}+\frac{1}{p_{1}}\frac{1}{\sqrt{n}}(S_{[n\beta]}(x_{2})-\frac{1}{p_{2}}n\beta)\bigg{)}.$

Taking into account that $\bar{S}_{[n\beta]}(x_{2})$ converges in probability
to $\beta/p_{2}$, we see that this vector converges weakly to

$\displaystyle\begin{bmatrix}1&0&0\\\
0&1&\frac{1}{p_{1}}\end{bmatrix}\begin{bmatrix}\frac{\sqrt{1-p_{1}}}{p_{1}}W_{\beta}^{(1)}\\\
\frac{\sqrt{1-p_{1}}}{p_{1}}W_{\beta/p_{2}}^{(1)}\\\
\frac{\sqrt{1-p_{2}}}{p_{2}}W_{\beta}^{(2)}\end{bmatrix}=\begin{bmatrix}\frac{\sqrt{1-p_{1}}}{p_{1}}&0&0\\\
0&\frac{\sqrt{1-p_{1}}}{p_{1}}&\frac{\sqrt{1-p_{2}}}{p_{1}p_{2}}\end{bmatrix}\begin{bmatrix}W_{\beta}^{(1)}\\\
W_{\beta/p_{2}}^{(1)}\\\ W_{\beta}^{(2)}\end{bmatrix}.$

This implies that $(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$
converges weakly to a zero-mean normal distribution with covariance matrix

$\displaystyle\begin{bmatrix}\frac{\sqrt{1-p_{1}}}{p_{1}}&0&0\\\
0&\frac{\sqrt{1-p_{1}}}{p_{1}}&\frac{\sqrt{1-p_{2}}}{p_{1}p_{2}}\end{bmatrix}\begin{bmatrix}\beta&\beta&0\\\
\beta&\frac{\beta}{p_{2}}&0\\\
0&0&\beta\end{bmatrix}\begin{bmatrix}\frac{\sqrt{1-p_{1}}}{p_{1}}&0\\\
0&\frac{\sqrt{1-p_{1}}}{p_{1}}\\\
0&\frac{\sqrt{1-p_{2}}}{p_{1}p_{2}}\end{bmatrix},$

which equals

$\displaystyle\begin{bmatrix}\beta\frac{1-p_{1}}{p_{1}^{2}}&\beta\frac{1-p_{1}}{p_{1}^{2}}\\\
\beta\frac{1-p_{1}}{p_{1}^{2}}&\beta\frac{1-p_{1}p_{2}}{p_{1}^{2}p_{2}^{2}}\end{bmatrix}=\begin{bmatrix}\beta\frac{F(x_{1})}{(1-f(x_{1}))^{2}}&\beta\frac{F(x_{1})}{(1-f(x_{1}))^{2}}\\\
\beta\frac{F(x_{1})}{(1-f(x_{1}))^{2}}&\beta\frac{F(x_{2})}{(1-f(x_{2}))^{2}}\end{bmatrix}.$

We conclude that $(\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{2}))$
converges weakly to a zero-mean normal distribution whose covariance matrix is
specified by (10).

The previous arguments establish the convergence of the two-dimensional
distributions of $\widetilde{B}_{n}$ on a dense subset of the real line. The
computations already reveal that the limit is Gaussian and provide the
corresponding covariance matrix. Convergence of the higher-dimensional
distribution can be derived along similar lines.

###  5.2 Tightness

We now turn to tightness of the process $\widetilde{B}_{n}$ in the M1
topology. For completeness, we first present characterizations of tightness
and then use these to establish the M1 tightness of $\widetilde{B}_{n}$.
Further background for tightness and the M1 topology can be found in [29, Ch.
11–13].

A sequence of stochastic processes $Z_{n}$ on $\mathbb{R}$ is C-tight if for
every interval $[-\tau,\tau]$ the following two conditions are met. First, for
every $\epsilon>0$ there exists a positive constant $c$ such that
$\mathbb{P}(\mathopen{}\mathclose{{}\left|Z_{n}(-\tau)}\right|>c)<\epsilon$
for all $n$. Second, for every $\epsilon>0$ and $\eta>0$ there exist
$\delta>0$ and $N$ such that

$\displaystyle\mathbb{P}\Big{(}u_{\tau}(Z_{n},\delta)>\eta\Big{)}<\epsilon$

for all $n\geq N$. The oscillation function $u_{\tau}$ is defined by taking

$\displaystyle
u_{\tau}(z,x,\delta)=\sup\\{\mathopen{}\mathclose{{}\left|z(x_{1})-z(x_{2})}\right|\,:\,-\tau\vee(x-\delta)\leq
x_{1}<x_{2}\leq(x+\delta)\wedge\tau\\}$

and

$\displaystyle u_{\tau}(z,\delta)=\sup_{-\tau\leq
x\leq\tau}u_{\tau}(z,x,\delta).$

The process $Z_{n}$ converges weakly to a continuous stochastic process $Z$ if
and only if its fdds converge weakly to the fdds of $Z$ and $Z_{n}$ is
C-tight.

The process $\widetilde{B}_{n}$ clearly does not converge to a continuous
stochastic process: its limit exhibits a jump at the point $x^{*}_{\beta}$.
This is apparent in Figure 5 as well as in the form of the covariance
structure in (10). In view of this, the M1 topology is the natural choice for
establishing weak convergence of $\widetilde{B}_{n}$.

A sequence of stochastic processes $Z_{n}$ on $\mathbb{R}$ converges to $Z$ in
the M1 topology if its fdds converge weakly to the fdds of $Z$ on a dense
subset of $\mathbb{R}$ and $Z_{n}$ is tight in the M1 topology. The sequence
$Z_{n}$ is tight in the M1 topology if for every element $\tau$ in a positive,
unbounded sequence $\\{\tau_{k}\\}_{k=1}^{\infty}$ the following two
conditions are met. First, for every $\epsilon>0$ there exists a positive
constant $c$ such that
$\mathbb{P}(\sup\\{\mathopen{}\mathclose{{}\left|Z_{n}(x)}\right|\,:\,-\tau\leq
x\leq\tau\\}>c)<\epsilon$ for all $n$. Second, for every $\epsilon>0$ and
$\eta>0$ there exist $\delta>0$ and $N$ such that

$\displaystyle\mathbb{P}\Big{(}w_{\tau}(Z_{n},\delta)>\eta\Big{)}<\epsilon$

for all $n\geq N$. The oscillation function $w_{\tau}$ is defined as follows.
We interpret $[a,b]$ as the interval $[a\wedge b,a\vee b]$ and let
$\mathopen{}\mathclose{{}\left|c-[a,b]}\right|$ denote the distance between a
point $c$ and the interval $[a,b]$. We then take

$\displaystyle
v_{\tau}(z,x,\delta)=\sup\\{\mathopen{}\mathclose{{}\left|z(x_{2})-[z(x_{1}),z(x_{3})]}\right|\,:\,-\tau\vee(x-\delta)\leq
x_{1}<x_{2}<x_{3}\leq(x+\delta)\wedge\tau\\},$

and

$\displaystyle v_{\tau}(z,\delta)=\sup_{-\tau\leq
x\leq\tau}v_{\tau}(z,x,\delta),$

which are used to define the oscillation function

$\displaystyle w_{\tau}(z,\delta)=u_{\tau}(z,-\tau,\delta)\vee
v_{\tau}(z,\delta)\vee u_{\tau}(z,\tau,\delta).$

We use these characterizations of tightness to prove that $\widetilde{B}_{n}$
is tight in the M1 topology.

To establish tightness of $\widetilde{B}_{n}$ in the M1 topology, we fix
$\tau>1+|x^{*}_{\beta}|$ as well as $\epsilon>0$ and $\eta>0$. As a first step
we note that the processes

$\displaystyle\widetilde{I}_{n}(x)=\sqrt{n}\Big{(}\frac{1}{n}I_{1+\beta
n}(x)-\beta\frac{F(x)}{1-F(x)}\Big{)}$

and

$\displaystyle\widetilde{H}_{n}(x)=\sqrt{n}\Big{(}\frac{1}{n}H_{n}(x)-(1-F(x))\Big{)}$

are both C-tight as a consequence of [29, Th. 11.6.5]. This implies that there
exist $\delta>0$ and $N$ such that the next three statements are true. First,
it implies that $\mathbb{P}(\sup\\{|\widetilde{B}_{n}(x)|\,:\,-\tau\leq
x\leq\tau\\}>c)<\epsilon$ for a sufficiently large constant, since
$|\widetilde{B}_{n}(x)|\leq|\widetilde{I}_{n}(x)|+|\widetilde{H}_{n}(x)|$.
Second, it implies that
$\mathbb{P}(u_{\tau}(\widetilde{B}_{n},-\tau,\delta)>\eta)<\epsilon/3$ and
$\mathbb{P}(u_{\tau}(\widetilde{B}_{n},\tau,\delta)>\eta)<\epsilon/3$. Third,
it implies that

$\displaystyle\frac{1}{n}I_{1+\beta n}(x_{1})$
$\displaystyle=\frac{1}{n}I_{1+\beta
n}(x_{2})+\Big{(}\beta\frac{F(x_{1})}{1-F(x_{1})}-\beta\frac{F(x_{2})}{1-F(x_{2})}\Big{)}+o\Big{(}\frac{\eta}{\sqrt{n}}\Big{)}$
$\displaystyle=\frac{1}{n}I_{1+\beta
n}(x_{2})+\Big{(}\frac{\beta}{1-F(x_{1})}-\frac{\beta}{1-F(x_{2})}\Big{)}+o\Big{(}\frac{\eta}{\sqrt{n}}\Big{)}$

and

$\displaystyle\frac{1}{n}H_{n}(x_{1})=\frac{1}{n}H_{n}(x_{2})+(F(x_{1})-F(x_{2}))+o\Big{(}\frac{\eta}{\sqrt{n}}\Big{)}$

for all $\mathopen{}\mathclose{{}\left|x_{1}-x_{2}}\right|<\delta$ on a set
$\Omega^{*}$ with probability at least $1-\epsilon/3$ for all $n\geq N$. Here
each $o(c)$ indicates a term whose absolute value is at most $c>0$.

We now analyze four separate cases for $x_{1}<x_{2}<x_{3}$ on the set
$\Omega^{*}$. We show that in each case

$\displaystyle\widetilde{B}_{n}(x_{1})+o(2\eta)\leq\widetilde{B}_{n}(x_{2})\leq\widetilde{B}_{n}(x_{3})+o(2\eta)$
(20)

or

$\displaystyle\widetilde{B}_{n}(x_{1})+o(2\eta)\geq\widetilde{B}_{n}(x_{2})\geq\widetilde{B}_{n}(x_{3})+o(2\eta),$
(21)

implying that
$|\widetilde{B}_{n}(x_{2})-[\widetilde{B}_{n}(x_{1}),\widetilde{B}_{n}(x_{3})]|\leq
2\eta$. The M1 tightness of $\widetilde{B}_{n}$ is a direct consequence of
this.

Case I: $\frac{1}{n}H_{n}(x_{3})>\beta$. In this case
$\frac{1}{n}H_{n}(x_{1})>\beta$ and $\frac{1}{n}H_{n}(x_{2})>\beta$, so

$\displaystyle
B_{n}(x_{1})=B_{n}(x_{2})+\Big{(}\frac{\beta}{1-F(x_{1})}-\frac{\beta}{1-F(x_{2})}\Big{)}+o\Big{(}\frac{\eta}{\sqrt{n}}\Big{)}$

and

$\displaystyle
B_{n}(x_{2})=B_{n}(x_{3})+\Big{(}\frac{\beta}{1-F(x_{2})}-\frac{\beta}{1-F(x_{3})}\Big{)}+o\Big{(}\frac{\eta}{\sqrt{n}}\Big{)}.$

This yields

$\displaystyle\widetilde{B}_{n}(x_{1})=\widetilde{B}_{n}(x_{2})+\sqrt{n}\Big{(}\frac{\beta}{1-F(x_{1})}-\frac{\beta}{1-F(x_{2})}+B(x_{2})-B(x_{1})\Big{)}+o(\eta).$
(22)

and

$\displaystyle\widetilde{B}_{n}(x_{2})=\widetilde{B}_{n}(x_{3})+\sqrt{n}\Big{(}\frac{\beta}{1-F(x_{2})}-\frac{\beta}{1-F(x_{3})}+B(x_{3})-B(x_{2})\Big{)}+o(\eta).$

Now note that

$\displaystyle\frac{\beta}{1-F(x_{i})}-\frac{\beta}{1-F(x_{i+1})}+B(x_{i+1})-B(x_{i})$
$\displaystyle{}=\min\Big{\\{}0,F(x_{i+1})-\beta\frac{F(x_{i+1})}{1-F(x_{i+1})}\Big{\\}}-\min\Big{\\{}0,F(x_{i})-\beta\frac{F(x_{i})}{1-F(x_{i})}\Big{\\}}\leq
0,$

since $\min\\{0,F(x)-\beta\frac{F(x)}{1-F(x)}\\}$ is nonincreasing. Therefore
(20) holds in this case.

Case II: $\frac{1}{n}H_{n}(x_{3})\leq\beta$ and
$\frac{1}{n}H_{n}(x_{2})>\beta$. In this case $\frac{1}{n}H_{n}(x_{1})>\beta$,
so (22) is satisfied. Based on the equality

$\displaystyle 1+\beta-\frac{1}{n}I_{1+\beta
n}(x_{2})-\frac{1}{n}H_{n}(x_{3})\geq 1+\beta-\frac{1}{n}I_{1+\beta
n}(x_{2})-\frac{1}{n}H_{n}(x_{2})\geq 0$

as well as analogous arguments as in the previous case, we see that

$\displaystyle\widetilde{B}_{n}(x_{3})\geq\widetilde{B}_{n}(x_{2})$
$\displaystyle+\sqrt{n}\max\Big{\\{}F(x_{3})-(1-\beta),\frac{\beta}{1-F(x_{3})}-1\Big{\\}}$
$\displaystyle{}-\sqrt{n}\max\Big{\\{}F(x_{2})-(1-\beta),\frac{\beta}{1-F(x_{2})}-1\Big{\\}}+o(2\eta).$

Since $\max\\{F(x)-(1-\beta),\frac{\beta}{1-F(x)}-1\\}$ is nondecreasing, we
conclude that $\widetilde{B}_{n}(x_{3})+o(2\eta)\geq\widetilde{B}_{n}(x_{2})$.
Therefore (20) holds in this case.

Case III: $\frac{1}{n}H_{n}(x_{2})\leq\beta$ and
$\frac{1}{n}H_{n}(x_{1})>\beta$. In this case
$1+\beta\leq\frac{1}{n}I_{1+\beta n}(x_{2})+\frac{1}{n}H_{n}(x_{1})$, which
leads to

$\displaystyle\widetilde{B}_{n}(x_{1})\geq\widetilde{B}_{n}(x_{2})$
$\displaystyle-\sqrt{n}\max\Big{\\{}F(x_{3})-(1-\beta),\frac{\beta}{1-F(x_{3})}-1\Big{\\}}$
$\displaystyle{}+\sqrt{n}\max\Big{\\{}F(x_{2})-(1-\beta),\frac{\beta}{1-F(x_{2})}-1\Big{\\}}+o(2\eta).$

Since $\max\\{F(x)-(1-\beta),\frac{\beta}{1-F(x)}-1\\}$ is nondecreasing, it
follows that $\widetilde{B}_{n}(x_{1})+o(2\eta)\geq\widetilde{B}_{n}(x_{2})$.
We also note that both $\frac{1}{n}H_{n}(x_{2})\leq\beta$ and
$\frac{1}{n}H_{n}(x_{3})\leq\beta$, so

$\displaystyle\widetilde{B}_{n}(x_{2})$
$\displaystyle=\widetilde{B}_{n}(x_{3})+\sqrt{n}\Big{(}\min\Big{\\{}\beta\frac{F(x_{2})}{1-F(x_{2})},0\Big{\\}}-\min\Big{\\{}\beta\frac{F(x_{1})}{1-F(x_{1})},0\Big{\\}}\Big{)}+o(\eta).$
(23)

Since $\min\\{\beta\frac{F(x)}{1-F(x)},0\\}$ is nondecreasing, this yields
$\widetilde{B}_{n}(x_{2})\geq\widetilde{B}_{n}(x_{3})+o(\eta)$. Therefore (21)
holds in this case.

Case IV: $\frac{1}{n}H_{n}(x_{1})\leq\beta$. In this case
$\frac{1}{n}H_{n}(x_{3})\leq\frac{1}{n}H_{n}(x_{2})\leq\frac{1}{n}H_{n}(x_{1})\leq\beta$,
implying that (23) is valid and remains so if we replace $x_{2}$ by $x_{1}$
and $x_{3}$ by $x_{2}$. Therefore (21) holds in this case.

##  6 Analyzing the bubble sort statistic

In this section we prove Theorem 3.4 and Lemma 3.3. The result of Theorem 3.2
combined with the fact that weak convergence of a process in the M1 topology
implies weak convergence of its supremum (cf. [29, Th. 13.4.1]), implies that
$\widehat{D}^{\beta}_{n}=\sup_{x\in\mathbb{R}}\sqrt{n}(\widehat{B}^{\beta}_{n}(x)-B^{\beta}_{0}(x))$
converges in distribution to $\sup_{x\in\mathbb{R}}Y^{\beta}_{0}(x)$. It
remains to show that $\sup_{x\in\mathbb{R}}Y_{0}^{\beta}(x)$ equals
${D}^{\beta}$ in distribution (Theorem 3.4), and to derive an expression for
the distribution of ${D}^{\beta}$ (Lemma 3.3).

We prove Theorem 3.4 in two steps. First, conditional on the jump of
$Y_{0}^{\beta}$ in $x^{*}_{\beta}$, we relate the distribution of
$\sup_{x\in\mathbb{R}}Y_{0}^{\beta}(x)$ to the distribution of generalized
Brownian bridges. Based on this, we then argue that
$\sup_{x\in\mathbb{R}}Y_{0}^{\beta}(x)$ indeed follows a generalized
Kolmogorov distribution.

Consider a right-continuous, centred Gaussian process
$\\{\check{Y}^{\beta}(t):t\in[0,1]\\}$ with covariance function

$\displaystyle\operatorname*{\mathbb{C}ov}(\check{Y}^{\beta}(t_{1}),\check{Y}^{\beta}(t_{2}))=\begin{cases}\frac{\beta
t_{1}}{(1-t_{1})^{2}}&\text{ if }t_{1}\leq t_{2}<1-\beta,\\\ \frac{\beta
t_{1}(1-t_{2})}{(1-t_{1})^{2}}&\text{ if }t_{1}<1-\beta<t_{2},\\\
t_{1}(1-t_{2})&\text{ if }1-\beta<t_{1}\leq t_{2}.\end{cases}$

Clearly $Y^{\beta}_{0}(t)=\check{Y}^{\beta}(F_{0}(t))$ and
$\sup_{t}|\check{Y}^{\beta}(t)|=\sup_{x}|Y^{\beta}_{0}(x)|$, since $F_{0}$ is
continuous. We define the function $\omega(t)=\beta t/(1-t)^{2}$, which is
strictly increasing for $0<t<1$. Additionally, we define the two stochastic
processes $\\{W_{1}(t):t\in[0,\tfrac{1-\beta}{\beta}]\\}$ and
$\\{W_{2}(t):t\in[0,\beta]\\}$, where

$\displaystyle W_{1}(t)$ $\displaystyle=\check{Y}^{\beta}(\omega^{-1}(t)),\
t\in[0,\tfrac{1-\beta}{\beta})\quad\text{and}\quad
W_{1}(\tfrac{1-\beta}{\beta})=\lim_{t\to(1-\beta)^{-}}\check{Y}^{\beta}(\omega^{-1}(t)),$
$\displaystyle W_{2}(t)$ $\displaystyle=\check{Y}^{\beta}(1-t),\
t\in[0,\beta]\quad\text{and}\quad
W_{2}(\beta)=\lim_{t\to(1-\beta)^{+}}\check{Y}^{\beta}(t).$

These processes can be thought of as the (time-scaled and reversed) left and
right parts of $\check{Y}^{\beta}$, viewed relative to the time point
$t=1-\beta$. Similar to before, we have

$\sup_{t}|\check{Y}^{\beta}(t)|=\sup_{t}|W_{1}(t)|~{}\vee~{}\sup_{t}|W_{2}(t)|$

and

$\displaystyle\operatorname*{\mathbb{C}ov}(W_{1}(t_{1}),W_{1}(t_{2}))$
$\displaystyle=t_{1}$ $\displaystyle\text{for}\quad 0\leq t_{1}\leq
t_{2}\leq\tfrac{1-\beta}{\beta},$
$\displaystyle\operatorname*{\mathbb{C}ov}(W_{1}(t_{1}),W_{2}(t_{2}))$
$\displaystyle=t_{1}t_{2}$ $\displaystyle\text{for}\quad
t_{1}\in[0,\tfrac{1-\beta}{\beta}],t_{2}\in[0,\beta],$
$\displaystyle\operatorname*{\mathbb{C}ov}(W_{2}(t_{1}),W_{2}(t_{2}))$
$\displaystyle=t_{1}(1-t_{2})$ $\displaystyle\text{for}\quad 0\leq t_{1}\leq
t_{2}\leq\beta.$

We proceed to study the distributions of $W_{1}$ and $W_{2}$ conditional on
$W_{2}(\beta)$. Since $W_{1}$ and $W_{2}$ are jointly Gaussian, their
conditional distribution is also jointly Gaussian. To determine this
conditional distribution, it suffices to determine the conditional means,
covariances, and cross-correlations.

First, we focus on the cross-correlation. For any
$t_{1}\in[0,\tfrac{1-\beta}{\beta}]$ and $t_{2}\in[0,\beta]$ we have

$\displaystyle\operatorname*{\mathbb{C}ov}(W_{1}(t_{1}),W_{2}(t_{2})\mid
W_{2}(\beta))$
$\displaystyle=\operatorname*{\mathbb{C}ov}(W_{1}(t_{1}),W_{2}(t_{2}))-\frac{\operatorname*{\mathbb{C}ov}(W_{1}(t_{1}),W_{2}(\beta))\cdot\operatorname*{\mathbb{C}ov}(W_{2}(t_{2}),W_{2}(\beta))}{\operatorname*{\mathbb{V}ar}(W_{2}(\beta))}$
$\displaystyle=t_{1}t_{2}-\frac{t_{1}\beta\cdot
t_{2}(1-\beta)}{\beta(1-\beta)}$ $\displaystyle=0,$

implying that $W_{1}$ and $W_{2}$ are independent conditional on
$W_{2}(\beta)$.

Second, we study the conditional distribution of $W_{1}$. For any
$t\in[0,\tfrac{1-\beta}{\beta}]$ we have

$\displaystyle\mathbb{E}(W_{1}(t)\mid W_{2}(\beta))$
$\displaystyle=\mathbb{E}W_{1}(t)-\frac{\operatorname*{\mathbb{C}ov}(W_{1}(t_{1}),W_{2}(\beta))}{\operatorname*{\mathbb{V}ar}(W_{2}(\beta))}\cdot\mathopen{}\mathclose{{}\left(W_{2}(\beta)-\mathbb{E}W_{2}(\beta)}\right)$
$\displaystyle=\frac{t\beta}{\beta(1-\beta)}\cdot W_{2}(\beta)$
$\displaystyle=\frac{t}{\frac{1-\beta}{\beta}}\cdot\frac{W_{2}(\beta)}{\beta}.$

For any $0\leq t_{1}\leq t_{2}\leq\tfrac{1-\beta}{\beta}$ we have

$\displaystyle\operatorname*{\mathbb{C}ov}(W_{1}(t_{1}),W_{1}(t_{2})\mid
W_{2}(\beta))$
$\displaystyle=\operatorname*{\mathbb{C}ov}(W_{1}(t_{1}),W_{1}(t_{2}))-\frac{\operatorname*{\mathbb{C}ov}(W_{1}(t_{1}),W_{2}(\beta))\cdot\operatorname*{\mathbb{C}ov}(W_{1}(t_{2}),W_{2}(\beta))}{\operatorname*{\mathbb{V}ar}(W_{2}(\beta))}$
$\displaystyle=t_{1}-\frac{t_{1}\beta\cdot t_{2}\beta}{\beta(1-\beta)}$
$\displaystyle=\frac{t_{1}\cdot(\tfrac{1-\beta}{\beta}-t_{2})}{\frac{1-\beta}{\beta}}.$

Similarly we obtain the conditional mean and covariance of $W_{2}$. For any
$t\in[0,\beta]$ and $0\leq t_{1}\leq t_{2}\leq\beta$ we have

$\displaystyle\mathbb{E}(W_{2}(t)\mid W_{2}(\beta))=\frac{t}{\beta}\cdot
W_{2}(\beta)$

and

$\displaystyle\operatorname*{\mathbb{C}ov}(W_{2}(t_{1}),W_{1}(t_{2})\mid
W_{2}(\beta))=\frac{t_{1}(\beta-t_{2})}{\beta}.$

The previous arguments demonstrate that, conditional on $W_{2}(\beta)$, the
processes $W_{1}$ and $W_{2}$ are independent. Using the formulas for the mean
and covariance of a generalized Brownian bridge. We conclude that

$\big{(}W_{1}\mid W_{2}(\beta)\big{)}\sim{\rm
BB}(\tfrac{1-\beta}{\beta},\tfrac{W_{2}(\beta)}{\beta})\quad\text{and}\quad\big{(}W_{2}\mid
W_{2}(\beta)\big{)}\sim{\rm BB}(\beta,W_{2}(\beta)).$ (24)

Further, since $W_{2}(\beta)$ is a normally distributed zero mean random
variable with variance $\beta(1-\beta)$, it is equal in distribution to
$\sqrt{\beta(1-\beta)}Z$, where $Z$ is a standard normal random variable.
Hence the representation in (24) is equivalent to (12) which is used to
construct $D^{\beta}$ as in (13). This concludes the proof of Theorem 3.4.

We close this section with a derivation of the cdf of $D^{\beta}$ as defined
in (7).

######  Proof 6.1 (Proof of Lemma 3.3).

Based on the conditional independence between $W_{1}$ and $W_{2}$ as well as
the definition of $\Psi(x;T,a)$ of (11), we obtain

$\displaystyle\mathbb{P}\Big{(}D^{\beta}\leq x\Big{)}$
$\displaystyle=\mathbb{E}\Big{[}\mathbb{P}\big{(}\sup_{t\in[0,(1-\beta)/\beta]}|W_{1}(t)|\leq
x,\sup_{t\in[0,\beta]}|W_{2}(t)|\leq x\mid Z\big{)}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\mathbb{P}\big{(}\sup_{t\in[0,(1-\beta)/\beta]}|W_{1}(t)|\leq
x\mid Z\big{)}~{}\mathbb{P}\big{(}\sup_{t\in[0,\beta]}|W_{2}(t)|\leq x\mid
Z\big{)}\Big{]}$
$\displaystyle=\int_{-\infty}^{\infty}\Psi\big{(}x;\tfrac{1-\beta}{\beta},\sqrt{\tfrac{1-\beta}{\beta}}\cdot
z\big{)}\Psi\big{(}x;\beta,\sqrt{\beta(1-\beta)}\cdot z\big{)}\phi(z)\,{\rm
d}z.$

The integrand in the last line is non-zero if and only if
$x>\sqrt{(1-\beta)/\beta}\cdot|z|$ and $x>\sqrt{(1-\beta)\beta}\cdot|z|$.
Finally, using that $\Psi(x;T,a)=\Psi(x;T,-a)$ for all $a\in\mathbb{R}$, we
obtain the desired result.

##  7 Numerical Examples

We now present two scenarios where the bubble sort based procedure of Section
2 is useful with the purpose of showing the potential value of nonparametric
testing with partial sorting. The first scenario deals with detecting partial
sorting of tabular data, and the second deals with detecting service time
dependent scheduling policies in queues. Both of these examples illustrate
that there are cases where partial sorting based procedures can be used to
detect lack of independence in a better manner than KS or WW.

In fact, in any situation where one considers the hypotheses (1), it is
possible to use partial sorting ($\beta<1$) in place of the standard $\beta=1$
KS test, or the WW test. However, in our experimentation with standard cases
involving simulated data from auto-regressive sequences, or data from hidden
Markov sequences, we did not always observe a clear benefit in partial
sorting. This may be due to the fact that the number of observations falling
on the empirical bubble sort curve (4) is typically significantly smaller than
those falling on the ecdf of the data. Nevertheless, as the following two
examples show, there certainly exist cases where partial sorting is
beneficial.

In both examples below, we assume that the distribution $F_{0}$ of (1) is
known and is not the subject of testing. Instead our focus is on testing for
lack of independence which is detected when $H_{0}$ is rejected. Further, in
both examples we consider a grid of sorting levels $\beta\in(0,1]$. As we see,
the performance of the test is heavily influenced by the choice of $\beta$.
Policies to optimally choose $\beta$, or perhaps to carry out more involved
procedures based on multiple $\beta$ values, are not within the scope of the
current paper. Our purpose here is simply to illustrate the potential value of
nonparametric testing with partial sorting.

In our experimentation with both examples, data is generated via simulations
in scenarios where $H_{0}$ does not hold. We then quantify performance by
fixing $\alpha=0.1$ and estimating the probability of rejecting $H_{0}$ via
repeated simulation experiments. The results are then plotted for
$\beta\in(0,1]$, where we keep in mind that $\beta=1$ is the KS test. We also
summarize results from the WW test on the same data for comparison.

### Example 1: Detection of magnitude based sorting of a hidden column in
tabular data

Consider the common setting of numerical tabular data, e.g. a data frame or a
spreadsheet, where different variables constitute different columns and
observations are in rows. With such data it is very common to sort the
observations by a given variable/column. This then modifies the order of other
columns accordingly.

Figure 7: Example 1: Detection of magnitude based sorting of a hidden column
in tabular data. Results for $n=1,000$, $\alpha=0.1$ (black line), and several
levels of $\rho$. The proportion of cases detected ($H_{0}$ rejected) is
plotted as a function of the sorting level $\beta$.

Assume now that we are presented with $x_{1},\ldots,x_{n}$, the variables of a
single column from the data and we wish to test if the data was sorted based
on an unobserved column, say $\widetilde{x}_{1},\ldots,\widetilde{x}_{n}$.
Here correlations between columns may leave “traces” of the sorting of the
unobserved $\widetilde{x}_{1},\ldots,\widetilde{x}_{n}$ in the observed
$x_{1},\ldots,x_{n}$. The KS test will not detect such sorting since it is
oblivious to the order of the data, but in many cases WW will do well.
However, assume further that $\widetilde{x}_{1},\ldots,\widetilde{x}_{n}$ is
sorted based on the magnitude $|\widetilde{x}_{i}|$. In this case, WW fails as
well because it is based on runs around the median. Nevertheless, as we
demonstrate, with a suitably chosen sorting level $\beta$, the partial bubble
sorting test procedure works well.

To experiment with this scenario we consider Gaussian random vectors
$(\widetilde{X},X)$ with zero mean and unit variance for both $\widetilde{X}$
and $X$. The correlation coefficient between $\widetilde{X}$ and $X$ is set at
$\rho$. We generate $n$ repeated observations from this random vector and then
sort the observations according to $|\widetilde{X}|$. Now taking the second
coordinate of the sorted data yields the resulting random vector
$X_{1},\ldots,X_{n}$ on which we execute the bubble sort testing procedure of
Section 2. In the extreme case of $\rho=0$ the data $X_{1},\ldots,X_{n}$ is
iid and $H_{0}$ holds. However, as the magnitude of the correlation, $|\rho|$
grows towards $1$, traces of the sorting of $|\widetilde{X}|$ may start to be
detected in the data.

Figure 7 presents numerical results from an experiment with such data. In this
case we considered samples of size $n=1,000$ and several correlation levels,
$\rho$, ranging from $\rho=0.5$ which is hardly detectable to $\rho=0.9$ which
is easily detectable. For each $\rho$, we ran $30,000$ Monte Carlo repetitions
where for each repetition we repeated over a grid of sorting levels $\beta$.
We also used the same data for the WW test and for all $\rho$ values were not
able to reject $H_{0}$ more than $10\%$ ($\alpha$) of the runs. Further, the
KS test ($\beta=1$) also fails. Interestingly, the bubble sort based procedure
works well for some $\beta$ values and not for others with a dip in the
ability to detect for $\beta\approx 0.3$. Our key point here though is that by
choosing a suitable sorting level $\beta$, the sorting of the of the hidden
column can be detected via the bubble sort test procedure but not via KS or
WW.

### Example 2: Service Time Dependent Scheduling in Queues

We now consider a queueing system with a single server which serves jobs one
by one when they are present, or idles otherwise. For general concepts of
queueing models, and statistics of queues see for example [1], or references
there-in. In our context we consider a finite set of job arrivals to the queue
and observe the sequence of service durations of processed jobs,
$X_{1},\ldots,X_{n}$, where $X_{1}$ is the service duration of the first job
served, $X_{2}$ is the service duration of the second job served, and so
fourth.

Apriori, before being served in the queue, the jobs arriving to the system
have iid service durations with a known distribution $F_{0}$. The scheduling
policy used to determine which job to serve next considers the available jobs
in the queue at any time at which the server is free. This policy can affect
the system behaviour and the nature of $X_{1},\ldots,X_{n}$. With scheduling
policies such as first come first serve, last come first serve, or random
order of service, the output sequence of processed job durations
$X_{1},\ldots,X_{n}$ is iid because these policies do not depended on the
actual service time of a job for scheduling. These policies and all other
policies which induce iid service durations of processed jobs constitute
$H_{0}$ in this example

Other policies may yield non-independent $X_{1},\ldots,X_{n}$ and all such
policies constitute $H_{1}$. One such policy is a policy where the server
always selects the job with the smallest duration from the queue. In cases
where there is occasional congestion in the system and multiple jobs wait in
the queue simultaneously, this policy yields non-independent sequences of
service times $X_{1},\ldots,X_{n}$. To see this consider the extreme example
in which all jobs arrive to the queue at the same time. In such a case the
resulting sequence is completely sorted. Such a scheduling policy can
generally be detected well by the WW test or the bubble sort testing procedure
(with $\beta<1)$, but not by the KS test.

Figure 8: Example 2: Service Time Dependent Scheduling in Queues. Results for
$n=100$, $\alpha=0.1$ (black line), and several levels of $\sigma$. The
proportion of cases detected ($H_{0}$ rejected) is plotted as a function of
the sorting level $\beta$.

With other scheduling policies in $H_{1}$ it can be more difficult to confirm
that the service requirements of the output process are not independent. We
now focus on an experiment with one such policy. The policy is that the server
randomly selects either the smallest or the largest job in the queue to go
into service. This implies that the service requirements of the output process
are dependent but it also yields sequences that are more difficult for the WW
test to detect. Nevertheless, as we show via a numerical example, for suitably
chosen sorting level $\beta$ the bubble sort based procedure exhibits solid
power.

In our example we resort to an arrival process where the arrival times jobs
are iid random variables which each follow the law of $e^{A}$, where $A$ has a
normal distribution with expectation $\log(n)$ and a variance $\sigma^{2}$. As
$\sigma\to 0$ we get that all jobs arrive simulataniously and as $\sigma$
grows we get a bigger spread between jobs implying lower congestion. We resort
to $n=100$ jobs so the mean of $A$ is at around $4.6$. We take service
durations of the jobs arriving to the queue as iid standard uniform random
variables. With these parameters, to obtain some congestion in the queue we
test levels of $\sigma$ at $0.01,0.05,0.075$, and $0.1$. To compare we also
consider $\sigma=0.15$ and $\sigma=0.2$ at which point queue congestion is
much less common.

For each $\sigma$, we ran $100,000$ Monte Carlo repetitions where for each
repetition we repeated over a grid of sorting levels $\beta$. As in Example 1
we set $\alpha=0.1$ and also carry out the WW test on the simulated data. The
WW rejects $H_{0}$ with estimated probabilities,
$0.209,0.203,0.200,0.184,0.123$ and $0.100$ for the corresponding $\sigma$
values ranging from $0.01$ to $0.2$. In contrast, Figure 8 presents the
estimated power of the test as a function of the sorting level and for
different values of $\sigma$. It is evident that in this example for
$\beta\approx 0.25$ the bubble sort based procedure performs better than WW as
long as there is some congestion in the queue ($\sigma$ is not too high).

##  8 Discussion and Concluding Remarks

We demonstrated that partial sorting can be used to design test statistics
that are sensitive to both the order and the underlying distribution of the
data. More specifically, we showed that applying a fixed number of bubble sort
iterations to a data set with iid elements gives rise to a limiting curve that
we called the bubble sort curve. This curve can be regarded as a
generalization of the cdf and shares similar properties. We then turned to the
uniform distance between the empirical bubble sort curve and the limiting
bubble sort curve. Analogous to the Kolmogorov–Smirnov test statistic, we used
this distance to define the bubble sort test statistic. We characterized its
asymptotic distribution as a generalized Kolmogorov–Smirnov distribution that
only depends on the sorting level. We showed in examples that the bubble sort
test statistics may outperform the classical test statistics if data are not
obtained from independent observations.

There are several straightforward extensions of the current work. We focused
in this paper on the one-sample bubble sort test statistic, but extending the
theory to a two-sample test statistic is elementary. Another extension
involves the left tail of the data. In the current setup, the largest data
points are guaranteed to be fully sorted. Roughly speaking, this means that
the right tail of the data is compared to the right tail of the underlying cdf
in the bubble sort test statistic. Something similar can be achieved for the
left tail by multiplying all data point by -1 or by changing the bubble sort
algorithm to sort from right to left instead of from left to right.

There are also several new research directions suggested by our results. Two
main questions stand out. The first question is what happens if we use a
different sorting algorithm. This may give rise to a different limiting curve
and a different test statistic with potentially better statistical properties.
We suspect that less efficient sorting algorithms may be beneficial here. The
second question is whether we can combine bubble sort test statistics with
different sorting levels $\beta_{1}<\beta_{2}<\dotsb<\beta_{m}$ for the same
data set. The examples showed that there is no single best sorting level
$\beta$, but combining different sorting levels may lead to a more powerful
test.

## References

  * [1] Azam Asanjarani, Yoni Nazarathy, and Peter Taylor, _A survey of parameter and state estimation in queues_ , Queueing Systems 97 (2021), no. 1, 39–80.
  * [2] Owen Astrachan, _Bubble sort: an archaeological algorithmic analysis_ , ACM Sigcse Bulletin 35 (2003), no. 1, 1–5.
  * [3] L Beghin and E Orsingher, _On the maximum of the generalized Brownian bridge_ , Lithuanian Mathematical Journal 39 (1999), no. 2, 157–167.
  * [4] Robert H Berk and Douglas H Jones, _Goodness-of-fit test statistics that dominate the kolmogorov statistics_ , Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete 47 (1979), no. 1, 47–59.
  * [5] K Bisewski, H. M. Jansen, and Y Nazarathy, _partial-sorting-non-parametric-stats github organization_ , https://github.com/partial-sorting-non-parametric-stats, 2020.
  * [6] Bogdan Ćmiel, Tadeusz Inglot, and Teresa Ledwina, _Intermediate efficiency of some weighted goodness-of-fit statistics_ , Journal of Nonparametric Statistics 32 (2020), no. 3, 667–703.
  * [7] M. D. Donsker, _Justification and extension of Doob’s heuristic approach to the Kolmogorov–Smirnov theorems_ , The Annals of Mathematical Statistics 23 (1952), no. 2, 277–281.
  * [8] J. L. Doob, _Heuristic approach to the Kolmogorov–Smirnov theorems_ , The Annals of Mathematical Statistics 20 (1949), no. 3, 393–403.
  * [9] W. Feller, _On the Kolmogorov–Smirnov limit theorems for empirical distributions_ , The Annals of Mathematical Statistics 19 (1948), no. 2, 177–189.
  * [10] Helmut Finner and Veronika Gontscharuk, _Two-sample kolmogorov–smirnov-type tests revisited: old and new tests in terms of local levels_ , The Annals of Statistics 46 (2018), no. 6A, 3014–3037.
  * [11] Jerome H Friedman and Lawrence C Rafsky, _Multivariate generalizations of the wald-wolfowitz and smirnov two-sample tests_ , The Annals of Statistics (1979), 697–717.
  * [12] Valery Glivenko, _Sulla determinazione empirica delle leggi di probabilita_ , Gion. Ist. Ital. Attauri. 4 (1933), 92–99.
  * [13] Veronika Gontscharuk and Helmut Finner, _Asymptotics of goodness-of-fit tests based on minimum p-value statistics_ , Communications in Statistics-Theory and Methods 46 (2017), no. 5, 2332–2342.
  * [14] Veronika Gontscharuk, Sandra Landwehr, and Helmut Finner, _The intermediates take it all: Asymptotics of higher criticism statistics and a powerful alternative based on equal local levels_ , Biometrical Journal 57 (2015), no. 1, 159–180.
  * [15] , _Goodness of fit tests in terms of local levels with special emphasis on higher criticism tests_ , Bernoulli 22 (2016), no. 3, 1331–1363.
  * [16] Leah Jager and Jon A Wellner, _Goodness-of-fit tests via phi-divergences_ , The Annals of Statistics 35 (2007), no. 5, 2018–2053.
  * [17] Arnold Janssen, _Global power functions of goodness of fit tests_ , The Annals of Statistics 28 (2000), no. 1, 239–253.
  * [18] Song-Hee Kim and Ward Whitt, _The power of alternative kolmogorov-smirnov tests based on transformations of the data_ , ACM Transactions on Modeling and Computer Simulation (TOMACS) 25 (2015), no. 4, 1–22.
  * [19] Donald Ervin Knuth, _The art of computer programming_ , vol. 3, Pearson Education, 1997.
  * [20] Andrey Kolmogorov, _Sulla determinazione empirica di una legge di distribuzione_ , Giorn. Inst. Ital. Attuari 4 (1933), 83–91.
  * [21] David M Mason and John H Schuenemeyer, _A modified kolmogorov-smirnov test sensitive to tail alternatives_ , The annals of Statistics (1983), 933–946.
  * [22] Amit Moscovich, Boaz Nadler, and Clifford Spiegelman, _On the exact Berk-Jones statistics and their $p$-value calculation_, Electronic Journal of Statistics 10 (2016), no. 2, 2329–2354.
  * [23] A Rényi, _On a group of problems in the theory of order statistics_ , Selected Transl. in Math. Statist. and Probab. Amer. Math. Soc 13 (1973), 289–298.
  * [24] N Smirnoff, _Sur les écarts de la courbe de distribution empirique_ , Matematicheskii Sbornik 48 (1939), no. 1, 3–26.
  * [25] Nikolai V Smirnov, _On the estimation of the discrepancy between empirical curves of distribution for two independent samples_ , Bull. Math. Univ. Moscou 2 (1939), no. 2, 3–14.
  * [26] Michael A Stephens, _An appreciation of kolmogorov’s 1933 paper_ , Tech. report, Stanford department of statistics, 1992.
  * [27] A. W. van der Vaart, _Asymptotic statistics_ , Cambridge University Press, New York, 1998.
  * [28] Abraham Wald and Jacob Wolfowitz, _On a test whether two samples are from the same population_ , The Annals of Mathematical Statistics 11 (1940), no. 2, 147–162.
  * [29] Ward Whitt, _Stochastic-process limits: an introduction to stochastic-process limits and their application to queues_ , Springer, New York, 2002.

KB’s research was funded by SNSF Grant 200021-196888. HMJ and YN were partly
funded by Australian Research Council (ARC) Discovery Project DP180101602.


# Towards Transformer-based Homogenization of Satellite Imagery for Landsat-8
and Sentinel-2

Venkatesh Thirugnana Sambandham
Department of Computer Science
Otto-von-Guericke University Magdeburg
<EMAIL_ADDRESS>
Konstantin Kirchheim
Department of Computer Science
Otto-von-Guericke University Magdeburg
<EMAIL_ADDRESS>
Sayan Mukhopadhaya
BASF Digital Farming GmbH
<EMAIL_ADDRESS>
Frank Ortmeier
Department of Computer Science
Otto-von-Guericke University Magdeburg
<EMAIL_ADDRESS>

###### Abstract

Landsat-8 (NASA) and Sentinel-2 (ESA) are two prominent multi-spectral imaging
satellite projects that provide publicly available data. The multi-spectral
imaging sensors of the satellites capture images of the earth’s surface in the
visible and infrared region of the electromagnetic spectrum. Since the
majority of the earth’s surface is constantly covered with clouds, which are
not transparent at these wavelengths, many images do not provide much
information. To increase the temporal availability of cloud-free images of a
certain area, one can combine the observations from multiple sources. However,
the sensors of satellites might differ in their properties, making the images
incompatible. This work provides a first glance at the possibility of using a
transformer-based model to reduce the spectral and spatial differences between
observations from both satellite projects. We compare the results to a model
based on a fully convolutional UNet architecture. Somewhat surprisingly, we
find that, while deep models outperform classical approaches, the UNet
significantly outperforms the transformer in our experiments.

##  1 Introduction

A large number of satellite clusters are constantly monitoring the physical
characteristics of the earth’s surface using active and passive
electromagnetic sensors. These satellites capture data across multiple bands
of the electromagnetic spectrum, which can be utilized in several important
applications, such as crop health monitoring [16, 14]. For these applications,
achieving a high temporal resolution for a particular area of the earth’s
surface is advantageous. Landsat-8 and Sentinel-2 are two prominent earth
observation satellite missions that provide publicly available data. A
comprehensive study on the revisiting time of Landsat-8, 9 as well as
Sentinel-2 A and B found that a combination of all 4 sensors could potentially
reduce the duration between two revisits to 2.3 days on average [10]. However,
not all of the available images are usable due to occlusions caused by clouds,
which are not transparent at the corresponding wavelengths. Furthermore, the
sensors have different properties and consequentially, several additional
challenges have to be considered:

  1. 1.

The spectral differences caused by different characteristics of the imaging
sensors used in each project have to be minimized.

  2. 2.

The differences in the spatial resolution of the image sources have to be
handled.

For example, consider fig. 1: loosely speaking, our goal is to generate image
with the quality of fig. 1(a) by a combination of fig. 1(b) and the
panchromatic band in fig. 1(c).

(a) Sentinel-2

(b) Landsat-8

(c) Panchromatic

Figure 1: A comparison of Sentinel-2 and Landsat-8 RGB composites of the same
region, as well as the corresponding panchromatic band from Landsat-8. We see
that the resolution of the Sentinel-2 is higher than the resolution of the
Landsat-8 image.

Approaches based on Deep Learning [18, 9] play a key role in several
applications in the field of remote sensing [21]. Extensive research was
conducted on the topic of deep learning based super-resolution (SR), with some
techniques showing exceptional improvements in the perceptual quality of low-
resolution images [1]. In particular, the UNet [15], which is based on an
encoder-decoder architecture, provides good results for super-resolution [13].

In this work, we explore the use of vision transformers for the harmonization
of images from Landsat-8 and Sentinel-2. We hypothesize that transformers can
outperform UNet-based models in the considered setting, due to the ability of
the incorporated attention mechanism to account for complex, non-local
interaction between variables.

##  2 Background

Reducing the spectral differences between the images acquired from multiple
sensors is a well-researched topic [19, 4, 17]. Most previous works derive an
explicit compensation function to convert the reflectances of one image source
to another.

Deep Learning based approaches have recently been applied to image super-
resolution and re-construction applications [2, 20]. Generally, deeper
Convolutional Neural Network (CNN) architectures (with residual connections
[5] and feature concatenation [6]) yield better results in common computer
vision tasks. Deep Learning based super-resolution architectures were also
applied to enhance the spatial resolution of low-resolution multi-spectral
images [3, 8, 7].

Motivated by their superior performance on Natural Language Processing tasks,
transformer-based models have recently been adapted for Computer Vision tasks.
Vision transformers (ViT) like the SwinIR [11] and the ESRT [12] have
successfully been applied to image reconstruction and super-resolution.

##  3 Models and Dataset

###  3.1 Data

We collect pairs of images from both satellite sources that were taken over
the Area Of Interest (AOI) on the same day. The AOI considered in this work
covers a large region of Germany. The areas used for training, validation, and
testing are disjoint from each other to avoid target leakage. The high-
resolution (HR) images from Sentinel-2 act as ground truth images; the
objective of the model is to homogenize the common bands from the low-
resolution (LR) images from Landsat-8 with respect to the HR source.

###  3.2 Preprocessing

After we fetch data from both sources, the images from various bands are
converted to reflectances, filtered for cloud coverage, and then corrected for
inter-sensor misregistration issues. The images are finally broken down into
smaller patches of size $1.28$ km$\times 1.28$ km. Each patch consists of
7-channels which include multi-spectral bands Blue, Green, Red, Near Infrared
(NIR), two Short Wave Infrared (SWIR), and the panchromatic image from the LR
source. After patching there were 63,740 patches for training, 22,344 patches
for testing, and 11,908 patches for validation. The LR images are bicubically
upsampled to the size of HR images to act as the input to the models.

###  3.3 Neural Network Models and Training

In this study, we compare two deep super-resolution architectures: a UNet [15,
13] and a slightly modified Efficient Super Resolution Transformer (ESRT)
[12]. The pixel-shuffle upsampling block in the ESRT architecture was removed,
since the images were already upsampled to the target size in the pre-
processing step. We train the models for 50 epochs, using the Adam optimizer
with a mini-batch size of 20, minimizing the $\ell_{1}$ norm between the
ground truth and the prediction. We use $1\times 10^{-5}$ as our learning
rate. Each model was trained on a single 32 GB NVIDIA V100 GPU. Since the GPU
could not hold 20 patches at a time while training the transformer-based
model, we used gradient accumulation over four batches of size five.

##  4 Results

We evaluate the models using a pixel-wise Normalized Root Mean Squared Error
(NRMSE) and perception-based Structural Similarity Index Measure (SSIM) as
described in table 1. The distribution of the calculated metrics from the test
data set is given as box plots in fig. 2. In the quantitative evaluation, we
observe that the UNet significantly outperformed the baseline upsampling and
the ESRT architecture with a mean SSIM of 0.91. Although the transformer
outperformed the baseline, it is not on par with the UNet in terms of both
metrics.

A qualitative comparison of the RGB composites of samples at different
performance levels (as measured by the UNets NRMSE) is given in fig. 3.
Significant improvement in the spatial quality can be observed in the images
generated by both neural networks. However, some artifacts appear in the
images generated by the ESRT, especially in the corners of the patches.
Another unsurprising observation is that models fail to reconstruct high-
frequency components of the image, like settlements, whose presence can cause
significant drops in performance metrics.

(a) NRMSE

(b) SSIM

Figure 2: Distribution of metrics calculated for different methods on the test set. Table 1: Mean and standard deviation of performance metrics on the test set. | SSIM | NRMSE
---|---|---
Bicubic | 0.8763 $\pm$ 0.043 | 0.1839 $\pm$ 0.056
UNet | 0.9114 $\pm$ 0.028 | 0.1460 $\pm$ 0.048
ESRT | 0.8898 $\pm$ 0.039 | 0.1705 $\pm$ 0.064
Figure 3: Results from sample patches across different levels of UNet NRMSE
metrics over the test set. Outliers often contain high-frequency components
that are not present in the panchromatic band.

##  5 Conclusion

We presented preliminary results for a transformer-based approach for
harmonizing multi-spectral images of the earths surface provided by the
Landsat-8 and Sentinel-2 satellites. The model constitutes a component in a
larger pipeline that enables higher temporal resolution of the observations of
a specific region, which can be utilized in several downstream applications,
for example, crop health monitoring. In first experiments, we somewhat
surprisingly find that, while deep models outperform classical approaches like
bilinear interpolation, the UNet outperforms the more complex transformer.

A possible explanation for the superior performance of the UNet is that our
initial hypothesis does not hold: non-local interactions might not be as
relevant for this application, and the UNets translation invariance, coupled
with a sufficiently large receptive field (i.e., sufficient depth),
constitutes a better inductive bias for the considered task.

## References

  * [1] Saeed Anwar, Salman Khan, and Nick Barnes.  A deep journey into super-resolution: A survey.  ACM Comput. Surv., 53(3), may 2020.
  * [2] Chao Dong, Chen Change Loy, Kaiming He, and Xiaoou Tang.  Image super-resolution using deep convolutional networks.  CoRR, abs/1501.00092, 2015.
  * [3] Mikel Galar, Ruben Sesma, Christian Ayala, and Carlos Aranda.  Super-resolution for sentinel-2 images.  International Archives of the Photogrammetry, Remote Sensing & Spatial Information Sciences, 2019.
  * [4] Alemu Gonsamo and Jing M. Chen.  Spectral response function comparability among 21 satellite sensors for vegetation monitoring.  IEEE Transactions on Geoscience and Remote Sensing, 51(3):1319–1335, 2013.
  * [5] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, page 770–778, 2016.
  * [6] Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger.  Densely connected convolutional networks.  Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, page 4700–4708, 2017.
  * [7] Sani M Isa, Suharjito, Gede Putera Kusuma, and Tjeng Wawan Cenggoro.  Supervised conversion from landsat-8 images to sentinel-2 images with deep learning.  European Journal of Remote Sensing, 54(1):182–208, 2021.
  * [8] Nicolas Latte and Philippe Lejeune.  Planetscope radiometric normalization and sentinel-2 super-resolution (2.5 m): A straightforward spectral-spatial fusion of multi-satellite multi-sensor images using residual convolutional neural networks.  Remote Sensing, 12(15):2366, 2020.
  * [9] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  Nature, 521(7553):436–444, 2015.
  * [10] Jian Li and Baozhang Chen.  Global revisit interval analysis of landsat-8 -9 and sentinel-2a -2b data for terrestrial monitoring.  Sensors, 20(22), 2020.
  * [11] Jingyun Liang, Jiezhang Cao, Guolei Sun, Kai Zhang, Luc Van Gool, and Radu Timofte.  Swinir: Image restoration using swin transformer.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 1833–1844, 2021.
  * [12] Zhisheng Lu, Juncheng Li, Hong Liu, Chaoyan Huang, Linlin Zhang, and Tieyong Zeng.  Transformer for single image super-resolution.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 457–466, 2022.
  * [13] Xiao-Jiao Mao, Chunhua Shen, and Yu-Bin Yang.  Image denoising using very deep fully convolutional encoder-decoder networks with symmetric skip connections.  CoRR, abs/1603.09056, 2016.
  * [14] Sayan Mukhopadhaya.  Land use and land cover change modelling using ca-markov case study: Deforestation analysis of doon valley.  Journal of Agroecology and Natural Resource Management, 3:1–5, 3 2016.
  * [15] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  International Conference on Medical image computing and computer-assisted intervention, pages 234–241, 2015.
  * [16] David P Roy, Michael A Wulder, Thomas R Loveland, Curtis E Woodcock, Richard G Allen, Martha C Anderson, Dennis Helder, James R Irons, David M Johnson, Robert Kennedy, et al.  Landsat-8: Science and product vision for terrestrial global change research.  Remote sensing of Environment, 145:154–172, 2014.
  * [17] Daniel Scheffler, David Frantz, and Karl Segl.  Spectral harmonization and red edge prediction of landsat-8 to sentinel-2 using land cover optimized multivariate regressors.  Remote Sensing of Environment, 241:111723, 2020.
  * [18] Jürgen Schmidhuber.  Deep learning in neural networks: An overview.  Neural networks, 61:85–117, 2015.
  * [19] Jose Luis Villaescusa-Nadal, Belen Franch, Jean-Claude Roger, Eric F. Vermote, Sergii Skakun, and Chris Justice.  Spectral adjustment model’s analysis and application to remote sensing data.  IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 12(3):961–972, 2019.
  * [20] Kai Zhang, Wangmeng Zuo, Yunjin Chen, Deyu Meng, and Lei Zhang.  Beyond a gaussian denoiser: Residual learning of deep CNN for image denoising.  CoRR, abs/1608.03981, 2016.
  * [21] Xiao Xiang Zhu, Devis Tuia, Lichao Mou, Gui-Song Xia, Liangpei Zhang, Feng Xu, and Friedrich Fraundorfer.  Deep learning in remote sensing: A comprehensive review and list of resources.  IEEE Geoscience and Remote Sensing Magazine, 5(4):8–36, 2017.