Category Archives: 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 !

Fractals and Kronecker product

A few years ago, I went to listen to Roger Nelsen who was giving a talk about copulas with fractal support. Roger is amazing when he gives a talk (I am also a huge fan of his books, and articles), and I really wanted to play with that concept (that he did publish later on, with Gregory Fredricks and José Antonio Rodriguez-Lallena). I did mention that idea in a paper, writen with Alessandro Juri, just to mention some cases where deriving fixed point theorems is not that simple (since the limit may not exist).

The idea in the initial article was to start with something quite simple, a the so-called transformation matrix, e.g.

https://latex.codecogs.com/gif.latex?T=\frac{1}{8}\left(\begin{matrix}1&%200%20&%201%20\\%200%20&%204%20&%200%20\\%201%20&%200&1\end{matrix}\right)
Here, in all areas with mass, we spread it uniformly (say), i.e. the support of https://latex.codecogs.com/gif.latex?T(C^\perp) is the one below, i.e. https://latex.codecogs.com/gif.latex?1/8th of the mass is located in each corner, and https://latex.codecogs.com/gif.latex?1/2 is in the center. So if we spread the mass to have a copula (with uniform margin,)we have to consider squares on intervals https://latex.codecogs.com/gif.latex?[0,1/4]https://latex.codecogs.com/gif.latex?[1/4,3/4] and https://latex.codecogs.com/gif.latex?[3/4,1],

Then the idea, then, is to consider https://latex.codecogs.com/gif.latex?T^2=\otimes^2T, where  https://latex.codecogs.com/gif.latex?\otimes^2T is the tensor product (also called Kronecker product) of https://latex.codecogs.com/gif.latex?T with itself. Here, the support of https://latex.codecogs.com/gif.latex?T^2(C^\perp) is

Then, consider https://latex.codecogs.com/gif.latex?T^3=\otimes^3T, where https://latex.codecogs.com/gif.latex?\otimes^3T is the tensor product of https://latex.codecogs.com/gif.latex?T with itself, three times. And the support of https://latex.codecogs.com/gif.latex?T^3(C^\perp) is

Etc. Here, it is computationally extremely simple to do it, using this Kronecker product. Recall that if https://latex.codecogs.com/gif.latex?%20%20%20%20%20\mathbf{A}=(a_{i,j}), then

https://latex.codecogs.com/gif.latex?%20%20%20%20%20\mathbf{A}\otimes\mathbf{B}%20=%20\begin{pmatrix}%20a_{11}%20\mathbf{B}%20&%20\cdots%20&%20a_{1n}\mathbf{B}%20\\%20\vdots%20&%20\ddots%20&%20\vdots%20\\%20a_{m1}%20\mathbf{B}%20&%20\cdots%20&%20a_{mn}%20\mathbf{B}%20\end{pmatrix}

So, we need a transformation matrix: consider the following https://latex.codecogs.com/gif.latex?4\times4 matrix,

> k=4
> M=matrix(c(1,0,0,1,
+            0,1,1,0,
+            0,1,1,0,
+            1,0,0,1),k,k)
> M
[,1] [,2] [,3] [,4]
[1,]    1    0    0    1
[2,]    0    1    1    0
[3,]    0    1    1    0
[4,]    1    0    0    1

Once we have it, we just consider the Kronecker product of this matrix with itself, which yields a https://latex.codecogs.com/gif.latex?4^2\times4^2 matrix,

> N=kronecker(M,M)
> N[,1:4]
[,1]  [,2] [,3] [,4]
[1,]     1    0    0    1
[2,]     0    1    1    0
[3,]     0    1    1    0
[4,]     1    0    0    1
[5,]     0    0    0    0
[6,]     0    0    0    0
[7,]     0    0    0    0
[8,]     0    0    0    0
[9,]     0    0    0    0
[10,]    0    0    0    0
[11,]    0    0    0    0
[12,]    0    0    0    0
[13,]    1    0    0    1
[14,]    0    1    1    0
[15,]    0    1    1    0
[16,]    1    0    0    1

And then, we continue,

> for(s in 1:3){N=kronecker(N,M)}

After only a couple of loops, we have a https://latex.codecogs.com/gif.latex?4^5\times4^5 matrix. And we can plot it simply to visualize the support,

