Digital resources in the Social Sciences and Humanities OpenEdition Our platforms OpenEdition Books OpenEdition Journals Hypotheses Calenda Libraries OpenEdition Freemium Follow us

Course on risk measures (in French)

The course on risk measure, in Luminy, starts at 16.00 on Monday (here). The slides can be found here,

Note that additional references can be downloaded on the internet, e.g. the short course on risk measures by Freddy Delbaen (here) or the article from the Encyclopedia of quantitative finance, by Hans Föllmer and Alexander Schied (there). See also here for the paper by Jean Marc Tallon, Johanna Etner and Meglena Jeleva, on decision theory under uncertainty.

Plus intéressants que les nobels: les ig-nobels

Il y a quelques semaines, Alessandro Pluchino, Andrea Rapisarda et Cesare Garofalo ont obtenu un des IG nobels (ici, probablement la plus haute distinction scientifique qui soit) pour avoir proposé une étude numérique du principe de Peter. Comme le notait Rue89, ils ont fait la “démonstration mathématique que la promotion aléatoire permet de gagner en efficacité dans une entreprise” ().

  • Le principe de Peter

Historiquement, tout vient d’un livre de Laurence Peter et Raymond Hull, paru en 1969, et intitulé précisément Peter principle. On y retrouve une formulation connu en management, avec le temps, tout poste sera occupé par un incompétent incapable d’en assumer la responsabilité“. Pour cela, on part de principes très simples, et que tout le monde admettra. Le premier est que quelqu’un de compétent à un poste est promu (à un niveau supérieur). Le second est que quelqu’un d’incompétent à un poste n’est pas promu à un niveau supérieur.

Aussi, un employé ne peut pas rester à un poste où il est compétent (puisqu’il est promu) et à l  suite des promotions, il devrait finir par atteindre un poste auquel il sera incompétent, il ne recevra plus de promotion, il restera donc indéfiniment à ce poste (pour lequel il est incompétent). Moralité, tous les postes finissent par être occupés par des employés incompétents, et à un instant donné on arrive à la formule lapidaire annonçant que la majorité du travail est effectuée par des salariés n’ayant pas encore atteint leur “seuil d’incompétence“.

Bon, c’est un peu simpliste, supposant qu’il existe une relation d’ordre, en particulier qu’une personne bonne à un niveau donné serait forcément bonne aux niveaux “inférieurs”. Il existe en particulier un compétence définie en une dimension unique, et donc avec un ordre naturel,

  • Le papier de Pluchino, Rapisarda et Garofalo

Le papier d’Alessandro Pluchino, Andrea Rapisarda et Cesare Garofalo (en ligne ici) essayent de tester numériquement ce principe, et le valident. Ils font en particulier un test très intéressant sur la stratégie des promotions. Ils montrent – sous quelques hypothèses de modélisation – que la stratégie visant à promouvoir le “meilleur” n’ait pas forcément optimale, et qu’on peut avoir intérêt à avoir une stratégie mixte en retenant parfois le plus mauvais… avec un seuil de mélange proche de 50%. Cette stratégie peut aussi se voir comme un stratégie simple: soit on offre une promotion au meilleur,  soit on offre une promotion par tirage au sort…. et comme c’est un peu ce qui se passe dans les tournois de soccer (si on pense que s’en remettre au tir au but c’est s’en remettre au hasard) on a un moyen simple pour tester cette théorie….

  • Étude sur les matchs de coupe du monde de soccer

J’ai regardé les résultats sur les coupes du monde allant de 1982 à 2010, en phase finale, soit 8+4+2=14 matchs par coupe, sauf en 1982, où je n’ai que deux matchs. Soit 100 matchs en tout. A chaque match, soit le gagnant s’est fait sur le score (le vainqueur gagne) soit sur le hasard (aux tirs au buts). Et pour tous ces matchs, je regarde si l’équipe qualifié a gagné, ou pas, le match suivant. J’obtiens la matrice suivante,

