Je tiens tout d'abord à exprimer ma gratitude envers mes directeurs de thèse. Olivier, pour ton amitié, et parce que tu as su me montrer que l'on pouvait chercher en s'amusant, même si on ne trouvait rien. A Pascal, pour m'avoir guidé ces trois années, pour ta rigueur et ton excellence sans qui cette thèse n'aurait pas été possible. Et enfin à François, qui a été mon bol d'air frais, et qui m'aura permis d'étudier des sujets auxquels je n'aurais jamais pensé. Tu es le phare de ces locaux de Palaiseau ! J'aurai beaucoup appris durant ces années, et je vous le dois en grande partie.

I thank Stéphanie Bidon and Maarten Weyn for accepting to review this thesis. Je remercie également Aurélien Bellet et Maurice Charbit d'avoir accepté de constituer mon jury. Je vous en suis très reconnaissant.

Je tiens aussi à exprimer mes remerciements à la société Sigfox, et en particulier à l'équipe Data. Renaud, Oliv', le Duc: votre goût pour la recherche, votre confiance et votre intérêt inébranlable pour l'avancée de mes travaux ont formé un terreau idéal pour le déroulement de cette CIFRE. J'espère que les algorithmes développés pourront un jour atteindre une précision de localisation au mètre. Je suis convaincu que votre équipe continuera à faire de grandes choses. Merci également à tous les autres collègues de Sigfox Paris, grâce à qui aller au bureau était un plaisir renouvelé chaque vendredi (ou les jours de foot et de petits déjeuner). Donc à Pourya, Antho, Aurélien, Antoine, Laura C., Laura D. et j'en oublie... j'espère sincèrement vous retrouver autour d'un verre, d'un pizza, d'une pâtisserie ou d'un ballon rond ! Certaines rencontres sont bonnes pour la santé. Parmi les camarades qui ont rendu ces années si plaisantes, Alex, Charlito, le Quent, Ricco, Vous avez été, tous à votre manière, une source impérissable de bon air (doit-on dire débonnaire).

Pierre, Mathus, voi due siete i miei due accoliti, non oso pensare come sarebbero stati quegli anni senza di voi. Grazie. À mes parents et à mes soeurs, qui m'auront toujours encouragé.

À Océane pour qui les mots m'ont toujours manqués, et à Zazou qui est venu illuminer ces derniers jours de doctorat.

Predicting the location of an emitter from the RSSI can be split in two categories of methods, which are referred to as likelihood-based ones and fingerprinting ones. The first ones, considered in this manuscript, first learn a likelihood model for the conditional probability of the RSSI vector given the position; then, the second stage consists to finding the position which best agrees with the observed RSSI vector, by means of a grid search for instance. In contrast, the second ones directly map the vector of RSSI into a location (radio-map), typically by means of a supervised learning algorithm. In this work, we first propose improvements of methods used for the RSSI-based geolocation problem. The first proposed technique relies on a semi-parametric Nadaraya-Watson estimator of the likelihood, followed by a maximum a posteriori estimator of the object's position. The second technique consists to learning a distance, constructed by means of a Gradient boosting regressor: a k-nearest neighbor algorithm is then used to estimate the position. The proposed methods are compared on two data sets originated from Sigfox network, and an indoor dataset from a three-story building. Experiments demonstrate the interest of the proposed methods, both in terms of location estimation performance, and ability to build radio maps. Results also show that the quality of the prediction is highly related to the chosen distance on the RSSI space. The metric learning problem is therefore a fundamental issue to improve RSSI-based geolocation technique.

Second, we introduce an original objective for learning a similarity between pairs of data points. In this manuscript, we propose to build the similarity by directly minimizing the regression error of an estimator. We thus obtain a task-driven learning objective. To minimize it, the similarity is chosen as a sum of regression trees and is sequentially learned by means of a modified version of XGBoost detailed in this document. This method benefits from the well-known qualities of XGBoost such as its efficiency and its scaling capabilities. Furthermore, our similarity, although non-parametric, does not require a storage of the size of the dataset. Finally, experiments show that our model outperforms other kernel regression models on several benchmark datasets.

Conditional Independence has been broadly used in the RSSI-based geolocation literature in order to decrease complexity of statistical models such as the ones presented in this manuscript. Testing CI is therefore critical for the performance of such estimators.

ABSTRACT

We introduce the weighted partial copula function for testing conditional independence. The proposed test procedure results from the following ingredients: (i) the test statistic is an explicit Cramér-von Mises transformation of the weighted partial copula, (ii) the regions of rejection are computed using a boot-strap procedure which mimics conditional independence by generating samples. Under CI, the weak convergence of the weighted partial copula process is established and endorses the soundness of our approach. Experiment finally demonstrate the competitive power of our approach compared to recent state-of-the-art methods.

Notation

Equal by definition "Je suis abasourdi par le nombre de personnes qui veulent 'connaître' l'univers alors qu'il est déjà suffisamment difficile de se repérer dans le quartier chinois de New York." Woody Allen

With the advances in wireless communications and low-power electronics, accurate position location may now be accomplished by a number of techniques which involve commercial wireless services. Emerging position location systems, when used in conjunction with mobile communications services, will lead to enhanced public safety and revolutionary products and services. The fundamental technical challenges and business motivations behind wireless position location systems are described, and promising techniques for solving the practical position location problem are treated.

Context of this Thesis

This thesis is the result of a CIFRE agreement (Convention Industrielle de Formation et de Recherche -Industrial Training and Research Agreement) between Télécom Paris Saclay and Sigfox, a French telecommunications operator created in 2009 by Christophe Fourtet and Ludovic Le Moan. Sigfox is specialised in Machine to Machine (M2M) via low-speed networks. It contributes to the Internet of Things (IoT) by enabling interconnection via a gateway. Its UNB ("Ultra narrow band") radio technology enables it to build a low-speed, energy-efficient cellular network. This type of network is deployed in the so-called Industrial Scientific and Medical (ISM) radio frequency bands, available worldwide without any license.

In the last few years, Internet of Things (IoT) has raised a great deal of attention in very diverse fields such as agriculture or health care. Experts agree that 30 billions objects will be part of the IoT by 2023 and 40% of these objects will need to be geolocated, e.g. for freight transport (Hatton). One of the most significant challenges for the field is the need for localization. Indeed, numerous applications of sensor networks need to track mobile objects, such as people, animals, cars, etc. To make these applications viable, device cost will need to be low (from a few dollars to a few cents depending on the application) and devices will need to last for years or even decades without battery replacement. Additionally, the network will need to organize without significant human moderation. Moreover, in order to enable connectivity for billions of devices, most of IoT dedicated networks are using long-range and low power communications. Therefore, the challenges for IoT networks are to achieve high scalability to handle massive number of devices, to achieve low cost and to have wide coverage while keeping low energy consumption. Devices which have these demands are difficult to be integrated into traditional cellular networks. That is why, LPWAN (Low Power Wide Area Network) dedicated technology, such as UNB (Ultra Narrow Band), developed and patented by Sigfox has emerged.

The knowledge of the geolocation of each device is a very valuable resource. Indeed, it allows Sigfox to provide this information to network users that leads to numerous applications such as in logistics or transport of merchandises, to monitor and track in smart buildings or even for proximity marketing and advertising in shopping malls. However, traditional localization techniques such as global positioning system (GPS) are thus not well-suited for the particular needs of the IoT industry. Providing a GPS on each device is cost and energy prohibitive for many applications, and furthermore not suited for indoor applications. All these requirements complicate greatly the localization of these objects.

Alternatively, range-based methods use measurements such as the two introduced here, the Time-of-Arrival (TOA) or Time-Difference-of-Arrival (TDoA) [START_REF] Ho | Solution and performance analysis of geolocation by tdoa[END_REF]; [START_REF] Cong | Hybrid tdoa/aoa mobile user location for wideband cdma cellular systems[END_REF], and the Received Signal Strength Indicator (RSSI) to estimate the distance between an emitting device and a receiving antenna. Other ranging methods are commonly used in the literature, such the Angle of Arrival (AOA) [START_REF] Niculescu | Ad hoc positioning system (aps) using aoa[END_REF], or the Frequency Difference of Arrival (FDoA) [START_REF] Amar | Localization of narrowband radio emitters based on doppler frequency shifts[END_REF]) but their study is beyond the scope of this thesis.

This chapter allows us to put the stakes of this thesis into perspective. Its purpose is to provide a better understanding of the singularity of the geolocation problem in a sensor network such as the Sigfox network.

Outline In this chapter, we address the problem of channel-based location estimation techniques. First, we discuss the general problem of geolocation in Appendix B.2.

To predict the location of an emitter, state-of-the-art channel-based approaches consist either in the estimation of the time delay from the transmitter to the receiver, from which it is possible to infer the distance between the emitter and the receiver, or to directly predict this latter distance from the observed power decay between the received power and the emission power. A introduction to the wireless channel, in particular the key parameters for model it, and a presentation of the discussed predictors are proposed in Section 1.3. Finally, the performance of both these approaches in the particular context of IoT communication are discussed in Section 1.4.

The estimation of the geographic location

Geolocation Principle

The terms geolocation and positioning are used to designate the real-world geographic estimation of the location of an object. This problem has entered our society with a massive diffusion. It is used in many applications such as navigation, communication, self-driving vehicles, connected objects and communicating cities, or more recently with issues such as controlling the amount of contamination in a population. It affects a wide variety of scientific fields such as Géolocalisation et Navigation par un Système de Satellites (GNSS), on which our economy increasingly depends. Experts agree that 30% of the gross domestic product would depend partly on GNSS by 2030, compared to 10% today.

The geolocation of an object refers to the (latitude, longitude)-coordinates, that is the position of the object on the Earth's surface. Sometimes, the word positioning is rather employed when it comes to identify the location of an object in a particular space such as a cell-phones in a shopping mall, or a robot in a building. In this thesis, we shall speak about geolocation when it comes to a localization on a global scale, while the term positioning is generally used in indoor and/or confined areas.

The radio frequency methods are used for most systems of geolocation. This family of methods, also called radiolocation methods make use of characteristics of received radio waves to predict the location of an emitting object. The examples are numerous. The widely used Global Positioning System (GPS) (see Figure 1.1) is based on the estimation of the Time-of-Arrival (TOA) of a signal at a satellite. When the time of transmission, the speed of propagation and the position of the satellite are known, the TOA leads to a very good estimator of the distance between the emitting object and the satellite.

Combining several TOA leads to the estimated position of the transmitter. The use of multiple receivers to locate a transmitter is known as multilateration (illustrated in Figure 1.2). In cellular telephony, the radiolocation is directly performed via the Base Stations (BS) of the cellular network by means of one or a combination of the following features:

• The TOA (or TDoA). In contrast with the GPS, these quantities are estimated w.r.t the BS of the cellular network.

• The Angle of Arrival (AOA) corresponds to the direction from which the signal is received. A practical way to determine the AOA is to consider that this direction is the one of the maximum signal strength during a complete rotation of the BS.

Combining several AOA results in the desired position estimation.

• The RSSI of great interest in this thesis. It corresponds to the power of received signal strength minus the emitting signal strength. It yields to a ranging system by means of the Log-distance path loss model described in Section 1.3.2 or by fingerprint-based methods (when the different locations of emission are known to exhibit very different power "signatures"). These methods are studied in detail in Chapter 4.

Network-Based Geolocation

This thesis focuses on network-based geolocation. These methods only use the network infrastructure. Among all the described geolocation methods, the latter are the cheapest and require the least of energy. Methods of this kind have met a tremendous success with the apparition of the Internet of Things (IoT) in the late 1990's. Essentially, the concept of IoT is to provide to any objects the ability to transfer data over a network without requiring human-to-human or human-to-computer interaction [START_REF] Rouse | Internet of things (iot)[END_REF]. Nowadays, the set of applications for IoT devices is dramatic: smart homes (Samuel, 2016), medical and healthcare applications [START_REF] Catarinucci | An iot-aware architecture for smart healthcare systems[END_REF], agriculture (Mekala and Viswanathan, 2017) or even transportation systems [START_REF] Zhou | Design and research of urban intelligent transportation system based on the internet of things[END_REF]. We next move to two concrete geolocation examples encountered in practice.

Logistic Tracking cargo and assets via the network can also result to great benefits in transportation systems. One can send specific alerts when remarkable events occur such as the arrival in a warehouse.

Ecology It has then been argued that the IoT will revolutionize the ecology area. First, in terms of energy management: the connectivity of significant numbers of energyconsuming devices (e.g. lamps, motors, pumps, etc.) can allow them to communicate with utilities not only to balance power generation but also optimize the energy consumption as a whole. Second, environmental monitoring applications of the IoT use sensors to help environmental protection by monitoring as instance air or water quality.

Other applications like earthquake or tsunami prediction systems can also be used to provide more effective aid.

A natural idea is to apply the radiolocation methods described above to the Sigfox IoT dedicated network. As stated before, the performance of such methods strongly relies on the network infrastructure. Without going into details here, we propose to give some elements that compromise their use and therefore motivate the subsequent works. First of all, Sigfox BS have a lack of directivity and thus cannot discriminate the incident wave directions. Therefore, the AOA-based methods are irrelevant here. In addition, Sigfox bases its communications on Ultra-Narrow-Band (UNB) Technology which allows to achieve both long range and extended battery life. Thus, each signal has a frequency band of 100 Hz width within the unlicensed frequency band (ranging from 868.0 to 868.6 MHz). This band is popular because it presents a good balance of range, building penetration and the ability to use small antennas. Nevertheless, it is well known that the performance of time-based methods using TOA or TDoA estimation strongly relies on the signal bandwidth. The Cramér-Rao Lower Bound (CRLB) on TOA estimators is thus often exhibited to quantify this effect. It states that the minimal variance of any unbiased estimator of the TOA is inversely proportional to B 3 × SNR, where B stands for the bandwidth and Signal to Noise Ration (SNR) is the well-known Signal-to-Noise ratio. This bound is then unfavourable to the use of such approaches and yields for instance to a standard deviation of range estimators at least equal to 20 km with Sigfox US infrastructure. Nevertheless, it also constitutes a first benchmark value to which our proposed geolocation methods are compared. It is shown [START_REF] Boucher | Analysis of discrete implementation of generalized cross correlator[END_REF] that another disfavourable element is that they are particularly memory demanding since their performance is directly related to the choice of the sampling interval.

RSSI-based Geolocation

The RSSI is a relevant feature for the geolocation task pursued. It is sufficiently explanatory without constituting a memory burden. The related literature is extremely vast and is often classified into two categories of methods: range-based and range-free. They essentially differ in the information used for the localization task. The first category of methods use measurements to predict the range between emitters and receivers (the coordinates of the BS are known in this case). Combining several estimated ranges allows the emitter position estimation by means of some mathematical methods such as trilateration [START_REF] Thomas | Revisiting trilateration for robot localization[END_REF]. There are thus always two phases: a ranging phase and a localization one. In contrast, the second family of methods do not base their predictions upon the range estimation. A simple example of such a method would estimates the position as the barycenter of the receiving BS. As far as the RSSI is the only measurements at our disposal, the range-based methods quickly show their limits [START_REF] Gayathri Chandrasekaran | Empirical evaluation of the limits on localization using signal strength[END_REF]. Indeed, they essentially found upon the log-distance path loss model of which the lack of realism is brought out by the investigation of Section 1.3.2.

In this thesis, we propose another classification of RSSI-based geolocation methods: likelihood-based and fingerprint-based methods. The first ones, that encapsulate the range-based methods already discussed, consist of learning (on a dataset) a model for the RSSI at a BS (denoted by X in the sequel) given the position (denoted by Y in the sequel). This learning phase is also referred to as the "calibration phase" in the literature. Once this model has been inferred, one can predict the emitter position as the one that best agrees the measured RSSI. The second ones are the fingerprint-based methods. The latter directly map the RSSI to the position by means of a function that has been previously learned on a dataset.

The recent advances in machine learning and its successes in a wide range of areas have driven the IoT community to apply these methods to RSSI-based geolocation. Let us now present the contributions we have developed to address this issue.

State-of-the-Art Gelocation Techniques

A formal introduction to the wireless channel and to the physical parameters involved in its modeling is necessary to understand the location estimation techniques developed in this chapter. Following Tse and Viswanath (2005), we start with the physical modeling of the wireless channel in terms of electro-magnetic waves and then define some important physical parameters. This will allow a better understanding of the foundations of the two geolocation methods discussed in this chapter.

Physical Modeling for Wireless Channel

First, let consider a fixed emitter. If the receiving antenna is sufficiently far away from this emitter (which we refer as the far field) the electric field and magnetic field at any given location are perpendicular both to each other and to the direction of propagation.

In response to a transmitted sinusoid cos 2πf t, the electric far field at time t can be written as follows:

E(f, t, (r, ϑ, ψ)) = α s (ϑ, ψ, f ) cos 2πf t -r c r , (1.1)
where, as shown in Figure 1.3, (r, ϑ, ψ) represents the polar coordinates of the point M at which the electric field is measured; r is then the distance from the emitter to M , (ϑ, ψ) are respectively the vertical and horizontal angles from the emitter to M . The constant c is the speed of light, and α s (ϑ, ψ, f ) is the radiation pattern of the sending antenna at frequency f in the direction (ϑ, ψ).

We observe that, as the distance r increases, the electric field decreases as r -1 and thus the power per square meter in the free space wave decreases as r -2 as illustrated in Figure 1.4. This attenuation is referred in the literature as the large-scale fading.

Nevertheless, we will see that this r -2 reduction of power with distance is not valid as soon as there are obstructions or reflections in the propagation space. Thereupon, two methods to estimate the unknown distance from the emitter emerge:

• First, a straightforward method that consists in mapping the received power to the desired distance. In free space, the signal power decays proportionally to r -2 . In real-world channels, multi-paths (due to reflectors) and shadowing are two major sources of environment-dependence in the measured power. These phenomena are generally modeled as random. As a consequence, the power decay is rather chosen proportionally to r -α , where α is the so-called path-loss exponent. Once this path-loss exponent has been estimated (through the calibration phase), it is possible to provide a range estimate from the measured power signal strength.

The procedure is described in the sequel.

• On the other hand, an intuitive method is to estimate the delay between the emission and the reception of the signal. To provide this Time Of Arrival (TOA) estimate, we suppose that the signal emanating from the source is known, and that the received signal is a superposition of an attenuated and delayed replica of the known signal plus noise. The remaining task is thus to estimate the arrival times of the replica and its attenuation amplitude. This is the well-known Time-of-Arrival estimation problem. There is a vast literature dealing with this problem. The most well-known estimators are the matched filter approach described in Ehrenberg et al. (1978), or the maximum likelihood (ML) estimator introduced in the work of [START_REF] Knapp | The generalized correlation method for estimation of time delay[END_REF]. This latter estimator is obtained as the peak of the cross-correlation function of both the emitted and the received signal. The Cramér-Rao Lower Bound (CRLB) (see Chapter 2) is taken as the standard of reference, because it should give a realistic indication of the attainable mean square error of the estimators. This problem is close to the problem of Time-Delay of Arrival (TDoA) or simply Time-Delay estimation. This latter aims to predict the time delay between the receptions of the signal at two different antennas. It allows to circumvent the imperfect knowledge of the emission time due to the lack of both precision and synchronisation of the emitting devices inner clocks. As far as location is concerned, the first method requires solving circles intersection equation while the second an hyperbolas intersection equation. We thus chose to present this problem here-below rather than the ToA problem.

In the sequel, we will discuss this two ranging-methods, as well as their theoretic reachable precision through their CRLB. For this purpose, we first introduce the singular context of IoT dedicated networks, particularly the Ultra Narrow Bound technology and its impact on the lower bound of the range estimate.

Location Estimation Techniques

The Time-Delay problem

The parameters of primary interest are the relative delays between receiving antennas which determine range (as well as the bearing) to the source. We assume hereafter that the transmission medium is homogeneous so that the signal wave-front is perfectly coherent over the receivers: the signal components received by various sensors are then delayed replicas of each other. The noise is considered to be additive, white, Gaussian, and incoherent from sensor to sensor.

Following the seminal work of [START_REF] Clifford | Coherence and time delay estimation[END_REF], we can derive a location estimation technique based on the time delay estimation. Let first consider the case of two receiving antennas at which the received signal can be modeled as follows:

x 1 (t) = s(t) + n 1 (t), (1.2) x 2 (t) = As(t -D) + n 2 (t), (1.3)
where it is assumed that s(t), n 1 (t), and n 2 (t) are real and jointly stationary random processes, s(t) is furthermore independent from n 1 (t) and n 2 (t).

The maximum likelihood estimator of the time delay maximizes the following crosscorrelation function (CCF).

R x 1 ,x 2 (τ ) x 1 (t)x 2 (t -τ )dt, (1.4)
We obtain the following Cross-Correlation (CC) estimator:

DCC arg max τ R x 1 ,x 2 (τ ), (1.5)
Carter and Knapp (1976) proposes a improved version of the CC estimator by introducing the so-called generalized cross-correlation (GCC) function. This function is obtained by applying pre-filters to amplify spectral components of the signal that have little noise and attenuate components with large noise. As such, the GCC requires knowledge (or estimates) of the signal and noise power spectra. There is a number of algorithms in the GCC family depending on the choice of the filters. Commonly used weighting functions include the constant weighting (in this case, the GCC becomes a frequency domain implementation of the cross-correlation method), the smoothed coherence transform (SCOT) [START_REF] Carter | The smoothed coherence transform[END_REF], the Hassab-Boucher transform [START_REF] Hassab | Performance of the generalized cross correlator in the presence of a strong spectral peak in the signal[END_REF] or the maximum-likelihood (ML) [START_REF] Knapp | The generalized correlation method for estimation of time delay[END_REF]. It is well known that the ML estimator obtained in the ideal propagation situation is optimal from a statistical point of view since the estimation variance can achieve the CRLB which has the rather simple form :

σ 2 CRLB = 3 8π 2 1 + 2SNR SNR 2 • 1 B 3 T , (1.6)
where B is the signal bandwidth, T is the observation time, and SNR is the Signalto-Noise Ratio defined as the ratio of the power of a signal (meaningful input) to the power of background noise (meaningless or unwanted input In other words, although the CRLB is a promising bound for the variance of our time delay estimator, it is not reachable in practice because it is obtained by too idealistic assumptions: the observation sample space has to be large enough; the environment should be free of reverberation; and the spectra of noise signals have to be known a priori. In practice, other problems may also arise. All the techniques described above measure the analog processing of the time delay parameter, while discrete signal processing methods are more commonly used. In the sequel, we briefly discuss the case when we only have at our disposal sampled signals. The problem remains the same: finding an estimator D of the true delay D using a finite set of samples of x 1 (t) and x 2 (t) using as in the continuous case, the peak of the CCF. A normal practice is to replace the CCF defined in Equation (1.5) by its time-averaged estimate:

Rx 1 ,x 2 (τ ) = 1 N N k=1 x 1 (kh)x 2 (kh -τ )dt, (1.11)
where h is the sampling interval, and (N -1)h is the estimation window width. The socalled Direct Correlator (DC) estimator of the time delay is defined as the maximizer of the latter quantity w.r.t. τ . It is thus necessary to interpolate the CCF in the neighborhood of the peak. A commonly used approach [START_REF] Boucher | Analysis of discrete implementation of generalized cross correlator[END_REF]) is to fit a parabola in the neighborhood of this maximum using 3 samples of Rx 1 ,x 2 (•). Thus, near its peak the CCF can be approximated by the following convex parabola:

Rx 1 ,x 2 (τ ) aτ 2 + bτ + c,
where (a, b, c) are the fitted parameters. The apex of the parabola is then used as a proxy of the time delay:

DCC = - b 2a .
(1.12)

It is shown that this methods yields a biased estimator of the time delay D [START_REF] Jacovitti | Discrete time techniques for time delay estimation[END_REF]. In general, this bias can be decomposed into a systematic bias, due to parabolic estimation; and another one resulting from the noise and a finite width of signals observation. The results of these works are enlightening in many ways. First because they state that in the range of SNR that are observed in the Sigfox network, the second bias is negligible with respect to the first one. The parabola misfit is therefor the keystone of these discrete methods. Moreover, the found expression for the bias is the product of two terms. A first term which depends linearly to the the sampling interval h. Let us recall that, for Sigfox, it is chosen as h = (1/2B) = 0.005 s. A second term is a function which depends only on the reference signal. Simulations of theoretical bias provided in these two works are shown in the [START_REF] Jacovitti | Discrete time techniques for time delay estimation[END_REF] where E[ DCC ] is computed for various D. As expected, the bias is the higher when the true TD lies between two samples points. The simulation results, obtained using the quantities observed in the Sigfox network invalidate these approaches because results to low time resolution.

The Received-Signal Strength and Log-distance Path Loss (PL) model

In telecommunications, Received Signal Strength Indicator (RSSI) is a measurement of the power present in a received radio signal. Wireless sensors communicate with neighboring sensors, and RSSIs can be measured by each receiver without any further bandwidth or energy requirements. Because RSSIs are relatively inexpensive and simple to implement in hardware, they are an important and popular topic of localization research. Typically, the ensemble mean received power in a real-world, obstructed channel decays proportional to r -α , where α is the "path-loss exponent", typically between 2 and 4 as acknowledged in Rappaport et al. (1996). Log distance Path Loss (PL) model is an extension to the Friis free space model introduced in Section 1.3.1. It is used to predict the propagation loss for a wide range of environments and to therefore overcome the limitation of the Friis free space model. The model encompasses random shadowing effects due to signal blockage by trees, buildings etc. This model is a parametric model of the RSSI given the distance d. The vector of parameters θ = (P 0 , α, d 0 ) is such that d 0 is some reference distance and the parameter P 0 represents the power in dBm at distance d 0 .