> image(N,col=c("white","blue"))

As we zoom in, we can visualize this fractal property,

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…

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.

Bayes is playing Russian roulette

There was (once again) a nice puzzle inhttp://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).

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.

25% ? sérieusemement, vous y croyez ?

(pour reprendre les statistiques postées par @guybirenbaum l’autre jour)

Réfléchissons deux minutes (ou un peu plus, c’est ce qu’on a fait ce midi avec @J_P_Boucher). On vous interroge sur un sondage: si vous croyez aux sondages, vous le dites. Mais si vous n’y croyez pas… vous répondez n’importe quoi (sinon c’est que vous y croyez, non ?). Bon essayons de formaliser ça afin de mieux comprendre… Soit https://latex.codecogs.com/gif.latex?X la réponse donnée, et https://latex.codecogs.com/gif.latex?Y la vérité (ce que pense vraiment la personne interrogée). On sait que

https://latex.codecogs.com/gif.latex?\mathbb{P}(X=C)=1-\mathbb{P}(X=CP)=\frac{1}{4}

Mais ce qu’on veut calculer, c’est https://latex.codecogs.com/gif.latex?p=%20\mathbb{P}(Y=C), la probabilité qu’une personne y croit vraiment aux sondages. Pour ça, on sait que

https://latex.codecogs.com/gif.latex?\mathbb{P}(X=C)=\mathbb{P}(X=C|Y=C)\cdot%20\mathbb{P}(Y=C)

https://latex.codecogs.com/gif.latex?+\mathbb{P}(X=C|Y=CP)\cdot%20\mathbb{P}(Y=CP)

Pour le premier terme, une personne qui y croit le dira (sinon c’est qu’elle n’y croit pas), donc https://latex.codecogs.com/gif.latex?%20\mathbb{P}(X=C|Y=C)=1. Par contre, une personne qui n’y croit pas peut dire n’importe quoi. Disons que https://latex.codecogs.com/gif.latex?%20\mathbb{P}(X=C|Y=CP) vaut https://latex.codecogs.com/gif.latex?%20\alpha. Dans ce cas

https://latex.codecogs.com/gif.latex?%20\frac{1}{4}=p+\alpha(1-p)

ou encore

https://latex.codecogs.com/gif.latex?%20p=\frac{1/4-\alpha}{1-\alpha}

Supposons par exemple qu’une petite proportion des personnes qui n’y croient pas disent y croire, disons 10%. Dans ce cas, la vraie probabilité qu’une personne croit au sondage serait plutôt 16%. Si https://latex.codecogs.com/gif.latex?%20\alpha était un peu plus grand (25%), on serait juste tombé sur la proportion des rigolos qui répondent n’importe quoi aux sondages, car personne n’y croit !

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),

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 ?

It is a simple rejection technique algorithm. It is possible to run some code to check the answer. The algorithm which return the taste of the last egg remaining is

> lastchocolate=function(dark=80,milk=20){
+ s=1
+ while(dark+milk>1){
+ if(s==1){
+ (eatnow=sample(c("D","M"),prob=c(dark,milk),size=1))
+ if(eatnow=="D"){dark=dark-1};
+ if(eatnow=="M"){milk=milk-1};
+ eatbefore=eatnow;s=2}
+ if(s==2){
+ if(dark+milk>1){
+ s=1;
+ eatnow=sample(c("D","M"),prob=c(dark,milk),size=1)
+ if(eatnow==eatbefore){s=2
+ eat=eatnow;
+ if(eatbefore=="D"){dark=dark-1};
+ if(eatbefore=="M"){milk=milk-1}}
+ }}
+ }
+ return(c(dark,milk))}

If we run it 2,000 times, we obtain

> set.seed(1)
> m=lastchocolate(dark=80,milk=20)
> for(s in 1:1999){m=cbind(m,lastchocolate(dark=80,milk=20))}
> apply(m,1,sum)
[1] 1022 978

So it looks like we have half chance to end up with a dark chocolate egg, and half chance to end up with a milk chocolate egg.

Let us prove that result… Let  denote the number of milk chocolate and the number of dark chocolate eggs, when we start. Consider an inductive proof of the fact that the probability has to be . The first step is when . Then

