Category Archives: Maps

Quel trajet minimal pour récupérer le maximum de bonbons à l’Halloween ? (partie 2)

Continuons la discussion de ce matin, avec l’exemple des capitales américaines. On peut trouver une telle liste liste en ligne sur wikipedia,

> capital=read.table(
+ "http://freakonometrics.free.fr/US-capital.csv",
+ header=TRUE,sep=";")
> capital=capital[-c(2,11),]

On va enlever Hawaï et l’Alaska, qui forcent à faire de trop grands détours. Dans TSP, sous R, il y a des listes de villes américaines, avec les coordonnées. En bricolant un peu, on va récupérer les coordonnées de la plupart des capitales,

> library(maps)
> library(maptools)
> library(sp)
> library(TSP)
> Lcap=paste(as.character(capital[,1]),
+ as.character(capital[,2]),sep=", ")
> M=as.matrix(USCA312)
> L=labels(M)[[1]]
> Lcap=Lcap[which(Lcap%in%L)]
> k=which(L%in%Lcap)
> M=M[k,]
> M=M[,k]
> listeUSA=TSP(M,L[k])

Et cette fois on est prêt,

> tour=solve_TSP(listeUSA, method = "nn")

on a lancé l’algorithme, et on peut faire un dessin,

> COORD=data.frame(USCA312_coords[k,])
> COORD=COORD[as.numeric(tour),] 
> map('state',fill=TRUE,col="yellow")
> lines(COORD$coords.x1,COORD$coords.x2,
+ lwd=3,col="blue")
> points(COORD$coords.x1,COORD$coords.x2,
+ pch = 19, cex = 1.2, col = "red")

pas mal, n’est pas ? Bon, mais comme d’habitude, rien ne garantie que l’on ait trouvé le trajet optimal… En fait, le soucis est que si on relance l’algorithme, on retombe sur un autre trajet,

voire encore un autre si on lancer une nouvelle fois l’algorithme,

C’est pénible, n’est ce pas ? En fait, si on lance 1,000 fois l’algorithme on obtient la distribution (des distances optimales) suivante

autrement dit, on a encore beaucoup de variabilité, avec une distance optimale qui peut varier d’environ 10% entre deux boucles d’algorithmes. Cela dit, cela reste faible par rapport à un trajet au hasard,

Et d’ailleurs, si on compare la distance obtenue en tirant les villes au hasard, on arrive à un taux de l’ordre de 27% (qui peut être comparé à ce que nous avions obtenu en tirant au hasard dans le carré unité)

>  HAZARD=OPTIMAL=RATIO=rep(NA,500)
>  for(s in 1:500){
+   HAZARD[s]=tsp.longueur(M,sample(1:nrow(M)))
+   tour=solve_TSP(listeUSA, method = "nn")
+   OPTIMAL[s]=tsp.longueur(M,as.numeric(tour))
+   RATIO[s]=HAZARD[s]/
+            OPTIMAL[s]
+ }

Bref, compte tenu de la (grande) dispersion des résultats optimaux obtenu, on pourrait dire que cet algorithme est assez mauvais. Même s’il est très rapide. Un autre avantage est que l’on peut spécifier la ville de départ, par exemple. Ici, si on souhaite partir de Floride, on peut obtenir le trajet suivant

> tour=solve_TSP(listeUSA, method = "nn",
+ control=list(start=41))

Passionnant n’est-ce pas ? Et encore, ça ne va pas m’aider à trouver le parcours optimal pour les enfants, car pour les enfants, il faut spécifier le point de départ (voire le point d’arrivée), il faut en plus rester sur la route (voire traverser aux passages pour piétons). On peut aussi vouloir faire une pause pipi au milieu de la promenade…. Bref, on a pas mal de contraintes ! Mais tout est expliqué dans le livre fabuleux de William Cook, In Pursuit of the Traveling Salesman, qui reprend toute l’histoire des algorithmes les plus fins pour résoudre ce joli problème de maths…

Quel trajet minimal pour récupérer le maximum de bonbons à l’Halloween ? (partie 1)

Bon, ben cette fois ça y est, c’est l’Halloween… et comme tous les ans, mon rôle devrait se limiter à creuser les citrouilles, décorer la maison… et attendre à la porte de la maison que les enfants des voisins viennent me dévaliser ma réserve de bonbons ! Je pourrais à la rigueur revendiquer un petit quelque chose dans le déguisement de mon fils (on a passé pas mal de temps à feuilleter les premiers tomes de Walking Deads, histoire de mieux saisir l’essence des personnages de morts-vivants, en plus de la Zombie Walk la fin de semaine passée, en centre-ville).

