Category Archives: Maths

Découpage équitable de la pizza avec des olives

mais non papa, tu comprends jamais rien ! ok, c’est bien si on a autant de pizza tous les deux… mais on veut aussi avoir autant d’olives tous les deux ! et autant de morceaux d’ananas !”. Ouais, bon, on n’est pas couché là !

On va quand même essayer de la couper cette pizza (avant qu’elle ne refroidisse) !

  • Partage de la pizza et des olives

Commençons par supposer que seules les olives soient importantes. La première idée que j’aurais serait de faire une régression médiane, qui va découper le nuage de points (les olives) en deux…

Sous R (oui, j’ai installé R sur la table de la cuisine), c’est assez facile à mettre en œuvre… Pour des raisons de symétrie, je vais faire les régressions de l’abscisse sur l’ordonnée, et réciproquement…

> library(quantreg)
> lq=rq(y ~ x, tau = 0.5, data=D)
> abline(lq,lwd=5,col="red")
> sum(y0>predict(lq))
[1] 22
> sum(y0<=predict(lq))
[1] 22
>
> lq=rq(x ~ y, tau = 0.5, data=D)
> abline(b=1/lq$coefficients[2],
+ a=-lq$coefficients[1]/lq$coefficients[2],lwd=5,col="red")
> sum(x0>predict(lq))
[1] 22
> sum(x0<=predict(lq))
[1] 22

Le soucis est que si la régression n’est pas contrainte, il y a une répartition 50-50 des olives, mais les parts ne sont pas équitables,

> library(quantreg)
> lq=rq(y ~ 0+x, tau = 0.5, data=D)
> abline(lq,lwd=5,col="red")
> sum(y0>predict(lq))
[1] 20
> sum(y0<=predict(lq))
[1] 24
>
> lq=rq(x ~ 0+y, tau = 0.5, data=D)
> abline(b=1/lq$coefficients,a=0,lwd=5,col="red")
> sum(x0>predict(lq))
[1] 25
> sum(x0<=predict(lq))
[1] 19

et si je supprime la constante les parts sont équitables (la régression passe par le centre la pizza), mais le partage des olives n’est pas équilibré…

On peut aussi calculer pour tous les angles possibles ce que donnerait des partages de la pizza (qui passent par le centre). C’est ce que représente la courbe rouge ci-dessous: l’extérieur du disque le ratio des olives dépasse 1, à l’intérieur il est inférieur à 1. Les découpages optimaux correspondent au cas où le ratio vaut précisément 1… ce qui est possible si on coupe comme indiqué ci-dessous, à droite,

Je pense qu’on doit pouvoir montrer qu’il est toujours possible de couper la pizza en deux, de manière équitable, avec autant d’olives de part et d’autre (s’il y a un nombre pair d’olives).
oui mais papa… tu nous avais dit que tu savais pas couper en passant par le centre…”. Ah oui, c’est vrai ! Damned…. bon, utilisons la méthode décrite dans le billet précédant (ici), a priori ça doit être possible de trouver un sommet pour lequel le partage des olives est équitable. En fait, on pourrait aussi tirer au hasard des points dans la pizza, et compter les olives en faisant la répartition de manière alternée,

> sommet=matrix(NA,500000,3)
> for(s in 1:nrow(sommet)){
+ a0=runif(10)*2-1; b0=runif(10)*2-1
+ if(a0^2+b0^2)<=1){
+ slope1=(y0-b0)/(x0-a0)
+ I1=((slope1>0)&(slope1<1))|((slope1<(-1)))
+ sommet[s,]=c(sum(I1),a0,b0)
}}

On a alors en rouge l’ensemble des points à partir desquels on peut faire des découpages en 8, avec des pointes à 45° (en commençant par couper horizontalement)

Et ça semble marcher !

oui mais papa… tu as oublié les ananas…”.

  • Partage de la pizza et des olives et des ananas

En utilisant la méthode précédente, rien de plus simple ! En plus de compter la répartition des olives, on va compter la répartition des ananas…! Sur le dessin à gauche on a les olives et les ananas. A droite, les points rouges correspondent à des répartitions équitables des olives, et les bleus des répartitions équitables des ananas… Reste à trouver des points à la fois rouges et bleus,

ce sont finalement les points violets ci-dessous ! On a alors une méthode pour couper la pizza en deux parts de même taille, avec autant d’olives et autant d’ananas de part et d’autres !

Découpage équitable de la pizza

Bon, maintenant que l’année scolaire est passée, je vais arrêter un peu les billets sérieux pour revenir pour mes mésaventures à la maison. Par exemple, hier soir, c’était pizza. Et mes enfants aiment partager à deux une pizza (sur laquelle on a rajouté plein de choses qu’ils aiment, genre des olives, ou des morceaux d’ananas). Mais il faut un découpage équitable. Les enfants ne veulent pas un partage juste, mais un partage moitié-moitié.

Classiquement, je coupe la pizza en deux (cf figure ci-contre, avec un côté en bleu pour mon fils, et un en rouge pour ma fille, histoire de simplifier un peu). Mais comme me l’a fait noté mon fils, si mon couteau ne passe pas par le centre de la pizza (on va dire le centre du cercle parce qu’il va falloir faire des maths… pour les recettes de cuisine, il faut changer de site, genre ici).
Donc pour faire simple, la proba d’avoir deux parts rigoureusement égale est très très faible… car il faut couper suivant une droite qui passe exactement au centre (sinon, il y a quelqu’un qui est défavorisé, alors on tire au sort… a prioric’est juste juste, mais ex-post, il y en a toujours un ou une qui trouve le système injuste et qui veut qu’on retire au hasard…)
Fort heureusement, il existe une technique simple, et rigoureusement juste, permettant d’avoir des parts équitables…
L’idée est la suivante: on commence par couper horizontalement, n’importe où. Ensuite, on coupe verticalement (à angle droit, en prenant son équerre) de manière à récupérer 4 parts. Au point d’intersection,  on coupe à 45° direction nord-est, et 45° direction nord-ouest…. Une fois le découpage effectué, on réparti les parts, à tour de rôle, les parts bleues pour mon fils, et celles en rouges pour ma fille,

Car le truc génial est que le découpage – en plus d’être fort élégant – est parfaitement équitable ! Comme les enfants sont toujours sceptiques, il faut bien sûr leur prouver que oui, ils en auront exactement pareil. Pour ça on reprend le couteau…

Si on recoupe les parts suivant les traits-ci-dessus, chaque part bleue pourra être superposée à une part rouge (la preuve mathématique repose sur des simple notions de symétrie…). Génial non ? On a même le théorème du mardi soir, à savoir que “si on coupe une pizza en 8 parts, avec des pointes à 45° partir d’un point choisi arbitrairement dans la pizza, alors la somme alternée des parts sera égale” (cf ici). En 1967, Upton a posé ce problème,

Larry Carter et Stan Wagon affirmaient que ce résultat resterait vrai si on coupait en 12 ou 16 parts (avec des angles égaux pour les pointes), mais pas pour 6, 10, 14 ou 18 par exemple. Le blog du chou romanesco (ici et ) prétend que c’est vrai pour les découpages en un nombres pairs. En fait, Rick Mabry et Paul Deiermann ont un théorème assez général ( pour la version complète) et qui généralise aussi en dimension 3 (pour une calzone en quelque sorte). Mais je sens que ça n’intéresse plus trop les enfants, qui ont déjà fini la pizza….