matchs gagnant perdant total
au mérite 46 34 80
au hasard 4 16 20
total 50 50 100

Si on fait un test d’indépendance du chi-deux sur ce tableau de contingence, on obtient qu’il n’y pas indépendance entre les variables,
> M=matrix(c(46,4,34,16),2,2)
> chisq.test(M)
Pearson’s Chi-squared test with Yates’ continuity correction
data:  M
X-squared = 7.5625, df = 1, p-value = 0.00596
Autrement dit, le fait d’être gagnant – ou perdant – à un match (entre les quarts et la finale) dépend fortement de la manière dont on s’est qualifié… Promouvoir en faisant un tirage au sort ne semble donc pas vraiment optimal, l’équipe qualifiée par hasard ne passant souvent pas l’étape suivante… Mais je suis intéressé par toute analyse sportive de ce principe, donc si quelqu’un a des pistes, je suis preneur ! (car les données ne manquent pas)

70% des salariés français touchent moins de 1 500 € par mois ?

Hier, @Louisa_A me demandait de commenter l’information comme quoi “70% des salariés français touchaient moins de 1 500 € par mois“. Tout d’abord, rappelons que selon l’INSEE, en 2008, le salaire médian net en France est de l’ordre de 1655 € par mois pour un temps plein (cf ici par exemple). Comme c’est la médiane (je renvoie ici pour une bonne blague du nouvel obs sur la médiane, comme quoi il ne faut pas toujours croire ce qu’on lit), c’est que 50% des français touchent ce salaire, ou plutôt toucheraient ce salaire s’il travaillaient à temps plein…. Le hic est que généralement l’insee normalise les données, et ne parle que de taux plein…

Mais heureusement l’insee met aussi a disposition des bases, en particulier la base dads2007, ici, avec un fichier de 25Mo (zippé) , contenant près de 2 millions d’observations, dont le temps travaillé en 2007 et le salaire (ou plutôt une tranche de salaire). Le code est plutôt simple (même si j’ai travaillé sur la première moitié de la base, suite à des soucis d’importation sur mon petit portable),

> library(foreign)
> base=read.dbf("D:\\r-data\\salaries07.dbf")
> TRANCHE=c(0,100,300,600,1000,1300,2000,3000,
+ 4000,5000,7000,9000,11000,12500,14000,16000,
+ 18000,20000,24000,30000,1000000000)
> salinf=TRANCHE[1:20]
> salsup=TRANCHE[2:21]

Bref,
on peut récupérer des salaires médians, par temps plein (afin de voir si nos données sont représentatives), en notant qu’un temps plein c’est 1820 heures dans l’année (cf ici)

> Xinf=salinf[base$TRNNETO+1]/base$NBHEUR_TOT*1820/12
> Xsup=salsup[base$TRNNETO+1]/base$NBHEUR_TOT*1820/12
> quantile(Xinf,.5,na.rm=TRUE)
 50% 
1484.745 
> quantile(Xsup,.5,na.rm=TRUE)
 50% 
1685.185

Comme nous ne sommes pas trop loin des 1655 indiqués, j’aurais tendance à faire confiance à nos chiffres (en se plaçant proches de la borne supérieure d’ailleurs). Mais si on travaille sur le revenu, et pas le revenu pour un temps plein, les chiffres sont un peu différents,

> Xinf=salinf[base$TRNNETO+1]/12
> Xsup=salsup[base$TRNNETO+1]/12
> quantile(Xinf,.5,na.rm=TRUE)
 50% 
1333.333 
> quantile(Xsup,.5,na.rm=TRUE)
 50% 
1500

Pour revenir à nos moutons, si on s’intéresse à la proportion des français qui touchent moins de 1500 € par mois, j’obtiens

> mean(Xinf<=1500)
[1] 0.6118895
> mean(Xsup<=1500)
[1] 0.5236287

