À proximité de l'aéroport d'Orly un ami se plaignait que les appareils ne respectaient pas toujours les couloirs aériens, connaissant la possibilité de recevoir la position des appareils je me suis amusé à créer une application qui relevait la présence d'un aéronef dans une zone définie à une altitude minimale définie. Cela a constitué un bon apprentissage des données flux sous Spark.

Vols aéronefs

Selon la définition de Air France :
L'ACARS (Aircraft Communication Addressing and Reporting System) est un système permettant l'échange d'informations (messages) entre l'avion et le sol sous forme numérique codée par liaison radio ou satellite.

Pour capter ce flux de données qui émettent sur 1090 MHz il est possible d'utiliser une simple clé usb qui reçoit la TNT et d'y ajouter une petite antenne spécifique.

Le site Flight24 de manière plus professionnel utilise essentiellement les flux de données de type ACARS.

Logiciel permettant le traitement du flux

Les données brutes de données de type ACARS ne sont pas directement utilisables, elles nécessitent un traitement complexe réalisé par le logiciel freeware dump1090.

Un exemple des données de sortie de dump1090

MSG,4,111,11111,400804,111111,2020/05/16,15:57:34.670,2020/05/16,15:57:34.665,,,390,326,,,64,,,,,0
MSG,6,111,11111,406091,111111,2020/05/16,15:57:34.719,2020/05/16,15:57:34.670,,,,,,,,2041,0,0,0,0
MSG,8,111,11111,800734,111111,2020/05/16,15:57:34.752,2020/05/16,15:57:34.732,,,,,,,,,,,,0
MSG,8,111,11111,AE4BE6,111111,2020/05/16,15:57:34.752,2020/05/16,15:57:34.732,,,,,,,,,,,,0
MSG,7,111,11111,800734,111111,2020/05/16,15:57:34.765,2020/05/16,15:57:34.734,ka,37000,,,,,,,,,,0
MSG,8,111,11111,AE4BE6,111111,2020/05/16,15:57:34.768,2020/05/16,15:57:34.735,,,,,,,,,,,,0
MSG,8,111,11111,400804,111111,2020/05/16,15:57:34.792,2020/05/16,15:57:34.737,me,,,,,,,,,,,0
MSG,8,111,11111,AE4BE6,111111,2020/05/16,15:57:34.792,2020/05/16,15:57:34.737,,,,,,,,,,,,0
MSG,8,111,11111,400804,111111,2020/05/16,15:57:34.800,2020/05/16,15:57:34.796

Explication des variables de dump1090


Sur ces 22 variables, nous allons en retenir 5 pour réaliser le filtrage:

• icao : Numéro d'identification de l'appareil indépendant de la compagnie
• callsign : Identifiant de l'appareil lié à la compagnie
• altitude : Altitude
• latitude : Latitude
• longitude : Longitude

Spark streaming

Dans les multiples possibilités offertes par Spark il est possible de traiter des données streaming en temps réels.
L'outil dump1090 permet de transmettre les données Acars en mode réseau, Spark peut récupérer ce type de flux (socketTextStream) et transformer ces données pour réaliser un traitement.


L'objectif initial est de définir une zone de contrôle et une altitude minimale à ne pas dépasser.

Pour créer le programme de traitement, il y aura trois étapes :

1. Réaliser le programme avec un jeu de données fixes afin de réaliser des tests

• La zone de test
  Elle est située au nord de Paris, latitude strictement supérieure à 48,903296 et strictement inférieure à 48,924444, longitude strictement supérieure à 2,331250 et strictement inférieure à 2,385399, l'altitude est strictement inférieure à 15000 pieds.
  
• Avec l'outil dump1090 dans un parc parisien j'ai réalisé la captation des données sur une période de 28 minutes, 103 435 lignes soit un peu moins de 4 000 lignes par minutes avec une antenne artisanale. Lien du fichier
• On traite les données afin qu'elles correspondent à la zone de test, on réalise différentes jointures afin de récupérer le icao et le callsign et avoir ces données sur une même ligne, car les données acars sont émises sous formes de packet de données différents.
  
  

2. Refaire le programme avec des données temps réels

• Le traitement des données est réalisé toutes les 5 minutes c'est-à-dire que pendant 5 minutes les données Acars s'écrivent dans Spark, le traitement et le filtrage sont réalisé après 5 minutes (StreamingContext(sc, Minutes(5)))
• Pour éviter qu'un échec de l'application entraine une perte des données, des checkpoints sont mis en places dans un environnement de données sécurisées de type hdfs Hadoop
• À noter que dans le programme il y a une simplification, la date du traitement du programme est la date de l'événement
• Le résultat du traitement est copier dans un disque hdfs sous le format suivant :
  (4CAA7B,AZA312,7425.0,48.91099,2.33134,2020/05/16,16:14:19.638)
  (440452,EJU89ZB,8775.0,48.9082,2.33269,2020/05/16,16:23:44.883)
  (440452,EJU89ZB,8775.0,48.90811,2.33149,2020/05/16,16:23:45.433)
  (392AEE,AFR1889,8300.0,48.90946,2.33734,2020/05/16,16:16:03.400)
  (392AEE,AFR1889,8300.0,48.90939,2.33636,2020/05/16,16:16:03.841)
  icao,callsign,altitude, longitude, latitude, date

3. Il faut créer un exécutable de l'application
Pour Spark sous scala nous utiliserons SBT.

• Lien de l'éxecutable en html
• Lien de l'éxecutable avec en plus la configuration SBT

En France chaque année il y a plus de 10 000 décès liés à un mauvais usage des médicaments c'est-à-dire "un mauvais dosage, mauvaise prise, non-respect du traitement prescrit, interaction entre plusieurs médicaments".

Dans cette étude nous nous intéresserons au problème lié à l'interaction entre plusieurs médicaments. Dans le cadre d'un mémoire de fin d'étude de médecine, j'ai apporté ma contribution dans la construction de modèles explicatifs. Cette contribution ne fut pas vaine, car l'étudiante est devenue médecin lors de sa soutenance avec les féliciations du jury.

Récupération de données

Cette étude met le focus sur les médecins. Nous cherchons à déterminer les facteurs qui permettent d'expliquer le niveau de connaissance des risques liés à la posologie.

À partir d'une base de données de médecin, un sondage en ligne a été réalisé et a permis de récupérer près de 5656 participants. N'étant pas propriétaire des données ces dernières ne sont pas en ligne.

