Freakonometrics

To content | To menu | To search

Wednesday, October 31 2012

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...

Tuesday, October 30 2012

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 à @, 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... la conjecture 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 !

Thursday, August 23 2012

Est-ce vraiment trop injuste ?

Depuis le début de la semaine, après avoir déposé les enfants au camp de jour du Musée des Beaux Arts, je viens à pieds à l'université. Chaque fois, je me dis que je pourrais prendre le bus, mais comme aucun ne vient, je commence en marchant, en me disant que si un bus passe, je le prendrais. Et tous les matins, j'arrive au bureau sans m’être fait dépassé par le moindre bus. Et bien sur, j'en ai croisé un paquet qui passaient en sens inverse...

C'est vraiment trop injuste ? En tant que statisticien, je dirais que non, tout simplement à cause du biais de sélection. Si on veut formaliser un peu, on va oublier l'aléa, et raisonner avec des bus qui passent de manière purement déterministe. On va aussi supposer qu'il y en a autant dans un sens que dans l'autre...

J'ai une distance http://freakonometrics.blog.free.fr/public/perso6/bus-01.gif à parcourir, on suppose que les bus sont espacés d'une distance http://freakonometrics.blog.free.fr/public/perso6/bus-02.gif, et qu'ils avancent à une vitesse http://freakonometrics.blog.free.fr/public/perso6/bus-03.gif. Autrement dit, les bus passent tous les http://freakonometrics.blog.free.fr/public/perso6/bus-04.gif secondes (si ma vitesse est exprimée en secondes).

Moi, j'avance à une vitesse http://freakonometrics.blog.free.fr/public/perso6/bus-07.gif avec http://freakonometrics.blog.free.fr/public/perso6/bus-06.gif (oui, on va supposer que je vais moins vite que le bus... ce qui n'est pas forcément une hypothèse faible aux heures de pointes, mais disons que le problème n'a de sens que si aller en bus me permet d'aller plus vite). Le temps que je vais mettre si je fais tous le trajet à pied est . Maintenant, comptons les bus qui passent en face. Je vais croiser tous ceux qui sont déjà sur ma portion de trajet, et il y en a http://freakonometrics.blog.free.fr/public/perso6/bus-08.gif. En plus, je croiserais tous ceux qui vont arriver à l'université, et qui n'y sont pas encore, soit http://freakonometrics.blog.free.fr/public/perso6/bus-09.gif, i.e. le temps qu'il me reste à marcher divisé par le temps qui s'écoule entre deux bus. On a alors un total de http://freakonometrics.blog.free.fr/public/perso6/bus-10.gif bus à croiser, en face.

De mon coté cette fois. Je peux compter tous ceux qui vont arriver à l'université, sauf que cette fois, je devrais enlever tous ceux qui sont déjà sur le trajet, et que je ne croiserais pas (car je vais moins vite qu'eux, par hypothèse). Le nombre de bus qui vont me doubler sera alors cette fois http://freakonometrics.blog.free.fr/public/perso6/bus-11.gif

Si je regarde le ratio du nombre de bus croisés sur le nombre de bus qui m'aura doublé, j'obtiens

http://freakonometrics.blog.free.fr/public/perso6/bus-13.gif

Les plus malins auront noté que ce ratio est aussi le ratio entre la vitesse des bus qui me croisent divisé par la vitesse des bus qui me doublent, quand je suis le référentiel (les vitesses sont exprimées par rapport à moi).

Par exemple, si je vais deux fois mois vite que le bus, il y aura trois fois plus de bus qui passent en sens inverse. Et sept fois plus si je vais juste 25% moins vite que le bus. Ce qui n'est pas absurde, quand on pense que je fais le trajet en 20 minutes, et le bus 15...

Saturday, June 30 2012

Simple and heuristic optimization

