arXiv:1611.00630v2 [] 27 Jun 2017

The accumulated persistence function, a new
useful functional summary statistic for
topological data analysis, with a view to brain
artery trees and spatial point process
applications
Christophe A.N. Biscio
Department of Mathematical Sciences, Aalborg University, Denmark
and
Jesper Møller
Department of Mathematical Sciences, Aalborg University, Denmark
June 28, 2017

Abstract
We start with a simple introduction to topological data analysis where the most
popular tool is called a persistent diagram. Briefly, a persistent diagram is a multiset of points in the plane describing the persistence of topological features of a
compact set when a scale parameter varies. Since statistical methods are difficult
to apply directly on persistence diagrams, various alternative functional summary
statistics have been suggested, but either they do not contain the full information of
the persistence diagram or they are two-dimensional functions. We suggest a new
functional summary statistic that is one-dimensional and hence easier to handle,
and which under mild conditions contains the full information of the persistence diagram. Its usefulness is illustrated in statistical settings concerned with point clouds
and brain artery trees. The appendix includes additional methods and examples,
together with technical details. The R-code used for all examples is available at
http://people.math.aau.dk/~christophe/Rcode.zip.

Keywords: clustering, confidence region, global rank envelope, functional boxplot, persistent homology, two-sample test.
1

1 Introduction
Statistical methods that make use of algebraic topological ideas to summarize and visualize complex data are called topological data analysis (TDA). In particular persistent
homology and a method called the persistence diagram are used to measure the persistence of topological features. As we expect many readers may not be familiar with
these concepts, Section 1.1 discusses two examples without going into technical details
though a few times it is unavoidable to refer to the terminology used in persistent homology. Section 1.2 discusses the use of the persistence diagram and related summary
statistics and motivates why in Section 1.3 a new functional summary statistic called
the accumulative persistence function (APF) is introduced. The remainder of the paper demonstrates the use of the APF in various statistical settings concerned with point
clouds and brain artery trees.

1.1 Examples of TDA
The mathematics underlying TDA uses technical definitions and results from persistent homology, see Fasy et al. (2014) and the references therein. This theory will not be
needed for the present paper. Instead we provide an understanding through examples
of the notion of persistence of k-dimensional topological features for a sequence of compact subsets Ct of the d-dimensional Euclidean space R d , where t ≥ 0, k = 0, 1, . . . , d − 1,
and either d = 2 (Section 1.1.1) or d = 3 (Section 1.1.2). Recalling that a set A ⊆ R d is

path-connected if any two points in A are connected by a curve in A, a 0-dimensional
topological feature of a compact set C ⊂ R d is a maximal path-connected subset of C,
also called a connected component of C. The meaning of a 1-dimensional topological feature is simply understood when d = 2, and we appeal to this in a remark at the end of
Section 1.1.1 when d = 3.

2

t=0

t = 0.5

t = 0.62

0.4

Birth time

0.8

1.0
0.5

APF(m)
0.0

APF0
APF1

0.0

0.4

k=0
k=1

0.0

Death time

0.8

1.5

t = 0.75

0.0

0.4

0.8

m

Figure 1: The four first panels show a simple example of spherically growing circles
centred at (-1,-1), (1,-1), and (0,1) and all with radii 0.5 when time is 0. The fifth panel
shows the persistence diagram for the connected components (k = 0) and the loops
(k = 1). The final panel shows the corresponding accumulated persistence functions.
1.1.1 A toy example
Let C ⊂ R2 be the union of the three circles depicted in the top-left panel of Figure 1. The
three circles are the 0-dimensional topological features (the connected components) of C,
as any curve that goes from a circle to another will be outside C. The complement R2 \ S
has four connected components, one of which is unbounded, whilst the three others are
the 1-dimensional topological features of C, also called the loops of C (the boundary of
each bounded connected component is a closed curve with no crossings; in this example
the closed curve is just a circle).
For t ≥ 0, let Ct be the subset of points in R2 within distance t from C. Thinking of t
as time, Ct results when each point on C grows as a disc with constant speed one. The
algebraic topology (technically speaking the Betti numbers) changes exactly at the times
t = 0, 0.5, 0.62, 0.75, see the first four panels of Figure 1: For each topological dimension
(k)

k = 0, 1, let ti

denote the time of the ith change. First, C0 = C has three connected
(0)

components and three loops as given above; we say that they are born at time t1

3

=

(1)

t1

= 0 (imaging there was nothing before time 0). Second, the loops disappear and

two connected components merge into one connected component; we say that the loops
(0)

(1)

and one of the connected components die at time t2 = t2 = 0.5; since the two merging
connected components were born at the same time, it is decided uniformly at random
which one should die respectively survive; the one which survives then represent the
(0)

(1)

new merged connected component. Third, at time t3 = t3 = 0.62, a new loop is born
and the two connected component merge into one which is represented by the oldest
(first born) connected component whilst the other connected component (the one which
(0)

was retained when the two merged at time t2

= 0.5) dies; the remaining connected

component "lives forever" after time 0.62 and it will be discarded in our analysis. Finally,
(1)

at time t4 = 0.75, the loop dies.
Hence, for each k = 0, 1, there is a multiset of points specifying the appearance and
disappearance of each k-dimensional topological feature as t grows. A scatter plot of
these points is called a persistence diagram, see Figure 1 (bottom-middle panel): For k = 0
(the connected components), the points are (0, 0.5) and (0, 0.62) (as (0, ∞) is discarded
from the diagram) with multiplicities 1 and 1, respectively; and for k = 1 (the loops),
the points are (0, 0.5) and (0.62, 0.75) with multiplicities 3 and 1, respectively. The term
"persistence" refers to that distant connected components and large loops are present
for a long time; which here of course just corresponds to the three circles/connected
components of C and their loops persist for long whilst the last appearing loop has a
short lifetime and hence is considered as "noise".
Usually in practice C is not known but a finite subset of points { x1 , . . . , x N } has been
collected as a sample on C, possibly with noise. Then we redefine Ct as the union of
closed discs of radius t and with centres given by the point cloud. Hence the connected
components of C0 are just the points x1 , . . . , x N , and C0 has no loops. For t > 0, it
is in general difficult to directly compute the connected components and loops of Ct ,
but a graph in R m (with m ≥ 2) can be constructed so that its connected components
correspond to those of Ct and moreover the triangles of the graph may be filled or not
in a way so that the loops of the obtained triangulation correspond to those of Ct .
4

Remark: Such a construction can also be created in the case where Ct is the union of
d-dimensional closed balls of radius t and with centres given by a finite point pattern

{ x1 , . . . , x N } ⊂ R d . The construction is a so-called simplicial complex such as the Čechcomplex, where m may be much larger than d, or the Delaunay-complex (or alphacomplex), where m = d, and a technical result (the Nerve Theorem) establishes that
it is possible to identify the topological features of Ct by the Čech or Delaunay-complex,
see e.g. Edelsbrunner and Harer (2010). It is unnecessary for this paper to understand
the precise definition of these notions, but as d = 2 or d = 3 is small in our examples, it
is computationally convenient to use the Delaunay-complex. When d = 3, we may still
think of a 1-dimensional topological feature as a loop, i.e. a closed curve with no crossings; again the simplicial complex is used for the "book keeping" when determining the
persistence of a loop. For example, a 2-dimensional sphere has no loops, and a torus in
R3 has two. Finally, when d ≥ 3, a k-dimensional topological feature is a k-dimensional

manifold (a closed surface if k = 2) that cannot "be filled in", but for this paper we omit
the precise definition since it is technical and not needed.

1.1.2 Persistent homology for brain artery trees

The left panel of Figure 2 shows an example of one of the 98 brain artery trees analysed in
Bendich et al. (2016). The data for each tree specifies a graph in R3 consisting of a dense
cloud of about 105 vertices (points) together with the edges (line segments) connecting
the neighbouring vertices; further details about the data are given in Section 2.2. As
in Bendich et al. (2016), for each tree we only consider the k-dimensional topological
features when k = 0 or k = 1, using different types of data and sets Ct as described
below. Below we consider the tree in Figure 2 and let B ⊂ R3 denote the union of its
edges.
Following Bendich et al. (2016), if k = 0, let Ct = {( x, y, z) ∈ B : z ≤ t} be the sub-level

set of the height function for the tree at level t ≥ 0 (assuming Ct is empty for t < 0). Thus
the 0-dimensional topological features at "time/level" t are the connected components of

5

25
20
15
5

10

Death time

80
60
40

Death time

20

0

0
0

20

40

60

80

Birth time

0

5

10

15

20

25

Birth time

Figure 2: A brain artery tree with the "water level" indicated (left panel) and the persistence diagrams of connected components (middle panel) and loops (right pane).
Ct . As illustrated in the left panel of Figure 2, instead of time we may think of t as "water
level": As the water level increases, connected components of the part of B surrounded
by water (the part in blue) may be born or die; we refer to this as sub-level persistence.
As in Section 1.1.1, we represent the births and deaths of the connected component in a
persistence diagram which is shown in Figure 2 (middle panel). The persistence of the
connected components for all brain artery trees will be studied in several examples later
on.
As in Bendich et al. (2016), if k = 1, we let B be represented by a point pattern C of 3000
points subsampled from B, and redefine Ct to be the union of balls of radii t ≥ 0 and
centres given by C (as considered in the remark at the end of Section 1.1.1). The loops
of Ct are then determined by the corresponding Delaunay-complex. The right panel of
Figure 2 shows the corresponding persistence diagram. The persistence of the loops for
all trees will be studied in the examples to follow.

1.2 Further background and objective
The persistence diagram is a popular graphical representation of the persistence of the
topological features of a sequence of compact sets Ct ⊂ R d , t ≥ 0. As exemplified

above it consists for each topological dimension k = 0, . . . , d − 1 of a multiset PDk of

points (bi , di ) with multiplicities ci , where bi and di is a pair of birth-death times for a
k-dimensional topological feature obtained as time t grows.

6

