Editerna

Aller au contenu | Aller au menu | Aller à la recherche

mardi, août 27 2013

Déménagement

Ce blog a changé d'adresse. Veuillez vous diriger vers la nouvelle version : http://editerna.free.fr/wp/

lundi, août 26 2013

Envoyer des e-mails avec R (sous Mac OS X)

Décidément, on peut vraiment faire pas mal de choses avec le logiciel R. La dernière fonctionnalité qui m'a laissé admiratif est la possibilité d'envoyer des e-mails en seulement quelques lignes de code. Je pense que ça peut s'avérer très appréciable à l'utilisation dans le cadre d'une mise en place de reporting. D'autres (je ne citerai personne !) y verront un avantage dans la possibilité de transmettre des documents personnalisés à une liste de contacts (par exemple, des sujets de TPs ou d'examens...).

Mon utilisation est toute autre. En ce moment, je dois exécuter une programme quotidiennement, et j'ai donc mis les mains dans le cambouis pour découvrir la planification de tâches sur Mac OS (avec le service cron). Le souci, c'est que ce genre de tâche peut devenir très envahissant dans une journée, s'il faut vérifier s'il n'y a eu aucun problème (comme il s'agit de récupération de données, l'environnement extérieur peut être soumis à une certaine variabilité). Aussi, j'aimerais recevoir un message contenant des sorties R pour m'assurer en un clin d'oeil du bon fonctionnement de l'opération.

On trouve sur la toile pas mal de pages indiquant comment procéder pour envoyer un e-mail avec R, avec ou sans pièce jointe. La solution qui m'a semblé la plus simple est d'utiliser le package sendmailR. Je me suis appuyé sur le contenu des réponses à une question sur stackoverflow.

Alors pourquoi, alors qu'il existe déjà pas mal de ressources fournissant des explications sur la procédure à mettre en oeuvre pour envoyer des e-mails avec R, rédiger une petite note ici ? Eh bien tout simplement parce que j'ai passé un peu de temps avant de réussir... et que je souhaite réunir les différentes informations nécessaires dans un même post.

La première chose à faire, est de configurer la machine pour envoyer des e-mails à partir de Localhost. Pour ce faire, il suffit de suivre à la lettre ce tutoriel : http://www.zenddeveloper.com/how-to-send-smtp-mails-with-postfix-mac-os-x-10-8/.

Il y a cependant un détail qui m'ennuie dans cette procédure (peut-être m'y prends-je de la mauvaise façon ?) : le mot de passe du compte Gmail utilisé est écrit en clair dans un fichier... Ainsi, avant l'étape 3 du tutoriel, il est possible d'ajouter quelques consignes : (i) supprimer le fichier sasl_passwd pour ne conserver que le fichier sasl_passwd.db. Ensuite, une petite modification des droits de lecture/écriture. Dans le terminal, il me semble que ça se fait avec :

$ chmod u+rw-x,g-rwx,o-rwx sasl_passwd.db

Avec ça, il est possible d'envoyer ses e-mails depuis le terminal, en passant par les serveurs d'envoi de Gmail. Plutôt pratique.

Enfin, il reste à écrire quelques lignes sous R, pour pouvoir effectuer la tâche depuis une session R.

library(sendmailR)
from <- sprintf("", Sys.info()[4])
to <- ""
subject <- paste("Rapport",format(Sys.time(),"%Y-%m-%d"),sep=" ")
body <- paste("Voici le rapport du ",format(Sys.time(),"%Y-%m-%d"),sep="")
attachmentPath <- "chemin/vers/fichier"
attachmentName <- "nomFichier.txt"
attachmentObject <- mime_part(x=attachmentPath,name=attachmentName)
bodyWithAttachment <- list(body,attachmentObject)
sendmail(from=from,to=to,subject=subject,msg=bodyWithAttachment)

Note : impossible d'afficher les premières lignes de code correctement... Les voici hors des balises :

from <- sprintf("<sendmailR@%s>", Sys.info()[4])

to <- "<prenom.nom@gmail.com>"

vendredi, juin 28 2013

Deuxièmes Rencontres R à Lyon

Le 27 et 28 juin 2013, se déroulaient les 2e Rencontres R, à Lyon. J'aimerais faire un petit bilan de ce que j'ai pu voir pendant ces deux journées.

En premier lieu, j'aimerais souligner la qualité de l'organisation. Tout était millimétré, l'accueil était très agréable et les membres de l'équipe organisatrice étaient très sympathiques. Dès l'arrivée, une malette (biodégradable) est remise à chaque participant, contenant un programme relié et en couleur, ainsi que quelques goodies provenant de Revolution Analytics.

Les lignes qui suivent sont plutôt un pense-bête personnel, mais je vais tenter de les rendre compréhensibles. Je vais me limiter ici à quelques présentations, surtout à celles qui me parlent le plus.

Session de tutoriels (26 juin 2013)

Avant les rencontres, j'ai assisté à deux tutoriels. Le premier était donné conjointement par Anne-Béatrice Dufour et Sylvain Mousset. Il était question de présenter l'intégration de code R dans des documents LaTeX, avec Sweave. Certains enseignants de l'Université Claude Bernard Lyon 1 ont réalisé un projet plutôt génial avec R+LaTeX, mais on reviendra sur ce point un peu plus tard.

Le second tutoriel était animé par Christophe Genolini, et avait pour but de familiariser les participants à la création de package avec R. Parmi les nombreuses ressources du site de Christophe, on retoure un document PDF pour créer ces fameux packages. S'il faut retenir une seule chose de cette session, il s'agit de l'accent qui a été mis sur la nécessité de créer des documentations complètes, qui favoriseront l'adoption des packages par les utilisateurs. Cette idée a été répétée à de nombreuses reprises durant les rencontres.

Première journée (27 juin 2013)

Après un petit mot d'ouverture par Stéphane Dray, dans une ambiance détendue, le premier orateur monte sur la scène. Il s'agit tout simplement de Hadley Wickham. Oui, à Lyon, ils ne font pas les choses à moitié !

Visualising big data in R

L'auteur du package ggplot2 ("djidjiplot tou" pour les intimes (@phnk)) vient présenter un nouveau package : bigvis. Ce dernier permet de travailler avec des jeux de données conséquents. Pour illustrer la présentation, Hadley s'appuie sur plus de 60 millions d'observations sur les distances et les durées de vols d'avions aux États-Unis. Les données sont stockées dans un fichier RDS. Pour les importer, il suffit de s'appuyer sur la fonction readRDS. L'idée, pour produire un graphique de la distance en fonction de la vitesse, est de procéder en trois étapes :

  1. condenser les données ;
  2. effectuer un lissage ;
  3. afficher les résultats sur un graphique.

Pour rendre la présentation un peu plus impressionnante, Hadley Wickham fait appel au package Shiny. Je suis scotché ! En un rien de temps (moins de 3 secondes), suite à un choix de fenêtre effectué en faisant glisser un curseur, la version condensée (fraîchement recalculée) des +60 millions de points est graphée sur la page.