Les plus curieux peuvent aussi consulter le court article ci-dessous,

Fin des débats sur la valeur de π

La rumeur avait agité pas mal de monde en août dernier lors du congrès international de mathématiques, à Hyderabad (même si à l’époque, ce sont surtout les médailles Fields qui avaient retenu toute l’attention, en France en tous les cas), mais finalement le communiqué de l’Union Internationale de Mathématiques (IMU) est finalement tombé hier soir: à compter du 1er juillet, π sera officiellement égal à 4.
Pour ceux qui ont suivi les débats (je n’ai eu que des bruits de couloir), Microsoft avait augmenté la pression ces derniers mois, quand le cap des 5000 milliards de décimales avait été franchi (ici). Comme l’avait dit Bill Gates avant le congrès à Hyderabad, bientôt la moitié de ma mémoire d’un processeur sera dédiée à stocker les décimales de π.

Et comme il l’avait souligné “la recherche sur les décimales de π allant plus vite que la recherche des améliorations de Windows, nous allons rapidement faire face à un choc informatique sans précédant”. Il avait comparé la situation au bug de l’an 2000 (en ajoutant que – comme il y a 11 ans – la transition se fera, selon lui, en douceur).
Dans le communiqué de l’IMU, il est mentionné que “π gardera son interprétation originelle « περίμετρος »” (périmètre en grec), et en particulier, l’IMU utilise la définition géométrique suivante:

http://freakonometrics.blog.free.fr/public/perso2/animation-pi.gif

Mais en quoi cela peut avoir un intérêt pour mon blog (hormis faire mon frimeur pour montrer que j’ai compris une démonstration géométrique). Tout simplement car π joue un rôle central en statistique et en probabilité (même si on a souvent tendance à l’oublier) au travers de la loi normale ! Pareil en risk management, y compris pour va valorisation des options, mais aussi le capital réglementaire calculé dans les accords de Bâle. Les banquiers se sont réjouit de l’annonce de l’IMU hier soir, car pour ceux qui l’auraient oublié, la loi normale intervient partout en finance . En particulier dans les calculs de quantiles. Rappelons que la densité s’écrit

https://freakonometrics.hypotheses.org/files/2015/12/dens-gauss-pi.gif (pour les utilisateurs de R, la version 2.12.3. sera lancée exceptionnellement plus tôt afin d’intégrer cette mise à jour – et je crois que MSExcel a prévu un addins qui sera bientôt en ligne, comme pour le passage à l’an 2000).
Si on regarde “à la main” ce que vont devenir les probabilités de dépassement de seuils, on obtient

> 1-pnorm(2)
[1] 0.02275013
> integrate(f=function(x){exp(-x^2/2)/sqrt(2*pi)},2,+Inf)
0.02275013 with absolute error < 1.5e-05
> integrate(f=function(x){exp(-x^2/2)/sqrt(2*4)},2,+Inf)
0.02016178 with absolute error < 1.3e-05

autrement dit, la probabilité de dépasser 2 va passer de 2.27% à 2.02%. Pour les quantiles à 99.5% (utilisées par les institutions financières), on a

> qnorm4(.995)
[1] 2.543701
> qnorm(.995)
[1] 2.575829

Autrement dit la baisse est finalement relativement faible. On peut alors s’attendre à une (très) légère baisse des fonds propres des banques et institutions financières (même si pour l’instant, les instances de Bâle ne se sont pas encore prononcées).

 

“Je comprends pas tout, mais j’aime bien”

http://freakonometrics.blog.free.fr/public/perso2/.2249021524.4_s.jpgPlusieurs fois, j’ai lu des commentaires sur le blog (ou sur le premier, à Rennes 1) commençant par quelque chose du genre “j’ai bien aimé même si j’ai pas tout compris“…
Étrangement, je pense à cette phrase régulièrement depuis quelques mois. En fait, en novembre dernier, j’avais acheté pour mon fils un livre écrit par Ivar Ekeland, Le chat au pays des nombres. S’il a bien aimé, c’est sans commune mesure avec la passion que voue ma fille pour ce livre ! Je dois lui lire au moins une fois par semaine ! Elle ne s’en lasse pas ! Et moi non plus d’ailleurs (je finirais par mieux connaitre ce livre que celui de théorie des jeux qu’il avait publié il y a un peu plus longtemps).

http://freakonometrics.blog.free.fr/public/perso2/.1296324_3197388_s.jpg http://freakonometrics.blog.free.fr/public/perso2/.1296324_3197389_s.jpg
Pourtant, le livre explique tout simplement que http://freakonometrics.blog.free.fr/public/perso2/ensembleQ.gif est dénombrable, en construisant une bijection entre http://freakonometrics.blog.free.fr/public/perso2/ensembleN.gif et http://freakonometrics.blog.free.fr/public/perso2/ensembleNN.gif, comme l’avait Cantor (et que l’infini, c’est quand même super grand).
http://freakonometrics.blog.free.fr/public/perso2/cantor_set.jpg