Mais cette année, pas de billet sur les zombies comme l’an passé, mais un vrai problème pratique: comment optimiser son trajet, quand on sait où sont les bonnes maisons, où on a de fortes chances d’avoir plein de bonbons (en particulier chez les parents des amis) ? C’est assez proche de l’idée de passer par un maximum de maisons, et ce à une heure raisonnable (histoire qu’il reste des bonbons à récupérer).

Pour tous ceux qui ont fait un peu de programmation dynamique, c’est un problème classique, connu sous le nom du “voyageur de commerce“. Ce problème est né en juillet 1954 avec “the shortest route for a traveling salesman” paru dans Newsweek ci-dessous (que l’on retrouve en détails dans le livre 50 Years of Integer Programming)

Enfin, disons qu’il a été popularisé à cette époque, car Euler évoquait déjà un problème similaire quand il travaillait sur le problème des sept ponts de Königsberg, dans “Solutio problematis ad geometriam situs pertinentis“,

(même si ce problème consiste à étudier l’existence d’un trajet sous des contraintes de passages par des nœuds, et pas vraiment de l’optimisation). On retrouve aussi un exemple en 1832, en Allemagne, dans un livre intitulé Der Handlungsreisende – wie er sein soll und was er zu thun hat, um Auftraege zu erhalten und eines gluecklichen Erfolgs in seinen Geschaeften gewiss zu sein – Von einem alten Commis-Voyageur où un trajet passant par 45 villes en Allemagne et en Suisse devait être optimisé

(le schéma ci-dessus montre celle décrite en 1832, pour un total de 1,285km, mais il est possible de faire une boucle de 1,248km).

En fait le problème a explosé en 1962 (à cette époque, il semble que les gens pouvaient encore penser que faire des maths pouvait être le fun), lors que Procter & Gamble ont offert un prix de 10,000$ à celui qui trouverait le chemin le plus court pour que Toody et Muldoom, les conducteurs de la voiture 54 dans une série populaire à l’époque, passent par 33 villes, aux États-Unis.

Cela dit, Carl Menger avait noté dès 1930 a noté que ce problème devrait etre difficile à résoudre. Facile à écrire, mais difficile à résoudre explicitement. Dans son livre, William Cooke rappelle qu’avec 33 villes, il y aurait

131,565,418,466,846,765,083,609,006,080,000,000

chemins possibles, dont il faudrait calculer la longueur. Et qu’avec l’ordinateur le plus puissant en 2009 (effectuant 1.5 million de milliards d’opérations à la seconde) il faudrait 28,000 milliards d’année pour mener à bien les calculs. Même si on essayait de tenir compte des progrès de l’informatique dans les mille prochaines années, ça risque de prendre du temps !..

En 1967, Jack Edmonds notait “I conjecture that there is no good algorithm for the traveling salesman problem” où “good” est une manière élégante pour dire un temps polynomial en , le nombre de villes. Ce qui s’appelle les problèmes de la classe . Il semblerait que le problème du voyageur de commerce soit un problème  et pas de la classe (cf le chapitre du livre de Sanjeev Arora et Boaz Barak, ou le chapitre 11 des notes de cours de Olivier Bournez, que j’ai pu découvrir grâce à @dmonniaux, ou formellement de la classe  et pas , i.e.  complet), voir aussi la discussion sur wikipedia),

Et depuis, l’institut Clay a proposé $1,000,000 à celui qui prouverait que les problèmes   ne sont pas de type  (ou qui prouveraient qu’ils le sont… laconjecture reste ouverte, semble-t-il). Bref, les enjeux financiers derrière ont considérablement augmenté (on pourra consulter Sipser (2007) pour la petite histoire)

Bon, et si l’histoire est intéressante, ce n’est pas une raison pour la lire en se croisant les bras… essayons de programmer un algorithme (simple) afin d’approcher une solution, car je ne pourrais me contenter de dire à mes enfants “c’est un problème  complet, papa a capitulé“.

Commençons par tirer 5 points au hasard dans le carré unité,

