New distance measures for classifying X-ray astronomy
data into stellar classes∗
Amparo Baı́lloa , Javier Cárcamoa and Konstantin Getmanb

arXiv:1603.06806v2 [] 19 Jan 2018

a

b

Departamento de Matemáticas, Universidad Autónoma de Madrid,

28049 Madrid (Spain)
Department of Astronomy and Astrophysics, Pennsylvania State University,
University Park PA 16802-6305 (U.S.A.)

Abstract
The classification of the X-ray sources into classes (such as extragalactic sources,
background stars, . . . ) is an essential task in astronomy. Typically, one of the classes
corresponds to extragalactic radiation, whose photon emission behaviour is well
characterized by a homogeneous Poisson process. We propose to use normalized
versions of the Wasserstein and Zolotarev distances to quantify the deviation of
the distribution of photon interarrival times from the exponential class. Our main
motivation is the analysis of a massive dataset from X-ray astronomy obtained by
the Chandra Orion Ultradeep Project (COUP). This project yielded a large catalog
of 1616 X-ray cosmic sources in the Orion Nebula region, with their series of photon
arrival times and associated energies. We consider the plug-in estimators of these
metrics, determine their asymptotic distributions, and illustrate their finite-sample
performance with a Monte Carlo study. We estimate these metrics for each COUP
source from three different classes. We conclude that our proposal provides a striking
amount of information on the nature of the photon emitting sources. Further,
these variables have the ability to identify X-ray sources wrongly catalogued before.
As an appealing conclusion, we show that some sources, previously classified as
extragalactic emissions, have a much higher probability of being young stars in Orion
Nebula. Keywords: Classification X-ray astronomy Wasserstein distance Zolotarev
metric Photon interarrival time Exponential distribution

Keywords: Classification; X-ray astronomy; Wasserstein distance; Zolotarev metric; photon
interarrival time; exponential distribution.

MSC : Primary 60K35; secondary 62G20, 62N05.
∗

Research by A.B. and J.C. was supported by the Spanish MEyC grants MTM2013-44045-P and
MTM2016-78751-P. K.G. acknowledges the support from the Chandra ACIS Team contract SV4-74018 (G.
Garmire & L. Townsley, PIs), issued by the Chandra X-ray Center, which is operated by the Smithsonian
Astrophysical Observatory on behalf of NASA under contract NAS8-03060.

1

1

Introduction and motivation

An important initial step in the analysis of stellar populations is the classification of
samples into different classes of sources (see [5]). The definition of the classes (foreground
stars, background stars, different types of pre-main-sequence stars, etc.) depends on the
research project, but it is always of interest to identify extragalactic sources (see [5]; [9]).
Frequently, the allocation has a degree of uncertainty, to the extent that some of the
astronomical sources might remain unclassified (see [11]) or even wrongly catalogued.
X-ray astronomy deals with the detection and observation of astrophysical objects by
means of the properties of their X-ray emissions. There are many astronomical sources
of X-rays, such as galaxy clusters, black holes or different types of stars. In X-ray
astronomy, classification of the data (that is, the X-ray sources) is accomplished using all
the information provided by source features such as its location and X-ray and infrared
properties (see [5]). As X-radiation is blocked by the atmosphere of Earth, cosmic X-ray
emissions can only be detected by space telescopes.
This article is motivated by a real dataset obtained as a result of Chandra Orion Ultradeep
Project (COUP). It was fulfilled with one of the “Great Observatories” of NASA, the
Chandra X-ray space telescope. Chandra was designed to observe X-ray emissions from
high-energy regions of the space such as supernovas, black holes or star clusters as the
Orion Nebula.
In this work, we focus on a massive collection of X-ray astronomical sources derived
from a 2003 exposure of Chandra to the Orion Nebula region ([12]). For each of the
sources captured by Chandra, the photon arrival times and associated energies were
collected during a nearly continuous observation period of almost 10 days. The majority
of these X-ray sources have been classified into one of three groups: lightly-obscured and
heavily-obscured low-mass young stars; and extragalactic sources. The X-ray classification
of young stellar objects in star forming regions is in general a complicated task, where
numerous source properties are used as features.
A very informative fact employed in this stellar classification is that, on the one hand,
extragalactic radiation usually has a constant photon emission rate. This can be illustrated
via the light curve of the astronomical source, a graph depicting its brightness (measured,
e.g., by the photon count rate) over the course of the observation period. The light curves of
COUP sources 111 (Figure 1 (a)) and 1304 (Figure 1 (b)), classified as extragalactic in [11]
are examples of a constant photon arrival rate. In this case, the point process constituted by
the photon arrival times is well-modeled by a homogeneous Poisson process. Thus, photon
interarrival times from an extragalactic source should be close to an exponential distribution.
On the other hand, young stars usually exhibit high-amplitude rapid variability ([29])
and their photon arrivals are generally affected by flares so the corresponding interarrival
distribution might deviate from the exponential one. As an example, Figure 1 (c) displays
the light curve of COUP source 89, corresponding to a young star, where we can see a
large flare at about 80 hours from the start of the observation period. In the astrophysical
literature, there are different proposals to quantify photon emission variability in stellar
X-ray sources (see [29]).

2

1.5
1.0
0.0

0.5

Count rate (counts/ks)

2.0

2.5

COUP source 111

0

50

100

150

200

250

300

200

250

300

200

250

300

Time (hrs)

(a)

50
40
30
20
0

10

Count rate (counts/ks)

60

COUP source 1304

0

50

100

150
Time (hrs)

(b)

15
10
0

5

Count rate (counts/ks)

20

25

COUP source 89

0

50

100

150
Time (hrs)

(c)
3
Figure 1: Light curves for COUP sources (a) 111 (extragalactic radiation); (b) 1304
(extragalactic radiation); (c) 89 (lightly obscured PMS star).

empirical
distribution
function Fn
of PIT

spatial
source
of X-rays

raw data:
photon arrival times
and
associated energies

photon
interarrival
times (PIT)
x1 , . . . , x n

distance of PIT
to exponential
d(Fn , Gµ̂ )

sample
mean µ̂
of PIT

Figure 2: Data analysis pipeline for each single cosmic X-ray source (input): the raw
data is the series of photon arrival times and their corresponding energies. PIT are
computed as differences between consecutive arrival times. The sample mean and the
empirical distribution function of PIT produce the output (the empirical distance of the
PIT distribution to the exponential one).
Here, we propose a new statistical methodology to quantify the deviation of a random
variable (namely, the photon interarrival time or PIT) from the exponential class. The
final aim is to generate a new input variable for the discriminant procedure distinguishing
among the source classes, and particulary extragalactic ones (see the data analysis pipeline
in Figure 2). We consider that a large estimated distance of the PIT to the exponential
class is an evidence that the corresponding source is not extragalactic. Specifically, we
use a normalized version of the so-called Wasserstein and Zolotarev ζ2 metrics, between
the photon interarrival times of each X-ray source and the exponential distribution. As
mentioned in [21, Section 15], the Zolotarev ζ2 -metric is appropriate when dealing with
exponential variables. Further, [22] argue that Wasserstein and Zolotarev distances
are more sensitive to extreme values than other probability metrics such as the usual
Kolmogorov distance. In general, it is often desirable to take into account extreme events
to compare distributions. This is specially relevant with data coming from astrophysical
studies (see [8]). We demonstrate that these distances can be used as informative variables
that detect groups or similitudes among the X-ray sources and help identify possible
outliers within a group. In this work the term “outlier” refers to a source whose distance to
the exponential distribution is substantially different from the distances of other members
of the same class. In fact, in the final analysis of the COUP data, we show that some of
the outlying COUP sources, initially classified as extragalactic radiation, could actually
be young stars in the Orion Nebula region.
The ideas in this paper are also potentially useful in other biological or physical problems
in which deviations from the exponential model need to be detected and quantified.
Alternatively, the proposed methodology allows assessing whether a homogeneous Poisson
process achieves a good approximation of an observed phenomenon.

4

This paper is structured as follows. In the next section, we describe in detail the COUP
dataset. In Section 3, we introduce the Wasserstein and Zolotarev distances and normalized
versions of them. We also consider the plug-in estimators of these metrics to be used
in practice and determine their asymptotic distributions. In Section 4, we carry out a
simulation study to assess the practical performance of our proposal with finite samples.
The COUP dataset is analyzed in depth in Section 5. Finally, the main conclusions are
collected in Section 6. In the Appendix, we include the proofs of the results stated in
Section 3 and other technical details.

2

Chandra Orion Ultradeep Project Dataset