In the majority of literature on TDA, including the analysis in Bendich et al. (2016) of
brain artery trees, long lifetimes are of main interest whereas short lifetimes are considered as topological noise. Short lifetimes are of interest in the study of complex structures such as branch polymers and fractals, see MacPherson and Schweinhart (2012);
and for brain artery trees Bendich et al. (2016) noticed in one case that "not-particularlyhigh persistence have the most distinguishing power in our specific application". In
our examples we demonstrate that short lifetimes will also be of key interest in many
situations, including when analysing the brain artery trees dataset from Bendich et al.
(2016).
Chazal et al. (2013) and Chen et al. (2015) note that it is difficult to apply statistical methodology to persistent diagrams. Alternative functional summary statistics have been suggested: Bubenik (2015) introduces a sequence of one-dimensional functions called the
persistent landscape, where his first function is denoted λ1 and is considered to be of
main interest, since it provides a measure of the dominant topological features, i.e. the
longest lifetimes; therefore we call λ1 the dominant function. Chazal et al. (2013) introduce the silhouette which is a weighted average of the functions of the persistent
landscape, where the weights control whether the focus is on topological features with
long or short lifetimes. Moreover, Chen et al. (2015) consider a kernel estimate of the
intensity function for the persistent diagram viewed as a point pattern. The dominant
function, the silhouette, and the intensity estimate are one-dimensional functions and
hence easier to handle than the persistence diagram, however, they provide selected
and not full information about the persistence diagram. In Section 1.3, we introduce
another one-dimensional functional summary statistic called the accumulative persistence function and discuss its advantages and how it differs from the existing functional
summary statistics.

7

1.3 The accumulated persistence function
For simplicity and specificity, for each topological dimension k = 0, 1, . . . , d − 1, we

always assume that the persistence diagram PDk = {(b1 , d1 , c1 ), . . . , (bn , dn , cn )} is such

that n < ∞ and 0 ≤ bi < di < ∞ for i = 1, . . . , n. This assumption will be satisfied

in our examples (at least with probability one). Often in the TDA literature, PDk is
transformed to the rotated and rescaled persistence diagram (RRPD) given by RRPDk =

{(m1 , l1 , c1 ), . . . , (mn , ln , cn )}, where mi = (bi + di )/2 is the meanage and li = di − bi is
the lifetime. This transformation is useful when defining our accumulative persistence
function (APF) by

n

APFk (m) =

∑ c i l i 1 ( m i ≤ m ),

i =1

m ≥ 0,

(1)

where 1(·) is the indicator function and we suppress in the notation that APFk is a function of RRPDk . The remainder of this section comments on this definition.
Formally speaking, when RRPDk is considered to be random, it is viewed as a finite
point process with multiplicities, see e.g. Daley and Vere-Jones (2003). It what follows it
will always be clear from the context whether PDk and RRPDk are considered as being
random or observed, and hence whether APFk is a deterministic or random function.
In the latter case, because APFk (m) is an accumulative function, its random fluctuations
typically increase as m increases.
Depending on the application, the jumps and/or the shape of APFk may be of interest
as demonstrated later in our examples. A large jump of APFk corresponds to a large
lifetime (long persistence). In the simple example shown in Figure 1, both jumps of
APF0 are large and indicate the three connected components (circles), whilst only the
first jump of APF1 is large and indicates the three original loops. For the more complicated examples considered in the following it may be hard to recognize the individual
jumps. In particular, as in the remark at the end of Section 1.1.1, suppose Ct is the
union of d-dimensional balls of radius t and with centres given by a finite point pattern

{ x1 , . . . , x N } ⊂ R d . Roughly speaking we may then have the following features as il-

lustrated later in Example 1. For small meanages m, jumps of APF0 (m) correspond to
8

balls that merge together for small values of t. Thus, if the point pattern is aggregated
(e.g. because of clustering), we expect that APF0 (m) has jumps and is hence large for
small meanages m, whilst if the point pattern is regular (typically because of inhibition
between the points), we expect the jumps of APF0 (m) to happen and to be large for modest values of m (as illustrated later in the middle panel of Figure 3 considering curves for
the Matérn cluster process and the determinantal point process). For large meanages,
jumps of APF0 (m) are most likely to happen in the case of aggregation. Accordingly, the
shape of APF0 can be very different for these two cases (as illustrated in the first panel of
Figure 3). Similar considerations lead us to expect different shapes of APF1 for different
types of point patterns; we expect that APF1 (m) is large respective small for the case of
aggregation respective regularity when m is small, and the opposite happens when m is
large (as illustrated in the last panel of Figure 3).
Clearly, RRPDk is in one-to-one correspondence to PDk . In turn, if all ci = 1 and the
mi are pairwise distinct, then there is a one-to-one correspondence between RRPDk
and its corresponding APFk . For k = 0, this one-to-one correspondence would easily be lost if we had used bi in place of mi in (1). We need to be careful with not
over-stating this possible one-to-one correspondence. For example, imagine we want
(1)

to compare two APFs with respect to Lq -norm (1 ≤ q ≤ ∞) and let PDk

(2)

and PDk

be the underlying persistence diagrams. However, when points (b, d) close to the diagonal are considered as topological noise (see Section 1.1.1), usually the so-called bot(1)

(2)

tleneck distance W∞ (PDk , PDk ) is used, see e.g. Fasy et al. (2014). Briefly, for ǫ > 0,
√
let N = {(b, d) : b ≤ d, l ≤ 2ǫ} be the set of points at distance 2ǫ of the diagonal
in the persistence diagram, and let S(b, d) = {( x, y) : | x − b| ≤ ǫ, |y − d| ≤ ǫ} be
the square with center (b, d), sides parallel to the b- and d-axes, and of side length 2ǫ.
(1)

(2)

(2)

Then W∞ (PDk , PDk ) ≤ ǫ if PDk

(bi , di ) a point of
(1)

(1)
PDk

has exactly one point in each square S(bi , di ), with

(repeating this condition ci times). However, small values of

(2)

W∞ (PDk , PDk ) does not correspond to closeness of the two corresponding APFs with
respect to Lq -norm.
Note that the dominant function, the silhouette, and the intensity estimate (see Sec-

9

tion 1.2) are in general not in a one-to-one correspondence with RRPDk . Like these
functions, APFk is a one-dimensional function, and so it is easier to handle than the
sequence of functions for the persistent landscape in Bubenik (2015) and the intensity
estimate in Chen et al. (2015) — e.g. confidence regions become easier to plot. Contrary
to the dominant function and the silhouette, the APF provides information about topological features without distinguishing between long and short lifetimes.

1.4 Outline
Our paper discusses various methods based on APFs in different contexts and illustrated by simulation studies related to spatial point process applications and by reanalysing the brain artery trees dataset previously analysed in Bendich et al. (2016). Section 2 specifies the setting for these examples. Sections 3, 4, and 5 consider the case of a
single APF, a sample of APFs, and two samples of APFs, respectively. Further examples
and details appear in Appendix A-F.

2 Datasets
2.1 Simulated data
In our simulation studies we consider a planar point cloud, i.e. a finite point pattern

{ x1 , . . . , x N } ⊂ R2 , and study as at the end of Section 1.1.1 how the topological features
of Ct , the union of closed discs of radii t and centred at x1 , . . . , x N , change as t grows.
Here { x1 , . . . , x N } will be a realisation of a point process X ⊂ R2 , where the count N is
finite. Thus PDk and RRPDk can be viewed as finite planar point processes (with multiplicities) and APFk as a random function. Note that N may be random, and conditional
on N, the points in X are not necessarily independent and identically distributed (IID).
This is a common situation in spatial statistics, e.g. if the focus is on the point process X
10

and the purpose is to assess the goodness of fit for a specified point process model of X
when { x1 , . . . , x N } is observed.

2.2 Brain artery trees dataset
The dataset in Bendich et al. (2016) comes from 98 brain artery trees which can be included within a cube of side length at most 206 mm; one tree is excluded "as the java/matlab
function crashed" (e-mail correspondence with Sean Skwerer). They want to capture
how the arteries bend through space and to detect age and gender effects. For k = 0,
sub-level persistence of the connected components of each tree represented by a union
of line segments is considered, cf. Section 1.1.2; then for all meanages, mi ≤ 137; and the
number of connected components is always below 3200. For k = 1, persistence of the

loops for the union of growing balls with centres at a point cloud representing the tree
is considered, cf. Section 1.1.2; the loops have a finite death time but some of them do
not die during the allocated time T = 25 (that is, Bendich et al. (2016) stop the growth
of balls when t > 25). Thus we shall only consider meanages mi ≤ 25; then the number
of loops is always below 2700. For each tree and k = 0, 1, most ci = 1 and sometimes
ci > 1.
For each tree and k = 0, 1, Bendich et al. (2016) use only the 100 largest lifetimes in their
analysis. Whereas their principal component analysis clearly reveal age effects, their
permutation test based on the mean lifetimes for the male and females subjects only
shows a clear difference when considering PD1 . Accordingly, when demonstrating the
usefulness of APF0 and APF1 , we will focus on the gender effect and consider the same
95 trees as in Bendich et al. (2016) (two transsexual subjects are excluded) obtained from
46 female subjects and 49 male subjects; in contrast to Bendich et al. (2016), we consider
all observed meanages and lifetimes. In accordance to the allocated time T = 25, we
need to redefine APF1 by
n

APF1 (m) =

∑ ci li 1(mi ≤ m, mi + li /2 ≤ T ),

i =1

11

m ≥ 0.

(2)

For simplicity we use the same notation APF1 in (1) and (2); although all methods and
results in this paper will be presented with the definition (1) in mind, they apply as well
when considering (2).
Finally, we write APFkF and APFkM to distinguish between APF’s for females and males,
respectively.

3 A single accumulated persistence function
There exists several constructions and results on confidence sets for persistence diagrams when the aim is to separate topological signal from noise, see Fasy et al. (2014),
Chazal et al. (2014), and the references therein. Appendix A and its accompanying Example 5 discuss the obvious idea of transforming such a confidence region into one for
an accumulate persistence function, where the potential problem is that the bottleneck
metric is used for persistence diagrams and this is not corresponding to closeness of
APFs, cf. Section 1.3. In this section we focus instead on spatial point process model
assessment using APFs or more traditional tools.
Suppose a realization of a finite spatial point process X0 has been observed and copies
X1 , . . . , Xr have been simulated under a claimed model for X0 so that the joint distribution of X0 , X1 , . . . , Xr should be exchangeable. That is, for any permutation (σ0 , . . . , σr )
of (0, . . . , r ), (Xσ0 , . . . , Xσr ) is claimed to be distributed as (X0 , . . . , Xr ); e.g. this is the
case if X0 , X1 , . . . , Xr are IID. This is a common situation for model assessment in spatial
point process analysis when a distribution for X0 has been specified (or estimated), see
e.g. Baddeley et al. (2015) and Møller and Waagepetersen (2016). Denote the APFk s for
X0 , . . . , Xr by A0 , . . . , Ar , respectively, and the null hypothesis that the joint distribution
of A0 , . . . , Ar is exchangeable by H0 . Adapting ideas from Myllymäki et al. (2016), we