Par la suite, Hadley nous touches deux mots sur Rcpp, qui permet d'écrire du code en C++ directement dans R. Plutôt génial (voir plus bas pour la présentation de Romain François).

L'analyse de données avec FactoMineR : les nouveautés

François Husson vient présenter les améliorations apportées au package FactoMineR. Premièrement, une gestion des données manquantes, à l'aide du package missMDA est proposée. D'un point de vue graphique, un argument fait son apparition, et il est plutôt sympathique : autolab. Finis les graphiques de projection illisibles, les légendes sont placées de manière à ce que leur recouvrement soit le moins grand possible.

rCarto, un package de cartographie statistique

Timothée Giraud nous montre le travail qu'il a réalisé pour produire de bien jolies cartes géographiques avec R, en travaillant notamment sur les légendes, ou de manière plus générale, sur leur habillage. Son package se nomme rCarto, et est disponible sur le CRAN.

What a statistician might want to know about human color vision, but was afraid to ask!

Le thème proposé par Kenneth Knoblauch laisse songeur. Il est question des couleurs utilisées pour produire des illustrations. Lorsque je produis un graphique, je ne pense jamais à ceux qui ont des troubles dans la vision des couleurs... Or, ils sont nombreux (8% selon Kenneth). Le package dichromat permet d'afficher les couleurs telles qu'elles sont perçues par les daltoniens. À méditer pour les prochaines sorties !

Kml3D: K-means pour données longitudinales jointes

Christophe Genolini vient nous parler des données longitudinales, et du package Kml3D. Il est vrai que c'est encore quelque chose que je ne connais quasiment pas. Je ne vais donc pas m'étendre sur le sujet. J'aimerais juste conserver une trace de ce que nous a proposé Christophe, qui bien qu'anecdotique, pourrait quand même s'avérer être amusant : la possibilité d'inclure un graphique 3D dans un document PDF, avec la fonction plot3dPdf.

R and the Cloud

Le deuxième conférencier invité, Karim Chine a fait un speech qui m'a particulièrement impressionné, mais qui me laisse perplexe dans le même temps. Il s'agit de présenter la relation possible entre R et des machines sur le cloud, notamment sur Amazon. En quelques secondes, Karim commande une machine sur le cloud, pour quelques centimes de dollars. C'est impressionnant, mais il n'y a rien d'extraordinaire (enfin... je pense que pour mes grand-parents, c'est une toute autre histoire !). Les serveurs d'Amazon sont apparemment aux États-Unis, la question de consommation d'énergie est occultée ici.

Là où j'ai été bluffé, c'est quand Karim nous a introduit sa plateforme Elastic-R, qui permet non seulement de commander des machines sur le cloud rapidement, mais surtout de faire de l'édition en direct, à plusieurs, de codes R. Un utilisateur (Romain, dont le nom m'a échappé) de l'amphi est invité à rejoindre la session créée par Karim, et à effectuer quelques modifications dans la console R. Mieux encore, Romain est invité à changer des valeurs dans une feuille de calcul Excel, qui interagit avec R ! Les modifications sont répercutées quasi instantanément sur les données, dans R ! Mais ça ne s'arrête pas là ! Karim reprend complètement la main, et ouvre un document Word, dans lequel il inscrit quelques lignes de code R, se terminant par un appel à la fonction plot. Comme par enchantement, le graphique correspondant s'affiche dans le document Word ! Cette librairie est prometteuse.

R2GUESS: a GPU-based R package for sparse Bayesian variable selection

Habib Saadi nous a présenté le package R R2GUESS qui ne devrait pas tarder à arriver sur le CRAN. Il permet de faire de la sélection de variables, en utilisant la puissance de calcul du GPU.

Intégration R et C++ avec Rcpp

Encore une partie de la conférence qui m'a pas mal plu. Romain François nous a présenté le package Rcpp, qui comme son nom le laisse deviner, permet d'écrire du code C++ directement dans R. L'objectif affiché ? Améliorer la rapidité des calculs. Les projets de Romain ont l'air tous mieux les uns que les autres...

R as a sound system

Pour finir la journée en beauté, avec une touche de poésie, nous avons eu le plaisir d'écouter Jérôme Sueur (conférencier invité), au sujet de la gestion des sons avec R, et plus particulièrement à l'aide des packages tuneR (ou audio) et seewave. Jérôme utilise R pour faire, entre autre, de la bioacoustique. R est utilisé pour analyser des sons. Une utilisation intéressante proposée par Jérôme est d'écouter des données. Oui, oui ! Écouter des données ! Dans un premier temps, une petite citation d'un homonyme de l'orateur, Jérôme Bonaldi s'affiche sur une slide :  C'est totalement inutile et donc rigoureusement indispensable ! (c'est plutôt bien choisi, pour une conférence à Lyon). Mais en réalité, écouter les données n'est pas aussi inutile que ce que l'on pourrait croire. Dès que je retrouve une connexion internet plus puissante, j'ajouterai un exemple ici. Grosso modo, l'écoute des données peut permettre d'identifier des points aberrants. Comme le fait remarquer Joel Gombin, La datasonification permet aux non voyants de "visualiser" des données.

Seconde journée (28 juin 2013)

Après une soirée autour d'un buffet dinatoire et accompagné de gens fort sympathiques, puis une bonne nuit de sommeil dans la résidence de l'INSA, il est temps d'attaquer la seconde journée des 2e rencontres R.

L'approche par comparaison de modèles avec R2STATS dans l'enseignement des statistiques en sciences humaines

Le premier conférencier invité de la journée est Yvonnick Noël. Il nous parle de l'interface R2STATS, qu'il a développé pour permettre à ses étudiants en psychologie de comparer des modèles statistiques sans passer par la rédaction de codes R. R2STATS permet d'estimer des GLM ou des GLMM, avec choix de distribution. La fenêtre de l'interface est organisée en six onglets : "fichiers", "données", "résultats", "graphiques" et "comparaisons". C'est assez intuitif, mais je trouve qu'on perd trop de liberté. Il est dommage de ne pas pouvoir récupérer le code R qui a servi à produire les sorties, pour le modifier (un peu à la manière de SAS qui permet de rappeler le code produit lors d'une procédure réalisée à l'aide du Wizzard).

Un site web d'enseignement R automatisé et à gestion partagée

Trois enseignants chercheurs (Anne-Béatrice Dufour, Daniel Chessel et Jean R. Lobry) ont mis en place un site web génial d'enseignements de statistique en biologie. Sur ce site, les enseignants peuvent déposer des cours, des fiches de TD, des sujets d'examens, etc. Simon Penel nous explique le fonctionnement de ce site.

Ce que je trouve génial dans ce qu'ils ont réalisé, c'est que les TD et les examens sont produits avec LaTeX et R via Sweave. Les sujets créés peut être différents chaque année, via la possibilité d'introduire des valeurs aléatoires grâce à R. L'enseignant réalise son sujet avec Sweave, et peut ainsi produire de manière automatique une feuille de composition ainsi qu'une de correction. Les années suivantes, la méthode des exercices ne change pas, mais les résultats si. Cela permet par exemple aux redoublants de ne pas avoir exactement le même sujet !