Among other things, the COUP analyzes X-flaring in pre-main-sequence (PMS) stars,
members of the Orion Nebula region that is composed of the rich revealed Orion Nebula
Cluster (ONC) and the filamentary molecular cloud called Orion Molecular Cloud 1 (OMC1). A PMS star is a premature star that has acquired all of its mass from its natal envelope
of interstellar dust and gas and contracts until it starts hydrogen burning (see, e.g., [23]).
These young stars have intense magnetic fields, detected through their X-ray emissions,
where plasma, confined in magnetic loops, is heated to X-ray emitting temperatures.
Detection of these high-energy emissions is only possible by space observatories, such as
the Chandra X-ray Observatory (see http://chandra.si.edu/). The nearest rich and
concentrated collection of PMS stars is in the ONC/OMC-1 star forming region.
In January 2003, the Chandra X-ray Observatory focused its Advanced CCD Imaging
Spectrometer on the ONC for a period of 13.2 days obtaining the deepest X-ray observation
ever taken of any star cluster (see [12]). This observation period was only interrupted
by the five passages of the Chandra spacecraft through the Van Allen radiation belts.
Graphically, this is reflected in the five gaps appearing in the light curves of Figure 1. The
results were an almost continuous observation of the photon arrival times and associated
energies for 1616 X-ray sources. COUP sources were compared with source positions from
previously existing catalogs, with the aim of physically associating the COUP sources
with already identified stars (whenever this was possible). The majority of these COUP
sources has been classified into one of three groups (see [9]):
– Lightly-obscured PMS sources: This class is constituted by 835 cool low-mass PMS
stars that are likely located in the ONC cluster. The term “lightly obscured” means
that the star X-ray emission is less absorbed by the material in the interstellar medium.
– Heavily-obscured PMS sources: This class corresponds to 559 low-mass PMS stellar
objects that are likely still embedded in the nascent OMC-1 cloud.
– Nonmembers: This group contains over 200 probable nonmembers of the Orion Nebula
star forming region. A large part are likely extragalactic sources and a few are foreground
stars or very faint sources without counterparts. For our analysis, in this group we only
consider the extragalactic X-ray emissions.
The original classification of these 1594 COUP sources is detailed in [11], who provide
an estimated probability of membership to the Orion cloud for each source. Thus, this
5

classification has a degree of uncertainty. The aim of this work is to introduce new
distance measures that contain relevant information for classifying X-ray astronomy data
into stellar classes. Specifically, we compute some probability metrics of the PIT to the
exponential class, as they exhibit different distributions depending on the nature of the
photon emission sources. These distances can be incorporated in any classification rule to
reduce the classification error.

3

Quantifying the discrepancy from extragalactic radiation

As mentioned in the introduction, PIT resulting from extragalactic radiation are usually
well-modeled by exponential distributions while, for the PMS stars classes, PIT distributions
might deviate from the exponential one. Our primary objective is to compare different
sources by computing some normalized probability metrics between their corresponding
PIT and the nearest exponential variable, and to highlight the classification power of these
features. This idea has been tackled before: for instance, in X-ray astronomy the usual
Kolmogorov distance to the exponential class is often used in data classification. However,
depending on the problem at hand, other distances could be more appropriate, e.g., to
take into account extreme values or to highlight a specific part of the distribution. Here,
we analyze in depth the performance of the Wasserstein and Zolotarev distances, more
sensitive to the behaviour in the tail of the distribution than others such as the usual
Kolmogorov and Cramér-von Mises metrics. In particular, a byproduct of the results in
this section is the possibility of testing for exponentiality (see the end of Subsection 3.2),
for which there are numerous proposals in the literature (see the reviews by [2] and [16]).
3.1

The choice of the metrics

A central question in the problem under consideration is the choice of a suitable metric to
measure how far the probability distribution of the positive random variable X of interest
(the PIT), with expectation µ > 0, is from the exponential variable Yµ with the same
mean. In this work, we focus on the family of integral probability metrics
dr (X, Yµ ) = sup {|Ef (X) − Ef (Yµ )| : f ∈ Fr } ,

r ∈ N,

(1)

where Fr is the class of real-valued functions f on R having r-th derivative f (r) a.e. and
such that |f (r) | ≤ 1 a.e. Observe that dr (X, Yµ ) is the maximum error in the expected
value within the class of smooth functions Fr due to the approximation of X by Yµ . For
notational convenience, if F and Gµ are the distribution functions of X and Yµ , respectively,
we indistinctly use dr (F, Gµ ) or dr (X, Yµ ). Note also that Gµ (x) = 1 − exp(−x/µ), for
x ≥ 0.
The case r = 1 in (1) has special relevance. By the Kantorovich–Rubinstein theorem, we
see that d1 ≡ ω is the famous L1 -Wasserstein distance. For r ≥ 2, dr ≡ ζr is the Zolotarev
metric of order r. For a general reference on these distances we refer to [21].
It can be proved that, when dr (X, Yµ ) < ∞, the moments of X and Yµ coincide up to order
r − 1. As a consequence, in general, when r ≥ 3 and X is not exponential, we have that
dr (X, Yµ ) = ∞. This follows from the fact that, for many distributions, equalities EX = µ
and EX 2 = 2 µ2 are too restrictive and they actually imply that X is exponential. For
6

instance, this happens for the variables with the HNBUE or HNWUE property, two large
families of random variables that include all the usual ageing and anti-ageing classes of
distributions as it follows from results on stochastic equality under convex domination (see,
e.g., [25, Theorem 3.A.42, p. 133]). Therefore, only the cases r = 1 and r = 2 make sense
for the discussed problem, that is, the Wasserstein distance) d1 ≡ ω and the Zolotarev
metric d2 ≡ ζ2 . These two metrics have easier-to-handle dual integral representations (see
[22]), given by
Z ∞
ω(F, Gµ ) =
|F (t) − Gµ (t)| dt,
(2)
0

Z

∞

Z

∞

(F (x) − Gµ (x)) dx dt.

ζ2 (F, Gµ ) =
0

(3)

t

The metrics ω and ζ2 have practical advantages for the problem at hand. First, as argued
in [22, p. 15], they are more sensitive to the differences in the probabilities corresponding
to extreme values than other common probability metrics such as the Kolmogorov distance,
κ(F, G) = supx |F (x) − G(x)|. Since the difference |F (x) − G(x)| converges to zero as x
tends to +∞ or −∞, the contribution of the terms corresponding to extreme events is
usually small. As a consequence, the differences in the tail behavior of X and Y will only
be reflected in κ(F, G) to a relatively small extent. However, representations (2) and (3)
show that extreme values have more weight in ω and ζ2 , as integrals of tail probabilities
appear in these distances. Additionally, Zolotarev ζ2 -metric is considered as the “natural
metric” when dealing with the exponential class (see [21, p. 340]).
It becomes soon apparent that these two metrics have the important drawback of not being
scale-independent, an essential requirement to compare X-ray sources with very different
photon emission rates (the average of photon emissions in each time interval). Even for
sources of the same nature, the distances ω and ζ2 to the exponential distribution could be
extremely different. This is clearly reflected in Figure 1, where COUP 111 and 1304, both
classified as extragalactic, would not be comparable without a suitable normalization.
To overcome this problem, next we introduce dimensionless versions of the ω and ζ2
distances. Let us observe that the Wasserstein ω and Zolotarev ζ2 metrics are homogeneous
of order 1 and 2, respectively. We define the normalized Wasserstein and normalized
Zolotarev metrics as
1
ω̄(X, Yµ ) = ω(X/µ, Yµ /µ) = ω(X, Yµ )
µ
and
ζ̄2 (X, Yµ ) = ζ2 (X/µ, Yµ /µ) =

1
ζ2 (X, Yµ ).
µ2

These are the two probability distances that should be used in practice (instead of their
unnormalized versions), as they are homogeneous of degree 0 and dimensionless.
3.2

Estimation and large sample behaviour

In practice, the distribution of the random variable X is usually unknown. Therefore, the
distances d(X, Yµ ) (for d = ω, ζ2 , ω̄ and ζ̄2 ) have to be estimated using a random sample
7

X1 , . . . , Xn from X. We propose to use the plug-in estimators obtained by replacing the
true distributionPF of X by the empirical distribution Fn of the sample X1 , . . . , Xn , that
is, Fn (t) = n−1 ni=1 I{Xi ≤t} , n ∈ N, t ≥ 0, where IA stands for thePindicator function of
the set A. Thus, we estimate d(F, Gµ ) by d(Fn , Gµ̂ ), where µ̂ = n1 ni=1 Xi is the sample
mean (and the maximum likelihood estimator of the rate parameter µ of an exponential
model).
Analyzing the asymptotic behavior of the empirical distances d(Fn , Gµ̂ ) is an important
issue to understand its performance and accuracy in practice, for instance, to assess
whether an exponential model provides a reasonably good approximation of X. Besides,
the asymptotic probability distribution potentially allows performing inference on d(X, Yµ ).
For the considered distances, we note that
1
d(Fn , Gµ̂ ) = d(F, Gµ ) + √ δn (d, F ),
n
where δn (d, F ) is the standardized version of the estimated distances
√
δn (d, F ) := n (d(Fn , Gµ̂ ) − d(F, Gµ )) , n ∈ N.