Using loose notation, given the distance d between the emitter and the receiving antenna, the RSSI X at this receiver is distributed as follows [START_REF] Mao | Path loss exponent estimation for wireless sensor network localization[END_REF]:

X[dBm] ∼ N P 0 -10α log d d 0 , σ 2 dB .
(1.13)

Estimating the range from the RSSI The estimated distance between the receiver and the emitter, can be thus estimated from the realization of the variable X. We recall that the log-likelihood of X given d is,

log p(X|θ, d) = c - X -P 0 + 10α log d d 0 2 2σ 2 dB , (1.14)
where c is independent of θ. The distance which maximizes the likelihood is

dML (X) = d 0 10 P 0 -X 10α . (1.15)
We can prove that the bias of the ML estimator is as follows:

E d 0 10 P 0 -X 10α = exp 10α √ 2σ dB log 10 2 C d.
In typical channels studied e.g. in Rappaport et al. (1996), C 1.2, and the range is then over estimated by a factor 20%. It is therefore preferable to use the bias-corrected estimator as follows:

dCML (X) = d 0 exp - 10α √ 2σ dB log 10 2 10 P 0 -X 10α .
(1.16)

Remark 1.1. Note that the variance of the bias-corrected estimator is proportional to the actual range. This result is however not surprising, and above all, informs on the need to have a very dense network in order to have a quality estimation of the emitter.

Introduction to Sigfox Radio System

Two of the most commonly used techniques for geolocation have been discussed above. This section provides a better understanding of the singularity of the geolocation problem in a sensor network such as the Sigfox network. The performance of the two described geolocation approaches in Sigfox network are studied.

Internet of Things, a new usage for the radio-communications

There is no consensus on the definition of the Internet of Things (IoT), although its initial use has been attributed to Kevin Ashton. In fact, many groups including academicians, researchers, developers and corporate people have their own. Yet, all the definitions have in common the idea that the first version of the Internet was about data created by people, while the next version is about data created by things. We can therefore cite [START_REF] Madakam | Internet of things (iot): A literature review[END_REF] and give the following definition: "An open and comprehensive network of objects that have the capacity to auto-organize, share information, data and resources". The connection between objects is not to be understood in terms of servers, computers or smartphones. In the context of IoT, sensors embedded in physical objects (from roadways to pacemakers) are linked through wireless networks. This implies that the implementation of IoT communication must be cost-effective and energy frugal since the object may not have been natively designed with an energy source. As a consequence, IoT comes with several constraints for the radio-communication system:

• The communication function of a connected device must have a marginal cost (when it is compared to the cost of the device itself). For ease of reference, the target price of a device for mass production should be less that 0.2 $.

• The volume of transmitted data by the connected devices must be low. Indeed, devices transmit information originated from sensors, alarms or GPS. trackers which communicate through small "packets". In Sigfox's network, the size of the messages must not exceed 12 bytes (see Section 1.4.2 for details).

• The volume of connected objects might be huge. Due to the low-cost and simplicity of use, it is possible to connect a large number of objects. This may result to a very high density of connections by square kilometer. Typically, in Sigfox network, the density of connected objects reaches 50, 000 per square kilometers (that is much higher than in cellular networks).

We thus note that in order to make use of IoT, many requirements come into play. The radio technology that can support all these requirements is encapsulated into the name of Low-Power Wide-Area Network (LPWAN 1 ), where we refer to [START_REF] Raza | Low power wide area networks: An overview[END_REF] for a in-depth study.

Ultra-narrowband benefits for LPWAN

While IoT is an expression to describe the usages, LPWAN rather encapsulates the underlying radio technology. The aim of this technology is to meet two contradictory needs: the low-energy consumption and the wide-area network. Therefore, LPWANs have the following characteristics:

• To reduce the infrastructure cost, the base stations of the network must be limited in number. The network counts a very large number of connected devices compared to the number of BS. Thus, one BS can receive numerous messages originated from multiple objects (sometimes at the same time).

• To reduce the energy cost of communication, the connected objects must transmit at very low power. They can therefore be endowed with small batteries (or even no battery at all thanks to energy harvesting). High sensitivity of the BS is thus a need to achieve a wide-area coverage with one-hop communication.

1 Machina Research, 2013

We now discuss briefly about Ultra-Narrow Band (UNB) and how this technology can meet the requirements described above.

UNB Technnology

Ultra Narrow-Band systems are those in which the channel has a very narrow bandwidth. This band is significantly smaller than the total available frequency resource, typically of few hundred Hertz. Essentially, Sigfox's technology is such that the receiver listens for a large spectrum range with a high dynamic range and the software is continuously looking for the signal of interest and tracks its center frequency, as a human radio operator would do (Chaudhari et al., 2020). In the Sigfox's network, the transmitted signal of UNB occupies a band of 100 Hz (in Europe, and 600 Hz in USA), inside a typical possible band of 192 kHz to 2 MHz. In Europe, the frequency band ranges from 868.0 MHz to 868.6 MHz and thus brings up to 6000 sub-channel of 100 Hz width each. Objects now can randomly select carrier frequency, without excessive collision rate. We assume that all objects have the same behavior and transmit a message of duration T seconds at a frequency f s . For the time-slotted case, any active node randomly selects a time-slot, since the temporal resource is separated by slots. Thus in the time domain, the transmissions either do not collide at all (when they choose different time-slots), or collide for the whole duration (when they choose the same time-slot).

All transmissions are performed within a dedicated band, which has band-width B Hz. Each transmission occupies a bandwidth b Hz which represents a small portion of the total channel bandwidth B Hz. From the base-station point of view, the total dedicated band contains, from time to time, transmitted ultra narrow signals at random carrier frequencies. For each detected transmission, the BS extracts the signal at the estimated frequency of interest, and decodes the packet. Such a detection and estimation can be done as described in a Sigfox patent (Artigue, 2017).

Campaign of validation of the PL model in Sigfox Network

We now conclude this section with a numerical investigation. We give numerous practical examples that invalidate the use of the log-distance path loss model in rural environment context of our concern. First, the log-normal distribution of the RSSI measurements is verified by examining the residuals X -(P 0 +10α log d d 0 ) using quantile-quantile plots shown in Figure 1.5. This figure shows that the observed data have a heavier tail than the one of the normal distribution: the quantiles move away from the diagonal outside a certain range. In that respect, the next numerical experiment investigate the lack of realism of the studied model. More specifically, the non-isotropy of the propagation model is highlighted. Indeed, these experiments show the high dependency of the RSSI upon the direction of the arrival signal. This is shown in Figure 1.6. This model of propagation, assuming that there is a free propagation space, and thus omitting the presence of obstacles in certain directions, is vain. A first extension of the model, is to encapsulate this dependence of the loss by making the path loss exponent α(ϑ, φ) direction-dependent. This approach is however interesting because it amounts to learning a path-loss exponent for every direction (ϑ, φ). Instead, a very general semiparametric model is introduced in Section 3.3.1. This model provides an estimation of the expected RSSI at a BS given the emission position (x, y). This enables to free ourselves from the knowledge of the positions of the BSs, which would have been needed if the dependence upon the direction between the emitter and the BS was modeled. This study has put into light how the context of IoT network jeopardizes the precision of the two presented methods of geolocation. It also identify more clearly the interests of using recent methods of machine learning to the problem of network-based geolocation. We therefore dedicate the future chapters to an intensive investigation of these approaches.

Geolocation as a Prediction Learning Task

The network is composed of d fixed BS, say (BS 1 , . . . , BS d ), whose respective coordinates (y 1 , . . . , y d ) in the complex plane are fixed but not necessarily known (unless it is specified).

We consider a connected device whose position Y is a random variable in some given subset Y, typically an open subset of R 2 . The device sends packets/messages which are collected by the neighboring BS. For a given message, each BS k (k = 1, . . . , d) computes a RSSI X k as the temporal mean of the received signal strength. The RSSI X k is typically real-valued in a certain subset X ⊂ R. However, some messages may not be detected by some BS, in which case we just set X k = NaN, where NaN stands for an unobserved value. Thus, for every k = 1, . . . , d, X k is a random variable in the set X X {NaN}. The aim of location estimate is then to predict the unknown position Y from the observation of the RSSI-vector

X (X 1 , . . . , X d ) .
It is the general prediction learning task: we wish to infer the relationship between the random variable X and a target random variable Y , taking values in sets X d and Y respectively. For a loss function : Y × Y → R + that defines a discrepancy on Y (in typical settings, (y, ŷ) = y -ŷ 2 ), the generic writing of this task is then to find a function h : X d → Y that verifies: (h(x i ), y i ).

h = arg min h∈H E (h(X), Y ) , ( 1
(1.18)

Equipped with this mathematical background, we now move to our application focus: the problem of network-based geolocation.

Contributions

The organization of this manuscript is as follows. Each chapter can be read independently.

Chapter 2 introduces the theoretical elements needed to the reading of this thesis. First, elements of statistical inference are introduced. Then, the notion of risk, and a lower bound on the so-called quadratic risk on an estimator is discussed. These notions will prove themselves useful to highlight that the parametric family of estimators are not well suited when applied to RSSI-based geolocation, and that we need to go beyond simple parametric estimation. This last point is therefore studied in this chapter as well.

Chapter 3 investigates machine learning approaches addressing the problem of geolocation. First, we review some classical learning methods to build a radio map. These methods are split in two categories, which we refer to as likelihoodbased methods and fingerprinting methods. Then, we provide a novel geolocation approach in each of these two categories. The first proposed technique relies on a semi-parametric Nadaraya-Watson (NW) estimator of the likelihood, followed by a maximum a posteriori (MAP) estimator of the object's position. The second technique consists in learning a proper metric on the dataset, constructed by means of a Gradient boosting regressor: a k-nearest neighbor algorithm is then used to estimate the position. The proposed methods are compared on two data sets originated from Sigfox network, and an indoor dataset performed in a threestory building. Experiments show the interest of the proposed methods, both in terms of location estimation performance, and ability to build radio maps.

Chapter 4 is dedicated to the similarity learning. Indeed, the choice of the similarity (or distance) is known to be critical for the performance of neighbors-based predictor such as k-NN regressor. That is why there is an increasing interest for optimizing distance and similarity functions. However, in most prior works, no link is made between the learned metrics and the estimator performance. In this chapter, we propose to build the metric by directly minimizing the regression error of our estimator, and thus obtain an ad-hoc learning objective. To minimize this objective, we propose a modified version of the eXtreme Gradient Boosting algorithm (XGBoost). The soundness of our approach is finally endorsed by conclusive numerical experiments on numerous datasets.

Chapter 5 introduces the weighted partial copula function for testing conditional independence (CI). The CI has been broadly used in the RSSI-based geolocation literature in order to decrease the statistical models complexity. Testing this assumption can leads to a better understanding of the likelihood based models proposed herein. The proposed test procedure results from the following ingredients: (i) the test statistic is an explicit Cramer-von Mises transformation of the weighted partial copula, (ii) the regions of rejection are computed using a bootstrap procedure which mimics conditional independence by generating samples from the product measure of the estimated conditional marginals. Under conditional independence, the weak convergence of the weighted partial copula process is established when the marginals are estimated using a smoothed local linear estimator.

Publications

The works presented in this manuscripts have resulted in several accepted publications and preprints, that are listed here in chronological order:

• Elgui, K., Bianchi, P., Portier In this chapter, we introduce the theoretical elements needed to the reading of this thesis. First, we properly define elements of statistical inference that are afterwards used for RSSI-based methods of geolocation. Then, the notion of risk, and a lower bound on the so-called quadratic risk on an estimator is discussed. These notions will prove themselves useful to highlight the irrelevance of the parametric family of estimators (when applied to RSSI-based geolocation), and the need to go beyond simple parametric estimation. This last point is therefore studied in this chapter as well.

Elements of Statistical Inference

Introduction

The problem of estimating statistical models depending on a finite number of parameters goes back to Fisher (Fisher, 1925) and has met with tremendous success ever since [START_REF] James | An introduction to statistical learning[END_REF].

Parametric approaches have the advantage of being very easy to compute, hence their predominance in the related literature. Furthermore, powerful tools are available for the statistical analysis of such estimators, such as a lower bound of their variance. This kind of knowledge, necessary to assess the quality of an estimator, is unfortunately difficult to obtain when going beyond parametric model, and often requires complex simulation methods. However, parametric models are limited in their expressive power, as they only provide an approximation, often imprecise, of the underlying statistical structure. Statistical models that explain the data in a more consistent way are often more complex; unknown elements in these models are usually, instead of scalar parameters, functions having certain smoothness properties. In this thesis, we shall indeed see that parametric estimation does not provide sufficiently satisfactory results when applied to RSSI-based location estimation. Hence, semi-parametric or non-parametric methods will often be preferred. This section is structured as follows.

Outline First, we set in Section 2.1.2 the framework of parametric estimation. This framework allows to introduce the Cramér-Rao Lower Bound of an estimator. This naturally leads to the introduction in Section 2.2 of the non-parametric estimation framework through the two main problems. Practical aspects of the regression problem are then investigated, and Tree-Based methods are thus introduced in Section 2.2 as a conclusion of this chapter.

Parametric Estimation

Statistical Model

We first consider a family P of probabilities over the space of observations (X , B(X )).

We refer to the two first chapters of Shorack and Shorack (2000) for a thorough definition of these mathematical objects. The family P is called the statistical model.

Then, we can define a parametrization of the model by attaching each P ∈ P with a parameter θ ∈ Θ, where Θ is any open subset of R d . One will then write P θ to designate the law thus labeled.

Definition 2.1 (Statistical model). We call a statistical model a family P of probabilities P over the space of observations (X , B(X )). If there exists a set Θ set such that

P = P θ , θ ∈ Θ ,
then, Θ is called the parameter space, and P a parametric model.

A model P θ , θ ∈ Θ is said to be dominated if any P θ ∈ P is absolutely continuous with respect to a common reference measure µ. In the sequel, we place ourselves in the case where the model is dominated. As a consequence, for all θ ∈ Θ, P θ always admits a probability density function p θ w.r.t the reference measure µ.

This fact allows to directly work on a family of densities rather than on a family of distributions. This leads to the introduction of the notion of likelihood, that we define as follows:

Definition 2.2 (Likelihood). The likelihood of the observation x is the mapping θ → p θ (x). This likelihood is also denoted p(x; θ).

The likelihood is a central element of a family of inference methods (e.g. maximum likelihood, or likelihood methods), for which the intuition is clear. Considering an observation x, the likelihood "assesses" how likely it is that this observation has been generated under P θ . The definition of the Maximum Likelihood estimator follows naturally:

Definition 2.3 (Maximum Likelihood). We consider a dominated statistical model P = P θ , θ ∈ Θ , p(•; θ) the probability density of P θ w.r.t. the reference measure µ. Suppose that we observe x drawn from P θ . A Maximum-Likelihood (ML) estimator of θ, is any estimator θ(x) that verifies: p(x; θ) ≥ sup t∈Θ p(x; t).

(2.1)

Applications of these estimators are considered in Chapter 3, in order to estimate the time-delay given the observation of both the emitted and received signals, or to estimate the range given the observation of the RSSI.

Risk and Cramér-Rao Lower Bound

We first recall the notion of risk for a parameter estimate. Then, lower bounds on the quadratic risk of unbiased estimators are established. In this thesis we focus on the estimation of parameters of interest that write as g(θ) for some function g : Θ → R.

An example is the Time Delay estimation in Section 1.3.2. Let now ĝ be an estimator of the parameter of g(θ). The risk evaluated at ĝ is defined as follows:

R(θ, ĝ) E θ g(θ), ĝ(X) , (2.2)
where : R × R → R + is some loss function. In typical settings, is the quadratic loss, and the corresponding risk is then called the quadratic risk, or Mean Squared Error (MSE). The MSE can be decomposed as follows:

MSE(θ, ĝ) E θ (g(θ) -ĝ(X)) 2 = E θ (ĝ(X) -E θ ĝ(X)) + (E θ ĝ(X) -g(θ)) 2 = E 2 θ ĝ(X) -g(θ) + E θ (E θ ĝ(X) -g(θ)) 2 b 2 (θ, ĝ) + Var θ ĝ(X) . (2.3)
This formula is known as bias-variance decomposition. Its simplicity is only provided by the use of the quadratic loss. The family of estimators that satisfy b(θ, ĝ) = 0 for all θ ∈ Θ is called the class of unbiased estimators.

In the sequel, we discuss the quadratic risk of such unbiased estimators via the notion of Fisher Information. We begin with the case where d = 1, that is θ ∈ Θ ⊂ R, to simplify the presentation. We consider an unbiased estimator of g(θ), i.e., a statistic ĝ : X n → R such that E[ĝ(X)] = g(θ) for all θ ∈ Θ. The Cramér-Rao bound gives a lower bound on the variance of ĝ, hence on its quadratic risk (since ĝ is unbiased). This theorem is defined under certain technical conditions on the statistical model at stake. These regularity conditions of the model are provided e.g. in [START_REF] Shao | Mathematical statistics: exercises and solutions[END_REF]. When these conditions are met, the model is said to be regular. We then define the Fisher information quantity [START_REF] Dembo | Information theoretic inequalities[END_REF]:

I(θ)      ∂ log p ∂θ (X; θ) 2      . (2.4)
Theorem 2.4. Let P = P θ , θ ∈ Θ be a regular model, g(θ) the parameter of interest and let ĝ(X) be a statistic such that E θ [ĝ(X)] = g(θ) and Var θ [ĝ(X)] < ∞ for all θ ∈ Θ. Suppose furthermore that 0 < I(θ) < ∞. Then,

Var θ [ĝ(X)] ≥ g (θ) 2 I(θ) . (2.5)
The proof is beyond the scope of this section, but we refer the interested reader to Kay (1993).

Corollary 2.5. Suppose that the conditions of Theorem 2.4 are satisfied and that ĝ is an unbiased and regular estimator of the parameter. Then,

Var θ [ĝ(X)] ≥ I -1 (θ), ∀θ ∈ Θ. (2.6)
This bound is called Cramér-Rao Lower Bound (CRLB) or Darmois-Fréchet bound.

This bound is applied to the case when θ is a scalar. The following proposition is an extension to the multi-dimensional case.

Proposition 2.6. Let ĝ be a statistic such that for every θ ∈ Θ, Var θ [ĝ(X)] < ∞, with 0 < I(θ) < +∞. We furthermore denote by g(θ

) = E θ [ĝ]. Then, θ → g(θ) is differentiable and: Var θ [ĝ(X)] ≥ ∇ θ g(θ) I(θ) -1 ∇ θ g(θ).
(2.7)

Non-Parametric Estimation

The problem of nonparametric estimation consists in estimating, from some observations, an unknown function belonging to a class of functions. We focus on two nonparametric estimation problems in this thesis: the estimation of a density and that of nonparametric regression Chapters 3 and 5.

An introduction to nonparametric methods is therefore proposed herein.

Kernel Density Estimation

Let X 1 , . . . , X n be n i.i.d. real-valued random variables whose common distribution is absolutely continuous w.r.t. the Lebesgue measure on R d . The problem of this section is to estimate the density of the common (unknown) distribution p. The estimate of p is thus a measurable application w.r.t. the observations (X 1 , . . . , X n ), pn :

x → pn (x; X 1 , . . . , X n ) = pn (x).
It happens that some prior can be known about this density function, and we can look for p in a parametric family. This case has been discussed in Section 2.1.2. But most of the time, this prior is not available and one may assume that p belongs to an hypothesis set of function H. This set can be the set of all Lipschitz continuous probability densities for instance [START_REF] Alexandre | Introduction to nonparametric estimation[END_REF]. In this thesis, we focus particularly on kernel density estimation methods. The kernel density estimate of p, also called the Parzen window estimate [START_REF] Parzen | On estimation of a probability density function and mode[END_REF], is a non-parametric estimate given by:

pn (x) = 1 n n i=1 K h (x -X i ), (2.8)
where

K h : R d → R + is defined as K h (x) = h -1 K(h -1 x) with h > 0, and K : R d → R + a kernel function, satisfying the condition R d K(u)du = 1.
(2.9) Popular kernels include for example:

• the rectangle kernel:

K(u) = 1 2 1 [-1,1] (|u|).
• the Gaussian kernel:

K(u) = 1 √ 2π exp -u 2 /2 .
• the Epanechnikov kernel:

K(u) = 3 4 (1 -u 2 )1 [-1,1] (|u|).
Multidimensional estimator can be easily obtained from the latter one-dimensional kernels. Consider for instance the case d = 2: we suppose that we have a sample of n i.i.d. pairs of r.v.

(U 1 , V 1 ), . . . , (U n , V n ), with a common density p(u, v) in R 2 . Then,
the kernel estimator of p(u, v) is given as follows:

pn (u, v) = 1 n n i=1 K h 1 (u -U i )K h 2 (v -V i ),
where h 1 and h 2 are the bandwidths associated to the "selection" of the neighbourhood of u and v respectively.

Nadaraya-Watson Estimator

We now study the case of non-parametric regression. We introduce the methods of construction of well-known Nadaraya-Watson estimators [START_REF] Elizbar | On estimating regression[END_REF], and to give guarantees on the risk of such estimators.

We consider a pair (X, Y ) of real-valued random variables, where Y is integrable:

E|Y | < ∞. We call regression function of Y given X the function m : R d → R : m(x) = E[Y |X = x].
(2.10)

We have at our disposal a sample (X i , Y i ) i=1...n of n, independent and identically distributed (i.i.d) pairs of random variable, with the same distribution as (X, Y ). The aim is to build from this sample an estimator of the regression function m.

Given a kernel K and a positive bandwidth

h > 0, for x ∈ R d , put K h (x) = h -1 K(h -1 x).
It follows easily that K h still verifies Equation (2.9). The Nadaraya-Watson estimator of the regression function m reads as follows:

m NW n (x) = n i=1 Y i K h (X i -x) n i=1 K h (X i -x) , if n i=1 K h (X i -x) = 0 , (2.11)
and m NW n (x) = 0 otherwise. To conclude this section, we briefly discuss the intuition beyond the expression of this estimator. We assume that the distribution of the random variable (X, Y ) admits a density p(x, y) w.r.t. the Lebesgue measure, and that p(x) = R p(x, y)dy > 0 a.s. We can write

E[Y |X = x] = R yp Y |X (y|x) dy = R yp(x, y) p(x, y)dy dy.
(2.12)

The following proposition shows that the Nadaraya-Watson estimator of the regression function is simply obtained by replacing in Equation (2.12) the unknown probability density functions p(x) and p(x, y) by their estimated versions pn (x) and pn (x, y) introduced in Section 2.2.1.

Proposition 2.7. Let pn (x) and pn (x, y) be the kernel estimates of p(x) and p(x, y) respectively as defined in Equation (2.8), then:

m NW n (x) = R y pn (x, y) pn (x) dy.
The proof can be found in [START_REF] Alexandre | Introduction to nonparametric estimation[END_REF].

Remark 2.8. The Nadaraya-Watson estimator m NW n satisfies the following property:

m NW n (x) = arg min a∈R n i=1 (Y i -a) 2 K h (X i -x). (2.13)
Thus m NW n is obtained by a local constant least squares approximation of the outputs Y i .

We are interested in the point-wise bias of such estimate E[m NW (x)] -m(x).

Remark 2.9. We draw to your attention that we chose to conserve the writing in its standard form, that is for the estimation of

E[Y |X = •].
However, in Section 3.3.1, Nadaraya-Watson estimators are in fact, employed to estimate X given the position Y . Indeed, the key question is whether it is possible to choose an appropriate kernel on X × X . One can either build an ad-hoc kernel on the space of the RSSI (see our proposed solution in Chapter 4), or leverage the fact that the choice of kernel is much easier on Y × Y. The latter approach implies the use of Maximum a Posteriori (MAP) estimator of Y , described in the next section.

Maximum A Posteriori estimators

We focus in this section on the Maximum A Posteriori (MAP) estimator of the target Y given X. In this thesis, this approach is at the center of the contribution proposed in Section 3.2.1.

First, we denote by p y , a prior on Y . In the case of a dominated model, one can determine the distribution a posteriori by explicitly writing its density. Suppose that P y , y ∈ R is a dominated model with P (dx|y) = p(x|y)µ(dx) and let ν be a dominant measure of the distribution p(•|x), and continue to note p Y the density, p(dy) = p(y)ν(dy). The joint density of the random vector (X, Y ) with respect to the product measure µ d ⊗ ν is then given by : p(x, y) = p Y (y)p(x|y).

(2.14)

Thus the density a posteriori is given (under the domination assumptions) as follows:

p(y|x) = p Y (y)p(x|y) R p(x|u)p Y (du) . (2.15)
This formalism has the advantage of making the density of the distribution of X|Y the one to be estimated. This approach has numerous benefits. It allows us to introduce propagation models as in Section 1.3.2. These propagation models are based on physical considerations and then provide relevant parametrization for the distributions of interest. Furthermore, assumptions (as the one in Section 3.2.1) are usually made and drastically decrease the algorithmic complexity of the parametric model inference.

We refer to the MAP estimator of Y the mapping Ŷ map defined as follows:

Ŷ map (x) = arg max z∈R p Y (z)p(x|z). (2.16)
One can also propose an estimator of the expectation of Y under the a posteriori distribution, noted Ŷ mmse :

Ŷ mmse (x) = R z p Y (z)p(x|z) R p(x|u)p Y (du) ν(dz).
(2.17)

Tree-Based Methods

Tree boosting is a highly effective and widely used machine learning method [START_REF] Jerome H Friedman | Stochastic gradient boosting[END_REF]; Maimon and Rokach, 2014). Among the machine learning methods used in practice, gradient tree boosting is one technique that shines in many applications. Recent well-optimized implementations of tree boosting (such as GBM or XGBoost) allow its incorporation into real-world production pipelines. This type of predictor is thus a relevant choice for practical application this thesis focuses on. We present hereafter the basics to define the eXtreme Gradient Boosting (XGBoost Chen and Guestrin ( 2016)) that will be essential for a good reading of this thesis, in particular the Chapters 3 and 4. The general framework is the regression problem in which an input random vector X ∈ X ⊂ R d is observed, and the goal is to predict a random response Y ∈ R by estimating the regression function m(x) = E Y |X = x . As usual, the objective is to use the data set Z n to built an estimate mn : X → R of the function m.

