Pierre-Yves Vandenbussche

Charles Teissèdre

Events Retrieval Using Enhanced Semantic Web Knowledge

Keywords: Events Retrieval, Semantic Web, Data Mashup

In this article, we present an experimental end user application to query DeRiVE 2011 challenge dataset in an innovative and intuitive manner. After enriching the dataset with external sources of information, it is indexed in a way that enables users to submit queries combining keywords, location and temporal anchor, in a single search field. The goal is to ease event retrieval providing a simple user interface to query and visualize events over time.

Introduction

While geolocation services have enjoyed strong progress, few initiatives take into consideration chrono-localization and temporal query processing for Information Retrieval over the Web [START_REF] Alonso | On the Value of Temporal Information in Information Retrieval[END_REF]. As Linked Open Data grows, things are changing, since more and more temporally anchored data is available. However, processing temporal data remains a challenge from (i) a modeling point of view, (ii) for data acquisition, (iii) as well as in terms of querying and navigating through it.

In this article, we address the last issue of querying and navigating through temporal data. We describe a system using the RDF data provided along with the DeRiVE 2011 challenge 1 . The dataset describes entertainment events related to music, such as advertisements for concerts or festivals. It also provides some information about agents involved in these events and about their location. The main objective of the system we present here is to hide data complexity and make it simple to query, providing a single search field as a first step in events retrieval. The goal is to make DeRiVE dataset temporally browsable. The considered use case consists in finding events occurring at a given period of time at a specific location.

After a brief overview of how temporal information is handled in the context of Information Retrieval over the Web, we will describe the way we processed the dataset to enrich it with external sources of information and to index it. We will then describe the final application to query and browse the dataset.

Temporal Information Retrieval over the Web

Retrieving temporal information over the Web of Content (i.e. HTML-based Web) and in the Web of Data (aka the Semantic Web) are two different issues, though they may converge on some points.

Temporal Search within the Web of Content. Major search engines currently offer few temporal search services. One such service is Google timeline feature 2 , which offers a way to visualize keywords frequency at different periods of time and to browse sentences where these keywords are associated with a date. However, temporal expressions are reduced to point in time with no duration extent, hence there is an important loss of information. Processing temporal information expressed in Web documents is a challenge from at least three different points of view: (i) modeling temporal references (models should be able to represent dates and intervals, but may also need to cope with approximate information (e.g. "by the end of the 13th century"), iterative occurrences (e.g. "every day from 10am to 8pm"), as well as deictics (e.g. "yesterday", "two months ago") and anaphorics (e.g. "the day before")); (ii) document annotation (it requires processing huge amount of documents with NLP techniques that necessarily have to deal with imperfect precision and recall rate) and (iii) relevancy ordering of the results from the temporal perspective (how to rank documents by relevance from the temporal perspective?).

Temporal Search within the Web of Data. While the modeling issue remains a difficulty, the acquisition process in this context is quite different, since the data to process is structured. Data acquisition however can be an issue as well. As for the querying process, the main querying language, SPARQL, allows filtering results in a timespan (i.e. intervals of well defined dates). This approach explains why generally only well defined temporal properties are effectively employed in LOD 3 .

The three Web sites that provided data for the challenge relies on this process: Upcoming Yahoo!, Last.fm and Eventful all propose similar approach to event retrieval. The main search scenario, with little variation depending on the Web site, follows this path: user has to provide a location, then a type of event (concert/festival), then eventually a musical genre, a date filter, etc. Such rich faceted search scenario is not made possible, though, with the DeRiVE challenge dataset, since no information on the type of event is provided.

Processing DeRiVE 2011 Dataset

The application we present here is an experimental retrieval engine with the goal to query and browse events temporally in the simplest way possible. It can be used both in the context of the Web of Content [START_REF] Teissèdre | Recherche d'information et temps linguistique : une heuristique pour calculer la pertinence des expressions calendaires[END_REF] and the Web of Data, as it relies on indexing process and NLP resources for temporal references extraction which can analyze either a query or Web documents. For the DeRiVE challenge, in order to get enough information to enable users to submit queries combining keywords, location and temporal information, we first had to enrich the dataset. The DeRiVE 2011 dataset is composed of 107.874 events and related knowledge. Knowledge is originating from Upcoming Yahoo! (12.15%), LastFm (53.04%) and Eventful (34.81%). It has been transformed by EventMedia [START_REF] Troncy | Linking events with media Proceedings of the 6th International Conference on Semantic Systems[END_REF]. The dataset is made of more than 1.800.000 statements. Temporal information consists in either single dates or intervals of dates.

Event geo-location augmentation. 98.794 events (91.58%) have latitude and longitude information. The first knowledge augmentation process concerns events' geolocation. It tries to fetch city, country and address information from coordinates, using Google and Yahoo! reverse geocoding API. In our application this geolocation information is used during query processing to cope with countries or cities. It is also used to propose a map visualization using Google maps API.

Event Image augmentation.

Images provide a simple way to ensure a pleasant way to experience event browsing. To associate images to events, we set up a strategy based on images information in the Semantic Web (via SPARQL queries on EventMedia and DBpedia) and on the Web (via Flickr API). As a result, at least one image was associated to 95.01% of events. SPARQL query example on EventMedia endpoint4 using event URI:

SELECT distinct ?imageURI ?image WHERE{ ?imageURI <http://linkedevents.org/ontology/illustrate> <http://data.linkedevents.org/event/dba9e034-fea0-4d01-ba4c-fb0515b89051>. {?imageURI <http://www.w3.org/ns/ma-ontlocator> ?image. } UNION{?imageURI <http://www.w3.org/ns/ma-ont#locator> ?image. } } Agent Information augmentation. Information about agents involved in an event is valuable for our application users. By enriching the dataset with Wikipedia links that point toward articles concerning these agents, users can further their search. We collected these links thanks to SPARQL queries on DBpedia endpoint. We have been able to find Wikipedia links for 25.22% of the agents. SPARQL query example on DBpedia endpoint5 using agent label:

SELECT distinct ?wikiLink WHERE { {?s <rdfs:label> "Bob Dylan"@en.} UNION{?s <http://xmlns.com/foaf/0.1/name> "Bob Dylan"@en.} {?s a <http://dbpedia.org/ontology/Person>.} UNION{?s a <http://dbpedia.org/ontology/Band>.} ?s <http://xmlns.com/foaf/0.1/page> ?wikiLink. }

An Experimental Temporal Search Engine to Retrieve Events

The system we have implemented is both a search engine and a tool to visualize and browse events6 . Temporal query relies on the search engine developed by [START_REF] Teissèdre | Recherche d'information et temps linguistique : une heuristique pour calculer la pertinence des expressions calendaires[END_REF]. The search engine is able to process queries with approximate temporal conditions like "around May 2007", even if this temporal expression does not exist in DeRiVE data. From the temporal perspective, event retrieval is based on an algorithm that calculates similarity scores between the temporal reference of the query and those that are associated to events. Based on Lucene and several modules to compute the dataset (see fig 1), the system can handle queries that may combine keywords, location and temporal information, such as "rock in London in August 2008" or "Bonn by the end of 2007". Temporal information, location information and event or agent description are indexed as different fields once the dataset is fully preprocessed. Queries are analyzed in such way that keywords, dates and location information are separated. Temporal data recognition in queries is performed thanks to an NLP module described in [START_REF] Teissèdre | Resources for Calendar Expressions Semantic Tagging and Temporal Navigation through Texts[END_REF]. The location extraction is performed thanks to a dictionary built during the indexing process: the dictionary contains cities and countries entities collected during the event geolocation enrichment process. Any other information that may appear in a query is considered as simple keywords, on which no semantic analysis is performed.

The events returned by the system are presented on a SIMILE timeline7 (see fig 2). The timeline on which results are displayed is fully browsable, which means that users can move over time: the system generates new queries on the fly as users move forward or backward in time.

Conclusion and Perspectives

The experimental application presented could be used as first step in events retrieval. Since the approach is generic and not tightly bound to DeRiVE dataset, it can be used in any other use case scenario where data is temporally anchored. If the dataset had contained information about musical genres, it could have been interesting to introduce faceted search with SolR tool, so that users could refine the results and eventually disambiguate query. Another interesting feature for possible improvement would be to synchronize a map for geolocation with the timeline, so as to present the results both in their temporal and geographic context.

Fig. 1 .

 1 Fig. 1. System's architecture.

Fig. 2 .

 2 Fig. 2. Screenshot of the UI for the query "rock in London in August 2008".

Dataset is available at: http://semanticweb.cs.vu.nl/derive2011/Challenge.html

URL for the query "revolution": http://bit.ly/relfGV

Despite Time Ontology[START_REF] Hobbs | An ontology of time for the semantic web[END_REF] capability to describe complex time knowledge representation, it is generally not used in all its' complexity.

URL: http://semantics.eurecom.fr/sparql

URL: http://dbpedia.org/sparql

The system can be tested at the following address: http://labs.mondeca.com/ChallengingTime/?locale=en&demo=eventMedia

URL: http://www.simile-widgets.org/timeline/

Acknowledgments. This project is partially granted by Datalift ANR project (ANR-10-CORD-009) and Chronolines ANR project (ANR-10-CORD-010).


Frank A Cowell
email: f.cowell@lse.ac.uk

Russell Davidson
email: russell.davidson@mcgill.ca

Emmanuel Flachaire
email: emmanuel.flachaire@univ-cezanne.fr

Goodness of Fit: an axiomatic approach

Keywords: Goodness of fit, axiomatic approach, measures of divergence, parametric bootstrap JEL codes: D63, C10

An axiomatic approach is used to develop a one-parameter family of measures of divergence between distributions. These measures can be used to perform goodness-of-fit tests with good statistical properties. Asymptotic theory shows that the test statistics have well-defined limiting distributions which are however analytically intractable. A parametric bootstrap procedure is proposed for implementation of the tests. The procedure is shown to work very well in a set of simulation experiments, and to compare favourably with other commonly used goodness-of-fit tests. By varying the parameter of the statistic, one can obtain information on how the distribution that generated a sample diverges from the target family of distributions when the true distribution does not belong to that family. An empirical application analyses a UK income data set.

Introduction

In this paper, we propose a one-parameter family of statistics that can be used to test whether an IID sample was drawn from a member of a parametric family of distributions. In this sense, the statistics can be used for a goodnessof-fit test. By varying the parameter of the family, denoted by α, a range of statistics is obtained and, when the null hypothesis that the observed data were indeed generated by a member of the family of distributions is false, the different statistics can provide valuable information about the nature of the divergence between the unknown true data-generating process (DGP) and the target family.

We do not seek just to add to the collection of convenient goodness-of-fit statistics based on the empirical distribution function (EDF) of the sample. Our approach here is to motivate the goodness-of-fit criterion in the same sort of way as is commonly done with other measurement problems in economics and econometrics. 1 The role of axiomatisation is central. We invoke a relatively small number of axioms to capture the idea of divergence of one distribution from another using an informational structure that is common in studies of income mobility. From this divergence concept one immediately obtains a class of goodness-of-fit measures that inherit the principles embodied in the axioms. As it happens, the measures in this class also have a natural and attractive interpretation.

In order to be used for testing purposes, the goodness-of-fit statistics should have a distribution under the null that is known or can be simulated. Asymptotic theory shows that the null distribution of the members of the family of statistics is independent of the parameter α, although that is certainly not true in finite samples. We show that the asymptotic distribution (as the sample size tends to infinity) exists, although it is not analytically tractable. However, its existence serves as an asymptotic justification for the use of a parametric bootstrap procedure for inference.

In addition to a set of simulation experiments designed to uncover the size and power properties of bootstrap tests based on our proposed family of statistics, we analyse a UK data set on households with below-average incomes, and show that we can derive a stronger conclusion by use of our tests than with most commonly used goodness-of-fit tests.

The paper is organised as follows. Section 2 sets out the formal framework and establishes a series of results that characterise the required class of measures. Section 3 derives the distribution of the members of this new class. Section 4 examines the performance of the goodness-of-fit criteria in practice, and section 5 concludes.

Axiomatic foundation

Our approach is to characterise a measure of aggregate discrepancy between two empirical income distributions and then to apply the discrepancy measure as a goodness-of-fit criterion.

Representation of the problem

We adopt a structure that is often applied in the income-mobility literature. Let there be an ordered set of n income classes; each class i is associated with income level x i where x i < x i+1 , i = 1, 2, ..., n -1. Let p i ≥ 0 be the size of class i, i = 1, 2, ..., n which could be an integer in the case of finite populations or a real number in the case of a continuum of persons. We will work with the associated cumulative mass u i = i j=1 p j , i = 1, 2, ..., n. The set of distributions is given by U

:= u| u ∈ R n + , u 1 ≤ u 2 ≤ ... ≤ u n .
The aggregate discrepancy measurement problem can be characterised as the relationship between two cumulative-mass vectors u, v ∈ U ; an alternative equivalent approach is to work with z : = (z 1 , z 2 , ..., z n ), where each z i is the ordered pair (u i , v i ), i = 1, ..., n and belongs to a set Z, which we will take to be a connected subset of R + × R + . The problem focuses on the discrepancies between the u-values and the v-values. To capture this we introduce a discrepancy function d :

Z → R such that d (z i ) is strictly increasing in |u i -v i |. Write the vector of discrepancies as d (z) := (d (z 1 ) , ..., d (z n )) .
The problem can then be approached in two steps.

1. We represent the problem as one of characterising a weak ordering2 on

Z n := Z × Z × ... × Z n .
where, for any z, z ∈ Z n the statement "z z " should be read as "the pairs in z constitute at least as good a fit according to as the pairs in z ." From we may derive the antisymmetric part and symmetric part ∼ of the ordering.3 2. We use the function representing to generate an aggregate discrepancy index.

In the first stage of step 1 we introduce some properties for , many of which correspond to those used in choice theory and in welfare economics.

Basic structure

Axiom 1 (Continuity) is continuous on Z n . Axiom 2 (Monotonicity) If z, z ∈ Z n differ only in their ith component then d (u i , v i ) < d (u i , v i ) ⇐⇒ z z .
For any z ∈ Z n denote by z (ζ, i) the member of Z n formed by replacing the ith component of z by ζ ∈ Z.

Axiom 3 (Independence) For z, z ∈ Z n such that: z ∼ z and z i = z i for some i then z (ζ, i) ∼ z (ζ, i) for all ζ ∈ [z i-1 , z i+1 ] ∩ z i-1 , z i+1 .
If z and z are equivalent in terms of overall discrepancy and the fit in class i is the same in the two cases then a local variation in component i simultaneously in z and z has no overall effect.

Axiom 4 (Perfect local fit) Let z, z ∈ Z n be such that, for some i and j,

u i = v i , u j = v j , u i = u i + δ, v i = v i + δ, u j = u j -δ, v j = v j -δ and, for all k = i, j, u k = u k , v k = v k . Then z ∼ z .
The principle states that if there is a perfect fit in two classes then moving u-mass and v-mass simultaneously from one class to the other has no effect on the overall discrepancy.

Theorem 1 Given Axioms 1 to 4 (a)

is representable by the continuous function given by

n i=1 φ i (z i ) , ∀z ∈ Z n (1)
where, for each i, φ i : Z → R is a continuous function that is strictly decreasing in |u i -v i | and (b)

φ i (u, u) = a i + b i u.
(2)

Proof. Axioms 1 to 4 imply that can be represented by a continuous function Φ : Z n → R that is increasing in |u i -v i |, i = 1, ..., n. Using Axiom 3 part (a) of the result follows from Theorem 5.3 of [START_REF] Fishburn | Utility Theory for Decision Making[END_REF]. Now take z and z in as specified in Axiom 4. Using (1) it is clear that z ∼ z if and only if

φ i (u i + δ, u i + δ) -φ i (u i , u i ) -φ j (u j + δ, u j + δ) + φ j (u j + δ, u j + δ) = 0
which can only be true if

φ i (u i + δ, u i + δ) -φ i (u i , u i ) = f (δ)
for arbitrary u i and δ. This is a standard Pexider equation and its solution implies (2).

Corollary 1 Since is an ordering it is also representable by

φ n i=1 φ i (z i ) (3)
where φ i is defined as in ( 1), ( 2) and φ : R → R, continuous and strictly monotonic increasing.

This additive structure means that we can proceed to evaluate the aggregate discrepancy problem one income class at a time. The following axiom imposes a very weak structural requirement, namely that the ordering remains unchanged by some uniform scale change to both u-values and vvalues simultaneously. As Theorem 2 shows it is enough to induce a rather specific structure on the function representing .

Axiom 5 (Population scale irrelevance) For any z, z ∈ Z n such that z ∼ z , tz ∼ tz for all t > 0.

Theorem 2 Given Axioms 1 to 5 is representable by

φ n i=1 u i h i u i v i (4)
where h i is a real-valued function.

Proof. Using the function Φ introduced in the proof of Theorem 1, Axiom 5 implies

Φ (z) = Φ (z ) Φ (tz) = Φ (tz )
and so, since this has to be true for arbitrary z, z we have

Φ (tz) Φ (z) = Φ (tz ) Φ (z ) = ψ (t)
where ψ is a continuous function R → R. Hence, using the φ i given in (1), we have for all z:

φ i (tz i ) = ψ (t) φ i (z i ) i = 1, ..., n.
or, equivalently

φ i (tu i , tv i ) = ψ (t) φ i (u i , v i ) , i = 1, ..., n. (5)
So, in view of [START_REF] Aczél | Functional Equations in Several Variables[END_REF], page 346 there must exist c ∈ R and a function h i : R + → R such that

φ i (u i , v i ) = u c i h i u i v i . (6)
From ( 2) and ( 6) it is clear that

φ i (u i , u i ) = u c i h i (1) = a i + b i u i , (7)
which implies c = 1 for non-trivial cases. Putting (6) with c = 1 into (3) gives the result. The function h i in Theorem 2 is arbitrary and it is useful to impose more structure. This is done in Section 2.3.

Mass discrepancy and goodness-of-fit

We now focus on the way in which one compares the (u, v) discrepancies in different parts of the distribution. The form of (4) suggests that discrepancy should be characterised terms of proportional differences:

d (z i ) = max u i v i , v i u i .
This is the form for d that we will assume from this point onwards. We also introduce:

Axiom 6 (Discrepancy scale irrelevance) Suppose there are z 0 , z 0 ∈ Z n such that z 0 ∼ z 0 . Then for all t > 0 and z, z such that d (z) = td (z 0 ) and

d (z ) = td (z 0 ): z ∼ z .
The principle states this. Suppose we have two distributional fits z 0 and z 0 that are regarded as equivalent under . Then scale up (or down) all the mass discrepancies in z 0 and z 0 by the same factor t. The resulting pair of distributional fits z and z will also be equivalent.4

Theorem 3 Given Axioms 1 to 6 is representable by

Φ (z) = φ n i=1 u α i v 1-α i ( 8
)
where α = 1 is a constant.5

Proof. Take the special case where, in distribution z 0 the discrepancy takes the same value r in all n classes.

If (u i , v i ) represents a typical compo- nent in z 0 then z 0 ∼ z 0 implies r = ψ n i=1 u i h i u i v i ( 9
)
where ψ is the solution in r to

n i=1 u i h i u i v i = n i=1 u i h i (r) . (10)
In (10) we can take the u i as fixed weights. Using Axiom 6 in (9) requires

tr = ψ n i=1 u i h i t u i v i , for all t > 0. (11)
Using (10) we have

n i=1 u i h i tψ n i=1 u i h i u i v i = n i=1 u i h i t u i v i (12)
Define variables q i and functions ψ i , i = 1, ..., n such that

q i := u i h i u i v i , (13)
u i v i = ψ i (q i ) . (14)
From ( 12) -( 14) we obtain

n i=1 u i h i tψ n i=1 q i = n i=1 u i h i (tψ i (q i )) . ( 15
)
Also define n + 1 functions θ i :

θ 0 (q, t) := n i=1 u i h i (tψ (q)) , (16)
θ i (q, t) := u i h i (tψ i (q)) , i = 1, ..., n. (17)
From ( 15) -( 17) we obtain

θ 0 n i=1 q i , t = n i=1 θ i (q i , t) ,
which implies (Aczél 1966, p. 142)

θ i (q, t) = b i (t) + B (t) q, i = 0, 1, ..., n where b 0 (t) = n i=1 b i (t) .
Therefore:

h i t u i v i = b i (t) u i + B (t) h i u i v i , i = 1, ..., n, (18)
the solution to which is

h i (v) = β i v α-1 + γ i , α = 1, β i log v + γ i α = 1, ( 19
)
where β i > 0 is an arbitrary positive number -see [START_REF] Eichhorn | Functional Equations in Economics. Reading Massachusetts[END_REF], Theorem 2.7.3. Substituting for h i (•) from ( 19) into (2) for the case where β i is the same for all i gives the result.

Aggregate discrepancy index

Theorem 3 provides some of the essential structure of an aggregate discrepancy index. To make further progress consider the behaviour of the index over the subset of Z consisting of all distributions with given total mass:

Z (ū, v) := z ∈Z| n i=1 z i = (ū, v) .
Clearly, for all z ∈Z (ū, v) the discrepancy index must take the form

Φ (z) = φ n i=1 u α i v 1-α i ; ū, v , (20)
where ū, v are parameters of the function φ that is the counterpart of φ in (8). It is reasonable to require that Φ (z) should take the value zero when z represents a "perfect fit," so that there is no discrepancy between the udistribution and the v-distribution; but there is a slight ambiguity about what the meaning of this requirement is.

A narrow interpretation of zero discrepancy is that

v i = u i , i = 1, .., n, in which case the form (20) requires that φ n i=1 u i ; ū, ū = 0, ( 21
)
in other words we have the restriction φ (ū; ū, ū) = 0. However, this restriction does not actually impose much additional structure. By contrast, suppose we take a broader interpretation of zero discrepancy, namely that if the v-distribution is obtained rescaling each component in the u-distribution by a factor k > 0 then there is no discrepancy; in other words suppose we say that the total mass is not relevant in the evaluation of discrepancy only the relative frequencies in each class. This interpretation requires that, if

v i = ku i , i = 1, .., n (where k = v/ū) then, from (20), we have φ k 1-α n i=1 u i ; ū, v = 0 (22) which implies φ ūα v1-α ; ū, v = 0. ( 23
)
This can only be true for all α if φ in (8) and φ in (20) can be written in the form

ψ n i=1 u i µ u α v i µ v 1-α , (24)
where

µ u := n -1 n i=1 u i and µ v := n -1 n i=1 v i .
A suitable cardinalisation of ( 24) gives the aggregate discrepancy measure

G α := 1 α(α -1) n i=1 u i µ u α v i µ v 1-α -1 , α ∈ R, α = 0, 1 (25)
where we have the following limiting forms for the cases α = 0 and α = 1, respectively

G 0 = - n i=1 v i µ v log u i µ u v i µ v , (26)
G 1 = n i=1 u i µ u log u i µ u v i µ v . ( 27
)
Expressions ( 25)-( 27) constitute a family of aggregate discrepancy measures where an individual family member is characterised by choice of α: a high positive α produces an index that is particularly sensitive to discrepancies where u exceeds v and a negative α yields an index that is sensitive discrepancies where v exceeds u.6

Goodness of fit

Our approach to the goodness-of-fit problem is to use the index constructed in section 2.4 to quantify the aggregate discrepancy between an empirical

1 α(α -1) Y F u (y) µ u α F v (y) µ v 1-α dy -1 , - Y F v (y) µ v log F u (y) µ u F v (y) µ v dy, and Y F u (y) µ u log F u (y) µ u F v (y) µ v dy.
Clearly there is a family resemblance to the [START_REF] Kullback | On Information and Sufficiency[END_REF] measure of relative entropy or divergence measure of

f 2 from f 1 Y f 1 log f 2 f 1 dy but with densities f replaced by cumulative distributions F .
distribution and a model. Given a set of observations {x 1 , x 2 , ..., x n } the empirical distribution function (EDF) is

F (x) = 1 n n i=1 ι x (i) ≤ x ,
where x (i) denotes the ith smallest observation and ι is an indicator function such that ι (S) = 1 if statement S is true and ι (S) = 0 otherwise. Denote the proposed model distribution by F (•; θ), where θ is a set of parameters and let

u i = F x (i) ; θ , i = 1, .., n v i = F x (i) = i n , i = 1, .., n;
then u i is a set of non-decreasing population proportions generated by the model from the n ordered observations. As before write µ u for the mean value of the u; observe that

µ v = 1 n n i=1 v i = n i=1 i n 2 = n + 1 2n .
Using ( 25)-( 27) we then find that we have a family of goodness-of-fit statistics

G α F, F = 1 α(α -1) n i=1 u i µ u α 2i n + 1 1-α -1 , (28)
where α ∈ R \ {0, 1} is a parameter. 7

Inference

If the parametric family F (•, θ) is replaced by a single distribution F , then the u i become just F (x (i) ), and therefore have the same distribution as the order statistics of a sample of size n drawn from the uniform U(0,1) distribution.

The statistic G α (F, F ) in ( 28) is random only through the u i , and so, for given α and n, it has a fixed distribution, independent of F . Further, as n → ∞, the distribution converges to a limiting distribution that does not depend on α.

7 In the cases α = 0 and α = 1 we have, respectively, that

G 0 F, F = - n i=1 2i n+1 log [n+1]ui 2iµ u and G 1 F, F = n i=1 ui µ u log [n+1]ui 2iµ u .
To see this, we make use of a result concerning the empirical quantile process; see [START_REF] Van Der Vaart | Weak Convergence and Empirical Processes[END_REF], example 3.9.24. Let F be a distribution function with continuous positive derivative f defined on a compact support. Let Fn be the empirical distribution function of an IID sample drawn from F , and let Q(p) = F -1 (p) and Qn (p) = F -1 n (p), p ∈ [0, 1], be the corresponding quantile functions. Since F is a discrete distribution, Qn (p) is just the order statistic indexed by np of the sample. Here x denotes the smallest integer not less than x. Then

√ n Qn (p) -Q(p) - B • F Q(p) f Q(p) , (29)
where the notation means that the left-hand side, considered as a stochastic process defined on [0, 1], converges weakly to the distribution of the righthand size, where f is the density of distribution F , and where B(p) is a standard Brownian bridge, that is, a Gaussian process with covariance function cov B(t), B(s) = min(s, t) -st.

The U(0,1) distribution certainly has compact support [0, 1], and its density is constant and equal to 1 on that interval. The result (29) in this case reduces to

√ n u np -p B(p). ( 30
)
We will be chiefly interested in the arguments t i defined as i/(n + 1), i = 1, . . . , n. Then we see that

√ n(u i -t i ) B(t i ). ( 31
)
This result expresses the asymptotic joint distribution of the uniform order statistics. Note that E(u i ) = t i . Write u i = t i +z i , where E(z i ) = 0. From (30), we see that the variance of

n 1/2 z i is t i (1 -t i ) plus a term that vanishes as n → ∞ 8 . Thus z i = O p (n -1/2 )
as n → ∞. We express the statistic G α (F, F ), under the null hypothesis that the u i do indeed have the joint distribution of the uniform order statistics, replacing u i by t i + z i and discarding terms that tend to 0 as n → ∞. We see that

G α (F, F ) = 1 α(α -1) 1 µ α u (1/2) 1-α n i=1 t i 1 + z i t i α -µ α u (1/2) 1-α . (32)
Now, by Taylor's theorem,

t i 1 + z i t i α = t i + αz i + α(α -1) 2
z 2 i t i + α(α -1)(α -2) 6 (θ i z i ) 3 t 2 i , (33)
where 0 ≤ θ i ≤ 1, i = 1, . . . , n, and so

n i=1 t i 1 + z i t i α = n 2 + nαz + α(α -1) 2 n i=1 z 2 i t i + o p (1), ( 34
)
where z is the mean of the z i , since it can be shown that the sum over i of the last term on the right-hand side of (33) if o p (1). Here, we have made use of the fact that n i=1 t i = (n + 1) -1 n i=1 i = n/2. By definition,

µ u = n -1 n i=1 u i = 1 2 + n -1 n i=1 z i = 1 2 + z. It follows that µ α u (1/2) 1-α = 1 2 (1 + 2z) α .
Using Taylor's theorem once more, we see that

µ α u (1/2) 1-α = 1 2 1 + 2αz + 2α(α -1)z 2 + 4α(α -1)(α -2) 3 (θ µ z) 3 , ( 35
)
with 0 ≤ θ µ ≤ 1. Now z is the estimation error made by estimating 1/2 by µ u , and so it is O p (n -1/2 ). The last term above is thus of order n -3/2 in probability. Putting together equations ( 34) and (35) gives

n i=1 t i 1 + z i t i α -µ α u 1 2 -α = α(α -1) 2 n i=1 z 2 i t i -2nz 2 + o p (1),
and so from (32) we arrive at the result

G α (F, F ) = n i=1 z 2 i t i -2nz 2 + o p (1). ( 36
)
It is striking that the leading-order term in (36) does not depend on α. For finite n, G α does of course depend on α. Simulation shows that, for n even as small as 10, the distributions of G α and of the leading term in (36) are very close indeed for α = 2, but that, for n even as large as 10,000, the distributions are noticeably different for values of α far enough removed from 2. The reason for this phenomenon is of course the factor of α -2 in the remainder terms in (33) and ( 35).

If the limiting asymptotic distribution of G α exists, it is the same as that of the approximation in (36), and, if the latter exists, it is the distribution of the limiting random variable obtained by replacing z i by n -1/2 B(t i ) (see ( 31)) and then letting n tend to infinity. For z first, we have

n 1/2 z = n -1/2 n i=1 z i = d n -1 n i=1 B(t i ) 1 0 B(t) dt. ( 37
)
Above, the symbol = d signifies equality in distribution, and the last step follows on noting that the second last expression is a Riemann sum that approximates the integral.

Similarly, we see that

n i=1 z 2 i /t i = d n -1 n i=1 B 2 (t i ) t i 1 0 B 2 (t) t dt. ( 38
)
From ( 37) and ( 38), we see that the limiting distribution of G α is that of

1 0 B 2 (t) t dt -2 1 0 B(t) dt 2 . ( 39
)
The denominator of t in the first integral may lead one to suppose that the integral may diverge with positive probability. However, notice that the expectation of the integral is

1 0 1 t EB 2 (t) dt = 1 0 (1 -t) dt = 1 2 .
A longer calculation shows that the second moment of the integral is also finite, so that the integral is finite in mean square, and so also in probability.

We conclude that the limiting distribution of G α exists, is independent of α, and is equal to the distribution of (39). We now turn to the more interesting case in which F does depend on a vector θ of parameters. The quantities u i are now given by u i = F (x (i) , θ), where θ is assumed to be a root-n consistent estimator of θ. If θ is the true parameter vector, then we can write x (i) = Q(v i , θ), where Q(•, θ) is the quantile function inverse to the cumulative distribution function (CDF) F (•, θ), and the v i have the distribution of the uniform order statistics. Then we have u i = F (Q(v i , θ), θ), and

µ u = n -1 n i=1 F Q(v i , θ), θ .
The statistic (28) becomes

G α = 1 α(α -1) 1 µ α u (1/2) 1-α n i=1 F Q(v i , θ), θ α t 1-α i -µ α u (1/2) 1-α .
Let p(x, θ) be the gradient of F with respect to θ, and define g(v, θ) to be p Q(v, θ), θ . As before, we let z i = v i -t i . Then a short Taylor expansion gives the approximation

F Q(v i , θ), θ = t i + z i + g (t i , θ))s(θ) + O p (n -1 ),
where s(θ) = θ -θ is the estimation error, and is of order n -1/2 . To leading order asymptotically, a calculation exactly like that leading to (36) gives

G α = n i=1 z i + g (t i , θ)s(θ) 2 t i -2 n -1/2 n i=1 z i + g (t i , θ)s(θ) 2 + o p (1).
(40) This asymptotic expression depends explicitly on θ, and also on the estimator θ that is used. In order to show that there does exist a limiting distribution for (40), we suppose that there exists a (vector) function h(x, θ) such that

n 1/2 ( θ -θ) = n 1/2 s(θ) = n -1/2 n i=1 h(x i , θ) + o p (1), ( 41
)
where E θ h(x, θ) = 0. The function h exists straightforwardly for most commonly used estimators, including maximum likelihood and least squares. Our sample is supposed to be IID, and so in (41) we can sum over the order statistics x (i) . Then a short Taylor expansion gives

n 1/2 s(θ) = n -1/2 n n=1 h Q(v i , θ), θ + o p (1) = n -1/2 n i-1 h Q(t i + z i , θ), θ + o p (1) = n -1/2 n i=1 h Q(t i , θ), θ + h Q(t i , θ), θ f Q(t i , θ), θ z i + o p (1), (42)
where f (x, θ) is the density that corresponds to F (x, θ) and h is the derivative of h with respect to its first argument. The integration over the whole real line means in fact integration over the support of the distribution F . Now, again by use of an argument based on a Riemann sum, we see that

n -1 n i=1 h Q(t i , θ), θ = 1 0 h Q(t, θ), θ dt + O(n -1 ) = ∞ -∞ h(x, θ) dF (x, θ) + O(n -1 ) = O(n -1 ), because E θ h(x, θ) = 0.
Thus the first term in the sum in ( 42) is O(n -1/2 ) and can be ignored for the asymptotic distribution. For the second term, we replace z i as before by n -1/2 B(t i ) to get

n 1/2 s(θ) 1 0 h Q(t, θ), θ f Q(t, θ), θ B(t) dt = ∞ -∞ h (x, θ)B F (x, θ) dx, ( 43
)
where for the last step we make the change of variables x = Q(t, θ), and note that dF (x, θ) = f (x, θ) dx.

Next consider the sum

n -1/2 n i=1 z i + g (t i , θ)s(θ)
that appears in (40). By the definition of g, g(t i , θ) = p Q(t i , θ), θ . Hence, with error of order n -1 , we have

n -1 n i=1 g(t i , θ) = n -1 n i=1 p Q(t i , θ), θ = 1 0 p Q(t, θ), θ dt = ∞ -∞ p(x, θ) dF (x, θ) ≡ P (θ).
Using (43), we have

n -1/2 n i=1 g (t i , θ)s(θ) P (θ) ∞ -∞ h (x, θ)B F (x, θ) dx,
and so

n -1/2 n i=1 z i +g (t i , θ)s(θ) 1 0 B(t) dt+P (θ) ∞ -∞ h (x, θ)B F (x, θ) dx. ( 44
)
Finally, we consider the first sum in (40). By arguments similar to those used above, we see that

n i=1 z i + g (t i , θ)s(θ) 2 t i 1 0 1 t B(t) + p Q(t, θ), θ × (45) ∞ -∞ h (x, θ)B F (x, θ) dx 2 dt.
So as to be sure that the integral converges with probability 1, we have to show that the non-random integrals

1 0 p Q(t, θ), θ t dt and 1 0 p 2 Q(t, θ), θ t dt are finite. Observe that 1 0 p Q(t, θ), θ t dt = ∞ -∞ p(x, θ) F (x, θ) dF (x, θ) = ∞ -∞ D θ log F (x, θ) dF (x, θ),
where D θ is the operator that takes the gradient of its operand with respect to θ. Similarly,

1 0 p 2 Q(t, θ), θ t dt = ∞ -∞ D θ log F (x, θ) 2 F (x, θ) dF (x, θ),
Clearly, it is enough to require that D θ log(F (x, θ) should be bounded for all x in the support of F (•, θ). It is worthy of note that this condition is not satisfied if varying θ causes the support of the distribution to change.

Under the condition just stated, the results ( 45) and ( 44) establish the existence of the limiting distribution of G α . In general, this distribution depends on the parameter vector θ, and so, in general, G α is not asymptotically pivotal with respect to the parametric family represented by the distributions F (•, θ). However, if the family can be interpreted as a location-scale family, then it is not difficult to check that, if θ is the maximum-likelihood estimator, the random quantity g (t i , θ)s(θ) does not in fact depend on θ. In addition, it turns out that the lognormal family also has this property. It would be interesting to see how common the property is, since, when it holds, the bootstrap benefits from an asymptotic refinement. But, even when it does not, the existence of the asymptotic distribution provides an asymptotic justification for the bootstrap.

It may be useful to give the details here of the bootstrap procedure used in the following section in order to perform goodness-of-fit tests, in the context both of simulations and of an application with real data. It is a parametric bootstrap procedure; see for instance [START_REF] Horowitz | Bootstrap Methods in Econometrics: Theory and Numerical Performance[END_REF] or [START_REF] Davidson | Bootstrap Methods in Econometrics[END_REF]. Estimates θ of the parameters of the family F (•, θ) are first obtained, after which the statistic of interest τ is computed, whether it is (28) for a chosen value of α or one of the other statistics studied in the next section. Bootstrap samples of the same size as the original data sample are drawn from the estimated distribution F (•, θ). Note that this is not a resampling procedure. For each of a suitable number B of bootstrap samples, parameter estimates θ * j , j = 1, . . . , B, are obtained using the same estimation procedure as with the original data, and the bootstrap statistic τ * j computed, also exactly as with the original data, but with F (•, θ * j ) as the target distribution. Then a bootstrap P value is obtained as the proportion of the τ * j that are more extreme than τ , that is, greater than τ for statistics like (28) which reject for large values. For well-known reasons -see [START_REF] Davison | Bootstrap Methods and their Application[END_REF] or [START_REF] Davidson | Bootstrap Tests: How Many Bootstraps[END_REF] -the number B should be chosen so that (B + 1)/100 is an integer. In the sequel, we set B = 999 unless otherwise stated. This computation of the P value can be used to test the fit of any parametric family of distributions.

Simulations and Application

We now turn to the way the new class of goodness-of-fit statistics performs in practice. In this section, we first study the finite sample properties of our G α test statistic and those of several standard measures: in particular we examine the comparative performance of the [START_REF] Anderson | Asymptotic Theory of Certain 'Goodness-of-Fit' Criteria based on Stochastic Processes[END_REF] statistic (AD), the Cramér-von-Mises statistic given by

CVM = ∞ -∞ F (x) -F (x, θ) 2 dF (x, θ), the Kolmogorov-Smirnov statistic KS = sup x | F (x) -F (x, θ)|,
and the Pearson chi-square (P) goodness-of-fit statistic

P = m i=1 (O i -E i ) 2 /E i ,
where O i is the observed percentage in the i th histogram interval, E i is the expected percentage in the i th histogram interval and m is the number of histogram intervals.9 Then, we provide an application using a UK data set on income distribution.

Tests for Normality

Consider the application of the G α statistic to the problem of providing a test for normality. It is clear from expression (28) that different members of the G α family will be sensitive to different types of divergence of the EDF of the sample data from the model F . We take as an example two cases in which the data come from a Beta distribution, and we attempt to test the hypothesis that the data are normally distributed. distribution is skewed to the left and Beta(2,5) is skewed to the right, while the normal distribution is of course unskewed. As can be deduced from ( 28), in the first case the G α statistic decreases as α increases, whereas in the second case it increases with α. 1, which shows normality tests with G α based on single samples of 1000 observations each drawn from the Beta(5,2) and from the Beta(2,5) distributions. Additional results are provided in the table with data generated by Student's t distribution with four degrees of freedom, denoted t(4). The t distribution is symmetric, and differs from the normal on account of kurtosis rather than skewness. The results in Table 1 for t(4) show that G α does not increase or decrease globally with α. However, as this example shows, the sensitivity to α provides information on the sort of divergence of the data distribution from normality. It is thus important to compare the finite-sample performance of G α with that of other standard goodness-of-fit tests.

Table 2 presents simulation results on the size and power of normality tests using Student's t and Gamma (Γ) distributions with several degrees of freedom, df = 2, 4, 6, . . . , 20. The t and Γ distributions provide two realistic examples that exhibit different types of departure from normality but tend to be closer to the normal as df increases. The values given in Table 2 are the percentages of rejections of the null H 0 : x ∼ Normal at 5% nominal level when the true distribution of x is F 0 , based on samples of 100 observations. Rejections are based on bootstrap P values for all tests, not just those that use G α . When F 0 is the standard normal distribution (first line), the results measure the Type I error of the tests, by giving the percentage of rejections of H 0 when it is true. For nominal level of 5%, we see that the Type I error is small. When F 0 is not the normal distribution (other lines of the Table ), the results show the power of the tests. The higher a value in the table, the better is the test at detecting departures from normality. As expected, results show that the power of all statistics considered increases as df decreases and the distribution is further from the normal distribution. Among the standard goodness-of-fit tests, Table 2 shows that the AD statistic is better at detecting most departures from the normal distribution (italic values). The CVM statistic is close, but KS and P have poorer power. Similar results are found in [START_REF] Stephens | Tests based on EDF statistics[END_REF]. Indeed, the Pearson chi-square test is usually not recommended as a goodness-of-fit test, on account of its inferior power properties.

Among the G α goodness-of-fit tests, Table 2 shows that the detection of greatest departure from the normal distribution is sensitive to the choice of α. We can see that, in most cases, the most powerful G α test performs better than the most powerful standard test (bold vs.italic values). In addition, it is clear that G α increases with α when the data are generated from the Gamma distribution. This is due to the fact that the Gamma distribution is skewed to the right.

Tests for other distributions

Table 3 presents simulation results on the power of tests for the lognormal distribution. 10 The values given in the table are the percentages of rejections of the null H 0 : x ∼ lognormal at level 5% when the true distribution of x is the Singh-Maddala distribution -see [START_REF] Singh | A Function for Size Distribution of Incomes[END_REF] -of which the CDF is

F SM (x) = 1 - 1 (1 + ax b ) c
with parameters a = 100, b = 2.8, and c = 1.7. We can see that the most powerful G α test (α = 1) performs better than the most powerful standard test (bold vs.italic values). The least powerful G α test (α = 5) performs similarly to the KS test. Table 4 presents simulation results on the power of tests for the Singh-Maddala distribution. The values given in the table are the percentage of rejections of the null H 0 : x ∼ SM at 5% when the true distribution of x is lognormal. We can see that the most powerful G α test (α = 5) performs better than the most powerful standard test (bold vs. italic values).

Note that the two experiments concern the divergence between Singh-Maddala and lognormal distributions, but in opposite directions. For this reason the G α tests are sensitive to α in opposite directions.

Application

Finally, as a practical example, we take the problem of modelling income distribution using the UK Households Below Average Incomes 2004Incomes -5 (2006)).

The application uses the "before housing costs" income concept, deflated and equivalised using the OECD equivalence scale, for the cohort of ages 21-45, couples with and without children, excluding households with selfemployed individuals. The variable used in the dataset is oe bhc. We exclude households with self-employed individuals as reported incomes are known to be misrepresented. The empirical distribution F consists of 3858 observations and has mean and standard deviation (398.28, 253.75).

We test the goodness-of-fit of a Singh-Maddala distribution, with 999 bootstrap samples used to compute the bootstrap P values. Table 5 presents the results with standard goodness-of-fit tests, Table 6 presents the results with G α tests. If we use standard goodness-of-fit statistics, we would not reject the Singh-Maddala distribution in most cases. Conversely, if we use G α goodnessof-fit statistics, we would reject the Singh-Maddala distribution in most cases. Our previous simulation study shows G α and AD have better finite sample properties. This leads us to conclude that the Singh-Maddala distribution is not a good fit, contrary to the conclusion from standard goodness-of-fit tests only.

Concluding Remarks

The family of goodness-of-fit tests presented in this paper has been seen to have excellent size and power properties as compared with other, commonly used, goodness-of-fit tests. It has the further advantage that the profile of the G α statistic as a function of α can provide valuable information about the nature of the departure from the target family of distributions, when that family is wrongly specified.

We have advocated the use of the parametric bootstrap for tests based on G α . The distributions of the limiting random variables (39) and (34) exist, as shown, but cannot be conveniently used without a simulation experiment that is at least as complicated as that involved in a bootstrapping procedure. In addition, there is no reason to suppose that the asymptotic distributions are as good an approximation to the finite-sample distribution under the null as the bootstrap distribution. We rely on the mere existence of the limiting distribution in order to justify use of the bootstrap. The same reasoning applies, of course, to the conventional goodness-of-fit tests studied in Section 4. They too give more reliable inference in conjunction with the parametric bootstrap.

Of course, the G α statistics for different values of α are correlated, and so it is not immediately obvious how to conduct a simple, powerful, test that works in all cases. It is clearly interesting to compute G α for various values of α, and so a solution to the problem would be to use as test statistic the maximum value of G α over some appropriate range of α. The simulation results in the previous section indicate that a range of α from -2 to 5 should be enough to provide ample power. It would probably be inadvisable to consider values of α outside this range, given that it is for α = 2 that the finite-sample distribution is best approximated by the limiting asymptotic distribution. That said, the maximum of G α has the same limiting distribution as G α for any fixed α, because there is only one limit distribution for all α. Thus bootstrapping the maximum over the chosen range is asymptotically justified, and seems likely to give a powerful test, to the extent that our simulation results can be generalised.

Figure 1 :

 1 Figure 1: Different types of divergence of the data distribution from the model

Table 1 :

 1 Normality tests with G α based on 1000 observations drawn from Beta and t distributions These observations are confirmed by the results of Table

	α	-2	-1	0	0.5	1	2	5	10
	B(5,2) 2.29 2.03 1.85 1.79 1.73 1.64	1.47	1.35
	B(2,5) 3.70 4.02 4.6 5.15 6.01 11.09 1.37e4 3.34e11
	t(4)	61.35 6.83 4.17 3.99 3.94 4.02	4.74	7.30

Table 2 :

 2 Normality tests: percentage of rejections of H 0 : x ∼ Normal, when the true distribution of x is F 0 . Sample size = 100, 5000 replications, 999 bootstraps.

			Standard tests				G α test with α =
	F 0	AD CVM KS	P	-2	-1	0	0.5	1	2	5
	N(0,1) 5.3	5.2	5.4	5.4	4.6	4.6	4.7	5.0	5.1	5.2	5.4
	t(20)	7.7	7.3	6.6	5.8 11.7 10.4 7.3	6.6	6.7	6.5	6.2
	t(18)	8.9	8.3	6.6	5.5 12.4 11.5 8.0	7.4	7.4	7.5	6.9
	t(16)	9.9	8.9	7.1	6.3 13.5 12.9 9.4	8.6	8.6	8.7	8.0
	t(14)	9.8	8.8	7.5	6.0 15.0 13.8 9.4	8.7	8.5	9.0	8.2
	t(12) 13.5	12.0	8.9	6.5 17.8 17.8 12.7 11.8 11.7 11.9 11.0
	t(10) 15.2	12.8 10.3 6.7 21.8 21.3 15.2 13.5 13.4 13.6 12.4
	t(8)	22.3	19.0 13.4 8.2 26.5 26.5 20.7 19.1 19.0 19.4 17.7
	t(6)	37.5	33.0 24.1 13.6 34.4 37.3 33.4 32.2 31.9 32.7 29.8
	t(4)	64.3	59.9 48.5 28.6 49.6 59.9 59.4 58.5 58.7 59.5 56.6
	t(2)	98.0	97.6 95.2 87.6 87.3 96.4 97.0 97.1 97.2 97.3 96.9
	Γ(20) 25.2	21.9 17.8 10.2 0.1	4.5 13.8 16.1 18.4 23.2 36.3
	Γ(18) 28.3	25.1 20.9 10.7 0.1	5.8 16.4 19.3 22.0 27.2 40.0
	Γ(16) 30.9	27.2 21.9 12.0 0.1	7.1 18.9 22.0 24.5 29.5 42.6
	Γ(14) 34.5	30.3 24.4 11.8 0.1	8.7 21.2 25.1 28.1 34.5 49.3
	Γ(12) 41.3	36.6 28.5 14.5 0.1	10.7 26.4 30.3 34.0 40.6 56.2
	Γ(10) 48.9	42.4 34.0 17.1 0.1	14.2 32.3 36.5 41.1 48.5 64.4
	Γ(8)	58.1	51.7 41.6 22.0 0.1	19.9 41.7 47.1 51.6 59.7 74.8
	Γ(6)	72.7	65.4 52.3 31.0 0.5	31.4 57.5 63.0 67.7 75.5 87.8
	Γ(4)	88.5	82.1 68.8 49.7 2.0	55.7 79.6 84.0 87.0 92.1 97.5
	Γ(2)	99.8	99.3 95.4 95.3 22.5 96.5 99.4 99.7 99.8 99.9 100

Table 3 :

 3 Lognormality tests: percentage of rejections of H

0 : x ∼ lognormal, when the true distribution of x is Singh-Maddala(100,2.8,1.7). 5000 replications, 499 bootstraps.

Table 4 :

 4 Singh-Maddala tests: percentage of rejections of H 0 : x ∼ SM, when the true distribution of x is lognormal(0,1). 1000 replications, 199 bootstraps.

		Standard tests				G α test with α =
	nobs AD CVM KS	P	-2	-1	0	0.5	1	2	5
	500 53.6	43.3 32.3 16.7 11.3 37.3 47.7 50.2 53.0 57.4 73.5
	600 65.8	52.6 37.4 20.1 18.6 51.3 60.1 62.4 64.5 68.4 83.3
	700 75.7	61.8 43.7 22.8 24.9 61.4 71.5 73.3 74.4 77.9 87.4
	800 82.3	69.3 53.1 27.6 37.9 72.5 79.3 80.6 82.6 85.8 93.6
	900 87.7	75.9 54.8 30.6 45.8 77.5 82.9 83.9 85.6 88.5 93.7
	1000 91.2	80.9 62.8 34.2 55.7 82.6 86.9 88.1 89.4 92.4 96.4

Table 5 :

 5 Standard goodness-of-fit tests: bootstrap P values, H 0 : x ∼ SM.

	α	-2	-1	0	0.5	1	2	5
	statistic 164	1.36 0.441 0.404 0.390 0.382 0.398
	p-value 0.003 0.002 0.003 0.006 0.007 0.010 0.013

Table 6 :

 6 G α goodness-of-fit tests: bootstrap P values, H 0 : x ∼ SM.

As examples of the axiomatic method seeSen (1976a) on national income,[START_REF] Sen | Poverty: An ordinal approach to measurement[END_REF] on poverty and[START_REF] Ebert | Measurement of inequality: an attempt at unification and generalization[END_REF] on inequality.

This implies that it has the minimal properties of completeness, reflexivity and transitivity.

For any z, z ∈ Z n "z z " means "[z z ] & [z z]"; "z ∼ z " means "[z z ] & [z z]".

Also note that Axiom 6 can be stated equivalently by requiring that, for a given z 0 , z 0 ∈ Z n such that z 0 ∼ z 0 , either (a) any z and z found by rescaling the u-components will be equivalent or (b) any z and z found by rescaling the v-components will be equivalent.

 5 The following proof draws on[START_REF] Ebert | Measurement of inequality: an attempt at unification and generalization[END_REF].

There is a natural extension to the case where one is dealing with a continuous distribution on support Y ⊆ R. Expressions (25) -(27) become, respectively:

In fact, the true variance of z i is t i (1 -t i )/(n + 2).

We use the standard tests as implemented with R, the number of intervals m is due to[START_REF] Moore | Tests of the chi-squared type[END_REF]. Note that G, AD, CVM and KS statistics are based on the empirical distribution function (EDF) and the P statistic is based on the density function.

Results under the null are close to the nominal level of 5%. For n = 50, we obtain rejection rates, for AD, CVM, KS, Pearson and G with α = -2, -1, 0, 0.5, 1, 2, 5 respectively, of 5.02, 4.78, 4.76, 4.86, 5.3, 5.06, 4.88, 4.6, 4.72, 5.18.


Jérémie Bigot
email: jeremie.bigot@math.univ-toulouse.fr

Fréchet means of curves for signal averaging and application to ECG data analysis

Keywords: Signal averaging, Mean shape, Fréchet means, Curve registration, Geometric variability, Deformable models, ECG data. AMS classifications: Primary 62G08, secondary 62P10

Signal averaging is the process that consists in computing a mean shape from a set of noisy signals. In the presence of geometric variability in time in the data, the usual Euclidean mean of the raw data yields a mean pattern that does not reflect the typical shape of the observed signals.

In this setting, it is necessary to use alignment techniques for a precise synchronization of the signals, and then to average the aligned data to obtain a consistent mean shape. In this paper, we study the numerical performances of Fréchet means of curves which are extensions of the usual Euclidean mean to spaces endowed with non-Euclidean metrics. This yields a new algorithm for signal averaging without a reference template. We apply this approach to the estimation of a mean heart cycle from ECG records.

Introduction

In many applications (biology, medicine, road traffic data) one observes a set of J signals that have a similar shape. This may lead to the assumption that such observations are random elements which vary around the same but unknown mean shape. Signal averaging is then the process that consists in computing a mean curve which reflects the typical shape of the observed signals. This procedure generally amounts to find an appropriate combination of the data to compute an average shape with a better signal-to-noise ratio. In many situations, the observed signals exhibit not only a classical source of random variation in amplitude, but also a less standard geometric source of variability in time. Due to this source of geometric variability, the usual Euclidean mean of the raw data may yield a mean curve that does not reflect the typical shape of the signals, as illustrated by the following application.

Signal averaging in ECG data analysis

An important application of signal averaging is the estimation of a mean heart cycle from electrocardiogram (ECG) records. An ECG signal corresponds to the recording of the heart electrical activity. It is a signal, recorded over time, that is composed of the succession of cycles of contraction and release of the heart muscle. Each recorded cycle is a curve composed of a characteristic P-wave, reflecting the atrial depolarization, that is followed by the so-called QRS complex which corresponds to the depolarization of the ventricles, and which ends with a T-wave reflecting the repolarization of the heart (see e.g. [START_REF] Guyton | Textbook of medical physiology[END_REF] for a precise description of an ECG recording). In this paper, we present results on two different data sets from the MIT-BIH database [GAG + 13]. In Figure 1(a), we display data from an ECG record of a healthy subject having no arrhythmia, and in Figure 1(b) we display data from an ECG record of a subject showing evidence of significant arrhythmia (note that, in all the figures showing ECG data, units on the vertical axis are in millivolts).

In the analysis of ECG data, it is generally assumed that the heart electrical activity repeats itself. Therefore, during an ECG record, one classically considers that the heart cycle of interest remains approximately the same with every beat, and that it is embedded in a random white noise with zero expectation that is uncorrelated with the mean shape to be estimated. After an appropriate segmentation of an ECG record, one observes a set of signals of the same length such that each of them contains a single QRS complex. The preliminary segmentation step is done by taking segments (of the same length) in the ECG record that are centered around the easily detectable maxima of the QRS complex of the beats. Identification of these maxima can be done using statistical methods to identity local extrema in noisy signals [START_REF] Bigot | Landmark-based registration of curves via the continuous wavelet transform[END_REF][START_REF] Gasser | Searching for Structure in Curve Sample[END_REF] or by applying appropriate digital filters to identify typical parts of the QRS complex [START_REF] Pan | A real-time QRS detection algorithm[END_REF]. In this paper, we used the approach in [START_REF] Bigot | Landmark-based registration of curves via the continuous wavelet transform[END_REF] to identify local maxima and to segment the ECG record into signals of the same length containing a single QRS complex. After this preliminary segmentation, one thus observes signals with approximately the same shape. In the normal case, we obtained J = 93 signals displayed in Figure 2, while in the case of cardiac arrhythmia we obtained J = 72 signals displayed in Figure 3.

To estimate the typical shape of a heart cycle and to improve the signal-to-noise ratio, one may use the Euclidean mean of these signals. In the normal case, this leads to satisfactory results as shown in Figure 4 since the average signal displayed in Figure 4(a) clearly reflects the typical shape of the heart beats displayed in Figure 4(b)-(g). However, in the case of cardiac arrhythmia, the electrical activity of the heart is more irregular. This can be seen in the shape of the heart beats displayed in Figure 3(b)-(g) which may vary significantly from one pulse to another. Due to noise and geometric variability in the measurements, a simple averaging step may cause a low-pass filtering effect that leads to a mean cycle that does not reflect the typical shape of the heart beats in the ECG record [START_REF] Rompelman | Coherent averaging technique: a tutorial review. part 1: Noise reduction and the equivalent filter[END_REF][START_REF] Rompelman | Coherent averaging technique: a tutorial review. part 2: Trigger jitter, overlapping responses and non-periodic stimulation[END_REF][START_REF] Laciar | Improved alignment method for noisy high-resolution ecg and holter records using multiscale cross-correlation[END_REF]. In the case of cardiac arrhythmia, the Euclidean mean causes a low-pass filtering effect in the shape of the QRS complex, as shown by the results displayed in Figure 5. Indeed, it can be seen in Figure 3 that, around the time point t ≈ 0.45 (which corresponds to the beginning of the QRS complex) there is rapid transition between a flat region and the peak of the R wave. For the Euclidean mean displayed in Figure 5(a), this transition is slower, and it can be seen that the average heart cycle in Figure 5(a) is a convolution by a smooth kernel of the signals in Figure 3(b)-(g) around the time point t ≈ 0.45. To obtain better results, it is necessary to use alignment techniques for a precise synchronization of the heart beats, and then to average the aligned data.

Deformable models

Assume for simplicity that the signals are observed on the time interval [0, 1], and that they can be extended outside [0, 1] by periodicity. An alignment technique consists in finding a time synchronization of a set of signals. To be more precise, define a deformation operator φ θ parametrized by θ ∈ R p as a smooth increasing function

φ θ : [0, 1] → R such that φ -1 θ (t) = φ -θ (t) for all t ∈ [0, 1].
In the paper, we shall consider the following families of deformation operators:

-Translation operators:

φ θ (t) = t -θ and φ -1 θ (t) = φ -θ (t) = t + θ, for θ ∈ R (p = 1) and all t ∈ [0, 1].
-Non-rigid operators:

φ θ : [0, 1] → [0, 1] is diffeomorphism of [0, 1] parametrized by some θ ∈ R p , i.e.
a smooth increasing function with φ θ (0) = 0 and φ θ (1) = 1 (a general method for constructing non-rigid deformation operators is described in Section 2.4). Given

f 1 , f 2 : [0, 1] → R and a family (φ θ ) θ∈R p of deformation operators, the problem of time synchronization of two signals is to find a θ ∈ R p such f 1 (φ θ (t)) ≈ f 2 (t) for all t ∈ [0, 1].
In ECG data analysis, the most widely used alignment technique is time synchronization using translation operators by temporal or multiscale cross-correlation, see [START_REF] Laciar | Improved alignment method for noisy high-resolution ecg and holter records using multiscale cross-correlation[END_REF][START_REF] Trigano | Semiparametric curve alignment and shift density estimation for biological data[END_REF] and references therein.

In the presence of variability in amplitude and shape in the data, it is reasonable to suppose that the signals at hand satisfy the following deformable (or perturbation) model:

Y ℓ j = f (φ θ * j (t ℓ )) + Z j (φ θ * j (t ℓ
)) + σε ℓ j , j = 1, . . . , J, and ℓ = 1, . . . , n,

where Y ℓ j denotes the ℓ-th observation for the j-th signal and t ℓ = ℓ n are equi-spaced design points in [0, 1]. The function f : [0, 1] → R in model (1.1) is the unknown mean shape of the signals. The ε ℓ j are supposed to be i.i.d. normal variables with zero expectation and variance 1. They represent a standard white noise in the measurements. The Z j 's are supposed to be independent of the ε ℓ j 's and are i.i.d. realizations of a smooth random process Z : [0, 1] → R with zero expectation. Therefore, the Z j 's represent smooth variations in amplitude of the signals around f . Finally, the θ * j 's are assumed to be i.i.d. random variables in R p with zero expectation. Hence, the random deformation operators φ θ * j represent geometric variability in time in the data.

Note that the deformable model (1.1) is well adapted to ECG data processing. The main type of perturbations related to the analysis of ECG data (see e.g. [START_REF] Laciar | Improved alignment method for noisy high-resolution ecg and holter records using multiscale cross-correlation[END_REF][START_REF] Trigano | Semiparametric curve alignment and shift density estimation for biological data[END_REF]) are the baseline wandering effect (a low-frequency signal), electromyographic (EMG) noise, and power-line interference which is an amplitude and frequency varying sinusoid. This source of additive noise is modeled in (1.1) by the zero-mean processes Z j . The physiological nature of the electrocardiographic signal also alters the recording from heartbeat to heartbeat in lag and duration. In the ECG signal, there are therefore variations in time of the heart cycle from one beat to another. This makes the heartbeats look shorter or longer in duration. After the segmentation of an ECG record into signals containing a single QRS complex, this local variability in time is modeled in (1.1) by the non-rigid deformation operators φ θ * j . Aligning heart beats using non-rigid deformation operators is an alternative to the cross-correlation method which is classical used in ECG data analysis. One of the goals of this paper is to show that using non-rigid operators to align signals obtained by the segmentation of an ECG record may significantly improve the shape of the mean heart cycle computed by averaging the data after a standard alignment using cross-validation.

Fréchet means of curves

The problem of estimating f and the deformation parameters θ * j in model (1.1) has been studied in [START_REF] Bigot | On the consistency of fréchet means in deformable models for curve and image analysis[END_REF] using the following procedure. First, for each j = 1, . . . , J, smooth the data (Y ℓ j ) n ℓ=1 to construct an estimator fj :

[0, 1] → R of f • φ θ * j .
In this paper, this smoothing step is done either by low-pass Fourier filtering or by wavelet thresholding. In a second step, estimate simultaneously the deformation parameters θ * j , j = 1, . . . , J by minimizing the following criterion ( θ1 , . . . , θJ ) = argmin

(θ 1 ,...,θ J )∈Θ 0 M (θ 1 , . . . , θ J ), (1.2)
where

M (θ 1 , . . . , θ J ) = 1 J J j=1 1 0 fj (φ -θ j (t)) - 1 J J j ′ =1 fj ′ (φ -θ j ′ (t)) 2 dt, (1.3)
and

Θ 0 = {(θ 1 , . . . , θ J ) ∈ (R p ) J , θ 1 + . . . + θ J = 0}.
Finally, in a third step take

f (t) = 1 J J j=1 fj (φ -θj (t)), (1.4)
as an estimator of the mean shape f . As explained in [START_REF] Bigot | On the consistency of fréchet means in deformable models for curve and image analysis[END_REF], the estimator f can be interpreted as a smoothed Fréchet mean of the observed signals. The Fréchet mean [START_REF] Fréchet | Les éléments aléatoires de nature quelconque dans un espace distancié[END_REF] is an extension of the usual Euclidean mean to spaces endowed with non-Euclidean metrics. We refer to [START_REF] Afsari | Riemannian l p center of mass: existence, uniqueness, and convexity[END_REF] and [START_REF] Huckemann | Intrinsic inference on the mean geodesic of planar shapes and tree discrimination by leaf growth[END_REF] for recent overviews of this notion and its application to the analysis of random variables taking their values in non-linear manifolds. Note that this procedure does not require the use of a reference template to compute the estimators θ1 , . . . , θJ of the deformation parameters. Instead, one can interpret 1

J J j ′ =1 fj ′ (φ -θj ′ (t))
as a template that is automatically computed during the minimization of (1.3), and onto which one searches to align all the curves fj . Therefore, this approach avoids the problem of taking one of the observed signals as a reference template which can lead to poor alignment results for data with a low signal-to-noise ratio. It should be also remarked that the constrained set Θ 0 (onto which the minimization of M (θ 1 , . . . , θ J ) is done) reflects the assumption that the deformation parameters θ * j in model (1.1) have zero expectation. The choice of this constraint is also related to identifiability issues in model (1.1), and we refer to [START_REF] Bigot | On the consistency of fréchet means in deformable models for curve and image analysis[END_REF] for a detailed discussion on that point.

In [START_REF] Bigot | On the consistency of fréchet means in deformable models for curve and image analysis[END_REF], the statistical properties of f and the θj 's have been studied in details in the asymptotic setting n → +∞ and/or J → +∞. However, the numerical performances of Fréchet means of curves have not been much studied in [START_REF] Bigot | On the consistency of fréchet means in deformable models for curve and image analysis[END_REF]. The goal of this paper is then two folds. First, we describe new algorithms to compute f and the θj 's in the case of translation and non-rigid operators. We also discuss the choice of the smoothing estimators fj . Secondly, we show the benefits of the methodology proposed in [START_REF] Bigot | On the consistency of fréchet means in deformable models for curve and image analysis[END_REF] for signal averaging of ECG data.

Previous work on signal averaging

In statistics, the problem of estimating the mean shape of a set of curves that differ by a time transformation is usually referred to as the curve registration problem. It has received a lot of attention in the statistical literature over the last two decades. Among the various methods that have been proposed, one can distinguish between landmark-based approaches see e.g. [START_REF] Gasser | Statistical tools to analyze data representing a sample of curves[END_REF], [START_REF] Bigot | Landmark-based registration of curves via the continuous wavelet transform[END_REF], and time synchronization to align a set of curves see e.g. [START_REF] Ramsay | Curve registration[END_REF], [START_REF] Wang | Alignment of curves by dynamic time warping[END_REF], [START_REF] Liu | Functional convex averaging and synchronization for time-warped random curves[END_REF], [START_REF] Kneip | Convergence and consistency results for self-modelling regression[END_REF], [START_REF] Trigano | Semiparametric curve alignment and shift density estimation for biological data[END_REF].

All these approaches are based either on the use of one the observed signals as a reference template (see e.g. [START_REF] Trigano | Semiparametric curve alignment and shift density estimation for biological data[END_REF] for curve alignment with translation operators), or on the Procrustean mean which is a standard algorithm to estimate a mean shape without a reference template. The Procrustean mean is based on an alternative scheme between estimation of deformation operators and averaging of back-transformed curves given estimated values of the deformation parameters. To be more precise, assume that Y j : [0, 1] → R denotes a linear interpolation of the data (Y ℓ j ) n ℓ=1 . Using our notations, Procrustean mean consists of an initialization step f0 = 1 J J j=1 Y j which is the Euclidean mean of the raw data taken as a first reference template. Then, at iteration 1 ≤ i ≤ i max , it computes for all 1 ≤ j ≤ J the estimators θj,i of the j-th deformation parameter as θj,i = arg min

θ∈R p 1 0 |Y j (φ θ j (t)) -fi-1 (t)| 2 dt
and then takes fi (t) = 1 J J j=1 Y j (φ θj,i (t))) as a new reference template. This procedure is repeated until the estimated reference template does not change. Usually the algorithm converges in a few steps.

However, in the case of a low signal-to-noise ratio, it is preferable to first smooth the data before alignment. In this paper, to minimize the criterion (1.3), we propose an alternative algorithm to Procrustean mean. We also show that the obtained smoothed Fréchet mean (1.4) yields better results. In the statistical literature, the criterion (1.2) has been first proposed by [START_REF] Gamboa | Semi-parametric estimation of shifts[END_REF] in the case of translation operators and then further studied by [START_REF] Vimond | Efficient estimation for a subclass of shape invariant models[END_REF] and [START_REF] Bigot | A deconvolution approach to estimation of a common shape in a shifted curves model[END_REF]. Its generalization to other deformation operators and the connection between minimizing (1.3) and the computation of Fréchet means of curves has been investigated in [START_REF] Bigot | On the consistency of fréchet means in deformable models for curve and image analysis[END_REF].

Organization of the paper

In Section 2 we describe a gradient descent algorithm to compute the θj 's. We detail this procedure in the case of translation and non-rigid operators. We also compare our approach with the Procrustean mean of curves. Numerical experiments on signal averaging for ECG data are presented in Section 3. We conclude the paper by a brief discussion on the numerical results and the benefits of our approach.

Methodology for mean pattern estimation 2.1 A gradient descent algorithm with an adaptive step

Denote by

∇M (θ) = ∂ ∂θ 1 M (θ), . . . , ∂ ∂θ J M (θ) ∈ (R p ) J ,
the gradient of the criterion (1.3) at θ = (θ 1 , . . . , θ J ) ∈ (R p ) J . A closed form expression of this gradient is given in Section 2.3 for translation operators, and in Section 2.4 for non-rigid operators.

Computation of the deformation parameters can then be done by minimizing the criterion M (θ) using the following gradient descent algorithm with the constraint that J j=1 θ j = 0:

Initialization : let θ 0 = 0, δ 0 = 1 ∇M (θ 0 ) , C(0) = M (θ 0
), and set i = 0.

Step 2 : let θ = θ i -δ i ∇M (θ i ) and let

θ new = θ1 -1 J J j=1 θj , . . . , θJ -1 J J j=1 θj . Let C(i + 1) = M (θ new ). While C(i + 1) > C(i) do δ i = δ i /κ, θ = θ i -δ i ∇M (θ i ), θ new =  θ 1 - 1 J J j=1 θj , . . . , θJ - 1 J J j=1 θj   ,
and set C(i + 1) = M (θ new ).

End while

Take θ i+1 = θ new .

Step 3 : if C(i) -C(i + 1) ≥ ρ(C(1) -C(i + 1)) then set i = i + 1 and return to Step 2, else stop the iterations, and take θ = θ i+1 .

In the above algorithm, ρ > 0 is a stopping parameter and κ > 1 is a parameter to control the choice of the adaptive step δ i .

Choice of the regularization parameter in the smoothing step

For the smoothing step, we present numerical results for low-pass Fourier filtering:

for t ∈ [0, 1] fj (t) = |k|≤ λj c (j) k e i2πkt , with c (j) k = 1 n n ℓ=1 Y j,ℓ e -i2πk ℓ
n , and where λj ∈ N is a regularization parameter (cut-off frequency). A possible data-based choice for λj is to use generalized cross validation (GCV), see e.g. [START_REF] Craven | Smoothing noisy data with spline functions: Estimating the correct degree of smoothing by the method of generalized cross-validation[END_REF].

wavelet smoothing by hard thresholding:

for t ∈ [0, 1] fj (t) = 2 m 0 k=0 α (j) m 0 ,k φ m 0 ,k (t) + m 1 m=m 0 2 m k=0 β (j) m,k 1 β (j) m,k ≥σ j √ 2 log(n) ψ m,k (t),
where

φ m 0 ,k (t) = 2 m 0 2 φ(2 m 0 t -k) and ψ m,k (t) = 2 m 2 ψ(2 m t -k)
are the usual scaling and wavelet basis functions at resolution levels 0 ≤ m 0 ≤ m ≤ m 1 and location k, α

(j) m 0 ,k , β (j)
m,k are respectively the empirical scaling and wavelet coefficients computed from the data (Y ℓ j ) n ℓ=1 (see e.g. [ABS01] for further details on wavelet thresholding). The universal threshold σj 2 log(n) depends on the estimation σj of the level of noise in the j-th signal. It is given by the median absolute deviation (MAD) of the empirical wavelet coefficients at the highest level of resolution m 1 (see e.g. [ABS01]).

The case of translation operators

When φ θ (t) = t -θ, it follows by simple computations that for j = 1, . . . , J

∂ ∂θ j M (θ) = - 2 J 1 0 ∂ ∂t fj (t + θ j )   1 J J j ′ =1 fj ′ (t + θ j ′ )   dt + 2 J | fj (1 + θ j )| 2 -| fj (θ j )| 2 .
Note that in the case of low-pass Fourier filtering, fj (t + 1) = fj (t) for any t ∈ R, and the gradient of M can be computed directly in the Fourier domain thanks to Parseval's relation

∂ ∂θ j M (θ) = - 2 J |k|≤λ ℜ   2iπkc (j) k e 2ikπθ j   1 J J j ′ =1 c (j ′ ) k e 2ikπθ j ′     ,
where ℜ(c) denotes the real part of a complex number c.

The case of non-rigid operators

To build a family (φ θ ) θ∈R p of parametric diffeomorphisms of [0, 1], we adapt to one-dimensional curves the approach proposed in [START_REF] Bigot | Statistical M-estimation and consistency in large deformable models for image warping[END_REF] to compute the mean pattern of a set of two-dimensional images. Let v : [0, 1] → R be a smooth parametric vector field given by a linear combination of p basis functions

{h k : [0, 1] → R, k = 1, . . . , p}, such that v(t) = p k=1 θ k h k (t) for t ∈ [0, 1],
where θ = (θ 1 , . . . , θ p ) ∈ R p is a set of reals coefficients. The function v is thus parametrized by the set of coefficients θ, and we write v = v θ to stress this dependency. In what follows, it will be assumed that the basis functions are continuously differentiable on [0, 1] and such that h k and ∂ ∂t h k vanish at t = 0 and t = 1. For the h k 's we took in our numerical experiments a set of p = 10 B-spline functions of degree 3 using equally-spaced knots on [0, 1]. Then, let t ∈ [0, 1] and for u ∈ [0, 1] consider the following ordinary differential equation (ODE)

∂ ∂u ψ(u, t) = v θ (ψ(u, t)) (2.1)
with initial condition ψ(0, t) = t. Then, it can be shown (see e.g. [START_REF] Younes | Shapes and diffeomorphisms[END_REF]) that for any u ∈ [0, 1] the solution of the above ODE is unique and such that

t → ψ θ (u, t) = t + u 0 v θ (ψ θ (s, t))ds
is a diffeomorphism of [0, 1] i.e. a smooth increasing function with ψ θ (u, 0) = 0 and ψ θ (u, 1) = 0. Then, we denote by

φ θ (t) = ψ θ (1, t)
the solution at u = 1 of the ODE (2.1). In this way, we finally obtain a diffeomorphism φ θ that is parametrized by the set of coefficients θ ∈ R p , and that is such that φ -1 θ (t) = φ -θ (t). Then, it can be also shown that the mapping θ → φ θ is differentiable, and thanks to Lemma 2.1 in [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF] it follows that

∂ ∂θ k φ θ (t) = ∂ ∂t φ θ (t) 1 0 ∂ ∂t ψ θ (u, t) -1 h k (ψ θ (u, t))du.
Now, one obtains that if φ θ is a diffeomorphism generated by the ODE (2.1), then for θ j = (θ

(1) j , . . . , θ (p) j ) ∈ R p ∂ ∂θ j M (θ) = ∂ ∂θ (1) j M (θ), . . . , ∂ ∂θ (p) j M (θ) , with for q = 1, . . . , p ∂ ∂θ (q) j M (θ) = - 2 J 1 0 ∂ ∂θ (q) j φ -θ j (t) ∂ ∂t fj (φ -θ j (t))   1 J J j ′ =1 fj ′ (φ -θ j ′ (t)) -fj (φ -θ j (t))   dt and ∂ ∂θ (q) j φ -θ j (t) = ∂ ∂t φ -θ j (t) 1 0 ∂ ∂t ψ -θ j (u, t) -1
h q (ψ -θ j (u, t))du.

Comparison between the smoothed Fréchet mean and the Procrustes mean

To illustrate the advantages of the smoothed Fréchet mean over the Procrustes mean of the raw data (without any smoothing), let us consider a set of J = 15 signals generated from the following deformable model using translation operators:

Y ℓ j = f (t ℓ -θ * j ) + Z j (t ℓ -θ * j ) + σε ℓ j , j = 1, . . . , J, and ℓ = 1, . . . , n, (2.2)
where n = 128, the θ * j 's are i.i.d. normal variables with zero mean and variance µ 2 = 0.004, the ε ℓ j 's are i.i.d. normal variables with zero mean and variance σ 2 , f is the signal displayed in Figure 6(a), and the Z j 's are i.i.d. copies of the Gaussian process Z defined by

Z(t) = a 0 + +∞ k=1 k -3/2 (a k c k (t) + b k s k (t)), t ∈ [0, 1],
where a 0 and (a k , b k ) k≥1 are i.i.d. normal variables with zero mean and variance σ 2 , c k (t) = √ 2 cos(2πkt) and s k (t) = √ 2 sin(2πkt). A sample of three signals out of J = 15 generated from model (2.2) is displayed in Figure 6(b)-(d). It can be seen that σ has been chosen to have a relatively low signal-to-noise ratio. The estimation of f by smoothed Fréchet mean using low-pass Fourier filtering with λj chosen by GCV for each j = 1, . . . , J is displayed in Figure 6(e), while the estimation obtained by Procrustes mean of the data without smoothing is shown in Figure 6(f). Clearly, the result obtained by the smoothed Fréchet mean is much more satisfactory. An important choice is thus the calibration of the regularization parameters λj . To illustrate the benefits of using the smoothed Fréchet mean with an automatic choice of the λj 's by GCV, we used M = 100 replications of model (2.2) with J = 15. For each repetition m = 1, . . . , M , we have computed an estimator f (m) P rocrustes of f using the Procrustes mean without smoothing, and an estimator f (m) F rechet of f using a smoothed Fréchet mean where each λj is chosen by denoising the j-th signal by GCV. In Figure 7, we display boxplots (for m = 1, . . . , M ) of the mean squared errors

M SE( f (m) F rechet ) := 1 n n ℓ=1 f (m) F rechet (t ℓ ) -f (t ℓ ) 2 and M SE( f (m) P rocrustes ) := 1 n n ℓ=1 f (m) P rocrustes (t ℓ ) -f (t ℓ ) 2 .
It can be seen that, over M = 100 replications, the mean squared error of the smoothed Fréchet mean is clearly much lower than the one obtained by Procrustes mean.

Application to ECG data analysis

Now, let us return to the analysis of the ECG records displayed in Figure 1. We propose to compare the results obtained by the smoothed Fréchet mean using either translation operators or non-rigid operators given by the family of parametric diffeomorphisms of [0, 1] described in Section 2.4. The smoothing of the J signals obtained after segmentation of the ECG records is done by wavelet thresholding with a data-based choice of the regularization parameters σj as explained in Section 2.2.

The computation of smoothed Fréchet means using translation operators in the normal and cardiac arrhythmia cases are displayed in Figure 8 and Figure 9. The results obtained using the Procrustes mean are very similar, and they are thus not reported. Note that using Fréchet or Procrustes means with translation operators is very similar to the alignment of signals using cross-correlation which is the widely used technique for signal averaging in ECG data analysis.

One can see that, for the normal case, the results are very satisfactory as the mean heart cycle displayed in Figure 8(a) is a good estimation of the typical shape of the signals displayed in Figure 8(b)-(g). For the case of cardiac arrhythmia, the results are not very satisfactory since there is still a low-pass filtering effect on the shape of the QRS complex (around the time point t ≈ 0.45) in the mean heart cycle displayed in Figure 9(a). In this case, the Fréchet mean does not represent very well the shape of the signals in Figure 9(b)-(g). This low-pass filtering effect is due to a much larger variability in shape of the signals (after segmentation) in the case of cardiac arrhythmia than in the normal case. Since the activity of the heart can be very irregular in cases of cardiac arrhythmia, modeling shape variability using only translation operators is not flexible enough. A more precise alignment to take into account a local variability in lag and duration of the heart beats is thus needed.

The computation of smoothed Fréchet means using non-rigid operators in the normal and cardiac arrhythmia cases are displayed in Figure 10 and Figure 11. In the normal case, the result is very similar to the one already obtained by the Fréchet mean using translation operators, as shown by the comparison of Figure 8 and Figure 10. In the case of cardiac arrhythmia, the results displayed in Figure 9 and Figure 11 clearly show the improvements obtained by using non-rigid deformation operators. The shape of the mean heart cycle in Figure 11(a) is a better estimation of the typical shape of the signals displayed in Figure 11(b)-(g). In particular, there is no low-pass filtering effect around the time point t ≈ 0.45 in the shape of the Fréchet mean.

Discussion and conclusion

We have presented a new algorithm for aligning heart beats extracted from an ECG record. Our approach is based on the notion of smoothed Fréchet means of curves using either translation or non-rigid deformation operators. In the case of a normal ECG obtained from a healthy subject, our approach yields a satisfactory mean heart cycle. Note that using Fréchet mean with translation operators is very similar to the computation of a mean heart cycle by aligning the data using crosscorrelation which is the widely used approach in ECG data analysis for signal averaging. When using non-rigid operators to align hearts beats having a high shape variability, with peaks showing an important variability in lag and duration from one pulse to another, our approach yields significant improvements over signal averaging by cross-correlation. The benefits of our procedure has been demonstrated for an ECG recording of a subject showing evidence of significant arrhythmia. We hope that the methods presented in this paper will stimulate further investigation into the development of better alignment procedures that take into account local shape variability in heart beats extracted from ECG records.

Figure 1 :

 1 Figure 1: (a) ECG recording during 60 seconds of a healthy subject having no arrhythmia, (b) ECG recording during 60 seconds of a subject showing evidence of significant arrhythmia.

Figure 2 :

 2 Figure 2: Normal case: six signals containing a single QRS complex out of J = 93 extracted from the ECG recording displayed in Figure 1(a). Units on the horizontal axis are arbitrary.

  0

Figure 4 :

 4 Figure 4: Normal case: (a) Euclidean mean of the J = 93 signals after segmentation of the ECG record. (b)-(g) Superposition of six signals containing a single QRS complex (dashed curves) with the Euclidean mean (solid curve).

Figure 5 :

 5 Figure 5: Case of cardiac arrhythmia: (a) Euclidean mean of the J = 72 signals after segmentation of the ECG record. (b)-(g) Superposition of six signals containing a single QRS complex (dashed curves) with the Euclidean mean (solid curve).

Figure 6 :

 6 Figure 6: (a) Unknown mean shape f , (b)-(d) Three noisy signals out of J = 15 sampled from model (2.2), (e) Estimation of f by smoothed Fréchet mean (solid curve) using low-pass Fourier filtering with λj chosen by GCV, (f) Estimation of f by Procrustes mean without smoothing (solid curve). The dashed curve in (e) and (f) is the mean shape f .

Figure 7 :

 7 Figure 7: Boxplots of M SE( f (m) F rechet ) (left-hand boxplot) and M SE( f (m) P rocrustes ) (right-hand boxplot) for m = 1, . . . , M .

Figure 8 :

 8 Figure 8: Normal case: (a) Fréchet mean using translation operators of the J = 93 signals after segmentation of the ECG record. (b)-(g) Superposition of six signals containing a single QRS complex (dashed curves) with the Fréchet mean using translation operators (solid curve).

  0

Figure 9 :

 9 Figure 9: Case of cardiac arrhythmia: (a) Fréchet mean using translation operators of the J = 72 signals after segmentation of the ECG record. (b)-(g) Superposition of six signals containing a single QRS complex (dashed curves) with the Fréchet mean using translation operators (solid curve).

Figure 10 :

 10 Figure 10: Normal case: (a) Fréchet mean using non-rigid operators of the J = 93 signals after segmentation of the ECG record. (b)-(g) Superposition of six signals containing a single QRS complex (dashed curves) with the Fréchet mean using non-rigid operators (solid curve).

Figure 11 :

 11 Figure 11: Case of cardiac arrhythmia: (a) Fréchet mean using non-rigid operators operators of the J = 72 signals after segmentation of the ECG record. (b)-(g) Superposition of six signals containing a single QRS complex (dashed curves) with the Fréchet mean using non-rigid operators (solid curve).

The author acknowledges the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-JCJC-SIMI1 DEMOS.