This week, at the Rmetrics conference, there has been an interesting discussion about heuristic optimization. The starting point was simple: in complex optimization problems (here we mean with a lot of local maxima, for instance), we do not necessarily need extremely advanced algorithms that do converge extremly fast, if we cannot ensure that they reach the optimum. Converging extremely fast, with a great numerical precision to some point (that is not the point we're looking for) is useless. And some algorithms might be much slower, but at least, it is much more likely to converge to the optimum. Wherever we start from.
We have experienced that with Mathieu, while we were looking for maximum likelihood of our MINAR process: genetic algorithm have performed extremely well. The idea is extremly simple, and natural. Let us consider as a starting point the following algorithm,

  1. Start from some 
  2. At step , draw a point  in a neighborhood of  
  • either  then  
  • or  then 
This is simple (if you do not enter into details about what such a neighborhood should be). But using that kind of algorithm, you might get trapped and attracted to some local optima if the neighborhood is not large enough. An alternative to this technique is the following: it might be interesting to change a bit more, and instead of changing when we have a maximum, we change if we have almost a maximum. Namely at step ,
  • either then  
  • or  then 
for some . To illustrate the idea, consider the following function
> f=function(x,y) { r <- sqrt(x^2+y^2);
+ 1.1^(x+y)*10 * sin(r)/r }
(on some bounded support). Here, by picking noise and  values arbitrary, we have obtained the following scenarios
> x0=15
> MX=matrix(NA,501,2)
> MX[1,]=runif(2,-x0,x0)
> k=.5
> for(s in 2:501){
+  bruit=rnorm(2)
+  X=MX[s-1,]+bruit*3
+  if(X[1]>x0){X[1]=x0}
+  if(X[1]<(-x0)){X[1]=-x0}
+  if(X[2]>x0){X[2]=x0}
+  if(X[2]<(-x0)){X[2]=-x0}
+  if(f(X[1],X[2])+k>f(MX[s-1,1],
+    MX[s-1,2])){MX[s,]=X}
+  if(f(X[1],X[2])+k<=f(MX[s-1,1],
+    MX[s-1,2])){MX[s,]=MX[s-1,]}
+}

 

It does not always converge towards the optimum,

and sometimes, we just missed it after being extremely unlucky

Note that if we run 10,000 scenarios (with different random noises and starting point), in 50% scenarios, we reach the maxima. Or at least, we are next to it, on top.
What if we compare with a standard optimization routine, like Nelder-Mead, or quasi gradient ?Since we look for the maxima on a restricted domain, we can use the following function,

> g=function(x) f(x[1],x[2])
> optim(X0, g,method="L-BFGS-B",
+ lower=-c(x0,x0),upper=c(x0,x0))$par
In that case, if we run the algorithm with 10,000 random starting point, this is where we end, below on the right (while the heuristic technique is on the left),

In only 15% of the scenarios, we have been able to reach the region where the maximum is.
So here, it looks like an heuristic method works extremelly well, if do not need to reach the maxima with a great precision. Which is usually the case actually.

Monday, May 7 2012

Bayes is playing Russian roulette

There was (once again) a nice puzzle in http://www.futilitycloset.com/. Bayes and a good friend are playing Russian roulette. The revolver has six chambers. He puts two bullets in two adjacent chambers, spin the cylinder, hold the gun to his friend's head, and pull the trigger. It clicks. So it is now Bayes's turn: he can choose either to spin the cylinder again or leave it as it is. Which is better? Hopefully, Bayes knows his theorem: if he does spin it, the probability of getting killed is 2 out of 6 (four empty chambers out of six), but if he does not, since his friend is still alive, then the hammer should be next to one of the four cylinders in red, below
So here, there is 3 chance out of 4 to survive, i.e. the probability of getting killed is 1 out of 4 (while it was 1 out of 3 when spinning). So Bayes should not spin. And as always, it is possible to see it is a more general result: more generally, in a revolver with http://freakonometrics.blog.free.fr/public/perso5/bullet01.gif chambers, it there are http://freakonometrics.blog.free.fr/public/perso5/bullet02.gif bullets in http://freakonometrics.blog.free.fr/public/perso5/bullet02.gif adjacent chambers,  if the first player survives, the probability of getting killed is k over http://freakonometrics.blog.free.fr/public/perso5/bullet01.gif, when spinning, while it would be 1 over http://freakonometrics.blog.free.fr/public/perso5/bullet03.gif if we don't. Not spinning is better if and only if

http://freakonometrics.blog.free.fr/public/perso5/bullet04.gif
i.e.
http://freakonometrics.blog.free.fr/public/perso5/bullet05.gif
So you'd better not spin, unless there was one bullet in the revolver, i.e. http://freakonometrics.blog.free.fr/public/perso5/bullet06.gif... or http://freakonometrics.blog.free.fr/public/perso5/bullet07.gif (in that case, it might not be a good idea actually to play the game).

Friday, April 27 2012

Proving tautological versus trivial results in mathematics

There is something that might be fun in mathematics, which is the connexion between trivial, tautological and difficult questions. Sometimes, things are so intuitive, that they seem to be obvious. But mathematicians aren't jedis, and they should not trust too much their intuition... I mean intuition is fine, but it is not a proof. It is like those standard results we learn in topology courses, e.g. "the closure of an open ball is not necessarily the closed ball". The other thing is that after a while, you try to prove something, until someone makes you realize that it is the definition... 
And this morning, while I was trying to make a coffee, @renaudjf came with a simple question (yes, it always starts like that). Consider the standard algorithm to generate a conditional random variable. Assume that  has a priori distribution , and that , given , has (conditional) distribution .
The standard idea is monte carlo simulation, to generate values of , is
  •  step 1: generate 
  •  step 2: given that generation of , generate 
"Can we prove that we actually generate from the (true, maybe hard to characterize) non-conditional distribution of  ? Or is it just trivial ?". After having those previous philosophical questions, we came to the point that if it was trivial, then we should be able to prove it. A standard way of writing the algorithm is to use the quantile based technique
  •   with ,
  •   with ,
For instance, to generate negative binomial distribution
n=1
theta=rgama(n,3,3)
X=rpois(n,lambda=theta)
Thus, let  where  and  are two independent random variables with a uniform distribution on the unit interval. Let us try to derive its distribution, i.e.
so
if we consider the following change of variate 
which is exactly the non-conditional distribution of .
And then, you're quite happy because you've been able to prove a trivial result ! But next time, I promise, we'll try to derive an amazing theorem that will change humanity... but next time only, first, let us prove trivial results.

Monday, April 16 2012

Pigeonholes and triangles

Once again, there was a nice maths puzzle on http://www.futilitycloset.com/ last week (but without further reference). The question was the following, "Five points are located in an equilateral triangle with 10-inch sides (or on its perimeter). What’s the maximum distance between the two closest points?" Actually, this is simply an application of Dirichlet's pigeonholes theorem, as mentioned in the answer of the puzzle, "Connect the midpoints of the triangle’s sides to make four smaller triangles. Because there are five points, two of them must fall within one of these triangles. The maximum distance between these two is 5 inches."

Thus, with Dirichlet's pigeonholes theorem, we know not only the maximal minimum distance, but also where points must be (on corners of inside triangles). Here, there might be two possibilities (with also the different shapes obtained using rotations),

Further, we also observe that this result is valid not only with five points, but with six. And if we go further, e.g. with nine points, we have the following

So actually, it is possible to have a simple conjecture: let denote the number of points, and let be so that

then the minimal distance is . Which can be visualized below,

Based on that pigeonhole theorem, I have the intuition that this result is valid (here we just count the number of inside-triangles), but can we check if this is correct or not ? One idea might be to draw points randomly, and so see where points might end, or at least the maximal distance obtained over millions of random draws... But standard monte carlo might take a while... so we can use two ideas. One idea is from quasi-monte carlo techniques: since we want points to be to be as separate as possible, we do not need to draw randomly in the triangle, but perhaps we can draw randomly points on some grid. A second idea from the latin hypercube technique (and that pigeonholes theorem): instead to generating points randomly in the triangle, perhaps we can draw them in specific regions. For instance, with five points, we know that 4 points have to be in those sub-triangles, and the additional point in any one of those triangle. And because of symmetry, with five points, we can claim that this additional point has to be specifically in one sub-triangle. A random sample with five points within the same sub-triangle will be useless (and a waste of computational time).

With the following code, we define a grid, for a triangle, either upward, or downward, starting from some point (on the left), with a given length, and a given number of subdivision.

TRIANGLES=function(xinf,yinf,l,n,updown="up"){
X=NA;Y=NA
for(i in n:1){
u=xinf+seq(0+(n-i)/2/(n-1),1-(n-i)/2/(n-1),length=i)*l
if(updown=="up") v=rep(yinf+(n-i)*sqrt(3)/2*l/(n-1),i)
if(updown=="down") v=rep(yinf-(n-i)*sqrt(3)/2*l/(n-1),i)
X=c(X,u);Y=c(Y,v)}
return(cbind(X[-1],Y[-1]))}

Here are grid with respectively 20 and 50 points on the lower side. It is then possible to define 4 grids, corresponding to the four sub-triangles,

k=3;st=6
firstgrid=TRIANGLES(0,0,(k-2)/(k-1),k-1)
secondgrid1=TRIANGLES(
firstgrid[1,1],firstgrid[1,2],1*(k-2)/(k-1),st)
secondgrid2=TRIANGLES(
firstgrid[2,1],firstgrid[1,2],1*(k-2)/(k-1),st)
secondgrid3=TRIANGLES(
firstgrid[3,1],firstgrid[3,2],1*(k-2)/(k-1),st)
secondgrid4=TRIANGLES(
firstgrid[3,1],firstgrid[3,2],1*(k-2)/(k-1),st,updown="down")

Then, we just draw randomly five points on that grid, in the four sub-triangles,

N=5
Dmax=0
setpointmax=matrix(0,4,2)
indice=c(1:4,sample(1:4,size=N-4,replace=FALSE))
tindice=table(indice)
indice1=sample(1:nrow(secondgrid1),size=
tindice[1],replace=FALSE)
indice2=sample(1:nrow(secondgrid2),size=
tindice[2],replace=FALSE)
indice3=sample(1:nrow(secondgrid3),size=
tindice[3],replace=FALSE)
indice4=sample(1:nrow(secondgrid4),size=
tindice[4],replace=FALSE)
setpoint=rbind(secondgrid1[indice1,],secondgrid2[indice2,],
secondgrid3[indice3,],secondgrid4[indice4,])

No, we can run a code, where we keep in mind locations of the five points each time we beak a record,

D=min(dist(setpoint,"euclidean"))
if(D>Dmax){Dmax=D
setpointmax=setpoint}

Here are some locations obtained after running the algorithm a few times (with five points)

On the graph below, we can visualize the time it takes before having a record, and the convergence towards 1/2 (which is the true value of the maximal distance)

The convergence is slow... extremely slow... However, we can run the same code for more than five points, e.g. seven points (actually, here sub-triangles are not used here, and it looks like we have been lucky here, since the convergence was rather fast),

Sunday, April 8 2012

Eating chocolate, an Easter problem

Assume that there are (say) 100 chocolate eggs in a basket, 20 are dark chocolate, while 80 are milk chocolate. Unfortunately, eggs are wrapped, and there is no way you can distinguish them. My daughter has the following algorithm for eating them (and she actually plans to eat all of them)

  1. if there are eggs in her basket, she picks one - at random - looks if it is either dark or milk chocolate, write it down on a piece of paper (just to remember how many of each kind are left), eat it, and move to strategy 2.
  2. if there are eggs in her basket, she picks one - at random - looks if it is either dark or milk chocolate, write it down on a piece of paper and:
  • if it is the same kind as the one she got before, then eat it, and go again to step 2.
  • if it is not the same kind as the one she got before, she wraps it back, and go again to step 1.

At the end, if there is only one egg left, the probability that it is a milk chocolate egg is exactly 1/2... Nice, isn't it ?

Continue reading...

Friday, April 6 2012

Nonconvexity, and playing indoor paintball

Following the two previous posts (here and there), on the number of people that don't get wet while playing with water pistols, consider now an indoor version, in a non-convex room (i.e. player behind wall are now, somehow, protected). In the previous posts, players where playing on a square field, and I briefly mentioned that if the field was a disk, results would have been (roughly) the same: so far, the shape of the field was not an issue. But what if the field is no longer convex,

library(sp)
plot(0:2,0:2,col="white",xlab="",ylab="")
MAP=Polygon(cbind(c(0,0,1,1,2,2,0),
c(0,2,2,1,1,0,0)))
polygon(MAP@coords,col="light blue")

and players hidden behind the wall cannot be reached (red lines above are impossible hits). As earlier, it is still possible to look at the closest neighbor, we just have to exclude pairs that can no longer hit each other.

And again, it is possible to plot safe zones in green.

Once again, it is possible to look more closely are those supposed-to-be "safe zones", i.e. by looking at the distribution of the location of players that were dry at the end of the game. With 11 players, we obtain


What about the distribution of the number of dry players, over a game ?
touch=function(x1,y1,x2,y2,n=251){
X=seq(x1,x2,length=n)
Y=seq(y1,y2,length=n)
sum(point.in.polygon(X,Y,MAP@coords[,1],
MAP@coords[,2], mode.checked=FALSE)==0)==0
}
 
NOTWETnc=function(n,p){
sx=runif(50)*2;sy=runif(50)*2
IN=which(point.in.polygon(sx,sy,MAP@coords[,1],
MAP@coords[,2], mode.checked=FALSE)==1)
Sx=sx[IN];Sy=sy[IN]
Sx=Sx[1:n];Sy=Sy[1:n]
IN=IN[1:n]
MI=matrix(NA,n,n)
for(i in 1:(n-1)){
for(j in (i+1):(n)){
MI[j,i]=MI[i,j]=touch(Sx[i],Sy[i],Sx[j],Sy[j])
}}
(d=as.matrix(dist(cbind(Sx,Sy),
method = "euclidean",upper=TRUE)))
diag(d)=999999
dpossible=d
dpossible[MI==FALSE]=999999
dmin=apply(dpossible,2,which.min)
#whonotwet=( (1:n) %notin% names(table(dmin)) )
notwet=n-length(table(dmin))
return(notwet)}
 
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)}
 
