Merci Laurie

Maxime, Clément, Ugo Manon Delphine

Pauline Thomas

Dianzhuo, Karim, Victorien, Moustapha, Mark, Chloé, Youssouph, Alexandrine, Maden, Wafaa Nathalie Aymeric

Elsa Rose

Julie Brice

Romain S D'epango

Yohann, Robin, Adeline, Ugo Merci Nazih

Clément, Sebastien Joseph Nicolas

Karzan Merci

BAYESIAN HIERARCHICAL MODELS FOR THE PREDICTION OF THE DRIVER FLOW AND PASSENGER WAITING TIMES IN A STOCHASTIC CARPOOLING SERVICE

Keywords: Hierarchical modelling, Gamma regression, GPS traces, MCMC, Multi-level moving average

REMERCIEMENTS

Mes premiers remerciements vont à mes deux directeurs de thèse Bertrand et Gérard. Vous avez été toujours disponibles et vous m'avez soutenu dès le début de ma thèse voire même avant son démarrage officiel. J'ai énormément appris grâce à nos échanges et vos encouragements ont toujours été précieux pour moi. J'ai eu l'immense honneur d'avoir été encadré par vous deux et j'espère que nos chemins vont se recroiser dans le futur.

Ensuite, je tiens à remercier Thomas et Arnaud les deux fondateurs d'Ecov, qui ont directement eu confiance en moi. Ce fut une expérience très instructive de travailler au sein de votre société. Merci de m'avoir donné cette chance.

Je remercie également Guillaume Saint-Pierre et Pierre Latouche pour m'avoir fait l'honneur de rapporter ma thèse. Merci pour vos remarques pertinentes

et Bénédicte pour tous les sujets administratifs entre Paris et Nantes.

Je remercie bien évidemment mes amis, votre soutien m'a toujours donné la force à accomplir ce travail. Merci Mous, Reda, Mounir, Andréas, Anne-Cécile, Selena, Laura, Alexis, Tassos, Aris, Dorothée, Jean-Christophe, Philippe A., Michel, Dimitri, Vago, Daphni, Nestor, Venia, Demba, Danie et Michel. Enfin, je remercie surtout ma famille. Merci Papa et Maman pour tout ce que vous avez fait et continuez à faire pour moi. Votre soutien et vos conseils ont toujours été les meilleurs. Merci aussi à mes deux grandes soeurs qui m'ont toujours supporté et conseillé de la meilleure manière. Vous avez toujours réussi à me réconforter pendant cette thèse (et même plus). Pour finir, je tiens à remercier ma conjointe Ioanna, merci de m'avoir soutenu pendant la thèse et bien avant. Tu as toujours su comment m'aider et je t'en remercie du fond de mon coeur.

CHAPTER 1

INTRODUCTION 1.1 Positionnement d'Ecov dans le marché du covoiturage

La pratique du covoiturage date des années 1950. En France, comme dans le monde entier, cette pratique s'est développée en parallèle avec l'internet. Ce n'est que plus récemment que les opérateurs de la mobilité ont commencé à collecter et exploiter les données sur la thématique du covoiturage. Un exemple important est le Registre National de Preuve de Covoiturage 1 qui centralise les données de plusieurs opérateurs de covoiturage en France. Tandis que le covoiturage de longues distances a largement fait ses preuves, avec BlaBlaCar (https://www.blablacar.fr/) comme acteur principal en France, le covoiturage du quotidien implique plusieurs acteurs et approches différentes.

La Figure 1.1 ci-dessous illustre les différentes catégories actuellement présentes en France. Dans ce sens, nous pouvons définir le covoiturage du quotidien comme un covoiturage correspondant à des trajets de courtes distances (moins de 200 km). Le passager informe les conducteurs potentiels de sa demande de covoiturage à l'aide de la console, demande qui est ensuite affichée sur le panneau électronique lumineux au bord de la voirie. Un conducteur peut récupérer le passager en toute sécurité sur la place de stationnement réservée à cette occasion. Photographie reproduite avec l'autorisation d'Ecov.

Introduction

Les arrêts connectés permettent aux passagers de faire leurs demandes de covoiturage via une console. Il est donc possible mais pas obligatoire d'utiliser un smartphone. Les conducteurs sont ensuite informés de cette demande grâce au panneau lumineux positionné sur la voirie et aussi via l'application mobile dédiée à cela. Cette approche est proche du système de bus et permet de faire le lien entre les conducteurs et les passagers.

Le réseau de covoiturage a ensuite évolué en proposant des lignes de covoiturage. Les passagers choisissent alors leur ligne de covoiturage déterminée par leur arrêts de départ et de destination. Ce changement de point de vue se reflète dans le Chapitre 2 et 3. Par exemple, la Figure 1.4 illustre le réseau de covoiturage du "Vexin", premier réseau implanté en Île-de-France, où les passagers pouvaient faire une demande de covoiturage sans se soucier de ligne de covoiturage avec arrêts déterminés.

Figure 1.4 -Carte du réseau de covoiturage du "Vexin". Les arrêts connectés où les passagers peuvent faire une demande de covoiturage sont indiqués par les cercles roses avec le logo du pouce.

La nouvelle structure, c'est à dire les lignes de covoiturage, est semblable aux lignes de bus des transports en commun. La Figure 1.5 représente le schéma Introduction des lignes de covoiturage qui sont développées dans la région de Lyon.

Figure 1.5 -Plan des lignes de covoiturage Lane. Schéma reproduit avec la permission d'Ecov.

Le passager utilisant une ligne de covoiturage a un choix précis de destinations parmi les arrêts desservis de cette ligne. Ainsi, cette dernière approche se démarque des concurrents, étant donné qu'il ne s'agit pas seulement d'une plateforme et/ou application mobile mais d'une solution complète combinant digital et infrastructure physique sur le territoire. Il est intéressant, afin d'avoir une image plus globale du marché, d'énumérer quelques acteurs du covoiturage qui travaillent sur la même problématique, à savoir le covoiturage de courte distance, cependant avec des approches différentes. Par exemple Klaxit (https://www.klaxit.com/), Karos (https://www. karos.fr/) et Mobicoop (https://www.mobicoop.fr/) proposent des plateformes 100% digitales destinées aux trajets domicile-travail. Ensuite, Citygo (https: //www.citygo.io/) et La Roue Verte (https://www.laroueverte.com/) se concentrent sur le covoiturage domicile-loisirs. La différence avec Ecov réside du fait qu'Ecov propose une solution qui s'adapte aussi bien au domicile-travail qu'au domicile-loisirs tout en combinant une infrastructure physique et digitale.

Besoin industriel et contexte statistique

Le présent travail de thèse est motivé par le besoin, d'une part de fiabiliser les informations disponibles pour les passagers (temps d'attente, flux disponible) au moment d'un covoiturage et d'autre part, de comprendre le comportement des conducteurs (heures de passage, arrêts desservis). L'enjeu majeur pour Ecov est d'être capable de faire une bonne prévision des temps d'attente. Comme nous l'avons vu précédemment, les passagers sont confrontés à un flux de conducteurs (non professionnels) qui sont sur leurs trajets du quotidien. On voit donc qu'il ne s'agit pas d'un simple calcul de temps d'attente grâce à une API ("Application Programming Interface"). Afin de répondre aux besoins industriels, cette thèse s'articule autour de cinq thématiques principales:

1. Les techniques de régression quantile pour la prédiction des temps d'attente.

2. La construction d'un processus de travail en vue d'utiliser les données issues du covoiturage.

3. La construction d'un modèle hiérarchique bayésien pour la prédiction des flux de trafic et des temps d'attente.

Introduction Nous reformulons cette question par "Quelle est la probabilité qu'un passager attende moins de y minutes avec α% de chance?". Les méthodes de régression standard s'intéressent à la moyenne de la variable cible et ne répondent pas à la question. Afin d'y répondre de la manière la plus complète possible, nous avons retenu la technique de régression quantile. On note α ∈ (0, 1) le niveau du quantile de Y et le α-quantile conditionnel de Y |X X X, défini par

q α (Y |X X X) = inf t∈R {F Y |X X X (t) ≥ α}, (1.1)
avec F Y |X X X la fonction de répartition de Y conditionnellement à X X X. Contrairement aux techniques de régression standard, où nous utilisons la méthode des moindres carrés, la régression quantile utilise une autre méthode de minimisation du quantile. Le problème de minimisation du quantile est défini par

q α = argmin a∈R E[ρ α (Y -a)],
dans le cas non-conditionnel et

q α (Y |X X X) = argmin f ∈L 1 (X X X) E[ρ α (Y -f )],
dans le cas conditionnel, où ρ α (y) = y (α-1 [y≤0] ) est la fonction de perte "Pinball" et L 1 (X X X) = {f : E[|f (X X X)|] < ∞}. La Figure 1.7 représente la fonction de perte "Pinball" pour différentes valeurs de α. Nous proposons une démonstration du Introduction problème de minimisation du quantile à la fois dans le cas non-conditionnel et conditionnel.

Figure 1.7 -Fonction de perte "Pinball" ρ α pour différentes valeurs de α. Bleu : α = 0, 95, orange α = 0, 75, vert α = 0, 5, rouge α = 0, 25.

Après avoir présenté le principe de régression quantile, nous nous arrêtons sur quatre approches possibles. La première et la plus ancienne est celle décrite dans [START_REF] Koenker | Quantile regression[END_REF]. Le quantile conditionnel est supposé linéaire tel que q α (Y |X X X) = f (X X X, β α ) avec f (X X X, β α ) = X X X β α . La deuxième est nommée "Quantile Regression Forest". Elle est introduite dans Meinshausen, 2006a et se base sur les Forêts Aléatoires de [START_REF] Breiman | Random forests[END_REF]. La différence essentielle entre cette approche et les Forêts aléatoires classiques est la pondération finale des arbres composant la forêt quantile. Dans cette même catégorie, nous présentons aussi la régression quantile par Forêts Aléatoires Généralisées dû à [START_REF] Athey | Generalized Random Forests[END_REF]. Cette fois-ci une toute nouvelle manière de construire les arbres aléatoires est développée et se base sur un critère de découpe prenant en compte la cible, le α-quantile de Y . Enfin, la dernière approche est le Gradient Boosting adapté à la régression quantile. Le principe consiste à appliquer la méthode classique de descente de gradient sur la fonction de perte "Pinball". De ce fait, on obtient en effet un estimateur du α-quantile.

Une fois ces notions introduites, nous détaillons les variables explicatives utilisées. Essentiellement, ces variables proviennent des demandes de covoiturages des passagers. Elles sont classées en trois catégories:

1. Les informations géographiques de la demande de covoiturage.

2. Les informations temporelles de la demande de covoiturage.

Les informations concernant les habitudes des passagers.

Par exemple, les couples origine-destination des demandes de passagers sont nombreuses pour la structure du réseau de covoiturage. Afin d'en tirer une information permettant d'améliorer les prédictions des modèles nous avons procédé a un partitionnement spatial des données. La Figure 1.8 montre la carte avec les sept groupes qui représentent les sept quartiers de destinations possibles. Ce partitionnement réduit le nombre de paires origine-destination observées de 252 à 66.

Figure 1.8 -Les quartiers de destination regroupés dans le réseau de covoiturage Vexin.

Introduction

Chapitre 3: Relaxing door-to-door matching reduces passenger waiting times: a workflow for the analysis of driver GPS traces in a stochastic carpooling service Ce chapitre a fait l'objet d'une publication, dans le journal "Transportation Engineering", co-écrite avec Tarn Duong, Safa Fennia et Constant Bridon. Il s'agit de trois membres de l'équipe Data/SIG d'Ecov.

À présent, nous travaillons sur les données issues des lignes de covoiturage. Les nouveautés par rapport au chapitre précédent sont; la collecte des traces GPS des conducteurs et la structure améliorée du service. Les traces GPS permettent de compléter les informations concernant les passagers, la géographie du territoire et la temporalité des demandes de covoiturage. Une trace GPS est une information spatio-temporelle très riche et l'utilisation de celle-ci permet une meilleure compréhension des habitudes des conducteurs. Ainsi, le traitement et l'analyse de ces données nous ont permis d'obtenir deux résultats essentiels pour le covoiturage:

• La correspondance porte-à-porte de trajectoires complètes de l'origine à la destination est un obstacle structurel pour la transformation du covoiturage en un service de transport en commun.

• Les arrêts de covoiturage fixes favorisent la rencontre entre les conducteurs et les passagers et diminuent les temps d'attente.

Pour illustrer les difficultés de la correspondance passager-conducteur dans l'espace et dans le temps pour les trajectoires porte-à-porte, nous pouvons la représenter par la partition d'un cube 3D divisé en sous-cubes plus petits, où l'axe x est la longitude, l'axe y la latitude et l'axe z le temps, comme le montre la Figure 1.9. Ainsi, pour que deux trajectoires correspondent dans l'espace et dans le temps dans un sens porte-à-porte, elles doivent partager le même souscube pour l'origine, et de même pour la destination. Ainsi, une correspondance porte-à-porte plus stricte entraîne une diminution du nombre de conducteurs disponibles pour partager leurs trajectoires avec des passagers. Ensuite, nous étudions les traces complètes des conducteurs et nous montrons l'intérêt des points de rencontre fixe. La Figure 1.10 illustre la classification hiérarchique spatiale faite sur 121 traces GPS de conducteurs. Nous voyons que le potentiel de covoiturage est bas lorsque nous nous basons sur la correspondance porte-à-porte entre les conducteurs et les passagers. Dans le cas de la correspondance porte-à-porte nous obtenons 9 clusters distincts dont certains représentent une trace unique. Au final, cette étude nous a permis de conclure que le potentiel de covoiturage, grâce aux arrêts de covoiturage fixe, a augmenté de 59%. Le processus de travail mis en place génère quatre résultats:

1. La visualisation des traces GPS sur des fonds de cartes.

2. Une fiche horaire pour des intervalles de tailles différents (à l'heure, la demi-heure, le quart d'heure) des flux de conducteurs par ligne.

Introduction 3. Une fiche horaire pour des intervalles de tailles différents (à l'heure, la demi-heure, le quart d'heure) de prédiction de temps d'attente.

4. L'évolution du temps d'attente en fonction du taux de participation des conducteurs sur les lignes de covoiturage.

Il est important de souligner que ce processus de travail proposé est la base de tous les travaux des chapitres suivants. Nous rappelons qu'en plus de la construction du processus de travail ce chapitre permet essentiellement de valider l'importance de pratiquer le covoiturage avec des points de rencontre fixes sur le territoire. Ceci permet effectivement de faire converger les conducteurs et les passagers afin de simplifier leur rencontre. En second lieu, nous avons aussi validé grâce aux données empiriques que ces points de rencontre fixes permettent de faire baisser les temps d'attente des passagers. Nous proposons aussi une estimation de l'ordre de grandeur de la baisse potentielle des temps d'attente en fonction du taux de participation des conducteurs. On donne une illustration de cette estimation à la Figure 1.12.

Figure 1.12 -Évolution du temps d'attente prédit des passagers en fonction du taux de participation des conducteurs sur la ligne de covoiturage Bourgoin > St-Priest, pendant les heures de fonctionnement du matin 06h30-09h00, du 2019-07-25 au 2020-02-17.

Chapitre 4 : Bayesian hierarchical models for the prediction of the driver flow and passenger waiting times in a stochastic carpooling service

Introduction

Nous souhaitons à présent construire un modèle pouvant apporter des informations plus précises concernant les flux de trafic sur une ligne de covoiturage ainsi que les temps d'attente qui en découlent. Les données disponibles sont de sources distinctes

• les traces GPS des conducteurs,

• les temps d'attente observés.

Les données concernant les conducteurs sont des séries temporelles représentant le flux de trafic quotidien sur la ligne de covoiturage. Nous représentons une telle série temporelle par la Figure 1.13, où un code couleur est utilisé pour différencier les différents types de jour. Cette classification des types de jours est importante. En effet, dans le domaine du covoiturage les habitudes des conducteurs sont influencées par le type de jour. Le premier étage s'intéresse à la modélisation des flux de trafic quotidien prenant en compte les différents types de jour. Il s'agit d'un modèle ayant une structure autorégressive d'ordre K ≥ 1. Afin de le définir, nous notons y i le flux de trafic quotidien pour le jour i. La modélisation du premier étage est définie par

y i = α DT(i) K k=1 η DT(i-k) y i-k + ε i ,
(1.2) avec ε i indépendantes et identiquement distribuées selon une loi Gaussienne centrée de variance σ 2 ε et DT une fonction qui caractérise les types de jours,

DT(i) =       
ORD si le jour i est ordinaire, SCH si le jour i est vacance scolaire, PWE si le jour i est férié ou weekend.

(1.3)

Le paramètre α DT est le coefficient pour le jour actuel i et η représente le coefficient de transition des types de jour qui précède le jour i sur le flux y i . Le second étage concerne la modélisation des temps d'attente. Nous distinguons à partir de maintenant les temps d'attente perçus et les pseudos temps d'attente. Le pseudo temps d'attente est défini comme le temps d'attente d'un passager sans la contrainte de queue d'attente. Notons n i le nombre de demandes des passagers pour le jour i et t i,1 < • • • < t i,n i les moments de la journée de ces demandes. Soit t i,j l'heure d'arrivée du conducteur pour la demande du passager au moment t i,j , i = 1, . . . , N et j = 1, . . . , n i . Le temps d'attente perçu w * i,j et le pseudo temps d'attente w i,j pour la demande du passager au moment t i,j sont définis par w * i,j = t i,j -t i,j , w i,j = t i,j -max(t i,j , t i,j-1 ).

La modélisation retenue pour ce second étage est la régression Gamma. Cette dernière permet d'exprimer une relation d'inverse proportionnalité entre les Introduction temps d'attente et les flux de trafic. En effet, cette relation est visible sur les données observées sur les différentes lignes. Soit la période de 24 heures d'une journée divisée en S intervalles égaux I 1 < • • • < I S , la fraction du flux quotidien de conducteurs y i sur chaque intervalle est y i β s avec β s ≥ 0 avec S s=1 β s = 1. Conditionnellement au flux de trafic y i , au temps de demande des passagers t i,j ∈ I s et β β β = (β 1 , . . . , β S ), nous supposons que les pseudos temps d'attente w i,j sont des variables aléatoires Gamma indépendantes de paramètres ν et β s y i , c'est-à-dire que Nous rappelons que nous nous sommes concentrés sur les pseudos temps d'attente w i,j , bien que les temps d'attente perçus w i,j soient plus pertinents, puisque ces derniers sont les vrais temps d'attente du point de vue du passager. Pour cela, en plus du modèle hiérarchique bayésien, nous avons proposé un cadre pour l'analyse future des temps d'attente perçus en nous basant sur des simulations poissonniennes des arrivées des passagers. Puisque la distribution prédictive a posteriori complète est disponible, nous pouvons l'intégrer par rapport à la distribution des arrivées passagers pour obtenir les pseudos temps d'attente. Ensuite, nous sommes alors en mesure de reconstruire les temps d'attente perçus à l'aide de l'équation suivante

w i,j |(y i , β β β, t i,j ∈ I s ) ∼ Γ(ν, β s y i ) (1.4) avec i = 1, . . . , N et j = 1, . . . ,
w * i,j+1 = w i,j+1 + [w i,j -ζ i |(w i,j > ζ i )] avec ζ i = t i,j+1 -t i,j .
La contribution essentielle du chapitre est la prédiction des flux quotidiens des conducteurs et des pseudos temps d'attente des passagers à l'aide d'un modèle hiérarchique bayésien de deux étages. La principale difficulté était de combiner deux sources (traces GPS des conducteurs et temps d'attente des passagers) ayant une temporalité différente. Une deuxième contribution est celle concernant la reconstruction des temps d'attente perçus à l'aide du cadre bayésien qui nous donne accès à la totalité de la loi prédictive a posteriori. Dans le cas des statistiques bayésiennes nous pouvons introduire de l'information sur les paramètres du modèle grâce à la construction d'un prior informatif. Dans ce chapitre nous développons un prior informatif afin d'améliorer les performances d'un modèle bayésien lors du lancement d'une nouvelle ligne de covoiturage. En notant L, le jeu de données long qui va servir à la construction du prior, et C le jeu de données court, où on va appliquer le transfert bayésien. Le prior informatif qu'on note π C pour le jeu de données C est construit à partir de la loi a posteriori π L (•|X L ) estimée sur le jeu de données L. On suppose que les lois a posteriori convergent vers leurs lois limites dans le cas du jeu de données C. En revanche, dans le cas de la situation d'historique court, les lois a posteriori n'ont pas eu le temps de converger. On suppose de ce fait que les données des nouvelles lignes (historique court) ressemblent aux données des lignes ayant un historique long. Les paramètres d'intérêts pour la construction sont respectivement θ L et θ C . La construction du prior informatif s'effectue en deux temps:

• On commence par ajuster le modèle non-informatif sur le jeu de données L.

• On construit ensuite le prior informatif qu'on va détailler ensuite pour le jeu de données L.

La loi informative est construite à partir d'une loi a priori hiérarchique normale qui s'appuie sur la moyenne µ L et et la variance

Σ L de la loi a posteriori π L (θ|X L ). Sa forme est θ|k, r ∼ N (kµ L , rΣ L ), k|m, l ∼ N (m(1, . . . , 1) , l Id).
Les lois a priori des hyperparamètres r, m et l sont définies par

r ∼ N (a r , b r ), m ∼ N (1, σ 2 m ), l ∼ N (1, σ 2 l ),
PE(δ) = 1 Ñ Ñ i=1 1{|q α ( wi ) -w i | < δ} (1.6)
avec q α ( wi ) est le quantile au niveau α de la loi prédictive a posteriori de la i ième observation et w i le i ième temps d'attente observé pour un échantillon de taille Ñ .

Introduction

Le quantile sélectionné pour comparer les scores PE est le quantile niveau 95% des lois a posteriori prédictives des temps d'attente. La Figure 1.24 compare les scores PE du quantile niveau 95% issues du modèle informatif et non-informatif pour le jeu de données B. Nous constatons que le score du modèle informatif est meilleure que celle du modèle non-informatif. Ce résultat montre que la loi a priori informative apporte une réelle information supplémentaire. • La récupération des données et leurs pré-traitements (dans le cas des données réelles) ou la simulation des données synthétiques (dans le cas des données simulées).

• L'entraînement du modèle hiérarchique bayésien sur les données récupérées ou simulées.

• La génération des lois a posteriori des paramètres du modèle ainsi que les lois prédictives a posteriori des temps d'attente.

• La génération des quantiles prédictifs a posteriori souhaités provenant de la loi prédictive a posteriori des temps d'attente.

• Le contrôle prédictif a posteriori du modèle hiérarchique bayésien.