>  n = 5
>  (x = matrix(runif(2*n),nr=n))
[,1]       [,2]
[1,] 0.4142621 0.35600725
[2,] 0.9477231 0.66938653
[3,] 0.2518656 0.19900873
[4,] 0.5314231 0.22954670
[5,] 0.5156522 0.09066741
et commençons par calculer la distance entre les points (i.e. une matrice, ou ici une matrice triangulaire inférieure)
>  (d = dist(x))
1         2         3         4
2 0.6186980
3 0.2258786 0.8399243
4 0.1723919 0.6056111 0.2812205
5 0.2840514 0.7222195 0.2851688 0.1397719
On va mettre ce triangle sous forme de matrice (pleine), ça sera plus simple à manipuler par la suite
>  (d = as.matrix(d))
1         2         3         4         5
1 0.0000000 0.6186980 0.2258786 0.1723919 0.2840514
2 0.6186980 0.0000000 0.8399243 0.6056111 0.7222195
3 0.2258786 0.8399243 0.0000000 0.2812205 0.2851688
4 0.1723919 0.6056111 0.2812205 0.0000000 0.1397719
5 0.2840514 0.7222195 0.2851688 0.1397719 0.0000000
Le principe, pour initialiser l’algorithme est de tirer au hasard un trajet, c’est à dire une permutation de nos points,
>  (o = sample(1:n))
[1] 2 4 1 5 3
Il faut alors calculer le 4 distances entre ces 5 villes (on ne prévoit pas de retour à la ville de départ, pour l’instant). Le plus simple est de jouer avec cette matrice de la réarranger,
>  d[o[1:(n-1)],][,o[2:n]]
4         1         5         3
2 0.6056111 0.6186980 0.7222195 0.8399243
4 0.0000000 0.1723919 0.1397719 0.2812205
1 0.1723919 0.0000000 0.2840514 0.2258786
5 0.1397719 0.2840514 0.0000000 0.2851688
et de prendre les 4 valeurs sur la diagonale
>  diag(d[o[1:(n-1)],][,o[2:n]])
[1] 0.6056111 0.1723919 0.2840514 0.2851688
Ce sont les 4 distances que l’on cherche, entre les 5 villes. La distance totale parcourue est alors
>  sum(diag(d[o[1:(n-1)],][,o[2:n]]))
[1] 1.347223
L’idée est ensuite de permuter deux des villes
>  (i=sample(1:n,2))
[1] 3 1
(ici la première et la troisième de la liste, mais on pourrait décider de fixer la première si on doit partir de là) et de voir si la nouvelle distance est plus courte, ou pas,
>  os=o
>  os[i]=o[rev(i)]
>  sum(diag(d[os[1:(n-1)],][,os[2:n]]))
[1] 1.785391
La distance étant plus longue, on ne va pas permuter. On pourrait tenter 1000 fois de permuter deux villes. Si on ne parvient pas à réduire la distance totale, c’est qu’on a atteint une valeur intéressante (à défaut de pouvoir prouver qu’elle est optimale). Et si on peut améliorer la distance, on permute. Le plus simple est donc de faire des petites fonctions. La première reprend le calcul de la distance totale, sur la diagonale de la matrice permutée, à matrice de distance donnée et à permutation donnée,
> tsp.longueur=function(matrice,ordres) {
+ n=length(ordres)
+ sum(diag(matrice[ordres[1:(n-1)],][,ordres[2:n]])) }
 La seconde est juste une boucle,
> tsp.optimal=function(d,N=1000) {
+ d=as.matrix(d)
+ n=ncol(d)
+ o=sample(1:n)
+ v=tsp.longueur(d,o)
+ k=0
+ while(k < N) {
+   i<-sample(1:n,2)
+   os=o
+   os[i]=o[rev(i)]
+   w=tsp.longueur(d,os)
+   if(w < v) {
+     v=w
+     o=os
+     k=0} else {k=k+1}}
+ list(ordre=o,longueur=v)}
On peut se lancer maintenant, avec par exemple 30 points simulés, sur le carré unité,
>  set.seed(1)
>  n=30
>  x=matrix(runif(2*n),nr=n)
>  o=sample(1:n)
>  plot(x,xlim=0:1,ylim=0:1,xlab="",ylab="")
>  lines(x[o,],col='blue')
>  r=tsp.optimal(dist(x))
>  os=r$ordre
>  lines(x[os,],col='blue',lwd=1.5)
avec la première trajectoire, à gauche, qui reste en trait fin sur le dessin de droite,

Pas mal, non ? Bon, en fait le soucis est que si on joue plusieurs fois avec notre fonction, on obtient toujours un graphique différent,

Autrement dit, je trouve (certes) des chemins pas trop long, mais je suis loin de trouver le chemin optimal. Et histoire de clotûrer le débat, notons qu’on peut faire un peu de simulation, afin de comparer la distance obtenue en faisant le trajet au hasard, et un trajet “optimal“,