(4)

Here we find conditions (as sharp as possible) on the random variable X so that δn (d, F )
converges in distribution as n → ∞, and determine its weak limit, δ∞ (d, F ). In this way,
we obtain that
√
d(Fn , Gµ̂ ) = d(F, Gµ ) + OP (1/ n), as n → ∞.
Though the detailed proofs of the asymptotic distribution of δn (d, F ) are collected in the
Appendix, we describe here in broad strokes the main ideas behind them. First, we note
that
δn (d, F ) = ρn (Xd,n , gd ),
(5)
where ρn : L1 × L1 → R is the functional defined by
√
√
ρn (f, g) := kf + ngk1 − nkgk1 ,

for f, g ∈ L1 ,

(6)

Xd,n are the stochastic processes given (for t ≥ 0) by
√
Xω,n (t) := n [(Fn (t) − Gµ̂ (t)) − (F (t) − Gµ (t))] ,


√ 1
1
Xω̄,n (t) := n (Fn (t) − Gµ̂ (t)) − (F (t) − Gµ (t)) ,
µ̂
µ
Z ∞

Z ∞
√
(Fn (x) − Gµ̂ (x)) dx −
(F (x) − Gµ (x)) dx ,
Xζ2 ,n (t) := n
t

Xζ̄2 ,n (t) :=

√



1
n 2
µ̂

(7)

t

Z

∞

t

1
(Fn (x) − Gµ̂ (x)) dx − 2
µ

Z

∞


(F (x) − Gµ (x)) dx ,

t

and gd are the (deterministic) functions defined by
1
(F (t) − Gµ (t)),
µ
Z ∞
1
gζ̄2 (t) := 2
(F (x) − Gµ (x)) dx.
µ t

gω (t) := F (t) − Gµ (t),
Z

gω̄ (t) :=

∞

(F (x) − Gµ (x)) dx,

gζ2 (t) :=
t

8

(8)
(9)

From (5), we see that establishing the (weak) convergence in L1 of the processes Xd,n in
(7), combined with the continuity of the linking functional in (6), immediately translates
into the convergence in distribution of δn (d, F ).
Before stating the main results, we need to introduce some definitions and notation. In
the sequel, BF := B ◦ F is the F -Brownian bridge, where B is a standard Brownian bridge
on [0, 1], that is, B is a centered Gaussian process with covariance function γ(s, t) =
min(s, t) − st and continuous paths, with probability 1.
We consider the Lorentz spaces of positive random variables defined by L2,1 := {X : Λ2,1 (X) < ∞}
and L4,2 := {X : Λ4,2 (X) < ∞}, where
Z ∞p
Z ∞ p
Λ2,1 (X) :=
t P(X > t) dt
P(X > t) dt and Λ4,2 (X) :=
0

0

(see [20, p. 279]). Conditions Λ2,1 (X) < ∞ and Λ4,2 (X) < ∞ are slightly stronger than
EX 2 < ∞ and EX 4 < ∞, respectively (see [15]). Finally, Lp := {X : EX p < ∞} (p > 0)
is the usual space of (positive) random variables with finite p-th moment.
The following two theorems characterize the asymptotic behavior of Xd,n in L1 . The results
are sharp in the sense that we obtain the exact integrability condition on X so that the
L1

processes converge in distribution in L1 as n → ∞. The symbol “−→w ” stands for the
weak convergence of a sequence of random processes in the space L1 as n → ∞ (see the
Appendix for the precise definition).
Theorem 1. Let X be a positive random variable with expectation µ > 0. If X ∈ L4/3 ,
the following assertions are equivalent:
(a) X ∈ L2,1 .
L1

(b) Xω,n −→w Xω,F , where Xω,F is a centered Gaussian process given by
Z
t −t/µ ∞
Xω,F (t) := BF (t) − 2 e
BF (s) ds, t ≥ 0.
µ
0

(10)

L1

(c) Xω̄,n −→w Xω̄,F , where Xω̄,F is a centered Gaussian process given by



Z ∞
t −t/µ
1
BF (t) + gω̄ (t) − 2 e
BF (s) ds , t ≥ 0,
Xω̄,F (t) :=
µ
µ
0
and the function gω̄ is defined in (8).
Theorem 2. Let X be a positive random variable with expectation µ > 0. The following
assertions are equivalent:
(a) X ∈ L4,2 .
L1

(b) Xζ2 ,n −→w Xζ2 ,F , where Xζ2 ,F is a centered Gaussian process given by


Z ∞
Z ∞
t
−t/µ
Xζ2 ,F (t) :=
BF (s) ds − 1 +
e
BF (s) ds, t ≥ 0.
µ
t
0
9

L1

(c) Xζ̄2 ,n −→w Xζ̄2 ,F , where Xζ̄2 ,F is a centered Gaussian process given by


Z ∞
Z ∞


1
t
−t/µ
Xζ̄2 ,F (t) := 2
BF (s) ds ,
e
BF (s) ds + 2 µ gζ̄2 (t) − 1 +
µ
µ
0
t
(11)
for t ≥ 0, and the function gζ̄2 is defined in (9).
Using (5) and Theorems 1 and 2, in the next theorem we derive the asymptotic distribution
of δn (d, F ). In the sequel “−→d ” stands for convergence in distribution as n → ∞, sgn(·)
denotes the sign function and Ac is the complement of the set A.
Theorem 3. Let X be a positive random variable with expectation µ > 0. For d = ω
or d = ω̄ (respectively, for d = ζ2 or d = ζ̄2 ), let us assume that X ∈ L2,1 (respectively,
X ∈ L4,2 ). Then, δn (d, F ) −→d δ∞ (d, F ), with
Z
Z
δ∞ (d, F ) :=
|Xd,F (t)| dt +
Xd,F (t) sgn(gd (t)) dt,
(12)
I(gd )c

I(gd )

where the processes Xd,F are defined in (10)-(11), the functions gd are given in (8)-(9),
and I(gd ) := {t ≥ 0 : gd (t) = 0}.
The next corollary, a direct consequence of Theorem 3, provides the asymptotic distribution
of δn (d, F ), when F is an exponential distribution function. In such √
a case, d(F, Gµ ) = 0
and the estimators behave as a random quantity at the order of OP (1/ n). It is interesting
to note that the limiting distribution of the normalized distances does not depend on the
unknown mean of the exponential distribution.
Corollary 1. For the processes Xd,F defined in (10)-(11), if X follows an exponential
distribution with mean µ, then
(a)
(b)
(c)
(d)

√
√
√
√

n ω(Fn , Gµ̂ ) −→d kXω,Gµ k1 = µ kXω,G1 k1 ;
n ω̄(Fn , Gµ̂ ) −→d kXω̄,G1 k1 = kXω,G1 k1 ;
n ζ2 (Fn , Gµ̂ ) −→d kXζ2 ,Gµ k1 = µ2 kXζ2 ,G1 k1 ;
n ζ̄2 (Fn , Gµ̂ ) −→d kXζ̄2 ,G1 k1 = kXζ2 ,G1 k1 .

The following corollary states that when X does not share any part of its distribution
function with the exponential one, the limiting distribution δ∞ (d, F ) in (12) is actually
normal, for all the considered distances.
Corollary 2. Let us assume that the conditions of Theorem 3 hold and let us further
assume that the set I(gd ) (defined in Theorem 3) has zero Lebesque measure. Then,
δ∞ (d, F ) defined in (12) has a zero mean normal distribution.
The previous results could be useful to compute (asymptotic) confidence intervals for
d(F, Gµ ). This can be implemented via the following bootstrap procedure: as F is usually
10