NSim=10000
Nnc=Vectorize(NOTWETnc)(n=rep(11,NSim))
Nc=Vectorize(NOTWET)(n=rep(11,NSim))
T=table(Nc)
Tn=table(Nnc)
plot(as.numeric(names(Tn)),
Tn/NSim,type="b",col="blue")
lines(as.numeric(names(T)),
T/NSim,type="b",col="red",pch=4)

On 11 players, we have the same distribution as the one on a square field. So convexity is not a key issue here...

Strange isn't it. And with an odd number of player, not only there is at least one dry player, but at least, half of the players (maybe minus one) have to be wet...

Monday, April 2 2012

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...

Tuesday, March 20 2012

Probabilité et contre-intuition

Je suis tombé un peu par hasard sur un vieux papier de Peter Enis, Enis (1973), qui revient sur une relation que l'on prend toujours pour acquis (moi le premier) à savoir http://freakonometrics.blog.free.fr/public/perso5/proba07.gif. Mais il faut des hypothèses complémentaires pour que cette relation soit vraie. Le fait que http://freakonometrics.blog.free.fr/public/perso5/proba09.gif existe par exemple ne suffit pas ! (je laisse les amateurs de résultats contre-intuitifs le temps de réfléchir à cette affirmation)

Continue reading...

Saturday, February 4 2012