out of chocolates, the probability to pick a milk chocolate egg is . Assume that it is  for all pairs  such that  and .
Assume that after some steps, there are  and  chocolates, with  (again, at least one egg has been eaten). The probability to have  is
Similarly, the probability to have  is
So, the probability that both are strictly positive is then
Then we can use our inductive assumption. Thus, the overall probability that the last egg is a milk one is
where the part on the left is  and the second one is . This probability is exactly one half (straightforward).

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…

Where hiding if you don’t want to get wet ?

Following the previous post, two additional remarks. Following a comment by@cosi, I have investigated quickly a binomial fit to the distribution of the number of people not getting wet, with a fixed number of players on the field. It looks like it should be a binomial distribution with a fixed probability (2/3) and with size parameter affine in the number of players. @guigui suggested some connexion with with “birds on a wire” problem (see e.g. http://www.cut-the-knot.org/)

n=p=rep(NA,20)
for(i in 1:20){
NSim=10000
N=Vectorize(NOTWET)(n=rep(3+2*i,NSim))
n[i]=mean(N)/(1-var(N)/mean(N))
p[i]=1-var(N)/mean(N)
}
plot(seq(5,43,by=2),n,col="red",type="b")

for the implied size parameter above, and below the implied probability parameter.

plot(seq(5,43,by=2),p,col="blue",type="b")

(as functions of the number of players). I’d be glad to get more details on that 2/3 probability.

Now, let us investigate another question sent by email: “Where should you hide if you don’t want to get wet ?” A first idea could be the following: given that some players are already on the field, where should I go if I do not want to get wet ? Below are some simulations for 7 or 25 players (already on the field). The red area is the area so that I will become someone’s target (perhaps even the target of two players…). The green area is the safe zone.

(with 7 players above, and 25 below)

It looks like, on the border, it might be safer than in the middle of the field. But we have to confirm that intuition… or at least see if that intuition is valid.

Based on what was done the other day, it is possible to look where people that got wet were located (instead of counting dry players as done in the previous function). So here, we simply look where non wet players were standing

NOTWET=function(n,p){
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)
whonotwet=( (1:n) %notin% names(table(dmin)) )
#plot(x[-whonotwet],y[-whonotwet],pch=19,col="blue",type="p")
#points(x[whonotwet],y[whonotwet],pch=19,col="red")
M=matrix(NA,p,p);u=seq(0,1,by=1/p)
for(i in 1:p){
for(j in 1:p){
M[i,j]=sum((x[whonotwet]>=u[i])&(x[whonotwet]<u[i+1])&
(y[whonotwet]>=u[j])&(y[whonotwet]<u[j+1]))
}}
return(M)}

based on function

"%notin%" <- function(x, y) x[!x %in% y]

On a given grid, we count people playing the game that ended dry (with might avoid boundary bias on nonparametric smooth estimator of distribution, as we’ll see later on). For instance with 11 players,

M11=matrix(0,25,25);
for(s in 1:100000){
M11=M11+NOTWET(11,25)
}

Then we can plot the distribution, on the field,

COL=rev(heat.colors(101)); p=25
u=seq(0,1,by=1/p)
plot(0:1,0:1,col="white",xlab="",ylab="")
for(i in 1:p){
for(j in 1:p){
polygon(c(u[i],u[i],u[i+1],u[i+1]),
 c(u[j],u[j+1],u[j+1],u[j]),border=NA,
col=COL[trunc(M11[i,j])/max(M11)*100+1])
}}

Red means a lot of non-wet people (i.e safer zones). Graphs below are with 7 and 11 players respectively (from the left to the right)

with the following distribution on the diagonal: corners are almost 4 times safer than the middle of the field, with 7 players,

Below are plotted distributions of locations of non-wet players when the total number of players was either 25 (on the left) and 101 (on the right)

with again on the diagonal

Hence, the border is rather safe, but next to the border, it is no safe any longer: is someone is standing right on the border, he will probably shoot at you: there is no one behind him ! This explains the stange behavior on the borders (and corners, thanks JP for the intuitive explanation).
But would it be completely different with a field shaped as a disk ?

using the previous technique of working on a fixed grid (or correcting for boundary bias, since the disk might cover only a fraction of the grid-square), or keeping coordinates of non-wet players, and using standard kernel-based estimator of the distribution (the light yellow circle outside the disk is simply due to bias of the kernel estimator on the border)