unknown, we substitute F for Fn in the limiting processes obtained in Theorems 1 and 2
(equations (10)-(11)). Next, we simulate a large number of trajectories of the processes
to obtain a Monte Carlo approximation of the asymptotic distribution δ∞ (d, F ) given in
Theorem 3. Finally, we use the Monte Carlo sample quantiles to construct the desired
interval. We also note that if we assume that the Lebesgue measure of the sets I(gd ) is zero,
then the procedure is simpler as Corollary 2 ensures that the limit distribution is a zero
mean normal distribution. Hence, in such a case it is enough to estimate the asymptotic
variance of the limit via Monte Carlo and use the quantiles of a normal distribution. This
latter interval is called standard normal interval in [7, p. 168]. However, let us observe
that asymptotic confidence intervals could be unprecise when the sample size is small.
As the distribution of δn (d, F ) for a fixed n could be extremely difficult to handle, an
interesting alternative is to construct a bootstrap confidence interval of d(F, Gµ ). In this
situation, the percentile interval proposed in [7, p. 170] is a reasonable choice.
As d(F, Gµ ) = 0 is equivalent to saying that F is exponential, the considered distances
can be additionally applied to goodness-of-fit tests for H0 : F is exponential.
√ As stated in
Corollary 1, for d = ω̄, ζ̄2 , the asymptotic distribution of the test statistic n d(Fn , Gµ̂ ) is
completely determined under H0 . In practice, this result allows us to derive an asymptotic
rejection region by Monte Carlo sampling from the asymptotic distribution. When the
sample size n is small, a better alternative is to use a parametric bootstrap procedure by
sampling from an exponential distribution with mean µ̂.

4

Simulations

The aim of this section is to analyze the finite-sample behaviour of the statistic δn (d, F )
given in (4), in particular to compare it with its asymptotic distribution δ∞ (d, F ) obtained
in Theorem 3. We only consider the normalized versions of the distances, d = ω̄ and
d = ζ̄2 . In any case, the simulations results for the unnormalized metrics are similar and
do not add relevant information.
To compute the normalized empirical distances ω̄(Fn , Gµ̂ ) = ω(Fn , Gµ̂ )/µ̂ and ζ̄2 (Fn , Gµ̂ ) =
ζ2 (Fn , Gµ̂ )/µ̂2 , we have used the following equalities
Z X(n)
ω(Fn , Gµ̂ ) = X(1) − µ̂ Gµ̂ (X(1) ) +
|Fn (x) − Gµ̂ (x)| dx + µ̂ e−X(n) /µ̂
(13)
X(1)

and
Z

∞

Z



0

0

X(n)

+

Z

X(1)

a2
2
!

t

−X(1) + µ̂ Gµ̂ (X(1) ) +

=2

dt + µ̂2 −

(Fn (x) − Gµ̂ (x)) dx

ζ2 (Fn , Gµ̂ ) = 2
Z

t

dt + µ̂2 −

(Fn (x) − Gµ̂ (x)) dx
X(1)

+

a2
, (14)
2

P
where X(1) ≤ · · · ≤ X(n) are the order statistics of the sample and a2 := ni=1 Xi2 /n.
Even though the integrals appearing in (13) and (14) can be expressed in terms of the
order statistics, from a computational viewpoint it is more convenient to approximate
them numerically: this was carried out by discretizing the integral on the equispaced grid
X(1) + k δ with k = 0, 1, 2, . . . , 20000 and δ = (X(n) − X(1) )/20000.
11

The asymptotic distribution δ∞ (d, F ) can be approximately sampled by generating trajectories of the Brownian bridge
R ∞ BF (t) on a bounded time interval [0, T ] and then approximating
the integrals such as 0 BF (t) dt by their discretization on [0, T ]. In our work, we have
chosen T as the smallest integer larger or equal to F −1 (1 − tol), where tol is a tolerance
limit equal to 10−6 . The integral discretization was carried out with an equispaced grid on
[0, T ] yielding 50000 subintervals.
Computing the normalized version, d = ω̄ or d = ζ̄2 , of the distances d(F, Gµ ) is equivalent
to computing the distance of the re-scaled variable X/µ to the exponential distribution
with mean µ = 1. As a consequence, in this Monte Carlo study the data-generating
distributions have expectation µ = 1. Specifically, we have considered:
- the exponential distribution with mean µ = 1;
- the Weibull distribution with shape parameter a > 0 and scale parameter λ = 1/Γ(1 +
a
1/a), with probability density f (x) = a/λ (x/λ)a−1 e−(x/λ) , for x > 0;
- the gamma distribution with shape parameter a > 0 and scale parameter λ = 1/a, with
density f (x) = aa /Γ(a) xa−1 e−ax , for x > 0.
For the Weibull distribution with mean 1 we have chosen the values a = 0.9 and a = 1.1
for the shape parameter. For the gamma distribution with mean 1, we have used shape
parameters a = 0.9 and a = 1.1 (see Figure 3).
In the simulations we have generated 10000 samples of size n = 100, 500, 1000 and 5000
from each of these distributions. For each Monte Carlo sample, we have computed the
statistic δn (d, F ) for d = ω̄ and d = ζ̄2 . The programming language used in this work is R
(www.R-project.org). The results of the simulations are summarized in Figures 4, 5, 6, 7
and 8. Each figure shows the evolution, as n increases, of the boxplots of these statistics
δn (d, F ) towards the boxplot of its asymptotic distribution δ∞ (d, F ) given in (12). This
latter boxplot is also based on 10000 samples from the corresponding limit distribution.
On the one hand, we observe that, the finite-sample behavior of δn (d, F ), both for d = ω̄
and for d = ζ̄, is stable for the exponential distribution. In particular, the quartiles and
median of δn (d, F ) are very similar for any n. For nonexponential distributions F , the
boxplot of δn (d, F ) resembles that of δ∞ (d, F ) for large sample sizes (n ≥ 1000). On the
other hand, we also observe that, the closer F is to the exponential distribution, the larger
n has to be for the distribution of δn (d, F ) to approach its limit. This is reasonable: when
F is not exponential, but close to it and the set I(gd ) has zero Lebesgue measure, the finite
sample behavior of δn (d, F ) is almost as if F were exponential, but the limit is actually
Gaussian (see Corollary 2).

5

Analysis of the COUP data

As mentioned in the introduction, photon interarrival times resulting from extragalactic
radiation are usually well-modeled by exponential distributions, while, for the classes of
PMS stars, PIT distributions might deviate from the exponential one. Hence, it is natural
to compare different sources by computing the distances between their corresponding
PIT and an exponential variable with the same mean. In this section, we compute the
12

0.2

0.4

0.6

0.8

exp(1)
Weibull(0.9)
Weibull(1.1)
Gamma(0.9)
Gamma(1.1)

0.5

1.0

1.5

2.0

2.5

3.0

x

●

6

3.0

Figure 3: Probability densities in the simulation study of Section 4.

●

●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●

●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●

●

●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

5000

∞

●

●
●

4

●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●

●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

0

0.5

1

1.0

2

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●

3

●
●

2.0

●
●

●

1.5

5

2.5

●

●

100

500

1000

100

n

500

1000

5000

∞

n

(a)

(b)

Figure 4: Boxplots of δn (d, F ) and its asymptotic distribution δ∞ (d, F ) for (a) d = ω̄ and
(b) d = ζ̄2 . The distribution F is exponential(1).

13

●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●

−1

−1

0

0

1

1

●

2

2

●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

−2

100

500

1000

●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●

●

●
●
●
●

5000

∞

−2

●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

100

500

1000

n

n

(a)

(b)

●
●
●
●
●
●

5000

∞

3

Figure 5: Boxplots of δn (d, F ) and its asymptotic distribution δ∞ (d, F ) for (a) d = ω̄ and
(b) d = ζ̄2 . The distribution F is Weibull with shape parameter a = 1.1.

●

●

●
●
●

●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●

●

●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

10

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●

●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●

0

−1

0

5

1

2

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●

−2

●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●

100

500

1000

5000

●
●
●
●
●
●
●

−5

●
●

∞

●

100

500

1000

n

n

(a)

(b)

5000

●
●
●
●
●
●
●
●
●
●
●
●

∞

Figure 6: Boxplots of δn (d, F ) and its asymptotic distribution δ∞ (d, F ) for (a) d = ω̄ and
(b) d = ζ̄2 . The distribution F is Weibull with shape parameter a = 0.9.

14

6

●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●

●
●
●

4

2

●
●
●

●

●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●

−1

0

0

2

1

●
●
●
●

●
●
●

−2

●
●

−2

●

●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●

●

100

500

1000

5000

∞

100

500

1000

n

n

(a)

(b)

5000

●
●
●
●
●

∞

8

Figure 7: Boxplots of δn (d, F ) and its asymptotic distribution δ∞ (d, F ) for (a) d = ω̄ and
(b) d = ζ̄2 . The distribution F is gamma with shape parameter a = 1.1.

●

●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

−2

−1

0

0

2

1

4

2

●
●
●
●
●
●
●

6

●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

100

500

1000

5000

●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

5000

∞

−4

−2

●
●

∞

100

500

1000

n

n

(a)

(b)