C’est à dire qu’à la louche, je pense que 55% des français touchent moins de 1500 € par mois. A condition que les chiffres soient en net. Si la citation parlait du brut, ça correspond en gros à du 1200 € par mois, et on a alors les proportions suivantes,

> mean(Xinf<=1500*.8)
[1] 0.4215607
> mean(Xsup<=1500*.8)
[1] 0.3143752

D’ailleurs un petite remarque sur ces histoires de revenu et de temps de travail. Si le temps travaillé par mois et le salaire horaire étaient des variables indépendantes, j’aurais pu faire des calculs rapides à partir de la distribution du temps de travail. Malheureusement, ces deux données ne sont pas indépendantes. Pour s’en convaincre, calculons le produit des espérances (salaire moyen horaire par temps moyen travaillé par mois) et l’espérance du produit (revenu moyen)

> Xinf=salinf[base$TRNNETO+1]/base$NBHEUR_TOT*1820/12
> Xsup=salsup[base$TRNNETO+1]/base$NBHEUR_TOT*1820/12
> mean(Xinf[Xinf<1000000000000],na.rm=TRUE)*
+ mean(base$NBHEUR_TOT/1820)
[1] 1609.669
> mean(salinf[base$TRNNETO+1])
[1] 16832.14

qui ne sont manifestement pas égaux. On a même un dépendance positive entre le salaire horaire et le temps de travail mensuel.  On retrouve ici le fait que les personnes qui effectuent les travaux avec les salaires horaires les plus bas sont aussi celles qui travaillent le moins. D’un point de vue technique, pour étudier les inégalités temps de travail/revenu horaire, c’est un peu compliqué à cause du produit. C’est pour ça qu’avec Stéphane Mussard on avait introduit la décomposition  des mesures d’inégalités pour des inégalités multiplicatives, ici. Mais on sort du cadre du papier, et ma réponse serait que non, je doute que la personne soit de bonne foi (la seule source que j’ai trouvé sur internet donnant ces chiffres était ici).

Comments on probabilities

The only thing I remember from courses I had in probability a few years ago is that we also have to clearly defined the event we want to calculate the probability. On the Freakonomics blog, last week, the Israeli lottery was mentioned (here, see also there where I mentioned that, and odds facts from the French lottery),

Yesterday, Andrew Gelman claimed (here) that there was a probability error… Well, since Andrew is really a statistician (and a good one… while I am barely an economist), I tried to do the maths…. and to understand where the error was coming from…

Since 6 numbers are drawn out of a pool of numbers from 1 to 37, the total number of combination at each lottery is

https://perso.univ-rennes1.fr/arthur.charpentier/latex/loto01.png
> (n=choose(37,6))
[1] 2324784

Over 8 lotteries (since there are two draws per week, we can assume there 8 draws per month), the probability of no identical draws is

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

Here is the R code for those who want to check, again,

> prod(n-0:7)/n^8
[1] 0.999988

Each month, the probability of “coincidence” (I define “coincidence” the event “over 8 draws, at least two times, we obtained the same 6-uplet” or more precisely (as mentioned here) “over one calendar month, at least two times, we obtained the same 6-uplet“) is

> (p=1-(prod(n-0:7)/n^8))
[1] 1.204407e-05

The occurrence of a coincidence each month as a Geometric distribution, with probability p. And it is classical, following Gumbel’s definition (here), to consider 1/p, called the “return period“, i.e. the number of months we have to wait until we observe a coincidence (i.e. a repetition in the same month), since for a geometric distribution

https://perso.univ-rennes1.fr/arthur.charpentier/latex/loto03.png
> 1/p/(12)
[1] 6919.034

Here, the (expected) return period is 6919 years.

From my point of view, this is “the incident of six numbers repeating themselves within a calendar month”, and this is an event of once in 6919.034 years. On the other hand the median of a geometric distribution is