NOTWET=function(n){
x=(runif(n*20)*2-1)*1
y=(runif(n*20)*2-1)*1
I=which((x^2+y^2<1))
x=x[I];y=y[I]
x=x[1:n];y=y[1:n]
(d=as.matrix(dist(cbind(x,y),
method = "euclidean",upper=TRUE)))
diag(d)=999999
dmin=apply(d,2,which.min)
whonotwet=( (1:n) %notin% names(table(dmin)) )
return(cbind(x[whonotwet],y[whonotwet]))
}

M=t(c(0,0))
for(s in 1:10000){
M=rbind(M11,NOTWET(25))
}
M=M[-1,]

library(ks)
HP=matrix(c(.001,0,0,.001),2,2)
K=kde(x=M11, H=HP)
image(K$eval.points[[1]],K$eval.points[[2]],K$estimate2,
col=rev(heat.colors(101)),xlim=c(-1,1),ylim=c(-1,1))

 

And note that the distribution of the number of players ending the game dry is the same, for a square field, or a disk,

NOTWET2=function(n){
x=(runif(n*20)*2-1)*1
y=(runif(n*20)*2-1)*1
I=which((x^2+y^2<1))
x=x[I];y=y[I]
x=x[1:n];y=y[1: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)}

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=100000
Nsquare=Vectorize(NOTWET)(n=rep(25,NSim))
Ndisk=Vectorize(NOTWET2)(n=rep(25,NSim))
Tsq=table(Nsquare)
Tdk=table(Ndisk)
plot(as.numeric(names(Tsq)),Tsq/NSim,
type="b",col="red")
lines(as.numeric(names(Tdk)),Tdk/NSim,
type="b",pch=4,col="blue")


But so far, it was still simple… I wonder what it might become if we consider a non-convex place, with walls, where player might hide…. Next time, a post on indoor paint-ball !

Sunday evening, stupid games…

This evening, while I was about to wash the dishes, I heard my elders starting a game (call them Him and Her)
Him: “I have picked – in my head – a number, lower than 50. Try to guess…”
Her: “No way, too difficult…”
Him: “You can try five different numbers…”
Her: “.,. um … No, no way…”
Me: “Wait… each time we suggest a number, you tell us if yours is either above, or below ?”
You can see me coming clearly, can’t you ? Using a simple subdivision rule, we have a fast algorithm (and indeed, if I have to choose between washing the dishes and playing with the kids…)
Him: “um…. ok”
Her: “Daddy, are you sure we will win ?”
Me: “Well… I cannot promise that we will win… but I am rather sure [sic] that we will win quite frequently: more gains than losses…” (I guess).
Her: “Great ! I am playing with daddy…”

Him: “um…. wait, is it one of you trick, again ? I don’t to play anymore… Do you want to see the books we’ve chosen at the library ?”
Her: “Sure…”
Me: “What ? no one wants to see if I was right ? that we have indeed more than 50% chances to win…”
Him and her: “No !”
The point of that story ? If we listen to kids, science will not go forward, trust me. But I am curious… I want to see if my intuition was correct. Actually, the intuition was based on the fact that

> 2^5
[1] 32 
> 2^6
[1] 64

so in 5 or 6 steps the algorithm of subdivision should converge. I guess… I mean, I do not know for sure, since 50 is not a power of 2, so it might be difficult, each time, to split in two: we have to deal only with integers here…
To be sure, let us substitute my laptop to my son… to pick up numbers, randomly (yes, sometimes I feel like I am Doctor Tenma, 天馬博士). The algorithm is simple: there are bounds, and at each stop I should suggest the middle of the interval. If the middle is not an integer, I suggest either the integer below or the integer above (with equal probabilities).

cutinhalf=function(a,b){
m=(a+b)/2
if(m %% 1 == 0){m=m}
if(m %% 1 != 0){m=sample(c(m-.5,m+.5),size=1)}
return(round(m))}

The following functions runs 10,000 simulations, and tells us how many times, out of 5 numbers suggested, we got the good one.