Figure 8: Boxplots of δn (d, F ) and its asymptotic distribution δ∞ (d, F ) for (a) d = ω̄ and
(b) d = ζ̄2 . The distribution F is gamma with shape parameter a = 0.9.

15

empirical normalized Wasserstein and Zolotarev distances, ω̄(Fn , Gµ̂ ) and ζ̄2 (Fn , Gµ̂ ), to
the exponential class for the sample of PIT of each X-ray source in the COUP dataset,
described in Sections 1 and 2. For the sake of comparison, we also include the usual
Kolmogorov metric κ(Fn , Gµ̂ ) in the analysis as it is commonly used by astrophysicists
in the classification of X-ray sources. The three types of sources of interest (namely,
lightly-obscured PMS stars, heavily-obscured PMS stars and extragalactic) are compared
via the values of these metrics.
For each of the X-ray sources in the three groups, we have computed the series of times
between consecutive photon detections, taking into account the five observation gaps due
to the passages through the high-radiation belts (see Section 2). We have also kept in
mind that the complicated Chandra ACIS X-ray background, which is the combination
of both celestial and instrumental backgrounds, tends to have peaks at energies below
0.5 keV and above 8 keV. To remove a significant portion of the X-ray background, we
have followed [12], who apply an energy filter to the data by cleaning X-ray events with
energies out of the total (0.5-8) keV band.
In order to obtain reasonably good estimates of the distances to the exponential distribution,
we have kept only the COUP series with at least 100 PIT. We have finally analyzed 1090
samples of PIT, of which 73, 644 and 373 had been previously classified as extragalactic
sources, lightly-obscured and heavily-obscured PMS stars, respectively. For each of these
1090 COUP sources, we have computed the distances κ(Fn , Gµ̂ ), ω̄(Fn , Gµ̂ ) and ζ̄2 (Fn , Gµ̂ )
of the empirical distribution of PIT to the exponential distribution with the same sample
mean.
The interesting issue is whether the distance of PIT to the exponential distribution depends
on the source class. The results are summarized in Figures 9 and 10. In Figure 9 we
have displayed the boxplots of log(d(Fn , Gµ̂ )), for d = κ, ω̄, ζ̄2 , separated according to
the three types of COUP sources. Outliers have been identified by the COUP source
number (see [12]). As expected, for the three distances we see that the distribution of PIT
due to extragalactic radiation is the nearest to the exponential
distribution. Moreover,
√
Figure 10 displays the empirical distribution functions of n d(Fn , Gµ̂ ) (d = ω̄, ζ̄2 ), for the
different COUP groups and the distribution function of δ∞ (d, G1 ) when the distribution
function is exponential with mean 1. The empirical distribution function corresponding to
the extragalactic group is again the nearest to the distribution function of δ∞ (d, G1 ), the
asymptotic distribution in Corollary 1. Additionally, Figure 9 shows that the normalized
Zolotarev distance separates better the three classes than the normalized Wasserstein
and the Kolmogorov-Smirnov metrics. Observe also that the PIT corresponding to
lightly-obscured PMS stars are noticeably nearer to the exponential than those of the
heavily-obscured group. A plausible explanation for this difference is the following: X-ray
emission in heavily-obscured PMS stars is more affected by the intervening interstellar
absorption, mainly the absorption from the gas in the molecular cloud and/or local
absorption in an envelope or a disk around a star. Other factor might also be that the
heavily absorbed sample is generally younger and more “diskier”, in the sense that it has a
higher fraction of stars that are still surrounded by their circumstellar disks; some of them
highly accreting. Thus, the age and the presence/absence of disks could play an additional
role, by affecting the X-ray production mechanisms of the strong X-ray flares in PMS stars

16

and in turn affecting their X-ray photon arrival patterns (see [13] and [14] for differences
in the strong X-ray flares between the diskless and disky-accreting stellar populations).
The normalized metrics ω̄ and ζ̄2 are also capable of detecting outlying COUP sources
(see Figure 9), which are interesting to interpret. The X-ray light-curves and the photon
arrival diagrams in the COUP atlas (see [12, Figure Set 12]), a numerical and graphical
summary of all the COUP sources, show one thing in common among different types of
outliers: the outlying sources nearest to the exponential distribution do not exhibit long
and powerful X-ray flares, while all outliers farthest from 0 do. In other words, some
outliers in Figure 9 could be in fact misclassified X-ray sources. For instance, the COUP
sources 751, 895, 961, and 1015, classified as extragalactic sources by [11] and having the
highest normalized Wasserstein and Zolotarev distances to the exponential distribution
among those sources of the extragalactic sample (Figure 9), are likely to be young stellar
objects in the Orion Nebula region.
To analyze the nature of the above-mentioned four outliers with more precision, we
decided to incorporate another informative variable with the potential of discriminating
among the three COUP classes specified above. Taking into account that, for each of
the photons belonging to a COUP source, we observe both the photon arrival time and
its associated energy, it is natural to choose some feature summarizing the energies in
the COUP series. The COUP atlas provides extra information (in the form of tabulated
quantities for each COUP source). Specifically, we have considered the quantity MedEn,
the median energy (in keV) of the source, due to its strong correlation with the absorbing
column density characterizing interstellar absorption (see [9]). To further clarify why
MedEn is an adequate choice, in Figure 11 we have plotted the empirical distribution
functions of the photon energies rescaled to [0,1] for all the COUP sources. We clearly see
that MedEn discriminates well between lightly obscured stars and extragalactic sources
whereas distances to the exponential class separate better extragalactic sources from
heavily obscured ones. Consequently, together these two quantities might separate well
the three considered COUP classes, though, on their own, each of these variables fails to
achieve a low misclassification rate. To check this, we have carried out three classification
procedures on these data: quadratic discriminant analysis, the k-nearest neighbours (k-NN)
rule and model-based discriminant analysis, respectively implemented in the R packages
MASS ([28]), class ([28]) and mclust ([24]). We can choose k via cross-validation (CV).
Using the package caret ([19]), we have run 10-, 5-fold and leave-one-out (LOO) CV to get
insight into the value of 5 ≤ k ≤ 63 yielding the largest accuracy. In the case of Zolotarev
distance ζ̄2 the optimal value of k is 5 for the three CV procedures. For the Wasserstein
metric ω̄, the optimal k’s are 35, 33 and 39 for 10-, 5-fold and LOO CV, respectively. For
the Kolmogorov metric κ, the optimal k’s are 19, 17 and 21 (or 23 and 25) for 10-, 5-fold
and LOO CV, respectively. Since, for each metric, the optimal values of k are similar,
we have chosen k = 5 for k-NN with the ζ̄2 metric, k = 35 with ω̄ and k = 19 in the
case of κ. The percentage of correct classifications with the three discriminant methods
and the three distances appear in Table 1: they are all nearly the same and remarkably
high (around 90%). The normalized Zolotarev distance of the PIT to the exponential
distribution achieves the best correct classification rate, no matter which procedure is used.
Although the improvement over the usual Kolmogorov metric may not seem significant,
this is probably due to the low proportion of extragalactic sources in the sample. To
17

-1
-2

895

-4

-3

log(κ(Fn , Gµ̂ ))

1015

514
383

Extragalactic

Lightly obscured

Heavily obscured

-2

895
751
1015
961

-4

-3

log(ω̄(Fn , Gµ̂ ))

-1

0

(a)

1008
939
514

8
111

1199

Extragalactic

Lightly obscured

Heavily obscured

(b)
797

224
308

2

160
424
333

0
-2

751
895
961
1015

383,939
1071

-4

log(ζ̄2 (Fn , Gµ̂ ))

79,101,90

1304
111

Extragalactic

1182
1143
1269
1199

Lightly obscured

Heavily obscured

(c)
18
Figure 9: Analysis of COUP data. Boxplots of log(d(Fn , Gµ̂ )) for (a) d = κ, (b) d = ω̄
and (c) d = ζ̄2 .

0.8

1.0

1.0
0.8

0.6

Extragalactic
Asymptotic
Lightly-obs.
Heavily-obs.

0.0

0.0

0.2

0.2

0.4

0.4

0.6

Extragalactic
Asymptotic
Lightly-obs.
Heavily-obs.

0

1

2

3

4

0

x
(a)

1

2

3

4

5

6

x
(b)