Les variables utilisées pour la réalisation du modèle

Variables à prédire

En fonction de questions techniques liées aux interactions médicamenteuses, deux notes ou scores sont calculés.
Le score de général va de 0 à 9 sachant que le meilleur score est de 9.
Le score de gravité va de 0 à 5 sachant que le meilleur score est de 5.

Variables explicatives

Choix du type de modèle

Il ne s’agit pas de trouver le modèle qui permet de prévoir si un score est de 1 ou de 5. Ici les variables à
expliquer respectent un ordre, il est plus important de rechercher les déterminants qui font croitre le score c’est pourquoi nous utiliserons des modèles de régression logistique de type polytomique avec des données ordinales avec une distribution de probabilité cumulative.

Soit la variable Y : score sge ou sgr, j : niveau du score
Une probabilité cumulative pour Y est la probabilité que Y soit au-delà ou en deçà d’un niveau j.
Dans notre cas nous étudierons le cas de la probabilité cumulative pour que la probabilité que Y soit au-delà d’un niveau j c’est-à-dire que le score augmente.

Pour faciliter l’étude des variables explicatives, on réalise les modèles avec l’hypothèse des odds
proportionnels. (voir Score Test for the Proportional Odds Assumption sortie SAS).

Pour la variable explicative score sge ou sgr, le modèle est le suivant :
$$ logit [P(Y \ge j)] = \exp^{(\alpha_{j}+\beta_{x})} \ avec \ j = 1, ...., j-1 $$
avec le paramètre β décrivant les effets de x sur log odds of response de la catégorie j. Puisque nous sommes
dans le cadre d’odds proportionnels nous admettons que l’effet de x est identique pour tous les j-1 logit
cumulatif.

Pourquoi j-1 ?

Exemple pour le premier score
6 modalités et 5 courbes car la probabilité d'une des courbes = 1

La forme de la courbe

$$ P(Y \ge j) = \frac{\exp^{(\alpha_{j}+\beta_{x})}}{1+\exp^{(\alpha_{j}+\beta_{x})}} $$

α j , valeur des intercepts qui définit la position des courbes logit relatives aux scores avec j étant le niveau du score.
Les coefficients intercepts sont différents pour chaque courbe, plus j augmente plus on se déplace vers la gauche car les coefficients intercepts baissent.

βx , paramètre des variables explicatives x
L'intérêt des odds proportionnels est que pour les courbes logit ont les mêmes paramètres βx,seul change la
valeur de α pour chaque j.

Les coefficients de la courbe sont identiques pour les 5 courbes car nous avons testé l’hypothèse d’égalité des pentes (Test for the Proportional Odds Assumption) sinon pour chaque courbe nous aurions des coefficients différents et il ne serait pas possible d'avoir des odds ratios pour l'ensemble des données.

Exemple d'une des sorties
P(y>=0) = 1 pas de courbe, probabilité = 1
P(y>=1) : intercept 0.09
P(y>=2) : intercept -1.27
P(y>=3) : intercept -2.34
P(y>=4) : intercept -3.67
P(y>=5) : intercept -5.47

Les modèles sélectionnés

Modélisation d'avoir un score de gravité élevé (le plus élevé est 5)

Les sorties sas pour le modèle du score de gravité

Modèle 7 Analysis of Maximum Likehood Estimates
• le fait d'être plus âgé augmente le logit de 0.0235 et ainsi augmente la probabilité d'avoir un score de gravité plus élevé
• le fait de connaître l'existence du thésaurus présent sur le site de l'ansm augmente le logit de 0.4816 et ainsi augmente la probabilité d'avoir un score de gravité plus élevé
• le fait d'avoir fait ses études dans une université non parisienne augmente le logit de 0.2134 et ainsi augmente la probabilité d'avoir un score de gravité plus élevé
• le fait d'utiliser l'outil ansm augmente le logit de 0.8162 et ainsi augmente la probabilité d'avoir un score de gravité plus élevé

Rapport de chance relative (odds ratio)
• l'odds ratio de l'expérience est de 1.024 soit supérieur à un c'est à dire que si l'expérience augmente d'une unité la probabilité d'avoir un score de gravité plus élevé augmente, l'effet est faible mais est présent.
• lorsqu'un individu connaît l'existence du thésaurus présent sur le site de l'ansm nous avons 1,619 fois plus de chance d'avoir un niveau de gravité plus élevé par rapport à un individu qui n'en a pas sa connaissance.
• lorsque les études ont été réalisées dans une université non parisienne nous avons 1.532 plus de chance d'avoir un niveau de gravité plus élevé par rapport à un individu qui a réalisé ses études dans une université parisienne.
• lorsqu'un individu utilise l'outil ansm nous avons 2,262 fois plus de chance d'avoir un niveau de gravité plus élevé par rapport à un individu qui n'utilise pas ce dernier.

Modélisation d'avoir un score général élevé (le plus élevé est 10)

Les sorties sas pour le modèle du score général

Modèle 8 Analysis of Maximum Likehood Estimates
• le fait de connaître l'existence du thésaurus présent sur le site de l'ansm augmente le logit de 0.5971 et ainsi augmente la probabilité d'avoir un score général plus élevé
• le fait d'être de sexe masculin augmente le logit de 0.3559 et ainsi augmente la probabilité d'avoir un score général plus élevé

Rapport de chance relatives (odds ratio)
• lorsqu'un individu connaît l'existence du thésaurus présent sur le site de l'ansm nous avons 1.817 fois plus de chance
d'avoir un niveau général plus élevé par rapport à un individu qui n'en a pas sa connaissance.
• lorsque l'individu est de sexe masculin nous avons 1.427 fois de plus de chance d'avoir un niveau général plus élevé par
rapport à un individu de sexe féminin.

Pour le modèle de score de gravité nous avons 60,1% de paires concordantes, 38,1% de paires discordantes et 1.8% de paires égales.

Pour le modèle de score de général nous avons 57,1% de paires concordantes, 39,7% de paires discordantes et 3.3% de paires égales.


Autres

Détail de la signification des pourcentages de paires

Code Sas

p180 6.2.1 Dans l'exemple du livre l'ordre est différent avec P(Y<=j), dans nos modèles on a modélisé la probabilité d'être à un niveau plus élevé