Regression Trees

The basics of a regression tree Decision tree can be both applied to classification and regression problem. A regression tree can be viewed as a tuple (q, ω) where q : R d → {1, . . . , L} is the structure of the tree that maps an input to a leaf index, ω ∈ R L is the vector containing the leaves weights and L is the number of leaves. For a observed input x, we use the decision rule given by q to map this input into a leaf and define the output as the weight of the leaf as illustrated in Figure 2.1.

We now briefly discuss the process of learning a regression tree. As mentioned above, a tree divide the feature space into non-overlapping regions (via its structure q). Each region is associated to a specific leaf. Let us denote these regions R 1 , . . . , R L . To every input that falls into a region, let say R j , the model prediction will be the mean of the target values for the training inputs in R j . Although these regions could have any shape in theory, it is preferred to give them a rectangular shape to ease the learning of the tree.

Recursive Binary Splitting

This greedy approach consists in successively splitting the predictor space from the top of the tree (the root) to the leaves. At each step, the algorithm chooses the direction of split j, and the cutpoint s such that dividing the predictor space into the regions {x | x j < s} and {x | x j ≥ s} leads to the greatest possible reduction of the regression error. That is, we consider all predictors directions 1, . . . , d, and all possible values of the cutpoint s for each of these directions, and then choose the direction and cutpoint such that the resulting tree has the lowest regression error. Next, we repeat this process, looking for the best direction and cutpoint in order to minimize the regression error within each of the resulting regions. The process continues until it reaches a stopping criterion; common stopping criteria are that the resulting regions contains a minimum number of observations, or that the tree contains a maximum number of leaves. input x

x 1 < 3 x 1 ≥ 3 x 3 ≤ 1.4 a 1 a 2 x 4 ≤ 7.1 a 4 a 3 Figure 2.1 -Generic
representation of a decision tree. The decision (output) of the tree given the input x is one of the weights (a i ) i=1,...,6 . For instance, the split at the top of the tree results in two large branches. The left-hand branch corresponds to x 1 < 3, and the right-hand branch corresponds to x 1 ≥ 3. The tree has four terminal nodes, or leaves. The number in each leaf is the mean of the response for the observations that fall there. As an illustration, the prediction of this tree given the input x = (7, 1, 3.1, 5, -1) is a 4 : x 1 = 7 ≥ 3, and x 4 = 5 ≤ 7.1.

Pros and Cons of Trees

Regression trees play a central role in this thesis and are at the center of the main contribution (see Chapter 4). It is therefore appropriate to recall here that they have a number of advantages over the more classical approaches developed in Chapter 3. They also suffer from a lack of robustness: a small change in the data can substantially modify the final estimator. To overcome this lack of robustness, we discuss in the next section the concept of boosting.

Trees can easily handle qualitative predictors without the need to create dummy variables (see the practical case herein).

Trees are highly interpretable and can be furthermore graphically displayed.

Trees can be highly non-robust, a small change in the data can cause a large change in the final estimated tree.

However, the aggregation of many decision trees allows to obtain a final estimator, called ensemble model, with a reduced sensitivity to a change of data [START_REF] Breiman | Stacked regressions[END_REF], and with better predictive performance at the price of losing interpretation. This aggregation procedure is discussed in the next section.

Boosting

There exists three methods for constructing ensembles of decision trees: bagging, boosting, and randomization. In this sequel, we mainly focus on the boosting method [START_REF] Freund | A short introduction to boosting[END_REF][START_REF] Freund | Experiments with a new boosting algorithm[END_REF]. This method is known to give better results as shown in [START_REF] Thomas | An experimental comparison of three methods for constructing ensembles of decision trees: Bagging, boosting, and randomization[END_REF] where the experiment section show that over a set of 33 tasks, the Adaboost algorithm gives the best results in most cases. However, this latter algorithm has already been proved to be sensitive to noise. Therefore, two of the most recent boosting algorithms (which perform better than Adaboost) are preferred hereinafter: eXtreme Gradient Boosting (XGBoost) (Chen and Guestrin, 2016), and Generalized Boosted Models (GBM) [START_REF] Ridgeway | Generalized boosted models: A guide to the gbm package[END_REF].

First, we define an ensemble tree model as a predictor consisting of a collection of T randomized regression trees. Let consider the t th tree, and denote by f t (x) the predicted value of this tree at the query point x. The trees are combined in an additive manner and thus form the forest estimate:

mT n (x) m n (x; f 1 , . . . , f T ) = T t=1 f t (x).
The idea of the boosting procedure is to learn the trees sequentially. Each new tree is learned with respect to the residual of the previous model rather that to the outcome Y as shown in Algorithm 2.2. This means that at each step a basis function that leads to the largest reduction of empirical risk is added into the estimator. In case of tree boosting, since a tree corresponds to split of the space of inputs into regions (see Section 2.3.1), the model is then improved in areas where it does not perform well. Note that in boosting, unlike in bagging, the construction of each tree strongly depends on the trees that have already been previously built. In general, boosting procedure can be tuned by means of parameters:

• The number of trees, noted T ∈ N * herein and corresponds to the parameter n_estimators in sklearn.ensemble.GradientBoostingRegressor class. Unlike bagging and random forests, boosting can overfit if T is chosen too big. We therefore use a validation set to select T .

• The shrinkage parameter ν > 0, a positive scalar. This parameter controls the rate at which boosting learns. At each stage, the contribution of each tree is shrunk by ν. This parameter is called learning_rate.

• The complexity of the tree, which can be controlled using different parameters such as the maximum depth of the trees, the number of leaves in a tree, or the minimum number of samples required to be at a leaf node denoted respectively by max_depth, max_leaf_nodes and min_samples_leaf.

In the literature [START_REF] James | An introduction to statistical learning[END_REF], it is often recommended to use small trees, with just a few terminal nodes. For this plot, we consider the well-known prediction problem in which we want to predict the housing values based on 13 features such as per capita crime rate by town (crim), the average number of rooms per dwelling (rm) or the nitrogen oxides concentration (nox). This dataset can be found in the sklearn.datasets1 .

In particular, simple stumps with a depth of one perform well if enough of them are included. This model can even outperform the depth-two model or a random forest. This highlights one difference between boosting and random forests: in boosting, because the growth of a particular tree takes into account the other trees that have already been grown, smaller trees are typically sufficient.

GBM

We begin this section on tree-based methods with a brief description of the recent Boosting implementation known as GBM. This implementation can be found in the gbm package (LightGBM). The gbm package takes the approach described in Friedman (2002). Some of the terminology somehow differs for generality. In addition, the gbm package implements boosting for models commonly used in statistics but not commonly associated with boosting. The gbm implementation of the boosting procedure is described in Algorithm 2.1

Algorithm 2.1 Boosting as implemented in gbm() Result: mn = t f t m n = 0, t = 1, and r i = y i ∀i ∈ {1, . . . , n}, a loss function L. while Stopping criterion is False do

Compute the negative gradient as the response.

r i = ∂ ∂m(x i ) L y i , m(x i ) m(x i )= mn(xi)
Randomly select p × k cases from the dataset.

Fit a regression tree f t w.r.t. those randomly selected observations. mn ← mn + f t t ← t + 1

XGBoost Algorithm

We now discuss in details the XGBoost algorithm, and conclude this section with a presentation of Gradient tree boosting also known as Gradient Boosted Regression Tree (GBRT) and specifically a scalable end-to-end tree boosting system called eXtreme Gradient Boosting (XGBoost). Boosting for regression trees is described in Algorithm 2.2.

Algorithm 2.2 Pseudo-code of Boosting for Regression Trees

Result: mn = t f t mn = 0, t = 1, and r i = y i ∀i ∈ {1, . . . , n}. while Stopping criterion is False do fit a tree f t to the training data

(x i , r i ) i=1,...,n . mn ← mn + f t ∀i ≤ n : r i ← r i -f t (x i ) t ← t + 1
The main contribution of GBRT is to use a second-order approximation of the objective in order to ease its minimization. Let be a regression loss function, typically the L 2 -norm, and Ω a regularizer to be defined later, and consider the following regularized objective:

L( mn ) n i=1 mn (x i ), y i + s Ω(f s ). (2.18) Let m(t-1)
n be the obtained predictor at the (t -1) th iteration:

m(t-1) n (x) = t-1 s=1 f s (x). (2.19)
We need to add f t to minimize the following objective.

L (t) (f ) n i=1 m(t-1) n (x i ) + f (x i ), y i + Ω(f ). (2.20)
To ease the minimization of the objective, second-order approximation originated from [START_REF] Friedman | Additive logistic regression: a statistical view of boosting (with discussion and a rejoinder by the authors)[END_REF] can be used:

L (t) (f ) i≤n m(t-1) n (x i ), y i + g i f (x i ) + 1 2 h i f 2 (x i ) + Ω(f ), (2.21)
where

g i = ∂ m(t-1) n (x i ) (•, y i ), and
h i = ∂ 2 m(t-1) n (x i ) (•, y i )
are the first and second order gradient of the loss function w.r.t. the first coordinate. Since the term ( m(t-1) n (x i ), y i ) is a constant of f , we can remove it to obtain a simplified version of the objective as follows:

L(t) (f ) = i≤n g i f (x i ) + 1 2 h i f 2 (x i ) + Ω(f ), (2.22)
It is possible to give a closed-form expression of the leaves weights ω for a fixed structure q of the predictor f . To this aim, we need to reindex the sum in Equation (2.22). Let l be an index of a leaf, (q, ω) the tuple representing the tree regressor f and I l {i ≤ n : q(x i ) = l} the set of indices that falls into the leaf l. Now, on the set of indices I l , the tree regressor f is constant and equal to ω l . Consequently, noting L the size of ω we have

L(t) (f ) = l≤L i∈I l g i ω l + 1 2 h i ω 2 l + λ ω 2 2 + γ|L|, (2.23)
where we expanded the regularizer Ω(f ) = λ ω 2 2 + γL. For a fixed structure q (L is constant), the optimal values of the leaves weights ω l (q) are thus given for all l ≤ L by :

ω l (q) = I l g i I l h i + λ . (2.24)
We can compute the value of the objective at the point q, ω (q) and use it as a scoring function to measure the quality of a tree structure q. To lighten notation, we will simply write ω l instead of ω l (q). It comes directly that this optimal value is given by

L(t) (q) = - 1 2 l≤L I l g i 2 I l h i + λ + γL. (2.25)
Now that we have a score to evaluate the quality of a tree structure (like the impurity score for evaluating decision trees (Shalev-Shwartz and Ben-David, 2014)) we can apply the greedy procedure introduced in Section 2.3.1 by starting from the root and iteratively adds branches to the tree.

Outline of the thesis

We would like to recall the place that the elements seen in this chapter will have in the sequel. Section 1.3.2 focuses on the problem of channel-based location estimation techniques, by means of two well-known techniques: the time delay estimation, and the RSSI-based range estimation. Both of these approaches involve a parameter estimation as well as the notion of CRLB. Chapter 3 provides a review of methods of geolocation. It also proposes a original classification of these methods which are referred to as fingerprinting methods and likelihood-based methods. While the first family of methods is quite self-contained, the last one relies on the estimation of the likelihood by means of Nadaraya-Watson estimates described above. In Chapter 5, we propose to assess the conditional independence between variables. The proposed test deeply relies on the conditional distributions estimation; we besides provide a practical example based on the Nadaraya-Watson estimates. Finally, in Chapter 4 we propose to adapt the boosting methods seen in Section 2.3 to metric learning. The learned metric has the particularity to be optimal w.r.t. the regression task pursued.

RSSI-based Methods for Location Estimation

"Adam, où es-tu?"

Genèse 3:9

In this chapter, we investigate machine learning approaches addressing the problem of geolocation. First, we review some classical learning methods to build a radio map. These methods are split in two categories, which we refer to as likelihood-based methods and fingerprinting methods. Then, we provide a novel geolocation approach in each of these two categories. The first proposed technique relies on a semi-parametric Nadaraya-Watson (NW) estimator of the likelihood, followed by a maximum a posteriori (MAP) estimator of the object's position. The second technique consists in learning a proper metric on the dataset, constructed by means of a Gradient boosting regressor: a k-nearest neighbor algorithm is then used to estimate the position. The proposed methods are compared on two data sets originated from Sigfox network, and an indoor dataset performed in a three-story building. Experiments show the interest of the proposed methods, both in terms of location estimation performance, and ability to build radio maps. This chapter covers the following publications:

• Elgui, K., Bianchi, P., Portier

Introduction

In In such data, the observed signals can be very noisy, especially in urban environment (RSSI based methods are often assisted with accelerometers, gyroscopes or Bluetooth beacons to improve their accuracy Yiu et al. ( 2017)). It may also happen that, due to range limitation or network sensitivity, some messages are not detected by some BSs.

In the network of interest for instance, the lowest observed value of the RSSI is -166 dBm. Below this value, the signal is not collected by the BSs, either because of a lack of sensitivity, or because of an SNR that is too low. In certain cases however, it can also be characteristic of a signal emitted from too far away. Despite the fact that experience has shown the benefits of taking into account the information of non-reception, this information is still, rarely used in the literature. The methods proposed in this chapter, in contrast, take into consideration both the non-isotropy and the information of non-reception. The detailed contributions of this chapter are the following.

Related

Contributions.

• We build a relevant metric for RSSI based geolocation. In order to achieve this, we propose to learn a non parametric metric function by means of a Gradient Boosting algorithm. The proposed method benefits both from the simplicity and robustness of k-NN and the regression performance of XGBoost, and more specifically from its ability to deal with complex and high dimensional data such as RSSI vectors.

The idea is to learn the metric used by the k-NN explicitly to enhance a predictor estimation [START_REF] Bellet | A survey on metric learning for feature vectors and structured data[END_REF]. Herein, the metric is built to compare two RSSI's vectors, such that the k-NN regressor provides the most appropriate neighbors for the geolocation. The main idea is to learn this metric d such that for a couple of RSSI's vectors (x, x ), d(x, x ) is a relevant predictor of the Euclidean distance between the two emitters' locations y -y 2 . We propose to learn d as a sum of T regression trees. Those trees are obtained through XGBoost algorithm (see Section 2.3). The benefits w.r.t. a classic k-NN regressor are twofold. First, the estimate benefits from the information of reception/non reception of the signal at a BS. Second, it improves the model by optimizing the metric explicitly for the task of geolocation. This leads to higher accuracy of the model as demonstrated in Section 3.4.

• We propose a method that takes advantage of the Boolean variable representing the reception/non reception of the signal at BSs. Our proposed model relies on a specific likelihood model of the RSSI's given the object's position. The expression of the likelihood is based on a model assumption of Naive Bayes type: given the emitter's position, the coordinates of the RSSI vector are assumed independent. It allows to build an estimator whose implementation is practical. The distribution of a RSSI at a given BS given the location of the emitter, is set as a Gaussian distribution. The mean and the standard deviation of this distribution are obtained by a non-parametric Nadaraya-Watson estimator Section 2.2.2. The final location is obtained using a Maximum-A-Posteriori (MAP) estimator. As will be shown in Section 3.4.1, one of the assets of the method is that it performs well, even on small training data sets. Finally, as a generative model, it allows to compute useful statistical information about our estimate such as density level sets and confidence regions unlike standard machine learning methods such as the k-NN estimator described below.

• We discuss the computational complexity as well as the memory cost of the methods founds in the related literature.

• We provide detailed experiments using real data originated from the Sigfox network, and an Indoor Positioning Dataset from Zsolt Tóth (2016).

Outline The rest of the chapter is organized as follows. Section 3.2 investigates several popular geolocation techniques of the literature. Then, in Section 3.3, we introduce the two proposed predictors. Finally, Section 3.4 is devoted to the numerical experiments and discussions.

Machine Learning Methods for Geolocation

In this section, we discuss different off-the-shelf predictors which can be used to solve the geolocation task introduced above.

Likelihood-based Methods

We refer to as Likelihood-based methods the methods which learn from the dataset a likelihood model p(x|y) for the conditional probability of the RSSI vector X given the position Y . These methods are presented in the Section 2.1.2.

One first learns from the observed data Z n a mapping p(x|y) which represents the conditional pdf of X|Y . One way is to introduce a parametric likelihood model as presented in Section 2.1.2, such as the path-loss model discussed at the end of this paragraph (and previously introduced in Section 1.3.2). The parameters of the model are then learned from the dataset. Non-parametric methods can be used as well (see Section 3.3). One of the main advantages is that some prior hypotheses on the form of the likelihood p(x|y) can be, based on physical considerations, easily introduced. In the sequel, we will assume as in Kaemarungsi [START_REF] Vincenty | Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations[END_REF], the parameters P 0,k , ν k respectively represent the power in dBm at distance d 0 and ν k is the so-called path-loss exponent.

The parameter vector θ = (θ 1 , . . . , θ d ) is then estimated from the dataset Z n using a standard maximum likelihood approach as described in Chapter 2. Once the parameter θ has been estimated, the second phase consists in replacing in Equation (3.2) the unknown distributions p k s by their estimates and thus obtain our estimator. 2009). We briefly present some popular learning algorithms to build this mapping.

Fingerprinting Methods

k-Nearest Neighbours

This method is used in Yiu et al. (2017); [START_REF] Patwari | Location estimation in sensor networks[END_REF] in the context of outdoor geolocation. As discussed in Section 1.5, some messages may not be detected by some BS, in which case we just set the corresponding RSSI value to NaN. In this section, we endow the space of RSSI's vectors with the Euclidean distance. For this aim, we need to replace all the NaN values. We suggest, as in [START_REF] Yiu | Wireless rssi fingerprinting localization[END_REF], to replace all the NaN values either by the lowest RSSI amongst all observed RSSI (in this paper -200 dBm will be used), or by an arbitrary value (the value -110 dBm is used in Yiu et al. ( 2017)).

For every D-dimensional RSSI vector X, we let (x (1) , y (1) ), . . . , (x (n) , Y (n) ) be a reordering of the dataset Z n such that X -

x (1) ≤ • • • ≤ X -x (n) .
The unknown position Y is finally estimated as follows:

Ŷ = 1 k k i=1 y (i) , (3.3)
where the integer k is an hyperparameter (see Hall 2016), the authors propose to use a Multi Layer Perceptron (MLP) approach to reduce the uncertainty in an indoor location estimation system. In the present section, we suggest to use a two hidden layers perceptron (see Fig. 3.2). The first hidden layer is composed by n h = 250 nodes. Each of these node computes a single output from multiple real-valued inputs by forming a linear combination according to its input weights, and puts this output through the logistic sigmoid activation function. The second hidden layer is similar to the first, but has no activation function. For every node l ≤ 250:

h (1) l (X) = σ   d α l,d X d + β l   = σ α l , X + β l , (3.4)
where σ is the logistic sigmoid function σ(x) = e x /(1 + e x ). For every node l ≤ 250:

h (2) l h (1) (X) = p γ l , p h (1) p (X) + δ l = σ γ l , h (1) (X) + δ l , (3.5)
where h (1) (X) h

(1) 1 (X), . . . , h

(dn) O(d|Y|) Yes k-NN O(1) O(dn) O(dnk) No Ensemble Trees O(dT n log n) O(2 q T ) O(qT ) Yes Perceptron O(nn h ) O(n h ) O(n h ) No Table 3
.1 -Criteria of the four presented methods. We denoted by d the dimension of the RSSI vector X, n h the number of neurons in the perceptron model. At last, T stands for the number of trees of the ensemble trees regressor, and q the depth of a single tree.

Discussion

In this part, we discuss technical aspects of the methods described above. This discussion is structured around four criteria:

• The training complexity of the model training.

• The memory cost of the predictor.

• The prediction complexity of the model.

• The ability to take into account the non-reception.

In Table 3.1, we illustrate the different criteria for the learning methods introduced in Section 3.2. We denote n the size of the training set, d the dimension of the RSSI vector X. In this table, we also denote n h the number of neurons in the perceptron model. Finally, T stands for the number of trees of the ensemble trees regressor, and q the depth of a single tree.

The likelihood based method is one of the cheapest method in term of both training complexity and memory cost. Indeed, this method only requires to stack into memory the conditional marginals p 1 , . . . , p d . Thus, the memory space of these methods depends on the number of parameters. Non-parametric methods can be used as well (see Section 3.3) to estimate this conditional marginals, in which case we have a space complexity of O(dn). However, the prediction is expensive because it requires to compute the arg max of the likelihood over the set Y.

Finally, while ensemble trees methods and perceptron offer lower both memory cost and prediction complexity, they require a long training phase. The time complexity of a method training is yet crucial. If this complexity is too high, one may not be able to retrain the model whenever the network changes. In that sense, the k-NN and the likelihood based methods should be preferred.

Proposed Geolocation Methods

In this section, we present two new geolocation approaches. The first relies on a semiparametric Nadaraya-Watson (NW) estimator of the likelihood, followed by a maximum a posteriori (MAP)estimator of the object's position. The second technique consists in learning a proper metric, constructed by means of a Gradient boosting regressor: a k-nearest neighbor algorithm is then used to estimate the position.

Semi-Parametric Likelihood-Based Method

In

π k (z) P(X k = NaN|Y = y) = E 1 {NaN} (X k )|Y = y
and we constrain the model by assuming that, given Y and given that X k = NaN, X k follows a Gaussian distribution whose mean and variance are respectively denoted by m k (y) and σ 2 k (y):

m k (y) E(X k |Y = y, X k = NaN) , σ 2 k (z) Var(X k |Y = y, X k = NaN) . We denote by Φ(x; m, σ 2 ) = ( √ 2πσ) -1 exp -(x-m) 2 2σ 2
the normal density of mean m and variance σ 2 . We summarize our model as follows:

1. X 1 , . . . , X d are independent given Y .

2. For every k,

P(X k ∈ dx|Y ) = π k (Y )δ NaN (dx) + 1 -π k (Y ) Φ x; m k (Y ), σ 2 k (Y ) dx. (3.7)
Based on this model, the likelihood p(x|y) is fully determined by the mappings π k , m k and σ 2 k for all k = 1, . . . , d. The remaining task is to estimate these quantities using our dataset Z n . To this end, we propose to use a non-parametric approach, and to replace these mappings with their Nadaraya-Watson estimates introduced in Chapter 2.

The Nadaraya-Watson Estimator for Non-Parametric regression. Let K : Y → R + be a kernel, i.e., non-negative, symmetric function integrating to one, and let h > 0 be a scalar (the so-called bandwidth). Define K h (y) = h -1 K(h -1 y) for all y ∈ Y.

The Nadaraya-Watson estimates are respectively given for every k ≤ d by:

πk (y) n i=1 1 NaN (x i,k )K h y i -y n i=1 K h y i -y , (3.8) mk (y) D k (y) -1 n i=1 1 R (x i,k )x i,k K h y i -y , (3.9) σ2 k (y) D k (y) -1 n i=1 1 R (x i,k )(x i,k -m k (y)) 2 K h y i -y , (3.10)
where D k (y) where I X {k = 1, . . . , d : X k = NaN} stands for the set of the receiving BS's, and where pY (y) stands for an estimation of the prior on Y , which we suggest to estimate from the dataset Z n through the kernel density estimator:

n i=1 1 R (x i,k )K h y i -
pY (y) = n -1 n i=1 K h y i -y .
The training step of this algorithm is very efficient. Thus, the model can be re-fitted very quickly. This can be useful if a BS is removed or shifted for example. This method enables us, in a very simple manner, to take into account the Boolean variable representing the reception/non reception of the signal at BS's. In addition, as a generative model, it allows to compute useful statistical guaranties on the location estimation such as confidence level sets (see Section 3.4).

With this choice of estimates, we have for every k ≤ d:

πk (y) = arg min a∈R n i=1 1 NaN (x i,k ) -a 2 K h y i -y , (3.12) mk (y) = arg min a∈R n i=1 x i,k -a 2 1 R (x i,k )K h y i -y . (3.13)
Proof The result follows immediately from taking the value that cancels out the gradient.

Thus, these estimates are obtained by a local constant least squares approximation of the target. The locality is determined by a kernel K h that downweights all the Y i that are not close to y whereas a plays the role of a local constant to be fitted. Because of the use of the conditional probability distribution of the components of X given Y instead of the one of Y given X, the choice of an appropriate kernel is easy.