Good old days, politics and mathematics

2012 will be, in several countries, a presidential election year. Some decades ago, candidates were not always supposed to spend months on some demagogic democratic debates, but had time to spend on more important problems. Like mathematics... For instance, a congressman, who became the 20th president of the United States, James Garfield, gave the following proof to the Pythagorean theorem (actually, he wrote that proof five years before he become President). The legend claims that he found this proof in 1876 during a mathematics discussion with some of the members of Congress... Those were good old days, when politicians were interested in mathematics, and sciences. The proof he suggested was the following

http://freakonometrics.blog.free.fr/public/perso5/PYTH-1.gif

i.e. http://freakonometrics.blog.free.fr/public/perso5/pyth02.gif since

http://freakonometrics.blog.free.fr/public/perso5/pythagore01.gif

I found that nice graph in Roger Nelsen's book. For more details, Klebe (1995), or on wikipedia. And for those who love proofs without words, look at the 96 geometrical proofs of the Pythagorean theorem mentioned on http://www.cut-the-knot.org/.

Credit: Charles Logan, played by in Gregory Itzin, in 24

Sunday, January 22 2012

Even odds

This evening, I found a nice probabilistic puzzle on http://www.futilitycloset.com/ "A bag contains 16 billiard balls, some white and some black. You draw two balls at the same time. It is equally likely that the two will be the same color as different colors. What is the proportion of colors within the bag?"
To be honest, I did not understood the answer on the blog, but if we write it down, we want to solve
http://freakonometrics.blog.free.fr/public/perso5/futil-01.gif
Let us count: if http://freakonometrics.blog.free.fr/public/perso5/futil-04.gif is the total number of balls, and if http://freakonometrics.blog.free.fr/public/perso5/futil-03.gif is the number of white http://freakonometrics.blog.free.fr/public/perso5/futil-05.gif balls then
http://freakonometrics.blog.free.fr/public/perso5/cccccwwwwwww.gif
I.e. we want to solve a polynomial equation (of order 2) in http://freakonometrics.blog.free.fr/public/perso5/futil-07.gif, or to be more precise, in  http://freakonometrics.blog.free.fr/public/perso5/futil-00.gif
http://freakonometrics.blog.free.fr/public/perso5/futil-06.gif
If http://freakonometrics.blog.free.fr/public/perso5/futil-04.gif is equal to 16, then http://freakonometrics.blog.free.fr/public/perso5/futil-03.gif is either 6 or 10. It can be visualized below
> balls=function(n=16){
+ NB=rep(NA,n)
+ for(k in 2:(n-2)){
+ NB[k]=(k*(k-1)+(n-k)*(n-k-1))
+ }
+ k=which(NB==n*(n-1)/2)
+ if(length(k)>0){
+ plot(1:n,NB,type="b")
+ abline(h=n*(n-1)/2,col="red")
+ points((1:n)[k],NB[k],pch=19,col="red")}
+ return((1:n)[k])}
> balls()
[1] 6 10