will discuss how to construct a goodness-of-fit test for H0 based on a so-called global
rank envelope for A0 ; their usefulness will be demonstrated in Example 1.

12

In functional data analysis, to measure how extreme A0 is in comparison to A1 , . . . , Ar , a
so-called depth function is used for ranking A0 , . . . , Ar , see e.g. López-Pintado and Romo
(2009). We suggest using a depth ordering called extreme rank in Myllymäki et al. (2016):
Let T > 0 be a user-specified parameter chosen such that it is the behaviour of A0 (m)
for 0 ≤ m ≤ T which is of interest. For l = 1, 2, . . ., define the l-th bounding curves of
A0 , . . . , Ar by
l
Alow
(m) = min l Ai (m)
i =0,...,r

and

l
Aupp
(m) = max l Ai (m),
i =0,...,r

0 ≤ m ≤ T,

where min l and max l denote the l-th smallest and largest values, respectively, and
where l ≤ r/2. Then, for i = 0, . . . , r, the extreme rank of Ai with respect to A0 , . . . , Ar
is

n
l
l
Ri = max l : Alow
(m) ≤ Ai (m) ≤ Aupp
(m )

o
for all m ∈ [0, T ] .

The larger Ri is, the deeper or more central Ai is among A0 , . . . , Ar .

Now, for a given α ∈ (0, 1), the extreme rank ordering is used to define the 100(1 − α)%-