https://perso.univ-rennes1.fr/arthur.charpentier/latex/loto04.png
> -log(2)/log(1-p)/(12)
[1] 4795.88

which means that we have 50% chance to get such a coincidence over 4796 years.

Of course, if instead of looking at a longer period, say 100 draws, i.e. one year (here I define “coincidence” the event “over 100 draws, at least two times, we obtained the same 6-uplet“), we have in red the expected return period, and in blue the median of the geometric distribution,

> M=E=rep(NA,100)
> for(i in 2:100){
+ p=1-exp((sum(log(n-0:(i-1)))-i*log(n)))
+ E[i]=1/p/(100/i)
+ M[i]=-log(2)/log(1-p)/(100/i)
+ }
> plot(1:100,E,ylim=c(0,10000),type="l",col="red",lwd=2)
> lines(1:100,M,col="blue",lwd=2)
> abline(v=8,lty=2)
> points(8,E[8],pch=19,col="red")
> points(8,M[8],pch=19,col="blue")

or below of a log-scaled version

As Xi’an did (here), assume now that there is a lottery over 100 countries. Here I define “coincidence” the event “over k lottery draws over 100 around the world, at least two times, we obtained the same 6-uplet“, and then the previous graph becomes (with on the x axis the level of k)

Here I have a 12% chance if we consider probability to have identical numbers over a month…

But here, we can have one 6-uplet in Israel, and the other one in Egypt, say… If we want to get the same 6-uplet in the same country, the graph is now

i.e. each month there is a chance over one thousand…

> i=8
> p=1-exp((sum(log(n-0:(i-1)))-i*log(n)))
> 1-(1-p)^100
[1] 0.001203689

Note: actually, Xi’an mentioned that the probability that this coincidence [of two identical draws over 188 draws] occurred in at least one out of 100 lotteries (there are hundreds of similar lotteries across the World) is 53%! And I got the same,

> 1-(1-P[188])^100
[1] 0.5305219

Lecture notes on risk and insurance

I just finished some lectures notes on risk and insurance. The notes, that can be downloaded [pdf], are in French, and will be used at the JES (Journées d’Etudes Statistiques), organised at the CIRM (mentioned here). Previous notes on risk measures [pdf] and copulas [pdf]. Again, all comments are welcome…

Names of villages, in France

Keith Briggs published a post here on names of English place name element distribution, which contains almost twenty maps like the one where names ends by -bourn,bourne,burn (here) or -head (there). Actually, it is possible (Robin mentioned that already here) to do similar things in France… Consider the dataset containing the 35,250 commune names (here), it is an xls file containing the official name, the latitude, and on the longitude. To start with something simple, it is possible also to look at village containing “saint” in it

There are a lot, and there is no obvious geographic trend. For some simple geographic trend, t is possible to see where are villages having a name ending with “sur mer” (meaning literally “on the sea”) below on the left. Obviously, we cannot find such places in the Alps. Similarly for names ending with “Seine” they are clearly on the Seine river, on the right

> ville=read.table("D:\\r-data\\ville.csv",sep=";",header=TRUE)
> nrow(ville)
[1] 35376
> ville$maj=as.character(ville$Nom.Ville)
> n=nchar(ville$maj)
> I=substr(ville$maj,pmax(0,n-8),n)
> Ind=I=="-sur-Mer "
> sum(Ind)
[1] 98
> library(maps)
> map('france', fill = FALSE)
> X=ville[Ind,]
> x=as.numeric(as.character(X$Longitude))
> y=as.numeric(as.character(X$Latitude))
> points(x,y,pch=19,col="blue",cex=.6)

In order to continue with some geographic pattern, consider the end of the names, such as “-gny” (below on the left, in red) or “-ac” (below on the right, in blue)