But more generally, we can seek other http://freakonometrics.blog.free.fr/public/perso5/futil-04.gif's and other pairs of solutions of such a problem. I am not good in arithmetic, so let us run some codes. And what we get is quite nice: if http://freakonometrics.blog.free.fr/public/perso5/futil-10.gif admits a pair of solutions, then http://freakonometrics.blog.free.fr/public/perso5/futil-10.gif is the squared of another integer, say http://freakonometrics.blog.free.fr/public/perso5/futil-13.gif. Further, the difference between http://freakonometrics.blog.free.fr/public/perso5/futil-11.gif and http://freakonometrics.blog.free.fr/public/perso5/futil-12.gif is precisely http://freakonometrics.blog.free.fr/public/perso5/futil-13.gif. And http://freakonometrics.blog.free.fr/public/perso5/futil-12.gif will be one of the answers when the total number of balls will be http://freakonometrics.blog.free.fr/public/perso5/futil-20.gif. Thus, recursively, it is extremely simple to get all possible answers. Below, we have http://freakonometrics.blog.free.fr/public/perso5/futil-04.gif, http://freakonometrics.blog.free.fr/public/perso5/futil-03.gif, http://freakonometrics.blog.free.fr/public/perso5/futil-21.gif and the difference between http://freakonometrics.blog.free.fr/public/perso5/futil-21.gif and http://freakonometrics.blog.free.fr/public/perso5/futil-03.gif,
> for(s in 4:1000){
+ b=balls(s)
+ if(length(b)>0) print(c(s,b,diff(b)))
+ }
[1] 9 3 6 3
[1] 16 6 10 4
[1] 25 10 15 5
[1] 36 15 21 6
[1] 49 21 28 7
[1] 64 28 36 8
[1] 81 36 45 9
[1] 100 45 55 10
[1] 121 55 66 11
[1] 144 66 78 12
[1] 169 78 91 13
[1] 196 91 105 14
[1] 225 105 120 15
[1] 256 120 136 16
[1] 289 136 153 17
[1] 324 153 171 18
[1] 361 171 190 19
[1] 400 190 210 20
[1] 441 210 231 21
[1] 484 231 253 22
[1] 529 253 276 23
[1] 576 276 300 24
[1] 625 300 325 25
[1] 676 325 351 26
[1] 729 351 378 27
[1] 784 378 406 28
[1] 841 406 435 29
[1] 900 435 465 30
[1] 961 465 496 31
Thus, given http://freakonometrics.blog.free.fr/public/perso5/futil-22.gif, consider an urn with http://freakonometrics.blog.free.fr/public/perso5/futil-23.gif balls. We draw two balls at the same time. It is equally likely that the two will be the same color as different colors. Then the number of colors within the bag are respectively
http://freakonometrics.blog.free.fr/public/perso5/futil-24.gif
Finally, observe that the http://freakonometrics.blog.free.fr/public/perso5/futil-11.gif's are well known, from Pascal's triangle,