Agreti introduction to categorical data

Les bases de données Nosql proposent une maintenance simplifiée et la possibilité de dimensionner la puissance de calcul en fonction de l’importance du flux.
Pour avoir un ordre d’idée de l’importance des flux des données temporelles, pour rappel une seconde est 1 000 millisecondes, une milliseconde est 1 000 microsecondes et une microseconde est 1 000 nanosecondes, c’est-à-dire que pour une seconde nous avons un million de microsecondes et qu’en une heure nous avons 3.6e+9 microsecondes.

D’après Wikipédia, la vitesse de transaction pour le trading haute fréquence était de 20 millisecondes à la fin des années 2010 pour passer en 2011 à 113 microsecondes, ce changement d’échelle nécessite une faible latence que les bases de données NoSql propose.

Format des dates

1. Présentation des données

Dans cette partie nous étudierons l'injection de données d'une base de plusieurs millions d'enregistrements.
https://archive.ics.uci.edu/ml/datasets/Heterogeneity+Activity+Recognition

Exemple d'un document du fichier csv
Index,Arrival_Time,Creation_Time,x,y,z,User,Model,Device,gt
0,1424696633909,1424696631914042029,0.013748169,-0.0006256103500000001,-
0.023376465,a,nexus4,nexus4_1,stand

Dans cette base nous utiliserons les données du gyroscope de téléphone portable, fichier
Phones_gyroscope.csv, taille de 1.3 Go.
Les données des gyroscopes sont x, y, z.
Pour chaque enregistrement nous avons la date en millisecondes de type Unix de la mesure,
le type d'appareil utilisé, le numéro de l'appareil, le nom de l'utilisation et le type d'activité :
"Biking", "Sitting", "Standing", "Walking", "Stair Up", "Stair down".

Ce type de données provenant de capteurs peut être utilisé pour modéliser le type d'activité c'est-à-dire qu'avec seulement les données des capteurs il serait possible de connaître le type d'activité en cours.
Dans notre cas nous utiliserons ce fichier de presque 14 millions d'enregistrements pour tester la vitesse d'indexation d'Elasticsearch et de influxdb.

2. Elasticsearch

Détail

3. Influxdb

Détail

Conclusion

L'environnement d'Elasticsearch permet relativement simplement de récupérer un flux de données et de créer un tableau de bord. J'ai été séduit par la possibilité de créer des index selon la date ce qui permet de supprimer facilement des données obsolètes sans bloquer la base de données avec des opérations de sauvegardes ce qui réduit le risque de corrompre ses données.

Le langage de requête d'Elasticsearch gagnerait à être plus lisible, l'utilisation du format json oblige à gérer un grand nombre d'accolades et augmente le risque d'erreur, à la différence du sql une requête complexe peut prendre plusieurs dizaines de lignes, les différentes documentations traitent peu de ce type de requête. Il faudrait créer un outil graphique pour rendre l'outil plus accessible, cependant on ne peut retirer à Elasticsearch la rapidité d'exécution des requêtes, souvent moins de 200 millisecondes avec les données que j'ai utilisées.

Au sujet de la vitesse d'indexation des données j'ai été impressionné par la vitesse d'indexation de influxdb par rapport à Elasticsearch cependant je ne peux conclure si un outil est mieux qu'un autre n'ayant testé que la partie indexation.

Pour des données massives, il faudrait en plus de la vitesse d'indexation tester la scalabilité des deux outils, selon une étude du Cern il y a quelques années elasticsearch brillait par sa capacité à la "scalabilité" par rapport à influxdb depuis il y a eu de nombreuses nouvelles versions de ces deux outils une nouvelle étude serait intéressante.

La meilleure manière de comprendre un algorithme est de le programmer, dans le cas suivant je vais implémenter un algorithme de type kmeans. Pour tester ce dernier j'utilise des données fournis par le département informatique de l'université de l'Est de la Finlande.

Mon choix se porte sur le dataset s1 : deux dimensions, 5 000 points avec 15 clusters.

Clustering basic benchmark

L'algorithme est réalisé sur matlab.

Description simplifiée de l'implémentation de l'algorithme

• définition de l'espace du nuage des points (max, min)

• tirage aléatoire de k points correspondant au k groupe que l'on recherche dans l'espace du nuage des points

Boucle

• calcul pour chaque point de la distance de ces derniers avec les k centres.

• rattachement des points au groupe dont le centre est le plus proche

• calcul du total T des distances entre chacun des points et leur centre

• définition les nouveaux centres de chaque groupe qui sont les barycentres ou centres de gravité de chacun de groupes. (moyenne de la position des points du nuage)

• Arrêt de boucle de l'algorithme lorsque le total T n'évoluent plus c'est à dire que chacun des groupes a un centre de gravité stable

Résultat
Nous avons 15 groupes. Pour avoir une aussi bonne classification, un script a été créé qui a lancé ce programme un très grand nombre de fois et conservé la classification qui avait le total T le plus faible.

Interprétation des résultats en utilisant une analyse graphique de type silhouette

La description de cette méthodologie est issue de l'article de Peter J. Rousseeuw, Silhouettes:a graphical aid to the interpretation and validation of cluster analysis.

On considère habituellement la distance euclidienne.
Soit i, un membre du cluster
Soit a(i), la moyenne de la distance du membre i du cluster par rapport à l'ensemble des membres du cluster auquel il appartient

Soit d(i,k) : la moyenne de la distance du membre i du cluster auquel i appartient et l'ensemble des membres du cluster k.

A partir de d(i,k) on calcule l'indice de voisinage b(i), on recherche la distance d(i,k) la plus faible qui nous donne une indication de la proximité du membre i par rapport à un autre cluster.

Avec ces deux moyennes on peut construire un indice de similarité s(i):

• Dans le premiere cas où a(i)<b(i), nous avons la distance moyenne mininal du membre i par rapport aux autres clusters qui est supérieure à la distance moyenne des membres auquel appartient i, c'est à dire que ce membre n'a pas de cluster auquel il serait plus optimal d'appartenir. Cette situation correspond à un bon classement et à un indice de similarité positif. La meilleure des situations étant le cas où l'indice de similarité est proche de 1.

• Dans le deuxième cas où a(i)=b(i), l'indice de similarité est nulle c'est à dire qu'il existe un ou plusieurs clusters où le membre i pourrait aussi appartenir.