We finally discuss broadly the Cross-validation procedure for bandwidth selection.

The choice of kernel bandwidth h is a key element for a kernel regression. Therefore, we suggest a K-Fold Cross-Validation selection of bandwidth h. Hereinafter, we provide a criterion which consists in assessing the regression performance of the reception given the position Y . We randomly divide the set of observations into K groups of equal size. These groups are denoted by

{I k } K k=1 . We denote by M SE (k) (h) = 1 |I k | i∈I k 1 NaN (x i )-πh (y i ) 2
the mean squared error computed on the observations on I k , where πh (y i ) stands for the kernel regressor of 1 NaN (X)|Y trained on {1, . . . , n}\I k . Finally, we chose h as h = arg min h In that sense, a "good" metric is a mapping d for which the empirical risk

1 K k M SE (k) (h).

A new metric for k-NN

R n (d) n i=1 n j=1 d(x i , x j ) -d v (y i , y j )
2 is small. The main trick, is to search for a mapping d minimizing R n (d) within a relevant hypothesis class. We suggest to search for d under the form

d(x, x ) T t=1 f t (Φ(x, x )) ,
where f 1 , . . . , f T is a collection of T regression trees, and where Φ :

X d × X d → R d × R d is given by: Φ(x, x ) |x -x | 1 2 (x + x ) . (3.14)
The first part of Φ(x, x ) encodes the relative position of the RSSI vectors and the second part their absolute position, as opposed to the implicit mapping of the Euclidean distance which only encodes relative information. In practice, the minimization of ER n (d) w.r.t. f 1 , . . . , f T is intractable. An alternative is to use a Random Forest or an XGboost regressor, which separately optimizes the T regression trees. In practice, the learning stage is thus as follows:

• Compute the pairwise features Φ(x i , x j ) for all couples (i, j) in the dataset. • Use a regression tree ensemble method to predict the labels d v (Y i , Y j ) based on the features Φ(X i , X j ) .

Note that the obtained mapping d, though symmetric, is not mathematically speaking a metric. This point is however irrelevant regarding the application of interest. Given the obtained metric and given an observed RSSI vector X, the k-NN estimate of Y is computed as in Equation (3.3).

Numerical experiments

Performance Analysis of Location Estimation Methods

Datasets

To compare the performance of the different methods, we have used two datasets.

UCI Dataset An indoor dataset from Zsolt Tóth (2016), composed of 1540 measurements performed in a three-story building. The measurements are made at the 32 base stations set in the building, and their positions are known. The measurements were recorded by the same kind of Android devices in order to reduce the effect of the variety of the hardware. The recording was preformed at weekend to reduce the noise of the environment. Sigfox Dataset The second one, is a private dataset originated from the Sigfox network. It is composed of n = 1.5 • 10 6 observations. These measurements are collected at the 164 base stations (see Figure 3.8).

These datasets are randomly split in a training set (90%) and a test set (10%). The training set was used to perform cross-validation (each fold containing 10% of the training set) in order to find the optimal parameters of our algorithms.

Results and discussions

To compute the errors we employ the Vincenty distance between estimated and actual location. Figure 3.4 shows the cumulative distribution function (cdf) of the prediction error for all the presented methods. Our proposed metric learning based k-NN turns out to outperform the other methods of the paper. It is particularly interesting to notice how the learned metric can improve a k-NN. Between these two methods, the optimal choice for the number of neighbours k is different. This parameter is chosen according a cross-validation (see Figure 3.3). Figure 3.3 shows the empirical risk of the Euclidean k-NN, and the proposed enhanced k-NN with respect to k. The main conclusion of this figure is the optimal value of k: 200 (resp. 25) neighbors minimize the empirical risk of the Euclidean k-NN (resp. the proposed k-NN). However, as k increases, the Euclidean k-NN inability to incorporate position information causes its performance to degrade. Our method, by comparison, shows fewer signs of degradation as k increases. In Figure 3.5 we display the locations of the k nearest neighbours for these two methods. As expected, the standard deviation of the k nearest neighbours around the true position is much lower for our proposed metric learning based k-NN than for the k-NN. By contrast, the Log Loss model is neither relevant for the noisy urban dataset nor for the UCI indoor one. The main assumption made by this model is that the signal is propagating through free space. This assumption cannot be valid in both of these datasets. Hence, the poor performance of the log-loss model.

Finally, notice that the presented likelihood-based method's location errors are relatively stable when the number of training points decreases (see Figure 3.6). This point is emphasized in [START_REF] Yiu | Wireless rssi fingerprinting localization[END_REF] where the impact of the density of the RSSI radio map is discussed.

Heat Map estimation

A major benefit of the Semi-Parametric Likelihood-Based method is that density level sets can be computed easily. Thanks to the statistical framework, this method is able to evaluate the probability density of Y |X at all y ∈ Y which can lead us to build density level sets in which Y is most likely to lie given the observation of X. This is shown in Figure 3.7 where further information is provided on the uncertainty of the estimation.

PoI Prediction Case Study and Metric Learning Motivation

Geolocation is usually described as a regression problem when the goal is to estimate the location (latitude, longitude). The advantage of this approach lies in its generality as it allows to deal with all situations. Still, in various business use cases and industrial applications, device's movements are constrained/limited and they can only move within a set of places. In IoT, logistic is a major application that fits well with this description: devices can only move within a set of a priori unknown specific locations such as warehouses, logistic hubs, etc.

By clustering the devices' locations in the original data set, the location diversity, using a positioning solution such as GPS or WIFI, reduces to a finite set of values. In this context, a Point of Interest (PoI) represents a cluster of true measured locations. A PoI must contain the spatial location of the cluster (usually defined as its center) but it can also be enriched with other features including the spatial extension of the cluster and the density of points within the cluster for example.

Note that a PoI's location can be known, partially known or unknown depending on the application and customer. In our pre-treatment, PoIs are supposed to be unknown and must be discovered. Their detection is a well studied problem that can be tackled by means of unsupervised machine learning techniques such as clustering methods. Density-based spatial clustering of applications with noise (DBSCAN [START_REF] Birant | St-dbscan: An algorithm for clustering spatial-temporal data[END_REF]) has proven to be efficient for this task as PoIs are defined as areas with high device or message density.

Attention should be paid to the fact that the set of PoIs may vary over time. It is thus necessary to maintain an up-to-date list of such PoIs by performing clustering steadily and continuously. This is illustrated in Figure 3.8 where some raw Sigfox devices locations are displayed in a map (a) and are then reduced to a few PoIs (b) after a clustering step.

Once the PoIs have been identified, the next step consists in predicting the unknown PoI belonging of a device, using a measured RSSI-vector. To achieve this goal, a data set made of fully supervised examples will be used: this data set is a collection of observed RSSI-vectors for various devices, each vector being mapped to one of the previously discovered PoIs. In this section only, a "good" metric (see Figure 3.10) is a mapping σ that minimizes the leave-one-out error:

R n (σ)

n i=1 i =i σ(Φ(x i , x i )), 1 y i =y i (3.15)
where is a certain loss function, typically the 0 -1 loss function (y, y ) = 1 y =y and Φ(x, x ) (|x -x |, 1 2 (x + x )) . One may search for a mapping σ that minimizes within any relevant hypothesis class such as logit model, decision tree, ensemble methods, support vector machines (SVM). Note that the obtained mapping σ, though symmetric, is not mathematically speaking a metric (no triangle inequality). This point is however irrelevant regarding the application of interest. Given the learned σn and an observed vector X, the plug-in predictor of the target Y is as following:

Ŷ = arg max j=1,...,J i∈I j σn Φ(x i , X) , (3.16)
where we recall that I j = {i ≤ n : y i = j}.

Similarity performance analysis

In this section, different classifier models are compared for the choice of the similarity function: a logit model, an Euclidean distance based model, a random forest and a Gradient Boosting. The Euclidean distance based model computes a similarity value 1/(1 + ||x i -x j || 2 ) for all couples of RSSIs vector (x i , x j ) in the data set. To analyse their performance, a dataset composed of n = 50000 messages received on the Sigfox network was considered. Instead of forming the n(n -1) possible couples as in Equation (3.15), we propose, for the sake of a lower computational cost, to build n + m random couples of RSSIs vectors (x i , x j ) with their corresponding labels 1 y i =y j . In our experiment, we considered n + m = 1000000 couples of messages.

The obtained dataset was split into a train dataset of size n and a test dataset of size m. The training set is used to learn the similarity function σn . On the other hand, the test set is employed to assess the methods' performance using their ROC curves as shown in Figure 3.9. For our dataset, the gradient boosting model shows the best performance in predicting whether or not two RSSI vectors x and x belong to the same PoI. In contrast, the Euclidean distance between two vectors of RSSI x and x fails to correctly represent the underlying distribution of the data. Indeed, close x and x (w.r.t. the Euclidean distance) does not correspond to the equality of labels. To illustrate this point, we computed the similarities' distributions for couples (x i , x j ) given whether or not y i = y j in Figure 3.10. It shows that similarity provided by classical k-NN when x and x share the same class (target = 1) are very likely to be null in comparison to the ones given by the gradient boosting model. Note that this figure also highlights that all classifiers are much more confident when predicting target 0 than target 1. This result was expected as it is fairly easy for any method to predict that many couples x and x do not belong to the same POI since their respective RSSI signatures are very different (such as couples of RSSIs with no common received BSs).

An other interesting metric to compare the methods is the lift (Figure 3.11). The lift is a measure of the performance of a classifier as having an enhanced response (with respect to the population as a whole), measured against a random choice targeting model. A targeting model is good if the response within the target is much better than the average for the population as a whole. The lift is simply the ratio of these values: target response divided by average response. Figure 3.11 shows for instance that the 10% highest similarities provided by the Gradient Boosting classifier correspond to couples that are 9.2 times more likely to share the same PoI than an average couple. This experiment shows that the Gradient Boosting based similarity outperforms the other methods in competition.

Euclidean Logistic regression

In the sequel, the similarity function is thus learned as a sum of T classification trees through Gradient Boosting algorithm.

Method

PoI prediction performance analysis

In this section, we focus on the original task that consists in predicting the device's PoI. To do so, we consider three different approaches: a classical k-NN (with Euclidean distance), a metric learning based k-NN (with Gradient Boosting based similarity) and a random forest classifier that provides state of the art performance on large feature space classification problems. To compare the performance of the methods in competition to predict the POI affiliation of messages, we have at our disposal the Sigfox dataset which is split in two, a training set Z n , and a testing set Z m . Note that to avoid overfitting, it is important that Z n , and Z m are built using different devices: some devices are used to learn and some other devices are used to predict.

For this problem, the training set was used to either learn a similarity function σn (as presented above) or directly learn the classification model (random forest in our case). On the other hand, the test subset Z m was employed to evaluate the algorithms' performance.

First, the Exact Match Ratio of the predicted PoI (see Equation (3.16)), that can be written as follows

1 m m i=1 1[y i = Ŷ (x i )] (3.17) is considered.
The results for the different methods are summed up in Table 3.2. This result shows that the gradient boosting metric learning method got the highest classification score.

In addition, the Vincenty distance between the predicted POI and the "true" position (given by a GPS or Wifi apparatus) was also computed. In Figure 3.12, the cumulative distribution function of the spatial error is compared for the different predictive methods.

This result also demonstrates that our proposed local non linear metric learning method outperforms the other considered methods on our data set.

Plotting the cumulative distribution also allows for a deeper interpretation of the observed performance. Indeed, the curves' shapes can be divided in three parts as follows:

• A steep rise up to 150 -200 meters that correspond to messages that can be correctly mapped to PoIs, which can also be interpreted as a measure the expressiveness of the algorithms. Among them, messages that originate from isolated PoIs having specific RSSIs' signatures are common to the different classifiers.

• A plateau up to 700 meters followed by a small jump for all the methods. This jump mostly corresponds to errors when PoIs are extremely dense (distance between PoIs lower than 1 km). The presence of this jump show the limitation of the tested Figure 3.12 -PoI prediction performance: cumulative distribution function of the error measured as the distance between the predicted PoI location and the "true" device location given by GPS or Wifi.

algorithms.

• A slow convergence before reaching a cumulative probability of 1 after several kilometers.

This analysis shows that a gradient boosting metric learning coupled with a k-NN model enables to outperform a classical Euclidean metric k-NN as well as state of the art classification predictors.

Conclusion

In this chapter, we investigated machine learning approaches addressing the problem of geolocation. We presented most popular methods found in the literature. We proposed two new techniques: one based on a likelihood model of the RSSI vector X given Y and an other based on a enhanced k-NN regressor. To compare these methods, 1.5M observations were collected from the Sigfox network, as well as 1540 measurements from the UCI dataset. Results have shown that our metric learning method has the highest accuracy on both of the datasets, with a 90th-percentile of errors of 5 km on Sigfox dataset, and of 7 m on UCI dataset. As for the semi-parametric method, in addition to being very easy to fit, it goes beyond the simple estimation by providing heat maps and level sets, making it a suitable method for industrial applications. Moreover, experiences have shown that the two proposed methods are the most accurate when we provide a reduced training set. They are robust when only a small training set is available. Furthermore, results show that the k-NN with the learned metric reaches the highest accuracy in the particular context of localization within PoIs, which therefore demonstrates the interest of the proposed method. Metric learning methods therefore deserves a great deal of attention. An interesting topic would be to build the metric by directly minimizing the regression error of the estimator. This would differ from the present chapter in which we provide a supervised objective for the learned metric.

Extreme Gradient Boosting for Similarity Learning

"Si, marchant dans une forêt aléatoire, tu rencontres deux fois le même arbre, c'est que tu es perdu."

Most of the fingerprint methods presented in Chapter 3 assume that positions corresponding to similar fingerprints are close to each others. By comparing the fingerprints through a distance (or similarity) function, the estimate is thus defined as the position of the nearest fingerprint found in the database. The quality of the prediction is then highly related to the chosen distance on the RSSI space. The metric learning problem is therefore a fundamental issue to improve RSSI-based geolocation techniques. This chapter addresses this problem. However, the results are general and go beyond the framework of the geolocation. The metric learning problem has been shown to improve regression methods that rely on distances or similarities. There is a surge of interest for optimizing distance and similarity functions. However, in most prior works, no link is made between the learned similarity and the estimator performance. In this chapter, we propose to build the similarity by directly minimizing the regression error of an estimator, and thus obtain an ad-hoc learning objective. To minimize this objective, we propose a modified version of the eXtreme Gradient Boosting algorithm (XGBoost) presented in Section 2.3. Experiments show that our model outperforms other kernel regression models on several benchmark datasets including one that proposes to locate emitting devices given the RSSI.

Introduction

The need for an appropriate way to measure distance or similarity between data points is ubiquitous in machine learning. The problem of learning such similarity functions is difficult and has attracted a lot of interest for the past twenty years, in particular with the seminal work of [START_REF] Xing | Distance metric learning with application to clustering with side-information[END_REF]. This led to the emergence of the metric or similarity learning field, that aims to automatically learn a function from the data to assess the similarity or the distance between pairs. As stated in [START_REF] Kulis | Metric learning: A survey[END_REF], this practice has besides proven useful when used in conjunction with nearest-neighbors methods or with any other regression method that relies on a distance such as the ones presented in Section 2.2.2. This effect is all the more noteworthy as the dimension of the data increases, as shown for instance in the work of [START_REF] Chopra | Learning a similarity metric discriminatively, with application to face verification[END_REF] that demonstrates the benefits of metric learning for the task of face recognition. In addition, the quality of the geolocation is highly related to the chosen distance on the RSSI space.

The metric learning problem is therefore a fundamental issue to improve RSSI-based geolocation techniques.

Nevertheless, in most previous works, the objective used to learn a similarity is generally not linked to the accuracy of the final predictor. In this chapter, we address this problem by proposing a general method to learn the similarity by directly minimizing the regression error of the final predictor. Consequently, the regression error is optimized with respect to the similarity function. This approach then requires the minimization of a particular learning objective. In previous works, such as Weinberger and Tesauro (2007); Keller et al. ( 2006), the choice of similarities is restricted to Mahalanobis distance metrics, which can be represented by symmetric positive semi-definite matrices. This parametrization allows a simple minimization stage by means of gradient based optimizers such as delta-bar-delta [START_REF] Robert | Increased rates of convergence through learning rate adaptation[END_REF].

In contrast, we propose here to learn a non-parametric similarity. More specifically, the similarity is chosen as a sum of regression trees and is sequentially learned by means of a modified version of XGBoost. This approach is shown to be well adapted to minimize the objective function of interest and furthermore benefits from the well-known qualities of XGBoost. The advantages of our methods are many:

• XGBoost is very efficient and scales well. It has been shown to be one of the most powerful and scalable supervised method for handling high-dimensional data.

• In contrast to local metric methods discussed subsequently, the storage requirement of our method is independent of the size of the input data.

• Our learned metric is non-parametric. This means that it can handle specific problem such as the XOR example discussed below.

As far as we know, learning a similarity that is used to improve regression by means of tree boosting methods is new.

Related literature. Methods that learn a similarity have been widely studied in recent years e.g. in Kar

(x, x ) (x -x ) M(x -x )
, with the Positive Semi-Definite (PSD) matrix M as parameter. This (pseudo)-distance corresponds to computing the Euclidean distance after a linear transformation of the data. These methods differ by the choice of the objective or by the regularizer used to learn M. We can cite as instance [START_REF] Davis | Information-theoretic metric learning[END_REF] in which a Burg divergence is used as regularization for M; in [START_REF] Xing | Distance metric learning with application to clustering with side-information[END_REF], the learning objective is to minimize the distance between points sharing the same class subject to the constraint that dissimilar points are separated (i.e the distance between those points is maximized). However, learning a linear metric such as a Mahanalobis distance can be too restrictive due to the following problems: first, because distance properties such as the triangle inequality must be violated (see [START_REF] Xiong | Random forests for metric learning with implicit pairwise position dependence[END_REF]). Second, because they demonstrate their limitation in simple case such as the XOR example as in [START_REF] Kulis | Metric learning: A survey[END_REF]. Kernel methods are used to extend linear methods to the nonlinear case, and this is achieved by writing algorithms in terms of inner products and replacing these inner products with kernel functions. The resulting distance is then a "kernelized" distance between pairs of examples. Numerous works employ such a kernelization of linear methods, in particular in computer vision. For example, the pyramid match kernel [START_REF] Grauman | The pyramid match kernel: Efficient learning with sets of features[END_REF], or face identification [START_REF] Guillaumin | Is that you? metric learning approaches for face identification[END_REF]. Although some algorithms have been shown to be kernelizable, a new formulation of the problem must generally be derived, and sometimes leads to insurmountable implementation problems [START_REF] Kulis | Metric learning: A survey[END_REF].

The work of [START_REF] Xiong | Random forests for metric learning with implicit pairwise position dependence[END_REF] applies Random Forests to learn a non-parametric metric called Random Forest Distance (RFD) in order to improve a classifier. This work emphasizes the benefits of using a nonlinear similarity in the classification task pursued. [START_REF] Kedem | Nonlinear metric learning[END_REF] proposes the GB-LMNN algorithm that applies gradient-boosting to learn non-linear mappings directly in function space. Neural network-based approaches offer the flexibility of learning arbitrarily complex nonlinear mappings [START_REF] Chopra | Learning a similarity metric discriminatively, with application to face verification[END_REF]. However, they often demand high computational expense, not only in parameter fitting but also in model selection and hyper-parameter tuning. Recently, metric learning has also shown its benefits in feature selection [START_REF] Navot | Nearest neighbor based feature selection for regression and its application to neural activity[END_REF] or to improve classification performance [START_REF] Goldberger | Neighbourhood components analysis[END_REF]. In this chapter, we focus on the problem of learning a task-specific similarity that improves regression.

Learning a metric that is later used in regression has been introduced in Goldberger et al. (2005) (Neighborhood Component Analysis (NCA)). In this work, authors propose to learn a global Mahanalobis distance measure to be used in a k-NN classifier by optimizing the expected classification leave-one-out error. The main innovation of this work is to replace the actual leave-one-out classification error, a discontinuous function of the matrix parameter M by a differentiable cost function based on stochastic ("soft") neighbour assignments in the transformed space. This approach has been extended for regression in Weinberger and Tesauro (2007). In this work, in addition to performing regression, the authors show that their algorithm can also be viewed as an algorithm for dimensionality reduction. Note also [START_REF] Noh | Generative local metric learning for kernel regression[END_REF] that analyzed the effect of the metric on the asymptotic bias of the kernel estimator, and proposed a Mahanalobis distance to alleviate the bias of the estimator. The main differences between the abovementioned articles is that our approach makes use of a non-linear similarity instead of learning a Mahalanobis matrix that is used in kernel regression. The work of [START_REF] Xiong | Random forests for metric learning with implicit pairwise position dependence[END_REF] proposes to learn a non-parametric metric called Random Forest Distance (RFD) in order to improve a classifier. The idea of using tree-based methods to build a similarity function is a central element of our contribution.

Outline. The chapter is organized as follows: in Section 4.2 we present our learning formulation for kernel regression. Then, we derive our proposed XGBoost-inspired algorithm to learn the similarity. We show that despite a complex learning objective the tree gradient boosting offers a suitable framework to minimize it. Finally, in Section 4.4 we compare our estimator to state-of-the-art regression models on both synthetic and real datasets.

Tree Boosting for Metric Learning

In this section, we propose to extend the XGBoost algorithm, detailed at length in Section 2.3 for a regression purpose, to the problem of metric learning. First of all, we provide the similarity learning objective corresponding to the unsupervised setting introduced in the introduction. Given the objective, we propose to chose the similarity function as a sum of regression trees. Each tree is then learned sequentially. This new objective induces drastic changes in the boosting algorithm that are presented in this section.

Problem Setting

Let (X, Y ) ∈ R d × R p be a random vector. We consider Z N = {x j , y j } N j=1 , N = m + n i.i.d. realizations of (X, Y ). By splitting the dataset Z N in two, we obtain Z n = {x j , y j } n j=1 to learn the similarity and Z m = { xj , ȳj } N j=n+1 to learn the regressor. Let σ : X × X → R be a mapping that we suppose symmetric and non-negative. We consider the following estimator of m

(x) = E Y |X = x : mσ (x) m j=1 σ( xj , x)ȳ j , (4.1)
The objective we minimize with respect to σ is the following regularized regression error of the predictor mσ on Z n :

L(σ) = n i=1 y i , mσ (x i ) + Ω(σ), (4.2)
where is a regression loss function and Ω a regularizer on σ to be defined later on. This approach is original since it is driven by the minimization of a regressor error.

Learning the similarity with XGBoost

We propose to learn σ as a sum of T regression trees f s , s = 1, . . . , T . As the similarity takes a couple (x, x ) as input, it is therefore necessary to make adjustments in the various functions involved in the definition of σ:

• Let Φ be a feature transform of a couple (x, x ) (see Section 3.4 for a practical example).

• Let f s be a regression tree, which can be viewed as a tuple (q, ω) where q : R d × R d → {1, . . . , L} is the structure of the tree that maps an input to a leaf index, ω ∈ R L is the vector containing the leaves weights and L is the number of leaves. For a given tuple (x, x ), we first map it into a feature space via Φ that is the input of f s . Then, we use the decision rule given by q to map this input into a leaf and define the output of the tree as the weight.

• The output σ(x, x ) is finally the sum of the T weights. This procedure is illustrated in Figure 4.1.

We propose, following the idea of XGBoost, to build σ sequentially and we denote by σ(t) the similarity at the t th iteration. We have for any pair

(x, x ) ∈ R d × R d : σ(t) (x, x ) = t s=1 f s Φ(x, x ) .
Consequently, the predictor of m(x) at the t th iteration becomes: To lighten notation, we will simply write m(t) instead of mσ (t) so that:

mσ (t) (x) m j=1 σ(t) ( xj , x)ȳ j . (4.3) Observation x Φ (x,x ) . . . f 1 f 2 f T

Sum of the T weigths

m(t) (x) = m(t-1) (x) + m j=1 f t Φ( xj , x) ȳj . (4.4)
The function f t added at the iteration t is defined as the minimizer of the following objective:

L (t) (f ) = n i=1 y i , m(t-1) (x i ) + f (x i ) ȳ + Ω f , (4.5) where f (x) = f Φ( x1 , x) , . . . , f Φ( xm , x)
, ȳ = ȳ1 , . . . , ȳm , is now a differentiable regression loss function defined on R p × R p ; and for a tree f with L leaves and weights ω, Ω(f ) = 1 2 λ ω 2 2 + γL. Using a second-order approximation of (y i , •) at the point m(t-1) (x i ) it follows that:

L (t) (f ) n i=1 y i , m(t-1) (x i ) + g i f (x i ) ȳ + 1 2 ȳ f (x i )H (y i , •) f (x i ) ȳ + Ω f ,
where g i ∇ m(t-1) (y i , •) ∈ R p and H (y i , •) ∈ R p×p are respectively the Jacobian and the Hessian of the function y i , • at m(t-1) (x i ). The term y i , m(t-1) (x i ) being a constant of f , we can remove it and define a simplified version of the objective as follows:

L(t) (f ) = n i=1 g i f (x i ) ȳ + 1 2 ȳ f (x i )H (y i , •) f (x i ) ȳ + Ω f , (4.6)
Proposition 4.1. Denoting as in Section 2.3 (q, ω) the tuple such that f t ( x, x) = ω q(Φ( x,x)) . We suppose furthermore that the current number of leaves in the structure q is L. First denoting by

I l i ≤ n, j ≤ m : q Φ( xj , x i ) = l , (4.7)
the pairs of indices for which the couple ( xj , x i ) falls into the leaf l and

J l,l i ≤ n, j ≤ m, k ≤ m : q Φ( xj , x i ) = l and q Φ( xk , x i ) = l (4.8)
the triplets (i, j, k) for which the pairs ( xj , x i ) and ( xk , x i ) fall respectively into the leaves l and l . Defining accordingly

G =    I 1 g i ȳj , . . . , I L g i ȳj    ∈ R L and (4.9) H =      J 1,1 ȳ j H i ȳk • • • J 1,L ȳ j H i ȳk . . . . . . . . . J L,1 ȳ j H i ȳk • • • J L,L ȳ j H i ȳk      ∈ R L×L . (4.10)
Then, we have the following:

L(t) (q, ω) = ω G + 1 2 ω H + λI L ω + γL. (4.11)
Proof Let us re-index the sum in Equation (4.6), such that the objective can easily be minimized with respect to the weights of the leaves ω of the tree f .

L(t) (f ) = L l=1 I l g i f (x i ) ȳj + 1 2 L l =1 J l,l ȳ j f (x i )H i f (x i ) ȳk + Ω f (4.12) = L l=1 I l g i ω l ȳj + 1 2 L l =1 J l,l ω l ω l ȳ j H i ȳk + 1 2 λ ω 2 2 + γL. (4.13)
Finally, with the definitions of Equation (4.9), and Equation (4.10), L(t) can then be rewritten in matrix form:

L(t) (q, ω) = ω G + 1 2 ω H + λI L ω + γL, (4.14)
.

The next lemma states that the so defined matrix H is positive semi-definite (PSD).

Lemma 4.2. The matrix H λ H + λI L is positive definite for any λ > 0. As a consequence, for a fixed structure q, the optimal values of the leaves weights ω are thus given by:

ω = -H + λI L -1 G. (4.15)

Proof

First, let us recall the existing relations between the sets that are considered. We assume that we are at the leaf and that the tree that is build contains L -1 leaves. This current leaf is denoted "parent" in the sequel. By definition, I parent is the set of pairs (i, j) ∈ {1, . . . , n} × {1, . . . , m} such that the couple (x i , x j ) falls into the "parent" leaf. We also have J parent, other that is the set of triplets (i, j, k) ∈ {1, . . . , n}×{1, . . . , m}×{1, . . . , m} such that (i, j) ∈ I parent and (i, k) ∈ I other . We suppose from now on that "parent" is split into "left child" and "right child". We therefore have the following relations: where the symbol stands for the disjoint union.

I parent = I left child I right child , ( 4
Let us denote by Y (l) ∈ R n×m×m the vector defined as follows:

Y (l) [i, j, k] =              ȳj if (i, j) ∈ I l ȳk if (i, k) ∈ I l √ ȳj ȳk if (i, j) ∈ I l and (i, k) ∈ I l 0 o.w. (4.19)
It follows that for any pairs of leaves indices (l, l ) we have:

H l,l = i≤n j≤m k≤m Y (l) [i, j, k]Y (l ) [i, k, j] = Y (l) , Y (l ) . (4.20)
The matrix in scrutiny is then a Gram matrix. It is therefore PSD. For any λ > 0, the matrix H λ is positive-definite. Let us note that the fact that the matrix Y is defined in terms of √ ȳj ȳk requires that the ȳj ≥ 0 for any j ≤ m, and ȳk ≥ 0 for any k ≤ m.

It follows from Proposition 4.1 and Lemma 4.2 that the optimal value of the objective L(t) for a given structure q has a simple expression that is provided in the following proposition.

Proposition 4.3. For a fixed structure q (L is constant), the optimal value of the objective function is given as follows:

L(t) (q, ω ) = - 1 2 G H + λI L -1 G + γL. (4.21)
Proof The optimal value of the objective is obtained by evaluative the value of L(t) (q, •) at point ω as defined in Equation (4.15). We plug (4.15) into (4.14) to obtain :

L(t) (q, ω ) = ω G + 1 2 ω H + λI L ω + γL (4.22) = -G H + λI L -1 G+ (4.23) 1 2 ( H + λI L -1 G) H + λI L H + λI L -1 I L G + γL (4.24) = - 1 2 G H + λI L -1 G + γL. (4.25)
This optimal value of the objective is used to measure the quality of the structure q. The procedure is the same employed in Section 2.3: we start from the root, and iteratively add branches to the tree by splitting the nodes until the maximal number of leaves is reached.

Remark Attention should be paid to the fact that, in contrast with the XGBoost algorithm presented in Section 2.3 the leaves weights (ω l ) l=1,...,L can not be computed separately. The value of a single ω l will influence the value of the other weights because of the correlation terms ω l ω l ȳ j H i ȳk . As a consequence, the computation can not be parallelized over nodes. The algorithm of the modified XGBoost is given in Section 4.3.

L 1 regularization on ω Consider the case when we add a L 1 regularization, that is Ω(f ) = λ ω 2 2 + µ ω 1 + γ|L|. This kind of regularization promotes sparsity so that most of the similarities between pairs of observations end up with a null weight. In this case, Equation (4.14) becomes:

L(t) (q, ω) = ω G + 1 2 ω H + λI L ω + µ ω 1 + γL. (4.26)
Denoting by X ∈ R L×L the unique PSD matrix verifying X X = H + λI L and ỹ ∈ R L a solution vector to X ỹ = G, a minimizer of Equation (4.26) is also a1 solution to:

ω ∈ arg min R L 1 2 ỹ -Xω 2 2 + µ ω 1 P(ω)
.

(4.27)

The optimization can be directly performed using solvers such as Celer [START_REF] Massias | Celer: a fast solver for the lasso with dual extrapolation[END_REF].

Implementation of the Metric Learning Algorithm

We now discuss in details the XGBoost metric algorithm and give the practical elements needed for its implementation. We provide as well the pseudo code and a discussion on the algorithm complexity. For simplicity of the presentation, we now assume that the loss function (y, y ) = y -y 2 2 .

Let us now suppose that we have built so far t -1 trees f 1 , . . . , f t-1 from the datasets Z n and Zm , so that at this point g i = 1 2 (y i -m (t-1) (x i )) and H i = 2 Id p for all i ≤ n; and that want to build the next one. Starting from the root, we have I root = {1, . . . , n} × {1, . . . , m}, and J root,root = {1, . . . , n} × {1, . . . , m} × {1, . . . , m}. The algorithm now goes through a greedy procedure to select jointly the best direction and value of a split.

For a given couple (direction, value of split), we obtain two children denoted subsequently by "left child" and "right child" and compute accordingly the sets I left , I right , J left,left , J left,right , J right,right , and the corresponding matrices G and H . Then, the score of the couple (direction, value of split) can be computed as in Equation (4.21). The couple with the higher reduction of the objective function is selected and the root leaf is split accordingly. The outcome of this procedure is then a stump with 2 new leaves. The latter procedure is repeated until a stopping criterion is reached as described in Algorithm 4.2.

Algorithm 4.1 Find best split of leaf l for t th tree.

Require: I l , {J l,l } l ∈leaves , score q score ← score q s ← [] for k ∈ {1, . . . , 2d}(Coordinates of Φ(x, x )) do for s ∈ Splits of Φ[k] do Compute I left , I right Compute G , H score ← score ∨ G H + λI L+1 -1
G end for end for return s with max score, score.

Practical cost of greedy procedure The main computation cost lies in finding the best split. In order to do so, a split finding algorithm enumerates over all the possible splits on all the dimensions of the feature space. When a split is proposed, two new leaves are added to a former leaf (a left and a right child), the sets I child and J l,child for l ∈ {1, . . . , L} and child ∈ {left, right} are built and the corresponding split score is computed. The cost of computing the score of a split is small, since the matrix is only of size of the current number of leaves, which is in practice lower that 8. Supposing that the tree is well-balanced, that is the number of observations is equi-distributed over the leaves. Each leaf contains consequently approximately n×m L pairs of observations. The complexity of choosing the best split is then O 2d × #splits × n × m × L 3 . The exact greedy algorithm is very powerful since it enumerates over all possible splitting points. However, it is impossible to efficiently do so when the data does not fit entirely into memory. To support effective gradient tree boosting, an approximate algorithm is needed. To summarize, the algorithm first proposes candidate splitting points according to percentiles of feature distribution. The algorithm then maps the continuous features into buckets split by these candidate points, aggregates the statistics and finds the best solution among proposals based on the aggregated statistics.

Algorithm 4.2 Pseudo-code: f t -builder Require: X n , Xm , , T, max N leaves , (g i , H i ) i=1...n q ← I Root f t ← 0 s q ← -1 I Root ← Z n × Zm J Root,Root ← Z n × Zm × Zm Current Leaves(f t ) := CL ← {Root} while |CL| ≤ max N leaves do for l ∈ CL do I l ← Compute according Equation (4.7)
J l,l ← Compute according Equation (4.8) split, g ← Find Best Split(I l , J l,l , s q ) end for NodeToSplit, split ← Node with higher g if s q < g then s q ← g split NodeToSplit CL ← CL + Children(NodeToSplit) end if end while return f t .

Numerical Experiments and Analysis

In this section we use synthetic and real data to compare our method (XGBoost Metric Learning, referred to as XGBML) against other methods referenced in the literature: a k-NN regressor (k-NN) using an Euclidean metric, a Nadaraya-Watson (NW) based on conventional kernels as in [START_REF] Hastie | The Elements of Statistical Learning[END_REF], and a Nadaraya-Watson with a learned local metric (NW+LLM, see [START_REF] Noh | Generative local metric learning for kernel regression[END_REF]).

Comparison with State-of-the-Art Metric Learning Methods

Results on synthetic data. In this section, we consider four synthetic datasets. First a classic swiss roll dataset SR in which we want to predict the radius. Then the two moons dataset TM in which we want to classify the point into two groups. The three blobs dataset TB used to classify point among each blob, and finally a half disk dataset HD for which the data are uniformly distributed over the unit half disk and we want to predict the angle. These four datasets are shown in Figure 4.2. For each dataset, the training sample size is n = 500, and the design sample size is m = 50. For the three experiments, our models achieves one of the lowest test errors. All results are provided in Table 4.1. Furthermore, the XGBoost Metric Learning provides interesting insights on the data. Figure 4.3 shows the similarities between a test point (indicated by a star) and the points of the design dataset. The width of the edges is increasing with the similarity between two points. We observe for instance that in Figure 4.3 that the highest similarities lie on the circle with radius of target-value. Furthermore, in Figure 4.4 we compute the Mean Squared Error (MSE) for different values of the ratio m/n. Our model achieves the lowest MSE for all ratios. Results on real data. In this section we consider two real world datasets. First, the classic housing dataset (Boston, available on scikit-learn2 ) in which we want to predict the house prices given 13 scalar features. Second, the Wireless Indoor Localization dataset (RSSI) which can be found on UCI website 3 . For both datasets, we proceed as follows. The dataset is split in 3: a design dataset Zm of size 50%, a training dataset Z n of size 40%, and a test dataset of size 10% used to compute the MSE. The results are shown in Table 4.2. Despite the higher dimensionality of these data, known for deteriorating kernel regressors performances, our method again exhibits the lowest test error for both datasets.

Benefits of Local Features on the Predictor Accuracy

We now highlight the benefits of the chosen pairwise position function Φ and in particular the improvements that the absolute position encoding have on the regression accuracy of our predictor. We build a predictor in which the features only encode for the relative information and compare it to our predictor. The models are compared The results are shown in Figure 4.5. Figure 4.5 shows that the addition of absolute position information drastically reduces the regression error of the estimator. Unsurprisingly, it also increases the impact of an additional tree on the error.

Conclusion and Perspectives

In this chapter, we proposed an original way to jointly learn a similarity function and the regressor based on this latter. Our approach consists in learning an similarity by directly minimizing the regression error. We further detailed at length an algorithm based on XGBoost to tackle this task. The proposed approach has been shown to be well adapted to minimize the objective function of interest and furthermore benefits from the wellknown qualities of XGBoost. Experiments demonstrated both the competitiveness and the generality of our approach, as well as the relevance of our learned metric. This work leaves exciting avenues for future works such as the study of our estimator consistency as instance. Analyzing the Rademacher complexity of the induced class of functions by invoking only Talagrand's lemma [START_REF] Ledoux | Probability in banach spaces[END_REF] results in a bound depending linearly on m. A more in-depth study is thus needed to refine this bound and would be the purpose of future work.

Conditional independence testing via weighted partial copulas

Taille et vocabulaire ne sont pas indépendants ; mais ils le sont conditionnellement à l'âge.

As stated in Chapter 3, conditional independence is closely related to the statistical model that we proposed for the couple of variable (X, Y ). Indeed, a direct consequence of the chosen model is that the coordinates of the variable X are independent given Y .

In this chapter, we introduce the weighted partial copula function for testing conditional independence. The proposed test procedure results from the following ingredients: (i) the test statistic is an explicit Cramér-von Mises transformation of the weighted partial copula, (ii) the regions of rejection are computed using a bootstrap procedure which mimics conditional independence by generating samples from the product measure of the estimated conditional marginals. Under conditional independence, the weak convergence of the weighted partial copula process is established when the marginals are estimated using a smoothed local linear estimator. Finally, an experimental section demonstrates that the proposed test has competitive power compared to recent stateof-the-art methods such as kernel-based test. This chapter covers the following publications:

• Elgui K., Bianchi P. & Portier F. Conditional independence testing via weighted partial copulas.arXiv e-prints, art. arXiv:2006.12839, June 2020.

Introduction

Let (Y 1 , Y 2 , X) be a triplet of real random variables. We say that Y 1 and Y 2 are conditionally independent given X if ∀(y 1 , y 2 , x) ∈ R 3 :

Pr(Y 1 ≤ y 1 , Y 2 ≤ y 2 | X = x) = Pr(Y 1 ≤ y 1 | X = x) Pr(Y 2 ≤ y 2 | X = x). (5.1)
This is denoted by Y 1 ⊥ ⊥ Y 2 | X and roughly speaking, it means that for a given value of X, the knowledge of Y 1 does not provide any further information on Y 2 (and vice versa). Determining conditional independence has become in the recent years a fundamental question in statistics and machine learning. For instance, it plays a key role in defining graphical models [START_REF] Koller | Probabilistic graphical models: principles and techniques[END_REF][START_REF] Francis | Learning graphical models with mercer kernels[END_REF]; see also [START_REF] Markowetz | Inferring cellular networks-a review[END_REF] for a study dedicated to cellular networks. Moreover the concept of conditional independence lies at the core of sufficient dimension reduction methods [START_REF] Li | Sufficient dimension reduction: Methods and applications with R[END_REF] and is useful to conduct variable selection in regression [START_REF] Lee | Variable selection via additive conditional independence[END_REF]. Finally, conditional independence is relevant in many application fields such as economy [START_REF] Huber | A test of the conditional independence assumption in sample selection models[END_REF] or psychometry [START_REF] Richard | Conditional independence in a clustered item test[END_REF]. In this chapter, we propose a new statistical tests to assess conditional independence.

The approach taken is related to the well-studied problem of (unconditional) independence testing, in which the most intuitive way to proceed is to compute a distance between the estimated joint distribution and the product of the estimated marginals [START_REF] Hoeffding | A non-parametric test of independence[END_REF]. Inspired by Kendall (1948), rank-based statistics have been extensively used in independence testing [START_REF] Hendrik | Asymptotic normality of nonparametric tests for independence[END_REF][START_REF] Ruschendorf | Asymptotic distributions of multivariate rank order statistics[END_REF][START_REF] Frits | Asymptotic normality of multivariate linear rank statistics in the non-iid case[END_REF]. Because rank-based statistics do not depend on the marginals, they have appeared as a key tool for modelling the joint distribution of random variables without being affected by their margins. This has led to the introduction of the copula function [START_REF] Deheuvels | An asymptotic decomposition for multivariate distribution-free tests of independence[END_REF], defined as the cumulative distribution function associated to the ranks. We refer to [START_REF] Fermanian | Weak convergence of empirical copula processes[END_REF]; [START_REF] Segers | Asymptotics of empirical copula processes under non-restrictive smoothness assumptions[END_REF] for recent studies on the estimation of the copula function. The copula function, which in principle measures the dependency between random variables, has been used with success in independence testing [START_REF] Genest | Test of independence and randomness based on the empirical copula process[END_REF][START_REF] Genest | Local efficiency of a cramér-von mises test of independence[END_REF]. Because the asymptotic distribution of the copula function is difficult to estimate, the related bootstrap estimate properties are of prime interest for inference [START_REF] Fermanian | Weak convergence of empirical copula processes[END_REF][START_REF] Rémillard | Testing for equality between two copulas[END_REF][START_REF] Bücher | A note on bootstrap approximations for the empirical copula process[END_REF].

The conditional copula of Y 1 and Y 2 given X is defined in the same way as the copula of Y 1 and Y 2 but uses the conditional distribution of Y 1 and Y 2 given X instead of the joint distribution of Y 1 and Y 2 . Compared to the copula, the conditional copula captures the conditional dependency between random variables and is thus useful to build conditional dependency measures [START_REF] Gijbels | Conditional copulas, association measures and their applications[END_REF]. Therefore, as in the case of independence testing, the conditional copula appears to be a relevant tool for building statistical test of conditional independence. This has been pointed out as a an "interesting open issue" in (Veraverbeke et al., 2011, Section 4).

In this work, a new statistical test procedure called the "weighted partial copula test" is investigated to assess conditional independence. The proposed approach follows from the use of an integrated criterion, the weighted partial copula, a function that equals 0 if and only if conditional independence holds. Given estimators of the conditional marginals of Y 1 and Y 2 given X, the empirical weighted partial copula is introduced to estimate the weighted partial copula and the test statistic results from a Cramér-von Mises transformation.

From a theoretical standpoint, the use of an "integrated" criterion enables to establish, in a general non-parametric framework, a convergence rate of order n -1/2 for the empirical weighted partial copula. More precisely, by using a smoothed local linear estimator for the conditional marginals, we obtain the weak convergence of the empirical weighted partial copula rescaled by n 1/2 . The rate of convergence n -1/2 , which is the same as the one derived in the (unconditional) independence test, is notable because conditional copula estimates are known [START_REF] Veraverbeke | Estimation of a conditional copula and association measures[END_REF] to converge at a slower rate, (nh d ) -1/2 where h is a smoothing parameter going to 0. Note finally that integrated criterion for testing has been frequently used in the conditional moment restrictions literature (see [START_REF] Lavergne | Smooth minimum distance estimation and testing with conditional estimating equations: uniform in bandwidth theory[END_REF] and the reference therein).

Inspired by the independence testing literature [START_REF] Beran | Nonparametric tests of independence between random vectors[END_REF][START_REF] Kojadinovic | Tests of independence among continuous random vectors based on cramér-von mises functionals of the empirical copula process[END_REF], the computation of the quantiles is made using a bootstrap procedure which generates bootstrap samples from the product of the marginal estimators to mimic the null hypothesis. Thanks to this bootstrap procedure, one is allowed to perform the weighted partial copula test using any marginal estimates as soon as one can generate from them.

Related literature. Testing for conditional independence has been considered only recently in the literature. Some of the existing approaches are based on comparing the (estimated) conditional distributions involved in the definition of conditional independence. The distributions can be compared using their conditional characteristic functions as in [START_REF] Su | A consistent characteristic function-based test for conditional independence[END_REF], their conditional densities as proposed in [START_REF] Su | A nonparametric hellinger metric test for conditional independence[END_REF], or their conditional copulas as studied in [START_REF] Bouezmarni | Nonparametric copula-based test for conditional independence with applications to granger causality[END_REF]. Unfortunately, the estimation of these conditional quantities are subjected to the well-known curse of dimensionality, i.e., the convergence rates are badly affected by the dimension of the conditioning variable. As a consequence, the power of the previous tests rapidly deteriorates if the conditioning variable has a large dimension. Note also that Bergsma (2010) uses partial copulas to derive the test statistic. Unfortunately, partial copulas fail to capture the whole conditional distribution and lead to detect a null hypothesis much larger than conditional independence.

Other approaches rely on the characterization of conditional independence using crosscovariance operators defined on reproducing kernel Hilbert spaces [START_REF] Fukumizu | Dimensionality reduction for supervised learning with reproducing kernel hilbert spaces[END_REF]. Extending the Hilbert-Schmidt independence criterion proposed in Gretton et al. 2017), the authors propose to train a classifier (e.g., XGBoost) to distinguish between two samples, one is the original sample, another one is a bootstrap sample generated in a way that reflects conditional independence. According to the accuracy of the trained classifier the test rejects, or not, conditional independence. This is further referred to as the classifier based conditional independence test (CCI-test).

Outline. In Section 5.2, we introduce the weighted partial copula test and provide implementation details including the mentioned bootstrap procedure. In Section 5.3, we state the main theorems (weak convergence results). In Section 5.4, the theory is illustrated by numerical experiments. Our approach is compared to the ones described in Zhang et al. (2012) when facing simulated datasets. Proofs are given in a supplementary material file.

The Weighted Partial Copula Test

Set-up and Definitions

Let f X,Y be the density function (with respect to the Lebesgue measure) of the random triplet

(X, Y ) = (X, Y 1 , Y 2 ) ∈ R d × R 2 . Let f X and S X = {x ∈ R : f X (x) > 0} denote
the density and the support of X, respectively. The conditional cumulative distribution function of Y given X = x is given by y → H(y | x) for x ∈ S X . The generalized inverse of a univariate distribution function F is defined as F -(u) = inf{y ∈ R : F (y) ≥ u}, for all u ∈ [0, 1], with the convention that inf ∅ = +∞. Since H( • |x) is a continuous bivariate cumulative distribution function, its copula is given by the function

C(u | x) = H F - 1 (u 1 |x), F - 2 (u 2 |x) | x , for u = (u 1 , u 2 ) ∈ [0, 1] 2 and x ∈ S X , where F 1 ( • |x) and F 2 ( • |x) are the margins of H( • |x).
We are interested in testing the null hypothesis that Y 1 and Y 2 are conditionally independent given X, that is,

H 0 : Y 1 ⊥ ⊥ Y 2 |X.
By definition (Dawid, 1979), H 0 is equivalent to H(y|x) = F 1 (y 1 |x)F 2 (y 2 |x), for every y ∈ R 2 and almost every x ∈ S X . Using the conditional copula introduced before, it follows that

H 0 ⇔ C(u | x) = u 1 u 2 ,
for every u ∈ [0, 1] 2 , and almost every x ∈ S X .

Let w : R d → R be a measurable function. The weighted partial copula is given by, for every u ∈ [0, 1] 2 and almost every t ∈ R,

W (u, t) = E (C(u | X) -u 1 u 2 )w(t -X) .
The proposed test follows from the observation, that H 0 is satisfied if and only if the function W is identically equal to 0 under a mild condition on w. This is presented in the following lemma whose proof is given in the supplementary material.

Lemma 5.1. Suppose that w : R d → R is integrable with respect to the Lebesgue measure and with a Fourier transform being non-zero almost everywhere, then H 0 is equivalent to W (u, t) = 0, for every u ∈ [0, 1] 2 and almost every t ∈ R.

The Test Statistic

In the following, we define a general estimator of W relying on some empirical copula construction that works for any estimate of the marginals F 1 and F 2 (see Section 5.2.5 for a typical example). That is, we first compute sample based observations of F j (Y j |X), j = 1, 2, by estimating each marginal F j . Those are usually called pseudo-observations. Second we define an estimate of W based on the ranks of the pseudo-observation. For the sake of generality, the estimator used for the conditional marginals is left unspecified in the subsequent development.

Let (X i , Y i1 , Y i2 ), for i ∈ {1, . . . , n}, be independent and identically distributed random vectors, with common distribution equal to that of (X, Y 1 , Y 2 ). We estimate the conditional margins in some way, producing random functions Fn,j ( • |x), j = 1, 2, and then we proceed with the pseudo-observations Ûij = Fn,j (Y ij |X i ). Let Ĝnj , for j ∈ {1, 2}, be the empirical distribution function of the pseudo-observations ( Û1j , . . . , Ûnj ), i.e. Ĝn,j (u) =

1 n n i=1 1 { Ûij ≤ u} , for u ∈ [0, 1].
From a conditioning argument, the weighted partial copula is given by

W (u, t) = E[(1 {F 1 (Y 1 |X) ≤u 1 } 1 {F 2 (Y 2 |X) ≤u 2 } -u 1 u 2 )w(t -X)].
(5.2)

The previous expression suggests the introduction of the following so-called empirical partial copula process, given by

Ŵn (u, t) = 1 n n i=1 1 { Ûi1 ≤ Ĝ- n,1 (u 1 )} 1 { Ûi2 ≤ Ĝ- n,2 (u 2 )} -u 1 u 2 w(t -X i ).
(5.

3)

The use of the transform G - n,1 and G - n,2 implies that Ŵn depends on Ûi1 and Ûi2 only through their ranks. Indeed, because Ĝn,j is a càd-làg function with jumps 1/n at each Ûij , it holds that Ûij ≤ Ĝn,j (u j ) is equivalent to ( Rij -1)/n < u j , where Rij = n Ĝn,j (U ij ) is the rank of U ij among the sample ( Û1j , . . . , Ûnj ). Hence, we have