Ensuite, nous nous intéressons à l'automatisation de l'utilisation du paquet bayesian-hierarchical-model. En effet, il est de grand intérêt de pouvoir automatiser le processus de prédiction des temps d'attente. Nous rappelons que les résultats obtenus par le modèle prédictif sont consommés par plusieurs services d'Ecov (site internet, application mobile, borne d'informations voyageurs, logiciel de tableau de bord). La solution retenue pour automatiser l'exécution du

Introduction

The carpooling services provided by Ecov1 differ from the usual carpooling services, since they are based on stochastic matching between driver(s) and passenger(s). To facilitate this matching, physical meeting points are constructed on the roadside. These physical meeting points are then equipped with electronic signs which display the carpooling requests made by the passengers to drivers on the road. The potential drivers thus can see a carpooling request in real-time and can decide to stop and embark the passenger(s) directly. This shifts away from a deterministic matching used by Uber-like and taxi-like services. The placement of these meeting points are decided with local government authorities in order to respond to the mobility requirements in each local area, and they constitute a carpooling network as illustrated in Figure 2.1. The objective of Ecov is to design a carpooling service for frequent, short distance journeys (from 10km to 40km roughly) in sparsely populated peri-urban and rural areas. Its first carpooling network was located in the northwestern of France, a region called "Vexin" as shown in Figure 2.2. In the early experimental phase of the implementation of these carpooling services, the passengers are even allowed to request a carpooling trip to any destination within the service area, so the carpooling network is even more complicated than that represented in the schematic diagram in Figure 2.1 where the connected edges indicate carpooling was provided only between the meeting points. Whilst this is highly convenient for the passengers, the complete freedom in the choice of the destination leads to an explosion of the number of possible origin-destination pairs which a reliable carpooling service must be assured.

A key component in the reliability of a carpooling (or indeed an transportation service) is the accuracy of the predicted waiting times of passengers, since Chapter 2 -Quantile Regression for waiting time prediction in a stochastic carpooling network the latter often prioritise the estimated time of arrival (ETA) when choosing between different modes of transport, along with the financial cost. For Ecov, as a transport provider, it is crucial to (i) quantify accurately the passenger waiting times and (ii) to reduce any excessive waiting times to a level suited to the market. These reliable waiting times are the foundation for building, 2.2. Quantile Regression maintaining and expanding the carpooling service.

We focus on the first of these requirements. Therefore, we must be able to answer a passenger's question "What is the waiting time for my carpooling request to my destination if I make the request now from this location?" as accurately as possible. Whilst it appears that this can be answered by supplying the mean of the observed waiting times, this can hide the large variability that is inherent in predicting waiting times for the arrival of the suitable carpooling driver without any prior matching. To supply accurate information to the passenger based on highly variable stochastic phenomena, Ecov wishes to provide the response to the passenger's question above in the form of a probabilistic statement of the type that the passenger will wait less than y minutes with α% chance.

Mathematically, this α ∈ (0, 1) defines the α-quantile of Y , where Y is the random variable of the carpooling request waiting time. The quantile of Y is more precisely defined by

q α (Y ) = inf t∈R {F Y (t) ≥ α}, (2.1)
where

F Y is the cumulative distribution function of Y , i.e. F Y (y) = P(Y ≤ y)
. The availability of auxiliary information, e.g. about the geographical and temporal characteristics of the carpooling request, can greatly assist in the prediction of the quantiles of the waiting times. This auxiliary information can be incorporated into a single, unified framework of Quantile Regression. The waiting time predictions from a Quantile Regression are more detailed than simple mean predictions, which can be utilised to reassure passengers about the uncertainties in the stochastic matching in an Ecov carpooling service. In Section 2.2 we focus on the definition of Quantile Regression as a minimisation problem, and then elaborate the solutions provided by the classical linear model approach, the Random Forests, the Generalised Random Forests approaches, and the gradient boosting approach. In Section 2.3 we describe in further detail the data collected from the Vexin carpooling network and then compare these different implementations of Quantile Regressions for these carpooling requests. We end with some concluding remarks and some details of the supporting mathematical results.

Quantile Regression

Classical regression relies on a well-known result that the mean of the response random variable Y is the minimiser of the expected squared error loss problem

E[Y ] = argmin a∈R E[(Y -a) 2 ].
Chapter 2 -Quantile Regression for waiting time prediction in a stochastic carpooling network That is, E[Y ] is the value that has smallest possible expected squared loss from the response Y . Then it is straightforward to pose the question of what the result would be if the squared distance were replaced by other distance measures. It is also well-known that for the absolute loss, then the result is the median of Y , i.e.

med[Y ] = argmin a∈R E[|Y -a|].
Since the median is the special case of the 0.5-quantile, then a further question to ask is which loss function would result in a generalisation of this median result to an arbitrary α-quantile, for any 0 ≤ α ≤ 1. This is achieved by the Pinball loss

ρ α (y) = y(α -1 [y≤0] ), (2.2)
where 1 [•] is the indicator function. The Pinball loss in Equation (2.2) is a convex loss function so it leads to a well-defined minimisation problem. In An equivalent definition for the quantile given by Equation (2.1) is given by the 2.2. Quantile Regression solution of the Pinball loss minimisation problem, as asserted in Lemma 1.

Lemma 1. Let Y be an integrable real-valued random variable whose distribution function F Y is strictly increasing. Let α ∈ (0, 1). Then the quantile q α is the unique value which satisfies

q α = argmin a∈R E[ρ α (Y -a)].
Whilst this is an established result, we provide its proof in Appendix 2.5.1.

In practice, the complete random variable Y is not known, and we only have observed instantiations of it, i.e. Y 1 , . . . , Y n which is a random sample from F Y . As a result of Lemma 1, we can define an estimator of the quantile by a principle of M-estimation (empirical risk minimisation). Thus, the empirical estimator of the α-quantile is

qα ∈ argmin a∈R 1 n n i=1 ρ α (Y i -a) .
Van der Vaart (2000, Chapter 21) establishes the weak consistency of this estimator under weaker conditions than those of Lemma 1. Koenker, Chesher, et al. (2005, Chapter 4) establishes the consistency of this estimator, i.e. ||q α -q α || → 0 as n → ∞ under the conditions of Lemma 1.

For the regression context, the random vector X X X, known as the explanatory variables or covariates, provides auxiliary information of the behaviour of the response variable Y . Analogous to the situation above where unconditional mean E[Y ] is the minimiser of an expected squared loss problem, then conditional mean E[Y |X X X] is the minimiser defined by

E[Y |X X X] = argmin f ∈L 2 (X X X) E[(Y -f ) 2 ],
where L 2 (X X X) = {f : E[f (X X X) 2 ] < ∞} is the space of all square integrable functions with respect to X X X. It is straightforward to replace Y by Y |X X X in Equation (2.1) to obtain the conditional α-quantile as

q α (Y |X X X) = inf t∈R {F Y |X X X (t) ≥ α}, (2.3) where F Y |X X X is the conditional distribution of Y |X X X.
However, this definition is not tractable for many reasons. Therefore we propose in Lemma 2 to provide an equivalent definition of Equation (2.3). As described in Armerin, 2014, it is a generalisation of the Pinball loss minimisation problem for the conditional α-quantile.

Lemma 2. Let Y be an integrable real-valued random variable whose conditional distribution F Y |X X X is almost surely increasing given a random vector X X X. Let α ∈ (0, 1).

Then the conditional quantile q α (Y |X X X) satisfies

q α (Y |X X X) = argmin f ∈L 1 (X X X) E[ρ α (Y -f )],
Chapter 2 -Quantile Regression for waiting time prediction in a stochastic carpooling network where L 1 (X X X) = {f : E[|f (X X X)|] < ∞} is the space of all integrable functions with respect to X X X.

The proof of Lemma 2 is given in Appendix 2.5.2. With this lemma, we are able to proceed to the definition of a Quantile Regression. Quantile Regression was introduced in [START_REF] Koenker | Regression quantiles[END_REF]and then deepened in Bassett Jr et al., 1982[START_REF] Koenker | Quantile regression[END_REF]. Just as classical regression methods are based on minimising the expected squared loss in order to estimate the conditional mean, the aim of Quantile Regression is to estimate conditional quantiles. Many authors have shown interest in Quantile Regression and have subsequently expanded on the standard techniques. For instance, the definition of the confidence intervals for quantiles regression estimates in [START_REF] Kocherginsky | Practical confidence intervals for regression quantiles[END_REF] Another is the adaptation for sequential time series forecasting proposed in [START_REF] Biau | Sequential quantile prediction of time series[END_REF]. More recently tree-based algorithms for Quantile Regression have been developed. For example, in [START_REF] Meinshausen | Quantile regression forests[END_REF][START_REF] Athey | Generalized Random Forests[END_REF] two different algorithms based on a Random Forests algorithm are proposed for Quantile Regression. Moreover, in Zheng, 2012 an adaptation of Gradient Boosting Machine for Quantile Regression is offered. We recall that in the regression context the aim is to give an estimate of q α (Y |X X X), that we denote qα (Y |X X X), by using the observed responses Y i given the explanatory variables X X X i ∈ R p with i = 1, . . . , n. This scalar represents an estimation of the quantile of waiting times. For instance, if we set α = 0.9 then q0.9 (Y |X X X) is the estimation of the 90% quantile. This value means that in 90% of the time, the observed waiting times does not exceed q0.9 (Y |X X X).

Linear Quantile Regression

The observed data set is denoted D n = {(X X X i , Y i ) : 1 ≤ i ≤ n}. Towards this goal of defining a regression model, we introduce a parametric definition of the quantile function, i.e. q α (Y |X X X) = f (X X X, β α ) where f is a known, non-random function satisfying the Lemma 2, and β α ∈ R p is the regression parameter. The objective is to estimate the β α parameter since it immediately yields an estimate of the quantile function as qα (Y |X X X) = f (X X X, βα ). The function f can take several forms, and we start with the historical approach, where it is a linear product of the covariate matrix and the regression parameter, that is f (X X X, β α ) = X X X β α . Following the Lemma 2, we obtain the empirical estimator of the regression parameter is

βα ∈ argmin β∈R p 1 n n i=1 ρ α (Y i -X X X i β) .
Contrary to the classical regression case, here there is not an explicit expression for the solution for the quantile regression, as pointed out in Koenker &

Quantile Regression

Hallock (2001). To solve this problem, efficient numerical methods like linear programming are often deployed. For example, [START_REF] Portnoy | The Gaussian hare and the Laplacian tortoise: computability of squared-error versus absolute-error estimators[END_REF] and [START_REF] Bosch | A convergent algorithm for quantile regression with smoothing splines[END_REF] use an interior point method to compute βα . Despite that βα does not have an explicit expression, much is known about it. Its asymptotic normality and its rate of convergence is detailed in Koenker, Chesher, et al. (2005, Chapter 4) under the following conditions: C1: The distribution function F Y is absolutely continuous at Y i , and its density is continuous and is uniformly bounded away from 0 and ∞ at the points q α (Y i |X X X i ).

C2: There exist positive definite matrices D 0 and D 1 (α) such that, as n → ∞,

• n -1 n i=1 X X X i X X X i → D 0 • n -1 n i=1 f Y (q α (Y i |X X X i ))X X X i X X X i → D 1 (α) • max i=1,...,n n -1/2 X X X i → 0,
where f Y is the density function of Y . When these two conditions are verified we obtain the asymptotic normality, i.e.

n 1/2 ( βα -β α ) N (0, α(1 -α)D -1 1 D 0 D -1 1 ).
The detailed proof of this result is given in Koenker, Chesher, et al. (2005, Chapter 4) and Van der Vaart (2000, Chapter 21). As for classical regression, [START_REF] Sherwood | Traffic flow forecasting: comparison of modeling approaches[END_REF] proposed an adaptive Lasso penalisation for Quantile Regression.

Let λ ≥ 0 be the penalisation parameter, then the Lasso minimisation to obtain the regression parameter estimate is expressed as

βα ∈ argmin β∈R p    1 n n i=1 ρ α (Y i -X X X i β) + λ p j=1 |β j |    .
Analogously Q. [START_REF] Li | Bayesian regularized quantile regression[END_REF] introduced the Elastic-Net penalty in a Bayesian framework: let λ 1 , λ 2 ≥ 0, then βα ∈ argmin

β∈R p    1 n n i=1 ρ α (Y i -X X X i β) + λ 1 p j=1 |β j | + λ 2 p j=1 β 2 j    .
Another Bayesian approach has been developed by [START_REF] Yu | Bayesian quantile regression[END_REF], and a Support Vector Machine algorithm by [START_REF] Crambes | Support vector machine quantile regression approach for functional data: Simulation and application studies[END_REF] We leave these penalisation methods and focus on Random Forest and Gradient Boosting methods in the following sections.

Chapter 2 -Quantile Regression for waiting time prediction in a stochastic carpooling network

Quantile Regression with Random Forests

The construction of Quantile Regression Forests is introduced in [START_REF] Meinshausen | Quantile regression forests[END_REF] and is an extension of the Random Forests algorithm [START_REF] Breiman | Random forests[END_REF]. In comparison to linear Quantile Regression, this approach is a non-parametric method to estimate the conditional quantiles. Quantile Regression Forests is an ensemble method that aggregates a collection of decision trees. This aggregation differs from the classical approach of Random Forests which aim to estimate the conditional expectation E[Y |X X X] via a weighted mean of the responses Y i , i = 1, . . . , n. Indeed, Quantile Regression Forests aim to estimate the full conditional distribution using this same weighted mean of the Y i , and hence the conditional quantiles q α (Y |X X X) follow as a corollary.

The mechanism of Quantile Regression Forests is borrowed from the Random Forests algorithm [START_REF] Breiman | Random forests[END_REF]. A Quantile Regression Forest is composed of several trees where each decision tree is trained on D * n which is resampled uniformly (with or without replacement) from the observed data D n . Each tree is a sequence of decisions which determines a partition of the resampled data

D * n . Let denote this partition C b 1 ∪ • • • ∪ C b M
for the b th tree. Each decision is a binary split that takes places in the middle of two consecutive data points in order to avoid possible ties. The splits are made recursively by maximising a CART-decision rule. For more details concerning the CART-decision rule and variable selection, see [START_REF] Biau | A random forest guided tour[END_REF]. The splitting stops when the current cell contains fewer points than a threshold called the nodesize ∈ {1, . . . , a n } with a n the number of re-sampled data points in each tree. Each partition class C b is also called a leaf of the tree, also denoted b for the th leaf of the b th tree with b = 1, . . . , B. Every leaf is a partition class that can also be interpreted as a neighbourhood of a subset of the observations. For a classical Random Forests Regression, for tree b, for a new data point x x x, a set of weights w i (x x x, b) are defined by

w i (x x x, b) = 1 {X X X i ∈ b (x x x)} | b (x x x)| . (2.4)
For the forest composed of B trees, the corresponding weights are obtained by summing over all individual trees

w i (x x x) = 1 B B b=1 w i (x x x, b), (2.5)
with n i=1 w i (x x x) = 1. Then the prediction by the forest (i.e. an estimate of the conditional expectation E[Y |X X X]) is the weighted sum over all observations of D * n . Let denote Ŷ (x x x) the prediction of the new data point x x x and is obtained by

Ŷ (x x x) = n i=1 w i (x x x)Y i . 2.2. Quantile Regression Figure 2.4 is an illustration of the two partitions C 1 1 ∪ • • • ∪ C 1 8 and C 2 1 ∪ • • • ∪ C 2 8
made from two decision trees trained on a subset {X X X 1 , . . . , X X X 19 }, and where the remaining observed data {X X X 20 , . . . , X X X 23 } are excluded from the construction of these trees. For simplicity we consider the subset to be 2-dimensional and denote by X X X (1) , X X X (2) the two variables used for the splitting. In the first tree on the left of Figure 2.4, the observations X X X 4 , X X X 10 , X X X 12 share the same neighbourhood C 1 6 . The new point x x x falls within this neighbourhood. According to Equation (2.4), with b = 1, the weights are w i (x x x, 1) = 1 3 for i = 4, 10, 12 and 0 for observations outside the neighbourhood C 1 6 . In the same way, for the second decision tree on the right of Figure 2.4, the new observation x x x falls within the neighbourhood C 2 7 which contains X X X 4 and X X X 18 . Then the weights are w i (x x x, 2) = 1 2 for i = 4, 18 and 0 for observations outside the neighbourhood C 2 7 . Consider a forest composed by these two decision trees, then the weights at x x x according Equation (2.5) are w 4 (x x x) = 5 12 , w 10 (x x x) = 1 6 , w 12 (x x x) = 1 6 and w 18 (x x x) = 1 4 . Then the final estimate of the conditional expectation at x x x by this forest is

Ŷ (x x x) = 5 12 Y 4 + 1 6 Y 10 + 1 6 Y 12 + 1 4 Y 18 .
* 19 = {X X X 1 , • • • , X X X 19 } (black circles).
The remaining observed data {X X X 20 , . . . , X X X 23 } (gray circles) are excluded from their construction. The resulting partition classes are

C 1 1 ∪ • • • ∪ C 1 8 and C 2 1 ∪ • • • ∪ C 2 8 . A new data point is x x x (black cross).
In the context of Quantile Regression Forests, we must estimate the full conditional distribution F Y |X X X and then the conditional quantile. Through the definition of the conditional distribution Y |X X X we obtain

F Y |X X X (y) = P(Y ≤ y|X X X) = E[1 {Y ≤y} |X X X].
Since this is a conditional expectation of an indicator function, the connection with classical Random Forests is more apparent. The estimation of F Y |X X X given x x x Chapter 2 -Quantile Regression for waiting time prediction in a stochastic carpooling network is obtained from the weighted mean over the observations of 1 {Y ≤y} instead of Y , that is

FY |x x x (y) = n i=1 w i (x x x)1 {Y i ≤y} .
(2.6)

The last step is to plug this estimator of the conditional distribution in Equation (2.3) to obtain the estimator of the conditional α-quantile as

qα (Y |x x x) = inf t∈R { FY |x x x (y) ≥ α}.
(2.7) Continuing with the same example in Figure 2.4 for a classical Random Forest, the estimate of the conditional distribution at the new data point x x x is FY |x x x (y) = 5 12 1 {Y 4 ≤y} + 1 6 1 {Y 10 ≤y} + 1 6 1 {Y 12 ≤y} + 1 4 1 {Y 18 ≤y} ), since the weights remain the same between the classical Random Forest and the Quantile Regression Forest. Moreover, under this approach, the prediction intervals for the new data point x x x are straightforward to compute. Denote α min for the lower level and α max the upper level, then we have directly

qα min (Y |x x x) = inf t∈R { FY |x x x (y) ≥ α min } qαmax (Y |x x x) = inf t∈R { FY |x x x (y) ≥ α max }.
The prediction interval Î(x x x) = [q α min (Y |x x x), qαmax (Y |x x x)] is well-defined. Indeed, since the distribution FY |X X X is a continuous piecewise and increasing function, it is clear that qα min (Y |x x x) ≤ qαmax (Y |x x x) whenever α min ≤ α max .

Quantile Regression with Generalized Random Forests

Generalized Random Forests, introduced in Athey et al., 2019, attempt to generalise Random Forests to several tasks: a more sophisticated non-parametric Quantile Regression, a conditional average partial effect estimation and a heterogeneous treatment effect estimation via instrumental variables. The general aim of Generalized Random Forests is to estimate the parameter function θ(•) by solving a local moment equation of the form

E[ψ θ(X X X i ),ν(X X X i ) (Y i )|X X X i = x x x] = 0,
(2.8)

where x x x is a test point and ψ θ,ν is the scoring function which is indexed by θ and ν (an optional nuisance parameter). Since we are focused on generalising the form of the non-parametric Quantile Regression, the optional nuisance parameter ν is not required, and in the sequel, we refer to

E[ψ θ(X X X i ) (Y i )|X X X i = x x x] = 0
as the local estimation problem. The required scoring function for a nonparametric α-Quantile Regression is ψ

θ(X X X i ) (Y i ) = α-1 {Y i ≤θ(X X X i )} .
We demonstrate in Appendix 2.5.3 that the scoring function ψ θ(•) (•) is indeed well adapted for an α-Quantile Regression. This demonstration is based on the fact that the scoring function ψ θ(•) (•) is intimately related to the Pinball loss ρ α (•). So the local estimation problem becomes

E[α -1 {Y i ≤θ(X X X i )} |X X X i = x x x] = 0.
(2.9) In practice, the objective of a Generalized Random Forest is to solve an empirical version of the local estimation problem defined by Equation (2.9). Similarly to Quantile Regression Forests, we use the weights w i obtained according to Equation (2.5) for the final aggregation of the forest. Thus the estimate of θ(•) is the solution of the empirical local estimation problem

θ(x x x) ∈ argmin θ(•) n i=1 w i (x x x)ψ θ(x x x) (Y i ) 2 ,
(2.10) with • 2 the Euclidean norm and the argmin is taken over the set of all measurable real-valued functions. Recall that these weights are a measure of the proximity between the re-sampled training data D * n and the new data point x x x. It is important to note here that the trees where these weights are computed for Generalized Random Forests are not built following the same process as Quantile Regression Forests. Indeed, a new tree-building procedure is detailed below, and the resulting trees are aggregated following Equation (2.10).

Generalized Random Forests offer a more convenient way to estimate quantiles based on the partitioning of the trees of a forest. A new splitting rule is proposed which is adapted to Quantile Regression, rather than the usual variance-based splitting rule, to minimise the heterogeneity within each child Chapter 2 -Quantile Regression for waiting time prediction in a stochastic carpooling network node, by maximising the heterogeneity between each child node. Each split starts at the parent node, denoted P , given the available data in this node D P = {(X X X i , Y i ) : i = 1, . . . , n P } with n P the number of available data points. We define θP the solution of the local estimation problem

θP ∈ argmin θ(•) {X X X i ∈D P } ψ θ(x x x) (Y i ) 2 .
(2.11)

Since that the splits are made in order to improve the quality of the estimate of θ(•), then this means that the splits should take into account the regression objective i.e. the α-quantile of Y . So we split the parent node P into two child nodes C 1 , C 2 which minimises the following error

err(C 1 , C 2 ) = 2 j=1 P[X X X ∈ C j |X X X ∈ P ]E[( θC j -θ(X X X)) 2 |X X X ∈ C j ],
where the θC j are the fit over the child nodes C j following Equation (2.11). The direct minimisation of this error is not possible since θ(x x x) is unknown. The alternative is to maximise the heterogeneity between the child nodes C 1 and C 2 .

The appropriate heterogeneity measure is the following delta criterion

∆(C 1 , C 2 ) = n C 1 n C 2 n 2 P ( θC 1 -θC 2 ) 2 .
(2.12) Athey et al. (2019, Section 3, Proposition 1) asserts that maximising

∆(C 1 , C 2 ) is equivalent to minimising err(C 1 , C 2 ).
This duality is illustrated in Figure 2.6, where we have two examples of splits on the parent nodes P and P containing the subset D P = D P = {X X X 1 , . . . , X X X 9 } (for simplicity the parent nodes are the same). The subset represented by the blue circles {X X X 2 , X X X 3 , X X X 6 , X X X 8 } should ideally be found after the split in a different child node to the black circles to achieve maximum heterogeneity between the child nodes. The split on the left, made on the parent node P , has a better separation than the split represented on the right on the parent node P . This leads to err(C 1 , C 2 ) < err(C 1 , C 2 ) and also its dual problem i.e.

∆(C 1 , C 2 ) > ∆(C 1 , C 2 ).
Whilst for Quantile Regression Forests, the heterogeneity between C 1 and C 2 is generally measured as the difference in the group variances, for the Generalized Random Forests this is measured by Equation (2.12). However, maximising this delta criterion ∆(C 1 , C 2 ) by brute force requires numerating all possible splits on parent node data D P . This is computationally infeasible since it also requires estimating θC 1 and θC 2 according Equation (2.11) for each split. A computationally feasible alternative is based on approximating this delta criterion ∆(C 1 , C 2 ), which in turn is based on approximations θC 1 and θC 2 of θ Figure 2.6 -Illustration of two splits on parents nodes P and P applied in the available data D P = D P = {X X X 1 , . . . , X X X 9 }. The blue circles should ideally be separated in one child node while the black circles in the other child node to achieve maximum heterogeneity between the child nodes.

in the child nodes. These last two quantities are the gradient approximations of θC 1 and θC 2 obtained by

θC j = θP - 1 |{i : X X X i ∈ C j }| {i:X X X i ∈C j } A -1 P ψ θP (Y i ),
where A P is a consistent estimate of the expectation of the gradient of ψ θ such that

E[A P ] = ∇E[ψ θ(X X X i ) (Y i )|X X X i ∈ P ] = ∇E[α -1 {Y i ≤θ(X X X i )} |X X X i ∈ P ].
The first (labelling) step consists of computing θP and the matrix A -1 P on the available data of the parent node D P following Equation (2.11), and from which we obtain

A -1 P ψ θP (X X X i ) = 1 {Y i > θP (X X X i ) } . The second step is a standard CART variance-based split which is applied to A -1 P ψ θP (X X X i )
. By construction, this split maximises

∆(C 1 , C 2 ) = 2 j=1 1 |{i : X X X i ∈ C j }|   {i:X X X i ∈C j } A -1 P ψ θP (X X X i )   2 . Athey et al. (2019, Proposition 2, p. 10) demonstrate that ∆(C 1 , C 2 ) and ∆(C 1 , C 2 )
are approximately equivalent under suitable conditions. So we are able to grow B individual trees according to maximisation of the approximate child heterogeneity outlined above. The final forest can now be built using the weights w i (x x x) that are computed as for classical Random Chapter 2 -Quantile Regression for waiting time prediction in a stochastic carpooling network Forest Regression, by applying Equations (2.4) and (2.5). The final estimate of the conditional α-quantile by a Generalized Random Forest at at new data point

x x x is qα (Y |x x x) = n i=1 w i (x x x)ψ θ(x x x) (Y i ),
(2.13)

with ψ θ(x x x) (Y i ) = α -1 {Y i ≤ θ(X X X i )} .

Quantile Regression with Gradient Boosting

The general aim of any boosting algorithm is to form an ensemble learner which is a linear combination of boosted weak learners. "Boosting" in this sense is to improve the accuracy of a weak learner such as a shallow decision tree. The boosted learner is built after a series of M iterations of an initial weak learner where it is improved at each iteration. Gradient Boosting is thus named because each iteration attempts to improve the learning process as guided by the gradient direction. It is a popular algorithm with several variants, all of which can be incorporated into Quantile Regression. Gradient Boosting is based on the work originally conducted on Adaboost [START_REF] Schapire | The strength of weak learnability[END_REF][START_REF] Freund | Boosting a weak learning algorithm by majority[END_REF][START_REF] Freund | A decision-theoretic generalization of online learning and an application to boosting[END_REF][START_REF] Schapire | Experiments with a New Boosting Algorithm[END_REF]. Since Breiman in 1997 [START_REF] Breiman | Arcing the edge[END_REF] made the fundamental observation that this procedure is a gradient-descent-type algorithm, the works of J. [START_REF] Friedman | Additive logistic regression: a statistical view of boosting (with discussion and a rejoinder by the authors)[END_REF][START_REF] Friedman | Greedy function approximation: a gradient boosting machine[END_REF][START_REF] Friedman | Stochastic gradient boosting[END_REF] formalised and baptised it with its current name of Gradient Boosting. These were refined by [START_REF] Bühlmann | Boosting algorithms: Regularization, prediction and model fitting[END_REF] and a new version named Extreme Gradient Boosting (XGBoost) was developed in [START_REF] Chen | XGBoost: A Scalable Tree Boosting System[END_REF]. More recent works have focused on the computational efficiency, where Biau, Cadre & Rouvière, 2019 proposed Accelerated Gradient Boosting (AGB), which in turn is based on Nesterov's accelerated gradient descent [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O(1/k 2 )[END_REF]. These computational efficiencies are not limited to AGB, and can also be applied to XGBoost. The incorporation of Gradient Boosting into Quantile Regression relies on the approach of Zheng, 2012 who replaced the traditional squared cost function, as for a classical regression, with the Pinball function ρ α (•). Thus the weak learners are boosted according the gradient of the Pinball function. Since Biau & Cadre (2020, Section 3) proved the convergence and consistency of Gradient Boosting algorithms for any welldefined cost function, then Quantile Regression with Gradient Boosting is guaranteed to inherit these desirable convergence properties.

Maintaining the previous notation, where

D n = {(X X X i , Y i ) : i = 1, .
. . , n} is the sample of observed data and f is the target function that we wish to approximate in the minimisation problem described in Lemma 2, the Gradient Boosting algorithm can be summarised in four steps. Firstly, at iteration m, for m ∈ 1, . . . , M , the steepest-descent step consists of computing the so-called "pseudo residuals" r r r [m] ∈ R p . The components of r r r [m] are the negative gradient of the Pinball function of the value of f obtained at iteration m -1, i.e.