Some pretend that “-ac” comes from Gaelic, and can be found in Celtic regions (here in Brittany). Obviously, there is also an origin in Occitany (south west of France). And this gave also in Oïl region the “-gny” (in North and North-East regions). Consider similarly end of the names, such as “-an” (below on the left, in red) or “-ey” (below on the right, in blue)

Still about the end of the names, it is also possible to look for village ending either with “-a” (below on the left, in red) or “-o” (below on the right, in blue)

We are now in the southern part of France…. “-a” in Corse and Pyrénées, while “-o”can be found in Corse, and in Brittany. For the beginning “ker-” or “lan-”  (below on the left, in red) or “castel-” (below on the right, in blue),

“ker-” appears in 18,000 location names (as mentioned here) but only in some village names. It is similar to “castel-” in the southern part of France.
To go a bit further, 40 years ago, Georges Brassens sang a song entitled “La ballade des gens qui sont nés quelque part“.

He says that people are usually extremely proud of their villages…. Actually, their are more people proud of living over something than under something: below are villages containing “sous” (i.e. under below on the left, in red) or “sur” (i.e. over below on the right, in blue)

On the other hand, villages containing “grand” or “grande” (i.e. tall or big below on the left, in red) or “petit” or “petite” (i.e. small below on the right, in blue) seem to be correlated: close to a city with “grand” there is a village with “petit” in it. For instance Virieu-le-Grand and Virieu-le-Petit, or  Essigny-le-Grand and Essigny-le-Petit.

And finally, I found surprising to see so many village containing “montagne” (ie mountain below on the left, in orange) or starting with “Mont” (below on the right, in purple) that are far from mountains,

You do not need to live close to some mountains to get mountains in it. Even in Brittany you can find dozen of villages starting with “Mont”….

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

Extracting information from a keyboard…

Yesterday, Baptiste published a post on “ethno-photography” (here). As he mentioned it, in Paris 8, they experience a real absence of serious cleaning of office equipment. He then shows the keyboard of the only computer they can use in the sociology department (for forty researchers),

Apart from the fact that everyone in France should be ashamed to see how much is spent in universities (which is the first information we have from that picture), we should also be able to guess in which langage people work in this department.
I considered three books (two in French, one in English) and I would like to see the frequency of each letter,

  • Mauss, manuel d’éthnographie (here), 1926
  • Durkheim, Livre II: Les croyances élémentaires in Les formes élémentaires de la vie religieuse (here), 1912
  • Ferri, Criminal Sociology (here), 1896

Those three books are in rich text format, I just changed it to get text files… Then, it is easy to count appearance of letters. E.g. for Mauss,

> library(corpora)
> textfile=scan("MAUSS-manuel.txt",
+ what="char", sep="\n")
Read 1550 items
> textfile<-tolower(textfile)
> M=NA
> for(i in 1:length(textfile)){
+ line=textfile[i]
+ M=c(M,strsplit(as.character(line),"")[[1]])
+ }
> T=table(M)
> T
M
    '     -           \t     !     "     %     &     (     )     ,     .     / 
 5308  1049 86589    44     3     3     2     2   370   391  6609  4909    12 
    :     ;     ?     @     ]     _     ~     ’     =     «     »     ¬     ° 
  819  1178   113     1     1     4     1    39     1   108   107   823     3 
    …     0     1     2     3     4     5     6     7     8     9     a     à 
    1    69   213    83    73    34    48    33    28    64   151 30559  1651 
    â     ä     b     c     ç     d     e     é     è     ê     ë     f     g 
  224     3  3562 14678   110 17713 63955 10354  1798  1000     5  4555  4911 
    h     i     î     ï     j     k     l     m     n     ñ     o     º     ô 
 4359 30851   226    47  1147   247 24792 12844 32525     6 25562     2   151 
    ö     œ     p     q     r     s     t     u     ù     û     ü     v     w 
   12    52 12696  4667 28237 37630 32945 25001   211    40     9  4787   164 
    x     y     z 
 1996  1222   343