• Dans le troisième cas où a(i)>b(i), nous avons la distance moyenne minimal du membre i par rapport aux autres clusters qui est inférieure à la distance moyenne des membres auqel appartient i, c'est à dire que ce membre n'est pas dans le cluster le plus optimal. Cette situation correspond à un mauvais classement et à un indice de similarité négatif. La plus mauvaise des situations étant le cas où l'indice de similarité est proche de -1.

Résulat répresentation graphique silhouette : fonction inclue dans matlab

A noter que pour chaque cluster, les indices de similarité des membres i sont classés par ordre décroissant.
Cette réprésentation graphique nous permet de constater que le nombre de cluster est conforme et que le classement est bonne qualité, la majorité des membres de chaque cluster a un indice de similarité supérieur à 0.75.

Code Matlab

Each year, Facebook using Kaggle platform to find their future data scientist.

The best competitor participants could pass a job interview. The interest of this contest is the members of kaggle couldn't exchange trick and models, its imply that your rating is the real rating and not a copy/paste program, method process.

We must find if a bidder is a human or a robot with bidding data provided by Facebook.
The score is an area under the ROC curve which quantifies the overall ability of the test to discriminate between human and robot.

At the end of the competition, my score is a 0.92, near to the best score of 0.94. My rank is 198 of 989, not bad.

Kaggle link

1. Data
We have two sets of data :

• Bidder dataset : bidder_id, robot flag
  Training dataset :
  2 013 distinct bidders with 103 robots (5,11%)

  Testing dataset :
  4 700 distinct bidders for testing dataset

• Bid dataset : bid_id, 7 656 335 transactions
  We have 7.6 millions of bid_id, related to bidder_id

I group Bid dataset by bidder_id, and create engineering features.


2. Feature engineering
Time
The time is not in hours/minutes/seconds, we have to explore this feature, to simplified time feature, I divide this time by 100 000 000 000 to control the number of time value.

If we plot time value against bid_id value, we find three not joined periods.

• Period 1: 96319 - 96455 : 136 distinct units of new time
• Period 2: 96955 - 97092 : 137 distinct units of new time
• Period 3: 97592 - 97728 : 136 distinct units of new time

I decide to create this new features :

By bidder_id, for each period :

• Number of bid
• Number of bid at same times
• Standard deviation between bid

Other features

Number of merchandise by bidder_id
Count distinct device by bidder_id
Number of country by bidder_id
Count distinct IP by bidder_id
Count distinct URL by bidder_id

3. Overfitting

Before the end of the contest, when we have found a model, we could send our prediction to kaggle which compute your score from only 30% of data, the final score is computed with 100% of data.

If for the same model with different submission we don’t have a stable score, maybe our model is unable to generalize.

My first models were clearly overfit : my training score was near 0.99, for gradient boosting and random forest with a testing score of 0.92, when I was testing on kaggle the score was near 0.89.

To limit overfitting, I have tried to reduce the difference between training score and testing score, and lower training score from 0.99 to 0.95-0.92.

I have made ensemble models for each models, and an ensemble models of ensemble models for the final prediction.


4. Models

Neural Network :
I am using the new python libraries, https://github.com/aigamedev/scikit-neuralnetwork .

3 layers :
Maxout : units=500, pieces=2 ==> Sigmoid: units=10 ==> Softmax
The validation size is 0.4

To have a validation error near to 0.16, we must have training score greater than 0.89. I have created an iterative process which select only these models, and made an ensemble of these models. It was a long processing, more than one day to create an ensemble of 5 models because 90% of models don't respect criteria.

Gradient boosting :
I am using scikit learn libraries, http://scikit-learn.org/stable/

My model is a gradient boosting with calibrated classifier cross validation.
"When performing classification you often want not only to predict the class label, but also obtain a probability of the respective label. This probability gives you some kind of confidence on the prediction. Some models can give you poor estimates of the class probabilities and some even do not not support probability prediction. The calibration module allows you to better calibrate the probabilities of a given model, or to add support for probability prediction"

GBM:
Number estimators : 500
Max depth : 6
Learning rate: 0.001
Max features : sqrt

I have created an ensemble of calibrated classifier. For each calibrated classifier, I choose a cross validation of 10 and I repeat the procedure of calibrated classifier 100 times to get more precise result.

5. Ensemble models

To find the weight of GBM and NN which maximize the area under the ROC curve, I have used a logistic regression.

0.25Neural Network + 0.85Gradient Boosting

Code Python :

• Data Management 1/2
• Data Management 2/2
• Neural Network
• GBM
• Final Model

In this article, I described a method how we can code CART algorithm in python language. I think it is a good exercise to build your own algorithm to increase your coding skills and tree knowledge. It's not the fastest algorithm implementation but it's enough to understand CART and object oriented programming.

Cart algorithm

In 1983, Breiman et al first described CART algorithm like classified and Regression Tree. It could be used for regression and classification of binary target.

It's a recursive algorithm, at each iteration it finds the best splitting of data which could increase the probability of predicting the target values.

Tree definition

A tree is composed of nodes. Each node could have one, two children or could be a leaf node, a node with no children.

In a list, I stored each object Node.

Tree=[node1, node2, node3, node4….]

Tree structure

To store links and paths between Nodes, I am using Breadth first Search layout.

If the node(x) has a left child the index of the child node is : node(2*x + 1)

If the node(x) has a right child the index of the child node is : node(2*x + 2)

class Node:

def init(self,t,L,R,D,S,V,M) :

self.t=t Index of Node

self.L=L Index of Left child

self.R=R Index of right child

self.D=D Depth of the Node

self.S=S Value of split

self.V=V Feature of split

self.M=M Subset array

self.X=X Execution Flag

Splitting criteria : Gini impurity gain

$$i (n) = 1 - \sum p^2_j $$ where

$$p_j ( t) = \frac{n_t ( t)}{n ( t)}$$

is the relative proportion of category j in Node(t)

For instance, if a data set has only one class target, its gini index is zero, is a purity data set.

We have two subsets : left and right induced by the value A and feature of the split.

$$ \Delta i ( n) = i ( n) - p_{left} i ( n_{left}) - p_{right} i ( n_{right}) $$

For each features, and for each split value A, we compute gini impurity gain and we choose the partition that maximize gini impurity gain.


Example of split calculus

In this example, we considered only one feature x which is continuous.