Génération automatique de documents pédagogiques avec R pour l'enseignement et l'évaluation des étudiants

Autant vous prévenir tout de suite, les vingt minutes de Sylvain Mousset m'ont fait rêver. Il est question de génération automatique de documents avec R. On reste donc dans le même esprit que la présentation précédente. Sylvain fait rêver une bonne partie de la salle (je pense) en nous montrant comment il réalise certaines évaluations, à l'aide de QCM. Pour ce faire, il utilise un logiciel appelé AMC (Auto Multiple Choice), développé par Alexis Bienvenüe. Le logiciel se charge de créer un questionnaire au format LaTeX, avec ordre des questions variables pour chaque étudiant, et ordre des propositions de réponses défini aléatoirement. Fini la triche lors des QCM ? François Briatte se pose des questions intéressantes (j'aimerais bien connaître les réponses) : est-ce que la mise en place de ce type de QCM a eu un impact sur les notes ? Sur la triche ? Quoi qu'il en soit, ce type de questionnaires est alléchant, d'autant plus qu'AMC ne s'arrête pas là, et propose de corriger les copies ! Formidable, non ? Plusieurs centaines de copies corrigées en trois heures...

Pistes de réflexion pour la mise en oeuvre d'un enseignement à distance sur le test du khi carré d'indépendance pour des étudiants en master de sciences de l'éducation

Mehdi Khaneboubi introduit les difficultés rencontrées pour l'enseignement de l'analyse de données avec R, à distance (avec des étudiants résidant dans plusieurs pays différents). L'avantage que propose R face à d'autres logiciels de statistique est sa présence sur de nombreuses plateformes, et même en ligne dans le cas où les étudiants se trouvent dans un cyber-café.

De la biologie à l'algèbre linéaire... en passant par R : expérimenter la notion de projection

Anne-Béatrice Dufour présente une manière d'introduire le sujet de la projection (dans le cadre de la régression linéaire) aux étudiants. J'ai trouvé son approche très intelligente, et assez intuitive. Le sujet d'étude doit interpeler les biologistes : étudier la relation entre la taille et l'empan (distance entre le bout du pouce et l'extrémité de l'auriculaire) de chaque main (data(survey) du package MASS). L'étudiant est amené à représenter graphiquement les mesures, pour à chercher de manière heuristique une droite qui permet de s'ajuster au mieux aux données, de manière à minimiser les écarts entre la droite et les observations.

Lightning talks 1

L'exercice est difficile pour les personnes devant présenter leur travail en 6 minutes. Toutes les 24 secondes, une diapositive défile automatiquement. Je me limite à deux intervenants ici.

Aspects de la cartographie thématique pour les sciences sociales avec R

Joël Gombin ouvre le bal avec la présentation de son package pour visualiser des données géographiques. L'idée de Joël est de proposer un aperçu des méthodes pour construire des cartes avec R, dans un domaine relatif aux sciences sociales. Les fonctions sont écrites en POO et disponibles sur son espace GitHub.

Prévision de consommation électrique avec R

Raphaël Nedellec expose des méthodes utilisées au sein du service R&D d'EDF pour faire des prévisions de consommation électrique. J'étais content de voir qu'un grand groupe comme EDF s'ouvre sur le logiciel R.

TraMineR : Une boîte à outils pour l'exportation et la visualisation de séquences

Gilbert Ritschard, dernier conférencier invité, nous montre le travail colossal qui a été réalisé autour de l'analyse de séquences. Après les prestations éclair lors des lightning talks, le contraste est très fort, puisque Gilbert est extrêmement posé, prend le temps d'expliquer chaque notion, et s'explique très clairement (je ne sous-entend aucunement que les intervenants des lightning talks ne s'exprimaient pas clairement, mais la rapidité avec laquelle s'enchaînent les diapositives entraîne un rythme effréné, d'où la présence du contraste). Le package TraMineR permet des représentations graphiques de séquences chronologiques. La conférence est illustrée avec des données de séquences de vies (études, chômage, emploi, etc.).

multcmm : fonction d'estimation de modèles mixtes à processus latent pour données longitudinales multivariées

Viviane Philipps présente la fonction multlcmm du package lcmm, qui permet d'estimer des modèles mixtes à processus latent, dans le cas où les données sont longitudinales multivariées. Je ne suis pas du tout familier avec ce genre de modèles, j'avoue ne pas avoir très bien suivi cette présentation.

Données tronquées sous R dans les modèles linéaires simples et à effets mixtes

Djeneba Thiam propose de se pencher sur des données censurées, avec effets aléatoires.

Prédiction d'un événement binaire à partir de données fonctionnelles : application aux bovins laitiers

Dans cette présentation, Cécile Sauder nous parle de ses travaux relatifs à la gestation des vaches. Il s'agit de prévoir l'état d'une vache (gestante ou non), 42 jours après son insémination, en fonction de ses courbes de lactation.

Génération de documents Word à partir de R : utilisation du package R2DOCS dans uneplateforme statistique en milieu industriel

David Gohel commence par fixer les bases : il veut lier R et Word, en réponse aux attentes de ses clients, qui  souhaitent travailler avec Word uniquement (entendre ici : pas avec Libre Office). Son package s'appelle R2DOCX, et est disponible sur GitHub. L'idée est de pouvoir faire du reporting avec R et Word, en permettant un respect de la charte graphique des clients. À l'avenir, David souhaite étendre son travail à la passerelle entre R et PowerPoint.

R-dyndoc une alternative à Sweave

Rémy Drouilhet procède à l'ultime présentation de ces 2e rencontres R. Il nous montre son projet, R-dyndoc, qui se veut être une alternative à Sweave. Rémy s'appuie sur le langage de programmation Ruby pour développer son outil. L'avantage de Dyndoc, serait la possibilité d'utiliser plus de librairies que Sweave. Pour un utilisateur lambda comme moi, je trouve que Sweave est plus abordable, dans la mesure où une fois qu'on sait se débrouiller un minimum avec LaTeX et R, il est facile de créer des documents avec Sweave. En revanche, avec Dyndoc, il faut encore changer dans la manière de raisonner pour construire les documents, et je ne suis pas prêt à investir du temps dans l'apprentissage de ce code.

Une autre alternative à Sweave, qui a été à plusieurs reprises mentionnée lors des rencontres, serait knitr, plus récent que Sweave. À creuser !

Voilà, le petit aperçu des 2e rencontres R s'arrête ici. Comme il y avait deux sessions en parallèle, il manque un bon nombre de présentations dans ce billet.

mercredi, juin 12 2013

Visualizing densities of spatial processes

In the past few weeks, Arthur Chaprentier (@freakonometrics) and I have been working again on an article untitled "Visualizing spatial processes using Ripley's correction: an application to bodily-injury car accident location".  The revised version of our work is available here: http://hal.archives-ouvertes.fr/hal-00725090.

Please note that this post is related to Arthur's: http://freakonometrics.hypotheses.org/7129