winning=function(lower=1,upper=50,tries=5,NS=100000){
SIM=rep(NA,NS)
for(simul in 1:NS){
interval=c(lower,upper)
(unknownnumber=sample(lower:upper,size=1))
success=FALSE
for(i in 1:tries){
picknumber=cutinhalf(interval[1],interval[2])
if(picknumber==unknownnumber){success=TRUE}
if(picknumber>unknownnumber){interval[2]=picknumber}
if(picknumber<unknownnumber){interval[1]=picknumber}
#print(c(unknownnumber,picknumber,success,interval))
};SIM[simul]=success};return(mean(SIM))}

It looks like the probability that we got the good number is higher than 60%,

> winning()
[1] 0.61801

Which is not bad. And if the upper limit was not 50, but something else, the probability of winning would have been the following.

VWN=function(n){winning(upper=n)}
V=Vectorize(VWN)(seq(25,100,by=5))
plot(seq(25,100,by=5),V,type="b",col="red",ylim=c(0,1))


Actually, after losing a couple of times, I am rather sure that my son would have to us that we can suggest only four numbers. In that case, the probability would have been close to 30%, as shown on the blue curve below (where four numbers only can be suggested)

Anyway, as intuited, with five possible suggestions, we were quite likely to win frequently. Actually with a probability of almost 2 out of 3…and 1 out of 3 if my son had decided to pick an number between 1 and 100, or only 4 possible suggestions… Those are quite large actually, when we think about it. It reminds me that McGyver story I mentioned a few months ago… Anyway, calculating probabilities is nice, but I still have to wash the dishes…

Maths can be cool (to impress your kids)

Just imagine that your kids need some help, to prepare fishes for April 1st, like

Her: “please, Daddy, help us to draw some fishes”

Me: “Sure, Daddy is a champion, actually, I do that everyday at work: drawing fishes – and more generally nice stuff – is exactly Daddy’s job”.

OK, no need to talk neither about Talbot’s curvesellipse negative pedal curvenor Burleigh’s ovals, unless you don’t want to scare them, e.g.

t=seq(0,2*pi,length=100) 
b=.8
c=sqrt(1-b^2)
x=cos(t)-c*sin(t)^2
y=(1-2*c^2+c*cos(t))*sin(t)/b
plot(x,y)

From now on, it is rather simple to draw fishes,

t=seq(0,2*pi,length=100)
y=cos(t)-sin(t)^2/sqrt(2)
x=cos(t)*sin(t)
plot(x,y,type="l",axes=FALSE,xlab="",ylab="")
polygon(c(-2,-2,2,2),
c(-2,2,2,-2),col="light blue",border=NA)
polygon(x,y,col="red",border=NA)
axis(1)
axis(2)
lines(x,y,type="l")

so we can easily get nice fishes,

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)

Supposons par exemple que http://freakonometrics.blog.free.fr/public/perso5/proba02.gif, conditionnellement à http://freakonometrics.blog.free.fr/public/perso5/proba03.gif, suive une loi normale centrée, d’écart-type http://freakonometrics.blog.free.fr/public/perso5/proba08.gif,

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

où http://freakonometrics.blog.free.fr/public/perso5/proba03.gif est une variable (strictement) positive. Presque sûrement http://freakonometrics.blog.free.fr/public/perso5/proba05.gifquelle que soit la loi de http://freakonometrics.blog.free.fr/public/perso5/proba03.gif. Et donc http://freakonometrics.blog.free.fr/public/perso5/proba04.gif. Mais si la loi de http://freakonometrics.blog.free.fr/public/perso5/proba03.gif présente des queues “trop épaisses” (sans pour autant être extrême pour autant d’ailleurs), il est fort possible que http://freakonometrics.blog.free.fr/public/perso5/proba02.gif ne soit plus intégrable ! Par exemple (c’est ce que propose Enis (1973)), on peut supposer que http://freakonometrics.blog.free.fr/public/perso5/proba03.gif suit une loi Gamma de paramètres 1/2 et 2, i.e. de densité

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

Dans ce cas, http://freakonometrics.blog.free.fr/public/perso5/proba02.gif suit une loi de Cauchy, qui n’est pas d’espérance finie. C’est pénible, hein ?

Bon, d’un autre côté, c’est très lié à la définition de l’espérance conditionnelle retenue. En l’occurrence, j’utilisais l’espérance de la loi conditionnelle. Mais je ferais peut être un billet un jour sur les diverses constructions de l’espérance conditionnelle.

Des calculs en base 12 à la magie des log-rendements