$$i ( n) = 1 - \left( \frac{10}{20} \right)^2 - \left( \frac{10}{20} \right)^2 = 0.5$$

$$i ( left) = 1 - \left( \frac{10}{13} \right)^2 - \left( \frac{3}{13} \right)^2 = 0.3550295$$

$$i ( right) = 1 - \left( \frac{10}{20} \right)^2 = 0$$

$$\Delta i ( n) = 0.5 - \left( \frac{13}{20} \right)0.355 = 0.2692308$$

If we have had more than one feature, we compute for each feature all the splits, and we choose the best gini impurity gain.

Code of the Cart Algorithm

GitHub link : Tree_Cart_clean.py

Function Build(Tr):

Tr : node of the tree

Ex: T[0] : node one of the list T(list of the all nodes)

Using result of ginisplit function to define or not children Node

Function Split(v,a):

v: feature, a:array

Ex: Split(2,a) : feature 2, array a

Compute maximum gini gain for one feature, it computes all the split possibility

Function Giniplit(matrix):

Matrix: Array of features

Agregate all the best split, feature, and find the best partitionning.

Function TT(depth): Main program

depth : depth of the tree, if you want a max depth of 3 you must type 2.

It's the maximum depth of root node.

Build Tree list and define node

Comparison between my Cart algorithm and sklearn Cart : Same data, same depth, same criterion

import numpy as np

import random


np.random.seed(42)

X=np.random.randint(10, size=(100, 4))

Y=np.random.randint(2, size=100)


a=np.column_stack((Y,X))

Tree from Sklearn

from sklearn import tree


clf = tree.DecisionTreeClassifier(criterion='gini',max_depth=3)

clf = clf.fit(X, Y)

tree.export_graphviz(clf,out_file='tree.dot')

Tree from my Cart algorithm

#Initialisation ot the Tree list 

Tree=[]


#Initialisation of Node[0] 

Tree.append(Node(0,"L","R",0,"S","V",a,"X"))


#2 for 3 depth 

TT(2)


#Print Tree

for index,node in enumerate(Tree):

print index,node.t,"*",node.L,node.R,"*","Depth:",node.D,node.S,len(node.M)

• Cart algorithm
• Gilles Louppe : Understanding Random Forests

To do :

- improve recursion loop in TT function to remove execution flag. This flag prevent algorithm to create children from a Node which has already been built.
- create predict function

- change the storage of array in object to rules to save memory

À partir des données historiques d'un patient (nombre et type de consultation, prise de médicaments, catégories de maladie...) à la période n, nous devons estimer le nombre de jours d'hospitalisation à la période n+1.

Pour cela nous disposons de données années 1 avec le nombre de jours d'hospitalisation années 2, de données années 2 avec le nombre de jours d'hospitalisation années 3 et enfin de données années 3 à partir desquelles on estimera le nombre de jours d'hospitalisation années 4.

J'utilise les données années 1 avec le nombre de jours d'hospitalisation années 2 en échantillon d'apprentissage.
J'utilise les données années 2 avec le nombre de jours d'hospitalisation années 3 en échantillon de test.

Le but est de minimiser le score suivant :

Soit : i : un individu, n : le total d'individu, pred : prévision du nombre de jours hospitalisés, act : le nombre de jours passés à l’hôpital.

Résultat : le meilleur modèle a un score de 0.461197, mon modèle arrive à un score de 0.470774.
En classement final, j'obtiens la 241e place sur un total de 1358 équipes participantes.

1.Action sur les données

• la distribution de DaysinHospital est asymétrique à droite, en réalisant une transformation log je limite cet effet

• pour réaliser un modèle, il faut que l'échantillon d'apprentissage, de test et sur lequel on applique le modèle soit proche or les variables relatives au délai de paiement ont des moyennes et des variances très différentes selon les échantillons, je décide de les exclure
  
  

2.Les meilleurs modèles

À la suite de nombreuses estimations et modèles, je sélectionne les 4 modèles suivants :

3.Ensemble modèle

La combinaison de différents modèles permet de réduire l'erreur c'est pourquoi à partir des quatre meilleurs modèles je construis un ensemble modèle.

Il faut trouver un modèle qui ne suridentifie pas l'échantillon d'apprentissage et ainsi ne pas fait augmenter l'erreur de l'échantillon de test lors de l'application du modèle sur ce dernier.

La régression linéaire donne une erreur de prévision très faible à partir de l'échantillon d'apprentissage (score : 0.28) mais fait augmenter l'erreur de l'échantillon du test (score : 0.51). Il y a donc ici sur-apprentissage.

Un ensemble modèle étant constitué de la combinaison de différents modèles, il me faut trouver la meilleure combinaison.

À travers un algorithme qui calcule l'erreur sur l'échantillon d'apprentissage et sur l'échantillon de test pour chacune des combinaisons, j'ai trouvé la combinaison suivante :

Le score de notre ensemble modèle est de 0.470774.

régression linéaire : bleu
gradient boosting : gris
random forest : rouge,
regularized greedy forest : vert
ensemble modèle : rouge pointillé

Lien de la compétition kaggle

Le fait que Bloomberg, entreprise sérieuse, propose depuis un peu plus d'une année un outil d'analyse des sentiments des tweets à destination des traders souligne le fait que l'utilisation de ces outils va se généraliser dans d'autres domaines. (Détail de la plateforme de Bloomberg)

Dans cet article, je présenterai les deux façons de récupérer les flux provenant de twitter en utilisant la librairie tweepy sous Python.

J'utiliserai pour illustrer les derniers jours de la campagne municipale parisienne opposant Nathalie Kosciusko-Morizet et Anne Hidalgo.

Dans un premier temps, il faut se déclarer développeur sur https://apps.twitter.com ainsi on pourra récupérer les différentes clefs qui nous permettront d'accéder aux flux. Je ne détaillerai pas plus la procédure qui est présente sur de nombreux forums. Un fois ces codes il faut accéder aux flux.

Il y a deux types de flux, ayant chacun des restrictions liées à l'api de twitter, le flux en direct et le flux différé.

1.Flux en direct

Nous récupérons l'ensemble des tweets pendant l'exécution du programme en appliquant des filtres qui peuvent être sur un mot clef, une langue... Pour le flux en direct, il est possible de suivre en théorie jusqu'à 400 thèmes en simultanée, la restriction est qu'il ne faut pas que l'ensemble des tweets que vous récupérez représente plus de 1 % de l'ensemble des tweets. Si ce cas se produit, il y aura une partie des tweets qui seront absents du résultat de votre requête.