r [m] i = - ∂ρ α (y i -f (X X X i )) ∂f (X X X i ) f =f [m-1] = [α1 {y i ≥f [m-1] (X X X i ) + (α -1)1 {Y i <f [m-1] (X X X i )} .
This steepest-descent is a greedy strategy since r r r [m] is the local direction in R p for which the Pinball loss is the most rapidly decreasing at f [m-1] . Secondly, a regression tree is fitted to the r

[m]

i . The terminal leaves of this regression tree are denoted j m , j = 1, . . . , j m , where j m is the number of leaves. Thirdly, the step size factor (or step length) η m is obtained as the solution to

η m = argmin η∈R ρ α (Y -(f [m-1] (X X X) -η)).
Fourthly, f [m] is updated following the rule f [m] = f [m-1] + η m r r r [m] .

That is, in iteration 1, we obtain f [1] , the first boosted approximation of the function f , and in iteration 2, the boosted approximation f [2] is obtained from r r r 1 , etc. If these four steps are iterated M times, then the resulting conditional quantile via this Gradient Boosting has an additive form, since we obtain

qα (Y |X X X) = f [M ] (X X X) = f [0] + M m=1
η m r r r [m] .

These steps are summarised in Algorithm 2.1.

Accuracy measures for Quantile Regression models

The output of a Quantile Regression is the prediction of an α-quantile. So the measure of the accuracy of a regression model is the based on the comparison of the empirical waiting times with the predicted α-quantile. The first metric we consider is the Pinball loss defined in Equation (2.2). For a Quantile Regression, the Pinball loss should be a small value for an accurate model. Let q(Y ) be the prediction for the α-quantile of an observation Y . The corresponding Pinball loss is

ρ α (Y -q(Y )) = (Y -q(Y ))(α -1 {Y ≤q(Y )} ). (2.14)
The interpretation of this Pinball loss is somewhat non-intuitive because it consists of a quantile-weighted difference between the predicted quantile q(Y ) and the observed value Y . Moreover, this weighting is different depending on whether the predicted quantile is above or below the observed value. For instance, if we perform a regression for the 75% quantile then if the predicted quantile is above the observed value, the weighting assigned by the Pinball loss is 0.25 and conversely, 0.75 if it is below. In other words, two models Chapter 2 -Quantile Regression for waiting time prediction in a stochastic carpooling network Algorithm 2.1: Quantile Regression with Gradient Boosting Input : Subset of observed data

D * n = {(X X X i , Y i ) : 1 ≤ i ≤ n}, desired quantile level α, number of iterations M Initialisation : Set f [0] (•) = 0 or α-quantile of D * n 1 for m=1 to M do 2 (1) Compute pseudo-residuals r [m] i = - ∂ρ α (y i -f (X X X i )) ∂f (X X X i ) f =f [m-1] = [α1 {Y i ≥f [m-1] (X X X i ) + (α -1)1 {Y i <f [m-1] (X X X i )}

3

(2) Fit regression tree to pseudo-residuals r

[m]
i , yielding terminal leaves j m with j = 1, . . . , j m

4

(3) for j = 1 to j m do 5 compute step size factor

η j,m = argmin η∈R      X X X i ∈ j m ρ α (Y i -(f [m-1] (X X X i ) -η))      6 end 7 (4) Update f [m] (X X X i ) = f [m-1] (X X X i ) + jm j=1 η j,m 1 {X X X i ∈ j m } 8 end Output : Estimated α-quantile function qα (Y |X X X) = f [M ] (X X X).
with different mean absolute errors may have the same Pinball loss. This is illustrated in subplots 1 and 4 in the toy example in Figure 2.7 where both of these models have the same Pinball loss. Whilst model 1 has larger distances between predicted and observed quantiles than model 4, the weight of the pinball losses in the model 4 are higher than those in model 1 because the former consistently underestimates. An alternative accuracy measure for Quantile Regression is called the Ramp loss and is defined by

1 n test ntest i=1 1 {Y i ≤q(Y i )}
which represents the empirical fraction of the quantile estimates which exceed the observed values Y i in the test dataset of size n test . This loss is introduced Biau & Patra, 2011 and used for sequential quantile prediction of time series. The Ramp loss for an α-Quantile Regression should be as close as possible to 1 -α to be an accurate model. For the regressions for the 75% quantile illustrated in models 1 and 3 of Figure 2.7, we have two predicted quantiles below their observed values. Given that n test = 8, then the Ramp losses are 2/8 = 0.25 in both cases. So based solely on the Ramp loss, models 1 and 3 have the same predictive performance. However a visual inspection reveals that this is not the interpretation that we wish to derive. In order to mitigate this, as well as the short-comings of the Pinball loss, we propose to examine both of these loss functions to inform our decisions. Consider again the example of the comparison of the models 1 and 4 in Figure 2.7. By considering the Ramp loss, we obtain directly that the better performing model is model 1. Likewise for the comparison between models 1 and 3. By additionally considering the Pinball loss, we can conclude that the model 3 is the most accurate model out of these four models. To compare the regression models, we use a k-fold cross validation procedure, which approximates the generalised error, i.e. the error that will be observed when the model is used in production. In summary this procedure consists in training the model on a systematic random sample of 90% of the observed data, and the predictions are carried out on the remaining 10% of the data (known as the test or prediction fold). These two steps are iterated for each 10% systematic sample as the test fold. For the kth test fold, the Pinball ρ α,k and Ramp errors R α,k are computed as the prediction errors. The generalised error is then obtained by averaging all the errors committed at each test fold. In addition to obtaining the generalised error we can also follow the evolution of the errors at each prediction fold and obtain information about the stability of Chapter 2 -Quantile Regression for waiting time prediction in a stochastic carpooling network the models with respect to their predictive performance. Figure 2.8 illustrates this cross-validation procedure.

Figure 2.8 -The k-fold cross validation procedure to calculate the generalised errors, Pinball and Ramp, for the competing models.

Comparison of Regression models for carpooling request waiting times

The main dataset is drawn from approximately 700 carpooling trips carried out in the Vexin carpooling network from 2016-01-18 to 2018-10-02, whose map is shown in Figure 2.2. The first data collection date corresponds to the service launch date. This service was in a highly experimental phase during this period in at least two aspects. Firstly because there was scarce existing data about the carpooling in this area. Secondly because Ecov's real-time carpooling service was at the time the world's first of its type. Recall that the carpooling network is composed of physical meeting points where the passengers can make their carpooling requests. These physical meeting points are accompanied of electronic panels placed in the roadside. When a passenger makes a request at the ticket machine, it is displayed on the panels. Drivers passing by are thus informed and can potentially respond in real-time. In this experimental

Comparison of Regression models for carpooling request waiting times

phase, the carpooling service allowed the passenger choose freely the desired destination from any location within the service area. Whilst this freedom makes the service highly convenient for passengers, it renders the problem of waiting time prediction extremely difficult because the combinatorial explosion in the number of possible origin-destination pairs.

Response variable: Waiting times

After a successful carpooling trip, i.e. a passenger is picked up and then dropped off by a driver, then we can define the waiting time as the difference between the time that the passenger embarks in the driver's car and the time of the initial carpooling request. The departure time is collected by the central data server from a manual notification sent by the passenger via the mobile phone application, whereas the carpooling request time is automatically sent from the ticket machine. A histogram of the waiting times is plot in Figure 2.9. These waiting times are skewed right in Figure 2.9 due the long right tail. Most of the waiting times are between 0 and 10 minutes, with a smaller peak around 15 minutes. The range is around 50 minutes which indicates a high variability in waiting times. For operational reasons, we consider that any values greater than 40 minutes to be outliers since they exceed Q 3 + 1.5•IQR where Q 3 is the third quartile and IQR is the inter-quartile range. After removing these outlier values, the range reduces to 40 minutes, which is remains stubbornly wide for a daily transportation service. In the box plot of waiting times per day in Figure 2.10, we observe that weekdays and Saturdays are similar with a median around 11 minutes, with Fridays being slightly lower. In France, Wednesday is a common day for working from home since there are also no school classes on Wednesday afternoon, this decrease in available drivers does not appear to lead to an increase in waiting times for passengers compared to other weekdays.

Chapter 2 -Quantile Regression for waiting time prediction in a stochastic carpooling network Also, Saturday mornings have school classes. Sundays, on the contrary, have significantly longer waiting times. This can be explained by the fact that on Sunday there are far fewer drivers than on other day of the week. For comparison, the average waiting time for a bus in the same region is 12 minutes according to the Global Transit Usage Report 2016 compiled by Moovit. 2 The comparison with public bus waiting times is more pertinent than with than other on-demand carpooling services such as Uber, since Ecov's carpooling service more closely resembles the former than the latter.

Explanatory variables

The explanatory variables or covariates employed to predict the waiting times are divided into three categories: (i) the geographical information, (ii) the temporal information, and (iii) the passenger information associated with the carpooling request. This information is necessary to match a passenger to a driver within spatio-temporal constraints. At this early experimental phase of the carpooling service, no direct information from the drivers is collected (we will see in the following chapters the effect of the availability this fourth category of auxiliary information from the drivers). Since the passenger has the freedom to indicate any location in the service area as a carpooling destination, there are 252 separate named destinations observed from the approximately 700 trips. The origins are always one of the 17 physical carpooling meeting points, as shown in Figure 2.2. We reduce the number of possible origin-destination pairs by performing a spatial k-means clustering based on the latitude and longitude of the destinations. Figure 2.11 shows the map with these seven clusters which represent the seven possible destination neighbourhoods. This clustering reduces the number of observed origin-destination pairs from 252 to 66. The reduction in the number of origin-destination pairs increases the number of carpooling trips within each of these pairs. Figure 2.12 shows the number of observations for the five most frequented origin-destination pairs before (left) and after (right) clustering. Whilst these are not the same five origin-destination pairs, we are still able to observe that, after clustering, we have around twice as many carpooling trips per origin-destination pair which will assist in the robustifying the regression predictions. Since driver information is not collected directly, there is no observed information about the distance and the duration of the carpooling trip. Since these values may have an important influence as explanatory variables in the regression models, then we approximate them using route-finding APIs (Application Programming Interfaces) via the latitude and longitude of the origins and the destinations.

The temporal information is derived from the timestamps with the carpooling request. From this information we are able to derive new explanatory variables, namely the hour of the day and the day of the week, and subsequently, the total number of carpooling requests per hour and per day of the week. These aggregated values add information about the overall trends in the driver population into the regression models, as shown in Figure 2.13. From the left plot, we see that the distribution of carpooling requests within a day is a bi-modal distribution with peaks in the morning and evening peak hours. So we add an indicator variable for the morning (6h00 -9h00) and evening peak hours (16h00 -19h00). From the right plot, there are fewer carpooling requests on Sundays and this leads us to add an indicator variable for holiday days network The passenger information is collected by Ecov via a standard passenger registration form. As a passenger uses the service over time, we obtain more information about their carpooling habits, such as the duration of a passenger's registration and the number of carpooling requests made. These measure the familiarity of the passenger with Ecov's carpooling services: internal studies conducted by Ecov have shown that a passenger who is familiar with the Ecov system is likely to depart more quickly than a first-time passenger. Whilst it is also possible to analyse more of the passenger's personal data (gender, age, etc), in order to comply with the privacy concerns in force in the European Union, namely the General Data Protection Regulation (GDPR), we chose not to do so.

Table 2.1 is a summary of all the geographical, temporal and passenger exploratory variables in the regression models to predict waiting times. The first column is the name of the variable, the second is a description and the third is the category (geographical, temporal or passenger).

Quantile Regression models comparison

We provide a comparison of the performance of the four regression algorithms presented in the previous sections, namely the linear Quantile Regression (QuantReg), the Quantile Regression Forests (QRF), the Generalized Random Forests (GRF) and Quantile Regression with Gradient Boosting (GBM), on these carpooling trips made in the network. The quantiles chosen for the comparison are the lower quartile 25% and the upper quartile 75%. We seek to test the models on quantiles of the distribution which are different from the median, though we do not focus on extreme quantiles, e.g. 5% or 95%. These extreme quantiles are difficult to estimate and require large amounts of data for robust prediction, which we do not have for the Ecov carpooling trips. See De Haan et al., 2007 for a complete overview of the analysis of extreme quantiles. Table 2.2 shows the Pinball losses for the four competing models for all test folds of the 10-fold cross validation. The upper rows corresponds to the accuracy of the Quantile Regressions for the 75% quantile and the lower rows for the 25% quantile. Overall, the GRF model has uniformly the lowest Pinball losses (around 3-4) for all test folds, compared to the three other models (around 5-9)

Chapter 2 -Quantile Regression for waiting time prediction in a stochastic carpooling network for both quantiles. This is confirmed when considering the overall generalised pinball loss, with the performance gain of the GRF more prominent for the 75% quantile than with the 25% quantile. In addition, the variability of all loss functions over all test folds is low, indicating that the predictions of the waiting times for all the regressions models are stable. Table 2.3 shows the Ramp losses for the four competing models for all test folds of the 10-fold cross validation. The upper rows corresponds to the accuracy of the Quantile Regressions for the 75% quantile and the lower rows for the 25% quantile. In comparison to the Pinball losses, the Ramp losses have larger variability between each of the 10 test folds within the same Regression model. So there is no uniformly best model in this respect. If we examine the overall generalised loss, then all models seem to be centred towards the optimal value 1 -α, indicating that all models overall are predicting the target quantile accurately, even though there is a large variability for the individual test folds. The model with the loss closest to 1 -α is QuantReg for α = 0.25, and GRF for α = 0.75.

α = 0.
Taking into account both the Pinball and Ramp losses, the GRF model seems to have most stable and accurate results of these considered Regression models. However this is only a relative comparison, since if we visually inspect the predicted quantiles of the waiting times for α = 0.25, 0.75 in Figure 2.14, not even the GRF has a sufficient accuracy in these waiting time predictions to warrant an operational deployment.

Conclusion

Quantile Regression is an interesting choice for predicting waiting times in a carpooling service, as it provides more detailed information than the mean waiting times provided by classical regression. The business objective was to harness this more detailed information to tailor the carpooling service to the customer behaviour and expectations. However, due to the complexity of the carpooling service as offered by Ecov during this early experimental phase, this business objective was not realised with this regression analysis, despite the strong mathematical framework underlying the quantile regression models. The operational complexity induced by allowing the passengers to select any arbitrary location in the service as a destination for a carpooling was not able to provide a reliable supply of carpooling drivers to their passengers. This creates a negative feedback where fewer passengers then continue to utilise the carpooling service. This ultimately leads to an insufficient number of completed carpooling trips (out of the total number of car trips made in the service area) required, as well as a lack of directly observed data of driver behaviour, for robust quantile regression analysis of the carpooling origin-destination pairs. From these lessons learnt in this early experimental phase, Ecov then pivoted their carpooling service to service a limited number of destinations, more similar to bus routes rather than Uber-like services, to aggregate sufficiently the passenger demand and driver supply to ensure a reliable matching between these two groups. Ecov also subsequently incorporated technology to collect the GPS traces of driver itineraries from their mobile phone application to better understand empirical driver behaviour. In the following chapters we network will analyse the data from these improved carpooling services and how they facilitate a better understanding of the carpooling potential in the service areas.

Appendix : Details of the supporting mathematical results

Quantiles as minimisers of the Pinball loss

Lemma 1 is proved by asserting that is equivalent to the following definition.

Definition 1. Let Y be a real-valued integrable random variable whose distribution F Y is strictly increasing. Then the quantile q α satisfies i

. ∀α ∈ ]0; 1[, F Y (q α ) ≥ α
ii. if q < q α then F Y (q) < α.

Appendix : Details of the supporting mathematical results

Proof of Definition 1. This proof is also given in [START_REF] Biau | Sequential quantile prediction of time series[END_REF]. We proceed by separating two cases, when q ≥ q α and q < q α . For q ≥ q α ,

E[ρ α (Y -q)] -E[ρ α (Y -q α )] = E[(Y -q)(α -1 [Y ≤q] )] -E[(Y -q α )(α -1 [Y ≤qα] )] = E[(Y -q)(α -(1 [Y ≤qα] + 1 [qα<Y ≤q ])) -(Y -q α )(α -1 [Y ≤qα] )] = E[(q α -q)(α -1 [Y ≤qα] )] -E[(Y -q)1 [qα<Y ≤q ].
Thus, the first term of right part

E[(q α -q)(α -1 [Y ≤qα] )] = (q α -q)E[α -1 [Y ≤qα] ] = (q α -q)(α -P[Y ≤ q α ]) = (q α -q)(α -F Y (q α )) ≥ 0.
Note that by definition q α -q ≤ 0 and by the Property 1 we have F Y (q α ) ≥ α.

For the second term of the right part, indicator

1 [qα<Y ≤q] implies that Y -q ≤ 0 and we obtain that -E[(Y -q)1 [qα<Y ≤q ] ≥ 0. For q < q α , E[ρ α (Y -q)] -E[ρ α (Y -q α )] = E[(Y -q)(α -1 [Y ≤q] )] -E[(Y -q α )(α -1 [Y ≤qα] )] = E[(Y -q)(α -1 [Y ≤q] ) -(Y -q α )(α -(1 [Y ≤q] + 1 [q<Y ≤qα] ))] = E[(q α -q)(α -1 [Y ≤q )] -E[(Y -q α )(α -1 [q<Y ≤qα] )] = (q α -q)(α -P[Y ≤ q]) -E[(Y -q α )(α -1 [q<Y ≤qα] )]
The first term of right part (q α -q) > 0 is strictly positive by definition and

α -P[Y ≤ q] = α -F Y (q) > 0.
We obtain Property 1 that F Y (q) < α. For the second term of the right part we have directly that -E[(Y -q α )(α -1 [q<Y ≤qα] )] ≥ 0. Thus, we obtain the desired result, i.e. for all q ∈ R, we have

E[ρ α (Y -q)] - E[ρ α (Y -q α )] ≥ 0.

Conditional quantiles as minimisers of a Pinball loss

Definition 2. Let Y be real-valued random variable whose conditional distribution

F Y |X X X is almost surely increasing. Then the conditional quantile q α (Y |X X X) is a random variable that satisfies i. q α (Y |X X X) is an element of L 1 (X X X), ii. P(Y < q α (Y |X X X)|X X X) ≤ α ≤ P(Y ≤ q α (Y |X X X)|X X X) almost surely.
Chapter 2 -Quantile Regression for waiting time prediction in a stochastic carpooling network Definition 3. Let f : R → R be a convex function and let

x 1 , x 2 ∈ R. Then f (x 2 ) = f (x 1 ) + x 2 x 1 D + f (t)dt,
where D + f is the right-hand derivative of f .

With these two definitions, we are able to state the proof of Lemma 2.

Proof Lemma 2. A more complete proof is given in Armerin, 2014. Let x 0 ∈ R, the Pinball function is convex, hence the right-hand derivative at point

(x -x 0 ) is D + ρ α (x -x 0 ) = α -1 {x≤x 0 } .
Then the Definition 3 is applied to Pinball loss function

ρ α (x 2 -x 0 ) = ρ α (x 1 -x 0 ) + x 2 x 1 D + ρ α (t -x 0 )dt. Let x 1 = q(X X X) ∈ L 1 (X X X)
, the space of all integrable functions with respect to X X X, and so by construction

x 2 = q α (Y |X X X) ∈ L 1 (X X X) for x 0 = Y . Then ρ α (q α (Y |X X X) -Y ) = ρ α (q(X X X) -Y ) + qα(Y |X X X) q(X X X) [α -1 {t≤Y } ]dt.
It follows that

E[ρ α (q α (Y |X X X) -Y )] = E[ρ α (q(X X X) -Y )] + E[ qα(Y |X X X) q(X X X) [α -1 {t≤Y } ]dt]. Write E[ qα(Y |X X X) q(X X X) [α -1 {t≤Y } ]dt] = E[ qα(Y |X X X) q(X X X) 1 {qα(Y |X X X)≥q(X X X)} [α -1 {t≤Y } ]dt] + E[ qα(Y |X X X) q(X X X) 1 {qα(Y |X X X)<q(X X X)} [α -1 {t≤Y } ]dt].
By definition of the definite integrals which states that the integral is bounded by the supremum and infimum of the integrated function. Then

E[ qα(Y |X X X) q(X X X) [α -1 {t≤Y } ]dt] ≤ E[1 {qα(Y |X X X)≥q(X X X)} (q α (Y |X X X) -q(X X X))(α -1 {qα(Y |X X X)≤Y } )] + E[ qα(Y |X X X) q(X X X) 1 {qα(Y |X X X)<q(X X X)} [α -1 {t≤Y } ]dt] ≤ E[1 {qα(Y |X X X)≥q(X X X)} (q α (Y |X X X) -q(X X X))(α -1 {qα(Y |X X X)≤Y } )] + E[1 {qα(Y |X X X)<q(X X X)} (q(X X X) -q α (Y |X X X))(-α + 1 {qα(Y |X X X)<Y } )].
Since q(X X X), q α (Y |X X X) ∈ L 1 (X X X), then by definition of the expectation of finite variables,

E[1 {qα(Y |X X X)≥q(X X X)} (q α (Y |X X X) -q(X X X))(α -1 {qα(Y |X X X)≤Y } )] = E[1 {qα(Y |X X X)≥q(X X X)} (q α (Y |X X X) -q(X X X))]E[(α -1 {qα(Y |X X X)≤Y } )|X X X]] = E[1 {qα(Y |X X X)≥q(X X X)} (q α (Y |X X X) -q(X X X))]E[(α -1 {qα(Y |X X X)≤Y } )|X X X]] = E[1 {qα(Y |X X X)≥q(X X X)} (q α (Y |X X X) -q(X X X))](α -P(q α (Y |X X X) ≤ Y )|X X X)] ≤ 0,
and

E[1 {qα(Y |X X X)<q(X X X)} (q(X X X) -q α (Y |X X X))(-α + 1 {qα(Y |X X X)<Y } )] = E[1 {qα(Y |X X X)<q(X X X)} (q(X X X) -q α (Y |X X X))]E[(-α + 1 {qα(Y |X X X)<Y } )|X X X]] = E[1 {qα(Y |X X X)<q(X X X)} (q(X X X) -q α (Y |X X X))]E[(-α + 1 {qα(Y |X X X)<Y } )|X X X]] = E[1 {qα(Y |X X X)<q(X X X)} (q(X X X) -q α (Y |X X X))](-α + P(q α (Y |X X X) < Y |X X X)] ≤ 0.
The previous two inequalities arise from Definition 2 of the conditional quantile.

Thus the result is E[ρ α (q α (Y |X X X) -Y )] ≤ E[ρ α (q(X X X) -Y )], for any function q(X X X) ∈ L 1 (X X X).

Pinball loss and its scoring function

We demonstrate the intimate link between the scoring function ψ θ and the Pinball loss ρ α . This relationship will justify the choice to use the ψ θ as scoring function in Equation (2.8) in order to estimate quantiles of the order α. Denote D -ρ α , D + ρ α as the left and right derivatives of ρ α . If h > 0 is a sequence of positive real numbers converging to 0, then the left derivative is

D -ρ α (y -θ(x x x)) = lim h→0 ρ α (y -θ(x x x) -h) -ρ α (y -θ(x x x)) h = lim h→0 (y -θ(x x x) -h)(α -1 {y≤θ(x x x)+h} ) -(y -θ(x x x))(α -1 {y≤θ(x x x)} ) h = lim h→0 (y -θ(x x x))(α -1 {y≤θ(x x x)+h} -α + 1 {y≤θ(x x x)} ) -h(α -1 {y≤θ(x x x)+h} ) h = lim h→0 (y -θ(x x x))(1 {y≤θ(x x x)+h} + 1 {y≤θ(x x x)} ) h -(α -1 {y≤θ(x x x)+h} ) = 1 {y≤θ(x x x)} ) -α.
By symmetry for the right derivative, we obtain

D + ρ α (y -θ(x x x)) = lim h→0 ρ(y -θ(x x x) + h) -ρ(y -θ(x x x)) h = α -1 {y≤θ(x x x)} .
We note directly here the connection between ψ θ and ρ α

D + ρ α (y -θ(x x x)) = ψ θ(x x x) (y) D -ρ α (y -θ(x x x)) = -ψ θ(x x x) (y).
Chapter 2 -Quantile Regression for waiting time prediction in a stochastic carpooling network

Solving the local moment condition in Equation (2.8) with the scoring function ψ θ(x x x) (y) = α -1 {y≤θ(x x x)} thus corresponds to finding the zeros of the derivatives of the Pinball loss. That is, we verify Lemma 2 that the quantile is a minimiser of Pinball loss.

CARPOOLING SERVICE

Carpooling has the potential to transform itself into a mass transportation mode by abandoning its adherence to deterministic passenger-driver matching for door-to-door journeys, and by adopting instead stochastic matching on a network of fixed meeting points. Stochastic matching is where a passenger sends out a carpooling request at a meeting point, and then waits for the arrival of a self-selected driver who is already travelling to the requested meeting point. Crucially there is no centrally dispatched driver. Moreover, the carpooling is assured only between the meeting points, so the onus is on the passengers to travel to/from them by their own means. Thus the success of a stochastic carpooling service relies on the convergence, with minimal perturbation to their existing travel patterns, to the meeting points which are highly frequented by both passengers and drivers. Due to the innovative nature of stochastic carpooling, existing offthe-shelf workflows are largely insufficient for this purpose. To fill the gap in the market, we introduce a novel workflow, comprising of a combination of data science and GIS (Geographic Information Systems), to analyse driver GPS traces. We implement it for an operational stochastic carpooling service in south-eastern France, and we demonstrate that relaxing door-todoor matching reduces passenger waiting times. Our workflow provides additional key operational indicators, namely the driver flow maps, the driver flow temporal profiles and the driver participation rates. Keywords : Data science, Stochastic Matching, GIS, Meeting point, Network

Chapter 3 -Relaxing door-to-door matching reduces passenger waiting times: a workflow for the analysis of driver GPS traces in a stochastic carpooling service

Introduction

Carpooling has seen an explosion of utilisation in recent years [START_REF] Furuhata | Ridesharing: The state-of-the-art and future directions[END_REF]. There are many underlying reasons for this, with concerns ranging from greenhouse gas emissions and air pollution to road congestion to land use, as well as economic costs [START_REF] Shaheen | Casual carpooling in the San Francisco Bay Area: Understanding user characteristics, behaviors, and motivations[END_REF]. It also attracts intense interest since carpooling is a crucial element of almost all development plans for smart cities [START_REF] Ghoseiri | Dynamic rideshare optimized matching problem[END_REF]. A broad definition of carpooling involves a driver sharing their journey with passengers. In this paper we employ a narrower definition. We additionally require that a non-professional driver would have undertaken their journey for their own reasons, regardless of whether the passengers would have been present or not. The driver may receive payment to offset the costs of the use of their vehicle, but the profit motive is non-existent or at least not their primary motivation [START_REF] Zhu | Understanding Motivations and Impacts of Ridesharing: Three Essays on Two French Ridesharing Platforms[END_REF]). Hence we do not consider taxi-like services (such as Uber, Lyft and Kapten etc.) to be carpooling services as they employ professional drivers who create a journey in response to a passenger request and are then paid the market rate for the service rendered.

Due to the altruistic nature of carpooling, service providers tend to be small, local, non-profit organisations. Though it does not preclude viable business models arising from carpooling with non-professional drivers: the market leader BlaBlaCar levies a commission fee for facilitating the matching of drivers and passengers [START_REF] Shaheen | Online and App-Based Carpooling in France: Analyzing Users and Practices -A Study of BlaBlaCar[END_REF]. This matching is managed by a centralised platform, which we call deterministic matching since a known driver is assigned in advance to collect the passenger. This deterministic matching is highly successful for infrequent, long distance, pre-reserved carpooling journeys, as witnessed by BlaBlaCar's status as a unicorn start-up company (a market capitalisation of at least 1000 million USD). Despite the success of deterministic passenger-driver matching in this market, attempts to export it other carpooling markets have not resulted in the same level of market penetration. This is most notable for frequent, short distance journeys (from 10 to 40 km roughly), which comprise the bulk of daily home-work commutes, and so carpooling remains a marginal practice in this market.

This paper focuses on short-distance, non-reserved carpooling, and it is what we refer to when employing 'carpooling' without any qualifiers. The advent of mass carpooling depends crucially on incentivising drivers and passengers to converge onto highly frequented meeting points (hotspots) along their door-to-door journeys [START_REF] Stiglic | The benefits of meeting points in ride-sharing systems[END_REF]. This type of incentivisation is well-established for a bus network where passengers embark/disembark only at the fixed bus stops. Thus mass carpooling requires a paradigm shift from considering carpooling as an exclusively private means of transport to a closer alignment to public transport models [START_REF] Cooper | Successfully changing individual travel behavior: Applying community-based social marketing to travel choice[END_REF].

Continuing with the public transport model, the meeting points are not 3.1. Introduction defined informally between passengers and drivers, but are decided in consultation with local government authorities so that they respond to the mobility requirements in the local area, taking into account various factors such as aggregated traffic flow, socioeconomic characteristics, pedestrian accessibility, local government regulations, etc. For our purposes, we consider that the identification of the meeting hotspots has been carried out beforehand. These meeting points are then connected to each other to define carpooling lines, which have massification potential, like traditional bus lines [START_REF] Stiglic | The benefits of meeting points in ride-sharing systems[END_REF][START_REF] Li | Modeling an enhanced ridesharing system with meet points and time windows[END_REF].

Like a bus service, no pre-reservations are required, as a passenger makes an ad hoc carpooling request at a meeting point, and this request for the desired destination is communicated to all passing drivers in real-time via an electronic sign on the side of a highly frequented main road. Unlike for deterministic passenger-driver matching mentioned above, a specific driver is not assigned to the passenger by a centralised platform, but the decision to collect a passenger at the meeting point is made spontaneously by a self-selected driver. Since the actual driver who collects the passenger is not known deterministically in advance, but is only known to be drawn from the population of drivers, this is known as stochastic matching. Due to the inherent variability of these driver arrivals, stochastic matching is only feasible when employed in conjunction with a network of highly frequented meeting points.

