Maudites Abeilles

Pendant les “vacances”, j’étais avec les enfants (le plus grand et la plus petite) au palais de la découverte, à Paris, à l’exposition temporaire “faîtes vos jeux“. Comme toujours, je suis surpris quand je compare les musées de sciences, et les musées d’art (j’avais fait un billet voilà quelques années en comparant le MoMaths et le Gugenheim, à New York City). Autant ma grande fille est autant intéressée par les tableaux que par les sciences, autant les deux autres préfèrent mille fois les seconds au premier. On a d’ailleurs passé plus de 4 heures au musée (dont une à faire toutes les animations de “faîtes vos jeux“), il a fallu y aller parce que j’avais une réunion de boulot (oui, les parisiens n’étaient pas en vacances, et j’avais une réunion impossible à déplacer). D’où les guillemets à “vacances” (le second jour de notre excursion, c’était une intervention de deux heures en fin d’après midi).

Mais je m’égare… A un des stands de l’animation, il avait le défi des abeilles,

Sur notre “damier” 8×8 (ici), on a 32 cases bleues et 32 rouges (ou oranges), et on cherche la probabilité qu’un chemin permette de passer d’un côté à l’autre. Sur le moment, le problème m’a pas mal pris la tête, et je me disais que la réponse devait être compliquée… en tout cas, que ça n’était pas aussi simple que ça…

En ligne, on peut trouver le guide de l’enseignant qui apporte des compléments aux divers ateliers, et qui propose, en l’occurrence, une solution

La réponse serait-ici 1/2. Mais la réponse ne me satisfait pas, en particulier la phrase “son absence à l’inverse impose un chemin orange” qui suppose qu’il existe forcément un chemin, soit rouge, soit bleu, pour passer d’un côté à l’autre, ce qui ne me semble pas trivial ! Quant

Par exemple sur le damier 4×4 ci-dessous,

il n’existe aucun chemin rouge allant  du bas vers le haut, et aucun chemin bleu non plus ! Et si on regarde, il existe (au moins) un chemin rouge allant de la gauche vers la droite, et un chemin bleu. Bref, la réponse 1/2 me semble douteuse….  Pour ceux qui veulent de vrai hasard (mais qui a autant de chance d’apparaître que celui au dessus), on a le graph ci-dessous

Mais maintenant essayons de la calculer cette probabilité. Par simulations au moins. Regardons un damier,

on peut le visualiser sous forme d’une matrice

ou en décalant un peu d’une matrice carrée 4×4 (ici)

Si on regarde ce que sont les voisins (afin de faire un chemin), sur le schéma initial, on cherchait les configuration suivantes

ce qui, dans notre matrice, correspond aux éléments ci-dessous,

Bref, on va pouvoir connecter des points d’un graph, et regarder si les 4 premiers éléments du graphs sont connectés aux 4 derniers

(le but est désormais de savoir si le nord et le sud connectent, peu importe la couleur). On peut commencer par créer le graph

d=4
set.seed(1)
C=sample(c(rep(1,d^2/2),rep(2,d^2/2)),size=d^2)
M=matrix(C,d,d)
> M
[,1] [,2] [,3] [,4]
[1,] 1 1 2 2
[2,] 1 2 1 1
[3,] 2 2 1 1
[4,] 2 2 1 2

ce qui correspond bien à notre matrice ci-dessous

On va alors créer le graph, c’est à dire chercher les noeuds qui communiquent

library(igraph)
B=NULL
for(i in 1:(d^2)){
x=rep(i,6)
y=c(i-d,i-d+1,i-1,i+1,i+d-1,i+d)
if(i%%d==1) y=c(i-d,i-d+1,NA,i+1,NA,i+d)
if(i%%d==0) y=c(i-d,NA,i-1,NA,i+d-1,i+d)
D=data.frame(x,y)
D=D[(D$y>=1)&(D$y<=d^2),] B=rbind(B,D[which(C[D$y]==C[D$x]),]) } B=B[(B[,2]>=1)&(B[,2]<=d^2),]
G=as.vector(t(B))
m=make_graph(G)

On peut ainsi voir ici que 1 communique avec 2, et 5

> B[B$x==1,]
x y
4 1 2
6 1 5

On peut d’ailleurs visualiser le graph, et voir qu’il est constitué de 4 graphs disjoints,

> plot(m)

(le coin 16 étant isolé, sans connection avec le reste – étrangement – il n’apparait pas sur le graph)

Pour trouver tous les points qui communiquent entre eux, on peut itérer un peu (voire beaucoup, peu importe)

G_iter=connect(make_graph(G), d^2)
connectpt=function(x) as.numeric(adjacent_vertices(G_iter,x)[[1]])

Ici, le nord et le sud communiquent s’il y a au moins un chemin (peu importe la longueur)

sum(unlist(Vectorize(connectpt)(1:d))%in%(d^2+1-(1:d)))>0

On peut alors une petite fonction, qui va générer une répartition des couleurs sur un damier carré (de longueur paire)