Un des gros avantages dans la culture anglo-saxonne, c’est que culturellement, les anglais ont été de grands amateurs du système duodécimal (en base 12). Par exemple, il y a 12 pouces dans un pied, ou le fait que jusqu’aux années 70, il y avait 12 pences dans un schilling… Et ce que j’apprécie, c’est que ce système pousse à manipuler les fractions (2, 3, 4, 6 étant des diviseurs de 12) alors que le système décimal pousse à utiliser les nombres à virgule. Et compte tenu de la proximité avec les États-Unis, même si ici le système métrique a été adopté, il reste des fonds de ces calculs en base 12.
Par exemple, à l’école à Montréal, les tables d’addition et de multiplication vont jusqu’à 12. Aussi, l’autre soir, je faisais réviser les tables (de multiplication) de 11 et 12 à mon fils, et on a passé du temps sur http://freakonometrics.blog.free.fr/public/perso5/12-18.gif et http://freakonometrics.blog.free.fr/public/perso5/12-19.gif. Car sa première réponse à combien font huit fois douze ? a été ben…. cent ?. J’ai bien entendu protesté (c’est un euphémisme) pour la forme, tout en sachant qu’il est naturel de croire que http://freakonometrics.blog.free.fr/public/perso5/12-01.gif, ou que http://freakonometrics.blog.free.fr/public/perso5/12-02.gif. N’est ce pas ce qu’on pense quand on apprend qu’un prix a augmenté de 20% avant de baisser de 20% le lendemain ? Car c’est la même erreur que l’on commet. Et mine de rien, raisonner avec des rendements, c’est pénible…. (tous ceux qui ont discuté taux d’intérêt avec un banquier le savent).
Pour rappel, on définie le rendement comme la variation relative, i.e.

http://freakonometrics.blog.free.fr/public/perso5/12-03.gif

L’intérêt de cette formule est évident

http://freakonometrics.blog.free.fr/public/perso5/12-04.gif

Et si on regarde ce qui se passe sur deux jours

http://freakonometrics.blog.free.fr/public/perso5/12-05.gif

Aussi, si un jour on a un rendement de http://freakonometrics.blog.free.fr/public/perso5/12-06.gif, et le lendemain un rendement de http://freakonometrics.blog.free.fr/public/perso5/12-07.gif

http://freakonometrics.blog.free.fr/public/perso5/12-08.gif

Autrement dit, si on subit une baisse de 10% puis une hausse de 10% (ou le contraire d’ailleurs), on ne revient pas à la position initiale, car http://freakonometrics.blog.free.fr/public/perso5/12-09.gif(ou ce que je voyais avec mon fils quand on notait que http://freakonometrics.blog.free.fr/public/perso5/12-10.gif). En fait, on perd de l’argent… Damned, c’est pénible ça, non ? D’autant plus que si on répète cette alternance (entre une baisse de 10% puis une hausse de 10%) sur un an, au final, on a perdu 85% de valeur… Autrement dit, avec des rendements d’espérance nulle, et symétriques, la volatilité à elle seule fait perdre de l’argent…
Si on s’intéresse maintenant aux log-rendement,

http://freakonometrics.blog.free.fr/public/perso5/12-11.gif

cette fois, on note que

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

(on retrouve très explicitement la croissance exponentielle des prix). En faisant des développements limités, i.e. si les rendements sont faibles, alors on observe que http://freakonometrics.blog.free.fr/public/perso5/12-13.gif. En effet

http://freakonometrics.blog.free.fr/public/perso5/12-15.gif

Aussi, si http://freakonometrics.blog.free.fr/public/perso5/12-14.gif vaut 0.05, on parlera (abusivement) de hausse de 5%. Si on regarde une évolution sur deux jours,

http://freakonometrics.blog.free.fr/public/perso5/12-16.gif

Autrement dit, si on subit une baisse de 10% (en log-rendement) puis une hausse de 10% (ou le contraire), on revient à la position initiale (et donc on a une variation nulle) car http://freakonometrics.blog.free.fr/public/perso5/12-20.gif.
C’est probablement pour ça que les financiers utilisent des log-rendements: ils n’ont pas (en toute rigueur) la même interprétation que la variation relative, mais au moins on peut faire les calculs de manière très intuitive.