Figure 10: Analysis of COUP data. Distribution √
function of δ∞ (d, G1 ) (as given in
Corollary 1) and empirical distribution functions of n d(Fn , Gµ̂ ) for (a) d = ω̄ and (b)
d = ζ̄2 .
gain a better insight of the advantages of the Zolotarev metric, in Table 2 we display the
confusion matrices of these classification procedures. If we look at the extragalactic class
(NM), using the Kolmogorov distance provides positive predictive values (PPVs) of 67%,
65% and 68%, for the quadratic, the k-NN and the mixture classification rules respectively.
With the Zolotarev metric the corresponding PPVs are 78%, 76% and 78%, respectively.
This significant increase in the classification accuracy was one of the motivations of this
work.
There could be concerns regarding the potential influence of the outlying sources and
their masking effect on the procedures employed in this section. For example, classical
discriminant methods (like the quadratic or model-based rules) are sensitive towards
outliers and may result in inaccurate parameter estimation or ill-posed problems. These
methods can be replaced by robustified classifiers that rely, e.g., on depth functions ([17])
or on robust estimates of the mixture parameters ([10]). However, robust versions of the
normalized distances of the PIT distributions to the exponential class are not adequate for
this problem. The reason is that outlying interarrival times are precisely the ones leading
to detection of flares, inconsistent with the usual behaviour of extragalactic radiation.
Of the three distances, κ, ω̄ and ζ̄2 , for the next diagram we have kept only the latter, as it
has the highest discriminant ability. Figure 12 displays the scatterplot of the logarithm of
MedEn in terms of the logarithm of the normalized Zolotarev distance to the exponential
distribution. At a glance, we see the three COUP classes separated in clusters (extragalactic:
high MedEn/low distance; heavily obscured: high energy/high distance; lightly obscured:
low energy/medium distance).

19

Table 1: Analysis of COUP data: Percentage of correct classifications based on the
logarithm of the median photon energy joint with the logarithm of the Kolmogorov
distance (first row), the normalized Wasserstein distance (second row) and the normalized
Zolotarev distance to the exponential distribution (last row).
Distance
Quadratic
Kolmogorov
88.99
Wasserstein
89.54
Zolotarev
90.18

k-NN
87.89
87.98
90.09

Model-based
89.08
88.99
90.64

Table 2: Analysis of COUP data: confusion matrices for classifications based on the
logarithm of the median photon energy and the logarithm of the Kolmogorov distance
(first row), the normalized Wasserstein distance (second row) and the normalized Zolotarev
distance (last row) to the exponential distribution. Numbers in bold font correspond to the
highest correct classification in each class. NM = non members, HO = Heavily obscured,
LO = Lightly obscured stars.
Classification rule
k-NN

20

Predicted

NM
HO
LO

NM
HO
LO

Predicted

NM
HO
LO

Actual
NM HO LO
50 26 1
23 292 27
0 55 616
Actual
NM HO LO
54 26 2
18 283 20
1 64 622
Actual
NM HO LO
57 15 3
16 317 33
0 41 608

Mixture

NM
HO
LO

Predicted

Predicted

NM
HO
LO

NM
HO
LO

Predicted

NM
HO
LO

Actual
NM HO LO
48 23 1
25 299 20
0 51 623
Actual
NM HO LO
55 21 2
18 297 18
0 55 624
Actual
NM HO LO
59 15 2
14 302 20
0 56 622

Predicted

Predicted

NM
HO
LO

Predicted

Quadratic

Predicted

Zolotarev Wasserstein Kolmogorov

Metric

NM
HO
LO

Actual
NM HO LO
48 22 1
25 308 28
0 43 615
Actual
NM HO LO
54 22 2
19 295 21
0 56 621
Actual
NM HO LO
60 15 2
13 307 21
0 51 621

Figure 11: Analysis of COUP data: Empirical distribution functions of the photon energies
rescaled to the interval [0,1].
Table 3: Analysis of COUP data: Posterior probabilities of membership for extragalactic
outliers.
Source
751
895
961
1015

Posterior probabilities (in percentage)
Extragalactic Heavily obs. Lightly obs.
0.0152
92.8997
7.08514
1.5713
98.4230
0.00576
0.6601
98.9164
0.42343
8.4486
91.5511
0.00003

In Figure 12, a simple visual inspection reveals that COUP sources 751, 895 and 961,
classified as extragalactic, have a higher probability of being heavily obscured stars. Indeed,
Table 3 displays the posterior probabilities of membership, derived from the quadratic
classification rule, for the misclassified extragalactic cases when these probabilities exceed
0.9. We see that all the sources in Table 3 are actually the largest outliers “detected”
by the normalized Zolotarev distance ζ̄2 in the extragalactic class (Figure 9(c)). This
emphasizes the information conveyed by ζ̄2 on the source class.

6

Conclusions

We have introduced normalized versions of two integral probability metrics, the Wasserstein
and Zolotarev distances, to quantify the discrepancy between photon interarrival times of
X-ray cosmic sources and the exponential distribution. The aim is to measure how different
the photon emitting source is from extragalactic radiation. The plug-in estimators of these
metrics show a good asymptotic behaviour. The analysis of more than one thousand X-ray
sources from the Chandra Orion Ultradeep Project with the proposed metrics reveals that
the information conveyed by photon interarrival times on the nature of each X-ray source
21

Source type
1.5

Extragalactic
Lightly obscured
Heavily obscured
1051

1.0

961
751

0.0

0.5

log(MedEn)

895

-4

-2

0

2

log(ζ̄2 (Fn , Gµ̂ ))

Figure 12: Analysis of COUP data: Scatterplot of the logarithm of the median energy
and the logarithm of the normalized Zolotarev distance to the exponential distribution.
The four outlying sources of the extragalactic class, identified by the Zolotarev distance,
are highlighted with filled circles.

22

is very well summarized by the normalized Zolotarev distance. We have further shown
that this metric, together only with the median energy of the X-ray source, yields a high
percentage of correct classifications of the sources into the classes previously provided by
astrophysicists. We remark that here we have only used two discriminating features while
the usual expert procedures in astronomy rely on many more. As a striking conclusion,
we have detected four sources, originally classified as extragalactic, which, as a matter of
fact, are very likely young stars in Orion Molecular Cloud 1.

Appendix
In this technical appendix, we collect the proofs of the results in Section 3.2. The main
ingredients are the following: first, we show that the sequences of stochastic processes
defined in (7) are equivalent in L1 to continuous functionals of the empirical process; then
we apply the central limit theorem (CLT) in suitable Banach spaces to find their weak
limits; finally, we use the continuity of the functional given in (6) to derive the asymptotic
distribution of δn (d, F ) in (4).
To begin, we recall that if the stochastic processes Pn and P take values in L1 a.s., it is said
that Pn converges in distribution to P in L1 if limn→∞ Ef (Pn ) = Ef (P), for all continuous
and bounded functions f : L1 −→ R. Note that if P and Pn are jointly measurable and
have almost all their trajectories in L1 , they can be identified with Borel-measurable
random elements in L1 . Therefore, the previous expectations are well-defined. In the
L1

following, we denote this weak convergence of probability measures in L1 by Pn −→w P.
An analogous definition can be given for the weak convergence in the weighted L1 space
defined by


Z
∞

W 1 :=

f ∈ L1 : kf kW 1 :=

(1 + t)|f (t)| dt < ∞ .
0

en be two stochastic processes with trajectories in L1 a.s. We say that Pn and
Let Pn and P
1
en are equivalent in L1 , denoted by Pn L∼ P
en , if kPn − P
en k1 −→P 0, where “−→P ” stands
P
1

L e
for convergence in probability. Roughly speaking, if Pn ∼ P
n , the two processes have the
1
L1
L1 e
1
e L
same asymptotic behavior in L because if Pn −→w P and Pn ∼ P
n , then Pn −→w P (see
for instance [?, Theorem 18.10]).

In
that is, En (t) :=
√ the sequel, En stands for the empirical process associated to X,
1
n(Fn (t) − F (t)), t ≥ 0, n ≥ 1. The asymptotic behavior of En in L and W 1 is collected
in the following lemma. Part (a) is a known result (see [?, Theorem 2.1]), while part (b)
can be found in [3, Lemma 2].
Lemma 1. We have that
L1

(a) En −→w BF if and only if X ∈ L2,1 .
W1

(b) En −→w BF if and only if X ∈ L4,2 .

23

Our first task is to find processes expressed as continuous functionals of En and equivalent
to Xd,n in L1 . To start, for t ≥ 0, we decompose Xd,n given in (7) as:
Xd,n = Ad,n + Bd,n + Cd,n ,

(15)

where
Aω,n := En ,

Bω,n :=

Aω̄,n := En /µ̂,
Z ∞
En ,
Aζ2 ,n (t) :=
t
Z ∞
1
Aζ̄2 ,n (t) := 2
En ,
µ̂ t

Bω̄,n :=

√
√

n(Gµ − Gµ̂ ),

Cω,n := 0,

n(Gµ − Gµ̂ )/µ̂,
Z ∞
√
(Gµ − Gµ̂ ),
Bζ2 ,n (t) := n
t
√ Z ∞
n
Bζ̄2 ,n (t) := 2
(Gµ − Gµ̂ ),
µ̂ t