simu=function(d){
C=sample(c(rep(1,d^2/2),rep(2,d^2/2)),size=d^2)
B=NULL
for(i in 1:(d^2)){
x=rep(i,6)
y=c(i-d,i-d+1,i-1,i+1,i+d-1,i+d)
if(i%%d==1) y=c(i-d,i-d+1,NA,i+1,NA,i+d)
if(i%%d==0) y=c(i-d,NA,i-1,NA,i+d-1,i+d)
D=data.frame(x,y)
D=D[(D$y>=1)&(D$y<=d^2),] B=rbind(B,D[which(C[D$y]==C[D$x]),]) } B=B[(B[,2]>=1)&(B[,2]<=d^2),] G=as.vector(t(B)) G_iter=connect(make_graph(G), d^2) connectpt=function(x) as.numeric(adjacent_vertices(G_iter,x)[[1]]) sum(unlist(Vectorize(connectpt)(1:d))%in%(d^2+1-(1:d)))>0}

On peut alors construire la fonction suivante, pour faire un peu de monte carlo

appr=function(d,nsimu) mean(Vectorize(simu)(rep(d,nsimu)))

Par exemple sur notre damier 4×4, on a environ 3 chances sur 4 d’avoir un chemin (peu importe sa couleur) qui permette de rejoindre le sud depuis le nord

> appr(4,10000)
[1] 0.7446

alors que sur le 8×8 du palais de la découverte, on a

> appr(8,10000)
[1] 0.7805

78% de chances (un peu plus). Ce qui semble confirmer qu’il n’y a pas forcément un chemin permettant d’aller du nord au sud.

On peut d’ailleurs faire un petit dessin pour noter que quand la dimension du damier augmente (je n’ai pas tenté ici le cas impair – avec une case de plus d’une certaine couleur), la probabilité de trouver au moins un chemin semblerait augmenter (même si ça reste à prouver, et en particulier, je serais curieux de voire vers quoi ça converge…)

 

Visualiser les résultats du premier tour

Juste quelques lignes de code, pour visualiser les résultats du premier tour en France. L’idée est de faire une carte assez minimaliste, avec des cercles centrés sur les centroïdes des départements. On commence par récupérer les données pour le fond de carte, un fichier 7z sur le site de l’ign.

download.file("https://wxs-telechargement.ign.fr/oikr5jryiph0iwhw36053ptm/telechargement/inspire/GEOFLA_THEME-DEPARTEMENTS_2016$GEOFLA_2-2_DEPARTEMENT_SHP_LAMB93_FXX_2016-06-28/file/GEOFLA_2-2_DEPARTEMENT_SHP_LAMB93_FXX_2016-06-28.7z",destfile = "dpt.7z")

On a dans ce fichier des informations sur les centroïdes

library(maptools)
library(maps)
departements<-readShapeSpatial("DEPARTEMENT.SHP")
plot(departements)
points(departements@data$X_CENTROID,departements@data$Y_CENTROID,pch=19,col="red")

Comme ça ne marche pas très bien, on va le refaire à la main, par exemple pour l’Ille-et-Vilaine,

pos=which(departements@data[,"CODE_DEPT"]==35)
Poly_35=departements[pos,]
plot(departements)
plot(Poly_35,col="yellow",add=TRUE)
departements@data[pos,c("X_CENTROID","Y_CENTROID")]
points(departements@data[pos,c("X_CENTROID","Y_CENTROID")],pch=19,col="red")
library(rgeos)
(ctd=gCentroid(Poly_35,byid=TRUE))
points(ctd,pch=19,col="blue")

Comme ça marche mieux, on va utiliser ces centroïdes.

ctd=as.data.frame(gCentroid(departements,byid=TRUE))
plot(departements)
points(ctd,pch=19,col="blue")

Maintenant, il nous faut les résultats des élections, par département. On peut aller scraper le site du ministère d’ intérieur. On a une page par département, donc c’est facile de parcourir. Par contre, dans l’adresse url, il faut le code région. Comme je suis un peu fainéant, au lieu de faire une table de correspondance, on teste tout, jusqu’à ce que ça marche. L’idée est d’aller cherche le nombre de voix obtenues par un des candidats.

candidat="M. Emmanuel MACRON"
library(XML)
voix=function(no=35){
testurl=FALSE
i=1
vect_reg=c("084","027","053","024","094","044","032","028","075","076","052",
"093","011","001","002","003","004")
region=NA
while((testurl==FALSE)&(i<=20)){
reg=vect_reg[i]
nodpt=paste("0",no,sep="")
# if(!is.na(as.numeric(no))){if(as.numeric(no)<10) nodpt=paste("00",no,sep="")}
url=paste("http://elections.interieur.gouv.fr/presidentielle-2017/",reg,"/",nodpt,"/index.html",sep="")
test=try(htmlParse(url),silent =TRUE)
if(!inherits(test, "try-error")){testurl=TRUE
region=reg}
i=i+1
}
tabs tab=tabs[[2]]
nb=tab[tab[,1]==candidat,"Voix"]
a<-unlist(strsplit(as.character(nb)," "))
as.numeric(paste(a,collapse=""))}

On peut alors tester

> voix(35)
[1] 84648

Comme ça semble marcher, on le fait pour tous les départements

liste_dpt=departements@data$CODE_DEPT
nbvoix=Vectorize(voix)(liste_dpt)

On peut alors visualiser sur une carte.

plot(departements,border=NA)
points(ctd,pch=19,col=rgb(1,0,0,.25),cex=nbvoix/50000)

Et on peut aussi tenter pour une autre candidate,

candidat="Mme Marine LE PEN"

et on obtient la carte suivante

plot(departements,border=NA)
points(ctd,pch=19,col=rgb(0,0,1,.25),cex=nbvoix/50000)