Et comme c’est un livre pour enfants, Ivar utilise l’approche proposée par David Hilbert (ici), à savoir l’hôtel infini. Chaque fraction (i.e. un élément de http://freakonometrics.blog.free.fr/public/perso2/ensembleQ.gif, en bleu au dessus) a une place dans une chambre de l’hôtel infini (dont le numéro de chambre est dans http://freakonometrics.blog.free.fr/public/perso2/ensembleN.gif, en rouge au dessus).
Le livre est génial, et je pousse tout le monde à le livre ! Mais je ne cesse de me demander ce que ma fille en retient, ce qu’elle comprend vraiment…  Car Ivar Ekeland est formidable: il arrive à parler de concept non triviaux à des enfants de 5 ans (pour les plus grands, je peux renvoyer ici pour une relecture du papier de John Nash sur l’interprétation d’une théorème de point fixe).
(Je pourrais aussi noter que ce livre me fait aussi m’interroger sur la notion d’impact factor des publications des chercheurs. Ivar est connu dans la communauté scientifique par bon nombre de papiers théoriques publiés dans les revues les plus prestigieuses, mais je doute qu’un chercheur ait passé autant de temps sur ses articles que ma fille sur son livre.)

A Million Random Digits: review of reviews

Recently on his blog (here), Robin mentioned an amazing book, called “A Million Random Digits” published by RAND corporation. The book was initially published in 1955, but RAND published a nice (and expensive) second edition.

A great thing is that on Amazon, there are several extremely interesting reviews of the book. E.g.

4.0 out of 5 stars Didn’t like the ending, February 10, 2009  By Damien Katz

Even though I didn’t really see it coming, the ending was kind of anti-climatic. But overall the book held my attention and I really liked the “10034 56429 234088” part. It’s nice to know I’m not the only one who feels that way.

5.0 out of 5 stars I found a typo, September 14, 2007  By fanfan

To whom do I write to report typographical errors? I noticed that the first “7” on the third line page 48 should be a “3”. The “7” that’s printed there now isn’t random. Other than that, this is really an excellent book.

5.0 out of 5 stars Superb and original plot, April 21, 2007  By Herr Tarquin Biskuitfaß

This one has a very unpredictable plot, sublime character development in a style that stubbornly defies any sort of development in its rare and iconoclastic brilliance, and is told remarkably with numbers instead of letters. Take, for example, this passage on page 202, “98783 24838 39793 80954”. I’m speechless. The symmetry is reminiscent of the I Ching, and it approaches a rare spiritual niveau lacking in American literature. It not only reads well, but it looks great too. I have a tattoo of page 214 on my arm, and I’m hoping to get 202 on my belly to celebrate my next birthday. It is an injustice that Rand Corporation has not received the Nobel Prize for Literature, nor even a Pulitzer.

3.0 out of 5 stars A serious reference work?, October 16, 2006  By BJ

For a supposedly serious reference work the omission of an index is a major impediment. I hope this will be corrected in the next edition.

1.0 out of 5 stars Not Nearly A Million, September 3, 2006  By Liron

This book does not even come close to delivering on its promise of one million random digits. My expectations were high after reading the first sentence, which contained ten unique digits. However, the author seems to have exhasted his creativity in this initial burst, because the other 99.999% of the book is filler in which those same ten digits are shamelessly reused!  If you are looking for a larger offering of numerals in various bases, I highly recommend “Peter Rabbit’s ABC and 123”.

3.0 out of 5 stars Wait for the audiobook version, October 19, 2006  By R. Rosini “Newtype”

While the printed version is good, I would have expected the publisher to have an audiobook version as well. A perfect companion for one’s Ipod.

5.0 out of 5 stars Wait for it…, February 10, 2009  By Cranky Yankee

It started off slow, single digit slow in the beginning but I stuck with it. I eventually learned all about the different numbers, 1,2,3,4,5,6,7,8,9 and 0 and their different combinations.  The author introduced them all a bit too quickly for my taste. I would have been perfectly happy with just 1,2,3,4 and 5 for the first 20,000 digits, but then again, I’m not a famous random-number author, am I?  After a while, patterns emerged and the true nature of the multiverse was revealed to me, and the jokes were kinda funny. I don’t want to spoil anything but you will LOVE the twist ending!  Like 4352204 said to 64231234, “2242 6575 0013 2829!”

Ok, I have to admit I tried to check a few of them (that’s my freaky part). For instance the first one is a fake: the two first numbers – for instance – never show up together (consecutively),

> DIGIT=read.table("
+ http://freakonometrics.blog.free.fr/public/data/digits.txt")
> DIGIT=DIGIT[,2:11]
> k=1
> I=apply(DIGIT[,1:2]==c(10034,56429),1,sum)==2
> for(k in 2:9){
+ I=cbind(I,apply(DIGIT[,k+0:1]==c(10034,56429),1,sum)==2)
+ }
> I0=which(apply(I,1,sum)>0)
> DIGIT[I0,]
 [1] V2  V3  V4  V5  V6  V7  V8  V9  V10 V11
<0 rows> (or 0-length row.names)

Nevertheless, I did have some fun reading those reviews. About the book, unfortunately I have to confess I stopped after 99998 appeared (the first time).

Will Rogers ou la magie des moyennes

Pour reprendre un exemple emprunté à@mathematicsprof, il existe un paradoxe intéressant (mais assez simpliste) sur les moyennes, appelé phénomène de Will Rogers. Le point de départ est une déclaration de Will Rogers, qui avait affirmé que “in the 1930’s, there was a mass migration of people from Oklahoma to California. As a result, the average I.Q. of both states…..increased”. Le paradoxe est connu en médecine sous le nom de “stage migration”, et se trouve mentionné dans bon nombre d’études académiques. Comment cela est-il possible ? En fait c’est assez simple… Considérons les deux groupes suivants

(1,2,3,4,5)(6,7,8,9,10)

Les moyennes sont respectivement de 3 et 8. Maintenant, si le 6 quitte le second groupe pour le premier,

(1,2,3,4,5,6)(7,8,9,10)

les moyennes deviennent alors 3,5 et 8,5. Formellement, si on considère deux groupes http://freakonometrics.blog.free.fr/public/perso2/will-1.gif et http://freakonometrics.blog.free.fr/public/perso2/will-2.gif, tels que http://freakonometrics.blog.free.fr/public/perso2/will-3.gif, si on transfère http://freakonometrics.blog.free.fr/public/perso2/will-6.gif, on détériore la moyenne des deux groupes, alors que si on transfère http://freakonometrics.blog.free.fr/public/perso2/will-8.gif, on améliore les deux….

Étonnant non (bon, pas tant que ça quand on voit la solution) ?

when Nuns or Hells Angels get in a plane

Today, at lunch, Matthieu told us a nice story (or call it a paradox if you like) about the probability to find you seat empty when you get in a place. 

  • a plane full of nuns

Assume that you are in the line to get in the airplane, you are the 100th in the line. The first one is scatter brained, he has his head in the clouds, and when he get in the airplane, he cannot remember where he should seat. His strategy is then extremely simple: he seats randomly in the plane. So he picks up randomly a seat, and he waits.

Then come 98 nuns (one by one). And nuns are extremely polite: if there is someone in their seat (the one that is on the ticket they have) then they do not complain, and pick up another seat randomly (among those available, of course). Then you arrive. The question is simple: what is the probability that someone is seated at your seat ?

Any idea…?

Maybe I should give more time to do the maths… and tell another story…

  • a plane full of Hells Angels

Consider almost the same problem as the one mentioned above. Except that now, it is not 98 nuns that are getting in the plane, but 98 Hells Angels. So the problem here is that Hells Angels are slightly less polite than nuns. When they find someone seating on the seat they should have, they do not shyly move to another seat, but they grunt and then our scatter brained man (who is actually seating in their seat) has to move somewhere else. And the question is the same: you are the 100th person to get in the plane, what is the probability that someone is seated at your seat ?Any idea….?

The important point is that the problem is exactly the same (at least from a mathematical point of view, maybe not for the stewardess, or from the guy who enter first in the plane). The point is that, at each time, there could be only one person (or less) seating in a seat which is not his or hers (in the sense that if we compare the list of the passenger at any time, and the list of seats taken, there should be only one – or less – difference). The difference in the two story is that in the first case, it will be a nun, while in the second one, it will be our shy guy.

  • Let us run simulations

If we do not see how to get that probability analytically, let us run some R code,

> set.seed(1)
> n=100; TEST=rep(NA,100000)
> for(s in 1:100000){
+ OCCUPIED=rep(FALSE,n)
+ OCCUPIED[sample(1:n,size=1)]=TRUE
+ for(j in 2:(n-1)){
+ FREE=which(OCCUPIED==FALSE)
+ if(OCCUPIED[j]==TRUE){OCCUPIED[sample(FREE,size=1)]=TRUE}
+ if(OCCUPIED[j]==FALSE){OCCUPIED[j]=TRUE}
+ }
+ TEST[s]=OCCUPIED[n]==TRUE
+ }
> mean(TEST)
[1] 0.49878

Here, we clearly see that the problem is the same (either with nuns or Hells Angels): we do not care about who will change his/her seat, but we just look at seats that are available… So the program is valid for the two problems (and the solution will then be the same). Another point is that the probability looks extremely simple: one over two !

  • an analytical expression

Consider the Hells Angels problem (for notations). Let http://freakonometrics.blog.free.fr/public/perso2/nonnes1.gif denote the probability that, at time http://freakonometrics.blog.free.fr/public/perso2/nonne6.gif, our shy guy is sitting in my seat. When he gets in the plane, the probability that he gets to my seat is

http://freakonometrics.blog.free.fr/public/perso2/nonne2.gif
Then, the probability that, after ith passenger’s entrance, our guy is sitting in my own seat is (since the initial proof was not correct, I remove it, see below for a nice proof) One can get that

http://freakonometrics.blog.free.fr/public/perso2/avion-ec-01.gif
So, we can get the probability that, when I get in, our guy is sitting in my own seat as
http://freakonometrics.blog.free.fr/public/perso2/avion-ec-07.gif

http://freakonometrics.blog.free.fr/public/perso2/avion-ec-08.gif

Hence, there is one chance out of two that my seat will be free… (which is what we got with Monte Carlo simulations).

But a faster proof is to observe that, in the Hells Angels case, our guy will be kicked out until he reaches either his seat, or mine. Since those two events are equiprobable, there is one chance out of two that he seats in my seat (and since no Hells Angel will seat in mine, only this first guy can). So the probability that someone is in my seat when I get in is one half.

Nice isn’t it ? And thanks Matthieu for the problem  (with his friend Claude’s solution with the Hells Angels, and Olivier and Renaud for their comments) !

What is the optimal strategy to marry the best one ?

Valentine’s day is a nice opportunity to post on hot and sexy topics… Well, it’s also an important day that I should not miss, probably as much as Saint Patrick’smy wife’s birthday. And as I mentioned last week (here), it is difficult to get the distribution of the age of marriage on the internet… So maybe we can build up a small model, to understand when do girls decide to get married… Consider a young girl who knows that he will not meet thousands of men willing to marry her (actually, one can consider the opposite point of view, with young man who can find only http://freakonometrics.hypotheses.org/files/2015/12/mariage01.png girls willing to marry him, the problem can be assumed as symmetric, especially if I do not want to get feminist leagues on my back).

Assume that http://freakonometrics.hypotheses.org/files/2015/12/mariage01.png men agree to marry her. Of course, among those http://freakonometrics.hypotheses.org/files/2015/12/mariage01.png men, our girl wants to marry the “best” one (assume that men can be ranked objectively). Of course, she cannot meet the “best” guy immediately, so men are met randomly, and after each “interview“, either she reject him (forever, we assume she cannot get back and admit she made a mistake), or agree to marry him. An important assumption is that rejected men cannot be recalled.

From a mathematical point of view, we need to find the optimal stopping time. Here, the problem is slightly different compared with that one (with optimal time to get a bonus) or this one (with the optimal time to sit in a bar and have a beer). Here, we do not give “grades” to guy. The only thing that is observed is their relative ranks. Our girl cannot know if she’s meting the best of all men (out of http://freakonometrics.hypotheses.org/files/2015/12/mariage01.png), but she knows if this one is better than the ones she already met. From a mathematical point of view, at time http://freakonometrics.hypotheses.org/files/2015/12/mariage02.png, she knows the relative rank of http://freakonometrics.hypotheses.org/files/2015/12/mariage02.png (compared with the first http://freakonometrics.hypotheses.org/files/2015/12/mariage04.png), not his absolute rank. We also assume that http://freakonometrics.hypotheses.org/files/2015/12/mariage01.png is known.

The optimal strategy is that she has to reject automatically the first http://freakonometrics.hypotheses.org/files/2015/12/mariage04.png (some kind of calibration period), and then, starting at time http://freakonometrics.hypotheses.org/files/2015/12/mariage02.png, she will marry the best over the ones she has already met.
So assume that our girl already met http://freakonometrics.hypotheses.org/files/2015/12/mariage04.png guys, and decided to reject all of them. So now she’s trying to see if the http://freakonometrics.hypotheses.org/files/2015/12/mariage02.png can be the optimal time to stop, and start looking seriously ….For an arbitrary cut-off http://freakonometrics.hypotheses.org/files/2015/12/mariage02.png, the probability that the best applicant will show up at some time http://freakonometrics.hypotheses.org/files/2015/12/wedd01.gif is http://freakonometrics.hypotheses.org/files/2015/12/wedding01.gif

http://freakonometrics.hypotheses.org/files/2015/12/wedding02.gif

i.e.

http://freakonometrics.hypotheses.org/files/2015/12/wdeeing03.gif

The http://freakonometrics.hypotheses.org/files/2015/12/wedd02.gif term is because there is only one “best” guy, and the http://freakonometrics.hypotheses.org/files/2015/12/wedd03.gifis the probability that he shows up at time http://freakonometrics.hypotheses.org/files/2015/12/wedd01.gif (this can be visualized below)

Thus, we can write

http://freakonometrics.hypotheses.org/files/2015/12/wedding04.gif

i.e.

http://freakonometrics.hypotheses.org/files/2015/12/wedding05.gif

Thus, since the minimum of http://freakonometrics.hypotheses.org/files/2015/12/mariage18.png is obtained when http://freakonometrics.hypotheses.org/files/2015/12/mariage19.png, which is the optimal time to stop (or here to start seeking), i.e. 36.7%.

Hence, the best strategy is to reject automatically the first http://freakonometrics.hypotheses.org/files/2015/12/mariage20.png=37% of the candidates (which is the maximum value of the function above), and then to select the first one (if possible) that is better than all previous candidates.

Consider the following Monte Carlo procedure: assume that she rejects – automatically – the first http://freakonometrics.hypotheses.org/files/2015/12/mariage02.png (we consider a loop with all possible values for http://freakonometrics.hypotheses.org/files/2015/12/mariage02.png) and then gets married with the first one who is the best one she’s seen during the calibration period (or overall, which is the same),

n=100
ns=1000000
MOY1=MOY2=rep(NA,n)
for(m in 2:(n-1)){
WHICH=rep(NA,ns); MARIAGE=rep(0,ns)
for(s in 1:ns){
Z=sample(1:n,size=n,replace=FALSE)
mx=max(Z[1:m])
STOP=FALSE
for(k in (m+1):n){
if((Z[k]>mx)&(STOP==FALSE)){
WHICH[s]=k
    STOP=TRUE
MARIAGE[s]=1
}
}
}
HIS=WHICH[is.na(WHICH)==FALSE]
TH=table(HIS)
MOY1[m]=mean(HIS)
MOY2[m]=mean(HIS)*mean(MARIAGE)
THH=rep(NA,100)
THH[as.numeric(names(TH))]=as.numeric(TH)/ns
}

If we run it over all possible http://freakonometrics.hypotheses.org/files/2015/12/mariage02.png we get

http://freakonometrics.hypotheses.org/files/2015/12/mariage-anim.gif

The “distribution” (in green) can be seen as the probability to marry the guy of level http://freakonometrics.hypotheses.org/files/2015/12/mariage06.png, given that the first http://freakonometrics.hypotheses.org/files/2015/12/mariage02.png were rejected. The sum is not one since there is a non null probability to marry no one. Actually, the probability to get married is the following

The more she waits, the smaller the probability of getting married. But on the other hand, the more she waits, the “better” the husband…. On the graph below is plotted the rank of the guy she marries, if she gets married (it was actually the vertical plain line in red on the animation)

So there is a trade-off. If not getting married gives a 0 satisfaction (lower than finally marrying anyone), and if marrying the guy with rank http://freakonometrics.hypotheses.org/files/2015/12/mariage06.png gives here satisfaction http://freakonometrics.hypotheses.org/files/2015/12/mariage06.png ,we have

(it was the vertical doted line in red on the animation). So it looks like it is optimal to test the first 35-38% men, and then to marry the best one she finds (if he is better than the best one she met during the “testing” procedure). So our previous analysis looks correct…

Now to go further, I have to admit that this model is known in academic literature as the secretary problem. In 1989, Thomas Ferguson wrote a nice paper inStatistical Science entitled who solved the secretary problem (here). Anthony Mucci published also an article in the Annals of Probability on possible extensions, in 1973 (here), or Thomas Lorenzen (there) in 1981. This problem is definitively an interesting one !

Tennis and risk management

As mentioned already here, while we were going to Québec City for the workshop, we had interesting discussions in the car, and Maciej mentioned an article recently published in The Actuary,

Hence, I wanted to discuss (extremely) rare event probabilities in tennis. The story is simple: in June 2010, at Wimbledon, Nicolas Mahut and John Isner have played the longest match ever. 980 points, 11 But first of all, we need a dataset. Thanks to Duncan Murdoch, I have been able to run a short code to build up a dataset:

CITIES=c("berlin","madrid","paris","rolandgarros","wimbledon","sydney",
"beijing","shanghai","singapore","tokyo","melbourne","melbourne-indoor")
YEARS=1970:2009
BASE0=data.frame(YEAR=NA,TRNMT=NA,LENGTH=NA,SETS=NA)
for(i in 1:length(CITIES)){
for(j in 1:length(YEARS)){
city=CITIES[i]
year=YEARS[j]
localization = paste("http://www.resultsfromtennis.com/",
year,"/atp/",city,".html",sep="")
essai = try(readLines(localization), silent=TRUE)
ERROR404=FALSE
if(inherits(essai, "try-error")){ERROR404=TRUE}
if(ERROR404==FALSE){
B=scan(localization,"character")
SETS=NA
LENGTH=NA
if(length(B)>270){
I=(substr(B,1,10)=="class=rez>")
sum(I)
X0=B[I]
X3=as.numeric(substr(X0,11,13))
X2=as.numeric(substr(X0,11,12))
X1=as.numeric(substr(X0,11,11))
X0=X3
X0[is.na(X3)==TRUE]=X2[is.na(X3)==TRUE]
X0[is.na(X2)==TRUE]=X1[is.na(X2)==TRUE]
JL=c(which(substr(B,1,9)=="class=nl>"),length(B))
IL=which(substr(B,1,10)=="class=rez>")
IC=cut(IL,JL)
base=data.frame(IC,X0)
LENGTH=as.numeric(tapply(X0,IC,sum))
SETS=as.numeric(tapply(X0,IC,length))/2}
BASE=data.frame(YEAR=year,TRNMT=city,LENGTH,SETS)
BASE0=rbind(BASE0,BASE)}}}
write.table(BASE0,"BASE-TENNIS-TOTAL.txt")

Here I consider only tournaments where players have to win 3 sets (and actually more tournaments than those in the code above), and I have something like a bit more than 72,000 matches,

> I=is.na(TENNIS$LENGTH)==FALSE
> BT=TENNIS[I,]
> nrow(BT)
[1] 72754
> maxr=function(x){max(x,na.rm=TRUE)}
> T=paste(BT$TRNMT,BT$YEAR)
> DUREE=tapply(BT$SETS,T,maxr)
> LISTE=names(DUREE[DUREE>3])
> BT=BT[T%in%LISTE,]

so, if we look briefly at matches over 35 years, we have the following boxplot (one boxplot per year),

The red line being the epic Isner-Mahut match in June 2010 (4-6, 6-3, 7-6, 6-7, 70-68, i.e. 183 games, here for the score card).

If we study theory (e.g. from Paul Newton and Kamran Aslam), a lot of results can be obtained for the expected value of the number of games, but if we want to study extremely rare events, we should generate Markov chains (with a lot of generation since the probability should be extremely small). But how many ? Consider below matches with more than 50 games,

The tail plot (over 50), i.e. the log-log Pareto plot indicates that it will be difficult to study tails,

and similarly with the Hill plot (assuming that tails are Pareto type….)

Anyway, if we want to study tails, we should consider a threshold high enough. For instance, with a threshold at 68 (we keep only 24 match), we have

> seuil=68+0.25
> GPD1=gpd(X,seuil,method = "ml")
> GPD2=gpd(X,seuil,method = "pwm")
>
> xi=GPD1$par.ests[1]
> mu=seuil
> beta=GPD1$par.ests[2]
> x=180
> P=exp((-1/xi)*log(1 + (xi * (x - mu))/beta))
> as.numeric((1-GPD1$p.less.thresh)*P)
[1] 5.621281e-09
>
> xi=GPD2$par.ests[1]
> mu=seuil
> beta=GPD2$par.ests[2]
> x=180
> P=exp((-1/xi)*log(1 + (xi * (x - mu))/beta))
> as.numeric((1-GPD2$p.less.thresh)*P)
[1] 3.027095e-09

I.e. the probability that one match last more than 183 games is 1 chance over a billion… With, say, 2500 match per year, that gives us a return period of 400 years. So yes, we might say that this way a rare event… So perhaps, generating several billions of chains, it should be possible to get a more precise estimation of the probability to play 183 games in a single match…

Finding roots of functions in actuarial science

The following simple code can be used to find roots of functions (based on the secant algorithm),

secant=function(fun, x0, x1, tolerence=1e-07, niter=500){
for ( i in 1:niter ) {
	x2 <- x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0))
	if (abs(fun(x2)) < tolerence)
		return(x2)
	x0 <- x1
	x1 <- x2
}}

It can be interesting in actuarial science, e.g. to find the actuarial rate so that to present values are equal. For instance, consider the following capital, given only if the insured is still alive (this example was initially considered here). We would like to find the rate so that the probable discounted value is 600,

> Lx=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",
+ header=TRUE,sep=";")
> capital=c(100,100,125,125,150,150)
> n=length(capital)
> x=0.035
> X=45
> f=function(x){
+ capital.act=capital*(1/(1+x))^(1:n)
+ PROBA=Lx[((Lx[,1]>X)*(Lx[,1]<=(X+n)))==1,2]/Lx[(Lx[,1]==X)==1,2]
+ return(sum(capital.act*PROBA))}
>
> f1=function(x){f(x)-600}
> secant(f1,0,0.1)
[1] 0.06022313
> f(0.06022313)
[1] 600

*

MacGyver et la théorie des groupes

L’autre jour, quelqu’un me félicitait pour mon blog, c’était une étudiante en maths, qui disait comprendre certaines choses en probas et en stats (entre autres sur les applications possibles)…. mais pour la théorie des groupes, elle ne voyait toujours pas d’intérêt…

J’étais bien entendu ravi que mon blog serve à quelque chose, mais en fait la théorie des groupes peut avoir des applications amusantes (bon, tout dépend aussi de ses centres d’amusement). Par exemple MacGyver pourrait avoir eu besoin de connaître des résultats sur les groupes pour sortir des prisonniers des geôles afghanes…

Considérons le gardien de prison sadique qui annonce le “jeu” suivant (Anatole, un collègue ici me l’avait suggéré l’autre jour en me mettant au défi de trouver une solution…. mais faut pas pousser, je suis plus fort que MacGyver ! bon, et il faut dire que le jeu est classique…on le trouve expliqué ici ou ):

Tous les prisonniers (disons qu’ils sont 20) se voient attribuer un numéro, entre 1 et 20. Le premier prisonnier est appelé et se retrouve face à 20 portes, derrière lesquelles sont cachés des numéros (allant de 1 à 20 là encore). Il a droit d’en ouvrir 10 (au maximum) et il doit retrouver son numéro. S’il n’y arrive pas,tous les prisonniers sont exécutés (lui compris). S’il y arrive, il n’est pas sauvé pour autant car le commandant appelle le second prisonnier, qui se retrouve à son tour devant les 20 portes, et doit trouver son numéro… Etc.

Les prisonniers commencent à paniquer, car tout le monde doit trouver son numéro pour être sauvé. En raisonnant rapidement, un des prisonniers explique qu’un seul prisonnier aurait 50% de chances d’y arriver. Avec deux, comme ils ont tous les deux 50% de chances, on tombe à 25%, etc. Avec n prisonniers, on est à 1/2n, qui est de l’ordre d’une chance sur un million avec 20 prisonniers…. Et qu’ils ont de la chance de n’être que 20 ! Avec 30 prisonniers, ils
auraient une chance sur un milliard de s’en sortir (ce qui est faible, on en conviendra).

Mais MacGyver, qui avait pris le cours Théorie des Groupes 101 à l’université (et qui avait R installé dans sa montre digitale) dit

  • pas du tout, j’ai une stratégie qui nous donne environ une chance sur trois de tous s’en sortir (vu comme ça c’est pas terrible, mais par rapport à 1 chance sur un milliard, c’est vachement bien),
  • et accessoirement, cette probabilité ne dépend pas du nombre de prisonnier (tant que l’on puisse ouvrir moitié moins de portes qu’il y a de prisonniers)

Bon, sur le moment, les autres prisonniers se fichent du fait que même s’ils étaient 1000 ils auraient une chance sur trois de survivre, ils veulent comprendre la stratégie de Mac… Comme le rappelle Mac, le commandant ne fait que choisir une permutation de https://perso.univ-rennes1.fr/arthur.charpentier/latex/prison99.png, même s’il est sadique (ou parce qu’il est sadique, les profs de maths sont sadiques, tout le monde sait ça). Et toute permutation peut se décomposer en un produit ds combinaisons cycliques (on parle de la décomposition en produit de cycles à supports disjoints, c’est un truc utilise je crois pour les amateurs de Rubik’s cube),

Par exemple, considérons la permutation suivante (c’est la première que j’ai eu sous R),

https://perso.univ-rennes1.fr/arthur.charpentier/latex/prison01.png

Elle se décompose en 5 permutations cycliques

https://perso.univ-rennes1.fr/arthur.charpentier/latex/prison02.png

La plus grande est la dernière, et elle permute sur 7 nombres. Autrement dit, si chacun utilise cette stratégie, il ouvrira au plus 7 portes. Le point important est qu’ils ne peuvent pas être pris dans une autre permutation que la leur, car ils ont pris leur porte, donc au pire, il tire le numéro de la première porte qu’ils ont ouvert… qui est le numéro qu’ils doivent tirer pour survivre ! Regardons sur notre exemple

> set.seed(1)
> portes = sample(1:n,size=n,replace=FALSE)
>  prisonnier=1
>  essai = c(prisonnier)
>  nessai = 1
>  i = porte[prisonnier]
>  while(nessai <= 50) {
+  essai = c(essai, i)
+  if(i == prisonnier) {
+  cat(“Le prisonnier”, prisonnier, “a ouvert”, nessai, “portes\n”)
+  break; }
+  else {
+  i = porte[i]}
+ n essai =n essai+1 }
Le prisonnier 1 a ouvert 4 portes
>  cat(“Le prisonnier”, prisonnier, “a ouvert les portes suivantes”, 
paste(essai[-length(essai)], collapse=” – “), “\n”)
Le prisonnier 1 a ouvert les portes suivantes 1 – 6 – 14 – 10

soit ici

Le prisonnier 1 a ouvert 4 portes
Le prisonnier 1 a ouvert les portes suivantes 1 – 6 – 14 – 10
Le prisonnier 2 a ouvert 7 portes
Le prisonnier 2 a ouvert les portes suivantes 2 – 8 – 9 – 19 – 18 – 17 – 12
Le prisonnier 3 a ouvert 2 portes
Le prisonnier 3 a ouvert les portes suivantes 3 – 11
Le prisonnier 4 a ouvert 5 portes
Le prisonnier 4 a ouvert les portes suivantes 4 – 16 – 7 – 15 – 5
Le prisonnier 5 a ouvert 5 portes
Le prisonnier 5 a ouvert les portes suivantes 5 – 4 – 16 – 7 – 15
Le prisonnier 6 a ouvert 4 portes
Le prisonnier 6 a ouvert les portes suivantes 6 – 14 – 10 – 1
Le prisonnier 7 a ouvert 5 portes
Le prisonnier 7 a ouvert les portes suivantes 7 – 15 – 5 – 4 – 16
Le prisonnier 8 a ouvert 7 portes
Le prisonnier 8 a ouvert les portes suivantes 8 – 9 – 19 – 18 – 17 – 12 – 2
Le prisonnier 9 a ouvert 7 portes
Le prisonnier 9 a ouvert les portes suivantes 9 – 19 – 18 – 17 – 12 – 2 – 8
Le prisonnier 10 a ouvert 4 portes
Le prisonnier 10 a ouvert les portes suivantes 10 – 1 – 6 – 14
Le prisonnier 11 a ouvert 2 portes
Le prisonnier 11 a ouvert les portes suivantes 11 – 3
Le prisonnier 12 a ouvert 7 portes
Le prisonnier 12 a ouvert les portes suivantes 12 – 2 – 8 – 9 – 19 – 18 – 17
Le prisonnier 13 a ouvert 2 portes
Le prisonnier 13 a ouvert les portes suivantes 13 – 20
Le prisonnier 14 a ouvert 4 portes
Le prisonnier 14 a ouvert les portes suivantes 14 – 10 – 1 – 6
Le prisonnier 15 a ouvert 5 portes
Le prisonnier 15 a ouvert les portes suivantes 15 – 5 – 4 – 16 – 7
Le prisonnier 16 a ouvert 5 portes
Le prisonnier 16 a ouvert les portes suivantes 16 – 7 – 15 – 5 – 4
Le prisonnier 17 a ouvert 7 portes
Le prisonnier 17 a ouvert les portes suivantes 17 – 12 – 2 – 8 – 9 – 19 – 18
Le prisonnier 18 a ouvert 7 portes
Le prisonnier 18 a ouvert les portes suivantes 18 – 17 – 12 – 2 – 8 – 9 – 19
Le prisonnier 19 a ouvert 7 portes
Le prisonnier 19 a ouvert les portes suivantes 19 – 18 – 17 – 12 – 2 – 8 – 9
Le prisonnier 20 a ouvert 2 portes
Le prisonnier 20 a ouvert les portes suivantes 20 – 13

En adaptant un code développé ici par Matt Asher, on peut d’ailleurs visualiser tout ça sur un petit dessin,
,

Bref, la grande question est “quelle est la probabilité que la plus grande sous permutation cyclique soit sur un ensemble de taille supérieure à 10” ? On peut faire un peu de Monte Carlo pour voir…. c’est ce que fait Matt. Mais sinon on peut aussi faire du dénombrement (ce qui est pareil que les probas). Le nombre de permutations contenant un cycle de longueur supérieure à 10 (strictement) se calcule de la manière suivante: soit k la taille du cycle, alors je peux choisir les https://perso.univ-rennes1.fr/arthur.charpentier/latex/prison00.png nombres de https://perso.univ-rennes1.fr/arthur.charpentier/latex/prison03.png manières. Il y a alors https://perso.univ-rennes1.fr/arthur.charpentier/latex/prison04.png permutations circulaires de ces k nombres, et on permute les https://perso.univ-rennes1.fr/arthur.charpentier/latex/prison05.png nombres qui reste comme on veut, ce qui donne https://perso.univ-rennes1.fr/arthur.charpentier/latex/prison06.png possibilités. Bref, j’en suis à

https://perso.univ-rennes1.fr/arthur.charpentier/latex/prison07.png

permutations ayant un cycle de longueur k, soit en tout,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/prison0_.png

Or comme on a en tout https://perso.univ-rennes1.fr/arthur.charpentier/latex/prison08.png permutations possibles, la probabilité de s’en sortir vivant est

https://perso.univ-rennes1.fr/arthur.charpentier/latex/prison09.png

Et là où Mac est fort, c’est qu’effectivement, avec n prisonniers

https://perso.univ-rennes1.fr/arthur.charpentier/latex/prison10.png

Soit environ 30% de chances de s’en sortir….
Moralité, la théorie des groupes peut être marrante et intéressante…. Bon, c’est la seule application à la vraie vie qui me vient à l’esprit, et à moins d’être enfermé en prison avec un geôlier sadique, c’est dur de trouver des applications. N’empêche que la solution est passionnante (mais MacGyver a toujours été mon héros pour ça…)

When should I optimally shoot at my son ? (part 2)

Following my previous post of yesterday, online here, assume now that I do not know if my son came when I turned my back at time, and missed me… Then the payoff function is the one propose by Vincent, i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel100.png

In that particular case,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel20.png

becomes

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel101.png

i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel102.png

If we draw those functions, on [0,1], the optimal value is solution of https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel103.png

i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel104.png (we focus only on solutions in [0,1]). Because here the game is symmetric, my son should also shoot at time https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel105.png
Thus, the payoff is then

 https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel106.png

Since we consider here a  zero-sum game, this cannot be a solution of the game. So the game does not have pure strategy solution.

Assume that now I have a mixed strategy, i.e. a distribution of the optimal time to shot. My strategy has distribution https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel110.png, with density https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel111.png (we assume here that the density exists, or we seek only solution that are differentiable). Assume further that there exists https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel112.png>0 such that the support of my optimal shooting time (the time to shoot is now a random variable) is (https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel112.png,1] (or [https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel112.png,1] since we assume that https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel110.png is differentiable). There is a discussion at the end of Vincent’s post where he needs that assumption, at the end. Actually, I think we can make it now, since we can rationally assume that
I will not shot at time 0 (even on a neighborhood of 0 since I have zero chance to hit my son).
The expected payoff function, assuming that my son shoots at time y is

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel1113.png

Since the zero-sum game is symmetric, again, the expected payoff should be zero. It comes that necessarily,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel120.png

if https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel121.png. Hence, if we differentiate (with respect to y), we have

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel122.png

and if we differentiate one more time, it comes

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel123.png

i.e. a general solution should be of the form https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel130.png.
Here, we have the same solution as the one considered in Vincent’s blog. His solution is obtained as follows (with slightly different expressions) conditional to https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel03.png, my expected payoff is

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel140.png

i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel131.png

With a simple integration by parts,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel132.png

where https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel134.png, i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel133.png

Thus,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel135.png

So, if we want to be indifferent to https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel03.png‘s strategy, https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel141.png, where

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel142.png

with https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel143.png,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel144.png

Consider solutions https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel150.png, then https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel151.png=1, i.e. either https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel152.png=1 and then https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel153.png is constant, or https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel152.png=-1. This means that https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel154.png is in proportional to https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel155.png.
If we substitute in the equation we had, initially, it comes that

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel160.png

i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel161.png

If we consider https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel03.png=a and https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel03.png=1, it comes that https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel112.png=1/3 while https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel162.png=1/4 (but we don’t really care about that normalizing constant).
It means that we should not start shooting before 1/3 of the tank is fulled. Actually, it makes sense, since

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel163.png

if https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel03.png<1/3 (while https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel164.png=0 if  https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel03.png>1/3).

When should I optimally shoot at my son ? (part 1)

when playing with water guns ! It is not a reverse oedipal fantasy, it is simply a (stochastic) game. Actually, that game was mentioned here, in the context of a duel with laser guns, and the answer is there. Vincent (alias @Vicnent) presented the problem this way: consider a single combat between two single warriors (a duel) with laser guns. The probability to kill the other one is proportional to time https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel01.png (duel started at time 0). Initially, the have no chance to kill the other one, and after one hour, they are certain to kill the other one. At time https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel01.png, they kill the other one with probability https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel01.png (and miss him with probability https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel04.png). What is the optimal time to shot ?

An alternative can be the following (yes, I prefer personal experience to science fiction): consider two warriors playing with water gun. We start with empty guns, and assume that both warriors fill their guns at two different tap water. The more water we have in guns, the more likely the other one will be wet. What is my optimal strategy to stop filling the gun, and start shooting at the other player ?
In order to formalize the game, assume that I will win 1 point if my son gets wet (and not me), I’ll get -1 if I am wet, and not him, and 0 if we are both wet (in Vincent’s game, we have the same payoff matrix, but I find odd to say that I have -1 if I die).
Let https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel02.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel03.png denote time where players shoot. Vincent derived the following expected payoff function: for me (I shoot at time https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel02.png) the payoff is

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel05.png

If we know that the other player shoot us, but missed us, it should be different: if I shot first, but missed my son, then he should wait until the end, and so, the payoff should be (thanks Jérôme)

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel06.png

Actually, we can consider a more general game, where probabilities are not proportional to https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel01.png, but functions https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel10.png for me (https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel13.png stands for dad), and https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel11.png for the other player (here my son https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel12.png). The expected payoff becomes,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel07.png

Thus my best strategy is obtained by solving

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel20.png

that is

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel24.png

Let

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel23.png

If https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel11.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel10.png are assumed to be continuous (we assume that water goes continuously into the tank of the gun), then there is https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel26.png such that

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel25.png

If https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel30.png, then https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel32.png and

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel34.png

and so https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel35.png. So

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel37.png

Similarly, if https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel31.png, then https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel33.png and

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel39.png

and so https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel40.png. So

https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel42.png

So finally, the optimal strategy is to shoot at the same time, https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel01.png.
In the game considered initially by Vincent, https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel43.png. Hence, https://perso.univ-rennes1.fr/arthur.charpentier/latex/duel01.png=1/2.
But it is also possible that I did not hear my son coming, shooting at me (and missing me). So in that case, the payoff is different, and closer to what Vincent proposed… but that will be in another post…

Mandelbrot, fractals and counterexamples in applied probability

Benoît Mandelbrot died yesterday. Like most of the blogs dealing with applied mathematics, it looks like I have to mention this event. Unfortunately, I don’t know much about fractals…

The first time I heard about Mandelbrot and chaos was when I have been working on fractional time series (see eghere). Murad Taqqu gave a very interesting short course in Paris, and I have been using it in two papers (actually one more should appear soon in Climate Change).

The second time was in Québec (city), five years ago, when Roger Nelsen gave a talk on copulas with fractal suppport (here). By that time, we were finishing our paper with Alessandro (in mathematical finance, and limiting theorems). I remember adding a Remark in the paper (that can be found here), since using that kind of copulas was a nice way to show that, without sufficient regularity conditions, the limit we were looking for had no sense.

A few months after, with Johan, in a paper on pitfalls on lower tail dependence for Archimedean copulas (here), we used again this fractal construction (here applied to Archimedean copulas) to find a nice counterexample to a (false) theorem on regular variation. Again, the goal was to understand how the dependence structure of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mandel01.png given https://perso.univ-rennes1.fr/arthur.charpentier/latex/mandel02.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mandel03.png changed when https://perso.univ-rennes1.fr/arthur.charpentier/latex/mandel04.png goes to 0. In the case of Archimedean copulas, the copula of the conditional pair is still Archimedean (with another generator, except for Clayton copula). The graph below show how https://perso.univ-rennes1.fr/arthur.charpentier/latex/mandel05.png changes, as https://perso.univ-rennes1.fr/arthur.charpentier/latex/mandel04.png decreases… Actually, I draw https://perso.univ-rennes1.fr/arthur.charpentier/latex/mandel06.png since Archimedean generators are not unique.

I have also decided to plot https://perso.univ-rennes1.fr/arthur.charpentier/latex/mandel07.png.

Here, we see that there is no way of talk about a possible limit for the conditional copula because of a fractal behavior in 0 of the generator (even if my fractal are not as nice as the one you can find on the internet….). So thanks Benoît for giving us a nice toy to build interesting counterexamples !

EM and mixture estimation

Following my previous post on optimization and mixtures (here), Nicolas told me that my idea was probably not the most clever one (there).
So, we get back to our simple mixture model,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM01.png

In order to describe how EM algorithm works, assume first that both https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png and  https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.pngare perfectly known, and the mixture parameter is the only one we care about.

  • The simple model, with only one parameter that is unknown

Here, the likelihood is

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM04.png

so that we write the log likelihood as

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM05.png

which might not be simple to maximize. Recall that the mixture model can interpreted through a latent variate https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png (that cannot be observed), taking value when https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM07.png is drawn from https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png, and 0 if it is drawn from https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png. More generally (especially in the case we want to extend our model to 3, 4, … mixtures), https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM08.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM09.png.
With that notation, the likelihood becomes

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM10.png

and the log likelihood

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM11.png

the term on the right is useless since we only care about p, here. From here, consider the following iterative procedure,
Assume that the mixture probability https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png is known, denoted https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png. Then I can predict the value of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png (i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM08.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM09.png) for all observations,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM14.png

So I can inject those values into my log likelihood, i.e. in

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM15.png

having maximum (no need to run numerical tools here)

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM16.png

that will be denoted https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM17.png. And I can iterate from here.
Formally, the first step is where we calculate an expected (E) value, where https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.pngis the best predictor of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM19.png given my observations (as well as my belief in https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png). Then comes a maximization (M) step, where using https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png, I can estimate probability https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png.

  • A more general framework, all parameters are now unkown

So far, it was simple, since we assumed that https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png and  https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png were perfectly known. Which is not reallistic. An there is not much to change to get a complete algorithm, to estimate https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM30.png. Recall that we had https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png which was the expected value of Z_{1,i}, i.e. it is a probability that observation i has been drawn from https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png.
If https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png, instead of being in the segment https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM31.png was in https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM32.png, then we could have considered mean and standard deviations of observations such that https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png=0, and similarly on the subset of observations such that https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png=1.
But we can’t. So what can be done is to consider https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png as the weight we should give to observation i when estimating parameters of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png, and similarly, 1-https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.pngwould be weights given to observation i when estimating parameters of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png.
So we set, as before

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM33.png

and then

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM34.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM35.png

and for the variance, well, it is a weighted mean again,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM36.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM37.png

and this is it.

  • Let us run the code on the same data as before

Here, the code is rather simple: let us start generating a sample
> X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z  = sample(c(1,2,2),size=n,replace=TRUE)
> X2=4+X20
> X = c(X1[Z==1],X2[Z==2])
then, given a vector of initial values (that I called https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png and then https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM99.png before),
>  s = c(0.5, mean(X)-1, var(X), mean(X)+1, var(X))
I define my function as,
>  em = function(X0,s) {
+  Ep = s[1]*dnorm(X0, s[2], sqrt(s[4]))/(s[1]*dnorm(X0, s[2], sqrt(s[4])) +
+  (1-s[1])*dnorm(X0, s[3], sqrt(s[5])))
+  s[1] = mean(Ep)
+  s[2] = sum(Ep*X0) / sum(Ep)
+  s[3] = sum((1-Ep)*X0) / sum(1-Ep)
+  s[4] = sum(Ep*(X0-s[2])^2) / sum(Ep)
+  s[5] = sum((1-Ep)*(X0-s[3])^2) / sum(1-Ep)
+  return(s)
+  }
Then I get https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png, or https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM99.png. So this is it ! We just need to iterate (here I stop after 200 iterations) since we can see that, actually, our algorithm converges quite fast,
> for(i in 2:200){
+ s=em(X,s)
+ }

Let us run the same procedure as before, i.e. I generate samples of size 200, where difference between means can be small (0) or large (4),

Ok, Nicolas, you were right, we’re doing much better ! Maybe we should also go for a Gibbs sampling procedure ?… next time, maybe….