Ŵn (u, t) = 1 n n i=1 1 {( Ri1 -1) < nu 1 } 1 {( Ri2 -1) < nu 2 } -u 1 u 2 w(t -X i ).
(5.4)

The test statistic is given by

Tn = [0,1] 2 ×R d
Ŵn (u, t) 2 dudt.

(5.5)

Remark 5.2. The test statistics Tn is of Cramér-von Mises type, as opposed to the Kolmogorov-Smirnov type (which would be defined taking the sup instead of integrating). In [START_REF] Genest | Test of independence and randomness based on the empirical copula process[END_REF] these two types of statistics are introduced in the context of (unconditional) independence testing. In our context, the Cramér-von Mises type is preferred over the Kolmogorov-Smirnov for practical reasons. Indeed, as shall be seen in the next section, a closed form exists for Tn .

Remark 5.3. The proposed estimate is of the same type as the copula estimator studied in Fermanian et al. ( 2004) and Portier and Segers (2018) but another definition might have been used. For instance, following the approach in Genest and Rémillard (2004), the statistic would be defined differently with

{ Rik ≤ (n+1)u k } in place of {( Rik -1) < nu k }, k = 1, 2.
Using one version or the other does not influence our results.

Computation of the Statistic

The following lemma provides a closed-form formula for the test statistics Tn .

Lemma 5.4. If w : R d → R is an integrable function, then

Tn = n -2 1≤i,j≤n M Ĝi , Ĝj w (X i -X j )
where Ĝi = ( Ri1 -1, Ri2 -1)/n, w = w w s with w s (x) = w(-x), and

M (u, v) = 1 -u 1 ∨ v 1 1 -u 2 ∨ v 2 - 1 4 1 -u 2 1 1 -u 2 2 + 1 4 1 -v 2 1 1 -v 2 2 + 1 9 .
Remark 5.5. The function w is left unspecified for the sake of generality. Examples include w(t) = exp(-t 2 ), w(t) = 1 |t|≤1 and other popular kernel functions such as the Epanechnikov kernel. In the simulations, we consider the Gaussian kernel as in this case, w remains Gaussian. In line with the result stated in Proposition 5.1, empirical evidences suggest that it does not have a leading role in the performance of the test. Another approach would have been to consider the (non-integrable) function w(t) = 1 {X≤t} as in [START_REF] Stute | Nonparametric model checks for regression[END_REF]. The same conclusion as in Lemma 5.1 remains valid in virtue of Proposition 16.10 in [START_REF] Billingsley | Probability and measure[END_REF].

Remark 5.6. Computing Tn requires n 2 operations. A sampling strategy might be to rather compute

1 |I × J| (i,j)∈(I×J) M ( Ri , Rj )w (X i -X j ),
with |I| = |J| denote random samples uniformly drawn in {1, . . . , n}.

Bootstrap Approximation

To compute the rejection regions of the test, we propose a bootstrap approach to generate new samples in a way that reflects the null hypothesis even when H 0 is not realized in the original sample. This has been notified as a guideline for bootstrap hypothesis testing in [START_REF] Hall | Two guidelines for bootstrap hypothesis testing[END_REF] and it enables, in practice, to control for the level of the test and to obtain a sufficiently large power.

The proposed bootstrap follows from the estimated conditional marginals of Y 1 |X and Y 2 |X, respectively Fn,1 and Fn,2 , and from the estimated distribution of X, denoted by Fn . First we choose X * uniformly over the (X i ) i=1,...,n , that is, X * ∼ Fn . Then we generate

Y 1 * ∼ Fn,1 (•|X * ), and
Y * 2 ∼ Fn,2 (•|X * ),
and execute the previous steps n times until obtaining an independent and identically distributed bootstrap sample of size n. We denote by (X * i , Y * i1 , Y * i2 ) i=1,...,n the obtained sample. We finally compute the test statistic based on this sample. We repeat this B times and obtain B realizations of the statistic under H 0 , denoted by (T * n,1 , . . . , T * n,B ). Now define the cumulative distribution function of the bootstrap statistics t → (1/B) B b=1 1 {T * n,1 ≤t} , and denote by ξ n (α) its quantile of level α ∈ (0, 1). The weighted partial copula test statistic with level α rejects H 0 as soon as Tn > ξ n (α).

A Generic Example using Nadaraya-Watson Estimator

In this section, the aim is to illustrate the proposed test procedure when using the classical Nadaraya-Watson estimator for the margins F j , j ∈ {1, 2} when d = 1.

Nadaraya-Watson estimator. Let K : R → [0, ∞) be the standard Gaussian density function on R. For x ∈ R and h > 0, put K h (x) = h -1 K(h -1 x). For j ∈ {1, 2}, the Nadaraya-Watson estimator of F j (•|x) is given by

Fn,j (y|x) = n i=1 1 {Y ij ≤y} K b j (x -X i ) n i=1 K b j (x -X i ) ,
(y ∈ R).

(5.6)

The choice of the bandwidths b 1 and b 2 is discussed below.

Post-nonlinear noise model. We consider Y 1 = cos(X + 1 ), Y 2 = cos(X +aY 1 + 2 ), where a ≥ 0 and X, 1 , 2 are independent standard Gaussian random variables with respective variances σ 2 X , σ 2 1 , σ 2 2 . When a = 0, H 0 is true, i.e., Y 1 ⊥ ⊥ Y 2 | X. On the other hand, H 0 is false when a > 0. The magnitude of a measures the distance to the null hypothesis. In the sequel we set σ 2 1 = σ 2 2 = 0.2 and σ 2 X = 1. Two illustrative samples are shown in Figure 5.1, one is drawn from the previous model with a = 0 and one other with a = 0.5.