√
Cω̄,n := gω̄ n(µ − µ̂)/µ̂,
Cζ2 ,n := 0,
Cζ̄2 ,n

√
n(µ2 − µ̂2 )
:= gζ̄2
,
µ̂2

with gω̄ and gζ̄2 defined in (8) and (9), respectively.
Lemma 2 provides equivalent expressions for the processes defined above.
Lemma 2. Let X be a positive random variable with mean µ > 0. For t ≥ 0, the following
assertions hold:
√
L e
e
n(µ̂ − µ)te−t/µ /µ2 .
(a) If X ∈ L4/3 , Bω,n ∼ B
ω,n , where Bω,n (t) :=
1

1

L e
(b) If X ∈ L2,1 , Aω̄,n ∼ A
ω̄,n := En /µ.

√
L e
e
(c) If X ∈ L4/3 , Bω̄,n ∼ B
n(µ̂ − µ)te−t/µ /µ3 .
ω̄,n , where Bω̄,n (t) :=
1

√
L e
(d) If X ∈ L4/3 , Cω̄,n ∼ C
n(µ − µ̂)gω̄ /µ.
ω̄,n :=
1

√
L e
e
(e) If X ∈ L4/3 , Bζ2 ,n ∼ B
n(µ̂ − µ)(1 + t/µ)e−t/µ .
ζ2 ,n , where Bζ2 ,n (t) :=
1

1

L e
e
(f ) If X ∈ L4,2 , Aζ̄2 ,n ∼ A
ζ̄2 ,n , where Aζ̄2 ,n (t) :=

R∞
t

En /µ2 .

√
L e
e
(g) If X ∈ L4/3 , Bζ̄2 ,n ∼ B
n(µ̂ − µ)(1 + t/µ)e−t/µ /µ2 .
ζ̄2 ,n , where Bζ̄2 ,n (t) :=
1

√
L e
n(µ − µ̂)2 gζ̄2 /µ.
(h) If X ∈ L2 , Cζ̄2 ,n ∼ C
ζ̄2 ,n :=
1

Proof. To show part (a), we use the mean value theorem twice to obtain
Z ∞
√
2
e ω,n k1 ≤ n(µ̂ − µ)
t|2 − t/ξt |e−t/ξt /ξt3 dt,
kBω,n − B

(16)

0

where ξt is a point between µ and µ̂. The integral in (16) is bounded by
Z ∞
max{µ, µ̂}2
t|2 − t/ξt |e−t/ξt /ξt3 dt ≤ 2(µ + µ̂)
→ 4/µ a.s.
min{µ, µ̂}4
0

24

(17)

Therefore, from (16)-(17), and by the Kolmogorov, Marcinkiewicz and Zygmund strong
law of large numbers (see, e.g., [18, Theorem 3.23]), we see that, whenever X ∈ L4/3 ,
e ω,n k1 → 0 a.s.
kBω,n − B
e ω̄,n k1 = kEn k1 |µ − µ̂|/(µµ̂). From Lemma 1 (a), we have
To see (b), we note that kAω̄,n − A
that kEn k1 −→d kBF k1 and the conclusion follows from the strong law of large numbers
and Slutsky’s theorem.
L1

Part (c) follows from (a), as it is straightforward to check that Bω̄,n ∼ Bω,n /µ, whenever
X ∈ L4/3 .
Part (d) is direct, whereas part (e) can be found in [3, Lemma 1]. The proof of part (f) is
similar to the one for (b) by using Lemma 1 (b).
1

e ζ̄ ,n L∼ Bζ ,n /µ2 . The
To show part (g), we observe that, from part (e), we have
that B
2
2
√
conclusion follows by checking that kBζ̄2 ,n − Bζ2 ,n /µ2 k1 = n(µ − µ̂)2 (1/µ + 1/µ̂)2 .
e ζ̄ ,n k1 = √n(µ − µ̂)2 (1/µ + 1/µ̂) ζ̄2 (X, Yµ ). As
Finally, it can be seen that kCζ̄2 ,n − C
2
ζ̄2 (X, Yµ ) < ∞ if and only if X ∈ L2 , we conclude that (h) is fulfilled.

The next corollary, which is a consequence of Lemma 2 and (15), shows that Xd,n are
equivalent in L1 to certain continuous functionals of the empirical process En .
Corollary 3. Let X be a positive random variable with mean µ > 0.
L1

(i) If X ∈ L4/3 , then Xω,n ∼ φω (En ), where φω : L1 → L1 is the linear operator
Z
t −t/µ ∞
φω (f, t) := f (t) − 2 e
f (x) dx, t ≥ 0.
µ
0
Moreover, kφω (f )k1 ≤ 2 kf k1 , and φω is therefore continuous.
L1

(ii) If X ∈ L2,1 , then Xω̄,n ∼ φω̄ (En ), where φω̄ : L1 → L1 is the linear operator



Z ∞
1
t −t/µ
φω̄ (f, t) :=
f (x) dx , t ≥ 0.
f (t) + gω̄ (t) − 2 e
µ
µ
0
Moreover, kφω̄ (f )k1 ≤ kf k1 (2 + ω̄(X, Yµ ))/µ, and φω̄ is therefore continuous.
L1

(iii) If X ∈ L4/3 , then Xζ2 ,n ∼ φζ2 (En ), where φζ2 : W 1 → L1 is the linear operator


Z ∞
Z ∞
t
−t/µ
φζ2 (f, t) :=
f (x) dx − 1 +
e
f (x) dx, t ≥ 0.
µ
t
0
Moreover, kφζ2 (f )k1 ≤ (1 + 2µ) kf kW 1 , and φζ2 is therefore continuous.
L1

(iv) If X ∈ L4,2 , then Xζ̄2 ,n ∼ φζ̄2 (En ), where φζ̄2 : W 1 → L1 is the linear operator
Z ∞



Z ∞

1
t
−t/µ
φζ̄2 (f, t) := 2
f (x) dx + 2 µ gζ̄2 (t) − 1 +
e
f (x) dx ,
µ
µ
t
0
for t ≥ 0. Moreover, kφζ̄2 (f )k1 ≤ kf kW 1 [1 + 2µ(1 + ζ̄2 (X, Yµ ))]/µ2 , and φζ̄2 is
therefore continuous.
25

We are now in condition to prove Theorems 1 and 2.
Proof of Theorem 1 Assume that (a) holds, i.e., X ∈ L2,1 . For d = ω or d = ω̄, from
L1

L1

Lemma 1 (a) and Corollary 3 (i) and (ii), we have that Xd,n ∼ φd (En ) −→w φd (BF ) = Xd,F ,
by the continuos mapping theorem. We conclude that (b) and (c) hold.
Conversely, let us assume that (b) is satisfied. By Corollary 3 (i), if X ∈ L4/3 , we obtain
L1

that φω (En ) −→P
w Xω,F . Observe now that φω (En ) can be rewritten as the normalized sum
n
−1/2
φω (En ) = n
i=1 Yω,i where Yω,1 , . . . , Yω,n are n independent copies of the zero mean
process
Yω (t) := P(X > t) − I{X>t} + (X − µ)te−t/µ /µ2 , t ≥ 0.
(18)
This means that the process Yω satisfies the CLT in L1 (and implies that Xω,F is a centered
Gaussian process), which is equivalent (see [1, p. 205]) to
Z ∞p
EYω (t)2 dt < ∞.
(19)
0

In particular, this implies that X ∈ L2 and denoting Z(t) := P(X > t) − I{X>t} (t ≥ 0),
from (18), (19), and by Minkowski inequality, we have that
Z ∞p
Z ∞p
2
EZ(t) dt ≤ σ +
EYω (t)2 dt < ∞,
0

0

where σ is the standard deviation of X. Last inequality amounts to (a) X ∈ L2,1 .
L1

For the proof that (c) implies (a), it is enough to note that Xω̄,n −→w Xω̄,F is equivalent
L1

to µ̂ Xω̄,n −→w µ Xω̄,F . The rest of the proof runs as with d = ω.



Proof of Theorem 2 To show that part (a) implies (b) and (c) it is enough to follow the
same steps as in the proof of the same implications in Theorem 1. We omit the details.
To finish, we will show that part (c) implies (a) (the remaining implication “(b) ⇒ (a)” is
simpler and similar). Let us assume that (c) is satisfied. In this situation, it is clear that
X ∈ L2 , as this integrability condition amounts to saying that the process Xζ̄2 ,n has its
L1

paths in L1 a.s. We have that µ̂2 Xζ̄2 ,n −→w µ2 Xζ̄2 ,F . Further, by Lemma 2, we conclude
L1