We provide an easy way to correct bias of density estimation of edges of regions, by investigating and extending Ripley's circumference method. An application to bodily-injury car accidents that happened in Finistère and Morbihan (two French regions) in 2008 is discussed in the paper, and some R codes are provided too. The purpose of this post is to explain the R code which enables to plot the results of the estimation. For those who are interested in the way to compute the estimation of the density, feel free to read the article on Hal. Please note that Arthur Charpentier wrote a nice article last year, explaining the method: http://freakonometrics.hypotheses.org/2416, but be aware that some tiny changes have been made to the functions.

There were 186 geolocated bodily-injury car accidents in 2008 in Finistère, and 180 in Morbihan. Their spatial distribution can be viewed on the maps below, with a standard kernel on the left and corrected one on the right.

Now, let's have a look at the R code, shall we? It has been tested on R 3.0.1 and R 2.15.3 too. First, you'll need to load the following packages: ggplot2 and ggmap.

Let us assume that the result of the density estimation is stored in the object called result

Let us send a query to Google Maps API to get a nice map.

theMap <- get_map(location = c(left = min(pol[, 1]), bottom = min(pol[,
2]), right = max(pol[, 1]), top = max(pol[, 2])), source = "google",
messaging = F, color = "bw")

Since we are using a function that relies on ggplot2 grammar, we need to put the data in the right format. To do so, we've created this little function :

getMelt <- function(smoothed) {
res <- melt(smoothed$ZNA)
res[, 1] <- smoothed$X[res[, 1]]
res[, 2] <- smoothed$Y[res[, 2]]
names(res) <- list("X", "Y", "ZNA")
return(res)
}

Let us apply this function to the list containing the results of the estimation :

smCont <- getMelt(result)

Now, we should define a color scale and some breaks.

breaks <- seq(min(result$ZNA, na.rm = TRUE) * 0.95, max(result$ZNA, na.rm = TRUE) * 1.05, length = 21)
col <- rev(heat.colors(20))

We should also change the labels we want to display in the legend :

theLabels <- round(breaks, 2)
indLabels <- floor(seq(1, length(theLabels), length.out = 5))
indLabels[length(indLabels)] <- length(theLabels)
theLabels <- as.character(theLabels[indLabels])
theLabels[theLabels == "0"] <- " 0.00 "

Finally, the map can be built :

P <- ggmap(theMap)
P <- P + geom_point(aes(x = X, y = Y, col = ZNA), alpha = 0.3,
data = smCont[!is.na(smCont$ZNA), ], na.rm = TRUE)
If one wants to add a contour :
P <- P + geom_contour(data = smCont[!is.na(smCont$ZNA), ], aes(x = X,
y = Y, z = ZNA), alpha = 0.5, colour = " white ")

If you display the map now, there is something missing: colors! So, the plot has to be updated.

P <- P + scale_colour_gradient(name = "", low = " yellow ", high = "red",
breaks = breaks[indLabels], limits = range(breaks), labels = theLabels)

We might want to remove some elements, such as the axis ticks and legends:

P <- P + theme(panel.grid.minor = element_line(colour = NA), panel.grid.minor = element_line(colour = NA),
panel.background = element_rect(fill = NA, colour = NA), axis.text.x = element_blank(),
axis.text.y = element_blank(), axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(), axis.title = element_blank(),
rect = element_blank())

And, to finish, let us draw the border of the polygon (assuming the polygon is named pol):

polDF <- data.frame(pol)
colnames(polDF) <- list("lon", "lat")
(P <- P + geom_polygon(data = polDF, mapping = (aes(x = lon, y = lat)),
colour = " black ", fill = NA))

To account for traffic, we've downloaded road sections from http://www.geofabrik.de/index.html, which provides files based on OpenStreetMap data. Each observation is a section of a road, and contains a few points identified by their geographical coordinates that allow to draw lines. We interpolate points between the extremities of road sections, weighting the distance separating each points according to the importance of the road. Hence, these points are used as a proxy for traffic.

splitroad <- function(listroad, h = 0.0025) {
pts = NULL
weights <- types.weights[match(unique(listroad$type), types.weights$type), "weight"]
for (i in 1:(length(listroad) - 1)) {
d = diag(as.matrix(dist(listroad[[i]]))[, 2:nrow(listroad[[i]])])
for (j in 1:(nrow(listroad[[i]]) - 1)) {
pts = rbind(pts, cbind(seq((listroad[[i]])[j, 1], (listroad[[i]])[j + 1, 1], length = weights * d[j]/h),
seq((listroad[[i]])[j, 2], (listroad[[i]])[j + 1, 2], length = weights * d[j]/h)))
}
}
return(pts)
}

The result can be seen on the maps below.

That's pretty much it. Oh, by the way, I will be presenting a poster on which Arthur and I have been working during the 2nd French R Meeting in Lyon, on the 27th-28th of June 2013. f you want to have a look, click on the image down below to view a PDF version (the document is written in French).

mardi, avril 9 2013

Charles BaudelR