Cross-validation selection of the bandwidth. The bandwidths b 1 and b 2 have a critical effect on the shape of the resulting estimates, and thus on the performance of our test. Indeed, these estimates of the margins Fj (y j |x) for j ∈ {1, 2} are used in the computation of the test statistic Tn as well as in the bootstrap procedure to simulate under the null (see Section 5. n (blue). The red area is the rejection region at 5%. Independently, Tn is computed 300 times to obtain the distribution of in gray. accordingly. We randomly divide the set of observations into K groups of nearly equal size. These groups are denoted by

{I k } k=1,...,K . Define MSE j,k (b) = (1/|I k |) i∈I k (Y ij - ĝ(-k) j,b (X i )) 2 , where ĝ(-k) j,b
stands for the Nadaraya-Watson estimate of the regression Y j |X computed on {1, . . . , n}\I k with bandwidth b. We choose b n,j as the minimizer of (1/K) K k=1 MSE j,k (b) over b. The success of the approach in distinguishing H 0 from its contrary is illustrated on Figure 5.2 considering the generic post-nonlinear noise model Section 5.4.3.

Remark 5.7. Though this example has been carried out using the Nadaraya-Watson estimate of the marginal distributions, other approaches to estimate the marginals can be used to conduct the weighted partial copula test. The only restriction on the employed marginal estimates comes from the bootstrap procedure in which the ability to generate according to the margins is required. For instance a k-nearest neighbours approach shall be considered in Section 5.4.

Weak Convergence

Smooth Estimator of the Margins

The theoretical results are provided in a general non-parametric setting, using a smoothed version of the local linear estimator [START_REF] Fan | Local polynomial modelling and its applications[END_REF] of the conditional marginals when d = 1 (see Remark 5.10 below). This estimate has been introduced in Portier and Segers (2018) and is a natural extension of the Nadaraya-Watson estimator (defined in the previous section). Such a non-parametric approach results in mild assumptions on the distribution of (X, Y 1 , Y 2 ). Let K : R → [0, ∞) and L : R → [0, ∞) be two kernel functions, i.e., non-negative, symmetric functions integrating to unity. Let (b n,j ) n≥1 and (h n,j ) n≥1 , for j = 1, 2, be four bandwidth sequences that tend to 0 as n → ∞. For (y, Y ) ∈ R 2 and h > 0, put

ϕ h (y, Y ) = y -∞ L h (t -Y ) dt.
(5.7)

with L h (y) = h -1 L(h -1 y). For j ∈ {1, 2}, we introduce the smoothed local linear estimator of F j (y j |x) defined by Fn,j (y j |x) = ân,j , (

where ân,j is the first component of the random pair (â n,j , bn,j ) = arg min

(a,b)∈R 2 n i=1 ϕ h n,j (y j , Y ij ) -a -b(X i -x) 2 K x -X i b n,j , (5.9)
where ϕ h in (5.7) serves to smooth the indicator function y → 1 {Y ≤y} . Kernels K and L do not have the same role: L is concerned with "smoothing" over Y 1 and Y 2 whereas K "localises" the variable X at x ∈ S X . For this reason, we purposefully use two different bandwidth sequences (b n,j ) n≥1 and (h n,j ) n≥1 . We shall see that the conditions on the bandwidth h n,j for the y-directions are weaker than the ones for the bandwidth b n,j for the x-direction. The assumptions related to the two kernels and bandwidth sequences are stated in ((Gc)) and ((Gd)) below. Note that if the previous optimization would be carried out only over a we would recover the Nadaraya-Watson estimator with a smoothed indicator function.

Weak Convergence of the Weigh Partial Copula

We rely on the following Hölder regularity class. Let δ ∈ (0, 1), k ∈ N, and M > 0 be scalars and let S ⊂ R be non-empty, open and convex. Let C k+δ,M (S) be the space of functions S → R that are k times differentiable and whose derivatives (including the zero-th derivative, that is, the function itself) are uniformly bounded by M and such that every mixed partial derivative of order l ≤ k, say f (l) , satisfies the Hölder condition

sup z =z f (l) (z) -f (l) (z) z -z δ ≤ M, (5.10)
where | • | in the denominator denotes the Euclidean norm. In particular, C 1,M (R) is the space of Lipschitz functions R → R bounded by M and with Lipschitz constant bounded by M .

We need the following four assumptions:

(Ga) The law P admits a density f X,Y on S X ×R 2 such that S X is a nonempty, bounded, open interval. For some M > 0 and δ > 0, the functions

F 1 ( • | • ) and F 2 ( • | • ) belong to C 3+δ,M (R × S X
) and the function f X belongs to C 2,M (S X ). There exists b > 0 such that f X (x) ≥ b for every x ∈ S X . For any j ∈ {1, 2} and any γ ∈ (0, 1/2), there exists b γ > 0 such that, for every y j ∈ [F - j (γ|x), F - j (1 -γ|x)] and every x ∈ S X , we have f j (y j |x) ≥ b γ .

(Gb) The function w : R → R is of bounded variation.

(Gc) The kernels K and L are bounded, non-negative, symmetric functions on R, supported on (-1, 1), and such that L(u) du = K(u) du = 1. The function L is continuously differentiable on R and its derivative is a bounded real function of bounded variation. The function K is twice continuously differentiable on R and its second-order derivative is a bounded real function of bounded variation.

(Gd) There exists α > 0 such that for any j = 1, 2, the bandwidth sequences b n,j > 0 and h n,j > 0 satisfy, as n → ∞,

nb 8 n,j → 0, nh 8 n,j → 0, b -1-α/2 n,j h 2 n,j → 0, nb 3+α n,j log b n,j → ∞, nb 1+α n,j h n,j log b n,j h n,j → ∞.
Let P denote the probability measure on the underlying probability space associated to the whole sequence (X i , Y i ) i=1,2,... . Let ∞ (T ) denote the space of bounded real functions on the set T , the space being equipped with the supremum distance. Define

U i1 = F 1 (Y i1 |X i ), U 2i = F 2 (Y i2 |X i )
, for any i = 1, . . . , n, and Wn (u, t) = Ẑn (u, t) -(f X w)(t)(u 1 Ẑn,2 (u 2 ) + u 2 Ẑn,1 (u 1 )),

(5.11)

for any u ∈ [0, 1] 2 , t ∈ R, with Ẑn (u, t) = n -1 n i=1 w(t -X i )(1 {U i1 ≤u 1 , U i2 ≤u 2 } -u 1 u 2 ) , (5.12)
and Ẑn,j (u j ) = n -1 n i=1 1 {U ij ≤u j } -u j . We now state our main result. Its proof is provided in the supplementary material.

Theorem 5.8. Assume that ((Ga)), ((Gb)), ((Gc)) and ((Gd)) hold. If H 0 holds, then for any γ ∈ (0, 1/2), we have when n → ∞ sup u∈[γ,1-γ] 2 , t∈R Ŵn (u, t) -Wn (u, t) = o P (n -1/2 ).

In addition, the process n 1/2 Ŵn (u, w)

u∈[γ,1-γ] 2 ,t∈R
converges weakly in ∞ ([γ, 1γ] 2 × R) to a certain Gaussian process.

Remark 5.9. Theorem 5.8 is a nontrivial extension of Theorem 2 in [START_REF] Portier | On the weak convergence of the empirical conditional copula under a simplifying assumption[END_REF]. By taking w t = 1 we would recover their result.

Remark 5.10. The approach employed follows from an approximation of the process W n by an oracle version of W n where the estimated marginals are replaced by the true ones.

In doing this, a crucial step consists in an embedding of some functions class involving estimated conditional quantiles into a Donsker class (van der Vaart and Wellner, 1996). First, because the estimated quantiles are difficult to control near the boundary of [0, 1], we need to restrict the proof to the interval [γ, 1 -γ] 2 . We believe that the extension to the whole interval is an interesting avenue for further research. Second, the regularity properties of the local linear estimate defined in (5.8) are essential to obtain that the resulting quantile functions are sufficiently smooth to be contained in a Donsker class. Third, as noticed in Portier and Segers (2018), the extension to higher dimensions, though feasible, is not straightforward and represents an interesting topic for future work. In the case d = 1, covered by Theorem 5.8, the rate of convergence, n -1/2 , is not affected by the size of the different bandwidths. We conjecture that this remains true in multiple dimensions with the same rate of order n -1/2 .

As a corollary of the previous weak convergence result, we obtain (invoking the continuous mapping theorem) the weak convergence, under H 0 , of a sightly modified version of Tn .

Corollary 5.11. Assume that ((Ga)), ((Gb)), ((Gc)) and ((Gd)) hold. If H 0 holds, then for any γ ∈ (0, 1/2) and any finite measure µ on R, we have that

n [γ,1-γ] 2 ×R Ŵn (u, t) 2 dudµ(t)
converges weakly to a tight non-negative random variable as n → ∞.

Numerical Experiments

In this section, we apply the proposed copula test on synthetic data to evaluate its performance based on the nominal level and the power of the test. We compare it with the KCI-test [START_REF] Zhang | Kernel-based conditional independence test and application in causal discovery[END_REF] presented in the related literature section. Since the level α is hard to set for the CCI-test of [START_REF] Sen | Model-powered conditional independence test[END_REF], this approach will only be considered when the proportions of correct decision will be computed (see Figure 5.6a).

In all the experiments, the function w is a Gaussian kernel given by w(t) = exp(-t 2 ) and the estimate of the marginals is the k-nearest neighbors version of (5.6) using the cross-validation approach of Section 5.2.5 to tune k.

We use four datasets, the linear model and the post-nonlinear noise model are considered. We also test our method on a probabilistic graphical model for testing causality detection. Finally, we apply our test in a practical setting, using the movie watching based brain development dataset [START_REF] Richardson | Development of the social brain from age three to twelve years[END_REF]. In all of our simulations we set α = 5% as the desired type-I error rate. Here, X ∈ R 2 , and we set

β 1 = β 2 = (1/ √ 2, 1 √ 2) T ∈ R 2 .
The shade of blue denote for the value of X T β.

Linear Model

Consider the joint distribution given by Y 1 = X T β 1 + 1 , Y 2 = X T β 2 + 2 , where X ∼ N (0, I d ), β 1 and β 2 are two constant vectors of [0, 1] d , and 1 , 2 are two standard Gaussian variables with Cov( 1 , 2 ) = a. When a = 0, H 0 is true. It is false when a > 0. We examine the effect of the constant a > 0, and the size of the dataset n on the type-II error rate. We also examine the type-I errors when the dimension of the variable X increases, in a setting where the null hypothesis H 0 holds. Two illustrative samples are shown in Figure 5.3, one is drawn from the previous model with a = 0 and one other with a = 0.5. We set

β 1 = β 2 = (1/ √ 2, 1 √ 2) T ∈ R 2 . Figure 5
.4 shows the attractive performance of our test compared to the KCI-test. Notably, we can see that in high dimensions, our test is more accurate with respect to the level set α than the KCI-test.

Causality Discovery

Hereinafter we consider a particular type of DAG called "Latent cause" model.

To draw samples from the alternative hypothesis, we break the conditional independence by adding an edge between the nodes Y 1 and Y 2 . The resulting graphs are shown in Figure 5.5 in dashed lines. For the "Latent cause" model of interest we have X ∼ N (0, 1), Y 1 |X ∼ N (X, 1), and Y 2 |X, Y 1 ∼ N (X + aY 1 , 1). It is easy to verify that H 0 is true when a = 0, and false otherwise. It can be seen in Figure 5.6 that for large sample size n, our test outperforms the ones in competition. Furthermore, our test is slightly more powerful than the KCI-test across a range of values of a but overall shows fairly similar performance.

Post-Nonlinear Noise

We apply the proposed test on the post-nonlinear noise model described in Section 5.2.5. We first examine the effect of the constant a > 0 on the probability of type-I and type-II error of our test. The results are shown in Figure 5.7. As expected, larger values of a yields lower type-II error probabilities. For every value a, we observe that the type-I error probability is closed to α. The performances of the tests are also compared when the sample size n changes. The role of n is critical and the results are shown in Figure 5.8. We note that the type-I error probability is again closed to α and that the type-II quickly vanishes when n increases. In this experiment, the proposed procedure outperforms the KCI-test.

Classification of Age Groups using Functional Connectivity

In this paragraph, we apply our test in a practical setting, using the movie watching based brain development dataset [START_REF] Richardson | Development of the social brain from age three to twelve years[END_REF] obtained from the OpenNeuro database1 . The dataset consists in 50 patients (10 adults and 40 children). The fMRI data consists in measuring the brain activity in 39 Region of Interest (ROI). For every patient, 168 measurements are provided for each ROI. We denote for j ∈ {1, . . . , 39} by X j the variable that represents the j th region signal value. Given two ROI j and j , we seek to test the null hypothesis that X j and X j are conditionally independent given X \{j,j } . The decisions given by our test allow us to obtain a map of connections between all the ROI, called connectome, given in Figure 5.9. In this figure, a line is drawn between two ROI whenever our test rejects the null for these two ROI. Here, due to the high dimension of the conditional variable, the margins are no longer estimated using a Gaussian kernel as in Section 5.2.5, but using a k-nearest neighbours approach.

L R L R

Mean Connectome over all adults Figure 5.10 -The blue circles represent the GPS coordinates of the BS that are used in this numerical experiment.

For a given x, the mapping Fn,j (y|x) is estimated for every y ∈ R as the proportion of samples i amongst the k-nearest nearest neighbours of x which satisfy Y ij ≤ y. The integer k is select by cross-validation.

As a sanity check, our connectome is used as an input feature of a classifier (Linear Support Vector Classifier (SVC)) in order to distinguish children from adults. We estimate the classification accuracy of our classifier using k-fold. The obtained accuracy is 97.4%. This result outperforms the standard correlation method (91.3%) and is close to the so-called tangent method (98.9%) which is known to be fitted for this task Dadi et al. (2019).

Copulas CI Test to Corroborate the Model of Equation (3.1)

The last numerical experiment is dedicated to corroborate the statistical model of Equation (3.1). Under this statistical model, the components X 1 , . . . , X d of the random vector X are independent conditionally to the position Y . In this section, we want to test that, for any couple of BS (i, j) ∈ {1, . . . , d} 2 , i = j, X i ⊥ ⊥ X j |Y. The datasets consists

Conclusion

In this chapter, we have developed test of conditional independence between Y 1 and Y 2 given X based on weighted partial copulas. First, under conditional independence, the weak convergence of the weighted partial copula process has been established under certain regularity conditions the marginals. We have shown that, empirically, our proposed test shows better performance, in terms of power, than recent state-of-the-art conditional independence tests such as the one based on a kernel embedding [START_REF] Zhang | Kernel-based conditional independence test and application in causal discovery[END_REF]. Furthermore, this test has been performed to corroborate one of the most commonly used assumption in RSSI-based geolocation literature, or more specifically, one consequence of this assumption. Numerical experiments have shown that, not only is this assumption not true for certain pairs of BS, but also that couples for whom the hypothesis was rejected share particular geographic dispositions with each others (geographic proximity or the fact that there are placed on a same road). An interesting topic would be to assess the improvement of prediction performance when these pairwise conditional independence are no longer assumed. Finally, the generality of the proposed approach makes it useful in many practical situations, and its soundness has been endorsed by conclusive numerical experiments.

Conclusion

In this thesis, we have proposed contributions to the RSSI-based geolocation problem in LPWAN.

First, we identified the limits of time-based approaches when applied to geolocation in the singular context of LPWAN; indeed, a preliminary investigation showed that the narrow band signals that are a requisite for IoT networks greatly jeopardize the performance of such geolocation methods. Then, we demonstrated that the well-known and widely used path-loss model is also not suited to address the task of geolocation because it fails to correctly represent the relation between the range to the emitter and the RSSI.

Then, we addressed this problem by means of recent techniques of machine learning, and presented most popular methods found in the literature. We furthermore proposed two new techniques. The first one relies on a semi-parametric Nadaraya-Watson (NW) estimator of the likelihood, followed by a MAP estimator of the object's position. The second one consists in learning a proper metric on the dataset, constructed by means of a Gradient boosting regressor. Contrary to prior works, the proposed methods take into consideration both the non-isotropy and the information of non-reception. This work emphasized how choosing a relevant distance on the RSSI space can greatly improve the quality of fingerprint methods estimation.

The quality of the prediction is indeed, highly related to the chosen distance on the RSSI space. The metric learning problem is therefore a fundamental issue to improve RSSIbased geolocation techniques. However, there are, as far as we know, few studies on the fingerprint problem from the perspective of metric learning. We therefore proposed an original way to jointly learn a similarity function and the regressor based on the latter. Our approach consists in learning a similarity by directly minimizing the regression error of the final predictor. We further detailed at length an algorithm based on XGBoost to tackle this task. The proposed approach has been shown to be perfectly adapted to minimize the objective function of interest. We then applied this method to popular regression problems such as the swiss roll dataset, and demonstrate that its range of applications go far beyond the field of geolocation.

Finally, independence conditional testing problem was tackled introducing the weighted partial copula function for testing conditional independence. The proposed test procedure is as follows: first, the test statistic is an explicit Cramér-von Mises transformation of the weighted partial copula. Second, the regions of rejection are computed using a boot-strap procedure which mimics conditional independence by generating samples from the product measure of the estimated conditional marginals. Furthermore, the weak convergence of the weighted partial copula process is established under the CI. Numerical experiments have shown that, not only is this assumption not true for certain pairs of BS, but also that couples for whom the hypothesis was rejected share particular geographic dispositions with each others (geographic proximity or the fact that there are placed on a same road). An interesting topic would be to assess the improvement of prediction performance when pair-wise conditional independence is no longer assumed.

([γ, 1 -γ] 2 × R).

A.2.2 Proof of Lemma A.1

Our proof is adapted from the proof of Theorem 1 in [START_REF] Portier | On the weak convergence of the empirical conditional copula under a simplifying assumption[END_REF]. For the sake of clarity, we start by recalling some of the results established in Portier and Segers (2018) that will be used further in our proof. Apart from this, the proof is self-consistent.

Fact 1. On a sequence of events whose probabilities tend to one, it holds that for every u j ∈ [γ, 1 -γ] and every (y j , x) ∈ R × S X , Fn,j (y j |x) ≤ u j ⇔ y j ≤ Fn,j (u j |x).

This is shown page 170 in [START_REF] Portier | On the weak convergence of the empirical conditional copula under a simplifying assumption[END_REF].

For u j ∈ [γ, 1 -γ], x ∈ S X , and j ∈ {1, 2}, define

∆n,j (u j |x) = F j F - n,j Ĝ- n,j (u j )| x | x - Ĝ(or)- n,j (u j ). (A.4)
Fact 2. We have for any j = 1, 2, sup

u j ∈[γ,1-γ] n 1/2 ∆n,j (u j |x) f X (x) dx = o P (1). (A.5)
This is shown page 171 in Portier and Segers (2018).

Fact 3. As established in (Portier and Segers, 2018, page 172), for each j = 1, 2,

x → F - n,j Ĝ- n,j (u j )|x : u j ∈ [γ, 1 -γ] ⊂ C 1+δ 1 ,M 1 (S X ), x → F - j Ĝ(or)- n,j (u j )|x : u j ∈ [γ, 1 -γ] ⊂ C 1+δ,M 2 (S X ),
with probability going to 1.

Therefore we apply Theorem

A.3 with ξ i = (X i , Y i ), X = S X × R 2 , V = [γ, 1 -γ] 2 × R
and E the space of measurable functions valued in R 4 and defined on S X × [γ, 1 -γ] 2 . Moreover, the quantities η 0 and ηn are given by, for every u ∈ [γ, 1 -γ] 2 and x ∈ S X ,

η 0 (u, x) = F -(u|x), F -(u|x) , ηn (u, x) = F - n Ĝ- n (u)|x , F -Ĝ(or)- n (u)|x . Identifying v ∈ V with (u, t) ∈ [γ, 1 -γ] 2 × R and η ∈ E with (η 1 , η 2 ), where η j , j ∈ {1, 2}, are valued in R 2 , the map m v,η : R 2 × S X → R is given by m v,η (y, x) = w t (x)(1 {y ≤ η 1 (u,x)} -1 {y ≤ η 2 (u,x)} ),
Finally, the space E 0 is the collection of those elements η

= (η 1 , η 2 ) in E such that {x → η 1 (u, x) : u ∈ [γ, 1 -γ] 2 } ⊂ C 1+δ 1 ,M 1 (S X ) 2 , {x → η 2 (u, x) : u ∈ [γ, 1 -γ] 2 } ⊂ C 1+δ,M 2 (S X ) 2 ,
where M 2 depends only on b γ and M . In the following we check each condition of Theorem A.3.

Verification of Condition (1) in Theorem A.3. Because the indicator function and w t are bounded, we have

|w t (x)| 1 {y ≤ F - n ( Ĝ- n (u)|x)} -1 {y ≤ F -( Ĝ(or)- n (u)|x)} 2 f X,Y (x, y) d(x, y) ≤ C w 2 j=1 sup u j ∈[γ,1-γ] 1 {y j ≤ F - n,j ( Ĝ- n,j (u j )|x)} -1 {y j ≤ F - j ( Ĝ(or)- n,j (u j )|x)} 2 f X,Y j (x, y j ) d(x, y j ),
so that we can focus on each margin separately. Recall that if the random variable U is uniformly distributed on (0, 1), then

E(1 {U ≤u 1 } -1 {U ≤u 2 } ) 2 = u 1 -u 2 . Writing ân,x (u j ) = F - n,j
Ĝn,j (u j )|x , we have

1 {y j ≤ ân,x(uj )} -1 {y j ≤ F - j ( Ĝ(or)- n,j (u j )|x)} 2 f X,Y j (x, y j ) d(x, y j ) = 1 {F j (y j |x) ≤ F j (ân,x(u j )|x)} -1 {F j (y j |x) ≤ Ĝ(or)- n,j (u j )} 2 f X,Y j (x, y j ) d(x, y j ) = S X F j (â n,x (u j )|x) - Ĝ(or)- n,j (u j ) f X (x) dx = S X ∆n,j (u j |x) f X (x) dx,
where ∆n,j has been defined in (A.4). Result 2 permits to conclude.

Verification of Condition (2) in Theorem A.3. This is given by Result 3.

Verification of Condition (3) in Theorem A.3. It is enough to show that the class of functions

(y, x) → w t (x)(1 {y ≤ g 1 (x)} -1 {y ≤ g 2 (x)} ) : t ∈ R, (g 1 , g 2 ) ∈ C 1+δ 1 ,M 1 (S X ) 2 × C 1+δ,M 2 (S X ) with R n,1 (u, t) = n -1/2 ( Ẑn ( Ĝ(or)- n (u), t) -Ẑn (u, t)), R n,2 (u) = ( Ĝ(or)- n,1 (u 1 ) -u 1 )( Ĝ(or)- n,2 (u 2 ) -u 2 ).
Now just recall the definition of Wn to obtain that

Ŵ (or) n (u, t) -Wn (u, t) = R n,1 (u, t) + P {w t }R n,2 (u) + P {w t } u 1 ρ n,2 (u 2 ) + u 2 ρ n,1 (u 1 ) , with ρ n,j (u j ) = ( Ĝ(or)- n,j (u j ) -u j ) + ( Ĝ(or) n,j (u j ) -u j ).
From Vervaat's Lemma (Segers, 2015, Lemma 4.3), we have that sup

u j ∈[0,1] |ρ n,j (u j )| = o P (n -1/2 ), sup u∈[0,1] 2 |R n,2 (u)| = O P (n -1 ).
Because the class of functions {(y, x) → w t (x)1 {F (y|x)≤u} : t ∈ R, u ∈ [0, 1] 2 } is Donsker (as demonstrated in Section A.2.1), the process Ẑn is asymptotically equicontinuous. This implies that

sup u∈[0,1] 2 |R n,1 (u)| = o P (n -1/2 ).
Consequently, each quantity in the above decomposition of Ŵ (or) n (u, t) -Wn (u, t) is o P (n -1/2 ), uniformly over u ∈ [0, 1] 2 and t ∈ R, and so comes the conclusion.

B

Résumé en français Ces dernières années, l'Internet des objets (IoT) a suscité une grande attention dans des domaines très divers tels que l'agriculture ou les soins de santé. Les experts s'accordent à dire que 30 milliards d'objets feront partie de l'IdO d'ici 2023 et que 40% de ces objets devront être géolocalisés, par exemple pour le transport de marchandises. L'un des défis les plus importants pour ce domaine est le besoin de localisation. En effet, de nombreuses applications des réseaux de capteurs doivent suivre des objets mobiles, tels que des personnes, des animaux, des voitures, etc. Pour que ces applications soient viables, le coût des dispositifs devra être faible (de quelques dollars à quelques centimes selon l'application) et les dispositifs devront durer des années, voire des décennies, sans remplacement de la batterie. En outre, le réseau devra s'organiser sans modération humaine importante. En outre, afin de permettre la connectivité de milliards de dispositifs, la plupart des réseaux dédiés à l'IdO utilisent des communications à longue portée et à faible puissance. L'IoT a suscité beaucoup d'intérêt ces dernières années. Il désigne les réseaux de dispositifs physiques dotés de capacités de communication. L'attente que tout soit connecté est à l'origine de cette tendance. On prévoit que, d'ici 2020, il y aura plus de 20 milliards d'objets communicants dans le monde (Hatton). Ces objets sont capables à la fois de collecter et de transférer des informations. Comme la maintenance fréquente des batteries doit être évitée en raison du nombre élevé de dispositifs prévus, une faible consommation d'énergie est également une exigence forte pour l'IdO. Par conséquent, les défis des réseaux IoT sont d'atteindre une grande évolutivité pour gérer un nombre massif de dispositifs, d'atteindre un faible coût et d'avoir une large couverture tout en gardant une faible consommation d'énergie. Les appareils qui répondent à ces exigences sont difficiles à intégrer dans les réseaux cellulaires traditionnels. C'est pourquoi une technologie dédiée LPWAN (Low Power Wide Area Network), telle que l'UNB (Ultra Narrow Band), développée et brevetée par Sigfox, est apparue. La connaissance de la géolocalisation de chaque appareil est une ressource très précieuse. En effet, elle permet à Sigfox de fournir cette information aux utilisateurs du réseau, ce qui débouche sur de nombreuses applications telles que la logistique ou le transport de marchandises, la surveillance et le suivi dans les bâtiments intelligents ou encore le marketing et la publicité de proximité dans les centres commerciaux. Cependant, les techniques de localisation traditionnelles telles que le système de positionnement global (GPS) ne sont donc pas bien adaptées aux besoins particuliers de l'industrie IoT. L'installation d'un GPS sur chaque appareil représente un coût et une consommation d'énergie prohibitifs pour de nombreuses applications, et n'est en outre pas adaptée aux applications intérieures. Toutes ces exigences compliquent grandement la localisation de ces objets.

Les méthodes alternatives basées sur la distance utilisent des mesures telles que les deux présentées ici, le temps d'arrivée (TOA) ou la différence de temps d'arrivée (TDoA) [START_REF] Ho | Solution and performance analysis of geolocation by tdoa[END_REF][START_REF] Cong | Hybrid tdoa/aoa mobile user location for wideband cdma cellular systems[END_REF], et l'indicateur de force du signal reçu (RSSI) pour estimer la distance entre un dispositif émetteur et une antenne réceptrice. D'autres méthodes de télémétrie sont couramment utilisées dans la littérature, comme l'angle d'arrivée (AOA) [START_REF] Niculescu | Ad hoc positioning system (aps) using aoa[END_REF], ou la différence de fréquence d'arrivée (FDoA) (Amar and Weiss, 2008) mais leur étude dépasse le cadre de cette thèse.

Nous sommes, comme indiqué dans l'introduction, intéressés par des scénarios où l'objet suivi est équipé d'un dispositif de communication, mais pas nécessairement d'un dispositif de positionnement tel que le GPS.

Ce chapitre nous permet de mettre en perspective les enjeux de cette thèse. Il a pour but de permettre de mieux comprendre la singularité du problème de géolocalisation dans un réseau de capteurs tel que le réseau Sigfox.

Généralement, les mesures de distance et d'angle utilisées pour la localisation sont affectées par des erreurs variant dans le temps et des erreurs statiques, dépendantes de l'environnement. Les erreurs temporelles (dues, par exemple, au bruit additif et aux interférences) peuvent être réduites en faisant la moyenne de plusieurs mesures dans le temps. Les erreurs liées à l'environnement sont le résultat de la disposition physique des objets (par exemple, les bâtiments, les arbres et les meubles) dans l'environnement particulier dans lequel le réseau de capteurs fonctionne. Comme l'environnement est imprévisible, ces erreurs sont imprévisibles et doivent être modélisées comme aléatoires. Toutefois, dans un environnement particulier, les objets sont essentiellement stationnaires et, par conséquent, pour un réseau de capteurs essentiellement stationnaires, les erreurs liées à l'environnement seront largement constantes dans le temps. La majorité des applications des réseaux de capteurs sans fil impliquent des capteurs essentiellement stationnaires. Comme un certain délai est acceptable dans ces applications, chaque paire de capteurs effectuera plusieurs mesures dans le temps et fera la moyenne des résultats afin de réduire l'impact des erreurs variant dans le temps. Tout d'abord, nous abordons le problème général de la géolocalisation dans Appendix B.2. Pour prédire la localisation d'un émetteur, les approches de pointe basées sur le canal consistent soit à estimer le délai entre l'émetteur et le récepteur, à partir duquel il est possible de déduire la distance entre l'émetteur et le récepteur; soit à prédire directement cette dernière distance à partir de la décroissance de puissance observée entre la puissance reçue et la puissance d'émission. Une introduction au canal sans fil, en particulier les paramètres clés pour le modéliser, et une présentation des prédicteurs discutés sont proposées dans la Section 1.3.

B.2 L'estimation de la localisation géographique

B.2.1 Principe de géolocalisation

Les termes géolocalisation et positionnement sont utilisés pour désigner l'estimation géographique dans le monde réel de l'emplacement d'un objet. Cette problématique est entrée dans notre société avec une diffusion massive. Il est utilisé dans de nombreuses applications telles que la navigation, la communication, les véhicules autopilotés, les objets connectés et les villes communicantes, ou plus récemment avec des problématiques telles que le contrôle de la contamination d'une population. Elle touche des domaines scientifiques très variés comme la Géolocalisation et Navigation par un Système de Satellites (GNSS), dont dépend de plus en plus notre économie. Les experts s'accordent à dire que 30% du produit intérieur brut dépendra en partie des GNSS d'ici 2030, contre 10% aujourd'hui.

La géolocalisation d'un objet fait référence aux coordonnées (latitude, longitude), c'està-dire à la position de l'objet sur la surface de la Terre. Parfois, le mot positionnement est plutôt employé lorsqu'il s'agit d'identifier l'emplacement d'un objet dans un espace particulier tel qu'un téléphone portable dans un centre commercial, ou un robot dans un bâtiment. Dans cette thèse, nous parlerons de géolocalisation lorsqu'il s'agit d'une localisation à l'échelle globale, alors que le terme positionnement est généralement utilisé dans des espaces intérieurs et/ou confinés.

Les méthodes radiofréquences sont utilisées pour la plupart des systèmes de géolocalisation. Cette famille de méthodes, également appelée méthodes de radiolocalisation, utilise les caractéristiques des ondes radio reçues pour prédire la localisation d'un objet émetteur. Les exemples sont nombreux. Le très répandu Global Positioning System (GPS) (voir Figure 1.1) est basé sur l'estimation du temps d'arrivée (TOA) d'un signal à un satellite. Lorsque l'heure de transmission, la vitesse de propagation et la position du satellite sont connues, le TOA conduit à un très bon estimateur de la distance entre l'objet émetteur et le satellite. La combinaison de plusieurs TOA conduit à l'estimation de la position de l'émetteur. L'utilisation de plusieurs récepteurs pour localiser un émetteur est connue sous le nom de multilatération. (illustré dans Figure 1.2). En téléphonie cellulaire, la radiolocalisation est effectuée directement par les stations de base (BS) du réseau cellulaire au moyen d'une ou plusieurs des caractéristiques suivantes :

• Le TOA (ou TDoA). Contrairement au GPS, ces quantités sont estimées par rapport à la station de base du réseau cellulaire.

• L'angle d'arrivée (AOA) correspond à la direction depuis laquelle le signal est reçu. Une façon pratique de déterminer l'AOA est de considérer que cette direction est celle où l'intensité du signal est maximale pendant une rotation complète de la station de base. En combinant plusieurs AOA, on obtient l'estimation de la position souhaitée.

• Le RSSI présente un grand intérêt dans cette thèse. Il correspond à la puissance du signal reçu moins la puissance du signal émis. Il donne un rendement à un système de télémétrie au moyen du modèle d'affaiblissement de chemin Log-distance décrit dans Section 1.3.2 ou par des méthodes basées sur les empreintes digitales (lorsque les différents lieux d'émission sont connus pour présenter des puissances très différentes signatures"). Ces méthodes sont étudiées en détail dans Chapter 4.

B.2.2 Géolocalisation basée sur le réseau

Cette thèse se concentre sur la géolocalisation basée sur le réseau. Ces méthodes utilisent uniquement l'infrastructure du réseau. Parmi toutes les méthodes de géolocalisation décrites, ces dernières sont les moins chères et nécessitent le moins d'énergie. Les méthodes de ce type ont rencontré un énorme succès avec l'apparition de l'Internet des objets (IoT) à la fin des années 1990. Essentiellement, le concept de l'IdO consiste à fournir à tout objet la capacité de transférer des données sur un réseau sans nécessiter d'interaction entre humains ou entre humains et ordinateurs. Aujourd'hui, l'ensemble des applications des dispositifs IoT est spectaculaire : maisons intelligentes (Samuel, 2016), applications médicales et de santé [START_REF] Catarinucci | An iot-aware architecture for smart healthcare systems[END_REF], agriculture (Mekala and Viswanathan, 2017) ou même systèmes de transport [START_REF] Zhou | Design and research of urban intelligent transportation system based on the internet of things[END_REF].

Nous passons ensuite à deux exemples concrets de géolocalisation rencontrés dans la pratique.

Logistique Le suivi des cargaisons et des biens via le réseau peut également présenter de grands avantages pour les systèmes de transport. On peut envoyer des alertes spécifiques lorsque des événements remarquables se produisent, comme l'arrivée dans un entrepôt.

Ecologie Il a ensuite été affirmé que l'IdO allait révolutionner le domaine de l'écologie. Tout d'abord, en termes de gestion de l'énergie : la connectivité d'un nombre important d'appareils consommateurs d'énergie (lampes, moteurs, pompes, etc.) peut leur permettre de communiquer avec les services publics non seulement pour équilibrer la production d'électricité mais aussi pour optimiser la consommation d'énergie dans son ensemble. Deuxièmement, les applications de surveillance environnementale de l'IdO utilisent des capteurs pour contribuer à la protection de l'environnement en surveillant, par exemple, la qualité de l'air ou de l'eau. D'autres applications, comme les systèmes de prévision des tremblements de terre ou des tsunamis, peuvent également être utilisées pour fournir une aide plus efficace.

B.3 Plan de la thèse et contributions

Dans cette thèse, nous proposons une autre classification des méthodes de géolocalisation basées sur le RSSI : les méthodes basées sur la vraisemblance et les méthodes basées sur les fingerprints. Les premières, qui englobent les méthodes "range based", consistent à apprendre (sur un jeu de données) un modèle pour le RSSI d'une station de base (notée X dans la suite) étant donné la position (notée Y dans la suite). Cette phase d'apprentissage est également appelée phase de calibrage dans la littérature. Une fois que ce modèle a été appris, on peut prédire la position de l'émetteur comme étant celle qui correspond au mieux au RSSI mesuré. Les secondes sont les méthodes basées sur les "fingerprints". Ces derniers font directement correspondre le RSSI à la position au moyen d'une fonction qui a été préalablement apprise sur un ensemble de données.

Les récentes avancées de l'apprentissage automatique et ses succès dans un large éventail de domaines ont poussé la communauté IoT à appliquer ces méthodes à la géolocalisation basée sur le RSSI. Présentons maintenant les contributions que nous avons développées pour répondre à cette problématique.

• Nous proposons d'abord des améliorations des méthodes utilisées pour le problème de géolocalisation basé sur le RSSI. La première technique proposée repose sur un estimateur semi-paramétrique de Nadaraya-Watson de la vraisemblance, suivi d'un estimateur du maximum a posteriori de la position de l'objet. La seconde technique consiste à apprendre une distance, construite au moyen d'un régresseur de type Gradient boosting : un algorithme de k-plus proches voisins est alors utilisé pour estimer la position. Les méthodes proposées sont comparées sur deux jeux de données provenant du réseau Sigfox, et sur un jeu de données intérieur provenant d'un bâtiment de trois étages. Les expériences démontrent l'intérêt des méthodes proposées, tant en termes de performance d'estimation de la position, que de capacité à établir des régions de confiance sur nos estimées. Les résultats montrent également que la qualité de la prédiction est fortement liée à la distance choisie sur l'espace RSSI. Le problème de l'apprentissage métrique est donc une question fondamentale pour améliorer la technique de géolocalisation basée sur le RSSI.

• Deuxièmement, nous introduisons un objectif original pour apprendre une similarité entre des paires de points de données. Dans ce manuscrit, nous proposons de construire la similarité en minimisant directement l'erreur de régression d'un estimateur. Nous obtenons ainsi un objectif d'apprentissage orienté vers la tâche. Pour le minimiser, la similarité est choisie comme une somme d'arbres de régression et est apprise séquentiellement au moyen d'une version modifiée de XGBoost détaillée dans ce document. Cette méthode bénéficie des qualités bien connues de XGBoost, telles que son efficacité et ses capacités de mise à l'échelle. De plus, notre similarité, bien que non-paramétrique, ne nécessite pas un stockage de la taille du jeu de données. Enfin, les expériences montrent que notre modèle surpasse les autres modèles de régression à noyau sur plusieurs jeux de données de référence.

• L'indépendance conditionnelle a été largement utilisée dans la littérature sur la géolocalisation basée sur le RSSI afin de réduire la complexité des modèles statistiques tels que ceux présentés dans ce manuscrit. Tester l'IC est donc essentiel pour la performance de tels estimateurs. Nous introduisons la fonction de copule partielle pondérée pour tester l'indépendance conditionnelle. La procédure de test proposée résulte des ingrédients suivants : (i) la statistique de test est une transformation explicite de Cramér-von Mises de la copule partielle pondérée, (ii) les régions de rejet sont calculées à l'aide d'une procédure bootstrap qui imite l'indépendance conditionnelle en générant des échantillons. Sous CI, la faible convergence du processus de la copule partielle pondérée est établie et confirme la solidité de notre approche. Les expériences démontrent enfin la compétitivité de notre approche par rapport aux méthodes récentes de l'état de l'art. First, we propose some methods for the RSSI-based geolocation problem. The observation is a vector of RSSI received at the various base stations. In particular, we introduce a semi-parametric Nadaraya-Watson estimator of the likelihood, followed by a maximum a posteriori estimator of the object's position. Experiments demonstrate the interest of the proposed method, both in terms of location estimation performance, and ability to build radio maps. An alternative approach is given by a k-nearest neighbors regressor which uses a suitable metric between RSSI vectors. Results also show that the quality of the prediction is highly related to the chosen metric. Therefore, we turn our attention to the metric learning problem. We intro-duce an original task-driven objective for learning a similarity between pairs of data points. The similarity is chosen as a sum of regression trees and is sequentially learned by means of a modified version of the so-called eXtreme Gradient Boosting algorithm (XG-Boost). The last part of the thesis is devoted to the introduction of a Conditional Independence (CI) hypothesis test. The motivation is related to the fact that for many estimators, the components of the RSSI vectors are assumed independent given the position. The contribution is however provided in a general statistical framework. We introduce the weighted partial copula function for testing conditional independence. The proposed test procedure results from the following ingredients : (i) the test statistic is an explicit Cram érvon Mises transformation of the weighted partial copula, (ii) the regions of rejection are computed using a boot-strap procedure which mimics conditional independence by generating samples. Under the null hypothesis, the weak convergence of the weighted partial copula process is established and endorses the soundness of our approach.

Institut Polytechnique de Paris 91120 Palaiseau, France

+ 0 0

 0 Set of integers from 1 to d included Y X Set of functions from X to Y R d 1 ×d 2 Set of real matrices of size d 1 by d 2 Id n Identity matrix in R n×nA i: i th row of matrix A A :j j th column of matrix A Tr A Trace of A ∈ R d×d Tr A = d i=1 A ii A Transpose of matrix A A † Moore-Penrose pseudo-inverse of matrix A • Euclidean norm on vectors and matrices S n ++ Positive definite matrices of size n × n S n Semipositive definite matrices of size n × n • S Mahalanobis matrix norm induced by S ∈ S n ++ A S = Tr(A SA) • 2 Spectral norm on matrices a ∨ b Maximum of real numbers a and b a ∧ b Minimum of real numbers a and b (a) + Positive part of a ∈ R a ∨ Vector

Figure 1

 1 Figure 1.1 -Principles of geolocation using GPS. This image comes from the Geolocation Wikipedia page.

Figure 1 . 2 -

 12 Figure 1.2 -Trilateration principle. Each circle represents all the possible locations of a mobile phone at a given distance (radius) of a cell tower. The aim of a trilateration algorithm is to calculate the (x, y) coordinates of the intersection point of the three circles.

Figure 1 . 3 -

 13 Figure 1.3 -Spherical coordinates (r, φ, ϕ): radial distance r, azimuthal angle φ, and polar angle ϕ.

Figure 1

 1 Figure 1.4 -Measured wideband path gain as a function of path length. Linear fit (-) with r -2 reduction.

Figure 1

 1 Figure 1.5 -Q-Q plot (w.r.t. normal distribution) of the residuals X -(P 0 + 10α log d d 0 ) shown for 3 BS in Paris.

Figure 1 . 6 -

 16 Figure1.6 -For each plots, we show the empirical mean of the RSSI at any point (r, ϑ), measured by the BSs placed at the center (0, 0). We notice that some directions (e.g. around 190 • for the left figure) show high values of RSSI and no loss of power in this directions. These behaviours correspond to major traffic routes.

Figure 1 . 7 -

 17 Figure 1.7 -Scatter plots of a sub-sample of the Sigfox training dataset. A red dot corresponds to an emitting position of a device in the Sigfox dataset.

  the considered IoT dedicated networks, with low power and bandwidth devices, localization is challenging. Standard methods such as famous range-based methods by means of TOA estimation as instance are not well suited for this application. Therefore, machine learning methods for localization have met a great deal of attention these past ten years. Besides they have shown their value for indoor scenarios as in Farjow et al. (2011); Wang et al. (2018). Nevertheless, the accuracy of such methods might be jeopardized in LPWAN[START_REF] Farid | Recent advances in wireless indoor localization techniques and system[END_REF]. In this chapter we focus on a baseline probabilistic RSSI-only localization algorithm. The main challenge comes from the large fluctuations of the observed RSSI values, for a given source location (see Fig.3.1).

Figure 3

 3 Figure 3.1 -Histograms from Sigfox network. Each plot shows histogram of the RSSI at a BS for three emitters of which the positions are fixed (the empirical standards deviations are resp. 13dBm, 7dBm and 13dBm). Given a fixed position of an emitter, the received RSSI's show large fluctuations.

Fingerprinting

  methods directly map the vector X into a location Ŷ (radio-map), typically by means of a supervised learning algorithm. There have been extensively studied Yiu et al. (2017); Brunato and Battiti (2005); Torres-Sospedra et al. (2015); Honkavirta et al. (

.

  Figure 3.2 -Figure of a two hidden layers perceptron

  the sequel, we propose a semi parametric likelihood-based method. Due to the NaNvalues, we modify the classical Gaussian Process model used Hähnel and Fox (2006) in order to leverage the information of non reception. As in Yiu et al. (2017); Hähnel and Fox (2006) we generate a likelihood model for signal strength measurements using a semi-parametric framework. We propose the following likelihood model for the conditional density p(x|y) of X given Y . Using the likelihood form presented in Section 3.2, it is sufficient to provide a model for the marginal conditional distributions p k (x k |y) of X k given Y , for every k = 1, . . . , d. Here, we recall that, unlike in Hähnel and Fox (2006); Yiu et al. (2017), X k is a random variable over the set R ∪ {NaN}. Densities are thus considered w.r.t. the reference measure λ + δ NaN where λ is the Lebesgue measure and δ NaN is the Dirac measure at the NaN-value. We define π k : Y → [0, 1] as

  y . Under standard technical conditions, πk , mk and σ2 k converge uniformly towards π k , m k and σ 2 k as n → ∞ and nh → ∞ (Tsybakov, 2013). Finally, the MAP location estimator can be written as: Ŷ = arg max y∈Y k∈I X (1 -πk (y)) log Φ(X k ; mk (y), σ2 k (y)) + k∈I c X (1 -πk (y)) + log pY (y) , (3.11)

  To enhanced the k-NN method, we can learn a appropriate metric Xie et al. (2016); Torres-Sospedra et al. (2015)) on X d . The impact of the distance on the prediction of k-NN regressor is widely discussed in Honkavirta et al. (2009), where several p-norms are compared. Nevertheless, there is no mention about learning such a distance. Let recall that NaN values are here replaced by a fixed real value as discussed in Section 3.2.2. Our contribution is to build a mapping d : X d × X d → [0, +∞) such that close RSSI (w.r.t. to the metric d) correspond to close object positions (w.r.t. to the Vincenty distance d v on Y).

Figure 3 . 3 -

 33 Figure 3.3 -Plots of RMSE (km) on Sigfox dataset vs. the number of neighbors k for the classic k-NN in (-), and the enhanced k-NN in (-). The RMSE is chosen as a criterion for the choice of k, as mentioned in Sec. 3.4.1

Figure 3 . 4 -

 34 Figure 3.4 -Comparisons of the presented methods in terms of their cdf of errors on the two datasets. left: Sigfox, right: UCI

Figure 3

 3 Figure 3.6 -The c.d.f. of errors on UCI dataset by the different presented methods assuming 90% (left) and 50% (right) of the database was available.

Figure 3 . 7 -

 37 Figure 3.7 -Heat map of the position Y |X = x for two different observations x. The axis are given in latitude and longitude scale. The red dots represent the true positions. The black dots represent the positions corresponding for the same observations of x in the test set.

Figure 3

 3 Figure 3.8 -a) Scatter plot of a sub-sample of the Sigfox training data set. A red bubble corresponds to an emitting position of a device in the Sigfox data set. b) PoIs obtained from the training data set displayed in a).