also known as triangular numbers,

http://freakonometrics.blog.free.fr/public/perso5/tr-pascal.gif

Maths can be magic, sometimes...

Wednesday, December 21 2011

Quel lancer faire à Serpents-Echelles ?

Hier soir, j'évoquais l'utilisation des chaînes de Markov au jeu Serpents et Échelles. Comme me le faisait remarquer Jean-Philippe, de manière assez étrange, les enfants sont toujours beaucoup plus content après avoir fait un 6 qu'après avoir fait un 1. Mais est-ce réellement la valeur optimale que le dé doit prendre ? Ça dépend bien sûr de la position. Par exemple, au premier lancer, on peut se demander ce que deviendrait le nombre de tours moyen nécessaires pour finir le jeu, conditionnellement aux 6 lancers possibles. On peut calculer, toujours conditionnellement aux 6 valeurs possibles de dés (et des positions où on va se retrouver), l'espérance du nombre de tours pour finir la partie,

esperance=function(h0){
INITIAL = as.numeric(which(M[h0+1,]>0))-1
ESPERANCE=rep(NA,length(INITIAL))
names(ESPERANCE)=INITIAL
for(k in 1:length(INITIAL)){
initial=rep(0,n+1); initial[INITIAL[k]]=1
distrib=initial%*%M
game=rep(NA,1000)
for(h in 1:length(game)){
game[h]=distrib[n+1]
distrib=distrib%*%M}
ESPERANCE[k]=sum(1-game)}
return(ESPERANCE)}