> RATIO=HASARD=OPTIMAL=matrix(NA,500,10)
> for(m in 1:10){
+ n=5*m
+ for(s in 1:500){
+  x=matrix(runif(2*n),nr=n)
+  r=tsp.opt(dist(x))
+  HAZARD[s,m]=tsp.longueur(as.matrix(dist(x)),1:n)
+  OPTIMAL[s,m]=r$longueur
+  RATIO[s,m]=HAZARD[s,m]/
+             OPTIMAL[s,m]
+ }
+}

Si on compare la longueur moyenne du trajet (sur 500 scénarios) avec un tour au hasard (en rouge) et ou tour optimisé (en bleu), on a

ce qui donne un ratio moyen (ou une moyenne des ratios, les deux sont représentés)

Autrement dit, avec une dizaine de points éparpillés (au hasard) dans le carré unité, un tour au hasard ne sera que deux fois plus long… étonnant non ?

Mais on doit pouvoir aller plus loin, parce que sous R, il y a des fonctions (et des packages) pour faire des algorithmes plus compliqués… à suivre avant la tournée des bonbons de demain soir !

Border bias and weighted kernels

With Ewen (aka @3wen), not only we have been playing on Twitter this month, we have also been working on kernel estimation for densities of spatial processes. Actually, it is only a part of what he was working on, but that part on kernel estimation has been the opportunity to write a short paper, that can now be downloaded on hal.

The problem with kernels is that kernel density estimators suffer a strong bias on borders. And with geographic data, it is not uncommon to have observations very close to the border (frontier, or ocean). With standard kernels, some weight is allocated outside the area: the density does not sum to one. And we should not look for a global correction, but for a local one. So we should use weighted kernel estimators (see on hal for more details). The problem that weights can be difficult to derive, when the shape of the support is a strange polygon. The idea is to use a property of product Gaussian kernels (with identical bandwidth) i.e. with the interpretation of having noisy observation, we can use the property of circular isodensity curve. And this can be related to Ripley (1977) circumferential correction. And the good point is that, with R, it is extremely simple to get the area of the intersection of two polygons. But we need to upload some R packages first,

require(maps)
require(sp)
require(snow)
require(ellipse)
require(ks)
require(gpclib)
require(rgeos)
require(fields)

To be more clear, let us illustrate that technique on a nice example. For instance, consider some bodiliy injury car accidents in France, in 2008 (that I cannot upload but I can upload a random sample),

base_cara=read.table(
"http://freakonometrics.blog.free.fr/public/base_fin_morb.txt",
sep=";",header=TRUE)

The border of the support of our distribution of car accidents will be the contour of the Finistère departement, that can be found in standard packages

geoloc=read.csv(
"http://freakonometrics.free.fr/popfr19752010.csv",
header=TRUE,sep=",",comment.char="",check.names=FALSE,
colClasses=c(rep("character",5),rep("numeric",38)))
geoloc=geoloc[,c("dep","com","com_nom",
"long","lat","pop_2008")]
geoloc$id=paste(sprintf("%02s",geoloc$dep),
sprintf("%03s",geoloc$com),sep="")
geoloc=geoloc[,c("com_nom","long","lat","pop_2008")]
head(geoloc)
france=map('france',namesonly=TRUE,
plot=FALSE)
francemap=map('france', fill=TRUE, col="transparent",
plot=FALSE)
detpartement_bzh=france[which(france%in%
c("Finistere","Morbihan","Ille-et-Vilaine",
"Cotes-Darmor"))]
bretagne=map('france',regions=detpartement_bzh,
fill=TRUE, col="transparent", plot=FALSE,exact=TRUE)
finistere=cbind(bretagne$x[321:678],bretagne$y[321:678])
FINISTERE=map('france',regions="Finistere", fill=TRUE,
col="transparent", plot=FALSE,exact=TRUE)
monFINISTERE=cbind(FINISTERE$x[c(8:414)],FINISTERE$y[c(8:414)])

Now we need simple functions,