Exemple de code pour le flux en direct

Dans le code ci-dessus nous récupérons le 'created_at', date de création et 'text', texte du tweet provenant du status.

Au status d'autres données sont rattachées.

Exemple provenant d'un forum

 'contributors': None,

 'truncated': False,

 'text': 'My Top Followers in 2010: @tkang1 @serin23 @uhrunland @aliassculptor @kor0307 @yunki62. Find yours @ http://mytopfollowersin2010.com',

 'in_reply_to_status_id': None,

 'id': 21041793667694593,

 '_api': <tweepy.api.api object="" at="" 0x6bebc50="">,

 'author': <tweepy.models.user object="" at="" 0x6c16610="">,

 'retweeted': False,

 'coordinates': None,

 'source': 'My Top Followers in 2010',

 'in_reply_to_screen_name': None,

 'id_str': '21041793667694593',

 'retweet_count': 0,

 'in_reply_to_user_id': None,

 'favorited': False,

 'retweeted_status': <tweepy.models.status object="" at="" 0xb2b5190="">,

 'source_url': 'http://mytopfollowersin2010.com',

 'user': <tweepy.models.user object="" at="" 0x6c16610="">,

 'geo': None,

 'in_reply_to_user_id_str': None,

 'created_at': datetime.datetime(2011, 1, 1, 3, 15, 29),

 'in_reply_to_status_id_str': None,

 'place': None

Il est à noter que le flux étant gratuit, il est possible que lors de la capture il y ait des coupures, il faut mettre en place une stratégie pour gérer les erreurs et éviter ainsi de perdre des tweets.

2.Flux différé

Nous récupérons les flux qui ont lieu dans le passé, il est possible de remonter jusqu'à une semaine. La restriction de l'api nous impose à respecter la fréquence de 180 résultats par tranche de 15 minutes, si cette fréquence n'est pas respectée la récupération du flux s'intérrompt du fait que vous avez dépassé les quotas.

Exemple de code pour le flux différé

3.Exemple de la campagne municipale parisienne

Pour illustrer la récupération des flux sous Twitter, j'ai récupéré les flux de la campagne électorale parisienne qui opposaient Anne Hidalgo et Nathalie Kosciusko-Morizet.

L'extraction couvre la période du 26 mars 2014 au 30 mars 2014 inclus. Les données récupérées étaient par minutes et ont été agrégées par heure.

La présence du mot clef "Hildalgo" conditionne la récupération de tweet relatif à Anne Hildago.

La présence du mot clef "NKM" conditionne la récupération de tweet relatif à Nathalie Kosciusko-Morizet.

Les retweets sont inclus.

Nous dénombrons trois pics qui correspondent à :

• le 26 mars à 18h : débat entre les deux candidates

• le 27 mars : dernier meeting des candidates

• le 30 mars au soir : proclamation des résultats

Hors de ces trois pics, on peut constater que les médias sociaux n'ont pas été utilisés à leur pleine mesure, nous ne sommes pas au niveau des Etats-Unis (10.3 millions de tweet pour 1h30 de débat pour la présidentielle américaine) toute proportion gardée.

Pour réaliser ce graphique, la version de ggplot pour Python a été utilisée.

Exemple du code

Critère de classement
5 - Si l'utilisateur a pris une chambre à cet hôtel
1 - L'utilisateur a cliqué pour avoir plus d'information sur l’hôtel
0 - L'utilisateur n'a ni cliqué ni pris une chambre à cet hôtel

Il faut donc pour chaque recherche proposer le meilleur classement soit pour les trois premiers résultats : 5,1,0

Résultat :
Le critère d'évaluation est le NDCG. Il fallait battre le NDCG d'Expedia qui est de 0.49748.
En conclusion, le meilleur modèle a un de NDCG de 0.54075, mon modèle obtient un NDCG de 0.47527.
En classement final, j'obtiens la 132e place sur un total de 337 équipes participantes. La taille de la mémoire nécessaire à l'apprentissage et à la prédiction, la lenteur de divers traitement sous R ont rendu plus ardue la recherche du meilleur modèle.

Pour les modèles de learning ranking (apprendre à classer), il y a trois types approches :

Pointwise
Le problème de ranking est transformé en un problème de régression, classification. La structure de la requête est ignorée.

Pairwise
Ici le problème de ranking est transformé en un problème classification, mais ici les différents résultats de la requête sont classés deux par deux. La structure de la requête dans son ensemble est ainsi ignorée. Cette méthodologie est lourde en calcul du fait du grand nombre de combinaison.(RankSVM)

Listwise
En comparaison avec l'approche pairwise qui minimise les erreurs de classification des documents des requêtes par deux, l'approche pairwise vise à minimiser l'erreur de ranking de la requête. On tient compte ici de la structure de la requête. (ListNet, Lambdarank)

Un excellent document pour commencer :
http://www.hangli-hl.com/uploads/3/4/4/6/34465961/learning_to_rank.pdf

NDCG
Normalized Discounted Cumulative Gain est la mesure de performance d'un système de recommandation. Elle varie de 0.0 à 1.0, où 1.0 représente le classement idéal. Cette mesure est communément utilisée pour évaluer la performance d'un moteur de recherche.

Exemple issu de wikipedia

1	3	0	N/A
2	2	1	2
3	3	1.585	1.892
4	0	2.0	0
5	1	2.322	0.431
6	2	2.584	0.774

SAS a une gestion très poussée de la mémoire, cet outil propriétaire permet de manipuler des bases de données de très grande taille et utilise la mémoire vive ainsi que le disque dur comme espace de travail, or cet outil n'est pas toujours disponible dans les petites structures et présente un coût non négligeable en licence.

R est un outil puissant de modélisation, mais dans le domaine du datamanagement il a un certain nombre de faiblesses : la taille des tables est limitée par la taille de la mémoire vive, les modules complémentaires de R qui facilitent le chargement de grosses tables ne permettent pas d'accéder à l'ensemble des librairies de modélisation.

Pour la partie modélisation, il semble difficile de s'astreindre des contraintes en mémoire vive, on utilisera volontiers des technologies cloud telles qu' Amazon EC2. Pour le datamanagement, il est possible de réaliser la plupart des manipulations de fichier avec des ressources très faibles avec PostgreSQL ou le langage de programmation PYTHON.