The effects of the double innovations of fixed meeting points and stochastic matching are only sparsely covered by the recent comprehensive review of general carpooling and taxi-like services over the past two decades (H. [START_REF] Wang | Ridesourcing systems: A framework and review[END_REF]. So there are few off-the-shelf workflows which are suitable for the analysis of the data arising from a stochastic carpooling service. We introduce a data science-GIS workflow which fills this gap in the market. Its main data source is the GPS traces, and its secondary sources are the meeting point locations, the origin-destination matrices, the route finder API and the base maps. Data wrangling/geoprocessing are then applied to these data sources, with the critical geoprocessing step being the topological simplification of the GPS traces onto the carpooling network. This topological simplification is essential to be able to mutualise GPS traces which share common arrival times at the meeting points. From these simplified GPS traces, we can produce the waiting times. The latter allow us to assert that stochastic matching at meeting points leads to reduced passenger waiting times in comparison to door-to-door matching. In addition to the waiting times, other outputs from this workflow are the driver flow maps, the driver flow profiles, and the driver participation rates. These additional outputs are obtained at low marginal cost but which are important elements for evidence-based decision making in a stochastic carpooling service.

In Section 3.2 we present the theoretical reasons why door-to-door matching is insufficient to ensure a regular carpooling service. In Section 3.3 we detail Chapter 3 -Relaxing door-to-door matching reduces passenger waiting times: a workflow for the analysis of driver GPS traces in a stochastic carpooling service our data science-GIS workflow for the analysis of GPS traces. In Section 3.4, we apply this workflow to an operational stochastic carpooling service to produce the passenger waiting times and the other outputs. We end with some concluding remarks.

Door-to-door matching is an obstacle to mass carpooling

As alluded to in the introduction, door-to-door matching of complete trajectories from the origin to the destination is a structural obstacle to the transformation of carpooling to a mass transit service. To illustrate the difficulties of passengerdriver matching in space and in time for door-to-door trajectories, we can represent it with partition of a 3D cube divided into smaller sub-cubes, where the x-axis is the longitude, the y-axis the latitude and the z-axis the time, as shown in Figure 3.1. On the left, there are 9 sub-cubes, where each sub-cube represents the origin/destination of a door-to-door trajectory. The blue sub-cube in the lower left represents all the trajectories whose origins are, say, within a 5 km radius around a residential neighbourhood between 07:00 and 09:00 on Tuesday, and the green sub-cube the trajectories whose destination are within a 5 km radius of the workplace between 08:00 and 10:00 on Tuesday. So for two trajectories to match spatio-temporally in a door-to-door sense, they must share the same sub-cube for the origin, and similarly for the destination: this condition is met by the 1 pair of green and blue sub-cubes among all possible 27 pairs of sub-cubes. On the right, the conditions for a door-to-door matching are stricter, say the origin is 1 km within the residential neighbourhood during 07:00 to 07:30, and the destination is 1 km within the workplace during 08:30 to 09:00. This represents 1 pair out of 125 pairs of sub-cubes. Thus stricter door-to-door matching leads to fewer drivers being available to share their trajectories with passengers.

To supplement the heuristic observations for door-to-door matching in Figure 3.1, we demonstrate that the probability that two users (i.e. a driver and a passenger) share the same origin and destination at the same time decreases rapidly as the spatio-temporal matching conditions become more stringent. For the sake of simplicity, we suppose that the origin and destination for a driver and a passenger are both represented by independent random variables which are uniform over all sub-cubes in Figure 3.1. Let U d O and U d D be the origin and destination of a driver, and likewise U p O , U p D for a passenger. These quantities are all uniform random variables U({1, . . . , n}) where n is the number of subcubes in Figure 3.1. Then the probability of a door-to-door match between the driver and passenger is p

(n) = P(U p O = U d O , U p D = U d D )
. Since an exact formula for this probability is difficult to obtain, we approximate it by a Monte Carlo

p(n) = 1 1000 1000 i=1 1 1 1{U p O,i = U d O,i , U p D,i = U d D,i }
where 1 1 1{•} is the indicator function. Figure 3.2 is the graph of the number of sub-cubes n versus the approximate probability of a door-to-door match p(n). If there is only 1 sub-cube (i.e. no spatio-temporal constraints) the probability of a match is 1. This probabilistic certainty decreases rapidly as the spatio-temporal constraints are added: for 27 sub-cubes, this probability is 0.6, and for 125 sub-cubes, it falls to 0.2. The previous analysis was based on the synthetic uniformly distributed origins and destinations. For a more realistic example, we analyse some data generated by an operational carpooling service. Our example is the 'Lane' carpooling service (lanemove.com) operated by Ecov, in conjunction with Instant System (instant-system.com), since May 2018 in the peri-urban regions around Lyon in south-eastern France. Our main data source is the GPS traces of drivers, which can be considered to be a form of crowd-sourced data collection [START_REF] Lee | Crowd-sourced carpool recommendation based on simple and efficient trajectory grouping[END_REF]. Passenger GPS traces are more difficult to obtain, and as we are not able to replicate exactly the synthetic example of passenger-driver matching above, so we use door-to-door matching of driver GPS traces to illustrate the diminishing probabilities. Since these GPS traces provide highly detailed spatiotemporal information, we are able to determine the number of strict door-to-workflow for the analysis of driver GPS traces in a stochastic carpooling service door matches which also pass by two meeting points, as well as the number of matches when door-to-door matching is relaxed. For an illustrative example in Figure 3. A hierarchical clustering with complete linkage was carried out on the spatial locations of these origins and destinations. The dissimilarity matrix used for this hierarchical clustering is composed of the Euclidean distance between the 4-vector comprising the (origin longitude, origin latitude, destination longitude, destination latitude) of each trajectory. This dissimilarity takes into account both the origin and the destination, but not the intermediate GPS points as these actual route taken is not critical for our purposes. We cut the dendrogram at h = 6000 to yield 9 spatial clusters. These clusters are represented with the different colours. So GPS traces with the same colour can be considered as door-to-door matches with each other.

The number of GPS traces per cluster is given in Table 3.1: as cluster 1 contains 75% of the mutualisable traces, this leaves the other 25% spread sparsely over the other 8 clusters, fragmenting the supply of the carpooling trajectories to passengers.

To quantify the augmentation of the carpooling potential by relaxing door-todoor matching, we compare the door-to-door cluster with the largest cardinality (76 traces) from Table 3.1 to the number of the trajectories (121 traces) which coincide with this carpooling line regardless of their true origin and destination. Table 3.1 -Spatio-temporal door-to-door matching fragments the number of mutualisable trajectories in the Bourgoin > St-Priest carpooling line, during its morning operating hours 06:30-09:00, from 2018-11-25 to 2018-12-01. The first line is the door-to-door cluster label and the second line is the number of traces in each cluster.

These counts are an empirical equivalent of Figure 3.1: the left corresponds to the 121 meeting point (i.e. relaxed door-to-door) matches, whereas the right the 76 door-to-door matches. Thus meeting point matching represents an increase of 45 traces or 59% of the carpooling driver potential due to relaxing door-to-door matching. Furthermore, [START_REF] Stiglic | The benefits of meeting points in ride-sharing systems[END_REF] and X. [START_REF] Li | Modeling an enhanced ridesharing system with meet points and time windows[END_REF] provide more complex synthetic models to affirm that meeting points are essential to the feasibility of the mass carpooling services, and assert that it is almost impossible for a carpooling service to be based on door-to-door spatio-temporal matching.

Whilst these examples demonstrate that incentivising drivers to converge to meeting points, rather than relying on door-to-door matching, increases the potential pool of mutualisable journeys, we have not yet demonstrated that Chapter 3 -Relaxing door-to-door matching reduces passenger waiting times: a workflow for the analysis of driver GPS traces in a stochastic carpooling service this leads to reduced waiting times. This would be straightforward for the synthetic examples but this is not the case for empirically observed drivers and passengers journeys. In the next section we introduce a general workflow which indeed allows us to confirm these reduced waiting times for empirical GPS traces.

Data science-GIS workflow for the analysis of GPS traces

The GPS traces analysis workflow is illustrated in Figure 3.4. The left rectangle of Figure 3.4 contains the main data sources: the GPS traces, the meeting point locations, the origin-destination matrices, the route finder API and the base maps. The first two are supplied in-house by the carpooling service provider, the origin-destination matrices are usually supplied by a third party which has carried out a mobility survey (e.g. a national statistical agency INSEE (2018)), the route finder API is provided by a GPS navigation operator (e.g. TomTom ( 2019)), and the base maps are accessed from a cartography provider (e.g. OpenStreetMap contributors ( 2019)). There are specialised data wrangling techniques specific to spatial databases, known collectively as geoprocessing, and these are carried out, in conjunction with traditional data wrangling, in the central rectangle. The critical geoprocessing task concerns the topological simplification of the GPS traces onto the carpooling network. Whilst GPS traces are a rich source of information of driver behaviour, they are voluminous and complex. Our approach is based on network analysis tools [START_REF] Guidotti | Never drive alone: boosting carpooling with network analysis[END_REF] and complexity reduction/harmonisation algorithms [START_REF] Douglas | Algorithms for the Reduction of the Number of Points Required to Represent a Digitized Line or its Caricature[END_REF]. This topological simplification is essential to be able to mutualise GPS traces which share common arrival times at the carpooling meeting points. Once these GPS traces are in a suitable format, we are able to produce the required outputs in the right rectangle, namely the predicted waiting times, the driver flow maps, the driver flow temporal profiles and the driver participation rates.

Data sources

Our primary data source are the driver GPS traces. A GPS trace is represented by an -sequence of triplets X X X = {(X i , Y i , T i )} i=1 where (X i , Y i ) are the longitude, latitude coordinates of the GPS sensor at the i th timestamp T i . We have n GPS traces X X X 1 , . . . , X X X n in the data collection period. The m meeting point locations are represented by their GPS coordinates M M M 1 , . . . , M M M m . The origin-destination matrix is such that its (j, k) th entry is the number of journeys from j th origin to the k th destination. In addition to the origin-destination matrix, we have the GPS coordinates of the origins and the destinations. Whilst it is common that they coincide, this is not required for our workflow. The base maps are graphics files of maps of the study area, which facilitate the fast and accurate map rendering at any desired scale.

Data wrangling/geoprocessing

From the m meeting points M M M 1 , . . . , M M M m , a directed graph is constructed where the meeting points are the nodes, and an edge is drawn between the two nodes if carpooling between these two meeting points is guaranteed by the service provider. Thus a carpooling line is represented by an acyclic sub-graph with at least two nodes.

The crucial data wrangling/geoprocessing process applied to the GPS traces is the topological simplification of GPS traces on a carpooling line. Around each of the m meeting points, a buffer zone of 1 km radius is drawn to obtain B(M M M 1 ), . . . , B(M M M m ). The intersection of the buffer zones and the GPS trace, X X X ∩ B(M M M 1 ), . . . , X X X ∩ B(M M M m ), is m sub-sequences of the GPS points of X X X. For those meeting points with non-empty intersections, we consider that the driver is able to collect a passenger at these points without onerous detours.

This spatial intersection only considers the spatial proximity of the driver to a passenger at a meeting point. For the carpooling to succeed, they also need to be in temporal proximity. Among the spatial intersections X X X ∩ B(M M M 1 ), . . . , X X X ∩ B(M M M m ), we examine the corresponding timestamps and retain only those in a suitably restrained time interval. If this reduced set of spatio-temporal intersections is non-empty then we proceed to the last data wrangling/geoprocessing step.

We compute the closest GPS points in X X X to the meeting points M M M j , as defined by

X X X M M M j = {(X k , Y k , T k ) : k = argmin 1≤i≤ ((X i , Y i ) -M M M j }, j = 1, . . . , m. From this closest point X X X M M M j ,
we can extract the corresponding timestamp to be an estimate of the arrival time at the meeting point M M M j . As an example, if the Chapter 3 -Relaxing door-to-door matching reduces passenger waiting times: a workflow for the analysis of driver GPS traces in a stochastic carpooling service meeting points M M M 1 , M M M 2 form the carpooling line {M M M 1 > M M M 2 }, and if the GPS trace X X X has well-defined estimated arrival times at M M M 1 and M M M 2 , then we are able to reduce the complexity of the GPS trace. That is the points of X X X can be reduced to the sequence of 4 points

X X X(M M M 1 , M M M 2 ) = {(X 1 , Y 1 , T 1 ) > X X X M M M 1 > X X X M M M 2 > (X , Y , T )}
where (X 1 , Y 1 , T 1 ) is the driver origin and (X , Y , T ) is the driver destination. With this simplified trace X X X(M M M 1 , M M M 2 ), we are still able to determine if the driver can fulfil a passenger request at a given time on the carpooling line {M M M 1 > M M M 2 }. The complex topology of X X X is thus simplified by retaining a small number of key derived indicators [START_REF] Lee | Crowd-sourced carpool recommendation based on simple and efficient trajectory grouping[END_REF].

We repeat these data wrangling/geoprocessing steps for all n GPS traces. The result is a reduced set of ñ GPS traces which correspond to the driver journeys which closely resemble the spatio-temporal characteristics of the likely passenger requests along the carpooling line.

Outputs

For the first output in the workflow in Figure 3.4, if we visualise the GPS traces of the reduced set of ñ meeting point matches with the base maps, then we obtain a map of the driver flow that matches to the passengers in the carpooling line, as in Figure 3.3. For the second output in the workflow, suppose that the initial time interval of interest is divided into n T sub-intervals τ j , j = 1, . . . , n T since we wish to quantify the driver flow at a higher temporal resolution. Computing the driver flows f (τ j ), j = 1, . . . , n T , is straightforward as it only requires an enumeration of the simplified GPS traces whose estimated arrival times fall within each sub-interval τ j . That is, the driver flow for the carpooling line

{M M M 1 > M M M 2 } during the time interval τ j is f (τ j , M M M 1 , M M M 2 ) = #{i : X X X i (M M M 1 , M M M 2 ) ∈ τ j , i = 1, . . . , ñ}. (3.1)
For the third output in the workflow, let W (t) be the waiting time until the driver arrival for a carpool request made at time t. For stochastic carpooling, since a specific driver is not dispatched to the given passenger, the problem is equivalent to the arrival time of the first driver from the population of available drivers. Assuming a Poissonian driver arrival process, the waiting time and the driver flow are inversely proportional to each other, W (t) ∝ len(τ j )/f ( τ j ) where t ∈ τ j and len(τ j ) is the length of the time interval τ j . For simplicity, we set the constant of proportionality to 1 as this corresponds to the assumption that all geolocated drivers are willing to respond to a carpooling request. It is a reasonable assumption that the vast majority of geolocated drivers are willing to pick up a passenger, according to unpublished evidence supplied by Ecov.

Data science-GIS workflow for the analysis of GPS traces

Thus for the carpooling line {M M M 1 > M M M 2 }, the passenger waiting times for the time interval τ j , j = 1, . . . , n T , are

W (τ j , M M M 1 , M M M 2 ) = len(τ j )/f (τ j , M M M 1 , M M M 2 ).
(3.2)

For the fourth output in the workflow, the driver participation rate is P = n 1 /n 0 where n 1 is the total number of the drivers who are motivated to carpool in response to a passenger request, and n 0 is the total numbers of drivers who undertake journeys in the same geographical region as the carpooling service.

Both n 1 and n 0 are difficult to define and to estimate precisely. We propose that ñ, calculated above as the number of drivers who share their geolocation, to be our proxy for n 1 , as the vast majority of carpooling journeys are assured by drivers who are willing to share their geolocation.

To enumerate all n 0 drivers in the same geographical region as the carpooling network is difficult since the GPS traces for all drivers are not available. Our proxy (ñ 0 ) is derived from inferring likely trajectories from the reference origindestination matrix. Usually this origin-destination matrix is provided at the county-level, but this is insufficiently detailed to decide if the drivers match with the meeting points on the carpooling lines. So we infer likely trajectories. These inferred likely trajectories are determined as the fastest route from the origins (county centroids) to the destinations (county centroids) by a route finder API. We employ a route finder API rather than an explicit model-based methodology, e.g. [START_REF] Tang | Estimating the most likely space-time paths, dwell times and path uncertainties from vehicle trajectory data: A time geographic method[END_REF], to infer these most likely routes. Model-based methods are the product of extensive theoretical and empirical work, and these tend to be difficult to access due to their proprietary nature. They also tend to be limited to dense urban regions, which are not the target regions for stochastic carpooling. Thus ñ0 is the sum of the driver flow from all origin-destination pairs whose likely trajectories coincide with the carpooling lines. The driver participation rate for a carpooling line

{M M M 1 > M M M 2 } is P = ñ/ñ 0 (3.3) where ñ = n T j=1 f (τ j , M M M 1 , M M M 2 ) from Equation (3.1
). Since there is no comparable door-to-door carpooling service operating concurrently with the meeting-point stochastic carpooling service, a direct comparison of empirical passenger waiting times is not possible. Instead, we propose an indirect comparison in three stages: (i) extract all driver GPS traces which connect two meeting points in a restrained time interval, as the meeting point matches, (ii) extract the largest hierarchical cluster of these GPS traces to serve as the door-to-door matches, and (iii) compute the driver flows using Equation (3.1) for both sets of matches, and then convert them using Equation (3.2) to passenger waiting time predictions.

In addition to the waiting times as an output, there are also the driver flow maps, the driver flow temporal profiles, and the driver participation rates. All Chapter 3 -Relaxing door-to-door matching reduces passenger waiting times: a workflow for the analysis of driver GPS traces in a stochastic carpooling service these outputs are useful in understanding the transport mix of the local area as well as the market penetration of the carpooling into the transport mix.

Case study of an operational stochastic carpooling service

Our case study focuses on the Lane carpooling service introduced earlier. Before we progress further into the data analysis of the driver GPS traces, we describe the operational details of this stochastic carpooling service. The physical meeting points require an integrated infrastructure to facilitate this real-time stochastic matching, as illustrated in Figure 3.5. The orange structure on the right functions like a bus shelter to provide protection from inclement weather whilst the passenger waits, and a prominent visual point of reference for drivers on the road. The passenger makes a carpooling request on the console (the green device with a horizontal yellow stripe). This request is displayed on the electronic sign on the roadside. In this configuration, the electronic sign is located close to the meeting point, but this can vary considerably according to the local geographical characteristics. A driver who wishes to embark the passenger in response to their request is able to do so safely in the reserved parking place.

Figure 3.5 -Configuration of a physical meeting point for the 'Lane' carpooling service. The orange structure is like a bus shelter. A passenger notifies potential drivers of their carpooling request using the console, which is then displayed on the roadside electronic sign. A driver can safely embark the passenger in the reserved parking place. Reproduced with permission from Ecov.

3.4. Case study of an operational stochastic carpooling service

Topological simplification of GPS traces on a carpooling line

The schematic diagram of the carpooling lines in the Lane network is shown on the left of Figure 3.6. The visual similarities of the schematic of this carpooling service with those associated with bus or train services is deliberately designed to induce the perception of carpooling as a form of public transport. There are 5 physical meeting points (Lyon Mermoz, St-Priest Parc Techno, Aéroport Lyon-St Exupéry, Villefontaine The Village, and Bourgoin La Grive Sortie 7), denoted by the circles with the stylised L, which function analogously to bus stops. According to mobility studies in this territory, the coloured lines connect the meeting points that have a sufficient driver flow between them to maintain a carpooling service with stochastic matching. These connected meeting points form a carpooling line, analogous to a bus line, where carpooling is only available between these meeting points. This carpooling network is represented as a directed graph, as shown on the right of Figure 3.6, where each node is a meeting point and the edge connects two nodes if they form segment of a carpooling line. For brevity, the node labels are abbreviated to the first letter, i.e. L = Lyon Mermoz, S = St-Priest Parc Techno, A = Aéroport Lyon-St Exupéry, V = Villefontaine The Village, and B = Bourgoin La Grive Sortie 7. We focus on the most frequented carpooling line, that is, the Bourgoin > St-Priest line (green line in Figure 3.6). The topology of the road network ensures that most of the journeys from Bourgoin to St-Priest pass also by the Villefontaine meeting point, that is, the B > S carpooling line includes both sub-graphs B > S and B > V > S. The period of data collection is 2019-07-25 (service launch date) to 2020-02-17 (last date for which consistent driver GPS traces are available), during the most frequented time period (the Chapter 3 -Relaxing door-to-door matching reduces passenger waiting times: a workflow for the analysis of driver GPS traces in a stochastic carpooling service morning operating hours 06:30-09:00).

A complete GPS trace X X X is displayed as the sequence of 530 blue circles in Figure 3.7. This GPS trace passes within 1 km of the B, V and S meeting point nodes, so its simplified topology consists of the 5-node sequence {origin > B > V > S > destination}, shown as the orange arrows. This simplified GPS trace represents a data compression rate of over 99% yet it retains the important information to decide the matching potential of this GPS driver trace with a passenger request on the Bourgoin > St-Priest carpooling line. In Table 3.2 is the average data compression for the ñ = 121 GPS traces on the Bourgoin > St-Priest line. The first column is the average number of GPS points in the complete driver traces #X X X, the second is the average number of GPS points of the simplified topologies # X X X, and the last column is the average data compression rate (1 -# X X X/#X X X).

Driver flow estimation

Of the ñ = 121 GPS traces that follow the Bourgoin > St-Priest carpooling line, 31 GPS traces have an arrival time at Bourgoin within 08:00 to 08:30, i.e., a close spatio-temporal match for a passenger request for a departure at the Bourgoin meeting point between 08:00 and 08:30 am, with a destination at the Table 3.2 -Data compression rate on the Bourgoin > St-Priest carpooling line, during the morning operating hours 06:30-09:00, for all driver GPS traces from 2019-07-25 to 2020-02-17. The first column is the average number of GPS points in the complete driver traces, the second is the average number of points of the simplified GPS traces, and the third column is the average data compression rate.

St-Priest meeting point. Of these 31 GPS traces, 17 of them are door-to-door matches (as defined as belonging to the largest door-to-door cluster of GPS traces in Table 3.1). These 17 traces are both door-to-door and meeting point matches and their simplified traces are displayed in Figure 3.8 as the orange diamonds/arrows. The simplified traces of the remaining 14 meeting point but not door-to-door matches are the blue diamonds/arrows. The latter represent an 82% increase in the number of drivers (from 17 to 31) who can potentially respond to a passenger carpooling request on the Bourgoin > St-Priest line. Chapter 3 -Relaxing door-to-door matching reduces passenger waiting times: a workflow for the analysis of driver GPS traces in a stochastic carpooling service Table 3.3 summarises the weekly evolution of the impact of meeting point matching over door-to-door matching. The first set of three columns focus on the entire morning operating hours (06:30-09:00) whereas the second set on the single 30 minute period (08:00-08:30) as this latter restricted period is a more realistic time frame that potential passengers are willing to wait for a driver to arrive. The first column contains the weekly aggregate number of door-to-door matches, the second the number of meeting point matches, and the third is the percentage increase due to meeting point matches, i.e. (# meeting point matches-# door-to-door matches)/# door-to-door matches. The number of door-to-door matches are enumerated from a similar hierarchical clustering to that in Table 3.1, and the number of meeting point matches are computed from Equation (3.1). This table demonstrates that the increase in the driver flow due to meeting point matching is maintained over the entire data collection period.

The simplified GPS traces also allow us to compute the driver flows for narrower time intervals than the 2.5 hour and 0.5 hour intervals in Table 3.3. Following Smith et al. (1997) and [START_REF] Mcshane | Traffic Engineering[END_REF] that 15 minutes intervals are a suitable choice because the variation in driver flows for shorter intervals is less stable, Table 3.4 displays the average driver flow for 15 minute intervals during 06:30 to 09:00. For robustness, we aggregate these driver flows over all weeks in the data collection period in applying Equation (3.1) since we have increased the intra-day temporal resolution.

Waiting time prediction

It is straightforward to convert these average driver flows in Table 3.4 into predicted waiting times using Equation (3.2). Suppose that a passenger makes a carpool request at 08:10 at the Bourgoin meeting point to travel to St-Priest. The expected waiting time is the length of the interval divided by the average driver flow in the interval 08:00-08:15, i.e. 7.5 minutes. Given that we have already established a highly detailed spatio-temporal profile of the average driver flow for a carpooling line in Table 3.4, then the predicted waiting times at the same temporal resolution are shown in Table 3.5.

For the Bourgoin > St-Priest carpooling line from 2019-07-25 to 2020-02-17, we observed roughly 1500 carpooling requests with a reliably recorded waiting time. Each box plot in Figure 3.9 displays the observed waiting times each 15 minute interval during the morning opening hours with at least one observed waiting time.

The accuracy of these predicted waiting times with respect to these observed ones is illustrated in Figure 3.10. Each box plot displays the RMSE (Root Mean Squared Error) between the observed and the predicted waiting times (from Table 3.5) for each 15 minute interval. For all 15 minute intervals, the median The first columns contains the number of door-to-door matches, the second the number of meeting point matches, and the third the percentage increase due to the meeting point matches.

RMSE is around 2 to 4 minutes which implies that the predicted waiting times fairly closely track the observed waiting times. This gives us confidence in our method for predicting waiting times in a stochastic carpooling service.

Since there is no comparable operational door-to-door carpooling service to Bourgoin > St-Priest stochastic carpooling line, then it is not possible to compare empirical waiting times from each type of carpooling. Our proxy is to compare the predicted waiting times for the door-to-door matches and the Chapter 3 -Relaxing door-to-door matching reduces passenger waiting times: a workflow for the analysis of driver GPS traces in a stochastic carpooling service meeting point matches. Since our method for waiting time prediction is fairly accurate for meeting point matches according to Figure 3.10, we reason that it will also yield accurate waiting times for door-to-door matches. In Table 3.6 are the predicted waiting times based on the weekly aggregate driver flows from Table 3.3 via Equation (3.2) for both door-to-door and meeting point matches. For all weeks, we observe a decrease in the predicted passenger waiting times.

From anecdotal evidence from Ecov, 15 minutes corresponds roughly to the maximum time that passengers are willing to wait for a driver to arrive since a pre-arranged meeting time has not made. This 15 minute threshold is exceeded by the door-to-door waiting times for most weeks, whereas the meeting point

Driver participation rate estimation

A key question for the service provider is what driver participation rate leads to passenger waiting times around 5 to 10 minutes, as observed in Figure 3.9? To respond to this question, we first need to enumerate the population of all drivers on a carpooling line. The county-level origin-destination matrix of home-work trajectories from the French official statistical agency (INSEE, 2018) is insufficiently detailed to decide if the drivers with these origins-destinations travel on the carpooling lines. So we infer likely trajectories, as determined as the fastest route by the TomTom route finder API (TomTom, 2019) starting on Tuesday 8am from the origins (county centroids) to the destinations (county centroids). A spatial intersection, similar to that carried out for the driver GPS traces, is computed to determine which trajectories pass within 1 km of the carpooling meeting points. These trajectories are shown in Figure 3.11. Note that there is no temporal information attached to this origin-destination matrix, but since they are home-work trajectories, we suppose they are effected in the morning peak hours which matches the time interval of the driver GPS traces.

If we then aggregate the corresponding driver flows in the origin-destination matrix, then we obtain ñ0 = 3821 drivers whose likely trajectories match the Bourgoin > St-Priest carpooling line. From Table 3.4, there is a daily average of ñ = 20 driver GPS traces between 06:30 and 09:30. This yields an estimated driver participation rate of P = ñ/ñ 0 = 0.52% from Equation (3.3). Even with this low driver participation rate, average passenger waiting times of 5-10 minutes are observed in Figure 3.9.

Chapter 3 -Relaxing door-to-door matching reduces passenger waiting times: a workflow for the analysis of driver GPS traces in a stochastic carpooling service The first columns contains the predicted waiting times for door-to-door matches, the second for meeting point matches, and the third the percentage decrease due to the meeting point matches.

If we were able to increase this low driver participation rate even modestly (to 1% or 5%), then the predicted passenger waiting times would fall substantially, as illustrated in Figure 3.12. In this case, these waiting times would be lower than those of bus lines and approach those of high frequency metro/subway train lines. The methods for increasing driver participation, as they lie largely outside of data science, are out of scope of this paper but are of intense interest to the service provider [START_REF] Zhu | More generous for small favour? Exploring the Role of Monetary and Pro-Social Incentives of Daily Ride Sharing Using a Field Experiment in Rural Île-de-France[END_REF][START_REF] Zhu | The limit of money in daily ridesharing: Evidence from a field experiment[END_REF][START_REF] Zhu | Understanding Motivations and Impacts of Ridesharing: Three Essays on Two French Ridesharing Platforms[END_REF].

Conclusions and future work

Stochastic real-time carpooling services differ from competing services which offer deterministic door-to-door matching for complete trajectories. Whilst the latter offer a high level of personal convenience in highly urbanised regions, door-to-door matching structurally inhibits mass adoption of carpooling, especially in peri-urban regions. Relaxing the strict door-to-door matching allows, Chapter 3 -Relaxing door-to-door matching reduces passenger waiting times: a workflow for the analysis of driver GPS traces in a stochastic carpooling service and subsequently implementing stochastic meeting point matching, allows for the mutualisation of high throughput road segments, and thus removes this obstacle. We introduced a novel data science-GIS workflow for a stochastic carpooling service. The crucial mathematical abstraction in this workflow is to reduce the complexity of driver GPS traces to a graph-based topology of the carpooling network. We illustrated this workflow on an operational stochastic carpooling service in a peri-urban region in south-eastern France. We provided quantitative justifications that the physical meeting points, by facilitating a critical mass of drivers and passengers drawn from a much larger geographical area, leads to passenger waiting times which are lower than those achieved by door-to-door matching. Our workflow is novel combination of two closely related, but historically separate, disciplines of data science and GIS into a single workflow. In addition to predicting the passenger waiting times, it is able to deliver outputs for the driver flow maps, driver flow spatio-temporal profiles, and driver participation rates. This workflow forms a solid prototype for other workflows to accompany the expansion of stochastic carpooling services to address the mobility requirements in neglected peri-urban regions in the future.

Introduction

Providing ecologically sustainable transportation that is accessible for all is one of the major challenges in the transition to post-carbon societies. A key component of the solution is carpooling services which cater to the mobility Chapter 4 -Bayesian hierarchical models for the prediction of the driver flow and passenger waiting times in a stochastic carpooling service requirements in marginalised peri-urban regions with sparser population and physical/digital infrastructure. These carpooling lines, which closely resemble traditional bus lines, connect the physical meeting points for drivers and passengers. The meeting points are placed strategically in highly frequented areas, which take into account various factors such as aggregated traffic flow, socioeconomic characteristics, pedestrian accessibility, local government regulations, etc. This concentrates the demand and the supply of carpooling so that they can reach a critical mass more quickly and more sustainably. These meeting points are where the passenger makes a carpooling request on an electronic console. Since a driver is not allocated in advance, this request is then displayed on a electronic sign on the roadside which informs all passing drivers of a passenger request to the specified destination. This is a real-time, stochastic matching between a passenger and a flow of potential drivers. This new type of passenger-driver matching, along with the aggregating effects of the highly frequented physical meeting points, enables carpooling to reach economical feasibility in peri-urban regions.

From a mathematical and technological point-of-view, it is vastly more difficult to provide a reliable waiting time of a driver arrival in stochastic matching than in deterministic matching. In the latter, a reliable waiting time requires only the tracking of a single assigned driver, whereas stochastic matching requires a more comprehensive understanding of the general driver flow. To assist in the construction of this understanding, the service operator can incentivise drivers to share their GPS locations in real-time. We focus on how to predict the driver flows and passenger waiting times from these GPS traces, which can be considered to be a form of crowd-sourced data collection [START_REF] Lee | Crowd-sourced carpool recommendation based on simple and efficient trajectory grouping[END_REF]. Due to the novelty of stochastic carpooling, there is a scarcity of empirically verified predictive models, apart from simple frequentist approaches [START_REF] Ray | Planning a real-time ridesharing network: critical mass and role of transfers[END_REF][START_REF] Papoutsis | Relaxing door-to-door matching reduces passenger waiting times: A workflow for the analysis of driver GPS traces in a stochastic carpooling service[END_REF]. We propose a more sophisticated Bayesian hierarchical approach where we first build predictive models of the potential driver flow from the observed GPS driver traces, and which are subsequently employed to predict the passenger waiting times. At the time when a passenger request is made, we model this instantaneous driver flow as a moving average of previous driver flows. Then we model the passenger waiting time as a regression with covariates based on the driver traffic flow from the first stage. Our objective is to construct a Bayesian two-stage hierarchical model which is able to predict well the daily driver flows and the hourly passenger waiting times.

The empirical data in this paper are extracted from the 'Lane' stochastic carpooling service (www.lanemove.com), operated by the carpooling provider Ecov (www.ecov.fr), in conjunction with Instant System (www.instant-system. com), in a peri-urban region in south-eastern France. See [START_REF] Papoutsis | Relaxing door-to-door matching reduces passenger waiting times: A workflow for the analysis of driver GPS traces in a stochastic carpooling service[END_REF] for more details on its set-up. The data collection period is the 382 days from 2018-05-15 (service launch) to 2019-05-31 (beginning of the following year's The passenger waiting times cover the period from 2019-07-25 to 2020-02-17. This range of dates is different to those for the driver GPS traces above since, due to operational technical difficulties, the passenger waiting times were not reliably recorded from 2018-05-15 until 2019-10-21, so these dates are excluded from the analysis. In Figure 4.2 are the (approximately) 1500 observed passenger waiting times aggregated for each day. The Lane service is guaranteed only for ordinary weekdays, and whilst the passengers and drivers are not prevented from using the service on other days, there are far fewer carpooling requests on school holiday weekdays and no requests on public holidays/weekends. In Section 4.2, we describe the two stages of the Bayesian hierarchical model for the daily driver traffic flow and the hourly passenger waiting times. In Section 4.3, we carry out a validation of the proposed model with simulated data. In Section 4.4 we apply it to empirical data drawn from an operational carpooling service, and compare the predictions of driver and passenger behaviour with those from frequentist and Bayesian non-hierarchical models. We end with some concluding remarks. passenger waiting times in a stochastic carpooling service

Bayesian hierarchical modelling of driver flow and passenger waiting times

As the driver flow and the passenger waiting time are fundamental quantities in transportation research, their estimation and prediction are the subject of a vast field of active research so we cite only a few references. Historically the simplest models for the driver flow are the moving window averages [START_REF] Stephanedes | Improved estimation of traffic flow for real-time control[END_REF]. More advanced methods draw from time series analysis, within a frequentist [START_REF] Ding | Traffic flow time series prediction based on statistics learning theory[END_REF] or a Bayesian framework [START_REF] Ghosh | Bayesian time-series model for short-term traffic flow forecasting[END_REF] have been posited. Established methods for waiting time prediction for stochastic carpooling tend to be frequentist approaches [START_REF] Ray | Planning a real-time ridesharing network: critical mass and role of transfers[END_REF][START_REF] Papoutsis | Relaxing door-to-door matching reduces passenger waiting times: A workflow for the analysis of driver GPS traces in a stochastic carpooling service[END_REF]. Due to the hierarchical relationship between the driver flows and the passenger waiting times in a stochastic carpooling service, it is natural to consider nested hierarchical models. A general introduction to hierarchical models is provided in [START_REF] Gelman | Multilevel (hierarchical) modeling: what it can and cannot do[END_REF][START_REF] Gelman | Multilevel (hierarchical) modeling: what it can and cannot do[END_REF][START_REF] Gelman | Data Analysis Using Regression and Multilevel/Hierarchical Models[END_REF]. Hierarchical models can be implemented with frequentist approaches, though we cite only industrial applications using Bayesian approaches here, e.g., image analysis learning (F.-F. [START_REF] Li | A Bayesian hierarchical model for learning natural scene categories[END_REF], football results prediction [START_REF] Baio | Bayesian hierarchical model for the prediction of football results[END_REF] and electricity load forecasting (S. [START_REF] Wang | A hierarchical Bayesian regression model for predicting summer residential electricity demand across the USA[END_REF]. Within the transport sector, examples include traffic accident prediction [START_REF] Deublein | Prediction of road accidents: a Bayesian hierarchical approach[END_REF] and traffic flow modelling [START_REF] Zammit | Bayesian hierarchical modelling of traffic flow -with application to Malta's road network[END_REF]. These latter approaches do not combine the driver traffic flow and the passenger waiting times and do not analyse data with differing time scales in the different stages in the hierarchical model, as we propose.

Our proposed Bayesian multi-level hierarchical model is composed of two nested stages, as illustrated in the flowchart in Figure 4.3. The input data (driver GPS traces) are preprocessed, as outlined in Appendix 4.7.1, so that they are suitable as input into the hierarchical models. The first model is a multi-level moving average model. It combines the robustness and simplicity of moving averages with the targeted adjustments of the multi-level coefficients [START_REF] Stephanedes | Improved estimation of traffic flow for real-time control[END_REF][START_REF] Ghosh | Bayesian time-series model for short-term traffic flow forecasting[END_REF]. The coefficient θ θ θ in this moving average model depends on the day types [START_REF] Kung | Exploring universal patterns in human home-work commuting from mobile phone data[END_REF][START_REF] Bao | Investigation of the traffic congestion during public holiday and the impact of the toll-exemption policy[END_REF],

Bayesian hierarchical modelling of driver flow and passenger waiting times

and so the number of components of θ θ θ depends on the number of different day types considered. The output from the first hierarchical model is the daily driver flow, which is the immediate input into the second hierarchical model. The latter is a Gamma regression, whose regression coefficient β β β has S components, with β s ∈ (0, 1) for s = 1, . . . , S, for each of the S time intervals of a 24-hour period. The role of β β β is to assign the daily traffic flow to these sub-daily time intervals. The output of this second hierarchical model is the temporal profile of the passenger waiting times w w w ∈ R S + for these sub-daily time intervals. The scarcity of the driver GPS traces allows us to model the driver flow robustly only at a daily level, whereas a higher temporal resolution of the output passenger waiting times is required for a carpooling service. Bayesian hierarchical models offer an intuitive treatment of these differing temporal resolutions within a single workflow.

Multi-level moving average for driver flows

From a visual inspection of the daily driver flows in Figure 4.1, a standard moving average which ignores the day types would be unable to account for the abrupt differences in the driver flow when consecutive days are of different day types. Let the day type function of day i be

DT(i) =       
ORD if day i is an ordinary workday SCH if day i is a school holiday PWE if day i is a public holiday or a weekend (4.1) Chapter 4 -Bayesian hierarchical models for the prediction of the driver flow and passenger waiting times in a stochastic carpooling service where i = 1, . . . , N . Thus a suitable K th order recurrence relation of the daily driver flow y i , for i ≥ K ≥ 1, satisfies

y i = α DT(i) K k=1 η DT(i-k) y i-k + ε i (4.2)
where α DT(•) is the coefficient for the current day i, η

DT(•) = 1{DT(•) = ORD} + η SCH 1{DT(•) = SCH}+η PWE 1{DT(•)
= PWE} are the coefficients for the past K driver flows, and {ε i } are a sequence of independent normal random variables N (0, σ 2 ε ). To ensure the identifiability of η DT(•) , without loss of generality, we set η ORD = 1 for all days.

The model in Equation ( 4.2) has a moving average structure of order K, but with two additional multi-level coefficients that make the average adaptive to the day types for the current day i and the previous K days. The multi-level coefficients η DT(•) allows us to model the current driver conditioned on the previous day types, whereas the multi-level coefficients α DT(•) re-scale these flows conditioned on the current day type. For example, if day i is a school holiday, then the right hand side of Equation (4.2) is

α SCH K k=1 [1{DT(i-k) = ORD}+η SCH 1{DT(i-k) = SCH}+η PWE 1{DT(i-k) = PWE}]y i-k . (4.3)
In the summand of Equation (4.3), the day type functions allow us to sum over the K previous days, even if they are of different types. If a previous day is a work day, then its contribution to the current driver flow is α SCH y i-k ; if a previous day is a school holiday then it is α SCH η SCH y i-k ; if a previous day is a public holiday/weekend then it is α SCH η PWE y i-k .

Our model in Equation (4.2) possesses a similar structure to an autoregressive model, though it does not strictly satisfy the definition of the latter. It cannot be defined with a back shift operator due to the action of the multi-level coefficients α DT(•) and η DT(•) , and the process {y i ∈ [0, ∞), i = 1, 2, . . .} is nonstationary due to the drift in the driver participation rate after the launch of the carpooling service.

Let θ θ θ = (α ORD , α SCH , α PWE , η SCH , η PWE , σ 2 ε ) and we fix η ORD = 1 identically, and have α ORD , α SCH , α PWE , η SCH , η PWE ∈ (0, 1) and σ 2 ε ∈ R + . Let the N days of observed daily driver flows be y i , i = 1, . . . , N , where N > K. Since the error variables are independent Gaussian random variables, then the conditional likelihood of y y y = (y K , y K+1 , . . . , y N ) is

L(y y y|θ θ θ) = 1 (2πσ 2 ε ) (N -K+1)/2 exp - 1 2σ 2 ε N i=K (y i -g i (θ θ θ)) 2
where

g i (θ θ θ) = α DT(i) K k=1 η DT(i-k) y i-k
. This conditional likelihood is formed by the product of the conditional densities of y i , given y i-K , . . . , y i-1 , for i = K + 1, . . . , N .

Bayesian hierarchical modelling of driver flow and passenger waiting times

In Bayesian analysis, the parameter of interest θ θ θ is treated as a random variable, and its prior distribution π represents our belief in its uncertainty. The posterior density π(θ θ θ|y y y) represents an update of the prior distribution by taking into account the observed data. In our case, we do not have access to existing knowledge that would provide an informative prior and thus we form a non-informative prior on θ θ θ, i.e. π(θ θ θ) ∝ σ -2 ε (Congdon, 2014, Chapter 1). This leads to the following posterior distribution

π(θ θ θ|y y y) = L(y y y|θ θ θ)π(θ θ θ)1(θ θ θ) Θ Θ L(y y y|θ θ θ)π(θ θ θ)dπ(θ θ θ) ∝ 1(θ θ θ) Θ σ N -K+1 ε exp - 1 2σ 2 ε N i=K (y i -g i (θ θ θ)) 2 (4.4)
where Θ = (0, 1) 5 × R + . For the inference on θ θ θ, Monte Carlo approximations are required since the posterior distribution (and its moments, quantiles etc.) cannot be calculated explicitly. The most widely used family of methods is the Monte Carlo Markov Chain (MCMC) which generates a Markov Chain {θ θ θ 0 , θ θ θ 1 , . . .} whose equilibrium distribution converges to the posterior distribution π(θ θ θ|y y y).

The next stage is to predict a driver flow ỹ in the future from the observed past driver flows y y y. Bayesian prediction is based on the posterior predictive distribution of ỹ|y y y. Its density p(ỹ|y y y) is given by p(ỹ|y y y) = Θ p(ỹ|θ θ θ, y y y)π(θ θ θ|y y y) dθ θ θ (4.5)

Since p(ỹ|y y y) is a compound probability distribution, we can simulate samples from this predictive distribution.

For the choice of an MCMC sampler, we use the NUT sampler [START_REF] Hoffman | The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo[END_REF], which is the default in the pyStan Python package (https: //pystan.readthedocs.io). This package is an interface to the state-of-art platform for Bayesian computations Stan (https://mc-stan.org). To carry out the integration and then a random draw from the posterior predictive distribution of daily driver flows p(ỹ|y y y) in Equation (4.4), we are only required to input the prior π(θ θ θ), the likelihood L(y y y|θ θ θ) and the recurrence relation which generates the vector of simulated driver flows y y y (Algorithm 4.2 in Appendix 4.7.2) into pyStan. The latter automatically simulates for the jth iteration, j = 1, . . . , J,

ỹ y y (j) =    
y (j,1) . . .

y (j,N )     ∼    
p(ỹ (j,1) |y y y) . . .

p(ỹ (j,N ) |y y y)     , (4.6)
and its final output is the sequence of posterior prediction vectors Ỹ = {ỹ y y (1) , . . . , ỹ y y (J) }.

Gamma regression for passenger waiting times

For simplicity, we assume that a passenger can only make one request at a time for themselves only at a carpooling meeting point, and the drivers can embark Chapter 4 -Bayesian hierarchical models for the prediction of the driver flow and passenger waiting times in a stochastic carpooling service only one passenger in their vehicle in the order that the passenger requests are made. For day i, let y i be the daily traffic flow, and that n i passengers make carpooling requests at times t i,1 < • • • < t i,n i . Let t i,j be the driver arrival time for the passenger request at time t i,j , i = 1, . . . , N and j = 1, . . . , n i . The perceived waiting time w * i,j and the pseudo waiting time w i,j for the passenger request at time t i,j are

w * i,j = t i,j -t i,j w i,j = t i,j -max(t i,j , t i,j-1 )
with the convention t i,0 = t i,1 for the first passenger on day i. Figure 4.4 illustrates the difference between the perceived and the pseudo waiting times for two passengers A, B who are both not the first passenger of the day. Passenger A arrives first and is the jth passenger of day i, and makes a carpooling request at time t i,j . Passenger B arrives immediately afterwards and is the (j + 1)th passenger with request time t i,j+1 . Suppose that there are at least two drivers en route to embark these passengers, and they have not received any passenger requests before passenger A's request. The first driver arrives at t i,j > t i,j+1 (i.e. after passenger B's request time) and the second driver at t i,j+1 > t i,j . The perceived waiting time for the passenger A is w * i,j = t i,j -t i,j (the blue brace in Figure 4.4) and for the passenger B it is w * i,j+1 = t i,j+1 -t i,j+1 (the green brace). The pseudo waiting time for passenger A is w i,j = w * i,j (the blue brace) since they are at the front of the queue, and for passenger B it is w i,j+1 = t i,j+1 -t i,j (the grey brace). The pseudo waiting time w i,j+1 for passenger B is the difference between their departure time and the departure time of the previous passenger A, and this is shorter than the perceived waiting time w * i,j+1 . From Figure 4.4, we observe that the perceived waiting times w * i,j and w * i,j+1

for passengers A and B overlap, whereas the pseudo waiting times w i,j and w i,j+1 do not overlap by construction. The overlapping nature of the interval processes that determine the perceived waiting times renders the problem of their unconditional prediction to be non-identifiable. Thus we focus on the pseudo waiting times, and we wish to formulate sub-daily predictions of them.

Let the 24 hour period of a day be divided into S equal intervals

I 1 < • • • < I S .
The fraction of the daily driver flow y i on each interval I s , s = 1, . . . , S is y i β s , where β s ≥ 0 and S s=1 β s = 1. Conditional on the traffic flow y i and the passenger request times t i,j ∈ I s , we suppose that the pseudo waiting times w i,j are independent Gamma random variables with parameters ν and β s y i , i.e.

w i,j |(y i , β β β, t i,j ∈ I s ) ∼ Γ(ν, β s y i ) (4.7)
for i = 1 . . . N and j = 1, . . . , n i . This Gamma regression model assumes that the pseudo waiting times depend on the time of day and on the daily driver flow. Furthermore, it ensures that the conditional mean pseudo waiting time is which is consistent with our intuition of the inverse relationship between the driver flow and the waiting time. Since β β β is constant for all i, then the model assumes that the relative proportions of the driver flow in the intervals I 1 , . . . , I S remain unchanged for all aggregate daily driver flows.

E[w i,j |(y i , β β β, t i,j ∈ I s )] = ν β s y i
A Dirichlet distribution is a natural choice as a prior distribution on the coefficients β β β: β β β ∼ Dir(S, α α α) where α α α = (α 1 , . . . , α S ) are the concentration parameters, since it imposes the constraint S s=1 β s = 1 on the coefficients. The corresponding Dirichlet density is p(

β β β) = S s=1 β αs-1 s /B(α α α) where B(α α α) = S s=1 Γ(α s ) Γ( S s=1 α s ) and Γ(x) = ∞ 0 u
x-1 e -u du. The β β β vector allows us to rebuild the temporal distribution of the sub-daily traffic flow from an aggregated daily driver flow.

Let t t t i = (t i,1 , . . . , t i,n i ) be the vector of the n i observed passenger carpooling request times for the day i ∈ {1, . . . , N }, and t t t = (t t t 1 , . . . , t t t N ) be all observed passenger carpooling request times. Likewise for the passenger pseudo waiting times w w w i for day i, and w w w for all days. Let y y y = (y 1 , . . . , y N ) be the observed driver flows for all days. It is reasonable to assume that the waiting times are mutually independent given (β β β, y y y, t t t). The conditional likelihood of the pseudo waiting times is thus given by the joint density of w w w given (β β β, y y y, t t t)

L(w w w|β β β, y y y, t t t)

= N i=1 p(w w w i |β β β, y y y, t t t) = N i=1 p(w w w i |β β β, y i , t t t i )
with p(w w w i |β β β, y i , t t t i ) = S s=1 {j:t i,j ∈Is} (β s y i ) ν w ν-1 i,j exp(-β s y i w i,j )/Γ(ν). Then we Chapter 4 -Bayesian hierarchical models for the prediction of the driver flow and passenger waiting times in a stochastic carpooling service obtain the posterior density of β β β, using a non-informative prior on β β β, as π(β β β|y y y, t t t, w w w)

∝ N i=1 S s=1 {j:t i,j ∈Is} (β s y i ) ν Γ(ν) w ν-1 i,j exp(-β s y i w i,j )1{β β β ∈ R S + }.
Let ws be the pseudo waiting time for a future day for a passenger who makes a carpooling request in the time interval To carry out the integration and then a random draw from the posterior predictive distribution of p( ws |ỹ, y y y, w w w) in Equation (4.8), we are only required to input the posterior predicted value of the driver flow ỹ (Equation (4.5)), the recurrence relation which generates the vector of simulated driver flows y y y (Algorithm 4.2 in Appendix 4.7.2), and the recurrence relation which generates the vector of simulated passenger pseudo waiting times w w w (Algorithm 4.3 in Appendix 4.7.2) into pyStan. The latter automatically simulates this N × S matrix distribution, for the jth iteration, j = 1, . . . , J, and its final output is the sequence of posterior prediction matrices W = { W W W (1) , . . . , W W W (J) )}.

W W W (j) =     w(j)

Model validation with simulated pseudo waiting times

We choose parameter values to produce simulated data which are comparable to those observed in the Lane carpooling service. We set the initial day i = 1 to be 2018-01-01, and the weekdays (ORD), school holidays (SCH) and public holidays/weekends (PWE) to be those observed in south-eastern France. For the simulation algorithms, the number of days is N = 365, the day types coefficients are θ θ θ = (0.333, 0.33, 0.331, 1, 1, 1), the autoregression order is K = 3, the error variance is σ 2 ε = 5, the 24 hour period is divided in S = 8 equal intervals of 3 hours, the first Gamma shape parameter is ν = 7, the Gamma regression 4.3. Model validation with simulated pseudo waiting times parameters are β β β = (0.012, 0.01, 0.011, 0.013, 0.018, 0.016, 0.017, 0.019), and the number of replicates (waiting times per time interval) is J = 10.

We generate one simulated data set of N = 365 days, each with one daily driver flow y i , i = 1, . . . , N (Algorithm 4.2), and J = 10 passenger pseudo waiting time N × S matrices W = {W W W (1) , . . . , W W W (J) } (Algorithm 4.

Wtest = { W W W (1) test , . . . , W W W (J)
test }. For brevity we have omitted the comparison of the driver flows and focus on the passenger waiting times for these simulated data: we make a more thorough comparison of both driver flows and passenger waiting times for the empirical data in the sequel.

From a passenger point of view, whilst the magnitude of the waiting times are important as a perception of the service quality, it is equally important that these posterior predicted waiting times be as close to the observed ones, whatever their magnitude. For example, suppose that a driver arrives 12 minutes after a passenger makes a carpooling request. In this case, a prediction of 15 minutes is better than 5 minutes since the former is closer to the observed waiting time than the latter (which is too optimistic). Therefore we propose the following metric to measure these discrepancies for a given threshold δ:

PE(W, W; δ) = 1 J Ñ S Ñ i=1 S s=1 J j=1 1{| wi,s -w (j) i,s | < δ} (4.9)
where wi,s = 1 J J j=1 w(j,i) s is the mean of the posterior predicted waiting times distribution for day i, and time interval I s .

We focus on the temporal profiles, over the S = 8 periods of a day, of the waiting times. In Figure 4.5 are the quantiles of the waiting times for all time intervals I s , s = 1, . . . , 8, for all days i = 1, . . . , Ñ in the test phase. The grey box plots are of the observations W test,i,s and the light, medium and dark purple circles superimposed over the box plots are the 50%, 75%, 95% quantiles of the posterior predictions Wtest,i,s . For operational purposes, short term prediction for the coming week is sufficient, and this is verified by the close agreement of the quantiles of the posterior predicted pseudo waiting times with their observed values for all Ñ = 5 prediction days. The advantage of an MCMC approach here is that we are able to reproduce the entire sampling distribution of the predicted waiting times, which is more comprehensive than point or interval predictions.

The PE metric from Equation (4.9), as a function of δ, illustrated in Figure 4.6, 4.4. Model validation with the Lane carpooling service

Model validation with the Lane carpooling service

Our objective is to employ the two-stage Bayesian hierarchical model to predict the daily driver flow distribution and the passenger pseudo waiting time distribution for the hourly intervals I s , s = 1, . . . , S, with S = 24 for the upcoming week. These are then compared to the observed driver flows and the pseudo waiting times from the same period.

Daily driver flows

We have approximately 5000 GPS traces for the 382 days from 2018-05-15 to 2019-05-31. We first apply the preprocessing, as outlined in Appendix 4.7.1, to convert the driver GPS traces into a format suitable for computing the daily driver flows y i . For the driver flow moving average model in Equation (4.2), the θ θ θ coefficient has a different value for each day type, since these day types are a key determinant of home-work daily commutes. This is verified empirically by the box plots of the daily driver flow in Figure 4.7. The daily driver flow for an ordinary weekday (ORD) approaches 150 trajectories, which is about double the driver flow on school holidays (SCH), and more than 20 times larger than on the public holiday/weekends (PWE). tions. These predictions are presented in Figure 4.9: the Bayesian hierarchical multi-level BHML in purple, the frequentist baseline model BASE in black, the Bayesian Prophet PROP in green, and as well as the observed daily driver flows in blue. The PROP predictions are mostly too low on week days and too high on weekends for all six test weeks in comparison to the observed driver flows, whilst the BHML appears to have marginally better prediction performance than the BASE.

In Figure 4.10 are the MSEs between the observed and predicted daily driver flows: the frequentist baseline model BASE in black, the Bayesian Prophet PROP in green, and the Bayesian hierarchical multi-level BHML in purple. Overall the BHML has the best prediction accuracy for all test week scenarios. PROP is the uniformly the worst of these three models for all test weeks. BASE is the best for the test scenario #1 (all ORD after PWE/SCH period) and #6 (all ORD) with almost zero prediction MSE, though the difference with BHML is not so large. These two test scenarios are where all days in the test week are the same day type. For the other test week scenarios #1, #3, #4, #6, BHML has the smallest prediction MSE, some times by a large margin. These test week scenarios include a day which is a different day type to the other days within the test week, which the BHML handles the best.

For the service operator, the sharp differences in the driver flow for different day types within the same week has operational repercussions. For example, since the driver flow is consistently low for all pubic holidays, the service operator must communicate to passengers that the service quality on a public holiday is not the same as that for an ordinary weekday. This is analogous to a public holiday schedule in lieu of a usual weekday schedule provided by a bus operator. On a more positive note, since the driver flow for weekdays on either side of the public holiday is similar to other weekdays further away, then the service operator can also communicate that this temporary reduction in service quality is limited to the public holiday itself and the usual weekday service 4.4. Model validation with the Lane carpooling service level can be assured on the preceding and following weekdays.

Temporal profiles of passenger pseudo waiting times

For the passenger pseudo waiting time Gamma regression, the β β β coefficient, which determines in the intra-day distribution of the waiting times, is considered to be constant for all days. In Figure 4.11 are the mean observed daily traffic flows for each weekday from the Lane carpooling service, where the day is divided into 15 minute intervals. Since the service operating hours are 06:00-09:00 and 16:00-19:00, there are few drivers outside them. Each dot in the figure is the mean number of drivers for each 15 minute interval for each week day from 2018-05-15 to 2019-05-31. Each week day has a similar shape so this gives some empirical justification for supposing a constant β β β for all days. For the service operator, this means that it can treat all non-public holiday weekdays as similar to each other. Before we examine the predictions from this Gamma regression model, we provide some heuristic justification of the model itself, namely concerning the choice of the Gamma distribution and the conditioning of the waiting times with respect to the driver flow. As expressed in Equation (4.7), the pseudo waiting times are represented by different Gamma distributions for each hourly interval. In Figure 4.12, we observe that each of these empirical distributions (i) resemble Gamma distributions and (ii) have a different rate parameter β s within each different hourly interval.

The other main assumption of the Gamma regression model in Equation (4.7) is that the mean pseudo waiting time is a decreasing function of the driver flow. In Figure 4.13, we have divided the observed daily driver flows into three categories: low (<100 vehicles per day), medium (100-200 vehicles per day), and high (>200 vehicles per day). Overall we observe that the mean waiting time is inversely proportional to the driver flow level. There are a maximum of around 40 observed waiting times per hourly interval per week, which are not sufficient to infer robustly their distribution within each interval. To remedy this data sparsity, we aggregate a moving window of test data so for time interval I s on day i, we combine its observed pseudo waiting times w test,i,s with those for the same time interval from the previous 5 weeks with the same day of week and same day type, i.e. {w test,i-k,s : DT(i -k) = DT(i), DN(i -k) = DN(i), k = 1, . . . , 35}. These days added to the test data are correspondingly removed from the training data. We aggregate the final 5 weeks to be a single test phase, so the Test scenario # 7 is composed of training weeks (2019-07-25 -2020-01-12) with 1289 observed training waiting times, and test weeks (2020-01-13 -2020-02-17) with 520 observed test waiting times. We make predictions for only the last test week (2020-02-10 -2020-02-17), so the number of prediction weekdays remains Ñ = 5.

In Figure 4.15 are the box plots of the observed pseudo waiting times and the quantiles for the posterior predictions, for hourly intervals for the Test scenario # 7. The observed pseudo waiting times are displayed as the grey box plots, and the 50%, 75%, 95% quantiles of the posterior predicted waiting times are the light, medium and dark purple circles. The advantage of the BHML is that we have the entire sampling distribution of the predicted waiting times, which is more comprehensive than point or interval predictions of the usual regression models. The median and upper quartile of the predicted pseudo waiting times tend to track those for the observed waiting times, especially for the 06:00-07:00, 17:00-18:00 and 18:00-19:00 intervals. From anecdotal evidence provided by Ecov, 15 minutes corresponds roughly to the maximum time that passengers are willing to wait for a driver to arrive if a pre-arranged meeting time has not been made. With the BHML predictions, we can assert that 95% of the waiting times for passenger requests do not exceed this 15 minutes threshold during most of the operating hours. Whilst this information could also be established with the empirical quantiles of the observed waiting times, the advantage of consistent waiting times can be provided to passengers. Since they are nonprofessional drivers, then non-monetary incentives are crucial in maintaining their participation in a stochastic carpooling service [START_REF] Zhu | More generous for small favour? Exploring the Role of Monetary and Pro-Social Incentives of Daily Ride Sharing Using a Field Experiment in Rural Île-de-France[END_REF][START_REF] Zhu | The limit of money in daily ridesharing: evidence from a field experiment[END_REF].

In contrast to the comparison of the predicted driver flows from the BHML to those from competing models in the previous section, the comparison of the Gamma regression with other possibilities in the second stage of the BHML is not considered here. According to [START_REF] Papoutsis | Relaxing door-to-door matching reduces passenger waiting times: A workflow for the analysis of driver GPS traces in a stochastic carpooling service[END_REF], the main predictor for the passenger waiting times is the driver flow, and so we conjecture that the choice of the passenger waiting times prediction model is of secondary importance with respect to the choice of the driver flow prediction model.

Temporal profiles for passenger perceived waiting times

We have focused on the pseudo waiting times w i,j , though the perceived waiting times w * i,j are more pertinent, since the latter are the true waiting times from the passenger point-of-view. The current set-up of the Lane carpooling service is not able to collect reliable perceived waiting times since the passenger arrivals are not reliably tracked. Also, as alluded to earlier, the unconditional prediction of the perceived waiting times is a non-identifiable problem. So a complete analysis of them is out of scope of this manuscript. Nonetheless, we can provide a framework for their future analysis based on conditioning on simulated passenger arrival processes. Since the complete posterior predictive distribution in Equation (4.8) is available, we can integrate it with respect to the passenger arrival distribution to obtain the pseudo waiting time distribution. Then we are able to reconstruct the perceived waiting times using

w * i,j+1 = w i,j+1 + [w i,j - ζ i |(w i,j > ζ i )] with ζ i = t i,j+1 -t i,j .
Let the passenger arrivals be a Poisson process, for t > 0, N (t) = max{n : n k=0 A k ≤ t} where N (0) = 0, A 0 = 0, and A k ∼ E(λ) and λ > 0 are independent exponential random variables. The parameter λ is the rate of the passenger arrivals, and it measures the mean number of arrivals over a unit of time. In our case, the unit of time is one hour and we focus on the opening hours of the service (I 1 = 6:00-7:00, I 2 = 7:00-8:00, I 3 = 8:00-9:00, I 4 = 16:00-17:00, I 5 = 17:00-18:00, I 6 = 18:00-19:00). The base passenger arrival rate is denoted by λ 1 , as shown in Table 4.2. The two other scenarios involve λ 2 , an increase by 50%, and λ 3 , an increase of 100%.

For the base passenger arrival scenario, there are few situations where the passenger requests overlap each other, and so there is little difference between the pseudo and the perceived waiting times. The second and third scenarios with increased rates of passenger arrivals lead to increased overlapping passenger requests. This is of intense interest to the service operator because it informs them how the passenger perceived waiting times respond to the increased passenger requests whilst maintaining the current driver flow. In Figure 4.17 passenger waiting times in a stochastic carpooling service

Morning intervals

Evening intervals 06:00-07:00 07:00-08:00 08:00-09:00 16:00-17:00 17:00-18:00 18:00-19:00 are the box plots for the pseudo and perceived waiting times for the three passenger arrival scenarios from Table 4.2 with a constant driver flow. The pseudo waiting times of the left panel are similar for the three scenarios since they do not account for overlapping passenger arrivals. On the other hand, we observe that the perceived waiting times on the right panel tend to increase as the number of passengers arriving increases. For λ 2 with a 50% increase in the passenger arrivals, the perceived waiting times remain acceptable for a stochastic carpooling service (median less than 15 minutes). However for λ 3 with a 100% increase in the passenger arrivals, the perceived waiting times exceed 15 minutes for many passengers. For the service operator to reduce the waiting times, the driver flow must be increased by increasing the driver participation rate.

λ

Conclusions

The main contribution of this paper is the prediction of the daily driver flows and the hourly passenger waiting times using a nested two-stage Bayesian hierarchical model. The first stage is a multi-level moving average model of the daily driver flows, where the multi-level coefficient depends on if the current day is a work day, a school holiday or a public holiday/weekend. The second stage is a Gamma regression where the covariates are the daily driver flows from the first stage, and the response variables are the hourly passenger waiting times. The predicted driver flows and passenger waiting times are robust going into the future, since we demonstrated that they are not due to over-fitting. Furthermore, since we analyse the data from an operational carpooling service, we are able to provide operational advice. For the service operator, the baseline frequentist model is the simplest to implement, and so may be sufficient under certain cost-benefit scenarios. However only the more complex BHML can be utilised for more in-depth data analysis, such as quantiles and confidence regions of passenger waiting times, and for forward planning with the effect of increased passenger requests on waiting times. We focused on modelling the driver arrival processes and assumed the passenger arrivals to be non-random in the first stage, and on pseudo waiting times in the second stage. One main advantage of the Bayesian hierarchical framework is that it is straightforward to generalise any of the models in the constituent stages (i) to allow the passenger arrivals to also be a random process, and (ii) to predict both the pseudo and perceived passenger waiting times. These perceived waiting times are of intense operational interest to stochastic carpooling service providers.

Chapter 4 -Bayesian hierarchical models for the prediction of the driver flow and passenger waiting times in a stochastic carpooling service GPS coordinates M M M 1 , . . . , M M M m . Around each of the m pick-up/drop-off locations M M M 1 , . . . , M M M m , a ball of 1 km radius is drawn to obtain B(M M M 1 ), . . . , B(M M M m ). The intersection of these balls and the GPS trace, X X X ∩ B(M M M 1 ), . . . , X X X ∩ B(M M M m ), is m sub-sequences of the GPS points of X X X. For those pick-up/drop-off locations with non-empty intersections, we consider that the driver is able to collect a passenger at these points without onerous detours.

This only considers the spatial proximity of the driver to a passenger at a pick-up/drop-off location. For the carpooling to succeed, they also need to be also in temporal proximity. Among the spatial intersections X X X ∩ B(M M M 1 ), . . . , X X X ∩ B(M M M m ), we examine the corresponding timestamps and retain only those in a suitably restrained time interval. If this set of spatio-temporal intersections is non-empty then we proceed to compute the closest GPS points in X X X to the pick-up/drop-off locations M M M j , as defined by

X X X M M M j = {(X k , Y k , T k ) : k = argmin 1≤i≤ ((X i , Y i ) -M M M j }, j = 1, . . . , m. From this closest point X X X M M M j ,
we extract the corresponding timestamp T k to be an estimate of the driver arrival time at M M M j .

As an example, suppose that there two pick-up/drop-off points M M M 1 , M M M 2 at which the GPS trace X X X has well-defined estimated arrival times. Then the points of X X X can be reduced to the sequence of 4 points X

X X = {(X 1 , Y 1 , T 1 ) > X X X M M M 1 > X X X M M M 2 > (X , Y , T )} where (X 1 , Y 1 , T 1 )
is the driver origin and (X , Y , T ) is the driver destination. With this simplified trace X X X, we are still able to determine if the driver can fulfil a passenger request at M M M 1 for a trip going to M M M 2 at time t. The complex topology of X X X is simplified by retaining a small number of key derived indicators [START_REF] Lee | Crowd-sourced carpool recommendation based on simple and efficient trajectory grouping[END_REF].

Simulation algorithms

Algorithm 4.1 simulates a driver flow for a single day with no day types. Equation (4.2) with no day types simplifies to y i = α K k=1 y i-k + ε i , which is a true autoregressive model. The inputs are the day i, the coefficient α, the autoregression order K, and the error variance σ 2 ε . The output is a single driver flow for day i. The repeat loop ensures that the simulated driver flow is strictly positive. To simulate a sequence of N driver flows, we initialise the values generated by Algorithm 4.1 for i = 1, . . . , K days, and then iterate Algorithm 4.1 sequentially for i = K + 1, . . . , N .

With Algorithm 4.1 defined, it is straightforward to define one with day types (i.e. Equation (4.2)) in Algorithm 4.2. The latter has similar inputs: the day i, the day type coefficients θ θ θ, the autoregression order K, the error variance σ 2 ε , and the vector coefficients θ θ θ. The output is the daily driver flow for day i, accounting for the day types before day i.

Algorithm 4.3 simulates the passenger pseudo waiting times in Equation (4.7) for a sequence of days. The inputs are the number of days N , the day type co-4.7. Appendix Algorithm 4.1: Daily driver flow without day types

1 procedure TRAFFICFLOW(i, α, K, σ 2 ε ) 2 if i <= K then 3 initialise y ←-N (30, σ 2 ε ) 4 else 5 repeat 6 y ←-N (α K k=1 TRAFFICFLOW(i -k, α, K, σ 2 ε ), σ 2 ε ) 7 until y > 0;
8 return: y driver flow for day i Algorithm 4.2: Daily driver flow with day types

1 procedure TRAFFICFLOWDT(i, θ θ θ, K, σ 2 ε ) 2 if DT(i) == ORD then 3 y ←-TRAFFICFLOW(i, α ORD , K, σ 2 ε ) 4 else 5 if DT(i) == SCH then 6 y ←-TRAFFICFLOW(i, α SCH η SCH , K, σ 2 ε ) 7 else 8 if DT(i) == PWE then 9 y ←-TRAFFICFLOW(i, α PWE η PWE , K, σ 2 ε )
efficients θ θ θ, the autoregression order K, the error variance σ 2 ε , the first shape parameter for the Gamma distribution ν, the S regression parameters β β β, and the number of replicates of the waiting times J. The output are J replicates of a pseudo waiting time for each time interval I s , s = 1, . . . , S, for each day i = 1, . . . , N . The TRAFFICFLOWDT procedure (Algorithm 4.2) is called outside of the replicates loop since all waiting times on a given day are simulated from the same daily driver flow.

An iteration of the nested loop in Algorithm 4.3 in the Appendix results in a single N × S matrix of pseudo waiting times drawn from the appropriate Gamma distributions

W W W (j) ∼     Γ(ν, β 1 y 1 ) . . . Γ(ν, β S y 1 ) . . . . . . Γ(ν, β 1 y N ) . . . Γ(ν, β S y N )     .
Chapter 4 -Bayesian hierarchical models for the prediction of the driver flow and passenger waiting times in a stochastic carpooling service Algorithm 4.3: Passenger pseudo waiting times

1 procedure WAITINGTIME(N, θ θ θ, K, σ 2 ε , ν, β β β, J) 2 S ←-LEN(β β β) 3 for i in 1:N do 4 Y [i] ←-TRAFFICFLOWDT(i, θ θ θ, K, σ 2 ε ) 5 for j in 1:J do 6 for i in 1:N do 7 for s in 1:S do 8 W W W (j) [i, s] ←-Γ(ν, β s Y [i])
9 return: W W W (1) , . . . , W W W (J) sequence of waiting time matrices for j = 1, . . . , J. These are collated into the sequence

W = {W W W (1) , . . . , W W W (J) } =            w (1) 1,1 . . . w (1) 1,S . . . . . . w (1) N,1 . . . w (1) N,S     , . . . ,     w (J) 1,1 . . . w (J) 1,S . . . . . . w (J) N,1 . . . w (J) N,S            .
As Equation (4.8) generates only a single posterior prediction ws for a time interval I s , we collate these ws for s = 1, . . . S into an S-vector, and in turn collate N of these S-vectors of posterior prediction distributions row-wise into a N × S matrix.

Competing models for driver flows

In addition to the multi-level moving average model for driver flows, we consider a baseline frequentist model and a Bayesian Prophet model. The baseline frequentist model has multi-levels like our model, but without the Bayesian moving average structure. To account for the school/public holidays, as proposed by [START_REF] Gould | Forecasting time series with multiple seasonal patterns[END_REF], if day i is not a school/public holiday then the average is calculated over all previous days with the same day of week as day i; and if day i is a school/public holiday, then the average is over all previous school/public holidays. That is, Notons θ ∈ R d le paramètre inconnu à estimer à partir des observations x 1 , . . . , x n . On note aussi f (x|θ) la densité de probabilité selon laquelle les observations sont échantillonnées. Contrairement aux statistiques fréquentistes, le paramètre θ est modélisé par une loi de probabilité π appelée loi a priori. Notons que la fonction f (x|θ) représente la densité de la loi des observations conditionnellement au paramètre θ. La règle de Bayes permet d'obtenir la loi du paramètre θ conditionnellement aux observations x. Il s'agit de la loi a posteriori définie par [START_REF] Cucchi | Ex-situ priors: A Bayesian hierarchical framework for defining informative prior distributions in hydrogeology[END_REF] propose un prior informatif dans le secteur de l'hydrogéologie. Cette fois le prior informatif est construit à partir de données provenant de divers sites hydrogéologiques. [START_REF] Egidi | Developments in Bayesian Hierarchical Models and Prior Specification with Application to Analysis of Soccer Data[END_REF] propose une autre application qui concerne l'industrie du sport et plus particulièrement le football et utilise un prior informatif afin de prédire les performances individuelles des joueurs. Ce dernier prior est basé sur les données historiques de chacun des membres de l'équipe sportive.

y i = 1 |T d (i)| k∈T d (i) y i-k 1{DT (i) = HOL}+ 1 |T HOL (i)| k∈T HOL (i) y i-k 1{DT (i) = HOL}+ε i
π(θ|x) = f (x|θ)π(θ) f (x|θ)π(θ)
Dans ce chapitre nous nous appuyons sur des approches similaires pour proposer à notre tour un prior informatif dans le domaine du covoiturage. Plus précisément, nous développons un prior informatif afin d'améliorer les performances d'un modèle bayésien lors du lancement d'une nouvelle ligne de covoiturage. Disposer d'un modèle prédictif performant dès l'ouverture d'une nouvelle ligne de covoiturage est un enjeu stratégique pour assurer son bon fonctionnement dès son lancement. Dans ce contexte, la difficulté principale rencontrée est évidemment le manque de données historiques. L'objet de ce chapitre est de développer des solutions par transfert bayésien afin de pallier à ce problème.

Méthodologie: construction du prior informatif

La construction du prior informatif dans le cas du covoiturage Ecov est motivé d'une part par la ressemblance des données d'un territoire à un autre, et d'autre part par le manque de données historiques à l'ouverture. La Figure 5.1 illustre la ressemblance entre les flux des conducteurs sur trois lignes distinctes. Nous remarquons que les flux quotidiens des trois lignes suivent des évolutions similaires au cours de la même période. Plus précisément, nous pouvons (5.3)

Les lois a priori des hyperparamètres r, m et l sont définies par

r ∼ N (a r , b r ), m ∼ N (1, σ 2 m ), l ∼ N (1, σ 2 l ),
avec a r et b r des réels positifs choisis afin que la loi a priori de r soit noninformative. Au final, le prior informatif est défini par π(θ, k, r, m, l) ∼ π(θ|k, r)π(k|m, l)π(m)π(r)π(l)π(σ 2 m )π(σ 2 l ), (5.4)

avec π(σ 2 l ) ∝ σ -2 l , π(σ m ) 2 ∝ σ -2 m , π(l) ∝ |σ -2 l | 1 2 exp - 1 2 σ -2 l (l -1) 2 , π(m) ∝ |σ -2 m | 1 2 exp - 1 2 σ -2 m (m -1) 2 , π(r) ∝ |b r | -1 2 exp - 1 2 b -1 r (r -a r ) 2 , π(k|m, l) ∝ |s| -d 2 exp - 1 2 l -1 d i=1 (k i -m) 2 , π(θ|k, r) ∝ |r| -d 2 exp - 1 2 (θ -Kµ L ) r -1 (Σ L )(θ -Kµ L ) .
La Figure 5.3 représente de manière schématique la méthode proposée pour la construction du prior informatif.

y i = K k=1 η DT(i-k),DT(i) y i-k + ε i (5.
w i,j |(y i , β β β, t i,j ∈ I s ) ∼ Γ(ν, β s y i ) pour i = 1 . . . N et j = 1, . . . , n i .
(5.7)

Le paramètre β β β est un vecteur de taille S. Il répartit les flux quotidiens sur les S différents intervalles de la journée. Le paramètre ν est le premier paramètre de la loi Gamma modélisant les temps d'attente. On note θ θ θ le vecteur de paramètre des deux étages du modèle. Il est défini θ θ θ = (η ORD, ORD , η SCH,SCH , η ORD,SCH , η SCH,ORD , ν, β 1 , . . . , β S , σ 2 ε ).

Dans la suite nous construisons un prior informatif uniquement pour les paramètres du deuxième étage qui concerne les temps d'attente. On note θ * le vecteur de taille S + 1 qui représente les paramètres de la loi Gamma, on obtient θ θ θ * = (ν, β 1 , . . . , β S ).

(5.8)

Résultats

Détails de la procédure d'évaluation

La construction du prior informatif est faite suivant l'équation 5.4. De plus, pour les lois a priori des hyperparamètres σ 2 l et σ 2 m nous avons utilisé des priors non-informatifs de Jeffreys comme énoncé dans l'introduction. Le paramètre θ θ θ * est estimé grâce à une méthode de Monte-Carlo par chaînes de Markov 5.4. Résultats (MCMC). Concernant l'échantillonneur MCMC, nous utilisons l'échantillonneur NUT [START_REF] Hoffman | The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo[END_REF], qui est implémenté dans le paquet Python pyStan (https://pystan.readthedocs.io). Nous avons utilisé des chaînes MCMC longues pour s'assurer que les lois ont pu converger correctement.

Nous avons évalué les résultats du modèle avec prior informatif et sans prior informatif dans un premier temps sur des données simulées puis sur les données réelles d'Ecov. Dans le cadre de l'évaluation de l'apport du prior informatif on discutera deux aspects i) la qualité de l'estimation des paramètres et ii) les performances prédictives. Plus précisément, la qualité de l'estimation des paramètres va être étudiée dans le cadre de données simulées pour lesquelles les vrais paramètres sont connus. En ce qui concerne les performances prédictives, nous allons évaluer l'apport du prior informatif sur les données réelles issues des lignes de covoiturage.

Jeu de données simulées Présentation des jeux de données simulées

Le protocole de simulation est celui présenté dans le Chapitre 4 et décrit en Appendice 4.7.2. Les valeurs de ν et β β β sont choisies de manière à obtenir des valeurs cohérentes avec les flux et les temps d'attente observés dans le covoiturage. Le jeu de données L correspond aux données simulées sur une période de 390 jours avec en moyenne une dizaine d'observations de temps d'attente par intervalles (six intervalles par jour). Le jeu de données court C correspond aux données simulées sur une période de 90 jours avec en moyenne trois observations de temps d'attente par intervalles. On note ν L ,β L les paramètres de simulation pour le jeu de données L et ν C ,β C pour le jeu de données C. On choisit

β C = λ 1 β L et ν C = λ 2 ν L avec λ 1 = 2 et λ 2 = 1.

Résultats : Estimation des paramètres du modèle

La qualité de l'estimation des paramètres du modèle en utilisant le prior informatif, par rapport à l'utilisation du prior non-informatif est visible sur les lois a posteriori des paramètres β C et ν C . La Figure 5.6 représente les lois a posteriori des paramètres du modèle informatif et non-informatif. Nous constatons que les distributions a posteriori dans le cas du modèle informatif sont moins étendues et plus piquées sur les vraies valeurs des paramètres recherchés. La Table 5.2 compare les écarts-types des paramètres du modèle pour les deux cas (informatif et non-informatif). Les distributions a posteriori des paramètres issues du modèle informatif sont plus précises; le prior informatif facilite l'estimation des paramètres du modèle. Nous rappelons que dans les deux cas, informatif et non-informatif, les lois a posteriori sont issues des chaînes MCMC de même taille, à savoir longue, permettant leurs convergences. De ce fait, nous pouvons nous assurer que dans le cas du modèle non-informatif les lois a posteriori ont pu converger vers leurs lois limites. Par la suite nous allons évaluer l'apport du prior informatif sur des données réelles d'Ecov. Nous allons évaluer la qualité de l'estimation des paramètres mais aussi les performances prédictives.

Jeu de données d'Ecov Présentation des jeux de données d'Ecov

Nous nous intéressons maintenant au transfert de prior bayésien pour les données de covoiturage Ecov. Nous allons appliquer la méthode à trois jeux de données qui correspondent à trois scénarios différents. Le premier jeu de 123 Chapter 5 -Construction d'un prior informatif par transfert bayésien : application au covoiturage données est le jeu de données long, il s'étale sur une plus grande période d'usage. Ce jeu de données est appelé A. Le deuxième jeu de données A correspond à un sous-échantillonnage du jeu de données A. Il s'étale sur une plus petite période et contient à peu près la moitié des observations du jeu de données A. Le troisième et dernier jeu de données B correspond lui aussi à un jeu de données court mais cette fois les observations proviennent d'une toute nouvelle ligne de covoiturage avec encore peu d'observations. Il s'agit de la situation d'historique court sur un nouveau territoire. Les détails des jeux de données sont résumés dans la Table 5.

Transfert bayésien

La méthodologie appliquée est la suivante :

1. On commence par ajuster le modèle non-informatif sur le jeu de données A. Chapter 5 -Construction d'un prior informatif par transfert bayésien : application au covoiturage

Résultats : Estimation des paramètres du modèle

Nous remarquons tout d'abord que les lois a posteriori du modèle non-informatif des jeux de données A et A se ressemblent fortement. Néanmoins, celles du jeu de données A sont plus étendues que celles du jeu de données A ce qui est naturel puisque A est une sous population de A. Les lois a posteriori des paramètres du modèle non-informatif pour le jeu de données B sont aussi moins précises et par ailleurs centrées sur des valeurs différentes. Un des objectifs de la construction du prior informatif est de diminuer l'incertitude autour des distributions a posteriori des paramètres, ce qui va entraîner une amélioration de la précision du modèle. La Figure 5.9 compare les distributions a posteriori des paramètres issues du modèle informatif et non-informatif pour le jeu de données A . De même, la Figure 5.10 représente la même comparaison mais pour le jeu de données B. Ces résultats confirment que le prior informatif est efficace. En effet, la totalité des lois a posteriori des paramètres du modèle informatif dans les deux cas (A et B) sont plus concentrées que celles du modèle non-informatif. On rappelle que, pour les jeux de données A et B, nous n'avons pas accès aux vraies valeurs a posteriori des paramètres. Sachant que les résultats sur les données simulées indiquent que le modèle informatif est plus efficace, on peut espérer que les estimations des répartitions depuis le prior informatif ont un biais comparable à celles provenant du modèle non-informatif. Nous allons voir dans la section suivante comment cette amélioration se traduit sur la performance des lois prédictives a posteriori.

Résultats : Performances prédictives

Le deuxième objectif de la construction du prior informatif est d'améliorer les performances prédictives du modèle. Rappelons qu'en statistique bayésienne nous pouvons estimer la loi prédictive a posteriori. Ceci nous donne la possibilité d'évaluer le gain du modèle informatif par rapport au non-informatif pour différents niveaux de quantiles de la loi prédictive. Il est en effet intéressant ici d'évaluer cet apport en terme de quantiles pour deux raisons. La première raison est que dans le domaine du covoiturage nous nous intéressons à prédire les quantiles de loi des temps d'attente (cf. Chapitre 2). La deuxième raison est que l'amélioration des performances prédictives est largement plus visible au niveau des valeurs des quantiles élevés. Les moyennes des lois a posteriori prédictives du modèle informatif et du modèle non-informatif sont similaires. En revanche, les quantiles des deux lois se distinguent car l'information apportée par le prior informatif rend la première plus précise. Afin de quantifier le gain du modèle informatif nous avons recours au score PE. Ce score, introduit dans cette thèse, est adapté pour le covoiturage. Il mesure le pourcentage de bonnes information supplémentaire. Effectivement, dans le cas du jeu de données B, nous observons que la totalité des prédictions du quantile posterior issues du modèle informatif, au niveau 95%, ne dépassent pas une erreur de plus de 9 minutes. Par contre, pour le même seuil de 9 minutes, le quantile posterior du modèle non-informatif atteint seulement 40% de ses prédictions. Cette nette amélioration illustre l'intérêt et l'efficacité du prior informatif sur un jeu de données court.

(δ) = 1 Ñ Ñ i=1 1{|q α ( wi ) -w i | < δ} (5.9) avec q α (
Une autre façon de résumer la loi prédictive a posteriori des temps d'attente est de calculer son intervalle de plus haute densité, HDI (Highest Density Interval). Un point se trouvant dans cet intervalle a une crédibilité plus élevée que tout point se trouvant à l'extérieur de ce même intervalle. Nous fixons donc un niveau de crédibilité à 94% et comparons les intervalles de plus haute densité des lois a posteriori prédictives du modèle informatif et non-informatif. De cette manière nous étudions les performances prédictives des deux modèles. Cette comparaison doit au préalable s'accompagner d'une étude sur les scores PE S 1 2 3 4 5 6 A 16% 14% 8% 12% 1% 1% B 48% 20% 30% 30% 22% 26% Table 5.4 -Pourcentages de diminution de l'intervalle de plus haute densité des lois prédictives a posteriori issues du modèle informatif sur les six créneaux du premier jour de prédiction.

Conclusions

Dans ce chapitre nous avons développé une méthode pour la construction d'une loi a priori informative par transfert bayésien dans le cas d'un jeu données court.

Chapter 5 -Construction d'un prior informatif par transfert bayésien : application au covoiturage Cette construction est faite à partir d'un jeu de données long L. Puis, le modèle informatif est utilisé sur un jeu de données court C, qui est supposé ressembler au jeu de données long. Nous avons montré, à l'aide de données simulées et réelles, que l'apport de cette loi a priori informative améliore la qualité de l'estimation des paramètres du modèle ainsi que les performances prédictives.

Cette méthode proposée repose sur l'incorporation des informations apprises par le modèle sur le jeu de données long L à la loi a priori du modèle sur le jeu de données court C. Une perspective possible des travaux de ce chapitre est l'utilisation d'autres données disponibles, externes à celles utilisées actuellement par le modèle, afin de mieux construire le prior informatif. En effet, Cucchi et al., 2019 ont proposé de mutualiser des données de plusieurs sites différents pour renforcer la qualité du prior. Nous pensons que cette perspective, en parallèle avec notre méthode, pourra améliorer davantage encore les performances du prior informatif. L'étape de la mise en production du modèle va donner la possibilité à la société Ecov de proposer un service de covoiturage plus complet. En effet, les informations concernant le flux de trafic et les temps d'attente vont alimenter les différents outils d'information et de synthèse disponibles de la société. Ensuite, chaque outil expose ces informations de manière à répondre aux différents besoins que ce soit des passagers ou des différentes équipes au sein d'Ecov. Dans les faits, un passager peut s'informer de son temps d'attente et du potentiel de covoiturage à plusieurs moments. Par exemple, la veille d'un covoiturage le passager peut consulter ces informations à l'aide du site internet et de l'application mobile d'Ecov. Également, au moment même du covoiturage le temps d'attente est affiché sur la borne d'informations des voyageurs (BIV) qui se trouve à l'arrêt de covoiturage. De plus, l'équipe "gestion relation utilisateur" d'Ecov qui accompagne les passagers lors d'un covoiturage a lui aussi besoin de connaître les temps d'attente et les flux de trafic potentiel. Dans ce cas, l'équipe "gestion relation utilisateur" consulte un logiciel de tableau de bord. C'est dans le but de répondre à ces besoins que nous avons développé un paquet en Python appelé bayesian-hierarchical-model. Ce paquet permet la collecte des données, le pré-traitement des données, l'apprentissage du modèle et l'exposition des résultats du modèle vers les différents canaux énumérés précédemment (site Chapter 6 -Valorisation Industrielle : Mise en production du modèle hiérarchique bayésien internet, application mobile, la BIV, logiciel de tableau de bord).

Le paquet bayesian-hierarchical-model a été conçu à partir des outils développés dans ce manuscrit et constitue un aboutissement de ces travaux de thèse. Nous allons décrire les différentes étapes et la structure du paquet ainsi que l'automatisation de son utilisation.

Description des données et du modèle à industrialiser

Rappelons

y i = K k=1 η DT(i-k),DT(i) y i-k + ε i (6.
i,j |(y i , β β β, t i,j ∈ I s ) ∼ Γ(ν, β s y i ) pour i = 1 . . . N et j = 1, . . . , n i . (6.3)
Le paramètre β β β est un vecteur de taille S et il répartit les flux quotidiens sur les S différents intervalles de la journée. Le paramètre ν est le premier paramètre de la loi Gamma modélisant les temps d'attente. des outils de gestion des différents paquets Python. Anaconda (Anaconda Software Distribution, 2020) permet la gestion de différents paquets proposant directement un environnement où la plupart des paquets populaires en science des données sont pré-installés. Un autre outil est la bibliothèque Matplotlib [START_REF] Hunter | Matplotlib: A 2D graphics environment[END_REF] pour visualiser les données sous formes de graphiques. Cette liste n'est pas exhaustive mais donne une vue globale sur les différents outils utiles lors du prototypage d'un modèle d'apprentissage statistique. Concernant l'industrialisation du modèle, des outils de gestion des versions du code tel que le Git [START_REF] Chacon | Pro git[END_REF] et GitLab 1 sont utilisés. Ils donnent la possibilité de garder un historique de l'évolution du code et de mieux partager le code entre collaborateurs. De plus, plusieurs bonnes pratiques existent afin de bien utiliser ce genre d'outils, il s'agit d'un sujet important pour la bonne organisation des équipes de développement. Ensuite, un autre outil incontournable est PyCharm 2 qui est un environnement de développement intégré pour le langage Python. Il permet de faciliter les bonnes pratiques de développement en offrant un environnement hautement configurable et contient une extension pratique de débogage graphique. Au moment du déploiement du modèle nous avons fait le choix d'utiliser une plateforme de planification de flux de travail qui s'appelle Apache Airflow 3 . La dernière étape à ne pas négliger est la documentation qui doit être claire et à jour, pour cela nous avons fait le choix d'utiliser Sphinx 4 qui est un générateur de documentation.

Il est important de noter que le prototypage du modèle est un travail de recherche où les aspects statistiques du modèle sont privilégiés ainsi que les performances prédictives. Tandis que pour l'industrialisation du modèle, on s'intéresse à la qualité du code, la rapidité d'exécution du modèle et l'automatisation de son entraînement et de son déploiement.

La structure du paquet

Le paquet développé s'appelle bayesian-hierarchical-model et est entièrement écrit en Python suivant les standards des paquets utilisés en apprentissage automatique. Nous avons emprunté la structure des modèles implémentés dans la bibliothèque Scikit-Learn [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF]. Cela signifie que nous avons développé une méthode nommée ".fit()" afin d'entraîner le modèle sur les données d'entraînement et une méthode ".predict()" permettant la prédiction du modèle sur les nouvelles données. Nous avons aussi implémenté plusieurs autres méthodes et fonctions concernant la récupération et le pré-traitement des données et la simulation des données jouets. La structure du paquet, illustrée Chapter 6 -Valorisation Industrielle : Mise en production du modèle hiérarchique bayésien

Après avoir détaillé la structure du paquet nous illustrons son utilisation par des exemples concrets.

Utilisation du paquet bayesian-hierarchical-model

Simulation de flux de trafic. L'utilisateur du paquet a la possibilité de générer des données simulées selon le protocole de simulation détaillé dans le Chapitre 4. La fonction développée pour cela s'appelle TrafficFlowGenerator. Les arguments nécessaires sont le nom de la colonne contenant les dates ("date_col"), la période concernée pour la simulation en donnant la date de début et de fin ("start" et "end") au format "aaaa-mm-jj", le pays concerné ("country") et la zone souhaitée ("zone") quand cela s'applique pour le pays choisi. Nous devons préciser le pays ainsi que la zone voulue pour permettre au simulateur de prendre en compte l'impact des types de jour sur le trafic. L'exemple suivant illustre l'utilisation de cette fonction pour la simulation de trafic de flux pour la période "2019-01-01" à "2020-01-01" en France dans la zone "A" (Poitiers, Bordeaux, Limoges, Clermont-Ferrand, Lyon, Grenoble, Dijon, Besançon). Les données obtenues sont illustrées par la • Le DataFrame 5 fitting_df contenant les informations de la période d'entraînement concernant les flux de trafic et les types de jour.

• Le DataFrame prediction_df contenant les informations de la période de prédiction concernant les flux de trafic et les types de jour.

• La matrice W contenant les temps d'attente observés pour la période d'entraînement.

• La matrice W_test contenant les temps d'attente observés pour la période de prédiction.

Ensuite, la phase de prédiction est lancée avec la méthode ".build_model()" qui déclenche l'entraînement et l'inférence du modèle. Nous obtenons par conséquent les distributions des lois a posteriori des paramètres du modèle et des lois prédictives a posteriori des temps d'attente. Dans l'exemple ci-dessous toutes ces lois sont enregistrées dans la variable "samples".

from b a y e s h i e r a r c h i c a l . python . b a y e s i a n _ h i e r a r c h i c a l _ m u l t i l e v e l import B a y e s i a n H i e r a r c h i c a l M u l t i l e v e l i t e r a t i o n s = 350 model = B a y e s i a n H i e r a r c h i c a l M u l t i l e v e l ( s t a r t = " 2018 -01 -01 " , w e e k _ t o _ p r e d i c t= " 2018 -11 -08 " , num_train_wt =20 , i t e r a t i o n s = i t e r a t i o n s , c h a i n s =1) Construction du modèle sur les données réelles Dans le cas de l'utilisation du modèle sur des données réelles, voir des données externes pour d'autres applications, il suffit de créer l'objet BayesianHierarchicalMultilevel() avec des arguments spécifiques. L'exemple suivant illustre le cas d'utilisation du modèle sur des données réelles d'Ecov. Les arguments concernant la simulation des données ("start","week_to_predict", "wt_num) sont remplacés par les DataFrame "fitting_df" et "prediction_df" et les matrices "W" et "W_test". L'utilisateur doit préparer ces données en respectant les formats qui sont décrits dans la documentation. Évidement, les arguments "prediction_df" et "W_test" sont facultatifs car en production nous n'avons pas accès aux données futures. L'exemple suivant illustre l'utilisation du modèle hiérarchique bayésien sur des données réelles d'Ecov qu'on a nommées "fitting_df_real", "prediction_df_real", "W_real" et "W_test_real".

Conclusions

Dans ce chapitre nous avons présenté la structure du paquet Python bayesianhierarchical-model ainsi que ses fonctions. Puis, nous avons exposé comment son automatisation est faite dans le but d'avoir des prédictions en temps réel. Nous rappelons que le paquet est développé en Python et en ce qui concerne l'inférence bayésienne nous avons utilisé PyStan. Une perspective possible est le changement de cette dernière technologie en la remplaçant par PyMC3 7 qui est un paquet pour la modélisation statistique bayésienne et l'apprentissage automatique probabiliste en python. Une deuxième perspective qui est apparue en fin de thèse est la définition d'un nouveau type de jour qui correspond aux périodes de confinement liée à la crise sanitaire. Les études préliminaires ont Le Chapitre 2 présente le principe de régression quantile ainsi que quatre approches possibles de celle-ci. Nous rappelons que notre problématique est formulée de la manière suivante "Quelle est la probabilité qu'un passager attende moins de y minutes avec α% de chance?". Les données utilisées dans ce chapitre sont celles provenant de la structure initiale du service de covoiturage, à savoir le réseau de covoiturage. Nous proposons une comparaison des performances des quatre algorithmes de régression quantile retenus.

Le Chapitre 3 détaille, quant à lui, un processus de travail combinant les outils des Systèmes d'Information Géographiques (SIG) et de la science des données afin d'exploiter les données issues du covoiturage. Nous rappelons qu'à partir de ce chapitre les données exploitées sont celles provenant de la nouvelle structure du service, à savoir les lignes de covoiturage. Les nouveautés par rapport au chapitre précédent sont la collecte des traces GPS des conducteurs et la structure améliorée du service. Suite à l'analyse de ces données nous obtenons deux résultats : i) que la correspondance porte-à-porte de trajectoires complètes de l'origine à la destination est un obstacle structurel pour la transformation du covoiturage en un service de transport en commun et ii) que les arrêts de covoiturage fixes favorisent la rencontre entre les conducteurs et les passagers et diminuent les temps d'attente.

Le Chapitre 4 est consacré à la construction d'un modèle hiérarchique bayésien permettant la prévision des flux de trafic et des temps d'attente. Nous disposons des données industrielles provenant de deux sources distinctes : la première source concerne les traces GPS des conducteurs et la seconde les temps d'attente observés. Le modèle est constitué de deux étages imbriqués : le premier étage s'intéresse à la modélisation des flux de trafic quotidien tandis que le second étage modélise les temps d'attente pour différents créneaux de la journée.

Le Chapitre 5 expose une méthode de construction de loi a priori informative pour le modèle hiérarchique bayésien en situation d'historique court. Nous construisons cette loi informative à partir d'un jeu de données long L et d'une loi a priori non-informative. Par la suite, le modèle avec prior informatif est appliqué sur un jeu de données court C supposé semblable au jeu de données L. Nous montrons à l'aide d'applications, sur des données simulées et réelles, que l'apport de la loi a priori construite améliore la qualité d'estimation des paramètres et les performances prédictives du modèle.

Enfin, le Chapitre 6 décrit les étapes nécessaires pour la mise en production du modèle hiérarchique bayésien. Le résultat concret du chapitre est la production d'un paquet Python bayesian-hierarchical-model pour la société Ecov.

Perspectives

Dans le Chapitre 3 le processus de travail constitue un prototype solide pour des flux de travail futurs de la société. Une perspective possible est l'industrialisation et la mise en production de ce processus de travail pour l'intégrer directement dans les outils de la société.

Comme nous l'avons expliqué dans le Chapitre 4, les données utilisées pour la construction du modèle hiérarchique sont issues des conducteurs et des passagers du service. Une première piste possible est l'inclusion de données externes. Des études ont montré que les données météorologiques peuvent améliorer la compréhension des flux de trafic. De plus, le modèle hiérarchique bayésien a la possibilité (avec des modifications mineures) de prendre en compte des données de comptage de trafic provenant de sources diverses. Ces comptages permettront d'améliorer à leur tour la connaissance des flux de trafic. Une dernière source de données à envisager est l'utilisation des informations des services "concurrents" pour un même territoire (la fréquence des transports en commun, les Véhicules de Tourisme avec Chauffeur disponibles, les autres applications de covoiturage). Cette dernière source de données a pour but d'améliorer la connaissance des temps d'attente. La seconde piste d'amélioration concerne directement la construction du modèle. Nous rappelons que le coefficient η DT(i-k) permet de modéliser l'impact des 1 ≤ k ≤ K jours qui précèdent le i ième jour. Une perspective envisageable est de modéliser de la même manière l'impact entre les différents créneaux de la journée. Le coefficient η DT(i-k) deviendrait une matrice représentant les interactions entre les différents créneaux de la journée.

La méthode proposée dans le Chapitre 5 repose sur l'incorporation des informations apprises par le modèle sur le jeu de données long L à la loi a priori du modèle sur le jeu de données court C. Une perspective possible des travaux de ce chapitre est l'utilisation d'autres données disponibles, externes à celles utilisées actuellement par le modèle, afin de mieux construire le prior informatif.

Un prolongement possible des travaux présentés dans le Chapitre 6 consiste à ajouter au paquet bayesian-hierarchical-model la méthode de construction du prior informatif par transfert bayésien introduit dans le Chapitre 5. Un autre prolongement envisageable est l'amélioration de l'exposition du modèle pour les clients interne et externe. Cela consiste à implémenter des fonctions permettant une meilleure visualisation des résultats obtenus par le modèle. Un dernier prolongement éventuel est de remplacer la librairie Python PyStan par une librairie concurrente nommée PyMC3.

Title: Potential and prediction of waiting times for carpooling in a territory

Keywords : carpooling, waiting times, quantile regression, bayesian hierarchical model, informative prior Abstract : This thesis focuses on the potential and prediction of carpooling waiting times in a territory using statistical learning methods. Five main themes are covered in this manuscript. The first presents quantile regression techniques to predict waiting times. The second details the construction of a workflow based on Geographic Information Systems (GIS) tools in order to fully leverage the carpooling data.

In a third part we develop a hierarchical bayesian model in order to predict traffic flows and waiting times. In the fourth part, we propose a methodology for constructing an informative prior by bayesian transfer to improve the prediction of waiting times for a short dataset situation. Lastly, the final theme focuses on the production and industrial exploitation of the bayesian hierarchical model.

Figure 1 . 1 -

 11 Figure 1.1 -Les acteurs industriels du covoiturage du quotidien, image provenant du site de la Fabrique des Mobilités https://wiki. lafabriquedesmobilites.fr/wiki/Le_covoiturage

Figure 1 . 3 -

 13 Figure 1.3 -Configuration d'un arrêt de covoiturage connecté. La structure orange ressemble à un abribus. Le passager informe les conducteurs potentiels de sa demande de covoiturage à l'aide de la console, demande qui est ensuite affichée sur le panneau électronique lumineux au bord de la voirie. Un conducteur peut récupérer le passager en toute sécurité sur la place de stationnement réservée à cette occasion. Photographie reproduite avec l'autorisation d'Ecov.

Figure 1 . 6 -

 16 Figure 1.6 -Histogramme des temps d'attente observés sur le réseau de covoiturage du "Vexin" pendant la période du 2016-01-18 au 2018-10-02.

Figure 1

 1 Figure 1.9 -Correspondance spatio-temporelle porte-à-porte, à gauche les conditions d'appariement sont moins strictes et à droite les conditions d'appariement sont restreintes. Le sous-cube bleu représente l'origine (quartier résidentiel), le vert la destination (lieu de travail), et les trajectoires qui partagent les mêmes sous-cubes d'origine et de destination sont considérées comme des correspondances porte-à-porte.

Figure 1 .

 1 Figure 1.10 -Traces GPS de plusieurs conducteurs provenant du service de covoiturage stochastique opérationnel "Lane" à Lyon. Les différents clusters de traces GPS représentent les mises en relation de porte-à-porte. Les losanges pleins représentent les origines et destinations. Les points de rencontre sont les cercles noirs pleins, notés B = Bourgoin, S = St-Priest.

Figure 1 .

 1 Figure 1.11 -Processus de travail Data Science-SIG pour l'analyse des traces GPS des conducteurs dans un service de covoiturage stochastique. (Gauche) Sources de données d'entrée spatio-temporelles. (Centre) Tâches de traitement des données et de géotraitement. (Droite) Résultats générés.

Figure 1 .

 1 Figure 1.13 -Flux quotidiens de conducteurs au sein du service de covoiturage Lane, du 2018-05-15 au 2019-05-31. Les jours de semaine ordinaires (ORD) sont en orange, les vacances scolaires (SCH) en vert et les jours fériés/week-end (PWE) en bleu.

Figure 1 .

 1 Figure 1.14 -Temps d'attente des passagers (en minutes) au sein du service de covoiturage Lane du 2019-10-22 au 2020-01-15. Les jours de semaines ordinaires (ORD) sont en orange, et les vacances scolaires (SCH) en vert.

Figure 1 .

 1 Figure 1.15 -Schéma du modèle hiérarchique bayésien pour la prédiction des flux de conducteurs et des temps d'attente des passagers. Les données d'entrée (traces GPS des conducteurs) sont en gris, les modèles hiérarchiques en vert, les paramètres du modèle en orange et les résultats du modèle en violet.

Figure 1 .

 1 Figure 1.16 -Les box-plot gris sont les pseudos temps d'attente simulés pour des intervalles de 3 heures, pour les Ñ jours de prédiction. Les prédictions des quantiles des pseudos temps d'attente niveau 50%, 75% et 95% sont respectivement les cercles violet clair, violet moyen clair et violet foncé.

Figure 1 .

 1 Figure 1.17 -Évolution de l'erreur MSE pour la prédiction du flux quotidien de conducteurs pour différentes périodes. Le modèle hiérarchique bayésien BHML est en violet, le modèle fréquentiste de base BASE en noir et le modèle bayésien Prophet PROP en vert.

Figure 1 .

 1 Figure 1.18 -Prédictions du flux quotidien de conducteurs pour six différentes périodes de test. Les flux journaliers de conducteurs observés sont en bleu. Le modèle hiérarchique bayésien BHML est en violet, le modèle fréquentiste de base BASE en noir et le modèle bayésien Prophet PROP en vert.

Figure 1 .

 1 Figure 1.19 -Les box-plot gris sont les pseudos temps d'attente observés pour des intervalles d'une heure. Les prédictions des quantiles des pseudos temps d'attente niveau 50%, 75% et 95% sont respectivement les cercles violet clair, violet moyen clair et violet foncé.

Introduction

  avec a r et b r des réels positifs fixés afin que la loi a priori de r soit noninformative. À terme, le prior informatif est défini parπ(θ, k, r, m, l) ∼ π(θ|k, r)π(k|m, l)π(m)π(r)π(l)π(σ 2 m )π(σ 2 l ) (1.5)avec les lois a priori des hyperparamètres σ 2 l et σ 2 m des priors non-informatifs de Jeffreys. Le schéma de la méthodologie du transfert bayésien par construction du prior est illustré à la Figure1.20.

Figure 1 .

 1 Figure 1.20 -Schéma de la méthode proposée pour la construction du prior informatif en situation d'historique court.

Figure 1 .

 1 Figure 1.21 -Histogrammes des lois a posteriori des paramètres du modèle informatif en orange et du modèle non-informatif en vert.

Figure 1 .

 1 Figure 1.22 -Histogrammes des lois a posteriori des paramètres du modèle de temps d'attente pour le jeu de données B. En bleu, les lois a posteriori avec un prior non-informatif et en rouge, les lois a posteriori avec un prior informatif.

Figure 1 .

 1 Figure 1.23 -Comparaison des lois a posteriori prédictives des temps d'attente pour le premier jour de prédiction pour les six créneaux du jeu de données B. La courbe en rouge est la loi prédictive lorsque le prior est informatif et en noir sa région HDI (Highest Density Interval). En violet, la région HDI de la loi prédictive des temps d'attente en non-informatif.

Figure 1 .

 1 Figure 1.24 -Comparaison de le score PE du quantile posterior 95% pour le jeu de données B. En bleu le score PE pour le modèle avec prior non-informatif et en rouge en avec prior informatif.

Introduction

  modèle est Airflow 2 . Il s'agit d'une plateforme de planification de flux de travail. La Figure1.26 illustre le schéma suivi pour l'intégration des données, leurs traitements, l'entraînement du modèle et l'exposition des résultats. Ce schéma fonctionne selon le concept d'ETL (Extract Load Transform). En particulier, nous avons segmenté les différentes tâches nécessaires en trois ETL dépendants les uns des autres.

Figure 1 .

 1 Figure 1.26 -Schéma de la structure des ETL mis en place afin d'exposer en production les prédictions des flux de trafic et des temps d'attente aux différents services d'Ecov.

Figure 2

 2 Figure 2.1 -Schematic of an Ecov carpooling network. The blue circles are the physical meeting points where the passenger can request a carpooling, and links between the circles indicate that a carpooling service is assured for this origin-destination pair.

Figure 2

 2 Figure 2.2 -Map of the "Vexin" carpooling network. The physical meeting points where passengers are able to make a carpooling request are indicated by the pink circles with a thumb logo.

  Figure 2.3 are illustrated the curves for the Pinball loss for α = 0.25, 0.5, 0.75, 0.95. Apart from the green curve (α = 0.5) which is symmetric, all the other Pinball losses are asymmetric about the origin (blue α = 0.95, orange α = 0.75, red α = 0.25). The induced asymmetry allows for y and -y, whilst they have the same absolute loss |y| and squared loss y 2 , to have different Pinball losses. Since the curves of the Pinball loss consist of straight lines which change direction at 0, they resemble the trajectories of a ball in a Pinball arcade machine, hence their name.From a mathematical point of view, the sign of y in the Pinball loss plays an important role, unlike the more commonly considered squared and absolute losses which are non-negative loss functions.

Figure 2 . 3 -

 23 Figure 2.3 -Pinball loss functions ρ α for different values of α. Blue: α = 0.95, orange α = 0.75, green α = 0.5, red α = 0.25.

Figure 2

 2 Figure 2.4 -Classical Random Forest with two decision trees trained on the same 2-dimensional re-sampled subset D* 19 = {X X X 1 , • • • , X X X 19 } (black circles).The remaining observed data {X X X 20 , . . . , X X X 23 } (gray circles) are excluded from their construction. The resulting partition classes areC 1 1 ∪ • • • ∪ C 1 8 and C 2 1 ∪ • • • ∪ C 2 8 .A new data point is x x x (black cross).

Figure 2 .

 2 Figure 2.5 illustrates the curves for the scoring function ψ 0 (y) for α = 0.25, 0.5, 0.75, 0.95.

Figure 2 . 5 -

 25 Figure 2.5 -Scoring function ψ 0 (y) for different values of α. Blue: α = 0.95, green α = 0.75, orange α = 0.5, purple α = 0.25.

Figure 2

 2 Figure 2.7 -Illustrative examples of loss metrics of Quantile Regressions on a test sample of size 8. The observed values are the black dots and the predicted 75% quantiles are the red crosses.

Figure 2

 2 Figure 2.9 -Histogram of the waiting times in the Vexin carpooling network during the period 2016-01-18 to 2018-10-02.

Figure 2 .

 2 Figure 2.10 -Box plot of the waiting times per day of the week in the Vexin carpooling network during the period 2016-01-18 to 2018-10-02.

Figure 2 .

 2 Figure 2.11 -Clustered destination neighbourhoods in the Vexin carpooling network.

Figure 2 .

 2 Figure 2.12 -Number of carpooling trips for the five most frequented origindestination pairs in the Vexin carpooling network. (Left) Non-clustered origindestinations. (Right) Clustered origin-destinations.

Figure 2 .

 2 Figure 2.13 -Aggregated counts of the number of carpooling requests in the Vexin carpooling network for 2016-01-18 to 2018-10-02. (Right) Per day of week. (Left) Per hour.

Figure 2 .

 2 Figure 2.14 -Predicted quantiles from the Regression models for the waiting times in the Vexin carpooling network. The predictions for α = 0.25 are in green, and for α = 0.75 in brown, for a random test fold from the cross validation for the models GBM (top left), QRF (top right), GRF (bottom left) and QuantReg (bottom right). The black circles are the observed waiting times.

Figure 3

 3 Figure 3.1 -Spatio-temporal door-to-door matching fragments the population of mutualisable trajectories. (Left) Relaxed matching conditions. (Right) Restricted matching conditions. Blue sub-cube represents the origin (residential neighbourhood), green the destination (workplace), and trajectories which share the same origin and destination sub-cubes are considered to be door-to-door matches.

Figure 3

 3 Figure 3.2 -Probability of door-to-door matches for uniformly distributed drivers and passengers, as a function of the number of sub-cube partition classes. Higher number of sub-cubes represent more stringent spatio-temporal matching conditions.

  3, we analyse 121 GPS traces of drivers who travelled from the Bourgoin meeting point (solid black circle labelled B) to the St-Priest meeting point (solid black circle labelled S) in the Lane carpooling service during the morning operating hours (06:30 to 09:00) for the work week 2019-11-25 to 2019-12-01.

Figure 3

 3 Figure 3.3 -Spatio-temporal door-to-door matching fragments the number of mutualisable trajectories in an operational stochastic carpooling service. The clusters of GPS traces of door-to-door matches are colour coded, with the GPS points as the solid circles, and the origins/destinations as the solid diamonds. The meeting points are the solid black circles, denoted B = Bourgoin, S = St-Priest.

Figure 3 . 4 -

 34 Figure 3.4 -Data science-GIS workflow for the analysis of driver GPS traces in stochastic carpooling service. (Left) Spatio-temporal input data sources. (Centre) Data wrangling and geoprocessing tasks. (Right) Generated outputs.

Figure 3

 3 Figure 3.6 -(Left) Schematic diagram of the Lane carpooling network, which resembles the geographically restrained trajectories of a public transport service. Reproduced with permission from Ecov. (Right) Carpooling network represented as a directed graph. Nodes are the meeting points, edges connect meeting points whenever a carpooling service between them is assured.

Figure 3

 3 Figure 3.7 -Topological simplification of a GPS trace in the Bourgoin > St-Priest carpooling line, during its morning operating hours 06:30-09:00, on 2019-11-28. The complete GPS trace are the 530 blue circles; the sequence of five nodes, as its simplified topology, are the orange arrows, and the orange diamonds are the origin, carpooling meeting points, destination nodes. The meeting points are S = St-Priest, V = Villefontaine, and B = Bourgoin.

Figure 3 . 8 -

 38 Figure 3.8 -Matching on meeting points increases the number of driver spatiotemporal matches in comparison to door-to-door matching, during a single 30 minute period (08:00-08:30), on 2019-11-28. The orange arrows are the 14 GPS traces which are both meeting point and door-to-door, and the blue arrows are the 17 GPS traces which are meeting point matches but not door-to-door match matches. The diamonds are the origin/destination points. The solid black circles are the meeting points: S = St-Priest, V = Villefontaine, and B = Bourgoin.

  5 -Waiting time predictions for a carpool request on the Bourgoin > St-Priest carpooling line, per 15 minute intervals, during the morning operating hours 06:30-09:00 for a typical day from 2019-09-02 to 2020-02-17.

Figure 3 . 9 -

 39 Figure 3.9 -Box plots of observed waiting times on the Bourgoin > St-Priest carpooling line, per 15 minute intervals during the morning operating hours 06:30-09:00, from 2019-09-02 to 2020-02-17.

Figure 3 .

 3 Figure 3.10 -RMSE between the observed and predicted waiting times on the Bourgoin > St-Priest carpooling line, per 15 minute intervals during the morning operating hours 06:30-09:00, from 2019-09-02 to 2020-02-17.

3. 5 .

 5 Figure 3.11 -Likely driver itineraries from the TomTom route finder API in the same geographical region as the Bourgoin > St-Priest carpooling line. The origins and destinations (county centroids) are the orange diamonds. The solid black circles are the meeting points: S = St-Priest, B = Bourgoin.

Figure 3 .

 3 Figure 3.12 -Evolution of the predicted passenger waiting time as a function of the driver participation rate in the Bourgoin > St-Priest carpooling line, during the morning operating hours 06:30-09:00, from 2019-07-25 to 2020-02-17.

  summer holiday season in France). The daily driver flows in the Lane network are presented in Figure4.1, where we enumerate each driver GPS trace, rather than each unique driver. The ordinary weekdays (ORD) are in orange, the school holidays (SCH) in green, and the public holidays/weekends (PWE) in blue. The classic temporal cycles of driver flow data indicate that a moving average is a relevant approach for prediction.

Figure 4

 4 Figure 4.1 -Daily driver flow in the Lane carpooling service, from 2018-05-15 to 2019-05-31. The ordinary weekdays (ORD) are in orange, the school holidays (SCH) in green and the public holidays/weekend days (PWE) in blue.

Figure 4

 4 Figure 4.2 -Passenger waiting times (in minutes) in the Lane carpooling service from 2019-10-22 to 2020-01-15. The ordinary work days (ORD) are in orange, and the school holidays (SCH) in green.

Figure 4 . 3 -

 43 Figure 4.3 -Flowchart of Bayesian hierarchical model for driver flow and passenger waiting time prediction. The input data (driver GPS traces) are in grey, the hierarchical models in green, the model parameters in orange, and the model outputs in purple.

Figure 4

 4 Figure4.4 -Perceived and pseudo waiting times for the case of two passengers at a carpooling meeting point. Passenger A is at the head of the queue so their perceived waiting time (blue brace) coincides with their pseudo waiting time. For passenger B, their pseudo waiting time (grey brace) is the difference between their departure time and the departure time of passenger A, which is shorter than their perceived waiting time (green brace).

  3), and the corresponding N × S posterior prediction matrices W = { W W W (1) , . . . , W W W (J) }. The data from these N = 365 days from 2018-01-01 to 2018-12-31 form the reference training data set. With the same parameters, we simulate a further Ñ = 5 days (2019-01-01 to 2019-01-05) as the oracle test data set of Ñ × S matrices W test = Furthermore, from the training data only, for these same extra Ñ days, we generate the corresponding Ñ × S posterior prediction matrices

Figure 4 . 7 -

 47 Figure 4.7 -Aggregate driver flow by day type, for the Lane carpooling service from 2018-05-15 to 2019-05-31. The ordinary weekdays (ORD) are in orange, the school holidays (SCH) in green and the public holidays/weekend days (PWE) in blue.

Figure 4

 4 Figure 4.8 -Evolution of the goodness-of-fit of the daily driver flow estimations over the training period (2018-05-15 to 2019-05-19, test scenario # 6). Goodnessof-fit is measured by the weekly aggregated estimation MSE. Bayesian hierarchical multi-level BHML is in purple, frequentist baseline BASE in black, and Bayesian Prophet PROP in green.

Figure 4 .

 4 Figure 4.11 -Mean driver flows for 15 minute interval for each weekday for the Lane carpooling service, from 2018-05-15 to 2019-05-31. Monday is in blue, Tuesday in orange, Wednesday in green, Thursday in pink, Friday in violet.

4. 4 .

 4 Model validation with the Lane carpooling service from 2019-07-25 to 2020-02-17. Although there are S = 24 hourly intervals, only those 6 which correspond to the service operating hours (06:00-09:00 and 16:00-19:00) contain any observed passenger waiting times.

Figure 4 .

 4 Figure 4.14 -Box plots of the weekly number of observed pseudo waiting times for each hourly interval for weekdays from 2019-07-25 to 2020-02-17.

Figure 4 .

 4 Figure 4.17 -Evolution of passenger waiting times as a function of increased passenger arrival rates. (Left) Pseudo waiting times. (Right) Perceived waiting times. λ 1 is the base passenger arrival rate, λ 2 is an increase by 50% and λ 3 is an increase of 100%.

  day i is an ordinary workday or a weekend day HOL if day i is a school or a public holiday; Dans de nombreuses situations en apprentissage statistique on dispose d'un modèle qui a été entraîné sur un premier jeu de données et on souhaite utiliser cet apprentissage pour un second jeu de données. Une situation typique est lorsque le second jeu de données est de petite taille et suit une distribution proche de celle du premier jeu de données. L'approche retenue s'appuie sur la construction d'un prior informatif pour l'ajustement d'un prédicteur provenant d'un modèle hiérarchique bayésien. Et plus spécifiquement celui du Chapitre 4. L'application visée ici, comme pour le reste de ce manuscrit, concerne le covoiturage. Nous rappelons tout d'abord quelques notions de base du paradigme bayésien dans le cas des modèles paramétriques.

  5.2. Méthodologie: construction du prior informatifaussi observer la manière dont les flux sont distribués dans une même journée. La Figure5.2 illustre le comportement des flux moyens par intervalles de 15 minutes pour les différents jours de la semaine.

Figure 5

 5 Figure 5.1 -Exemple des flux quotidiens sur trois lignes de covoiturage différentes.

Figure 5

 5 Figure 5.2 -Flux moyens des conducteurs par intervalles de 15 minutes pour chaque jour de semaine sur la période du 2018-05-15 au 2019-05-31. Le lundi est en bleu, le mardi en orange, le mercredi en vert, le jeudi en rose, le vendredi en violet.

5. 3 .Figure 5

 35 Figure 5.3 -Schéma de la méthode proposée pour la construction du prior informatif en situation d'historique court.

  Figure 5.5 -À gauche sont représentés les temps d'attente simulés du jeu de données long L et à droite les temps d'attente simulés du jeu de données court C.

Figure 5

 5 Figure 5.6 -Distributions des lois a posteriori des paramètres du modèle informatif en orange et du modèle non-informatif en vert.

  Récapitulatif des jeux de données pour les trois différents scénarios.

Figure 5

 5 Figure 5.7 -Distributions des temps d'attente observés pour les trois différents jeux de données. En bleu, le jeu de données A, en vert, le jeu de données A et en jaune le jeu de données B.

  , on construit le prior informatif comme détaillé en Section 5.2 pour le jeu de données B.L'étape 2 a pour but de vérifier si l'apport du prior informatif est correct car le jeu de données A est un sous-échantillon du jeu de données A tandis que la dernière étape correspond à l'enjeu essentiel de la construction du prior informatif. Effectivement, l'amélioration des performances du modèle sur un nouveau jeu de données, en l'occurrence le jeu de données B à partir d'un autre jeu de données long, est stratégique d'un point de vue industriel. Nous rappelons que les paramètres du modèle concernés par la construction du prior informatif sont θ θ θ * = (ν, β 1 , . . . , β 6 ). Ce sont les paramètres du deuxième étage du modèle hiérarchique bayésien détaillé en Section 5.3 modélisant les temps d'attente. Nous avons aussi ajusté le modèle non-informatif sur les deux jeux de données A et B de façon à pouvoir comparer l'intérêt du prior informatif. Les lois a posteriori obtenues avec prior non-informatif sont illustrées par la Figure5.8.

Figure 5 . 8 -

 58 Figure 5.8 -Distributions a posteriori des paramètres du modèle non-informatif concernant les temps d'attente pour chaque jeu de données. En bleu, le jeu de données A, en jaune, le jeu de données B et en vert, le jeu de données A .

Figure 5

 5 Figure 5.9 -Histogrammes des lois a posteriori des paramètres du modèle de temps d'attente pour le jeu de données A . En bleu, les lois a posteriori avec un prior non-informatif et en rouge, les lois a posteriori avec un prior informatif.

Figure 5 .

 5 Figure 5.10 -Histogrammes des lois a posteriori des paramètres du modèle de temps d'attente pour le jeu de données B. En bleu, les lois a posteriori avec un prior non-informatif et en rouge, les lois a posteriori avec un prior informatif.

  Figure 5.11 -Comparaison du score PE du quantile posterior 95% pour le jeu de données A . En bleu le score PE pour le modèle avec prior non-informatif et en rouge en avec prior informatif.

Figure 5 .

 5 Figure 5.12 -Comparaison du score PE du quantile posterior 95% pour le jeu de données B. En bleu le score PE pour le modèle avec prior non-informatif et en rouge en avec prior informatif.

Figure 5 .

 5 Figure 5.13 -Comparaison des lois a posteriori prédictives des temps d'attente pour le premier jour de prédiction pour les six créneaux du dataset A . La courbe en rouge est la loi prédictive lorsque le prior est informatif et en noir sa région HDI (Highest Density Interval). En violet, la région HDI de la loi prédictive des temps d'attente en non-informatif.

Figure 5 .

 5 Figure 5.14 -Comparaison des lois a posteriori prédictives des temps d'attente pour le premier jour de prédiction pour les six créneaux du dataset B. La courbe en rouge est la loi prédictive lorsque le prior est informatif et en noir sa région HDI (Highest Density Interval). En violet, la région HDI de la loi prédictive des temps d'attente en non-informatif.

  Ecov propose aux collectivités plusieurs solutions de mobilité comme par exemple la solution Lane https://www.lanemove.com/. Il s'agit d'un service de covoiturage structuré en ligne se situant au Sud-Est de Lyon. Afin d'améliorer l'expérience des passagers, la connaissance des flux de trafic et des temps d'attente est essentielle. Nous rappelons que tout au long de cette thèse nous avons développé des outils et des modèles dans ce but. L'objet de ce dernier chapitre est de montrer les différentes étapes et enjeux pour mettre en production le modèle hiérarchique bayésien développé dans le Chapitre 4 permettant la prédiction des flux de trafic et de temps d'attente.

  Le passage d'un modèle prédictif expérimental vers une version stable et industrielle de celle-ci nécessite plusieurs pratiques et outils informatiques. Durant la phase de conceptualisation du modèle, nous construisons un modèle prototype, dans ce cas, il fonctionne localement sur un Jupyter notebook[START_REF] Kluyver | Jupyter Notebooks-a publishing format for reproducible computational workflows[END_REF]. Ensuite, pour la phase d'industrialisation du modèle nous avons développé un paquet en Python autoportant de manière à ce qu'il soit exécuté automatiquement à des moments clés (tous les matins, à chaque demande de covoiturage). Nous avons cartographié à la Figure6.1 les outils informatiques principaux nécessaires pour les deux phases.

Figure 6

 6 Figure 6.1 -Outils principaux utilisés pour le prototypage du modèle (à gauche) et pour son industrialisation (à droite).

  4. Le fichier "setup.py" est indispensable pour la bonne construction et distribution d'un paquet Python. En effet, il est responsable de la transmission des instructions indispensables à la construction du paquet tels que le nom, la version et les dépendances. 5. Le dossier tests contient les fichiers Python avec les tests unitaires afin de vérifier le bon fonctionnement des différentes méthodes et fonctions. L'exécution du modèle est aussi testée sur des données artificielles dans le Jupyter Notebook "Tests.ipynb".

Figure 6

 6 Figure 6.2 -La structure du paquet bayesian-hierarchical-model, vue sous forme d'arborescence affichant les dossiers et fichiers présent du paquet.

  Figure 6.3 -Exemple de trafic de flux généré par la fonction TrafficFlowGenerator du paquet bayesian-hierarchical-model.

  f i t t i n g _ d f , p r e d i c t i o n _ d f ,W, W_test = model .d a t a _ p r e p a r a t i o n ( ) samples = model . build_model ( ) 5. Un DataFrame est une structure de données bidimensionnelle, les données sont alignées de façon tabulaire en lignes et en colonnes.

  Figure 6.5 -Exemple du contrôle prédictif a posteriori des temps d'attente simulés avec la distribution prédictive a posteriori des temps d'attente via la fonction graphical_ppc_hist() du paquet bayesian-hierarchical-model.

Figure 6 . 6 -

 66 Figure 6.6 -Exemple du contrôle prédictif a posteriori des temps d'attente simulés avec les quantiles de la distribution prédictive a posteriori des temps d'attente via la fonction graphical_ppc_boxplots() du paquet bayesian-hierarchicalmodel.

Figure 6 . 7 -

 67 Figure 6.7 -Exemple du graphe du score PE produit par la fonction plot_pe() du paquet bayesian-hierarchical-model sur des données simulées.

Figure 6 . 8 -

 68 Figure 6.8 -Schéma de la structure des ETL mis en place afin d'exposer en production les prédictions des flux de trafic et des temps d'attente aux différents services d'Ecov.

  7. https://docs.pymc.io/ CONCLUSION Synthèse L'objectif de cette thèse est de développer des outils méthodologiques destinés au domaine du covoiturage, et plus précisément, à la prévision des temps d'attente et au potentiel du covoiturage sur un territoire. Le premier thème, à savoir les temps d'attente, est abordé dans les Chapitres 2 et 4. Le thème du potentiel de covoiturage est couvert par le Chapitre 3. Le Chapitre 5 présente une méthodologie pour la construction d'un prior informatif afin d'améliorer les estimations des temps d'attente dans une situation de jeu de données court. Enfin, le Chapitre 6 détaille la mise en production et l'exploitation industrielle du modèle hiérarchique bayésien.

Titre:

  Potentiel et prévision des temps d'attente pour le covoiturage sur un territoire Mot clés : covoiturage, temps d'attente, régression quantile, modèle hiérarchique bayésien, prior informatif Resumé : Cette thèse s'intéresse au potentiel et à la prévision des temps d'attente concernant le covoiturage sur un territoire en utilisant des méthodes d'apprentissage statistique. Cinq thèmes principaux sont abordés dans le présent manuscrit. Le premier présente des techniques de régression quantile afin de prédire des temps d'attente. Le deuxième détaille la construction d'un processus de travail empruntant des outils des Systèmes d'Information Géographique (SIG) afin d'exploiter pleinement les données issues du covoiturage. Dans un troisième temps nous construisons un modèle hiérarchique bayésien en vue de prédire des flux de trafic et des temps d'attente. En quatrième partie nous proposons une méthode de construction d'une loi a priori informative par transfert bayésien dans le but d'améliorer les prédictions des temps d'attente pour une situation de jeu de données court. Enfin, le dernier thème se concentre sur la mise en production et l'exploitation industrielle du modèle hiérarchique bayésien.

Introduction Chapitre 5 : Construction d'un prior informatif par transfert bayésien: application au covoiturage À

  présent, nous souhaitons transférer l'apprentissage du modèle bayésien hiérarchique présenté dans le Chapitre 4 d'une ligne de covoiturage vers une nouvelle ligne. Ce transfert bayésien va permettre d'améliorer les performances du prédicteur sur le nouveau jeu de données qui est de petite taille. Disposer d'un modèle prédictif performant dès l'ouverture d'une nouvelle ligne de covoiturage est un enjeu stratégique pour assurer dès son départ son bon fonctionnement.

Table 2

 2 2.3. Comparison of Regression models for carpooling request waiting times

	Explanatory Variable	Description	Category
	origin_id	Identifier of physical meeting point (origin)	Geographical
	destination_cluster	Identifier of destination cluster Geographical
	origin_latitude	Latitude of origin	Geographical
	origin_longitude	Longitude of origin	Geographical
	destination_latitude	Latitude of destination	Geographical
	destination_longitude	Longitude of destination	Geographical
	travel_time	Duration of carpooling trip	Geographical
	travel_distance	Distance of carpooling trip	Geographical
	day_of_week	Day of week	Temporal
	time_of_day	Time of day (continuous)	Temporal
	is_peak_hour	Peak hour indicator	Temporal
	is_vacation_period	Holiday indicator	Temporal
	service_seniority	# days since registration	Passenger
	number_of_requests	# carpooling requests since registration	Passenger
	validated_requests	# completed carpooling trips since registration	Passenger

.1 -Summary table of the explanatory variables in the regression models.

Table 2

 2

	75					kth test fold				All
	Model	1	2	3	4	5	6	7	8	9	10
	QuantReg 8.17 7.97 7.43 12.05 8.70 10.46 9.84 8.20 10.07 8.98 9.19
	GBM	7.17 8.44 9.11 7.36 7.79 8.26 8.38 9.37 9.29 9.36 8.45
	QRF	8.29 6.25 7.54 8.80 8.93 9.145 7.10 9.17 9.45 8.75 8.34
	GRF	4.38 4.32 4.62 4.70 4.53 4.02 5.40 4.54 4.16 3.87 4.45
	α = 0.25
	QuantReg 4.43 5.40 4.79 4.69 4.38 5.35 4.02 5.43 4.87 4.64 4.80
	GBM	4.45 4.62 5.22 5.13 4.88 4.15 4.24 5.23 4.12 4.85 4.69
	QRF	5.11 4.61 4.73 5.21 5.45 4.56 4.04 5.89 4.99 4.45 4.90
	GRF	3.21 3.26 3.85 3.38 3.63 2.78 3.55 3.50 3.43 2.94 3.35

.2 -Pinball losses for all 10 test folds of the 10-fold cross validation for the 75% and 25% Quantile Regressions. The first column is the Regression model names, the middle columns are for the 10 test folds, and the last column is the overall generalised pinball loss with the lowest value in bold.

  Table2.3 -Ramp losses for all 10 test folds of the 10-fold cross validation for the 75% and 25% Quantile Regressions. The first column is the Regression model names, the middle columns are for the 10 test folds, and the last column is the overall generalised ramp loss with the closest value to 1 -α in bold.

										2.4. Conclusion
	α = 0.75					kth test fold				All
	Model	1	2	3	4	5	6	7	8	9	10
	QuantReg 29.23 35.38 35.94 12.50 23.44 23.44 12.50 37.50 21.87 26.56 25.84
	GBM	30.77 23.08 29.69 32.81 35.94 23.44 23.44 26.56 23.44 18.75 26.79
	QRF	32.31 38.46 29.69 26.56 29.69 26.56 29.69 29.69 25.00 28.12 29.58
	GRF	23.33 21.67 36.67 25.00 31.67 13.33 23.33 21.67 21.67 20.00 23.83
	α = 0.25
	QuantReg 75.38 76.92 81.25 64.06 75.00 67.19 68.75 79.69 71.87 76.56 73.67
	GBM	86.15 66.15 62.50 89.06 71.87 73.44 75.00 62.50 71.87 73.44 73.20
	QRF	69.23 72.31 67.19 68.75 62.50 68.75 71.87 60.94 75.00 71.87 68.84
	GRF	85.00 66.67 88.33 70.00 76.67 66.67 68.33 68.33 78.33 75.00 74.33

  t ∈ I s . If we observe a new daily driver flow ỹ on this future day, then the posterior predictive distribution of the waiting time ws is

		1
	p( ws |ỹ, y y y, w w w) =	0	p( ws |ỹ, β s )π(β s |y y y, w w w, t t t) dβ s ,	(4.8)

which is then collated into an S-vector (p( w1 |ỹ, y y y, w w w), . . . , p( wS |ỹ, y y y, w w w)) for all time intervals I 1 , . . . , I S .

Table 4 .

 4 2 -Poisson passenger arrival rates per hourly intervals. λ 1 is the base passenger arrival rate, λ 2 is an increase by 50% and λ 3 is an increase of 100%.

	1	8	6	4	6	4	4
	λ 2	12	9	6	9	6	6
	λ 3	16	12	8	12	8	8

  dθ .En statistique bayésienne, le prior permet d'ajouter à une modélisation une connaissance a prior, telle que la connaissance du métier ou la connaissance du terrain apportée par exemple par un expert. Le choix de cette loi a priori est important car il influence directement l'inférence du modèle. Lorsque nous n'avons pas de connaissance a priori nous avons recours à un prior dit noninformatif. Plusieurs priors non-informatif peuvent être construits, et nous Chapter 5 -Construction d'un prior informatif par transfert bayésien : application au covoiturage utilisons ici celuide Jeffreys, 1946. Ce prior a l'avantage d'être invariant lors de la reparamétrisation des coordonnées du paramètre. Il est défini par l'information de Fisher du modèle. Dans la communauté statistique, la définition des priors informatifs fait parfois l'objet de débats. En effet, les intuitions et les hypothèses portant sur un prior sont subjectifs. Il s'agit d'une question cruciale car un prior informatif bien construit permettra souvent d'éviter des modélisations aberrantes entraînant des conclusions erronées. Cette situation est d'autant plus sensible dans un cas de jeu de données de petite taille(voir Vanpaemel, 2011).La construction des priors informatifs est un sujet d'intérêt pour les applications industrielles. Dans le domaine de l'énergie,[START_REF] Launay | Construction of an informative hierarchical prior for a small sample with the help of historical data and application to electricity load forecasting[END_REF] ont construit un prior informatif pour améliorer les prédictions de consommation d'énergie des particuliers. Le prior développé dans cet article est construit sur un jeu de données long afin d'être ultérieurement appliqué sur un jeu de données court. De même,

	π(θ) ∝ det	1 2 I(θ)	(5.1)
	avec I,

  Le paramètre η représente le coefficient de transition du type de jour qui précède le jour i sur le flux y i .Le deuxième étage modélise les temps d'attente par créneaux dans la journée. Il est strictement identique au modèle présenté au Chapitre 4. On rappelle qu'il s'agit d'un modèle de régression Gamma qui permet de modéliser les flux journaliers en fonction des temps d'attente par créneaux. Pour le jour i, nous avons n i temps d'attente qu'on note w i,1 , . . . , w i,n i pour des demandes de covoiturages effectuées aux instants t 1,i < • • • < t 1,n i . L'intervalle s de la journée est noté I s . On suppose que les temps d'attente vérifient

				5)
	avec ε l'erreur et DT une fonction qui caractérise les types de jours,
	DT(i) =	  	ORD si le jour i est un jour ordinaire, SCH si le jour i est un jour férié scolaire.	(5.6)

  5 pour perturber les paramètres initiaux. Le récapitulatif des paramètres de simulation pour chaque jeu de données est détaillé par la Table5.1.

	Jeu de données ν	β 1	β 2	β 3	β 4	β 5	β 6 Nombre de jours simulés
	L	10 0.12 0.1 0.11 0.13 0.18 0.16	390
	C	15 0.24 0.2 0.22 0.26 0.36 0.32	90
	Table 5.1 -Tableau récapitulatif des paramètres de simulation pour le jeu de
	données long L et court C.

Table 5 .

 5 2 -Comparaison des écarts-types des lois a posteriori des paramètres du modèle informatif et non-informatif.

  wi ) est le quantile au niveau α de la loi prédictive a posteriori de la i ième observation et w i le i ième temps d'attente observé pour un échantillon de taille Ñ . Le niveau du quantile choisi pour notre comparaison est α = 95%. Ce niveau avait été retenu lors de l'industrialisation du modèle parce qu'il permettait d'assurer un bon niveau de service à nos utilisateurs. La Figure5.11 compare les scores PE du quantile niveau 95% issues du modèle informatif et non-informatif pour le jeu de données A . La Figure5.12 représente la même comparaison mais pour le jeu de données B. Nous constatons que dans le cas

du transfert de A vers A comme dans le cas du transfert A vers B le score PE du modèle informatif est meilleure. Ce résultat est d'autant plus visible sur le jeu de données B ce qui montre que la loi a priori informative apporte une réelle covoiturage

  tout d'abord que les données utilisées par le modèle proviennent de deux sources distinctes i) les traces GPS des conducteurs et ii) les temps d'attente observés des passagers. Nous rappelons aussi que le modèle développé pour la prédiction des flux de trafic et des temps d'attente est composé de deux étages imbriqués. Comme pour le Chapitre 5, nous travaillons sur un modèle dérivé de celui présenté au Chapitre 4. Le premier étage modélise les flux quotidiens sur une ligne de covoiturage. Ce premier étage est défini par le modèle

  Le paramètre η représente le coefficient de transition du type de jour qui précède le jour i sur le flux y i .Le deuxième étage modélise les temps d'attente par créneaux dans la journée. Il est strictement identique au modèle présenté au Chapitre 4. On rappelle qu'il s'agit d'un modèle de régression Gamma modélisant les flux journaliers en fonction des temps d'attente par créneaux. Pour le jour i, nous avons n i temps d'attente qu'on note w i,1 , . . . , w i,n i pour des demandes de covoiturages effectuées aux instants t 1,i < • • • < t 1,n i . L'intervalle s de la journée est noté I s . On suppose que les temps d'attente vérifient w

	  	ORD si le jour i est un jour ordinaire, SCH si le jour i est un jour férié scolaire.	(6.2)

1) avec ε l'erreur et DT une fonction qui caractérise les types de jours, DT(i) =

  6.3. Industrialisation du modèle hiérarchique bayésienpar son arborescence en dossier à la Figure6.2, est divisée en cinq parties.1. Les fichiers "README.rst" et "README_files" contiennent des informations importantes sur le projet pour les futurs utilisateurs. Par convention, il est placé au tout premier niveau de l'arborescence pour que les utilisateurs trouvent immédiatement les informations.2. Le dossier bayeshierarchical contient les fichiers Python et Stan qui correspondent aux différents modules du modèle. Le fichier "__init__.py" permet l'interprétation des répertoires en tant que modules, il est toujours chargé en premier dans un module.3. Le dossier doc contient tous les fichiers textes qui constituent la documentation du paquet. Ce répertoire est généré automatiquement par l'outil Sphinx.

Nous comparons ensuite les performances de ces quatre algorithmes de régression quantile sur les données d'Ecov. En raison de la complexité du service de covoiturage tel qu'il est proposé par Ecov au cours de cette phase expérimentale et grâce à cette analyse de régression, Ecov a ensuite pivoté son service de covoiturage pour desservir un nombre limité de destinations, plus proches des lignes de bus. Ecov a également intégré une technologie permettant de collecter les traces GPS des itinéraires des conducteurs à partir de leur application mobile afin de mieux comprendre leurs comportements. Dans les chapitres suivants, nous analyserons les données de ce service amélioré et la manière dont il facilite une meilleure compréhension du potentiel de covoiturage.

Ce chapitre fait l'objet d'une publication en cours co-écrite avec Anne Philippe (Professeure à l'Université de Nantes), Bertrand Michel (Professeur à l'Université de Nantes et l'École Centrale de Nantes) et Tarn Duong (Docteur en Statistique et membre de l'équipe Data/SIG d'Ecov).

Nous souhaitons à présent exploiter de manière industrielle les travaux exposés dans les chapitres précédents. En d'autres termes, l'objet de ce chapitre est de présenter les différentes étapes pour la mise en production du modèle hiérarchique bayésien introduit dans le Chapitre 4.La première étape est la construction du paquet Python nommé bayesianhierarchical-model reprenant les outils et les algorithmes développés dans ce manuscrit. Le développement de ce paquet a nécessité d'assurer le passage d'un modèle prédictif expérimental vers une version stable et industrielle de celui-ci. À ce titre, nous avons utilisé différents outils informatiques, certains d'entre eux sont dédiés pour la phase de prototypage du modèle et d'autres pour la phase d'industrialisation. Nous avons cartographié à la Figure1.25 les outils informatiques choisis durant cette thèse pour les deux phases.

https://airflow.apache.org/

www.ecov.fr

2.2. Quantile Regression

https://moovit.com/blog/50-million-global-transit-report/

On construit ensuite le prior informatif comme détaillé en Section 5.2 pour le jeu de données A .

https://airflow.apache.org/

Acknowledgements

The authors thank Ecov for providing the data sets of driver GPS traces and passenger waiting times. The authors also thank Bertrand Michel from the Central Engineering School of Nantes and Gérard Biau from Sorbonne University for their feedback.

Acknowledgements

The authors thank Ecov for providing the data sets of the driver GPS traces and the passenger waiting times. The authors also thank Safa Fennia, Madeleine Zuber, Flavien Sindou and Constant Bridon from Ecov, and Gérard Biau from Sorbonne University for their feedback.

Chapter 4 -Bayesian hierarchical models for the prediction of the driver flow and passenger waiting times in a stochastic carpooling service Figure 4.5 -Predictions of simulated pseudo waiting times for 3-hourly intervals, for all Ñ prediction days. The observed waiting times are the grey box plots, and the 50%, 75%, 95% quantiles of the posterior predicted waiting times are the light, medium and dark purple circles.

during both the training phase PE(W, W; δ) (blue curve) and the test phase PE(W test , Wtest ; δ) (red curve). The test predictions are more accurate than the training predictions for small values of δ < 4 minutes since red PE curve is below the blue PE curve in this interval. This reverses for δ greater than 4 minutes, and after 12 minutes, both curves level off at 1. Thus the posterior predictions from our proposed Bayesian hierarchical model can have robust prediction performance according to this PE metric. Chapter 4 -Bayesian hierarchical models for the prediction of the driver flow and passenger waiting times in a stochastic carpooling service weeks, and the fourth column is the number of training days (N ). For these training-test scenarios, in addition to our proposed Bayesian hierarchical multilevel (BHML) predictions, we compute predictions from a baseline frequentist model (BASE), and a Bayesian Prophet model (PROP). The details of these competing models are described in Appendix 4.7.3. We input the daily driver flows into the first hierarchical model from the Bayesian hierarchical multi-level model BHML to produce the posterior predicted daily driver flows ỹi , as well the corresponding predictions/estimations from the frequentist baseline model BASE and the Bayesian Prophet model PROP. For the Test scenario # 6, the training phase covers the dates 2018-05-15 to 2019-05-19. In Figure 4.8 is the evolution of the goodness-of-fit of the three different models for daily driver flow estimation (leaving out the first week 2018-05-15 to 2018-05-21 which serves as the 'burn-in' period). The goodness-offit is measured by the MSE of the estimated and the observed daily driver flows, aggregated per week. Visually the BHML tends to have the best goodness-offit (smallest MSE) for most weeks. The sum of these weekly MSEs are BASE: 421.9, PROP: 816.9, BHML: 297.2, which confirms our visual impression that the BHML achieves the best overall estimation accuracy.

Therefore we can be confident that the Bayesian hierarchical multi-level moving average model has good estimation accuracy/goodness-of-fit, but this good performance does not necessarily translate to prediction [START_REF] Makridakis | Forecasting in social settings: the state of the art[END_REF]. So for each scenario described in Table 4.1, we compute the BHML, BASE and PROP models for the training phase, and then the days of the test phase are input into each these training models to yield the daily driver flow predic-Chapter 4 -Bayesian hierarchical models for the prediction of the driver flow and passenger waiting times in a stochastic carpooling service Now that we have verified that a Gamma regression model is suitable for the data observed in the Lane carpooling service, we proceed with the BHML to form predictions of the passenger pseudo waiting times. Since there are insufficient passenger carpooling requests to robustly compute observed hourly waiting time profiles over an entire day for the school holidays (SCH) and the public holidays/weekends (PWE), we restrict ourselves to forming predictions for the weekdays (ORD). In Figure 4.14 are the box plots of the weekly number of observed pseudo waiting times for each hourly interval for the weekdays passenger waiting times in a stochastic carpooling service the BHML is that it gives a more solid basis that this performance will continue into the future. For the service operator, the BHML implies that the key factor in determining the passenger waiting time (i.e. the output from the second stage) is the driver flow (i.e. the output from the first stage). It is thus imperative that a consistent level of driver participation in the carpooling service is maintained so that

Appendix

GPS traces pre-processing

A GPS trace is an -sequence of triplets X X X = {(X i , Y i , T i )} i=1 where (X i , Y i ) are the longitude, latitude coordinates of the GPS sensor at the i th timestamp T i . The m pick-up/drop-off locations in the carpooling network are represented by their

is the set of days with the same day of week before day i; T HOL (i) = {k : k < i, DT (i -k) = HOL} is the set of school/public holidays before day i; and DN is the day of week number function,

The Bayesian Prophet model, devised by [START_REF] Taylor | Forecasting at scale[END_REF]Facebook Core Data Science Group, 2019, is an additive model with three components:

where g(i) is the trend, s(i) is the seasonality, and h(i) is the holiday effect. The linear trend is

where k is the growth rate, m is the offset, a a a is the change point indicator, δ δ δ is the growth rate adjustment, and γ γ γ is the piece-wise continuity adjustment to ensure that g is continuous. The seasonality component is a Fourier decomposition s(i) = L =1 [α cos(2π i/P ) + β sin(2π i/P )] where (α , β ) are the Fourier coefficients, L is the number of Fourier coefficients and P is the period (in days). The holiday effect is h(i) = h h h(i) κ κ κ where, say, h h h(i) = (1{DT(i) = SCH}, 1{DT(i) = PWE}) is the vector of indicator variables of the type of holiday of day i, and κ κ κ is the weight vector, usually equal to the all-ones vector. [START_REF] Taylor | Forecasting at scale[END_REF] provide the details for the construction of the change point function a a a(t) and the continuity adjustment parameter γ γ γ. These authors set the number of Fourier coefficients to be L = 10 for yearly cycles and L = 3 for weekly cycles. What remains is to estimate the trend growth rate k, the offset m, the growth rate adjustments δ δ δ and the Fourier coefficients α α α.

Chapter 6 -Valorisation Industrielle : Mise en production du modèle hiérarchique bayésien

• T : Le nettoyage et la correction des données brutes, par exemple les temps d'attente négatifs et la conversion au bon fuseau horaire des horodatages. Une fois les données corrigées les deux tables sont fusionnées en une seule.

• L : L'enregistrement de cette table sur le data_warehouse d'Ecov.

ETL2:

• E : La récupération de la table produite par l'ETL1.

• T : La transformation des données afin de coïncider au format de données d'entrée du modèle comme illustré par la Figure 6.4 et l'entraînement du modèle avec la méthode ".build_model()".

• L : L'enregistrement des lois a posteriori des paramètres du modèle dans le data_warehouse d'Ecov.

ETL3:

• E : La récupération des lois a posteriori des paramètres du modèle par l'ETL2.

• T : La génération des lois a posteriori prédictives des flux de trafic et des temps d'attente basés sur les lois a posteriori des paramètres du modèle. Ensuite le calcul des quantiles prédictifs a posteriori qui nous intéresse.