cercle=function(n=200,centre=c(0,0),rayon)
{theta=seq(0,2*pi,length=100)
m=cbind(cos(theta),sin(theta))*rayon
m[,1]=m[,1]+centre[1]
m[,2]=m[,2]+centre[2]
names(m)=c("x","y")
return(m)}
poids=function(x,h,POL)
{leCercle=cercle(centre=x,rayon=5/pi*h)
POLcercle=as(leCercle, "gpc.poly")
return(area.poly(intersect(POL,POLcercle))/
area.poly(POLcercle))}
lissage = function(U,polygone,optimal=TRUE,h=.1)
{n=nrow(U)
IND=which(is.na(U[,1])==FALSE)
U=U[IND,]
if(optimal==TRUE) {H=Hpi(U,binned=FALSE);
H=matrix(c(sqrt(H[1,1]*H[2,2]),0,0,
sqrt(H[1,1]*H[2,2])),2,2)}
if(optimal==FALSE){H= matrix(c(h,0,0,h),2,2)

before defining our weights.

poidsU=function(i,U,h,POL)
{x=U[i,]
poids(x,h,POL)}
OMEGA=parLapply(cl,1:n,poidsU,U=U,h=sqrt(H[1,1]),
POL=as(polygone, "gpc.poly"))
OMEGA=do.call("c",OMEGA)
stopCluster(cl)
}else
{OMEGA=lapply(1:n,poidsU,U=U,h=sqrt(H[1,1]),
POL=as(polygone, "gpc.poly"))
OMEGA=do.call("c",OMEGA)}

Note that it is possible to parallelize if there are a lot of observations,

if(n>=500)
{cl <- makeCluster(4,type="SOCK")
worker.init <- function(packages)
{for(p in packages){library(p, character.only=T)}
NULL}
clusterCall(cl, worker.init, c("gpclib","sp"))
clusterExport(cl,c("cercle","poids"))

Then, we can use standard bivariate kernel smoothing functions, but with the weights we just calculated, using a simple technique that can be related to one suggested in Ripley (1977),

fhat=kde(U,H,w=1/OMEGA,xmin=c(min(polygone[,1]),
min(polygone[,2])),xmax=c(max(polygone[,1]),
max(polygone[,2])))
fhat$estimate=fhat$estimate*sum(1/OMEGA)/n
vx=unlist(fhat$eval.points[1])
vy=unlist(fhat$eval.points[2])
VX = cbind(rep(vx,each=length(vy)))
VY = cbind(rep(vy,length(vx)))
VXY=cbind(VX,VY)
Ind=matrix(point.in.polygon(VX,VY, polygone[,1],
polygone[,2]),length(vy),length(vx))
f0=fhat
f0$estimate[t(Ind)==0]=NA
return(list(
X=fhat$eval.points[[1]],
Y=fhat$eval.points[[2]],
Z=fhat$estimate,
ZNA=f0$estimate,
H=fhat$H,
W=fhat$W))}
lissage_without_c = function(U,polygone,optimal=TRUE,h=.1)
{n=nrow(U)
IND=which(is.na(U[,1])==FALSE)
U=U[IND,]
if(optimal==TRUE) {H=Hpi(U,binned=FALSE);
H=matrix(c(sqrt(H[1,1]*H[2,2]),0,0,sqrt(H[1,1]*H[2,2])),2,2)}
if(optimal==FALSE){H= matrix(c(h,0,0,h),2,2)}
fhat=kde(U,H,xmin=c(min(polygone[,1]),
min(polygone[,2])),xmax=c(max(polygone[,1]),
max(polygone[,2])))
vx=unlist(fhat$eval.points[1])
vy=unlist(fhat$eval.points[2])
VX = cbind(rep(vx,each=length(vy)))
VY = cbind(rep(vy,length(vx)))
VXY=cbind(VX,VY)
Ind=matrix(point.in.polygon(VX,VY, polygone[,1],
polygone[,2]),length(vy),length(vx))
f0=fhat
f0$estimate[t(Ind)==0]=NA
return(list(
X=fhat$eval.points[[1]],
Y=fhat$eval.points[[2]],
Z=fhat$estimate,
ZNA=f0$estimate,
H=fhat$H,
W=fhat$W))}

So, now we can play with those functions,

base_cara_FINISTERE=base_cara[which(point.in.polygon(
base_cara$long,base_cara$lat,monFINISTERE[,1],
monFINISTERE[,2])==1),]
coord=cbind(as.numeric(base_cara_FINISTERE$long),
as.numeric(base_cara_FINISTERE$lat))
nrow(coord)
map(francemap)
lissage_FIN_withoutc=lissage_without_c(coord,
monFINISTERE,optimal=TRUE)
lissage_FIN=lissage(coord,monFINISTERE,
optimal=TRUE)
lesBreaks_sans_pop=range(c(
range(lissage_FIN_withoutc$Z),
range(lissage_FIN$Z)))
lesBreaks_sans_pop=seq(min(lesBreaks_sans_pop)*.95,
max(lesBreaks_sans_pop)*1.05,length=21)

plot_article=function(lissage,breaks,
polygone,coord){
par(mar=c(3,1,3,1))
image.plot(lissage$X,lissage$Y,(lissage$ZNA),
xlim=range(polygone[,1]),ylim=range(polygone[,2]),
breaks=breaks, col=rev(heat.colors(20)),xlab="",
ylab="",xaxt="n",yaxt="n",bty="n",zlim=range(breaks),
horizontal=TRUE)
contour(lissage$X,lissage$Y,lissage$ZNA,add=TRUE,
col="grey")
points(coord[,1],coord[,2],pch=19,cex=.1,
col="dodger blue")
polygon(polygone,lwd=2,)}

plot_article(lissage_FIN_withoutc,breaks=
lesBreaks_sans_pop,polygone=monFINISTERE,
coord=coord)

plot_article(lissage_FIN,breaks=
lesBreaks_sans_pop,polygone=monFINISTERE,
coord=coord)

If we look at the graphs, we have the following densities of car accident, with a standard kernel on the left, and our proposal on the right (with local weight adjustment when the estimation is done next to the border of the region of interest),

Similarly, in Morbihan,

With those modified kernels, hot spots appear much more clearly. For more details, the paper is online on hal.

French dataset: population and GPS coordinates

A short post today based on recent work by @3wen (Ewen Gallic, graduate Student in Rennes, spending a year in Montreal). Since we were working on a detailed French dataset (per commune), we needed a dataset containing a list allcommunes, with population and location. GPS coordinates were extracted from Google, using the following php file, inspired by http://www.andrew-kirkpatrick.com/ on Google geocoding api with php webpage. Population was interpolated from INSEE’s datasets, i.e. http://www.insee.fr/ (since data are over a 35 year period, from 1975 to 2010, changes have been taken into account as carefully are possible – e.g. merges and splits of cities – based on thatdescription). A spline model has been used for all cities (with three degrees of freedom, and null and negative interpolation became one, since we’ll be using loglinear models afterwards). Names are from that dataset, still on INSEE’s website, http://www.insee.fr/.

A zipped file can be downloaded here popfr19752010.zip, but it is also possible to use the code below (it is a 24Mo dataset). Since it was hard to find such a dataset online (different files can be found, but we found none with population and location), we have decided to upload that dataset. Please let us know if there are problems with those data…

> base=read.csv(
+ "http://freakonometrics.free.fr/popfr19752010.csv",
+ header=TRUE)

Using that code, it is possible to locate all the communes in France (metropolitan), for instance

> library(maps)
> map("france")
> points(base$long,base$lat,cex=.1,col="red",pch=19)
> points(base$long,base$lat,cex=2*base$pop_2010/
+ max(base$pop_2010),col="blue",pch=19)

Several additional lines of code on that dataset (and also others) will be uploaded, soon.

Cette oeuvre est mise à disposition sous licence Paternité – Partage à l’Identique 3.0 non transposé. Pour voir une copie de cette licence, visitez http://creativecommons.org/. Date : 24 mai 2012, par Ewen GALLIC. Sources : INSEE, API Google Maps v3 et GeoHack (coordonnées GPS), propres calculs (estimation de population à partir des données INSEE).

  • reg : code region INSEE (character)
  • dep : code departement INSEE (character, corse 201 et 202 au lieu de 2A et 2B)
  • com : code commune INSEE (character)
  • article : article du nom de la commune (character)
  • com_nom : nom de la commune (character)
  • long : longitude (numeric)
  • lat : latitude (numeric)
  • pop_i : estimation de la population à la date i (ramenée à 1 si <=0), i=1975,…,2010 (numeric)

Playing with fire (or water)

A few days ago,http://www.futilitycloset.com/published a short post based on the fourth problem of the 1987 Canadian Mathematical Olympiad (from on a problem from the 6th All Soviet Union Mathematical Competition in Voronezh, 1966). The problem is simple (as always). It is about water pistol duels (with an odd number of players)

The answer is nice, an can be read on the blog.

What puzzled me in this problem is the following: if we know, for sure, that at least one player won’t get wet, we don’t know exactly how many of them won’t get wet (assuming that if they shoot at the closest, they hit him for sure) ? It is simple to run simulations, e.g. assuming that players are uniformly distributed over a square,

NOTWET=function(n){
x=runif(n)
y=runif(n)
(d=as.matrix(dist(cbind(x,y), method = "euclidean",upper=TRUE)))
diag(d)=999999
dmin=apply(d,2,which.min)
notwet=n-length(table(dmin))
return(notwet)}

It is then rather simple to get the distribution of the number of player that did not get wet,

N25=Vectorize(NOTWET)(n=rep(25,NSim))
T=table(N25)
plot(as.numeric(names(T)),T/NSim,type="b")

The graph for different values for the total number of players is the following (based on 25,000 simulations)

If we investigate further, say with 51 players, we have a distribution for the total number of players that did not get wet which looks exactly like the Gaussian distribution,

NSim=25000
N51=Vectorize(NOTWET)(n=rep(51,NSim))
T=table(N51)
plot(as.numeric(names(T)),T/NSim,type="b",col="blue")
u=seq(0,51,by=.1)
lines(u,dnorm(u,mean(N51),sd(N51)),col="red",lty=2)

If anyone has an intuition (not to say a proof) for that, I’d be glad to hear it…

Maps with R, and polygon boundaries

With R, it is extremely easy to draw maps. Let us start with something simple, like French regions. Baptiste mentioned on his blog that shapefiles can be downloaded from http://ign.fr/ website. Hence, if you extract the zip file, it is possible to get claims frequency per region (as done in the course ACT2040),

> library(maptools)
> library(maps)
> departements<-readShapeSpatial("DEPARTEMENT.SHP")
> region<-tapply(baseFREQ[,"nbre"],
+ as.factor(baseFREQ[,"region"]),sum)/
+ tapply(baseFREQ[,"exposition"],
+ as.factor(baseFREQ[,"region"]),sum)
> depFREQ=rep(NA,nrow(departements))
> names(depFREQ)=as.character(
+ departements$CODE_REG)
> for(nom in names(region)){
+ depFREQ[names(depFREQ)==nom] =
+ region[nom]}
> plot(departements,col=gray((depFREQ-.05)*20))
> legend(166963,6561753,legend=seq(1,0,by=-.1)/20+.05,
+ fill=gray(seq(1,0,by=-.1)),cex=1.25, bty="n")

Another application is on earthquakes. It is possible to use shapefiles of tectonic plates contour, and to relate earthquakes to plates. Shapefiles can be found onhttp://www.colorado.edu/ (here).

http://freakonometrics.blog.free.fr/public/perso4/plate-tekto.gif

First, we can extract the shapes of the tectonic plates

> plates = readShapePoly("plates.shp",
+ proj4string=CRS("+proj=longlat"))
> PP=SpatialPolygons2PolySet(plates)

Consider Montreal,

> montreal=c(-73.600,45.500)

Given that specific location, it is possible to use the following code to get the associated plate,

> PLATE.loc=function(pt){
+ K=NA
+ for(k in 1:17){
+ c=point.in.polygon(pt[1], pt[2],
+ PP[PP$PID==k,c("X")],PP[PP$PID==k,c("Y")],
+ mode.checked=FALSE)
+ if(c>0){K=k}
+ }
+ return(K)}
> abline(v=montreal[1],col="red")
> abline(h=montreal[2],col="red")
> PLATE.loc(montreal)
[1] 1

and then to plot the associated tectonic plate very easily

> PLATE=function(k0){
+ library(maps)
+ map("world")
+ polygon(PP[PP$PID==k0,c("X")],PP[PP$PID==k0,c("Y")],
+ col="red")
+ for(k in (1:17)[-k0]){polygon(PP[PP$PID==k,c("X")],
+ PP[PP$PID==k,c("Y")],col="light blue")}
+ map("world",add=TRUE)}
> PLATE(PLATE.loc(montreal))

Those code were used in the paper written with Mathieu, and that will be presented on January 30th at the Geotop seminar.

Une région géographique n’est pas une variable continue

En relisant les devoirs maisons, je me suis rendu compte que certains avaient tenté de regrouper les régions (géographiques) par régions homogènes. Sauf que les régions étaient codées par un numéro (selon la codification officielle). Par exemple, dans une des bases, nous avions des assurés dans 4 zones géographiques, à savoir la région 82 (région Rhône-Alpes en rouge) la région 54 (région Poitou-Charentes en vert) la région 73 (région Midi-Pyrénées en bleu) et enfin la région 41 (région Lorraine en mauve).

> unique(baseFREQ$region) 
[1] 82 54 73 41

Une idée intéressante pour regrouper les régions pouvait être d’utiliser les arbres. Les régions étant des couleurs (on le voit bien sur la carte) et pas des variables quantitatives, il est normal de travailler sur des facteurs. D’ailleurs le code pour faire la carte est le suivant,

> library(maps) 
>  france<-map(database="france") 
>  dpt=c("Ain","Ardeche","Drome","Isere","Loire ","Rhone",  
+ "Savoie","Haute-Savoie","Charente","Charente-Maritime", 
+ "Deux-Sevres","Vienne","Ariege","Aveyron","Haute-Garonne",  
+ "Gers","Lot","Hautes-Pyrenees","Tarn","Tarn-et-Garonne", 
+ "Meurthe-et-Moselle","Meuse","Moselle","Vosges") 
>  couleur=c(rep(2,8),rep(3,4),rep(4,8),rep(6,4))  
>  match=match.map(france,dpt) 
>  color=couleur[match] 
>  map(database="france", fill=TRUE, col=color)

L’arbre sur les régions en tant que facteurs donne le découpage suivant

>  baseFREQ$fregion=as.factor(baseFREQ$region) 
>  ARBRE1=tree(nombre~fregion,data=baseFREQ,split="gini")  
>  plot(ARBRE1) 
>  text(ARBRE1)

Bon, R a la mauvaise idée de recoder les classes (mais il garde l’ordre, i.e. a correspond à la région 41, b à 54, c à 73 et d à 82). Visuellement, on retient qu’il est possible de considérer deux grandes régions, AC (i.e. 41 et 73) et BD (i.e. 54 et 72). L’intérêt des arbres sur des variables qualitatives, des facteurs, c’est que tous les regroupements sont possibles. En revanche, si on fait un arbre sur la région qui est lue en tant que nombre (quantitatif), on obtient

>  ARBRE2=tree(nombre~region,data=baseFREQ,split="gini") 
>  plot(ARBRE2) 
>  text(ARBRE2)

Il est alors impossible de regrouper dans une même classe deux régions séparées par un nombre, i.e. on ne peut regrouper 41 et 82 dans la même classe. R suggère de distinguer peut être trois régions, à savoir 82 (à droite), puis 73 (au centre) et enfin de mettre éventuellement 41 et 54 ensemble. Ce qui n’est pas la stratégie optimale quand on regroupe des facteurs.

Dis moi comment tu t’appelles…

A chaque fois que je prends le métro à Paris, ou que je traverse la gare Montparnasse en descendant du train, je suis surpris par les affiches de publicités qui nous harcèlent. Mais ce qui me surprendra le plus, ce sont les compagnes où on nous présente un Stéphane, un Mathieu, voire un Jean-Pierre, qui doivent nous rendre les affiches plus “humaines” (enfin, je pense). Les gens qui font du marketing ne sont pas forcément stupides (même si personnellement, me montrer la gentille bouille de Stéphane ne va pas me faire aimer quelqu’un qui verbalise des usagers à deux-roues), et je me suis demandé si le choix du prénom (oui, je suis persuadé que le brave Stéphane de l’affiche ne s’appelle pas Stéphane… depuis quand le monde de la pub’ coïnciderait avec le monde réel ? même mes enfants savent que ce n’est pas le cas) était vraiment pertinent…. Bref, j’ai demandé au réseau Quetelet un accès à la base des prénoms…

Pour mon nouvel ami Stéphane, on note effectivement que les Stéphane ont, en moyenne, un peu plus de 36 ans (comme sur la photo), et si on regarde les quantiles, la moitié ont entre 32 ans (quantile à 25%) et 41 ans (quantile à 75%). Bref, on cible bien les trentennaires.

 

Maintenant, on  peut se demander un peu plus précisément qui est ciblé par cette campagne, on peut regarder d’où viennent les Stéphane,

On retrouve bien un flash pendant une petite dizaine d’années, plutôt dans la partie nord de la France… Bref, la campagne semble très ciblée….

Pour les Mathieu, on veut nous vendre des téléphones portables, je crois…Et effectivement, on est sur un âge moyen beaucoup plus faible, de l’ordre de 25 ans. En revanche, un peu plus étalé que tout à l’heure, avec la moitié ayant entre 15 et 28 ans.


Ici la distribution est plus étalée sur la France, dans le temps,

Bon, le dernier, c’est Jean-Pierre…. alors Jean-Pierre, on ne peut pas trop le voir. On voit qu’il captive des enfants1 en leur racontant son expérience, ce qui correspondrait à une personne un peu âgée, avec de l’expérience. Alors effectivement, si on regarde dans l’historique, les Jean-Pierre ont 55 ans en moyenne, avec 50% entre 50 et 62 ans.


Et là aussi, on retrouve une population proche des Mathieu quant à sa répartition géographique (des pays de cheminots ?).

1 au delà de la fascination exercée par notre ami Jean-Pierre sur ces élèves, pour ma part, ce sont les élèves qui me fascinent… une classe presque exclusivement constituée de gauchers ! peut être pour accentuer le fait qu’il ne faut pas croire que la pub’, c’est comme dans la vraie vie, je ne sais pas… Remarquez, je n’ai rien contre les gauchers sur les photos ! mais c’est un peu comme les auvergnats: un ça va, mais c’est quand il y en a trop que ça devient gênant….