1.PostgreSQL

PostgreSQL est un système de gestion de base de données gratuit et très puissant, l'installation est très simple et permet de charger et manipuler des tables de grandes tailles.

Par le passé ayant déjà utilisé des bases de données sur des outils limités tels que Access je constate que la mise en place de PostgreSQL aurait permis à de nombreuses structures de gagner en rapidité et en optimisation des traitements.

Lors d'installation, il faut créer un utilisateur

Lancement de postgres : sudo -i -u postgres puis psql.

Commande base

Lister les tables : \d

Lister les colonnes d'une table : \d nom_de_la_table

Effacer table : drop table nom_de_la_table;

Création de la table et de ces variables

CREATE TABLE TRAIN (ID BIGINT PRIMARY KEY,TITLE TEXT,BODY TEXT,TAGS TEXT);

Importer un fichier csv dans la table TRAIN

COPY TRAIN FROM '/home/alfard/Documents/Kaggle/Facebook2/Train.csv' csv header;

Une fois la table importée on peut réaliser les traitements en utilisant le langage SQL.

Exporter une table en fichier csv

COPY matable TO '/root/pg/extraction.csv' WITH CSV;

2.Python

Python est un langage de plus en plus en vogue dans le monde de la data, il est simple à apprendre, dispose de plus en plus de librairies spécialisées et a reçu récemment un soutient financier de près de 3 millions de dollars de la DARPA.

Il est possible avec Python de charger un fichier ligne et par ligne, d'y appliquer des traitements et de l'enregistrer ligne par ligne, la consommation en mémoire vive est alors très faible.

#A noter sur python, on compte à partir de 0 et non de 1

import pandas as pd
import numpy as np
import csv


f = open('/home/alfard/Documents/Kaggle/Facebook2/Train.csv',"rb")
fw= open('/home/alfard/Documents/Kaggle/Facebook2/TrainClean.csv',"w")

 
fileopen = csv.reader(f,delimiter=',', quotechar='"')
filewrite = csv.writer(fw,delimiter=',', quotechar='"')

#Mise en place d'un compteur pour connaître le nombre de ligne traitées 
p=0
for row in fileopen:
    #Exemple de traitement de caractère
    for i in range(len(stopwords)):
        
        row[1] = row[1].lower()
        row[1] = row[1].replace(stopwords[i],' ')
        row[1] = " ".join(row[1].split())
        

        #Ecriture du résultat du traitement dans le fichier TrainClean.csv       
    filewrite.writerow([row[1],row[3]])
    p=p+1
    print p
        # row[1],row[3] signifie que je copie la 2e et le 4e élement de la ligne traitée,
        # à la fin de traitement, l'ensemble des éléments de la ligne constituera une colonne 
        # CAD la 2e colonne et la 4e colonne 
        

#Obligatoire pour réaliser l'écriture du csv
f.close()
fw.close()

À la suite du concours Heritage Health Prize j'ai constaté que sous R le moteur de calcul matriciel n'était pas rapide (je ne fais pas référence aux calculs des librairies spécialisées) ce qui n'est pas vraiment compatible avec une utilisation professionnelle.
Il existe pourtant des solutions qui permettent de faire du calcul matriciel avec un rapport coût/performance intéressant. J'en ai testé quelques-unes.

Le modèle ensemble est une combinaison de modèles. Il faut déterminer les meilleures combinaisons c'est à dire celles qui permettent de minimiser l'erreur de l'échantillon d'apprentissage et de test.

Pour chaque modèle, nous avons une prévision, le poids de ces prévisions dans l'ensemble modèle est recherché de la façon suivante :

Boucle
Prévision de l'ensemble modèle
(
prévision modèle 1poids de 1+
prévision modèle 2poids de 2+
prévision modèle 3poids de 3+
prévision modèle 4poids de 4
)
/100

Poids de 1+Poids de 2+Poids de 3+Poids de 4 = 100

Le poids de 1,2,3 et 4 varie de 0 à 100, soit près de 156 849 combinaisons.

Pour chacune de ces combinaisons, je calcule l'erreur de l'échantillon d'apprentissage et de test.

Julia : http://julialang.org/
Julia est un nouveau langage qui a été crée par le MIT qui a pour objectif de permettre le calcul matriciel avec un code proche de Matlab et de performance proche du C/Fortran. Il est relativement simple de transcrire du code Matlab en Julia. Inutilisable professionnellement à l'heure actuelle car trop jeune cependant il a un énorme potentiel.

Armadillo : http://arma.sourceforge.net/
Armadillo est une librairie qui permet de faire du calcul matriciel en C++ avec un code proche de Matlab. Elle présente d'excellente performance.

Ne disposant pas de la librairie IML je n'ai pas pu réaliser de test sous SAS

Conclusion : Il existe toujours des solutions pour accélérer un traitement.

Les données doivent devenir le pétrole du 21ie siècles. À travers de nombreux cours et de concours de modélisation une communauté de data scientist est en train de se constituer. Le big data est un terme à la mode alors que déjà les statisticiens font du big data depuis un certain temps.

Pourquoi une telle effervescence ?

• des outils libres ont fait leurs apparitions et ainsi ont permis de diffuser le data mining, alors qu'auparavant seuls des outils onéreux étaient disponibles.
• il est possible de louer de la puissance de calcul en faisant appel à des services de cloud comme Azure ou Amazon EC2 pour quelques euros

La contrepartie à ce développement est de mettre de coté l'importance des statistiques pour la réalisation et la compréhension des modèles.

Je ferais une analogie avec le monde du jeu d'échec.

En 1988, Gary Kasparov a perdu une partie face à Deep Blue qui était capable de calculer 300 millions de coups par seconde. Pour caricaturer le statisticien comme le jour d'échec connait les raccourcis pour trouver le meilleur modèle en limitant les calculs, tandis que l'informaticien alias Deep Blue connait la puissance brute pour trouver le meilleur modèle hors pour devenir un bon data scientist il faut combiner ces deux aspects, il faut être à la fois un bon statisticien et un bon informaticien.

Or aujourd'hui, il n'y a de moins en moins une analyse en amont des données, du travail d'enquêtes pour comprendre ces dernières. J'ai vu de nombreux concours de modélisations où la connaissance des données n'était pas un préalable obligatoire :