Then, we can summarize in to see proportion of standard 26 letters, and we have, for Mauss,

and for Durkeim,

If we compare the two, we have almost the same proportions,

If we look at our book in English now, we have

i.e., if we compare with Mauss for instance

So we have much more E in French than in English, but still, people writing in English use a lot the E. So looking at the E should not give us any clue…. But we can see that in English, the H is as common as the L, or the C. Not in French, where L is much more frequent than the H. But on the picture, the C is more clear than the H. We can also look at the U, which is common in French, not in English… Here, on the keyboard, it is perfectly clear… so I guess people use it frequently.
So I would say that they write more in French than the write in English, on that computer.
Actually, the same idea has been used a long time ago on calculators to see that Benford’s law works: some numbers are really used (as well as the legend pretends that some pages in logarithm books were never used….), see here orthere. So Baptiste, if one day the keyboard is cleaned up, please send me another picture after a few weeks to see if things have changed….
An for those who cannot imagine how it is to work in some universities in France, just look at his blog (here). Pictures are unbelievable….Good luck Baptiste….

Lottery, and martingales

I recently got a comment on a post I published one year ago, here, about the fact that in September 2009, on the 6th and the 10th, the 6 same numbers came out at the lottery, in Bulgaria (but  I do not understand the question: the author of the comment ask about the order the numbers came out…)
Xi’an published also a post on that topic, there, since last week, the same thing happened in Israel.
All that reminded me a discussion I had with a colleague about another post (here) where I mentioned that I found a strange distribution of numbers in the French lottery (the old one actually). For those who want to check, all historical events are here, in a zip file. My colleague was wondering if I found the martingale to win the lottery…

First, I do not like that term, since martingale is something different from a mathematical point of view… Second, let us look if it would have been possible to make some money… (free lunch ?)

> loto=read.table("D:\\loto.csv",dec=",",header=TRUE,sep=";")
> ntirage=nrow(loto)
> loto=loto[51:ntirage,]
> ntirage=nrow(loto)
>   N=as.matrix(loto[,c("boule_1","boule_2","boule_3","boule_4","boule_5","boule_6")])
> n=as.vector(N)
> length(n)
[1] 28848
> (TN=table(n))
n
1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20
607 576 571 618 579 598 608 582 588 590 562 577 577 580 591 630 558 567 594 608
21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40
578 562 579 583 574 589 602 572 550 598 604 582 545 646 597 618 599 636 609 588
41  42  43  44  45  46  47  48  49
576 589 577 585 618 596 560 571 604

So, it might look nice, but we have to compare that distribution with the one we should have with “independent” draws. It is not possible to look at a discrete uniform distribution: the six numbers are not independent. Each day, the 49 balls are back in the urn, but within a day, we do not have independent draws (it is a sample without replacement of balls). Hence, with 4808 lottery draws, each number cannot be obtained more than 4808 times. So, let us use monte carlo techniques to  look at the theoretical distribution,

> M=matrix(NA,49,1000)
> for(s in 1:1000){
+ B=NA
+ for(i in 1:ntirage){B=c(B,sample(1:49,size=6,replace=FALSE))}
+ B=B[-1]
+ M[,s]=sort(table(B))
+ }
> q50=function(x){quantile(x,.5)}
> Q50=apply(M,1,q50)
> lines(1:49,Q50,col="red",lwd=2)
> q10=function(x){quantile(x,.1)}
> Q10=apply(M,1,q10)
> q90=function(x){quantile(x,.9)}
> Q90=apply(M,1,q90)
> polygon(c(1:49,49:1),c(Q10,rev(Q90)),col="light blue",border=NA)
> lines(1:49,Q10,col="red",lty=2)
> lines(1:49,Q90,col="red",lty=2)
> lines(1:49,Q50,col="red",lwd=2)
> points(1:49,sort(TN),pch=19,type="b")