that µ̂2 Xζ̄2 ,n ∼ Zn , where
∞

√
En + h(t) n(µ − µ̂),

t ≥ 0,

h(t) := 2µ gζ̄2 (t) − (1 + t/µ)e−t/µ ,

t ≥ 0.

Z
Zn (t) :=
t

with
(20)
P
The process Zn can be rewritten as the normalized sum Zn = n−1/2 ni=1 Zζ̄2 ,i , where
Zζ̄2 ,1 , . . . , Zζ̄2 ,n are n independent copies of the zero mean process
Zζ̄2 (t) := E(X − t)+ − (X − t)+ + h(t) (µ − X).

26

Therefore, the process Zζ̄2 satisfies the CLT in L1 (and we see that Xζ̄2 ,F is a centered
Gaussian process). Using again [1, p. 205], we obtain that
Z ∞q
EZζ̄2 (t)2 dt < ∞.
(21)
0

Finally, by Minkowski inequality and Fubini theorem, we have that
Z ∞q
Z ∞q
2
2
E(X − t)+ dt ≤ EX /2 + σkhk1 +
EZζ̄2 (t)2 dt.
0

0

From R(20),
we can check that khk1 ≤ 2µ(1 + ζ̄2 (X, Yµ )) < ∞, and from (21), we conclude
∞p
that 0
E(X − t)2+ dt < ∞. This implies that X ∈ L4,2 because t2 P(X > 2t) ≤
2
E(X − t)+ .

Proof of Theorem 3 We have that δn (d, F ) = ρn (Xd,n , gd ), with ρn defined in (6). The
continuity of ρn (in L1 ) was analyzed in [6, Lemma 4], where it was shown that if fn → f
in L1 and g ∈ L1 , then, for I(g) := {t ≥ 0 : g(t) = 0},
Z
Z
lim ρn (fn , g) = ρ(f, g) :=
|f | +
f sgn(g).
(22)
n→∞

I(g)

I(g)c

Therefore, from Theorems 1 and 2, (22), and the extended continuous mapping theorem (see
[?, Theorem 1.11.1]), we conclude that δn (d, F ) = ρn (Xd,n , gd ) −→d ρ(Xd,F , gd ) = δ∞ (d, F ),
as n → ∞.

Proof of Corollary 2 From Theorem 3, we have that
Z ∞
δ∞ (d, F ) =
Xd,F (t) sgn(gd (t)) dt.
0

As Xd,F is a centered Gaussian process and gd is nonrandom, we conclude that δ∞ (d, F ) is
normally distributed.


Acknowledgements
The authors are grateful to three reviewers and the associate editor for their insightful
comments which have improved the presentation of the paper.

References
[1] Araujo, A. and Giné, E. (1980). The Central Limit Theorem for Real and Banach
Valued Random Variables. Wiley.
[2] Ascher, S. (1990). A survey of tests for exponentiality. Communications in StatisticsTheory and Methods, 19(5), 1811–1825.
[3] Baı́llo, A., Cárcamo, J. and Nieto, E. (2015). A test for convex dominance with respect
to the exponential class based on an L1 distance. IEEE Transactions on Reliability,
64, 71–82.
27

[4] del Barrio, E., Giné, E. and Matrán, C. (1999). Central limit theorems for the
Wasserstein distance between the empirical and the true distributions. The Annals of
Probability, 27, 1009–1071.
[5] Broos, P.S., Getman, K.V., Povich, M.S., Townsley, L.K., Feigelson, E.D. and Garmire,
G.P. (2011). A naive Bayes source classifier for X-ray sources. The Astrophysical
Journal Supplement Series 194, 4.
[6] Cárcamo, J. (2017). Integrated empirical processes in Lp with applications to estimate
probability metrics. Bernoulli, 23(4B), 3412–3436.
[7] Efron, B., and Tibshirani, R.J. (1993). An introduction to the bootstrap. CRC press.
[8] Feigelson, E.D. and Babu, G.J. (2012). Modern Statistical Methods for Astronomy.
With R Applications. Cambridge University Press.
[9] Feigelson, E.D., Getman, K., Townsley, L., Garmire, G., Preibisch, T., Grosso, N.,
Montmerle, T., Muench, A. and McCaughrean, M. (2005). Global X-ray properties of
the Orion Nebula region. The Astrophysical Journal Supplement Series, 160, 379–389.
[10] Garcı́a-Escudero, L.A., Gordaliza, A. and Mayo-Iscar, A. (2014). A constrained robust
proposal for mixture modeling avoiding spurious solutions. Advances in Data Analysis
and Classification, 8, 27–43.
[11] Getman, K.V., Feigelson, E.D., Grosso, N., McCaughrean, M.J., Micela, G., Broos, P.
Garmire, G. and Townsley, L. (2005a). Membership of the Orion Nebula population
from the Chandra Orion Ultradeep Project. The Astrophysical Journal Supplement
Series, 160, 353–378.
[12] Getman, K.V., Flaccomio, E., Broos, P.S., Grosso, N., Tsujimoto, M., Townsley, L.,
Garmire, G.P., Kastner, J., Li, J., Harnden, F.R., Wolk, S., Murray, S.S., Lada, C.J.,
Muench, A.A., McCaughrean, M.J., Meeus, G., Damiani, F., Micela, G., Sciortino,
S., Bally, J., Hillenbrand, L.A., Herbst, W., Preibisch, T. and Feigelson, E. (2005b).
Chandra Orion Ultradeep Project: Observations and source lists. The Astrophysical
Journal Supplement Series, 160, 319–352.
[13] Getman, K.V., Feigelson, E.D., Broos, P.S., Micela, G. and Garmire, G.P. (2008a).
X-ray flares in Orion young stars. I. Flare characteristics. The Astrophysical Journal,
688, 418–436.
[14] Getman, K.V., Feigelson, E.D., Micela, G., Jardine, M.M., Gregory, S.G. and Garmire,
G.P. (2008b). X-ray flares in Orion young stars. II. Flares, magnetospheres, and
protoplanetary disks. The Astrophysical Journal, 688, 437–455.
[15] Grafakos, L. (2014). Classical Fourier Analysis. Springer.
[16] Henze, N., and Meintanis, S.G. (2005). Recent and classical tests for exponentiality:
a partial review with comparisons. Metrika, 61(1), 29–45.

28

[17] Hubert, M., Rousseeuw, P. and Segaert, P. (2017). Multivariate and functional
classification using depth and distance. Advances in Data Analysis and Classification,
11, 445–466.
[18] Kallemberg, O. (1997). Foundations of Modern Probability. Springer.
[19] Kuhn, M. Contributions from J. Wing, S. Weston, A. Williams, C. Keefer,
A. Engelhardt, T. Cooper, Z. Mayer, B. Kenkel, the R Core Team, M. Benesty, R. Lescarbeau, A. Ziem, L. Scrucca, Y. Tang, C. Candan and T. Hunt.
(2016). caret: Classification and Regression Training. R package version 6.0-72.
https://CRAN.R-project.org/package=caret
[20] Ledoux, M. and Talagrand, M. (2002). Probability in Banach Spaces. Springer.
[21] Rachev, S.T., Klebanov, L., Stoyanov, S.V. and Fabozzi, F. (2013). The Methods of
Distances in the Theory of Probability and Statistics. Springer.
[22] Rachev, S.T., Stoyanov, S.V. and Fabozzi, F. (2011). A Probability Metrics Approach
to Financial Risk Measures. Wiley-Blackwell.
[23] Schulz, N.S. (2012). The Formation and Early Evolution of Stars: From Dust to Stars
and Planets. Second edition. Springer.
[24] Scrucca L., Fop M., Murphy T.B.and Raftery A.E. (2016). mclust 5: Clustering,
classification and density estimation using Gaussian finite mixture models. R Journal,
8, 289–317.
[25] Shaked, M. and Shanthikumar, J.G. (2006). Stochastic Orders. Springer.
[26] van der Vaart, A. (1998). Asymptotic Statistics. Cambridge University Press.
[27] van der Vaart, A.W. and Wellner, J.A. (1996). Weak Convergence and Empirical
Processes: with Applications to Statistics. Springer.
[28] Venables, W.N., and Ripley, B.D. (2002). Modern Applied Statistics with S. Springer,
New York.
[29] Wolk, S.J., Harnden, F.R., Jr, Flaccomio, E., Micela, G., Favata, F., Shang, H.
and Feigelson, E.D. (2005). Stellar activity on the young suns of Orion: COUP
observations of K5-7 pre-main-sequence stars. The Astrophysical Journal Supplement
Series 160, 423–449.

29