lα
lα
global rank envelope as the band delimited by the curves Alow
and Aupp
where
)
(
r
1
1 ( Ri < l ) ≤ α .
lα = max l :
r + 1 i∑
=0

Under H0 , with probability at least 1 − α,
lα
lα
Alow
(m) ≤ A0 (m) ≤ Aupp
(m )

for all m ∈ [0, T ],

(3)

see Myllymäki et al. (2016). Therefore, the 100(1 − α)%-global rank envelope is specifying a conservative statistical test called the extreme rank envelope test and which
accepts H0 at level 100α% if (3) is satisfied or equivalently if
r
1
1( Ri < R0 ) > α,
r + 1 i∑
=0

(4)

cf. Myllymäki et al. (2016). A plot of the extreme rank envelope allows a graphical interpretation of the extreme rank envelope test and may in case of rejection suggest an
alternative model for X0 .

13

There exist alternatives to the extreme rank envelope test, in particular a liberal extreme
rank envelope test and a so-called global scaled maximum absolute difference envelope, see Myllymäki et al. (2016). It is also possible to combine several extreme rank
envelopes, for instance by combining APF0 and APF1 , see Mrkvička et al. (2016). In the
following example we focus on (3)-(4) and briefly remark on results obtained by combining APF0 and APF1 .

Example 1 (simulation study).

Recall that a homogeneous Poisson process is a model

for complete spatial randomness (CSR), see e.g. Møller and Waagepetersen (2004) and
the simulation in the first panel of Figure 3. Consider APFs A0 , A1 , . . . , Ar corresponding
to independent point processes X0 , X1 , . . . , Xr defined on a unit square and where Xi for
i > 0 is CSR with a given intensity ρ (the mean number of points). Suppose X0 is claimed
to be CSR with intensity ρ, however, the model for X0 is given by one of the following
four point process models, which we refer to as the true model:
(a) CSR; hence the true model agrees with the claimed model.
(b) A Baddeley-Silverman cell process; this has the same second-order moment properties as under CSR, see Baddeley and Silverman (1984). Though from a mathematical point of view, it is a cluster process, simulated realisations will exhibit both
aggregation and regularity at different scales, see the second panel of Figure 3.
(c) A Matérn cluster process; this is a model for clustering where each cluster is a
homogenous Poisson process within a disc and the centers of the discs are not
observed and constitute a stationary Poisson process, see Matérn (1986),
Møller and Waagepetersen (2004), and the third panel of Figure 3.
(d) A most repulsive Bessel-type determinantal point process (DPP); this is a model for
regularity, see Lavancier et al. (2015), Biscio and Lavancier (2016), and the fourth
panel of Figure 3.
We let ρ = 100 or 400. This specifies completely the models in (a) and (d), whereas
14

Figure 3: Simulated point patterns for a homogeneous Poisson process (first panel), a
Baddeley-Silverman cell process (second panel), a Matérn cluster process (third panel),
and a most repulsive Bessel-type DPP (fourth panel).
the remaining parameters in the cases (b)-(c) are defined to be the same as those used
in Robins and Turner (2016). In all cases of Figure 3, ρ = 400. Finally, following the
recommendation in Myllymäki et al. (2016), we let r = 2499.
For each value of ρ = 100 or 400, we simulate each point process in (a)-(d) with the Rpackage spatstat. Then, for each dimension k = 0 or 1, we compute the extreme rank
envelopes and extreme rank envelope tests with the R-package spptest. We repeat
all this 500 times. Table 1 shows for each case (a)-(d) the percentage of rejection of the
hypothesis that X0 is a homogeneous Poisson process with known intensity ρ. In case of
CSR, the type one error of the test is small except when k = 0 and ρ = 100. As expected
in case of (b)-(d), the power of the test is increased when ρ is increased. For both the
Baddeley-Silverman process and the DPP, when k = 0 and/or ρ = 400, the power is
high and even 100% in two cases. For the Matérn cluster process, the power is 100%
when both ρ = 100 and 400; this is also the case when instead the radius of a cluster
becomes 10 times larger and hence it is not so easy to distinguish the clusters as in the
third panel of Figure 3. When we combine the extreme rank envelopes for APF0 and
APF1 , the results are better or close to the best results obtained when considering only
one extreme rank envelope.
Figure 4 illustrates for one of the 500 repetitions and for each dimension k = 0 and
k = 1 the deviation of APFk from the extreme rank envelope obtained when the true
model is not CSR. For each of the three non-CSR models, APFk is outside the extreme
15

CSR

DPP

Matérn cluster

Baddeley-Silverman

ρ = 100

ρ = 400

ρ = 100

ρ = 400

ρ = 100

ρ = 400

ρ = 100

ρ = 400

APF0

3.6

4

77.4

100

100

100

45.6

99.6

APF1

3.8

4.6

28.2

57.8

100

100

65.8

100

APF0 , APF1

4.8

3.6

82.4

100

100

100

60.8

100

Table 1: Percentage of point patterns for which the 95%-extreme rank envelope test rejects the hypothesis of CSR (a homogeneous Poisson process on the unit square with
intensity ρ = 100 or ρ = 400) when the true model is either CSR or one of three alternative point process models.
rank envelope, in particular when k = 0 and both the meanage and lifetime are small,
cf. the middle panel. This means that small lifetimes are not noise but of particular
importance, cf. the discussion in Section 1.2. Using an obvious notation, for small m,
MC
we may expect that APFDPP
(m) < APFCSR
0
0 (m ) < APF0 (m ) which is in agreement

with the middle panel. For large m, we may expect that APFDPP
(m) > APFCSR
0
0 (m ) and
MC
APFCSR
0 (m ) > APF0 (m ), but only the last relation is detected by the extreme rank
CSR
envelope in the left panel. Similarly, we may expect APFMC
1 (m ) > APF1 (m ) for small
CSR
m, whereas APFMC
1 (m ) < APF1 (m ) for large m, and both cases are detected in the

right panel. Note that for the Baddeley-Silverman cell process and k = 0, 1, APFBS
k
has a rather similar behaviour as APFDPP
, i.e. like a regular point process and probably
k
because clustering is a rare phenomena.
A similar simulation study is discussed in Robins and Turner (2016) for the models in
(a)-(c), but notice that they fix the number of points to be 100 and they use a testing procedure based on the persistent homology rank function, which in contrast to our onedimensional APF is a two-dimensional function and is not summarizing all the topological features represented in a persistent diagram. Robins and Turner (2016) show that
a test for CSR based on the persistent homology rank function is useful as compared
to various tests implemented in spatstat and which only concern first and secondorder moment properties. Their method is in particular useful, when the true model is
a Baddeley-Silverman cell process with the same first and second-order moment properties as under CSR. Comparing Figure 4 in Robins and Turner (2016) with the results
16

1.5

Determinantal

0.5

1.0

2.5
1.5

APF1(m)

Baddeley−Silverman

1.0

APF0(m)

Matern cluster

Baddeley−Silverman

0.5

Matern cluster

Extreme rank envelope bounds

Matern cluster
Determinantal

2.0

8
6
4

APF0(m)

2

Extreme rank envelope bounds

Extreme rank envelope bounds

0.1

0.2

0.3

0.4

0.0

0

Baddeley−Silverman

0.0

0.0

Determinantal

0.000

0.010

m

0.020

m

0.030

0.0

0.1

0.2

0.3

0.4

0.5

0.6

m

Figure 4: 95%-extreme rank envelope for APFi when i = 0 (left panel and the enlargement shown in the middle panel) or i = 1 (right panel) together with the curves for the
three non-CSR models (Baddeley-Silverman cell process, Matérn cluster process, and
Bessel-type DPP). The envelope is obtained from 2499 realisations of a CSR model on
the unit square and with intensity 100.
in Table 1 when the true model is a Baddeley-Silverman cell process and ρ = 100, the
extreme rank envelope test seems less powerful than the test they suggest. On the other
hand, Robins and Turner (2016) observe that the latter test performs poorly when the
true model is a Strauss process (a model for inhibition) or a Matérn cluster process; as
noticed for the Matérn cluster process, we obtain a perfect power when using the extreme rang envelope test.

4 A single sample of accumulated persistence functions
4.1 Functional boxplot
This section discusses the use of a functional boxplot (Sun and Genton, 2011) for a sample A1 , . . . , Ar of APFk s those joint distribution is exchangeable. The plot provides a
representation of the variation of the curves given by A1 , . . . , Ar around the most central curve, and it can be used for outlier detection, i.e. detection of curves that are too
extreme with respect to the others in the sample. This is illustrated in Example 2 for the

17

brain artery trees dataset and in Appendix B and its accompanying Example 6 concerning a simulation study.
The functional boxplot is based on an ordering of the APFk s obtained using a so-called
depth function. For specificity we make the standard choice called the modified band
depth function (MBD), cf. López-Pintado and Romo (2009) and Sun and Genton (2011):
For a user-specified parameter T > 0 and h, i, j = 1, . . . , r with i < j, define



Bh,i,j = m ∈ [0, T ] : min Ai (m) , A j (m) ≤ Ah (m) ≤ max Ai (m) , A j (m)

,

and denote the Lebesgue measure on [0, T ] by |·|. Then the MBD of Ah with respect to
A1 , . . . , Ar is
MBDr ( Ah ) =

2
Bh,i,j .
r (r − 1) 1≤∑
i < j ≤r

(5)

This is the average proportion of Ah on [0, T ] between all possible pairs of A1 , . . . , Ar .
Thus, the larger the value of the MBD of a curve is, the more central or deeper it is in
the sample. We call the region delimited by the 50% most central curves the central
envelope. It is often assumed that a curve outside the central envelope inflated by 1.5
times the range of the central envelope is an outlier or abnormal curve — this is just
a generalisation of a similar criterion for the boxplot of a sample of real numbers —
and the range may be changed if it is more suitable for the application at hand, see the
discussion in Sun and Genton (2011) and Example 2 below.

Example 2 (brain artery trees). For the brain artery trees dataset (Section 2.2), Figure 5
shows the functional boxplots of APFk s for females (first and third panels) respective
males (second and fourth panels) when k = 0 (first and second panels) and k = 1 (third
and fourth panels): The most central curve is plotted in black, the central envelope in
purple, and the upper and lower bounds obtained from all the curves except the outliers
in dark blue. Comparing the two left panels (concerned with connected components),
the shape of the central envelope is clearly different for females and males, in particular
on the interval [40, 60], and the upper and lower bounds of the non-outlier are closer
18

20

40

60

80 100

m

0

20

40

60

80 100

1000
200
0

0

m

600

APF1(m) − Male

1000
600

APF1(m) − Female

0

200

1500
1000
0

500

APF0(m) − Male

1500
1000
500

APF0(m) − Female

0
0

5

10

15

20

25

0

5

10

m

15

20

25

m

Figure 5: Functional boxplots of APFs for females and males obtained from the brain
artery trees dataset: APF0F (first panel), APF0M (second panel), APF1F (third panel), APF1M
(fourth panel). The dashed lines show the outliers detected by the 1.5 criterion.

Figure 6: Brain artery tree of a male subject with APF0M and APF1M detected as outliers
by the 1.5 criterion.
to the central region for females, in particular on the interval [0, 50]. For the two right
panels (concerned with loops), the main difference is observed on the interval [15, 25]
where the central envelope is larger for females than for males.
The dashed lines in Figure 5 show the APFs detected as outliers by the 1.5 criterion, that
is 6 APF0F s (first panel), 3 APF1F s (third panel), 6 APF0M s (second panel), and 4 APF1M s
(fourth panel). For the females, only for one point pattern both APF0F and APF1F are
outliers, where APF1F is the steep dashed line in the bottom-left panel; and for the males,
only for two point patterns both APF0M and APF1M are outliers, where in one case APF1M
is the steep dashed line in the bottom-right panel. For this case, Figure 6 reveals an
obvious issue: A large part on the right of the corresponding tree is missing!
Examples 3 and 4 discuss to what extent our analysis of the brain artery trees will be
sensitive to whether we include or exclude the detected outliers.
19

4.2 Confidence region for the mean function
This section considers an asymptotic confidence region for the mean function of a sample A1 , . . . , Ar of IID APFk s. We assume that D1 , . . . , Dr are the underlying IID RRPDk s
for the sample so that with probability one, there exists an upper bound T < ∞ on the
death times and there exists an upper bound nmax < ∞ on the number of k-dimensional
topological features. Note that the state space for such RRPDk s is
n

Dk,T,nmax = {{(m1 , l1 , c1 ), . . . , (mn , ln , cn )} :

∑ ci ≤ nmax , mi + li /2 ≤ T, i = 1, . . . , n}

i =1

and only the existence and not the actual values of nmax and T play a role when applying our method below. For example, in the settings (i)-(ii) of Section 2.1 it suffices
to assume that X is included in a bounded region of R2 and that the number of points
N is bounded by a constant; this follows from the two versions of the Nerve Theorem
presented in Fasy et al. (2014) and Edelsbrunner and Harer (2010), respectively.
We adapt an empirical bootstrap procedure (see e.g. van der Vaart and Wellner (1996))
which in Chazal et al. (2013) is used for a confidence region for the mean of the dominant function of the persistent landscape and which in our case works as follows. For
0 ≤ m ≤ T, the mean function is given by µ(m) = E { A1 (m)} and estimated by the

empirical mean function Ar (m) =

1
r

∑ri=1 Ai (m). Let A1∗ , . . . , Ar∗ be independent uni-

form draws with replacement from the set { A1 , . . . , Ar } and set Ar∗ = 1r ∑ri=1 Ai∗ and
√
θ ∗ = supm∈[0,T ] r Ar (m) − Ar∗ (m) . For a given integer B > 0, independently repeat
this procedure B times to obtain θ1∗ , . . . , θB∗ . Then, for 0 < α < 1, the 100(1 − α)%quantile in the distribution of θ ∗ is estimated by
q̂αB = inf{q ≥ 0 :

1 B
1(θi∗ > q) ≤ α}.
B i∑
=1

The following theorem is verified in Appendix F.
Theorem 4.1. Let the situation be as described above. For large values of r and B, the functions
√
Ar ± q̂αB / r provide the bounds for an asymptotic conservative 100(1 − α)%-confidence region
20

20

40

60

80 100

400
200

APF1(m)
0

Confidence region
bounds − Female

0

0

Confidence region
bounds − Male

600

800

1000

1400
1000
600

APF0(m)

200

Confidence region
bounds − Female

Confidence region
bounds − Male

0

m

5

10

15

20

25

m

Figure 7: Bootstrap confidence regions for the mean APFkM and the mean APFkF when
k = 0 (left panel) and k = 1 (right panel).
for the mean APF, that is


√
√
lim lim P µ(m) ∈ [ Ar (m) − q̂αB / r, Ar (m) + q̂αB / r ] for all m ∈ [0, T ] ≥ 1 − α.

r →∞ B→∞

Example 3 (brain artery trees).

The brain artery trees are all contained in a bounded

region and presented by a bounded number of points, so it is obvious that T and nmax
exist for k = 0, 1. To establish confidence regions for the mean of the APFkM s respective
APFkF s, we apply the bootstrap procedure with B = 1000. The result is shown in Figure 7
when all 95 trees are considered: In the left panel, k = 0 and approximatively half of
each confidence region overlap with the other confidence region; it is not clear if there
is a difference between genders. In the right panel, k = 1 and the difference is more
pronounced, in particular on the interval [15, 25]. Similar results and conclusions are
obtained if we exclude the APFs detected as outliers in Example 2. Of course we should
supply with a statistical test to assess the gender effect and such a test is established in
Section 5 and applied in Example 4.
Appendix C provides the additional Example 7 for a simulated dataset along with a
discussion on the geometrical interpretation of the confidence region obtained.

21

5 Two samples of accumulated persistence functions
This section concerns a two-sample test for comparison of two samples of APFs. Appendix E presents both a clustering method (Appendix E.1, including Example 9) and
a unsupervised classification method (Appendix E.2, including Example 10) for two or
more samples.
Consider two samples of independent RRPDk s D1 , . . . , Dr1 and E1 , . . . , Er2 , where each
Di (i = 1, . . . , r1 ) has distribution PD and each E j has distribution PE (j = 1, . . . , r2 ),
and suppose we want to test the null hypothesis H0 : PD = PE = P. Here, the common

distribution P is unknown and as in Section 4.2 we assume it is concentrated on Dk,T,nmax

for some integer nmax > 0 and number T > 0. Below, we adapt a two-sample test
statistic studied in Præstgaard (1995) and van der Vaart and Wellner (1996).
Let r = r1 + r2 . Let A1 , . . . , Ar be the APFk corresponding to ( D1 , . . . , Dr1 , E1 , . . . , Er2 ),
and denote by Ar1 and Ar2 the empirical means of A1 , . . . , Ar1 and Ar1 +1 , . . . , Ar1 +r2 ,
respectively. Let I = [ T1 , T2 ] be a user-specified interval with 0 ≤ T1 < T2 ≤ T and used
for defining a two-sample test statistic by
r
r1 r2
KSr1 ,r2 =
sup Ar1 (m) − Ar2 (m) ,
r m∈ I
where large values are critical for H0 . This may be rewritten as
r
r
r
r 2 r1
r 1 r2
r1 r2
KSr1 ,r2 = sup
GD ( m ) −
GE ( m ) +
E { A D − A E } (m ) ,
r
r
r
m∈ I
r1
=
where GD

√

(6)

(7)



√
r1 Ar1 − E { A D } and GEr2 = r2 Ar2 − E { A E } . By Lemma F.2 in Apr

pendix F and by the independence of the samples, GD1 and GEr2 converge in distribution

to two independent zero-mean Gaussian processes on I, denoted GD and GE , respectively. Assume that r1 /r → λ ∈ (0, 1) as r → ∞. Under H0 , in the sense of convergence
in distribution,
lim KSr1 ,r2 = sup

r →∞

m∈ I

√

1 − λGD (m) −
22

√

λGE (m) ,

(8)

where

√
1 − λGD − λGE follows the same distribution as GD . If H0 is not true and

E A1 − Ar1 +1 (m) > 0, then KSr1 ,r2 → ∞ as r → ∞, see van der Vaart and Wellner

√

supm∈ I

(1996). Therefore, for 0 < α < 1 and letting qα = inf{q : P(supm∈ I | GD (m)| > q) ≤ α},
the asymptotic test that rejects H0 if KSr1 ,r2 ≤ qα is of level 100α% and of power 100%.

As qα depends on the unknown distribution P, we estimate qα by a bootstrap method:
Let A1∗ , . . . , Ar∗ be independent uniform draws with replacement from { A1 , . . . , Ar }. For

0 ≤ m ≤ T, define the empirical mean functions Ar∗1 (m) =
1
r2

1
r1

r1
∗
∗
∑i =
1 Ai (m ) and Ar2 (m ) =

r +r

2
∗
1
∑i =
r1 +1 Ai (m ), and compute

∗

θ =

r

r1 r2
sup A∗ (m) − Ar∗2 (m) .
r m ∈ I r1

(9)

For a given integer B > 0, independently repeat this procedure B times to obtain θ1∗ , . . . , θB∗ .
Then we estimate qα by the 100(1 − α)%-quantile of the empirical distribution of θ1∗ , . . . , θB∗ ,
that is
q̂αB

1 B
= inf{q ≥ 0 : ∑ 1(θi∗ > q) ≤ α}.
B i =1

The next theorem is a direct application of Theorem 3.7.7 in van der Vaart and Wellner
(1996) noticing that the APFk s are uniformly bounded by Tnmax and they form a socalled Donsker class, see Lemma F.2 and its proof in Appendix F.
Theorem 5.1. Let the situation be as described above. If r → ∞ such that r1 /r → λ with
λ ∈ (0, 1), then under H0



lim lim P KSr1 ,r2 >

r →∞ B→∞

q̂αB



= α,


whilst if H0 is not true and supm∈ I E A1 − Ar1 +1 (m) > 0, then


lim lim P KSr1 ,r2 > q̂αB = 1.

r →∞ B→∞

Therefore, the test that rejects H0 if KSr1 ,r2 > q̂αB is of asymptotic level 100α% and power
100%. As remarked in van der Vaart and Wellner (1996), by their Theorem 3.7.2 it is
23

possible to present a permutation two-sample test so that the critical value q̂αB for the
bootstrap two-sample test has the same asymptotic properties as the critical value for
the permutation test.
Other two-sample test statistics than (6) can be constructed by considering other measurable functions of Ar1 − Ar2 , e.g. we may consider the two-sample test statistic
Mr1 ,r2 =

Z

I

Ar1 (m) − Ar2 (m) dm.

(10)

Then by similar arguments as above but redefining θ ∗ in (9) by
r
Z
r1 r2
∗
A∗ (m) − Ar∗2 (m) dm,
θ =
r m ∈ I r1
the test that rejects H0 if Mr1 ,r2 > q̂αB is of asymptotic level 100α% and power 100%.

Example 4 (brain artery trees).

To distinguish between male and female subjects of

the brain artery trees dataset, we use the two-sample test statistic KSr1 ,r2 under three
different settings:
(A) For k = 0, 1, we let PD′k be the subset of PDk corresponding to the 100 largest
lifetimes. Then D1 , . . . , D46 and E1 , . . . , E49 are the RRPDk s obtained from the PD′k s
associated to female and male subjects, respectively. This is the setting used in
Bendich et al. (2016).
(B) For k = 0, 1, we consider all lifetimes and let D1 , . . . , D46 and E1 , . . . , E49 be the
RRPDk s associated to female and male subjects, respectively.
(C) The samples are as in setting (B) except that we exclude the RRPDk s where the
corresponding APFk was detected as an outlier in Example 2. Hence, r1 = 40 and
r2 = 43 if k = 0, and r1 = 43 and r2 = 45 if k = 1.
Bendich et al. (2016) perform a permutation test based on the mean lifetimes for the male
and female subjects and conclude that gender effect is recognized when considering PD1
24

APF0

APF1

I = [0, 137]

I = [0, 60]

I = [0, 25]

I = [15, 25]

Setting (A)

5.26

3.26

3.18

2.72

Setting (B)

7.67

3.64

20.06

1.83

Setting (C)

4.55

2.61

0.92

0.85

Table 2: Estimated p-values given in percentage of the two-sample test based on KSr1 ,r2
used with APF0 and APF1 on different intervals I to distinguish between male and female subjects under settings (A), (B), and (C) described in Example 4.
(p-value = 3%) but not PD0 (p-value = 10%). For comparison, under each setting (A)(C), we perform the two-sample test for k = 0, 1, different intervals I, and B = 10000.
In each case, we estimate the p-value, i.e. the smallest α such that the two-sample test
with significance level 100α% does not reject H0 , by p̂ =

1
B

∑iB=1 1(θi∗ > KSr1 ,r2 ). Table 2

shows the results. Under each setting (A)-(C), using APF0 we have a smaller p-value
than in Bendich et al. (2016) if I = [0, 137] and an even larger p-value if I = [0, 60];
and for k = 1 under setting (B), our p-value is about seven times larger than the pvalue in Bendich et al. (2016) if I = [0, 25], and else it is similar or smaller. For k = 1
and I = [0, 25], the large difference between our p-values under settings (B) and (C)
indicates that the presence of outliers violates the result of Theorem 5.1 and care should
hence be taken. In our opinion we can better trust the results without outliers, where
in contrast to Bendich et al. (2016) we see a clear gender effect when considering the
connected components. Notice also that in agreement with the discussion of Figure 5
in Example 2, for each setting A, B, and C and each dimension k = 0, 1, the p-values in
Table 2 are smallest when considering the smaller interval I = [0, 60] or I = [15, 25].
Appendix D provides an additional Example 8 illustrating the use of two-sample test in
a simulation study.

25

Acknowledgements

Supported by The Danish Council for Independent Research | Natural Sciences, grant
7014-00074B, "Statistics for point processes in space and beyond", and by the "Centre for
Stochastic Geometry and Advanced Bioimaging", funded by grant 8721 from the Villum Foundation. Helpful discussions with Lisbeth Fajstrup on persistence homology is
acknowledged. In connection to the brain artery trees dataset we thank James Stephen
Marron and Sean Skwerer for helpful discussions and the CASILab at The University
of North Carolina at Chapel Hill for providing the data distributed by the MIDAS Data
Server at Kitware, Inc. We are grateful to the editors and the referees for useful comments.

Appendix
Appendix A-F contain complements and additional examples to Sections 3-5. Our setting and notation are as follows. All the examples are based on a simulated point pattern

{ x1 , . . . , x N } ⊂ R2 as described in Section 2.1, with x1 , . . . , x N being IID points where
N is a fixed positive integer. As in Section 1.1.1, our setting corresponds to applications
typically considered in TDA where the aim is to obtain topological information about a
compact set C ⊂ R2 which is unobserved and where possibly noise appears: For specificity, we let xi = yi + ǫi , i = 1, . . . , N, where y1 , . . . , y N are IID points with support
C, the noise ǫ1 , . . . , ǫ N are IID and independent of y1 , . . . , y N , and ǫi follows the restriction to the square [−10σ, 10σ]2 of a bivariate zero-mean normal distribution with IID
coordinates and standard deviation σ ≥ 0 (if σ = 0 there is no noise). We denote this

distribution for ǫi by N2 (σ) (the restriction to [−10σ, 10σ]2 is only imposed for technical

reasons and is not of practical importance). We let Ct be the union of closed discs of radii
t and centred at x1 , . . . , x N , and we study how the topological features of Ct changes as
t ≥ 0 grows. For this we use the Delaunay-complex mentioned in Section 1.1.1. Finally,

we denote by C((a, b), r ) the circle with center (a, b) and radius r.
26

A Transforming confidence regions for persistence diagrams
used for separating topological signal from noise
As noted in Section 3 there exists several constructions and results on confidence sets
for persistence diagrams when the aim is to separate topological signal from noise,
see Fasy et al. (2014), Chazal et al. (2014), and the references therein. We avoid presenting
the technical description of these constructions and results, which depend on different
choices of complexes (or more precisely so-called filtrations). For specificity, in this appendix we just consider the Delaunay-complex and discuss the transformation of such
a confidence region into one for an accumulate persistence function.
We use the following notation. As in the aforementioned references, consider the persistence diagram PDk for an unobserved compact manifold C ⊂ R2 and obtained as in

Section 1.1.1 by considering the persistence as t ≥ 0 grows of k-dimensional topological
features of the set consisting of all points in R2 within distance t from C. Note that PDk

is considered as being non-random and unknown; of course in our simulation study
c k,N be the
presented in Example 5 below we only pretend that PDk is unknown. Let PD

random persistence diagram obtained as in Section 1.1.1 from IID points x1 , . . . , x N with
√
support C. Let N = {(b, d) : b ≤ d, l ≤ 2c N } be the set of points at distance 2c N of the

diagonal in the persistence diagram. Let S(b, d) = {( x, y) : | x − b| ≤ c N , |y − d| ≤ c N }

be the square with center (b, d), sides parallel to the b- and d-axes, and of side length
2c N . Finally, let α ∈ (0, 1).
Fasy et al. (2014) and Chazal et al. (2014) suggest various ways of constructing a bound
c N > 0 so that an asymptotic conservative 100(1 − α)%-confidence region for PDk with
respect to the bottleneck distance W∞ is given by


c
lim inf P W∞ (PDk , PDk,N ) ≤ c N ≥ 1 − α
N →∞

(11)

where W∞ is the bottleneck distance defined in Section 1.3. The confidence region given
by (11) consists of those persistence diagrams PDk which have exactly one point in each

c k,N , and have an arbitrary number of points in
square S(bi , di ), with (bi , di ) a point of PD
27

c k,N falling in N as noise and the
the set N . Fasy et al. (2014) consider the points of PD
remaining points as representing a significant topological feature of C.

Using (11) an asymptotic conservative 100(1 − α)%-confidence region for the APFk corresponding to PDk is immediately obtained. This region will be bounded by two funcbmin and A
bmax specified by PD
c k,N and c N . Due to the accumulating nature of APFk ,
tions A
k,N
k,N

the span between the bounds is an increasing function of the meanage. When using the

Delaunay-complex, Chazal et al. (2014) show that the span decreases as N increases; this
is illustrated in Example 5 below.

Example 5 (simulation study).

Let C = C((−1.5, 0), 1) ∪ C((1.5, 0), 0.8) and suppose

each point xi is uniformly distributed on C. Figure 8 shows C and an example of a simulated point pattern with N = 300 points. We use the bootstrap method implemented in
the R-package TDA and presented in Chazal et al. (2014) to compute the 95%-confidence
region for PD1 when N = 300, see the top-left panel of Figure 9, where the two squares
above the diagonal correspond to the two loops in C and the other squares correspond
to topological noise. Thereby 95%-confidence regions for RRPD1 (top-right panel) and
APF1 (bottom-left panel) are obtained. The confidence region for APF1 decreases as N
increases as demonstrated in the bottom panels where N is increased from 300 to 500.
As noticed in Section 1.3, we must be careful when using results based on the bottleneck
metric, because small values of the bottleneck metric does not correspond to closeness
of the two corresponding APFs: Although close persistence diagram with respect to the
bottleneck distance imply that the two corresponding APFs are close with respect to the
Lq -norm (1 ≤ q ≤ ∞), the converse is not true. Hence, it is possible that an APF is in the
confidence region plotted in Figure 9 but that the corresponding persistence diagram is
very different from the truth.

28

Figure 8: The set C in Example 5 (left panel) and a simulated point pattern of N = 300
independent and uniformly distributed points on C (right panel).

B Additional example related to Section 4.1 "Functional
boxplot"
The functional boxplot described in Section 4.1 can be used as an exploratory tool for the
curves given by a sample of APFk s. It provides a representation of the most central curve
and the variation around this. It can also be used for outliers detection as illustrated in
the following example.

Example 6 (simulation study).

We consider a sample of 65 independent APFk s, where

the joint distribution of the first 50 APFk s is exchangeable, whereas the last 15 play the
role of outliers. We suppose each APFk corresponds to a point process of 100 IID points,
where each point xi follows one of the following distributions P1 , . . . , P4 .
• P1 (unit circle): xi is a uniform point on C((0, 0), 1) perturbed by N2 (0.1)-noise.
• P2 (Gaussian mixture): Let yi follow N2 (0.2), then xi = yi with probability 0.5, and
xi = yi + (1.5, 0.5) otherwise.
• P3 (two circles): xi is a uniform point on C((−1, −1), 1) ∪ C((1, 1), 0.5) perturbed
by N2 ((0, 0), 0.1)-noise.

• P4 (circle of radius 0.7): xi is a uniform point on C((0, 0), 0.7) perturbed by N2 (0.1)noise.

29

1.2
0.6
0.2

0.4

Lifetime

0.8

1.0

1.0
0.8
0.6

Death time

0.4

0.0

0.2
0.0
0.0

0.2

0.4

0.6

0.8

1.0

−0.2

0.4

0.6

0.8

5

APF1

0

0

1

2

3
2

3

4

Confidence region bounds

APF1

APF1(m)

4

Confidence region bounds

1

APF1(m)

0.2

Meanage

5

Birth time

0.0

0.5

1.0

1.5

0.0

0.5

m

1.0

1.5

m

Figure 9: 95%-confidence regions obtained by the bootstrap method for PD1 (top-left
panel) and its corresponding RRPD1 (top-right panel) when C and x1 , . . . , x300 are as
in Figure 8. The bottom-left panel shows the corresponding 95%-confidence region for
APF1 . The bottom-right panel shows the 95%-confidence region for APF1 when a larger
point cloud with 300 points is used.
We let the first 50 point processes be obtained from P1 (the distribution for non-outliers),
the next 5 from P2 , the following 5 from P3 , and the final 5 from P4 . Figure 10 shows a
simulated realization of each of the four types of point processes.
Figure 11 shows the functional boxplots when considering APF0 (left panel) and APF1
(right panel). The curves detected as outliers and corresponding to the distributions
P2 , P3 , and P4 are plotted in red, blue, and green, respectively. In both panels the outliers
detected by the 1.5 criterion agree with the true outliers.
In the left panel, each curve has an accumulation of small jumps between m = 0 and
m ≈ 0.1, corresponding to the moments where the points associated to each circle are

connected by the growing discs in the sequence {Ct }t≥0 . The curves corresponding to

realisations of P3 have a jump at m ≈ 0.38 which corresponds to the moment where
30

P1 −
P2 −
P3 −
P4 −

Circle r=1
Gaussian mixture
Two circles
Circle r=0.7

Figure 10: Simulated realizations of the four types of point processes, each consisting of
100 IID points with distribution either P1 (black dots), P2 (blue crosses), P3 (red triangles), or P4 (rotated green crosses).
the points associated to the two circles used when defining P3 are connected by the
growing discs in the sequence {Ct }t≥0 . The points following the distribution P4 are
generally closer to each other than the ones following the distribution P1 as the radius
of the underlying circle is smaller. This corresponds to more but smaller jumps in APF0
for small meanages, and hence the curves of APF0 are lower when they correspond to
realisations of P1 than to realisations of P4 ; and as expected, for large meanages, the
curves of APF0 are larger when they correspond to realisations of P1 than to realisations
of P4 . Note that if we redefine P4 so that the N2 (0.1)-noise is replaced by N2 (0.07)-noise,
then the curves would be the same up to rescaling.
In the right panel, we observe clear jumps in all APF1 s obtained from P1 , P3 , and P4 .
These jumps correspond to the first time that the loops of the circles in P1 , P3 , and P4 are
covered by the union of growing discs in the sequence {Ct }t≥0 . Once again, if we have

used N2 (0.07)-noise in place of N2 (0.1)-noise in the definition of P4 , the curves would
be the same up to rescaling.

If we repeat everything but with the distribution P4 redefined so that C((0, 0), 0.7) is

replaced by C((0, 0), 0.8), then the support of P4 is closer to that of P1 and it becomes
harder in the case of APF0 to detect the outliers with distribution P4 (we omit the corresponding plot); thus further simulations for determining a stronger criterion would be
31

1.5
1.0

APF1(m)

0.5

6
4

Gaussian mixture
Two circles
Circle r=0.7

0

0.0

2

APF0(m)

8

Gaussian mixture
Two circles
Circle r=0.7

0.0

0.2

0.4

0.6

0.8

0.0

m

0.5

1.0

1.5

m

Figure 11: Functional boxplots of 65 APFs based on the topological features of dimension 0 (left panel) and 1 (right panel). In each panel, 50, 5, 5, and 5 APFs are obtained
from the Delaunay-complex of 100 IID points from the distribution P1 , P2 , P3 , and P4 ,
respectively. The APFs detected as outliers are plotted in red, blue, and green in the case
of P2 , P3 , and P4 , respectively.
needed.

C Additional example related to Section 4.2 "Confidence
region for the mean function"
This appendix provides yet an example to illustrate the bootstrap method in Section 4.2
for obtaining a confidence region for the mean function of a sample of IID APFk s.

Example 7 (simulation study).

Consider 50 IID copies of a point process consisting

of 100 independent and uniformly distributed points on the union of three circles with
radius 0.25 and centred at (−1, −1), (0, 1), and (1, −1), respectively (these circles were
also considered in the example of Section 1.1.1). A simulated realization of the point
process is shown in the left panel of Figure 12, and the next two panels show simulated
confidence regions for APF0 and APF1 , respectively, when the bootstrap procedure with
B = 1000 is used. In the middle panel, between m = 0 and m = 0.2, there is an accu32

1.2
0.1

0.2

0.3

0.4

m

0.5

0.6
0.4
0.2

APF1(m)
0.0

Confidence
region bounds

0.0

0

Empirical mean

0.8

1.0

6
5
4
3

APF0(m)

2
1

Confidence
region bounds

Empirical mean

0.0

0.2

0.4

0.6

0.8

1.0

m

Figure 12: A simulation of 100 independent and uniformly distributed points on the
union of three circles (dashed lines) with the same radius r = 0.5 and centred at

(−1, −1), (0, 1), and (1, −1) (left panel). The 95%-confidence regions for the mean APF0
(middle panel) and the mean APF1 (right panel) are based on 50 IID simulations.
mulation of small jumps corresponding to the moment when each circle is covered by
the union of growing discs from the sequence {Ct }t≥0 ; we interpret these small jumps
as topological noise. The jump at m ≈ 0.25 corresponds to the moment when the circles

centred at (−1, −1) and (1, −1) are connected by the growing discs, and the jump at

m ≈ 0.3 to when all three circles are connected by the growing discs. In the right panel,

at m ≈ 0.3 there is an accumulation of small jumps corresponding to the moment when
the three circles are connected by the growing discs and they form a loop at m = 0.25
in Figure 1. The disappearance of this loop at m = 0.69 in Figure 1 corresponds to the
jump at m ≈ 0.7 in Figure 12.

D

Additional example to Section 5 "Two samples of accumulated persistence functions"

Section 5 considered two samples of independent RRPDk s D1 , . . . , Dr1 and E1 , . . . , Er2 ,
where each Di (i = 1, . . . , r1 ) has distribution PD and each E j has distribution PE (j =
1, . . . , r2 ). Then we studied a bootstrap two-sample test to asses the null hypothesis H0 :
PD = PE , e.g. in connection to the brain artery trees. An additional example showing
33

Circle r=1
Circle r=0.95

Figure 13: A simulation of 100 independent and uniformly distributed points on the
circle centred at (0, 0) with radius 1 and perturbed by N2 ((0, 0), 0.2)-noise (red dots),
together with 100 independent and uniformly distributed points on the circle centred at

(0, 0) with radius 0.95 and perturbed by N2 ((0, 0), 0.2)-noise (blue crosses).
the performance of the test is presented below.

Example 8 (simulation study).

Let PD be the distribution of a RRPDk obtained from

100 independent and uniformly distributed points on C((0, 0), 1) perturbed by N2 (0.2)noise, and define PE in a similar way but with a circle of radius 0.95. A simulated
realisation of each point process is shown in Figure 13; it seems difficult to recognize
that the underlying circles are different. Let us consider the two-sample test statistics
(6) and (10) with I = [0, 3], r1 = r2 = 50, and α = 0.05. Over 500 simulations of the
two samples of RRPDk we obtain the following percentage of rejection: For Mr1 ,r2 , 5.2%
if k = 0, and 24.2% if k = 1. For KSr1 ,r2 much better results are observed, namely 73.8%
if k = 0, and 93.8% if k = 1, where this high percentage is mainly caused by the largest
lifetime of a loop.

34

E Further methods for two or more samples of accumulated persistence functions
E.1 Clustering
Suppose A1 , . . . , Ar are APFk s which we want to label into K < r groups by using a
method of clustering (or unsupervised classification). Such methods are studied many
places in the literature for functional data, see the survey in Jacques and Preda (2014). In
particular, Chazal et al. (2009), Chen et al. (2015), and Robins and Turner (2016) consider
clustering in connection to RRPDk s. Whereas the RRPDk s are two-dimensional functions, it becomes easy to use clustering for the one-dimensional APFk s as illustrated in
Example 9 below.
For simplicity we just consider the standard technique known as the K-means clustering
algorithm (Hartigan and Wong (1979)). For more complicated applications than considered in Example 9 the EM-algorithm may be needed for the K-means clustering algorithm. As noticed by a referee, to avoid the use of the EM-algorithm we can modify (9)
or (10) and thereby construct a distance/similarity matrix for different APFs which is
used to perform hierarchical clustering. However, for Example 9 the results using hierarchical clustering (omitted here) were not better than with the K-means algorithm.
Assume that A1 , . . . , Ar are pairwise different and square-integrable functions on [0, T ],
where T is a user-specified parameter. For example, if RRPDk ∈ Dk,T,nmax (see Section 4.2), then APFk ∈ L2 ([0, T ]). The K-means clustering algorithm works as follows.

• Chose uniformly at random a subset of K functions from { A1 , . . . , Ar }; call these
functions centres and label them by 1, . . . , K.
• Assign each non-selected APFk the label i if it is closer to the centre of label i than
to any other centre with respect to the L2 -distance on L2 ([0, T ]).

35

• In each group, reassign the centre by the mean curve of the group (this may not be
an APFk of the sample).
• Iterate these steps until the assignment of centres does not change.
The algorithm is known to be convergent, however, it may have several drawbacks as
discussed in Hartigan and Wong (1979) and Bottou and Bengio (1995).

Example 9 (simulation study).

Consider K = 3 groups, each consisting of 50 APF0 s

and associated to point processes consisting of 100 IID points, where each point xi follows one of the following distributions P1 , P2 , and P3 for groups 1, 2, and 3, respectively.
• P1 (unit circle): xi is a uniform point on C((0, 0), 1) perturbed by N2 (0.1)-noise.
• P2 (two circles): xi is a uniform point on C((−1, −1), 0.5) ∪ C((1, 1), 0.5) perturbed
by N2 (0.1)-noise.

• P3 (circle of radius 0.8): xi is a uniform point on C((0, 0), 0.8) perturbed by N2 (0.1)noise.
We start by simulating a realization of each of the 3 × 50 = 150 point processes. The left
panel of Figure 14 shows one realization of each type of point process; it seems difficult
to distinguish the underlying circles for groups 1 and 3, but the three APF0 s associated
to these three point patterns are in fact assigned to their right groups. The right panel
of Figure 14 shows the result of the K-means clustering algorithm. Here we are using
the R-function “kmeans” for the K-means algorithm and it takes only a few seconds
when evaluating each Ai (m) at 2500 equidistant values of m between 0 and T = 0.5.
As expected we see more overlap between the curves of the APF0 s assigned to groups 1
and 3.
We next repeat 500 times the simulation of the 150 point processes. A clear distinction
between the groups is obtained by the K-means algorithm applied for connected com-

36

ponents: The percentage of wrongly assigned APF0 s among the 500 × 3 × 50 = 75000
APF0 s has an average of 4.5% and a standard deviation of 1.6%. The assignment error
is in fact mostly caused by incorrect labelling of APF0 s associated to P1 or P3 . This is
expected as the underlying circles used in the definitions of P1 and P3 are rather close,
whereas the underlying set in the definition of P2 is different with two connected components as represented by the jump at m ≈ 0.4 in the middle panel of Figure 14.
Even better results are obtained when considering loops instead of connected components: The percentage of wrongly assigned APF1 s among the 75000 APF1 s has an average of 1.6% and a standard deviation of 1.0%. This is mainly due to the sets underlying
P1 , P2 , and P3 which have distinctive loops that results in clear distinct jumps in the
APF1 s as seen in the right panel of Figure 14.

0.2

0.4

m

0.8
0.4

APF1(m)
0.0

0.0

4
3
2

P1 − Circle r=1
P2 − Two circles
P3 − Circle r=0.8

1

P1 − Circle r=1
P2 − Two circles
P3 − Circle r=0.8

0

APF0(m)

5

6

1.2

P1 − Circle r=1
P2 − Two circles
P3 − Circle r=0.8

0.0

0.2

0.4

m

Figure 14: Left panel: Simulated example of the three point processes, each consisting of
100 IID points drawn from the distribution P1 (black dots), P2 (red triangles), or P3 (blue
crosses). Middle panel: The 150 APF0 s obtained from the simulation of the 150 point
processes associated to P1 , P2 , or P3 , where the colouring in black, red, or blue specifies
whether the K-means algorithm assigns an APF0 to the group associated to P1 , P2 , or P3 .
Right panel: As the middle panel but for the 150 APF1 s.

E.2 Supervised classification
Suppose we want to assign an APFk to a training set of K different groups G1 , . . . , GK ,

where Gi is a sample of ri independent APFk s Ai1 , . . . , Airi . For this purpose supervised
classification methods for functional data may be adapted.
37

We just consider a particular method by López-Pintado et al. (2010): Suppose α ∈ [0, 1]

and we believe that at least 100(1 − α)% of the APFk s in each group are IID, whereas
the remaining APFk s in each group follow a different distribution and are considered as
outliers (see Section 4.1). For a user-specified parameter T > 0 and i = 1, . . . , K, define
α

the 100α%-trimmed mean Ai with respect to Gi as the mean function on [0, T ] of the

100(1 − α)% APFk s in Gi with the largest MBDri , see (5). Assuming ∪iK=1 Gi ⊂ L2 ([0, T ]),

an APFk A ∈ L2 ([0, T ]) is assigned to Gi if

α

i = argmin k A j − Ak,

(12)

j∈{1,...,K }

where k · k denotes the L2 -distance. Here, the trimmed mean is used for robustness and
allows a control over the curves we may like to omit because of outliers, but e.g. the
median could have been used instead.

Example 10 (simulation study).

Consider the following distributions P1 , . . . , P4 for a

point xi .
• P1 (unit circle): xi is a uniform point on C((0, 0), 1) which is perturbed by N2 (0.1)noise.
• P1′ (two circles, radii 1 and 0.5): xi is a uniform point on C((0, 0), 1) ∪ C((1.5, 1.5), 0.5)
and perturbed by N2 (0.1)-noise.

• P2 (circle of radius 0.8): xi is a uniform point on C((0, 0), 0.8) which is perturbed
by N2 (0.1)-noise.

• P2′ (two circles, radii 0.8 and 0.5): xi is a uniform point on C((0, 0), 0.8) ∪ C((1.5, 1.5), 0.5)
and perturbed by N2 (0.1)-noise.

For k = 0, 1 we consider the following simulation study: K = 2 and r1 = r2 = 50; G1
consists of 45 APFk s associated to simulations of point processes consisting of 100 IID
points with distribution P1 (the non-outliers) and 5 APFk s obtained in the same way but
from P1′ (the outliers); G2 is specified in the same way as G1 but replacing P1 and P1′ with
38

0.8

6

0.6

5
4

0.4

3
0.00

0.10

0.20

0.2

A0.2
2

A0.2
1
A0.2
2

0.0

0

2
1

A0.2
1

0.30

0.0

m

0.2

0.4

0.6

0.8

1.0

m

Figure 15: Top panels: The 20%-trimmed mean functions with respect to G1 and G2
when considering APF0 s (left) and APF1 s (right) obtained from the Delaunay-complex
and based on 100 IID points following the distribution P1 (solid curve) or P2 (dotted
curve). Bottom panels: Examples of point patterns with associated APF0 s assigned to
the wrong group, together with the circles of radius 0.8 and 1.
P2 and P2′ , respectively; and we have correctly specified that α = 0.2. Then we simulate
100 APFk s associated to P1 and 100 APFk s associated to P2 , i.e. they are all non-outliers.
Finally, we use (12) to assign each of these 200 APFk s to either G1 or G2 .
0.2

The top panels in Figure 15 show the 20%-trimmed means A1

0.2

and A2

when k = 0

(left) and k = 1 (right). The difference between the 20%-trimmed means is clearest when
k = 1 and so we expect that the assignment error is lower in that case. In fact wrong
assignments happen mainly when the support of P1 or P2 is not well covered by the
point pattern as illustrated in the bottom panels.

39

Repeating this simulation study 500 times, the percentage of APF0 s wrongly assigned
among the 500 repetitions has a mean of 6.7% and a standard deviation of 1.7%, whereas
for the APF1 s the mean is 0.24% and the standard deviation is 0.43%. To investigate how
the results depend on the radius of the smallest circle, we repeat everything but with
radius 0.9 in place of 0.8 when defining the distributions P2 and P2′ . Then for the APF0 s,
the proportion of wrong assignments has a mean of 23.2% and a standard deviation of
2.9%, and for the APF1 s, a mean of 5.7% and a standard deviation of 1.9%. Similar to
Example 9, the error is lowest when k = 1 and this is due to the largest lifetime of a loop.

F Proof of Theorem 4.1
The proof of Theorem 4.1 follows along similar lines as in Chazal et al. (2013) as soon as
we have verified Lemma F.2 below. Note that the proof of Lemma F.2 is not covered by
the approach in Chazal et al. (2013).
We first need to recall the following definition, where BT denotes the topological space

of bounded real valued Borel functions defined on [0, T ] and its topology is induced by
the uniform norm.

Definition F.1. A sequence { Xr }r=1,2,... of random elements in BT converges in distribution

to a random element X in BT if for any bounded continuous function f : BT 7→ R, E f (Xr )

converges to E f (X ) as r → ∞.

Lemma F.2. Let the situation be as in Section 4.2. As r → ∞,

√


r Ar − µ converges in dis-

tribution towards a zero-mean Gaussian process on [0, T ] with covariance function c(m, m′ ) =

Cov ( A1 (m), A1 (m′ )), m, m′ ∈ [0, T ].
Proof. We need some notation and to recall some concepts of empirical process theory.
For D ∈ DTk,nmax , denote A D the APFk of D. Let F = { f m : 0 ≤ m ≤ T } be the class of

functions f m : DTk,nmax 7→ [0, ∞) given by f m ( D ) = A D (m). To see the connection with
40

empirical process theory, we consider
Gr ( f m ) =

√

r

1 r
f m ( Di ) − µ ( m )
r i∑
=1

!

as an empirical process. Denote k · k the L2 norm on F with respect to the distribution

of D1 , i.e. k f m (·)k2 = E A D1 (m)2 . For u, v ∈ F , the bracket [u, v] is the set of all
functions f ∈ F with u ≤ f ≤ v. For any ǫ > 0, N[] (ǫ, F , k · k) is the smallest integer
J

J ≥ 1 such that F ⊂ ∪ j=1 [u j , v j ] for some functions u1 , . . . , u J and v1 , . . . , v J in F with
r


R1
kv j − u j k ≤ ǫ for j = 1, . . . , J. We show below that 0 log N[] (ǫ, F , k · k) dǫ is

finite. Then, by Theorem 19.5 in van der Vaart (2000), F is a so-called Donsker class

which implies the convergence in distribution of Gr ( f m ) to a Gaussian process as in the
statement of Lemma F.2.

For any sequence −∞ = t1 < . . . < t J = ∞ with J ≥ 2, for j = 1, . . . , J − 1, and
for D = {(m1 , l1 , c1 ), . . . , (mn , ln , cn )} ∈ DTk,nmax , let u j ( D ) = ∑in=1 ci li 1(mi ≤ t j ) and

v j ( D ) = ∑in=1 ci li 1(mi < t j+1 ) (if n = 0, then D is empty and we set u j ( D ) = v j ( D ) = 0).

Then, for any m ∈ [0, T ], there exists a j = j(m) such that u j ( D ) ≤ f m ( D ) ≤ v j ( D ), i.e.
J −1

f m ( D ) ∈ [u j , v j ]. Consequently, F ⊂ ∪ j=1 [u j , v j ].

We prove now that for any ǫ ∈ (0, 1), the sequence {t j }1≤ j≤ J can be chosen such that for
j = 1, . . . , J − 1, we have kv j − u j k ≤ ǫ. Write D1 = {( M1 , L1 , C1 ), . . . , ( M N , L N , CN )},

where N is random and should not to be confused with N in Sections 2.1 and A (if
N = 0, then D1 is empty). Let n ∈ {1, . . . , nmax } and conditioned on N = n, let I be

uniformly selected from {1, . . . , n}. Then
E

n

v j ( D1 ) − u j ( D1 )

2

1 ( N = n)

o

(

1
Ci Li 1 Mi ∈ t j , t j+1
n i∑
=1

 2
≤ T 2 n4max E 1 ( N = n) 1 M I ∈ t j , t j+1


≤ T 2 n4max P M I ∈ t j , t j+1 | N = n ,

= n2 E 1 ( N = n )

as n ≤ nmax , Ci ≤ nmax , and Li ≤ T. Further,
E

n

n

v j ( D1 ) − u j ( D1 )
41

2

o

1 ( N = 0) = 0.



)2

Hence


E v j ( D1 ) − u j ( D1 )

2

nmax

=

∑

E

n =0

n

≤ T 2 n5max

v j ( D1 ) − u j ( D1 )
max

n =1,...,nmax

2

1 ( N = n)

o



P M I ∈ t j , t j +1 | N = n .

(13)

Moreover, by Lemma F.3 below, there exists a finite sequence {tn,j }1≤ j≤ Jn such that

P( M I ∈ (tn,j , tn,j+1)| N = n) ≤ ǫ2 / T 2 n5max and Jn ≤ 2 + T 2 n5max /ǫ2 . Thus, by choosing

{ t j } 1≤ j ≤ J =

[

n =1,...,nmax

{tn,j }1≤ j≤ Jn ,

we have J ≤ 2nmax + T 2 n6max /ǫ2 and
max

n =1,...,nmax



P M I ∈ t j , t j +1 | N = n ≤

ǫ2
.
T 2 n5max

Hence by (13), kv j − u j k ≤ ǫ, and so by definition, N[] (ǫ, F , k · k) ≤ 2nmax + T 2 n6max /ǫ2 .
Therefore
Z 1r
0





log N[] (ǫ, F , k · k) dǫ ≤

Z 1q
0

log (2nmax + T 2 n6max /ǫ2 ) dǫ < ∞.

This completes the proof.

Proof of Theorem 4.1. By the Donsker property established in the proof of Lemma F.2

√
√
and Theorem 2.4 in Gine and Zinn (1990), r ( Ar − Ar∗ ) and r Ar − µ converge in
√
distribution to the same process as r → ∞, so the quantile of supm∈[0,T ] r Ar (m) − Ar∗ (m)
√
converges to the quantile of supm∈[0,T ] r Ar (m) − µ(m) . Therefore, q̂αB provides the

bounds for the asymptotic 100(1 − α)%-confidence region stated in Theorem 4.1.

Lemma F.3. Let X be a positive random variable. For any ǫ ∈ (0, 1), there exists a finite

sequence −∞ = t1 < . . . < t J = ∞ such that J ≤ 2 + 1/ǫ and for j = 1, . . . , J − 1,

P X ∈ (t j , t j+1 ) ≤ ǫ.

Proof. Denote by F the cumulative distribution function of X, by F(t−) the left-sided
limit of F at t ∈ R, and by F−1 the generalised inverse of F, i.e. F−1 (y) = inf{ x ∈
42

R : F( x ) ≥ y} for y ∈ R. We verify the lemma with J = 2 + ⌊1/ǫ⌋, t J = ∞, and
t j = F−1 (( j − 1)ǫ) for j = 1, . . . , J − 1. Then, for j = 1, . . . , J − 2,
P X ∈ t j , t j +1
Finally,
P X ∈ t J −1 , t J









= F F−1 ( jǫ)− − F F−1 (( j − 1) ǫ) ≤ jǫ − ( j − 1) ǫ = ǫ.



= P(X > F−1 (( J − 2)ǫ)) = 1 − F F−1 (( J − 2)ǫ) ≤ 1 − ⌊1/ǫ⌋ǫ < ǫ.

References
Baddeley, A., Rubak, E. & Turner, R. (2015). Spatial Point Patterns: Methodology and Applications with R. Chapman and Hall/CRC Press.
Baddeley, A. J. & Silverman, B. W. (1984). A cautionary example on the use of secondorder methods for analyzing point patterns. Biometrics 40, 1089–1093.
Bendich, P., Marron, J., Miller, E., Pieloch, A. & Skwerer, S. (2016). Persistent homology
analysis of brain artery trees. The Annals of Applied Statistics 10, 198–218.
Biscio, C. A. N. & Lavancier, F. (2016). Quantifying repulsiveness of determinantal point
processes. Bernoulli 22, 2001–2028.
Bottou, L. & Bengio, Y. (1995). Convergence properties of the k-means algorithms. In:
Advances in Neural Information Processing Systems, volume 7, MIT Press, 585–592.
Bubenik, P. (2015). Statistical topological data analysis using persistence landscapes.
Journal of Machine Learning Research 16, 77–102.
Chazal, F., Cohen-Steiner, D., Guibas, L. J., Mémoli, F. & Oudot, S. Y. (2009). GromovHausdorff stable signatures for shapes using persistence. Computer Graphics Forum 28,
1393–1403.

43

Chazal, F., Fasy, B., Lecci, F., Rinaldo, A., Singh, A. & Wasserman, L. (2013). On the bootstrap for persistence diagrams and landscapes. Modeling and Analysis of Information
Systems 20, 111–120.
Chazal, F., Fasy, B. T., Lecci, F., Michel, B., Rinaldo, A. & Wasserman, L. (2014). Robust topological inference: Distance to a measure and kernel distance. Available on
arXiv:1412.7197.
Chen, Y.-C., Wang, D., Rinaldo, A. & Wasserman, L. (2015). Statistical analysis of persistence intensity functions. Available on arXiv: 1510.02502.
Daley, D. J. & Vere-Jones, D. (2003). An Introduction to the Theory of Point Processes. Volume
I: Elementary Theory and Methods. Springer-Verlag, New York, 2nd edition.
Edelsbrunner, H. & Harer, J. L. (2010). Computational Topology. American Mathematical
Society, Providence, RI.
Fasy, B., Lecci, F., Rinaldo, A., Wasserman, L., Balakrishnan, S. & Singh, A. (2014). Confidence sets for persistence diagrams. The Annals of Statistics 42, 2301–2339.
Gine, E. & Zinn, J. (1990). Bootstrapping general empirical measures. The Annals of Probability 18, 851–869.
Hartigan, J. A. & Wong, M. A. (1979). Algorithm AS 136: A k-means clustering algorithm. Journal of the Royal Statistical Society: Series C (Applied Statistics) 28, 100–108.
Jacques, J. & Preda, C. (2014). Functional data clustering: a survey. Advances in Data
Analysis and Classification 8, 231–255.
Lavancier, F., Møller, J. & Rubak, E. (2015). Determinantal point process models and statistical inference. Journal of the Royal Statistical Society: Series B (Statistical Methodology)
77, 853–877.
López-Pintado, S. & Romo, J. (2009). On the concept of depth for functional data. Journal
of the American Statistical Association 104, 718–734.

44

López-Pintado, S., Romo, J. & Torrente, A. (2010). Robust depth-based tools for the analysis of gene expression data. Biostatistics 11, 254–264.
MacPherson, R. & Schweinhart, B. (2012). Measuring shape with topology. Journal of
Mathematical Physics 53, 073516.
Matérn, B. (1986). Spatial Variation. Lecture Notes in Statistics 36, Springer-Verlag, Berlin.
Møller, J. & Waagepetersen, R. P. (2004). Statistical Inference and Simulation for Spatial Point
Processes. Chapman and Hall/CRC, Boca Raton.
Møller, J. & Waagepetersen, R. P. (2016). Some recent developments in statistics for spatial point patterns. Annual Review of Statistics and Its Application To appear.
Mrkvička, T., Myllymäki, M. & Hahn, U. (2016). Multiple Monte Carlo testing, with
applications in spatial point processes. Statistics and Computing , 1–17.
Myllymäki, M., Mrkvička, T., Grabarnik, P., Seijo, H. & Hahn, U. (2016). Global envelope
tests for spatial processes. Journal of the Royal Statistical Society: Series B (Statistical
Methodology) Available on arXiv:1307.0239.
Præstgaard, J. T. (1995). Permutation and bootstrap Kolmogorov-Smirnov tests for the
equality of two distributions. Scandinavian Journal of Statistics 22, 305–322.
Robins, V. & Turner, K. (2016). Principal component analysis of persistent homology
rank functions with case studies of spatial point patterns, sphere packing and colloids.
Physica D: Nonlinear Phenomena 334, 99–117.
Sun, Y. & Genton, M. G. (2011). Functional boxplots. Journal of Computational and Graphical Statistics 20, 316–334.
van der Vaart, A. W. (2000). Asymptotic Statistics. Cambridge University Press, Cambridge.
van der Vaart, A. W. & Wellner, J. A. (1996). Weak Convergence and Empirical Processes.
Springer Series in Statistics, Springer-Verlag, New York.

45