(à partir du code mis en ligne hier), e.g. pour le premier lancer,

> esperance(0)
1        2        3        5        6       14
32.16499 31.99947 31.82348 31.47954 31.36441 29.83020 

où la valeur indiquée est la position où l'on pourrait se retrouver (la dernière correspond à un 4). On note que le meilleur premier lancer possible est celui qui nous amène sur la première échelle, i.e. la 4ème case,

Si on regarde case par case, on a les "meilleurs" lancers de dés suivant (sans forcément l'unicité) en fonction de sa position sur le plateau

(avec en rouge les serpents et en bleu les échelles). Amusant non ?

Tuesday, November 22 2011

les garçons ont-ils plus de sœurs que les filles ?

Ce qui est pénible avec les paradoxes, en tant que scientifique, c'est qu'ils nous obligent à remettre en cause notre intuition, sans cesse. Car mine de rien, on fonctionne tous à l'intuition, c'est un peu comme ça qu'on monte un programme de recherche: on a l'intuition que l'on pourra obtenir des résultats, mais tout reste à prouver.

Prenons la question suivante "les garçons on-t-ils plus de sœurs que les filles ?" en moyenne. N'importe qui avec un peu de bon sens répondrait "oui", l'idée naturelle est que la fille ne peut pas se compter elle même comme sœur, donc les filles devraient en compter un peu moins... Et comme toujours la formulation est importante. On ne demande pas, par exemple "les garçons ont-t-ils plus chance d'avoir au moins une sœur que les filles ?" (auquel cas la réponse serait différente).

Maintenant réfléchissons un peu, et faisons un petit modèle où on suppose que la taille de la fratrie ne dépend pas des sexes de la fratrie. Et comptons, pour mieux voir ce qui se passe. Et commençons par un cas simple, par exemple avec trois enfants,

  • (G,G,G) - personne n'a de sœur
  • (G,G,F), (G,F,G) ou (F,G,G) - la fille n'a pas de sœur, par contre, il y a deux garçons avec une sœur
  • (G,F,F), (F,G,F) ou (F,F,G) - dans chaque cas, il y a deux filles avec une sœur, et un garçon avec deux sœurs,
  • (F,F,F) - ce qui fait trois filles avec chacune deux sœurs.
Si on récapitule pour les garçons, on a 3 cas avec 2 garçons ayant 1 sœur et 3 cas avec 1 garçon ayant 2 sœurs, soit (à un facteur de normalisation près) un nombre espéré de sœurs de 12. Et pour les filles, on a 3 cas avec 2 filles ayant 1 sœur et 1 cas avec 3 filles ayant 2 sœurs, soit (au même facteur de normalisation près) un nombre espéré de sœur de 12.
Posons http://freakonometrics.blog.free.fr/public/perso4/soeur07.gif la taille de la fratrie, que l'on va supposer fixé. Essayons de formaliser les calculs qu'on vient de faire. S'il y a http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif filles, et http://freakonometrics.blog.free.fr/public/perso4/soeur11.gif garçons,
  • http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif filles et http://freakonometrics.blog.free.fr/public/perso4/soeur11.gif garçons - chacun des http://freakonometrics.blog.free.fr/public/perso4/soeur11.gif garçons a http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif sœurs, alors que chacune des http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif filles a http://freakonometrics.blog.free.fr/public/perso4/soeur13.gif sœurs,
Bref, dans le calcul de l’espérance pour les garçons, on doit avoir des http://freakonometrics.blog.free.fr/public/perso4/soeur12.gif et pour les filles, http://freakonometrics.blog.free.fr/public/perso4/soeur14.gif. Et en faisant intervenir le nombre de combinaisons possibles de http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif filles, et http://freakonometrics.blog.free.fr/public/perso4/soeur11.gif garçons, et par qui était notre facteur de proportionnalité de tout à l'heure (i.e. le nombre de possibilités d'agencement des sexes dans les familles avec http://freakonometrics.blog.free.fr/public/perso4/soeur07.gif enfants, quand on se focalise sur le sexe de tous, sauf un, qui est choisi au hasard, et qui est la personne qui compte ses soeurs), et en notant http://freakonometrics.blog.free.fr/public/perso4/soeur06.gif l'espérance du nombre de sœurs chez les garçons, et http://freakonometrics.blog.free.fr/public/perso4/soeur05.gif chez les filles, on a (merci à Pierre Antoine qui m'a fait remarqué une erreur dans mes calculs initiaux)

alors que

Un peu de calculs permet de montrer que ces deux quantités sont en fait égales... et valent
ce qui a du sens car il reste http://freakonometrics.blog.free.fr/public/perso4/soeur07.gif-1 personnes dans la fratrie (je peux renvoyer au livre de Ruma Falk, Understanding Probability and Statistics: A Book of Problems, page 46, qui passe par les fonctions génératrices pour une démonstration propre). Amusant, non ? En tous les cas suffisamment pour attirer mon attention alors que j'essaye de récupérer d'un soucis de santé, en dormant dans le lit de ma (grande) fille au milieu des peluches et des poupées...

- page 1 of 4