Figure 3

 3 Figure 3.9 -ROC curves measured for different classification model based metric.

Figure 3 Figure 3

 33 Figure 3.10 -Distribution of the similarity values predicted by the different classification models.In blue (resp. in green), all the couples (x i , x j ) such that y i = y j (resp. y i = y j ) are displayed.

Figure 4 . 1 -

 41 Figure 4.1 -Illustration of the similarity prediction σ (T ) (x, x ).

Figure 4 . 2 -

 42 Figure 4.2 -Scatter plots of the three synthetic datasets. The colormap stands for the target value y.

Figure 4 . 3 -

 43 Figure 4.3 -Graph representation (edges) of the similarities between the point of interest (star) and the points of the design dataset (dots).

Figure 4

 4 Figure 4.5 -MSE computed on test dataset with respect to the number of weak regressors. In black (resp. red) the model uses relative (resp. absolute and relative) information.

  (2008), Zhang et al. (2012) defines a kernel-based conditional independence test (KCItest) by estimating the cross-covariance operator. A surge of recent research (Doran et al., 2014; Runge, 2017; Sen et al., 2017) has focused on testing conditional independence using permutation-based tests. The seminal work of Candès et al. (2018) had led to many conditional independence tests depending on the availability of an approximation to the distribution of Y 1 |X, such as the conditional permutation test (CPT) proposed in ?. In Sen et al. (

Figure 5 Figure 5

 55 Figure 5.1 -On the left (resp. right) 2500 realizations of (X, Y 1 , Y 2 ) from the postnonlinear noise model when a = 0 (resp. a = 0,5). The shade of blue denotes for the value of X.

Figure 5

 5 Figure5.3 -On the left (resp. right) 10000 realizations of (X, Y 1 , Y 2 ) drawn from the model in 5.4.1 when a = 0 (resp. a = 0.5). Here, X ∈ R 2 , and we setβ 1 = β 2 = (1/ √ 2, 1 √ 2) T ∈ R 2 .The shade of blue denote for the value of X T β.

Figure 5 . 5 -

 55 Figure 5.4 -Simulation results for the linear model. Figures 5.4a, 5.4b show the probability of acceptance (i.e. the type II error rate), plotted against the constant a and n.Figure 5.4c shows the probability of rejection (type I error) against d. The plots show the average probabilities with standard error bars.

  Figure 5.6 -Simulation results for the causal inference datasets.

Figure 5 . 7 -

 57 Figure 5.7 -Simulation results for the post-nonlinear noise model. Figures 5.7a, 5.7b and 5.7c, show respectively the ROC AUC score, the probability of acceptance (i.e. the type II error rate), and the type I error rate plotted against the constant a with standard error bars.

Figure 5 . 9 -

 59 Figure 5.9 -Mean connectome provided by our test over all adults.

Figure 5 .

 5 Figure 5.11 -four couples of BS for which the conditional independence given the position Y does not hold.

  Tamás Zsolt Tóth. Miskolc iis hybrid ips: Dataset for hybrid indoor positioning. In 26st International Conference on Radioelektronika, pages 408-412. IEEE, 2016. pages 46, 53 Titre : Apprentissage statistique pour la g éolocalisation d'objets connect és Mots cl és : Apprentissage de Similarit é, G éolocalisation par le r éseau, Ind épendance Conditionnelle, M éthodes de Boosting R ésum é : La g éolocalisation par le r éseau a suscit é beaucoup d'attention ces derni ères ann ées. Dans un contexte o ù les signaux sont à bandes étroites, par exemple pour l'Internet des Objets, les techniques de g éolocalisation bas ées sur le RSSI se distinguent Nous proposons tout d'abord quelques m éthodes pour le probl ème de la g éolocalisation bas ée sur le RSSI. En particulier, nous introduisons un estimateur semi-param étrique de Nadaraya-Watson de la vraisemblance, suivi d'un estimateur de maximum à post ériori de la position de l'objet. Les exp ériences d émontrent l'int ér êt de la m éthode propos ée en termes de performance d'estimation. Une approche alternative est donn ée par une r égression de type k-NN qui utilise une m étrique appropri ée entre les vecteurs de RSSI. Nous nous int éressons donc au probl ème de l'apprentissage de similarit é et nous introduisons un objectif sp écifiquemet choisi pour am éliorer la g éolocalisation. La fonction de similarit é est choisie comme une somme d'arbres de r égression et est apprise s équentiellement au moyen d'une ver-sion modifi ée de l'algorithme eXtreme Gradient Boosting (XGBoost). La derni ère partie de la th èse est consacr ée à l'introduction d'un test d'hypoth èse d'ind épendance conditionnelle (IC). En effet, pour de nombreux estimateurs, les composantes des vecteurs RSSI sont suppos ées ind épendantes sachant la position. La contribution est cependant fournie dans un cadre statistique g én éral. Nous introduisons la fonction de copule partielle pond ér ée pour tester l'ind épendance conditionnelle. La proc édure de test propos ée r ésulte des él éments suivants : (i) la statistique de test est une transformation de Cram ér-von Mises de la copule partielle pond ér ée, (ii) les r égions de rejet sont calcul ées à l'aide d'une proc édure de "boot-strap" qui imite l'ind épendance conditionnelle en g én érant des échantillons. Sous l'hypoth èse nulle, la faible convergence du processus de la copule partielle pond ér ée est établie et confirme le bien-fond é de notre approche. Title : Contributions to RSSI-based Geolocation Keywords : Metric Learning, Network-Based Geolocation, Conditional Independence, Boosting Methods Abstract : The Network-Based Geolocation has raised a great deal of attention in the context of the Internet of Things. In many situations, connected objects with low-consumption should be geolocated without the use of GPS or GSM. Geolocation techniques based on the Received Signal Strength Indicator (RSSI) stands out, because other location techniques may fail in the context of urban environments and/or narrow band signals.

  However, the simulation results of[START_REF] Hassab | A quantitative study of optimum and sub-optimum filters in the generalized correlator[END_REF];[START_REF] Scarbrough | On the simulation of a class of time delay estimation algorithms[END_REF] show that ML estimator actual performance can be much worse for a given SNR and observation time. More specifically, when the bound is plotted as a function of signal-

	σ CRLB, Sigfox US 8.46 • 10 -5 s,	(1.7)
	σ CRLB, Sigfox Eu 5.07 • 10 -4 s	(1.8)
	Multiplying by speed of light gives the following ranging error variances:
	σ Sigfox US 25km,	(1.9)
	σ Sigfox Eu 150km.	(1.10)

). In practice, evaluating this expression on the values observed in the Sigfox network in the U.S and Europe (see also

[START_REF] Sallouha | Localization in long-range ultra narrow band iot networks using rssi[END_REF]

) leads to to-noise ratio one observes a distinct threshold. The work of

[START_REF] Chow | Delay estimation using narrow-band processes[END_REF]

states that below this threshold, the true bound can exceed the CRLB by large factors of the order of f s /B 2 .

  Let us also mention[START_REF] Jedari | Wi-fi based indoor location positioning employing random forest classifier[END_REF] which compares the performance of k-NN with random forest to predict the indoor location of using RSSI-based fingerprinting method. This work emphasizes the importance of the distance function used by k-NN. Besides, different distance functions such as Euclidean, Manhattan and Minkowski are compared in terms of accuracy of the resulting predictor. In[START_REF] Xie | An improved k-nearestneighbor indoor localization method based on spearman distance[END_REF], the authors exploit the Spearman rank correlation of RSSI measurements to improve k-NN. However, the authors of these prior works do not study the choice of the optimal metric. In contract, we propose to learn a (non parametric) metric function, by means of a Gradient Boosting algorithm. This metric is learned such that, when evaluated on a pair of RSSI vectors, it gives a good approximation of the Vincenty distance between the two emitters.Another approach is to model the signal strength as a function of the position. In the literature, these methods are generally referred to as range-based methods because they use a signal propagation model that maps an observed RSSI value to a distance (range) estimate. The localization procedure then tries to find the location which best agrees with the measured signal strengths. These methods then require the introduction of a likelihood model, which can be parametric or semi-parametric. We therefore refer to these approaches as likelihood-based methods. A commonly used parametric model is the so-called log distance path-loss model Wang et al. (2012); Xu et al. (2014). This model is a radio propagation model that provides a mapping between the received signal strength and the distance between the emitter and the receiver. The main advantage of parametric models is that they only require the calibration of few parameters, which helps the training phase of the model and do not require a large training set. However, they often fail to represent arbitrary distributions Yiu et al. (2017), which could deteriorate the geolocation accuracy (we refer to Zhou et al. (2018); Botteron (2003) for

literature. Machine learning techniques have been extensively considered for RSSI based localization, for instance to track customers

[START_REF] Oosterlinck | Bluetooth tracking of humans in an indoor environment: An application to shopping mall visits[END_REF]

, or in application to autonomous vehicles

[START_REF] Campbell | Sensor technology in autonomous vehicles : A review[END_REF]

. The first fingerprint-based method was proposed in

[START_REF] Bahl | Radar: An in-building rf-based user location and tracking system[END_REF]

. In this prior work, k-Nearest Neighbor (k-NN) is used to tackle the task of predicting the location given RSSI with RADAR, a system for locating and tracking users inside a building.

In recent years, neural networks have also been widely used for RSSI-based localization, in particular for indoor use cases

Iqbal et

al. (2018); Ahmadi and Bouallegue (2015); Zhang et al. (2016). For instance, Zhang et al. (2016) proposes to deal with the high variability of the RSSI by using a four-layer deep neural network structure to extract relevant features. in-depth studies of the geolocation accuracy of log-loss model). To remedy this, semiparametric models have been introduced. The latter benefit from the advantages of parametric models while still being able to represent arbitrary distributions. To build a likelihood in a semi-parametric fashion, many methods have been employed such as kernel methods Mirowski et al. (2012); Mahfouz et al. (2013, 2015) and Gaussian Processes (GPs) Schwaighofer et al. (2004); Hähnel and Fox (2006) or more recently in Bisio et al. (2017) with Smart 2 . First, as far as we know, in the related literature, few of works (Piórkowski and Grossglauser, 2006) have payed attention on the variable representing the reception/non reception of the signal at BSs discussed above. In most works, an arbitrary low RSSI value is attributed whenever a signal is not received as in Dashti et al. (2015); Janssen et al. (2020). Moreover, almost all of those found in the literature give a statistical relationship between the RSSI and the distance between emitter and receiver through a propagation model involving physical considerations. This implies that they therefore disregard the non-isotropy of the environment.

  , . . . , x d ) and where p 1 , . . . , p d are the conditional marginals to be learned from the training set. This implies that the components X 1 , . . . , X d of the random vector X are independent conditionally to Y .Assume now that a new message arises from the unknown position Y with an observed RSSI vector x. we can then make the choice to use the MAP estimator presented in Section 2.2.

		d
			p k (x k |y) ,	(3.1)
		k=1
	where x = (x 1 ŶMAP (X) arg max y∈Y	p(y|X)
		d
	= arg max y∈Y	k=1	log p

and Krishnamurthy (2004); Mazuelas et al. (2009); Li (2006); Bshara et al. (2010) that the conditional pdf admits the following expression: p(x|y) = k (X k |y) + log p Y (y).

(3.2)

The remaining task is to provide a model for p k 's, which can be learned from the training set. This step is commonly referred to as the calibration step and is described in the preliminary Chapter 2.

Let us now provide a practical example with the broadly used log-loss (or path-loss) parametric model

(Yiu et

al., 2017; Bshara et al., 2010; Barsocchi et al., 2009) described in Section 1.3.2. The conditional distribution p k (x k |y) of X k |Y is supposed to have the form p θ k (x k |y) where θ k = (P 0,k , ν k , σ 2 k ) is a triplet of parameters, p θ k ( . |y) is a Gaussian distribution of variance σ 2 k and mean P 0,k -10ν k log 10 d v (y, y k )/d 0 . Here, d 0 is some reference distance and d v stands for the Vincenty distance

  et al. (2008) for a discussion on the choice of k).

	Ensemble Trees Methods

Two important classes of ensemble methods are bagging methods such as random forests

[START_REF] Breiman | Random forests[END_REF]

, and boosting methods such as Gradient Tree boosting

[START_REF] Jerome H Friedman | Stochastic gradient boosting[END_REF]

. A Random Forest model has been applied as a classifier for a indoor-context geolocation

[START_REF] Jedari | Wi-fi based indoor location positioning employing random forest classifier[END_REF]

in which it gets better accuracy than a k-NN based method. In Li et al. (

2018

) used both TDOA and RSSI as input to indoor location estimation. The authors proposed a combination of a k-NN (to remove outliers) and a Random Forest Regressor.

Multi-Layers Perceptron

In

[START_REF] Ahmad | Modular multilayer perceptron for wlan based localization[END_REF]

;

Dai et al. (

  .16) J parent,other = J left child,other J right child,other if other = parent (4.17) J parent,parent = J left child,left child J right child,right child

	J left child,right child J right child,left child	o.w,	(4.18)

  Ŵn (u, t) -Wn (u, t) = o P (n -1/2 ). Lemma, the process { Ŵn (u, t)} u∈[γ,1-γ] 2 , t∈R and { Wn (u, t)} u∈[γ,1-γ] 2 , t∈R have the same weak limit in ∞ ([γ, 1 -γ] 2 × R). Now note that {x → w t (x) : t ∈ R} is a Euclidean or VC class with constant envelop C w = sup x∈R |w(x)| (Nolan and Pollard, 1987, Lemma 22, (ii)), i.e., the covering numbers are polynomials. Moreover, the class of indicator functions is also Euclidean (van der Vaart and Wellner, 1996, Example 2.5.4). This implies that both classes have finite entropy integrals and therefore are Donsker (van der Vaart and Wellner, 1996, Chapter 2.1, equation (2.1.7)). Using the preservation of the Donsker property through products and sums (van der Vaart and Wellner, 1996, Example 2.10.7 and 2.10.8), the class {(y, x) → w t (x)1 {F (y|x)≤u} : t ∈ R, u ∈ [0, 1] 2 } is Donsker. As a result, the process { Wn (u, t)} u∈[γ,1-γ] 2 , t∈R converges weakly to a tight Gaussian process in ∞

	Based on Lemma A.1 and A.2, we deduce that
	sup
	u∈[γ,1-γ] 2 , t∈R
	Invoking Slutsky's

  B.1 Contexte de la thèseCette thèse est le résultat d'une convention CIFRE (Convention Industrielle de Formation et de Recherche) entre Télécom Paris Saclay et Sigfox, un opérateur de télécommunications français créé en 2009 par Christophe Fourtet et Ludovic Le Moan. Sigfox est spécialisé dans le Machine to Machine (M2M) via des réseaux bas débit. Il contribue à l'Internet des objets (IoT) en permettant l'interconnexion via une passerelle. Sa technologie radio UNB ("Ultra narrow band") lui permet de construire un réseau cellulaire à bas débit et économe en énergie. Ce type de réseau est déployé dans les bandes de fréquences radio dites industrielles, scientifiques et médicales (ISM), disponibles dans le monde entier sans aucune licence.

  Une idée naturelle est d'appliquer les méthodes de radiolocalisation décrites ci-dessus au réseau Sigfox. Comme indiqué précédemment, les performances de ces méthodes dépendent fortement de l'infrastructure du réseau. Sans entrer dans les détails ici, nous proposons de donner quelques éléments qui compromettent leur utilisation et motivent donc les travaux ultérieurs. Tout d'abord, les stations de base Sigfox manquent de directivité et ne peuvent donc pas discriminer les directions des ondes incidentes. Par conséquent, les méthodes basées sur l'AOA ne sont pas pertinentes ici. De plus, Sigfox base ses communications sur la technologie dite Ultra-Narrow-Band (UNB) qui permet d'atteindre à la fois une longue portée et une autonomie prolongée. Ainsi, chaque signal possède une bande de fréquence d'une largeur de 100 Hz dans la bande de fréquence sans licence (allant de 868.0 à 868.6 MHz). Cette bande est populaire car elle présente un bon équilibre entre la portée, la pénétration des bâtiments et la possibilité d'utiliser de petites antennes. Néanmoins, il est bien connu que les performances des méthodes temporelles utilisant l'estimation du TOA ou du TDoA dépendent fortement de la largeur de bande du signal. La limite inférieure de Cramér-Rao (CRLB) sur les estimateurs TOA est donc souvent employée pour quantifier cet effet. Elle stipule que la variance minimale de tout estimateur sans biais du TOA est inversement proportionnelle à B 3 × SNR, où B représente la largeur de bande et Signal to Noise Ratio (SNR) est le rapport signal/bruit bien connu. Cette limite est donc défavorable à l'utilisation de telles approches et conduit par exemple à un écart-type des estimateurs de portée au moins égal à 20 km avec l'infrastructure Sigfox US. Néanmoins, elle constitue également une première valeur de référence à laquelle doivent être comparées les méthodes de géolocalisation que nous proposons dans la suite. On montre[START_REF] Boucher | Analysis of discrete implementation of generalized cross correlator[END_REF] qu'un autre élément défavorable est qu'elles sont particulièrement gourmandes en mémoire puisque leurs performances sont directement liées au choix de l'intervalle d'échantillonnage.B.2.3 Géolocalisation basée sur le RSSILe RSSI est une caractéristique pertinente pour la géolocalisation. Il est suffisamment explicatif sans constituer une charge mémoire. La littérature est extrêmement vaste et les méthodes sont souvent classée sous deux catégories: les méthodes "range-based" et les méthodes "free range". La première catégorie de méthodes utilise des mesures pour prédire la distance entre émetteurs et récepteurs (les coordonnées de la station de base sont connues dans ce cas). La combinaison de plusieurs portées estimées permet d'estimer la position de l'émetteur au moyen de certaines méthodes telles que la trilatération :[START_REF] Thomas | Revisiting trilateration for robot localization[END_REF]. Il y a donc toujours deux phases : une phase de télémétrie et une phase de localisation. En revanche, la deuxième famille de méthodes ne base pas ses prédictions sur l'estimation de la distance. Un exemple simple d'une telle méthode serait d'estimer la position comme le barycentre de la station de base réceptrice. Dans la mesure où le RSSI est la seule mesure dont nous disposons, les méthodes basées sur la distance montrent rapidement leurs limites[START_REF] Gayathri Chandrasekaran | Empirical evaluation of the limits on localization using signal strength[END_REF]. En effet, elles reposent essentiellement sur le modèle log-distance path loss dont le manque de réalisme est mis en évidence dans Section 1.3.2.

sklearn Boston Dataset

recall that the solution might not be unique

scikit-learn Boston Dataset

Accession number ds000228.

Remerciements

A

Proofs of Chapter 6

A.1 Proofs of the basic lemmas of Section 5.2

A.1.1 Proof of Lemma 5.1

The "only if" part is obvious. Suppose that the function W = 0 and let u ∈ [0, 1] 2 . Define g(x) = C(u | x) -u 1 u 2 )f X (x). We have (g w) = 0, a.e. on R d , where stands for the standard convolution product with respect to the Lebesgue measure. Applying the Fourier transform gives that F(g)F(w) = 0 which, by assumption, yields F(g) = 0. By the Fourier inversion theorem we obtain that g = 0 a.e. on R d . That is for any u ∈ [0, 1] 2 and any x ∈ S X , C(u | x) = u 1 u 2 .

A.1.2 Proof of Lemma 5.4

Write

Using the notation Ĝi , we have

First, let compute the first term of the right hand side. Let notice that the value of the integrand is 1 if u 1 > Ĝi1 ∨ Ĝj1 and u 2 > Ĝu2 ∨ Ĝj2 and 0 otherwise. Thus we obtain for this term:

Now let derive the second integral term of the right hand side, the third term will follow directly.

[0,1] 2

By combining (A.1) and (A.2) we obtain the desired result.

A.2 Proof of Theorem 5.8

We use notation from empirical process theory. Let P n = n -1 n i=1 δ (X i ,Y i ) denote the empirical measure. For a function f and a probability measure Q, write Qf = f dQ. The empirical process is

For any pair of cumulative distribution functions F 1 and

A.2.1 Sketch of the proof

We introduce an oracle copula estimator, defined as the empirical copula based on the unobservable random pairs

be its generalized inverse. The oracle estimator of W is then

with w t (x) = w(t -x). A crucial result is that the processes Ŵn and Ŵ (or) n are asymptotically equivalent as stated in the following lemma.

Lemma A.1. Assume that (G(Ga)), (G(Gb)), (G(Gc)) and (G(Gd)) hold. If H 0 holds, then for any γ ∈ (0, 1/2), we have when n → ∞,

Using the notation from empirical process theory, introduced below, we have

A second crucial result is the following one, where it is shown that Ŵ (or) n is asymptotically equivalent to Wn . Lemma A.2. Assume that (G(Ga)), (G(Gb)), (G(Gc)) and (G(Gd)) hold. If H 0 holds, we have when n → ∞,

We are based on Theorem 2.1 stated in (van der Vaart and Wellner, 2007) and reported below; for a proof see for instance van der Vaart and Wellner (1996, Lemma 3.3.5). Let ξ 1 , ξ 2 , . . . be independent and identically distributed random elements of a measurable space (X , A) and with common distribution equal to P . Let P denote the probability measure on the probability space on which the sequence ξ 1 , ξ 2 , . . . is defined. Let G ξ,n be the empirical process associated to the sample ξ 1 , . . . , ξ n . Let E and V be sets and let {m v,η : v ∈ V, η ∈ E} be a collection of real-valued, measurable functions on X .

Theorem A.3 (Theorem 2.1 in (van der Vaart and Wellner, 2007)). Let ηn be random elements in E. Suppose there exist η 0 ∈ E and E 0 ⊂ E such that the following three conditions hold:

1.

Then it holds that

The empirical process notation allows us to write

where w t (x) = w(t -x). To establish (A.3), we rely on the decomposition

Let γ ∈ (0, 1/2). The proof consists in showing that the empirical process term Ân,1 (u, t) goes to zero, uniformly over (u, t) ∈ [γ, 1 -γ] 2 × R, in probability (first step) and that the bias term Ân,2 (u, t) goes to zero, uniformly over (u, t) ∈ [γ, 1-γ] 2 ×R, in probability (second step). Assumption H 0 , will be crucial for treating the bias term in the second step.

First step: We show that

By Result 1, it holds that (with a slight abuse of notation) Under H 0 , we have, for every u ∈ [0, 1] 2 , t ∈ R,

Consequently, using the bound F 2 ≤ 1 and Ĝ(or)-

∆n,j (u|x) f X (x) dx

It remains to use Result 2 to obtain the conclusion.

A.2.3 Proof of Lemma A.2

Recall the definition of Ŵ (or) n (u, t) and Ẑn (u, t) that are given in Section A.2.1 and that under H 0 , in virtue of Lemma 5.1, P {w t (1 {F ≤u} -u 1 u 2 )} = 0. Notice that Ŵ (or) n (u, t) -Ẑn (u, t) = P n {w t (1 {F ≤ Ĝ(or)- n (u)} -1 {F ≤u} )} = n -1/2 G n {w t (1 {F ≤ Ĝ(or)- n (u)} -1 {F ≤u} )} + P {w t (1 {F ≤ Ĝ(or)- n (u)} -1 {F ≤u} )} = R n,1 (u, t) + Ĝ(or)n,1 (u 1 ) Ĝ(or)n,2 (u 2 ) -u 1 u 2 P w t = R n,1 (u, t) + R n,2 (u) + u 1 ( Ĝ(or)n,2 (u 2 ) -u 2 ) + u 2 ( Ĝ(or)n,1 (u 1 ) -u 1 ) P w t