Looking at the graph, it looks like some numbers appeared too frequently, especially the ones that did not appear frequently (bottom left). So, since I have removed the last 50 draws, let us see if we could have used that information, somehow…

> nb=names(sort(TN))
> loto=read.table("D:\\loto.csv",dec=",",header=TRUE,sep=";")
> loto=loto[1:50,]
> N=as.matrix(loto[,c("boule_1","boule_2","boule_3","boule_4","boule_5","boule_6")])
> n=as.vector(N)
> TN=table(n)
> TN[nb]
> barplot(TN[nb])

Unfortunately, numbers that came out too frequently over 4800 draws did not appear that frequently of the last 50. Playing top number might not have been a great strategy.

(numbers that came out frequently are on the right, while those we did not see much are on the left)… What about worst numbers: if I had decided to play the 6 that did not come out very frequently (we’ve seen earlier that they should have appeared even less, actually), would it have been interesting ? As we can see, our top 2 numbers were numbers that did not appear frequently earlier (29 and 47 appears respectively 10 and 11 times over 50 draws)….
Over 50 draws of 6 balls, the expected frequency of 6 given number is around 36.7,..

> S=rep(NA,10000)
> for(s in 1:10000){
+ B=NA
+ for(i in 1:50){B=c(B,sample(1:49,size=6,replace=FALSE))}
+ B=B[-1]
+ S[s]=sum(B%in%(1:6))
+ }
> mean(S)
[1] 36.7694

But here for the top 6, we have

> z=TN[nb]
> sum(rev(z)[1:6])
[1] 29

i.e. the top 6 appeared 29 times over 50 draw of 6 balls (which looks low) and for the worst 6, it is a bit higher,

>  sum(z[1:6])
[1] 38

If we look at the theoretical density of the frequency of 6 given number, we have

i.e. our worst 6 is a nice average (in green) while top 6 did not appear frequently this time (here in blue) ! So we could not have used that information….
Anyway, if some of you are interesting using statistics to get a free lunch, with the nouveau loto, I did not see any strange pattern (data can be downloaded here).

I am terribly sorry, but I cannot help anyone winning at the French Lottery….

Prétentieux, va !

Depuis quelques semaines, mes billets sont presque exclusivement en anglais… En fait, même si je fais cours en français, et que la plupart de mes collègues ici parlent français, je me rends compte qu’on parle naturellement anglais…. en tous les cas plus naturellement qu’en France. Les séminaires sont toujours en anglais, dès qu’on discute à plus de 3, la probabilité qu’une personne ne parle pas français est tellement grande qu’on parle le plus souvent en anglais….

Finalement, comme toujours quand je ne suis pas en France, je trouve naturel d’écrire en anglais… (ok, d’essayer d’écrire quelque chose qui ressemble à de l’anglais). Mais hier soir, ma femme m’a fait remarqué que c’était prétentieux d’écrire en anglais… Et effectivement, quand je traîne sur les blogs francophones, presque tous sont exclusivement en français…. Bref, je me remets à douter.
J’ai vu l’autre jour que David Monniaux s’en prenait sur son blog à un édito de Jacques Julliard (ici), et d’ailleurs David publie de plus en plus en anglais ces derniers temps (en tous les cas dès que l’on quitte les sujets franco-français). J’ai vu aussi que de plus en plus de doctorants blogueurs se lancent en anglais…
Le débat reste ouvert, et je suis preneur de toutes les remarques….

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…

Lecture notes on risk measures

I just finished some lectures notes on risk measures. The notes, that can be downloaded [pdf], are in French, and will be used at the JES (Journées d’Etudes Statistiques), organised at the CIRM (mentioned here). The short course will follow a short introductionary course by Jean Marc Tallon. He recently published on his website a survey on decision theory under uncertainty (here) which is almost what he will talk about in Marseille. All comments are welcome…

"sendo l'intento mio scrivere cosa utile a chi la intende…"