• la KDD Cup 2009 sponsorisée par Orange pour expliquer le comportement du consommateur avait la majorité des variables non expliquées
• le concours Kaggle pour prévoir la durée d'hospitalisation avec un ensemble de prix s'élevant à 3 millions de dollars a un certain nombre de variables dont l'articulation est obscure pour l'ensemble des participants

Cette méconnaissance des données risque de transformer le travail de modélisation en un simple concours d'algorithme et de puissance de calcul. Il faut utiliser la puissance de calcul à bon escient, et ne pas négliger le travail d'études statistiques des données en amont.

En 2010 des courriers ont été envoyés dans les foyers américains pour réaliser le recensement. Afin d'optimiser sa communication le bureau du recensement recherche le meilleur modèle qui sera à même de prévoir le taux de retour des plis, pour cela nous avons à notre disposition le taux de retour réel de 2010, ainsi qu'un ensemble de données sociodémographique. Les données sont au niveau des "block group" qui regroupent jusqu'à 25 000 habitants. Il a été possible d'inclure des données externes après validation par le bureau du recensement américain.

L'évaluation se fait par la moyenne absolue des erreurs pondérées par le poids des populations (weighted mean absolute error).

Résultat : le meilleur modèle a un WMAE de 2.54476, mon modèle arrive à un WMAE de 3.27624.
En classement final j'obtiens la 98e place sur un total de 244 équipes participantes.

2: Avec un nombre si élevé de variables il est difficile de choisir les données qui appartiendront au modèle, la première approche qui j'avais mis en place fut de supprimer les valeurs qui étaient trop corrélées entre elles, cela fut fastidieux sans des résultats probants.
Dans mon esprit un modèle ne doit pas avoir trop de variables explicatives, car le but d'un modèle c'est de pouvoir déterminer les variables les plus déterminantes.
J'ai pour la première fois utilisés la sélection de variables en utilisant un random forest, cette technique semble utilisée dans de nombreuses situations, des utilisateurs de SAS Miner ne disposant pas de random forest utilisent R dans un premier temps pour essayer de déterminer les variables les plus explicatives.
Par un processus itératif, j'isole 40 variables.

3: Pour tenter d'améliorer le modèle, j'ai inclus les données géographiques (longitude, latitude) ainsi que les taux de chômage de chacune des zones, je n'ai pas d'amélioration significative du modèle.

4: Sur mes données j'applique les trois modèles suivants : random forest, gradient boosting ainsi que neural network.
Le gradient boosting donne les meilleures prévisions suivies des random forest et du neural network. Ne disposant pas de puissance de calcul nécessaire, je n'ai pas pu réaliser le meilleur des tunings. L'utilisation du cloud amazon EC2 devrait à l'avenir me permettre d'éviter cet écueil.

5: En combinant mes trois modèles j'ai significativement amélioré mes estimations, j'ai réalisé un ensemble modèle linéaire, une simple régression du taux de retour réel versus les taux de retour estimé par les trois modèles. L'utilisation la plus connue du modèle ensemble fut lors du concours Netfix.
Papier de référence au sujet de la combinaison de prédictions

Lien de la compétition kaggle

Pour cela nous disposons d'un historique de deux années comportant 67 millions de clics et 27 millions de recherches, ainsi d'un catalogue des jeux Xbox avec leur référence produit respective.

Résultat : la meilleure méthodologie arrive à prévoir en moyenne 78 % de clics sur un produit, ma méthodologie arrive à en prévoir en moyenne 64 %. En classement final, j'obtiens la 65e place sur un total de 97 équipes participantes.

Pour chaque requête, il y a des probabilités différentes de choisir différents produits, je conserve la probabilité de choix la plus élevée que je mets dans un dictionnaire que j'appelle dictionnaire historique.
Le cas des requêtes existantes a été ainsi traité.

Si une requête n'existe pas dans l'historique, j'utilise alors un dictionnaire que j'appelle dictionnaire nouveau. Ce dictionnaire est constitué de la manière suivante :
je procède à un rapprochement textuel c'est-à-dire si je réalise la requête suivante :"tetris", alors le moteur va rechercher dans le catalogue des produits un nom de produits qui se rapproche le plus du nom "tetris" en utilisant la fonction SAS compged (voir ci-dessous). Je n'ai pas réalisé volontairement de filtre sur le catalogue des produits pour n'intégrer que les produits xbox ce qui fait réduire les résultats de ma prévision.

Normalisation : suppression des accents, des majuscules, des espaces

La fonction compged

Elle permet de calculer sous la forme d'un score la distance généralisée entre deux chaines de caractères. Plus le score GED est faible plus il y a correspondance des deux chaînes de caractères.

Cette fonction est très utile, je l'ai déjà utilisée dans un programme pour détecter des emails crées par des robots, les emails ayant un corps identique auquel une série de chiffres (exemple : toto1@yahoo.fr, toto2@yahoo.fr, toto3@yahoo.fr).

Lien de la compétition kaggle

L’intérêt des compétitions de Kaggle est qu'une fois l'épreuve terminée il est possible de découvrir les meilleures solutions. La méthodologie gagnante sera bientôt accessible dans le détail, car le bureau du recensement a décidé de partager le modèle avec l'ensemble de la communauté, en attendant il ressort plusieurs pistes qui ont permis aux compétiteurs de s'élever en haut du classement

• à partir des variables de départ de nouvelles variables ont été créées à la suite d'étude des relations des variables entres-elles souvent en réalisant de simples régressions, à noter qu'un module de R : Mars permet de rechercher ces interactions (à voir)

• Les modèles n'ont pas été simplement lancés en testant différents paramètres, mais un gros travail de "fitting" a été réalisé c'est-à-dire que pour chacun des paramètres des modèles des intervalles ont été définis et le programme de recherche du modèle a tourné le temps que les modifications de ces paramètres permettent de réduire l'erreur du modèle. Cette méthode a donné de bons résultats, mais a été gourmande en calcul : pour trouver le nombre d'arbres optimal en utilisant le gradient boosting cela a nécessité sur un ordinateur personnel près de trois journées de calcul. Il est donc nécessaire de savoir programmer pour trouver les meilleurs paramètres.

• il semble obligatoire de réaliser un ensemble model c'est-à-dire en utilisant différents modèles pour réaliser la prévision