Aujourd'hui est le 192e anniversaire de Charles Baudelaire, comme me l'a gentiment rappelé le Doodle de Google de ce jour.
Or, ce week-end, je me suis amusé à refaire un TD de première année en Fac d'info, que nous avais proposé notre prof, Olivier Ridoux. L'objectif était, à partir d'un texte en entrée, de fournir un nouveau texte qui "sonnait" comme le premier. Il est vrai qu'à la lecture, le résultat paraît presque lisible, même s'il ne veut absolument rien dire. Je me suis inspiré de la méthode pour tenter de reproduire ce qu'il nous était demandé. Pour être honnête, j'ai eu du mal à comprendre la syntaxe utilisée... Le langage de programmation que nous devions utiliser était Scheme. Pourtant, me direz-vous, d'après la page Wikipedia de R (j'étais assez surpris, mais finalement pas tant que ça), la sémantique lexicale de R est inspirée de Scheme !

Alors c'est bien beau tout ça, mais au final, comment est-ce que ça marche ? L'idée est simple ! Il faut se baser sur les probabilités d'apparition de chaque lettre. Ce qui reste encore flou pour moi, c'est à quel point faut-il intégrer l'historique dans le calcul de la probabilité ? Faut-il regarder la probabilité d'apparition de chaque caractère uniquement ? Regarder cette probabilité sachant la lettre qui précède ? Sachant les deux lettres qui précèdent ? etc.

Le TP prenait appui sur les travaux de Shannon, sur la théorie de l'information.
On considère une variable aléatoire \(X\) qui prend ses valeurs dans un ensemble fini, avec des probabilités d'apparition. Ce qu'on cherche à obtenir, c'est la probabilité d'apparition d'une lettre, en sachant la séquence de \(k\) lettres qui précèdent. On parle de n-gram. Pour évaluer cette probabilité, on considère qu'il suffit de connaître les \(n-1\) éléments précédents. On a donc :

\(\mathbb{P}(x_i \mid x_1, x_2, \ldots, x_{i-1}) = \mathbb{P}(x_i \mid x_{i-(n-1)}, x_{i-(n-2)}, \ldots, x_{i-1}) \).

Pour les besoins, j'ai récupéré 20 poèmes de Baudelaire sur ce site : http://baudelaire.litteratura.com/?rub=oeuvre&srub=pov&id=11.
Au total, j'ai un échantillon de 20794 caractères (espaces comprises).

La première fonction set à calculer les n-grams :

calculerNGrams <- function(n,texte){
# Table des frequences
freq.letters <- data.frame(unclass(textcnt(texte,n=n,split=NULL)))
colnames(freq.letters) <- list("occ")
freq.letters$letter <- rownames(freq.letters)
rownames(freq.letters) <- NULL
# On enleve le caractere "_" de la chaine
aRetirer <- grep("_", freq.letters$letter)
if(length(aRetirer)>0) freq.letters <- freq.letters[-aRetirer,]
return(freq.letters)
}

La fonction textcnt du package tau est bien pratique pour compter les occurrences de chaque séquence dans le texte. On définit la longueur de la séquence maximale. Comme ce que je veux, ce sont uniquement les séquences de n caractères, il faut effectuer un petit filtrage :


calculerNGrams.reduits <- function(n, freq.letters){
# On ne conserve que les n-grams
freq.letters <- freq.letters[which(nchar(freq.letters$letter)==n),]
# La frequence
freq.letters$freq <- freq.letters$occ / sum(freq.letters$occ)
# On ordonne le tout
freq.letters <- freq.letters[order(freq.letters$occ,decreasing=FALSE),]
return(freq.letters)
}

Les occurrences d'apparition permettent de calculer les fréquences normalisées.

Ensuite, une petite fonction pour tirer une lettre aléatoirement en pondérant la probabilité d'apparition de chaque lettre par la fréquence observée dans l'ensemble du texte :


tirerUneLettre <- function(ngram) sample(size=1,x=ngram$letter,prob=ngram$freq)

Une dernière fonction d'aide au calcul est créée :


# Afin de generer une lettre supplementaire
genererLettreSupple <- function(resul,n){
i=n
prefixe <- substr(resul,nchar(resul)-(n-2),nchar(resul))
lesiGrams <- get(paste("les",i,"Grams",sep=""))
# Si on ne trouve pas le prefixe dans les n-grams, on va cherche dans le (n-1)-grams
while(i>1 & !length(which(substr(lesiGrams$letter,1,i-1)==prefixe))>0){
lesiGrams <- get(paste("les",i-1,"Grams",sep=""))
i=i-1
prefixe <- substr(resul,nchar(resul)-(i-2),nchar(resul))
}
lesiGrams.temp <- lesiGrams[which(substr(lesiGrams$letter,1,(i-1))==prefixe),]
lettreSup <- tirerUneLettre(ngram= lesiGrams.temp)
lettreSup <- substr(lettreSup,nchar(lettreSup),nchar(lettreSup))
return(lettreSup)
}

À partir d'une chaine de caractère (resul), et d'un nombre de lettre à utiliser pour l'historique (n), une nouvelle lettre est proposée, en se basant sur la probabilité d'apparition.

Enfin, une dernière fonction pour générer le texte :


genererTexte <- function(texte,taille,n){
# Creation de la table des sequences
freq.letters <- calculerNGrams(n,texte)
# Creation des tables de frequence pour chaque n-gram.
creerNgram <- function(i){
assign(paste("les",i,"Grams",sep=""), calculerNGrams.reduits(i, freq.letters = freq.letters),envir=.GlobalEnv)
return(NULL)
}
lapply(1:n, creerNgram)
# INITIALISATION
resul <- tirerUneLettre(ngram=les1Grams)
for(i in 2:n) resul <- paste(resul,genererLettreSupple(resul,i),sep="")
# Fin de l'initialisation
while(nchar(resul)<=taille){
resul <- paste(resul,genererLettreSupple(resul,n=n),sep="")
}
return(resul)
}

En fournissant le texte d'entrée (les poèmes de Baudelaire), la taille du texte souhaité et la taille de l'historique sur lequel baser les calculs, la fonction retourne un joli petit texte !

Dans un premier temps, on utilise la fonction calculerNGrams pour obtenir l'apparition de chaque séquence de 1, 2, ..., n caractères. Ensuite, les tables de fréquences sont créées pour chacune des n longueurs de séquence. Pour obtenir le résultat, il est nécessaire de construire une chaine atteignant une longueur n avant de pouvoir se baser sur l'historique complet. Au lieu de tirer aléatoirement n caractères, on se base sur l'historique présent, même s'il n'est pas assez grand.

L'inconvénient de ce code, c'est qu'il utilise des boucles, et de fait, les calculs sont assez longs... Je vous laisse malgré tout admirer la poésie que R nous retourne :


cat(genererTexte(1000,n=5,texte=baudelaire))

l'on ne ses bellit le chagrins

que leur riches et mauvais je neigeuses légers le ciel, le carnaval où l'empêcher.

 au-desses bois, de pieds des azur, front le viande son des temps et formes maternels, le de ces dominations, il eut d'avoir

le fils mêlent de sève et ses de serein,

nous de les hiver tourse à la mers, avait l'un soleils, des pitié :

- ô douloureux des vous communesseins,

lassez béni, poète et la plus méchanter de ces farces esprit quand il est perdue.

dès lors les vallées,

vers les sainte vérité,

et que, pare lutin,

et son desse,

occupentions et vieille soirées,

tandis qu'il noyée au poètes,

vers les huées

aux la prophétique

le sillant volons, comme dans amers,

son homme,

comme dans l'armure, un agace longer au pied d'un fœtus que et tout un s'exhalant des de ce marcheront qui se comique n'ai-je donne mystique nous repentirs s'exhaleurs ardent.

il jourse à ce cœurs soleils.

les habite sereine les éternels, ô frères

qui chasse, l'invisibles de sa rate des horreur fascin

samedi, avril 6 2013

Carte R interactive

Depuis un moment, j'avais envie de créer une version un peu animée des cartes R. Il y a bien la possibilité d'exporter les cartes au format gif, mais il manque le côté où on devient acteur...

En cherchant sur la toile, on tombe rapidement sur le tutoriel de Sylvain Prévost qui est plutôt bien expliqué. En parcourant les commentaires, un certain Romain propose son lien avec une petite amélioration : l'affichage de résultats.

Ce qui me plaît dans la méthode employée pour faire cette carte, c'est que ça se base sur un fichier SVG. Or, bien évidemment, R propose d'exporter les graphiques au format SVG... Aussi, je vous propose aujourd'hui, à partir de R et d'un peu de bidouille, de créer une visualisation interactive (un peu) de la carte de France. Mes compétences en JavaScript étant nulles, le résultat n'est pas exceptionnel, mais il est facilement reproductible. Merci à @YoruNoHikage pour ses conseils avisés en javascript et en css.

Voici ce que ça donne : http://editerna.free.fr/blog/cartes/carte_france/.

Un aperçu statique est présenté ci-dessous (impossibilité d'insérer un script dans un billet...).

Ce dont on a besoin :

  • le logiciel R, avec les packages :
    • maps,
    • maptools,
    • RColorBrewer ;
La première chose à faire, est de créer un dossier dans lequel figurent :
Le fichier "creationJS.R", comme vous vous en doutez, va servir à créer le fichier JavaScript pour pouvoir afficher la carte. Il contient deux fonctions :
  1. creer.fichier.js, qui prend en entrée 6 paramètres. Le premier, 
    file, est le nom du fichier .js que l'on souhaite créer. Le second, valeurs, renseigne les valeurs que l'on souhaite assigner à chaque département. Le troisième, informations.sup, qui contient d'autres informations que l'on souhaite afficher en dehors de l'infobulle, dans la partie droite. Le quatrième, couleurs, indique les couleurs pour chaque département. Le cinquième, legende, renseigne les labels de la légende. Enfin, legende.lesCouleurs donne les couleurs des éléments de la légende ;
  2. creer.html, qui prend 3 arguments. Le premier, file.html, renseigne le nom du fichier HTML à créer. Le second, file.css indique le nom du fichier CSS qui contiendra les propriétés de styles. Le dernier, titre, est simplement pour donner un titre à la page.
Je ne rentre pas dans le détail de ces fonctions. Si vous avez des questions, n'hésitez pas. Je préfère m'attarder sur le deuxième fichier R, celui qu'il faut modifier pour obtenir une nouvelle carte.
Dans un premier temps, on se place dans le répertoire qui contient les deux fichiers R et le dossier js :

setwd("leDossierRacine")
source("creationJS.R")

Puis, on charge la carte de France


require(maps)
france <- map("france",fill=TRUE,plot=FALSE)

À présent, il faut récupérer les noms des départements. Ils vont servir à l'affichage, mais également comme nom de variable. C'est pourquoi je retire tout ce qui suit après les ":" dans les noms proposés par R.


lesDepartements.unique <- unique(sapply(strsplit(france$names,":"),function(x) x[1]))
lesDepartements <- sapply(strsplit(france$names,":"),function(x) x[1])

On note qu'il y a deux variables. En effet, on va avoir besoin de conserver les doublons, parce que derrière des noms similaires se cachent des polygones différents. Par exemple, pour toutes les îles.

Pour fournir un exemple facile à reproduire, tirons des valeurs pour chaque département selon une loi Uniforme de paramètres 0 et 30. Faisons de même pour les informations supplémentaire.

Dans la version présentée en introduction, je me suis amusé avec du texte généré à partir des poèmes de Baudelaire. On verra ça dans un prochain billet.


valeurs.temp <- round(runif(n=length(lesDepartements.unique),min=0,max=30))
informations.sup.temp <- round(runif(n=length(lesDepartements.unique),min=0,max=30))

On ordonne le tout dans une data frame :


df.departements.temp <- data.frame(dep.nom=lesDepartements.unique, valeur=valeurs.temp, informations.sup= informations.sup.temp,stringsAsFactors=FALSE)

Pour que les doublons aient droit, eux aussi, à leur valeur, ne les oublions pas :


valeurs <- df.departements.temp[match(lesDepartements, lesDepartements.unique),"valeur"]
informations.sup <- df.departements.temp[match(lesDepartements, lesDepartements.unique),"informations.sup"]

Maintenant, penchons-nous sur la légende. On choisit 5 couleurs et on scinde l'éventail des valeurs possibles en classes de même longueur.


nb.colors <- 5
require(RColorBrewer)
legende.lesCouleurs <- rev(brewer.pal(nb.colors,"YlOrRd"))
ratio <- (max(valeurs)-min(valeurs))/nb.colors
brks <- round(seq(min(valeurs),max(valeurs),by=ratio),digits=0)
couleurs <- legende.lesCouleurs[findInterval(valeurs,brks,all.inside=TRUE)]
require(maptools)
legende <- leglabs(brks,under="inférieur à",over="supérieur à")

Dans une dernière étape, on appelle les fonctions qui vont créer le script js, le fichier html et la feuille de style.


creer.fichier.js(file="js/script.js",valeurs, informations.sup,couleurs,legende, legende.lesCouleurs)
creer.html(file.html="index.html",file.css="style_carte.css",titre="Carte de la France")

Et voilà ! Il ne reste plus qu'à lancer le fichier "index.html".

jeudi, avril 4 2013

Points dans un cercle

Aujourd'hui, j'ai dû déterminer quels étaient les points géographiques (des stations météorologiques) se situant dans un cercle de centre et de rayon donnés. L'année dernière, j'avais tenté quelque chose, mais le résultat ne me plaisait pas vraiment. Je suis tombé sur une fonction qui m'avait échappé : gcDestination, du pakage maptools. Elle permet de retourner les coordonnées d'un point selon une distance et un angle donné.

Trois packages sont nécessaire pour pouvoir exécuter le code qui suit : 

library(maps)
library(maptools)
library(rgeos)

Pour illustrer le propos, penchons-nous sur la Bretagne.

Donnons-nous 10 points en Bretagne :

longitudes <- c(-2.373247, -2.523580, -3.112924, -1.474192, -4.151028, -4.367183, -2.143235, -1.528023, -3.733716, -3.311343)
latitudes <- c(47.69491, 48.30560, 48.39880, 48.24052, 48.03038, 48.53544, 48.49801, 48.01217, 47.87922, 48.75708)
lesPoints <- data.frame(lon=longitudes,lat=latitudes)

À présent, récupérons la carte de la France. J'utilise les méthodes que j'ai apprises sur le tutoriel de Baptiste Coulmont (section 2) et sur le billet d'Arthur Charpentier.

# La carte de la France
france <- map('france', plot=FALSE,fill=TRUE)
detpartement_bzh <- france$names[c(25,30,32,41)] #bretagne
# La carte de la Bretagne
bretagne <- map('france',regions=detpartement_bzh, fill=TRUE, col="transparent", plot=FALSE)
# Les noms des departements
# c.f. http://stackoverflow.com/questions/13316185/r-convert-zipcode-or-lat-long-to-county
IDs <-sapply(strsplit(bretagne$names,":"),function(x) x[1])
# Convertir en objet SP
bretagne.sp <- map2SpatialPolygons(bretagne,IDs=IDs,proj4string=CRS("+proj=longlat +datum=wgs84"))

On va repérer des points particuliers pour chaque département. C'est l'occasion d'introduire la fonction gCentroid, du package rgeos, qui permet de trouver le barycentre d'un polygone.

barycentre <- gCentroid(bretagne.sp,byid=TRUE)

Le rendu visuel est le suivant :

Barycentre Bretagne

La fonction qui vient sert à créer 360 points autour d'une référence.

# A partir d'un point (unPoint) repere par ses coordonnees (longitude, latitude), cree un polygone circulaire de rayon donne (distance)
cercleDistant <- function(unPoint,distance){
# Creation de 360 points distincts sur le cercle de centre
# unPoint et de rayon distance
do.call("rbind",
lapply(0:360,function(bearing){
res <- gcDestination(
lon=unPoint[1],
lat=unPoint[2],
bearing=bearing,
dist=distance,
dist.units="km",
model="WGS84")
rownames(res) <- NULL
return(data.frame(res))
})
)
}

Comme l'idée n'est de donner un exemple, on isole le barycentre des Côtes d'Armor :

# L'indice du barycentre des Cotes d'Armor
indice.armor <- which(rownames(coordinates(barycentre))%in%"Cotes-Darmor")
# Les coordonnees du barycentre des Cotes d'Armor
barycentre.armor <- coordinates(barycentre[indice.armor])

Créons un cercle de rayon 40km autour de ce point :

# Le cercle de centre barycentre.armor et de rayon 40km
cercle.armor <- cercleDistant(barycentre.armor,distance=40)

Il ne reste plus qu'à repérer les points de la data.frame initiale (celle qui contient les 10 couples de coordonnées) qui sont dans ce cercle :

# Indice des points qui sont dans le cercle cercle.armor
indices.points.armor <- which(point.in.polygon(lesPoints$lon,lesPoints$lat,cercle.armor$long,cercle.armor$lat)==1)
# Les points dans le cercle
points.cercle.armor <- lesPoints[indices.points.armor,]

Et pour finir, un petit affichage graphique :

map(bretagne)
points(barycentre.armor,col="red",pch=19)
points(cercle.armor,col="dodger blue",cex=.5,pch=1)
points(lesPoints,col="gold",pch=19)
points(points.cercle.armor,col="purple",pch=19)

jeudi, mars 14 2013

Gestion des dates sous R

Un pense bête pour la gestion des dates sous R.

Dans le cas où les dates sont stockées dans une data.frame au format character, il suffit d'utiliser la fonction strptime en indiquant le bon format de la date. Il se peut que les options de la machine viennent empêcher le bon fonctionnement de cette fonction. Aussi, il est possible de changer les paramètres d'internationalisation. Les lignes de code qui suivent parlent d'elles-mêmes.


Sys.getlocale() # Nos parametres
Sys.setlocale('LC_ALL','C') # utiliser les parametres par defaut pour le langage C
USDates <- data.frame(lesDates=c("2012-07-08 00:41:01",
"2011-08-08 12:34:28",
"2012-09-15 01:17:03",
"2012-07-14 22:25:18"),stringsAsFactors=FALSE)
USDates$lesDates2 <- strptime(USDates$lesDates,"%Y-%m-%d %H:%M:%S")
sort(USDates$lesDates2,decreasing=TRUE)
(USDates$lesDates2)
strftime(USDates$lesDates2,"%Y") #extraire l'annee
strftime(USDates$lesDates2,"%m") #extraire le mois
strftime(USDates$lesDates2,"%d") #extraire le jour
Sys.setlocale(locale = "") # restaurer les parametres locaux

mercredi, octobre 31 2012

Sélection de modèles à l'aide des critères AIC et BIC

Dans le cadre d'un projet pour un cours de prévisions, je suis amené à sélectionner différents modèles pour prévoir des valeurs futures par la suite. Dans le cours, on nous parle d'une procédure sur Winrats, qui donne les AIC et les BIC pour tous les modèles ARMA possibles, en ayant au préalable fixé les valeurs maximales des retards.

Soyons clairs, il est entendu que cette méthode de choix de modèle est vraiment naïve. Mais elle peut être un début. Sous R, il existe la fonction auto.arima (package forecast) qui retourne — dans des temps qui m'épatent — le "meilleur" modèle, en fonction du critère qu'on passe en argument. Mais le résultat obtenu ne me plait pas vraiment. J'aimerais voir s'il existe d'autres modèles possibles, donnant des critères de choix très proches, avec des modèles plus simples.

Avant de se lancer dans le code, un petit rappel sur les critères AIC (introduit par Akaike, 1974) et BIC (Schwarz, 1978).

\(AIC=\color{green}{2k} - 2\log(L)\) avec \(k\) le nombre de paramètres et \(L\) le maximum de la fonction de log-vraisemblance.

\(BIC=\color{green}{\log(n)k}-2\log(L)\) avec \(n\) le nombre d'observations.

On remarque que lorsque \(n>\exp(2)\), le critère AIC donne une pénalité liée au nombre de paramètres moins importante que le BIC. Les modèle choisis selon le BIC ont donc tendance à être plus parcimonieux.

Avant de décortiquer le code, voici le genre de résultat que je souhaite obtenir :

$AIC
    Q=0       Q=1         Q=2       Q=3       Q=4       Q=5      
P=0 "1789.42" "1630.21"   "1525.85" "1493.52" "1479.2"  "1461.4" 
P=1 "1486.53" "1446.86"   "1442.9"  "1443.33" "1445.04" "1446.65"
P=2 "1439.53" "1441.52"   "1443.52" "1445.2"  "1440.85" "1441.14"
P=3 "1441.52" "*1436.52*" "1437.29" "1438.82" "1440.69" "1442.59"
P=4 "1443.52" "1437.59"   "1443.02" "1440.78" "1442.77" "1445.74"
P=5 "1444.75" "1438.68"   "1440.67" "1442.07" "1440.65" "1442.73"

$BIC
    Q=0         Q=1       Q=2       Q=3       Q=4       Q=5      
P=0 "1793.64"   "1638.64" "1538.5"  "1510.38" "1500.27" "1486.69"
P=1 "1494.96"   "1459.5"  "1459.76" "1464.4"  "1470.33" "1476.15"
P=2 "*1452.17*" "1458.38" "1464.59" "1470.48" "1470.35" "1474.86"
P=3 "1458.38"   "1457.59" "1462.58" "1468.32" "1474.4"  "1480.52"
P=4 "1464.59"   "1462.88" "1472.53" "1474.5"  "1480.7"  "1487.89"
P=5 "1470.03"   "1468.18" "1474.38" "1480.01" "1482.8"  "1489.09"

Alors, on a besoin de deux packages : forecast et gregmisc.

L'idée est d'estimer chaque modèle, en s'étant donné les retards P (auto-régression) et Q (moyenne-mobile) maximum, puis de calculer les deux critères AIC et BIC.

criteres <- function(data,P,D,Q){
mod <- try(arima(as.numeric(data),order=c(P,D,Q),include.mean=FALSE,method="ML"))
if(is.null(names(mod))){
print(paste("Probleme pour ARIMA(",P,",",D,",",Q,")",sep=""))
return(c(mod=paste("ARIMA(",P,D,Q,")",sep=""),aic=NA,bic=NA))
}else{
n <- length(mod$residuals)
npars <- length(mod$coef)+1
loglik <- mod$loglik
aic <- -2*loglik+2*npars
bic <- -2*loglik+log(n)*npars
return(c(mod=paste("ARIMA(",P,D,Q,")",sep=""),aic=aic,bic=bic))
}
}

On va donc utiliser cette fonction pour chacun des modèles. Il est donc nécessaire de les énumérer. Il doit exister une fonction miracle pour créer ça un peu plus proprement. J'avais pensé utiliser les combinaisons de deux dans un vecteur, mais il manque les modèles ARMA(1,1), ARMA(2,2,), etc.


# Retourne les valeurs des criteres AIC et BIC pour une serie donnee
AICBIC=function(serie,maxP,maxQ,round=2){
matriceModeles=matrix(ncol=2)
for(i in 0:maxP){
for(j in 0:maxQ){
matriceModeles=rbind(matriceModeles,c(i,j))
}
}
matriceModeles=matriceModeles[-1,]
res=apply(matriceModeles,1,function(x) criteres(serie,P=x[1],D=0,Q=x[2]))
lesAIC=round(creerMatriceResul(maxP,maxQ,res,"aic"),round)
lesAIC[which(lesAIC==min(lesAIC,na.rm=TRUE),arr.ind=TRUE)]=paste("*",lesAIC[which(lesAIC==min(lesAIC,na.rm=TRUE),arr.ind=TRUE)],"*",sep="")
lesBIC=round(creerMatriceResul(maxP,maxQ,res,"bic"),round)
lesBIC[which(lesBIC==min(lesBIC,na.rm=TRUE),arr.ind=TRUE)]=paste("*",lesBIC[which(lesBIC==min(lesBIC,na.rm=TRUE),arr.ind=TRUE)],"*",sep="")
return(list(AIC=lesAIC,BIC=lesBIC))
}

Avec ça, on a tout ce qu'il faut. Maintenant, un petit exemple. Je simule deux ARMA différents. Le premier est un AR(2) et le second un ARMA(2,3).

set.seed(1)
n=100
AR2=arima.sim(model=list(ar=c(.5,.3)),n)
ARMA23=arima.sim(model=list(ar=c(.5,.3),ma=c(0.8,0.1,-0.1)),n)
AR(2) et ARMA(2,3)
Si on regarde les sorties de la fonction d'autocorrélation et d'autocorrélation partielle, on peut effectivement penser à un AR(2) pour le premier. On a bien une décroissance vers zéro de la fonction d'autocorrélation, et on observe une chute à zéro au 3e retard.
ACF et PACF AR(2)
Pour la second série, on rejette l'idée d'un processus auto-régressif pur, avec l'alternance de signe des autocorrélations partielles.
ACF et PACF ARMA(2,3)

Là où la fonction auto.arima est très rapide, ma fonction est quant à elle beaucoup plus gourmande en ressources. En fixant des ordres de lags max à 5, voici les résultats :

Vrai modèle auto.arima AICBIC
AR(2) AIC BIC AIC BIC
ARMA(1,2) ARMA(2) ARMA(3,1) AR(2)
ARMA(2,3) AIC BIC AIC BIC
ARMA(1,1) ARMA(1,1) ARMA(4,1) ARMA(1,1)

Ce petit exemple me montre qu'avec 500 observations, on a du mal à trouver le bon modèle avec les critères AIC et BIC... Maintenant je vois un peu mieux pourquoi un autre prof riait aux grands éclats quand on lui parlait de ces critères pour différencier nos modèles.

Pour ceux qui sont malgré tout intéressés, voici le code R.

samedi, octobre 20 2012

Montreal border

After reading a really nice post about geo-marked bike accidents (Mapping Bike Accidents in R), I wanted to reuse the map for something else. Unfortunately, the shapefile used to create it does not give the entire border, but an aggregation of borough instead.

So, after a couple hour looking on the web for a simpler SHP file, I gave up and finally decided to create it myself. To do so, I used the SHP file given on the Mapping Bike Accidents post.

Montreal

First, let's load data:


mtl<-readShapePoly('montreal_borough_borders.shp')

The following lines are not elegant at all, but I didn't know any better way to do the job.

So, basically, what is done here is quite simple:

  1. find every borough that has any border delimiting Montreal ;
  2. keep only the points I need

Since the object mtl is a spatial polygon data frame, I have to use this kind of syntax to get to the coordinates of the boroughs :

mtl@polygons[[1]]@Polygons[[1]]@coords

Here are some functions I need to find beginning and end points of the borough that are on the external frontier :

# plots the first borough in green, the second in blue, and each points from lIntersect
raccords=function(borough1,borough2,lIntersect){
plot(mtl)
lines(borough1,col="light green")
lines(borough2,col="dodger blue")
points(lIntersect,col=heat.colors(nrow(lIntersect)))
}
# Returns the position of lIntersect in borough
findPosition=function(borough,intersectPoint){
which(borough[,1]==intersectPoint[1] & borough[,2]==intersectPoint[2])
}
# Help function to establish union between two bourough
raccords2=function(borough1,borough2){
plot(rbind(borough1,borough2),type="l")
lines(borough1,col="light green")
lines(borough2,lwd=2,col="dodger blue")
}

The idea is simple: the first point of a given borough should be the last of another one. This is what we try to find with the function raccords.

coords1=mtl@polygons[[29]]@Polygons[[1]]@coords
coords2=mtl@polygons[[15]]@Polygons[[1]]@coords
lIntersect=matrix(intersect(coords1,coords2),ncol=2)

The last instruction gives us 21 points. Let's plot this.

intersection of two boroughs

The point we want is the one on the West. Since colors are picked from a list created by the functionheat.colors(.), red ones arrive on top of the list. So the point we want it the first of the list. We only have to find it in both boroughs :

end1=findPosition(coords1,lIntersect[1,]) #end of coords1
beg2=findPosition(coords2,lIntersect[1,]) #beginning of coords2

For some districts, I encountered some issues.

intersection boroughs 3 4

As you can see on the map, there is no intersection near the river (on the right). This is why I used the other raccord function, which helps finding the point in an heuristic way.

Once I found a beginning and an ending point for each district, I checked for possible error :

check=function(i){
plot(mtl)
coords=get(paste("coords",i,sep=""))
beg=get(paste("beg",i,sep=""))
end=get(paste("end",i,sep=""))
points(coords[beg:end,],col="red",cex=.1)
}
par(mfrow=c(4,4))
lapply(1:16,check)

The output stresses some issues. For some boroughs, what I expected to be an ending point was actually a beginning point. Hence, I had to use two different ways to finally create my polygon of Montreal.

assignCoord=function(i){
coords=get(paste("coords",i,sep=""))
beg=get(paste("beg",i,sep=""))
end=get(paste("end",i,sep=""))
assign(paste("co",i,sep="_"),coords[beg:end,],envir=globalenv())
}
lapply(c(1:2,4,7,10:15),assignCoord)
assignCoord2=function(i){
coords=get(paste("coords",i,sep=""))
beg=get(paste("beg",i,sep=""))
end=get(paste("end",i,sep=""))
assign(paste("co",i,sep="_"),coords[c(beg:nrow(coords),1:end),],envir=globalenv())
}
lapply(c(3,5,6,8,9,16),assignCoord2)
montreal=co_1
for(i in 2:17){
montreal=rbind(montreal,get(paste("co",i,sep="_")))
}

So here it is, my polygon of Montreal! If you want to check my R code, feel free.

Montreal as a Polygon

- page 1 de 2