[Geekery] Le mariage gay, une “loi détruisant notre société” ?

Comme certains ont pu le dire à l’époque des débats sur le mariage gay, en 2012, il y avait une crainte que cette nouvelle loi amène l’apocalypse et la fin de la civilisation en France. Il semblerait que, deux ans et quelques jours après, l’apocalypse ne soit pas encore arrivée, et que seuls quelques dizaines de milliers de mariage aient été célébrés. Mais on n’en est qu’à deux ans après l’adoption de la loi : que pourrait-il bien se passer à plus long terme ?

Pour éviter de recevoir trop de mails d’insultes et autres réjouissances, je vais vous rassurer de suite : il ne s’agira pas dans cet article de s’intéresser à l’impact sur les enfants, la société, ou que sais-je, du fait d’avoir autorisé deux hommes ou deux femmes à se marier. Non, ce qui m’a frappé dans cette loi, c’est son article 11. Celui-ci dit :

« En cas de désaccord entre les parents, signalé par l’un d’eux à l’officier de l’état civil, au plus tard au jour de la déclaration de naissance ou après la naissance, lors de l’établissement simultané de la filiation, l’enfant prend leurs deux noms, dans la limite du premier nom de famille pour chacun d’eux, accolés selon l’ordre alphabétique. »

Je résume (sans aucune mauvaise foi) : blablabla, le nom des enfants sera celui des parents accolés dans l’ordre alphabétique. C’est bien, l’ordre alphabétique, non ? C’est neutre, c’est pour l’égalité homme-femme, homme-homme et femme-femme, tout ça. Mais ce qui me perturbe, personnellement, c’est que si monsieur A épouse madame B (+1 point hétérocentrisme), leurs enfants s’appelleront A-B. Et si ils épousent les enfants de C et D, soient les C-D, eh bien leurs enfants à eux (les petits-enfants de monsieur A et madame B) auront A-C comme patronyme. Vous avez bien vu, les noms de famille B et D ont disparu. Pourquoi ? Parce qu’ils sont plus loin dans l’ordre alphabétique. Le mariage va donc bien détruire notre civilisation, en s’attaquant à nos patronymes ! Peut-on essayer d’évaluer la vitesse de cette décadence ?

Quelle sont les noms de famille les plus courants ?

La question du nom de famille le plus courant n’est pas compliquée : c’est Martin, avec une très large avance. Certains sites comme celui-ci permettent de le savoir, bien que leurs sources ne soient pas très clairement précisées. Mais ce qui nous intéresse ici est une distribution plus générale des noms de famille en France. Pour cela, nous allons utiliser les données du Bac 2014, en faisant des statistiques sur les noms des candidats. Cela peut créer quelques biais, parce que cela concerne une seule génération et non pas la population toute entière, et parce que toute la génération en question ne se présente pas au bac (qu’il soit général, professionnel ou technologique). Néanmoins, c’est une des sources les plus solides pour ce genre de questions, d’autant plus que leur âge uniforme (et leur absence d’enfants, a priori) va permettre de les considérer comme la première génération. La distribution, rapportée à 10.000 personnes, est la suivante :

distr_init

On a représenté la distribution des noms de famille suivant les initiales. On remarque déjà que c’est la lettre B qui est la plus courante (vous connaissez beaucoup de personnes dont le nom de famille commence par B ? C’est normal), et qu’il y a de fortes disparités entre les différentes lettres : pas grand monde n’a un patronyme commençant par Q, U ou X, du moins en France.

Avec un peu de mauvaise foi

Maintenant que nous connaissons la situation initiale de la répartition des noms de famille en France, du moins de leurs initiales, il est temps de s’intéresser à l’impact qu’aurait le fait d’attribuer aux enfants les deux noms de leurs parents 1 et 2 (-1 point hétérocentrisme), ordonnés dans l’ordre alphabétique, dans la limite d’un nom par parent. Pour cela nous allons utiliser un modèle très simple de générations : parmi les N individus de notre population, les N/2 premiers vont se marier avec les N/2 autres, au hasard (du moins au hasard selon l’initiale du nom de famille, d’autres facteurs pouvant jouer). Ensuite, on va supposer que chacun des couples aura deux enfants : leurs noms de familles commenceront par l’initiale du nom des deux parents qui est le premier dans l’ordre alphabétique. Evidemment, comme le mariage gay est autorisé, on ne s’intéressera pas au sexe des marié(e)s : ils auront leurs enfants de la façon qui leur convient, cela n’est pas la question soulevée ici.

À quoi va ressembler la distribution des noms de familles dans la génération suivante, pour les bacheliers de 2040 ? Les simulations donnent le résultat suivant :

distr_g1

On remarque déjà une belle modification de la distribution des noms de famille. Beaucoup plus de personnes ont un nom qui commence par une lettre de début de l’alphabet, et c’est évidemment l’inverse pour les noms dans le milieu et la fin de l’alphabet. Mais la modification est encore mineure : que se passera t-il dans 5 générations, soit dans 125 ans environ ?

distr_g5

Impressionnant. La distribution n’a plus rien à voir : tous les patronymes commencent par A, B, voire C pour quelques chanceux. La totalité de la culture présente dans nos patronymes a été anéantie. C’est véritablement l’apocalypse annoncée.

Soyons plus réalistes

J’ai un peu joué sur les mots dans la partie précédente. L’article de loi ne force personne à donner comme nom de famille à ses enfants les deux noms de ses parents dans l’ordre alphabétique, elle donne seulement un cadre pour toute situation où il y aurait désaccord ou problème. Il s’agit donc d’appliquer cette règle uniquement dans ces cas. Nous évaluons à environ 10% des cas cette règle : que ce soit par erreur, par désaccord, pour éviter tout problème dans le couple, il nous semble crédible qu’un couple sur 10 environ suive la règle édictée par la loi. Pour les 90% autres, il y a un accord entre les parents : on peut raisonnablement penser que cet accord se fait la moitié du temps en faveur du parent 1, et la moitié du temps en faveur du parent 2 (-1 point hétérocentrisme, j’espère finir en négatif). Dans ce cas, l’évolution d’une génération à l’autre est évidemment moins marquée que dans le cas précédent. Mais que se passe t-il à long terme ?

prop_AB_Reste

Ce (magnifique) graphe s’intéresse à l’évolution générale des initiales des noms de familles, entre ceux commençant par A, qui commencent à environ 4%, pour monter très haut génération après génération, ceux par B qui connaissent un pic aux alentours de la 25ème génération, tandis que toutes les autres lettres tendent à disparaître à plus ou moins long terme. Plus précisément, ce graphe dynamique montre l’évolution d’une génération à l’autre de la distribution des noms de famille :

Même avec des hypothèses plus réalistes, il n’y a que le délai avant l’apocalypse qui change. De manière plus précise, nous pouvons estimer à partir de combien de générations il n’y aura plus de noms commençant par chacune des lettres (sauf A, qui ne disparaît jamais, et B qui persiste de nombreuses générations) :

C D E F G H I J K L M N
100 98 74 82 87 77 57 71 70 81 81 63
O P Q R S T U V W X Y Z
56 72 45 70 70 64 36 62 48 13 39 46

Légende : Sous chaque lettre est indiquée le nombre de générations nécessaires avant la disparition de la totalité des noms de famille commençant par cette lettre. Par exemple, il faut attendre 46 générations pour que tous les patronymes commençant par Z disparaissent.

Les noms en X sont les premiers à disparaître, après 13 générations, soit un peu plus de 3 siècles, ce qui nous laisse tout de même le temps de voir venir. Mais, en parlant de noms de famille commençant par X, vous en connaissez beaucoup ? Xavier, d’accord, mais à part celui-là ? La plupart d’entre eux sont des noms d’origine étrangère, et donc viennent de l’immigration. Ce qui est précisément le facteur que nous avons oublié.

L’immigration va-t-elle nous sauver ?

Il serait plus correct de dire l’un des deux facteurs : pour être tout à fait exact, il faudrait mentionner que pas tout le monde ne se marie et n’a des enfants, mais ce comportement n’a aucune raison d’être lié au nom de famille, du moins pas directement. On peut donc supposer que 15% des gens n’ont pas d’enfants (selon le Figaro, c’est un homme sur 5, mais à 50 ans ce qui laisse un peu de marge – on peut aussi aller jeter un oeil sur fivethirtyeight (en) sur la question aux États-Unis).

Mais revenons à l’immigration. Comme nous n’avons aucune envie de nous lancer dans le débat sur le nombre d’immigrés arrivant ou repartant, laissons ça aux spécialistes. Ce qui nous intéresse ici est une étude à long terme, et un modèle fréquemment utilisé en démographie pour cela est celui de population stationnaire, où la taille n’évolue plus. On va donc supposer que les 15% d’enfants qui ne sont pas nés d’après l’hypothèse précédente sont remplacés par une immigration entrante de même taille.

Reste une question de taille : si l’on s’intéresse aux noms de famille des personnes venant vivre en France, comment connaître leur répartition ? On pourrait imaginer récupérer des données sur l’ensemble des pays du monde, et les pondérer selon le pourcentage d’arrivants par pays, qui n’est pas très bien connu, en prenant en compte des évolutions pour chacune des générations successives. Mais c’est très compliqué, et cela n’aurait concrètement rien apporté de tangible. Nous allons nous limiter à deux hypothèses qui sont aussi peu crédibles l’une que l’autre, mais qui ont le bénéfice d’être simples à mettre en oeuvre.

Tout d’abord, supposons que la totalité des pays du monde a la même distribution des noms de famille que la France. Dans ce cas, les noms des immigrants suivent la même répartition que celle de la France avant l’application de la loi, c’est à dire celle des bacheliers 2014. L’évolution est alors la suivante :

imm_FR

Une autre hypothèse est de dire que comme les noms de familles peuvent provenir de pays, de cultures et même d’alphabets très différents, il y a autant de chance que le nom de famille commence par chacune des lettres de l’alphabet latin. Autrement dit, la répartition des initiales des noms des arrivants est uniforme. L’évolution est alors la suivante :

imm_unif

Dans les deux cas, la situation ne semble plus se dégrader à moyen terme : les noms de familles commençant par toutes les lettres hormis les toutes les premières ne vont pas disparaître d’ici quelques siècles, bien qu’il y ait une légère hausse des noms commençant par A et B dans le cas d’une répartition des noms mondiale similaire à celle de la France. Cas très improbable, on le rappelle.

Il semblerait donc que l’immigration nous sauvera de l’apocalypse causée par l’adoption du mariage gay. Quel retournement de situation, n’est-ce pas ?

image source : Wikimedia, Le mariage de Sigebert Ier et de Brunehaut

[Sampling] L’algorithme CURIOS pour les nuls

La problématique de la baisse continue des taux de réponse aux enquêtes depuis les années 1950 conduit les méthodologues d’enquête à innover et à adapter les techniques d’enquête aux nouvelles données (voir par exemple la présentation de Carl-Erik Särndal aux JMS de 2012). Certaines de ces méthodes sont basées sur les techniques de priorisation, qui consistent à “prioriser” (relancer par téléphone, sur-représenter dans un échantillon, encourager par des cadeaux…) des individus parmi la population afin d’équilibrer les échantillons, c’est à dire d’éviter de se retrouver avec un échantillon uniquement composé de ceux qui auront eu la “gentilesse” de répondre à l’enquête spontanément, groupe qui ne présente pas nécessairement la diversité voulue pour l’analyse des résultats de l’enquête. Ces techniques de priorisation sont principalement utilisées dans le cadre d’enquêtes par téléphone (dites CATI). Nous souhaitons ici présenter une méthode applicable aux enquêtes CAPI, c’est à dire en face-à-face, avec un enquêteur se rendant au domicile de la personne interrogée, ce qui suppose une phase de reconnaissance et ne permet pas une priorisation à la volée.

L’algorithme CURIOS (Curios Uses Representativity Indicators to Optimize Samples) est utilisable dans le cadre d’une enquête en plusieurs vagues. À la fin de la première vague, on peut trouver les groupes à prioriser, et l’algorithme permet de tirer un échantillon de seconde vague qui réponde à des objectifs précis d’équivalence, de qualité et de spécificité. Expliquons tout cela plus en détail.

Cet article fait suite à la présentation du 02 avril 2015 aux JMS (Journées de Méthodologie Statistique) à Paris, et en particulier à la présentation (lien à venir) liée aux deux articles suivants : Algorithme Curios et méthode de ‘priorisation’ pour les enquêtes en face à face. Application à l’enquête Patrimoine 2014 et l’utilisation des R-indicateurs pour « prioriser » la collecte des enquêtes Ménages : une application à l’enquête Patrimoine 2010.

L’échantillonnage usuel

On suppose que l’on s’intéresse à une population de 100 individus, représentés comme suit :

population

Usuellement, la plupart des enquêtes sont réalisées en une seule phase : on tire un échantillon d’une taille convenable (c’est à dire suffisante pour permettre une bonne précision des résultats, mais pas trop importante pour limiter les coûts et la logistique de l’opération) selon un plan de sondage adapté et on réalise l’enquête auprès des individus sélectionnés. Par exemple, on peut obtenir l’échantillon suivant :

echt

Détection des groupes priorisables

Une fois que ces individus ont été enquêtés, certains auront répondu et d’autres non. On dispose à propos de chacun d’entre eux d’informations socio-démographiques présentes dans la base de sondage qui permettent usuellement d’étudier et de traiter le phénomène de non-réponse. Ici, on souhaite identifier les groupes qui sont sous-représentés parmi les répondants, c’est à dire des groupes dont le comportement de réponse diffère à la baisse. Pour cela, on utilise des indicateurs de représentativité tels que les R-indicateurs (Schouten, 2009) ou d’autres : variance due à la non-réponse, etc. L’échantillon peut donc être séparé en deux groupes :

select

où les individus entourés de bleu sont non-répondants, et les individus rouges sont ceux qui sont à prioriser, si une deuxième vague de l’enquête devrait avoir lieu. Pourquoi tous les individus non-répondants ne sont pas rouges ? Par exemple, on peut imaginer que les individus rouges sont ceux vivant en milieu rural et les noirs en milieu urbain, ou plus spécifiquement que les rouges sont les individus vivant en milieu rural dans des logements de plus de 70m² : la séparation entre rouges et noirs doit reposer sur des critères disponibles dans la base de sondage. Ces critères sont déterminés à l’aide des indicateurs de répresentativité.

Que fait CURIOS alors ?

Comme nous l’avons déjà mentionné, l’algorithme CURIOS nécessite que le sondage soit réalisé en deux phases. L’échantillon classique ne conviendrait pas. Il faut donc tirer dans une première phase une partie seulement des individus. Pour que cela soit plus visuel sur l’exemple, nous en tirons la moitié, mais il est plus recommandé de réaliser environ 75% de la collecte avant de lancer une deuxième vague. Concrètement, à la fin de la première phase on peut réaliser la procédure de priorisation et détecter les individus rouges :

select

On remarque bien ici que cette recommandation de 75% de la collecte effectuée en première phase permet de s’assurer que les individus que l’on identifiera comme sous-représentés le sont bien structurellement et non pas seulement à cause d’un artefact statistique. On peut alors appliquer la séparation en deux groupes, noirs et rouges, à la population entière (moins le premier échantillon, évidemment) :

select

On a bien remarqué que les rouges étaient moins bons répondants, et qu’il fallait les prioriser. La logique veut donc que dans l’échantillon de seconde vague le nombre de rouges soit plus important que prévu. Mais combien de rouges faut-il exactement ? L’algorithme CURIOS permet de quantifier cela précisement. On obtient par exemple :

select

Soit, en combinant les échantillons de première et de deuxième vague :

select

Cela peut impliquer que le nombre de répondants est plus faible que dans le cas de l’échantillon standard : il y a plus de rouges, et ils répondent moins bien. Mais cela n’est pas un problème car la qualité de l’échantillon est assurée, et l’exploitation des résultats sera possible.

Quelques simulations

Pour s’en convaincre, nous avons réalisé des simulations sur une enquête INSEE réalisée en 2010 sur le patrimoine des français. Les variables d’intérêt de l’enquête sont le patrimoine moyen des français (dit brut), le patrimoine net qui correspond au patrimoine brut moins les dettes, et la ventilation en patrimoine financier (comptes…), immobilier et professionnel (pour les entrepeneurs et commerçants). On veut comparer la précision de l’estimation de ces différents indicateurs dans l’enquête qui a eu lieu en 2010 en présence ou non d’une procédure de priorisation.

Pour cela, nous avons simulé un phénomène de non-réponse supplémentaire, et décidé d’affecter certains enquêteurs à la réalisation de quelques enquêtes en plus dans les zones les plus touchées. Dans un premier temps, cette affectation était faite aléatoirement : les enquêteurs pouvaient aller voir n’importe quel logement de la zone, aucun critère n’était donné. Dans un second temps, cette affectation était limitée aux ménages sous-représentés : les enquêteurs ne pouvaient enquêter que ces logements et n’avaient pas la possibilité d’en interroger d’autres. Les résultats obtenus sont les suivants :

simus

On remarque bien que pour la plupart des indicateurs (surtout les deux principaux, patrimoine brut et net) la méthode de priorisation semble améliorer la précision des estimations, de manière bien plus importante que la méthode aléatoire. Il est important de noter que cette amélioration de la précision se fait à nombre de fiches-adresses enquêtées constant voire inférieur, et donc à nombre de répondants inférieur en raison de la concentration des enquêteurs sur des populations par nature moins bons répondants.

[Geekery] How many stations of the Paris Métro could you pass through in alphabetical order?

I spend a lot of time in the Parisian subway, as I use it everyday to go to work, and sometimes to meet some friends of mine. I usually spend my time listening to music and asking myself lots of questions about its organisation : how the time frequencies of subways are decided ? How can empty trains become completly full within minutes ? I guess there are people working on that, but there are some sillier questions I doubt anyone is paid to look for an answer to.

One of these came recently into my mind when I was on the “ligne 13”, between “Montparnasse Bienvenue” and “Porte de Vanves”. I realized that I passed through 4 stations which were in alphabetical order.

So now I wonder : using only the subway, what is the maximal number of stations in alphabetical order you could go through? I already know it is at least 4 because I do that every morning or so, but could we go up to 5, 6, 7, or maybe more ? Let’s do some research!

map

5 stations in a row…

I’m a very lazy person, so I’m not going to do it by looking at all the lines in the Paris Metro, but I’d rather ask gently my computer to do it for me. First, I need to gather some data about the stations’ name on every line of the metro. These informations are available at different places, such as Wikipedia and “data.gouv”, which is the official French governement website on open data.

All left to do is an easy computation of the maximum number of stations in alphabetical order for each line, going back and forth. The algorithm used is thereby (in pseudocode):

maximum_number <- function(Line) {
    for Station in Line {
    if (Station > Station_prec) {Counter <- Counter + 1}
    else {Maximum = max(Maximum,Counter) ; Counter = 1}	
     }
  Maximum = max(Maximum,Counter);
  return(Maximum)
}

The results obtained are compiled in this graph:

graph1

The answer to my question seems to be 5 stations: "ligne 2" (Belleville / Couronnes / Ménilmontant / Père Lachaise /Philippe Auguste), "ligne 5" (Bobigny-Pablo Picasso / Bobigny-Pantin-Raymond Queneau / Église de Pantin / Hoche / Porte de Pantin) or "ligne 12" (Falguière / Montparnasse-Bienvenüe / Notre-Dame-Des-Champs / Rennes / Sèvres-Babylone).

... or maybe more?

Ok, so I was right 4 wasn't the best I could do. But I feel like I'm forgetting something... Of course, the connections between the lines! The Paris Metro is organized around huge hubs such as "Montparnasse-Bienvenüe" which belong to plenty of lines, so I guess using these big stations might get us to better than 5, right?

As you may know, the organization of the subway is quite complex. To simplify my study, I will consider that two lines are connected at a station iff the station's name is the same on the two lines. For instance, on the map previously shown, I consider "La Chapelle" and "Gare du Nord" not to form a connection, but "Gare du Nord" to be one on its own, between Line 4 and 5.

The algorithm I use is a more complex than the previous one. It requires to browse through all the subway network in order to find the higher number of stations in alphabetical order. For those with some computing background who are used to recursive functions, here is my pseudocode :

 max_length <- function(Station) {
	Connections = connection_list(Station)
        Values = void()
	for (Line in Connections) 
		{ Next_Station = Line[id_Station + 1]
                  Values = c(Values, max_length(Next_Station) }
	return(max(Values))
}

And the result is... 6! Well, I though that using the connections might lead to an higher number of stations, but whatever. There are numerous ways to achieve the maximum number:

  • Belleville (2) / Colonel Fabien (2) / Jaurès (2->5) / Laumière (5) / Ourcq (5) / Porte de Pantin (5)
  • Duroc (13) / Montparnasse Bienvenue (13->12) / Notre-Dame des-Champs (12) / Rennes (12) / Sèvres Babylone (12->10) / Vaneau (10)
  • Gaité (13) / Montparnasse Bienvenue (13->12) / Notre-Dame des-Champs (12) / Rennes (12) / Sèvres Babylone (12->10) / Vaneau (10)
  • Edgar Quinet (6) / Montparnasse Bienvenue (6->12) / Notre-Dame des-Champs (12) / Rennes (12) / Sèvres Babylone (12->10) / Vaneau (10)
  • Falguière (12) / Montparnasse Bienvenue (12) / Notre-Dame des-Champs (12) / Rennes (12) / Sèvres Babylone (12->10) / Vaneau (10)

Well, there are at least two completly different ways to achieve the record of 6 stations in a row: the first one, in the northeastern part of Paris, and the 4 others located near Montparnasse. Fun fact: before 1942, Montparnasse and Bienvenüe were two different stations (see there), and "Porte de Pantin" station wasn't built (see there): therefore the only possible circuit was :

  • Falguière (12) / Montparnasse Bienvenue (12) / Notre-Dame des-Champs (12) / Rennes (12) / Sèvres Babylone (12->10) / Vaneau (10)

But I guess there were more urgent things to deal with at this time.

Let's simulate some graphs!

I'm now wondering what would happen if I did this kind of analysis for all the subway network all around the world. But as I mentionned before, I'm kind of lazy, so I won't try to gather all the data about all these networks. Another approach is to try to know if the Paris Metro network is particular, or if it behaves like an average graph. So I'll try to simulate some networks, compute the maximum number of ordered nodes, and compare the mean value obtained to 6.

A graph is a representation of a set of objects, denoted the vertices, which are the stations here, and the relations between them. Each time two vertices A and B are in relation, there is a link, denoted an edge, going from A to B. In our subway example, it will mean that you can go directly from station A to station B without any stop. We will consider it as a symmetric graph, which means that when you can go from A to B you can also go back to A from B, even if in the Paris Métro there is some absurd lines (for instance, the West part of the "ligne 10").

Unfortunately, simulating a subway network is quite hard (citation needed). As I only have small knowledge of graph theory, I'll try the simplest way to simulate a graph, which happens to be the Erdos-Renyi method: we first generate 300 vertices (because there is about 300 stations in Paris), numbering them from 1 to 300. Then, for each possible link between two stations, we roll a (odd) dice to determine if there is an edge or not. The dice is designed in order to obtain an average of 371 edges, because there is 371 links in the parisian network.

Using R, I simulated 1000 graphs with this method. Then for each of them I apply an algorithm, similar to the one previously presented, to get the maximum number of ordered vertices, using the label given to them. I obtain the following results:

graph2

It seems that 6 is at the bottom of the distribution of the maximum number of ordered stations (average is about 7.8), which may lead us to think that we are unlucky in Paris. But we have to keep in mind a few things :

  • First, the model used here isn't at all a good model for subway network. So this comparison is a little fallacious, and we shouldn't jump to any conclusion.
  • Someting else that may explain the difference between the simulations and Paris Métro is that the names of the stations aren't really uniformly distributed: there are some spatial correlations between the names, for instance all the "Porte de" are located near Paris' exits, and "Bobigny-Pablo Picasso" and "Bobigny-Pantin-Raymond Queneau" are obviously close stations. This might lead to less randomness and biased results.

I guess I have to do the job for all the other subways. I think I will start with this one!

Image titre : Paris Métro Entrance, Abbesses

[Games] Les mots les plus probables à “Des chiffres et des lettres”

Pour changer un peu, aujourd’hui, un post écrit en français ! De temps à autres pendant mes congés, je prends le temps de regarder la célèbre émission “Des chiffres et des Lettres”. Je m’assois devant ma télévision avec une bonne tasse de thé bien chaud, du papier et un crayon pour compter les points, et je tente de terminer avec le score le moins éloigné possible du vainqueur du jour. Sur la partie “Le compte est bon” (les chiffres, donc), je crois m’en tirer honorablement. “Le mot le plus long” (les lettres) est pour moi un exercice beaucoup plus difficile ; il n’est pas rare que je n’aie rien à proposer, quand les candidats ont de leur côté trouvé le même “8 lettres” dont j’ignore la signification !

Bien que je n’aie que trop rarement l’occasion de regarder l’émission, il me semble que l’on retrouve assez souvent les mêmes mots parmi les solutions. De plus, certains mots rares semblent connus de beaucoup de candidats : si je comprends bien, il s’agit quand on s’entraîne de procéder un peu comme au Scrabble, en apprenant à reconnaître des types de combinaisons, et les listes de mots associés. Je me suis donc demandé si l’on pourrait construire des probabilités de tirage pour chaque mot, et en déduire la liste des mots les plus probables.

En cherchant sur les interwebs, je suis tombé sur beaucoup de solveurs “Le mot le plus long” programmés en javascript, quelques fanblogs avec des conseils pour joueurs en herbe, mais rien qui ne ressemble à un début de réponse à la question que je me posais (pour une analyse de la partie “Le compte est bon”, je vous renvoie à l’excellent post sur le blog Datagenetics).

Le problème

Il s’agit donc de dresser une liste des mots valides dans le jeu, et de leur attribuer à chacun une probabilité de tirage. Intuitivement, cette probabilité de tirage doit prendre en compte la fréquence des lettres tirées par l’ordinateur de l’émission, ainsi que la fréquence d’occurrence des voyelles relativement aux consonnes.

En effet, si vous avez bien compris le principe de l’émission (sinon, rassurez-vous, vous n’êtes pas les seuls), chaque candidat choisit à tour de rôle “consonne” ou “voyelle”, et l’ordinateur tire ensuite une lettre correspondant au choix du candidat. Les choix des candidats doivent donc logiquement influencer la probabilité que l’on attribue à chaque mot.

Afin de ne pas trop surcharger ce post, je n’écrirai pas ici le code utilisé pour chaque étape. Mon programme est écrit en R, et je le mettrai probablement à disposition sur mon gitorious prochainement.

Les données

Pour commencer, nous avons bien entendu besoin d’un dictionnaire adapté, c’est-à-dire contenant uniquement les mots valides selon la règle du “mot le plus long”. Le dictionnaire Lexique contient des informations qui permettent d’identifier la nature grammaticale des mots, ce qui va nous permettre de retirer de la liste les mots interdits par le règlement. Le fichier contient également quelques informations utiles sur chaque mot (sa fréquence d’usage, par exemple), ce qui devrait nous permettre de faire quelques statistiques intéressantes. La dernière mise à jour du dictionnaire Lexique semble dater de 2007. Même si des mots apparaissent et disparaissent chaque année du dictionnaire, un intervalle de 7 ans me semble largement acceptable pour répondre à mes besoins.

Il faut ensuite trouver la probabilité d’appartition de chaque lettre, conditionnellement au choix “Voyelle” ou “Consonne”. Le règlement de l’émission n’est pas très explicite à ce sujet, et après une petite recherche google, il semble que personne ne connaisse réellement les probabilités de sélection des lettres. Qu’à cela ne tienne, nous allons les estimer à l’aide de la loi des grands nombres et d’une base de données ! Jacques Colard tient à jour sur sa page une base de données recensant tous les tirages depuis 1997. Ce travail très impressionnant nous sera également utile pour analyser les choix “voyelle” ou “consonne” des candidats. Le dernier changement dans les règles de jeu date de fin 2011. Afin d’éviter tout problème, mon étude porte sur toutes les émissions diffusées entre janvier 2012 et décembre 2014. Le nombre de tirages dans la base (4671) me semble largement suffisant pour obtenir une bonne précision sur les estimations qui m’intéressent.

Le modèle de probabilité

Nous allons expliciter notre raisonnement à l’aide d’un exemple. Supposons que le jeu se joue avec 4 lettres, et que les candidats ont choisi “Consonne” – “Voyelle” – “Consonne” – “Voyelle” (que je note par le suite le pattern CVCV). Calculons la probabilité de sortie du mot “CAVE” avec ce choix. Commençons par les voyelles : A et E. Pour trouver la probabilité de sélection étant donné la probabilité d’apparition de A et E parmi les voyelles, on peut fonctionner avec un arbre de probabilités (la méthode qui est enseignée au lycée, il me semble) :

arbre_proba_horizontal.svg

En suivant le chemin A, puis E, on obtient :

\Pr(A) \cdot \Pr(E)

Mais nous aurions pu tout aussi bien choisir d’abord “E”, ensuite “A”, ce qui nous aurait donné un autre chemin à suivre dans l’arbre de probabilités. On a finalement :

\Pr(A,E) = 2 \cdot \Pr(A) \cdot \Pr(E)

On peut procéder exactement de même avec les consonnes, ce qui donne finalement :

\Pr(C,A,V,E) = (2 \cdot \Pr(A) \cdot \Pr(E)) \cdot (2 \cdot \Pr(C) \cdot \Pr(V))

Plus généralement, pour tout mot, la probabilité de présence étant donnée le pattern s’écrit :

\Pr(mot | pattern) = NP(voyelles) \cdot \underset{lettre \in Voyelles}{\displaystyle{\prod}} \Pr(lettre | voyelle) \cdot NP(consonne) \cdot \underset{lettre \in Consonnes}{\displaystyle{\prod}} \Pr(lettre | consonne)

avec :

NP(mot) = \text{nombre de permutations des lettres du mot} = \dfrac{(nombre~lettres)!}{\underset{lettres~repetees}{\displaystyle{\prod}} (occurences~lettre)!}

Je ne détaille pas ici le calcul du nombre de permutations des lettres du mot, qui est un problème classique de combinatoire. Je précise simplement qu’il faut penser à prendre en compte le fait que permuter deux lettres identiques ne constitue pas une nouvelle permutation du mot (vous pouvez vous rendre ici pour une explication détaillée).

Il faut ensuite trouver à partir de la probabilité conditionnelle au pattern la probabilité générale sur tous les patterns. On considère que tous les patterns sont des événements disjoints, ce qui permet d’écrire :

\Pr(mot) = \sum_{pattern} \Pr(mot | pattern)

Estimation des probabilités de tirages des lettres

B C D F G H J K L M N P Q R S T V W X Z
N 777 1664 1075 647 960 666 104 109 2094 1248 2901 1196 232 3562 4152 2846 479 32 197 98
% 0.031 0.0665 0.0429 0.0258 0.0383 0.0266 0.0042 0.0044 0.0836 0.0498 0.1159 0.0478 0.0093 0.1423 0.1658 0.1137 0.0191 0.0013 0.0079 0.0039
A E I O U Y
N 4097 8320 3849 2821 2267 317
% 0.1891 0.3839 0.1776 0.1302 0.1046 0.0146

La fréquence d’apparition de chaque voyelle (et de chaque consonne) est utilisée comme d’estimateur de la probabilité de tirage de chaque lettre. Comme attendu, les lettres rares (Y, Z, etc.) sont moins probables que les lettres courantes (A, E, S, etc.).

Estimation des probabilités des patterns

Mon idée était que les choix des candidats étaient différents suivant les lettres déjà tombées (imaginez un tirage où vous avez déjà obtenu Z,H,Y,U – je suppose que votre choix futur ne sera pas nécessairement le même que si vous avez obtenu E,R,S,A), mais les lettres dans la base semblent apparaître toujours dans le même ordre (à nombre de voyelles égal). La base a déjà dû être ordonnée, et je ne peux donc pas différencier les patterns par ordre d’apparition des voyelles et consonnes (bien que cela ne me semble pas être une approximation trop dramatique !). On s’aperçoit par contre que les fréquences ne sont pas du tout symétriques en nombre de voyelles et de consonnes :

CVCVCVCCCC CVCVCVCVCC CVCVCVCVCV VCVCVCVCVV VCVCVCVVVV VCVCVVVVVV
N 72.00 1805.00 2538.00 248.00 7.00 1.00
% 0.02 0.39 0.54 0.05 0.00 0.00

Note : cette façon de procéder suppose que les tirages sont effectués avec remise, ce qui ne semble pas vrai d’après le règlement de l’émission. On réalise donc en réalité une autre approximation. Cette approximation est potentiellement plus grave pour les lettres rares. En effet, supposons que l’on tire sans remise parmi 100 lettres, contenant 30 “E” et 2 “Z”. Si la première lettre tirée est un “E”, la probabilité de tirer un “E” en deuxième lettre est de 29/99 et la probabilité avec remise (celle que nous utilisons), est de 30/100, qui lui est pratiquement égal. Par contre, si l’on a tiré un “Z” en première lettre, la probabilité de tirer un “Z” en seconde lettre est de 1/99, qui est pratiquement deux fois inférieur à 2/100, la probabilité avec remise. Néanmoins, le nombre de mots avec deux “Z” dans le dictionnaire me semble suffisamment limité pour que l’on ne s’inquiète pas trop de cette approximation !

2,5% de chances qu’un tirage contienne un 10 lettres

Les tirages où il existe un 10 lettres semblent assez peu communs dans l’émission. Dans une des dernières émissions que j’ai pu suivre, les présentateurs de l’émission ont signalé la rareté de l’événement en félicitant un candidat qui avait trouvé un mot de 10 lettres. Essayons d’estimer la probabilité de cet événement en utilisant notre table. On peut utiliser la formule suivante : Si les An sont des événements deux à deux incompatibles, on a :

\Pr(\displaystyle{\cup}_n A_n ) = \displaystyle{\sum}_n \Pr(A_n)

En considérant l’union de tous les mots de 10 lettres, sont-ils deux-à-deux incompatibles ? Non, car des anagrammes peuvent coexister au sein d’un même tirage. Il suffit donc de retirer les anagrammes de notre dictionnaire, afin de pouvoir sommer les probabilités de chaque mot de 10 lettres pour obtenir la probabilité qu’il existe un 10 lettres dans le tirage. On obtient :

\Pr(existe~un~10~lettres) \approx 2.49 \%

Une vérification sur la période 2012 – 2014 montre qu’il y a eu 120 tirages avec existence d’au moins un 10 lettres sur 4671 tirages dans notre base. Ce qui donne une fréquence de 2.57% : pas trop mal !

Pour donner un ordre de grandeur de cette probabilité : en partant d’une émission diffusée à un instant t, il y a presque 40% de chances qu’aucun 10 lettres ne sorte avant 40 tirages, c’est-à-dire presque 7 émissions.

Les mots les plus probables (par nombre de lettres)

odalisque.jpg
Enfin, nous y voici ! Nous sommes désormais capables d’attribuer une probabilité à chaque mot valable pour “Le mot le plus long”. Le mot les plus probable est IE (qui est référencé comme nom dans le dictionnaire Lexique, mais n’est peut-être pas accepté par le dictionnaire de référence, étant une abréviation). Le deuxième mot le plus probable est “EU”, participe passé du verbe “avoir”.

Quelques mots ont une probabilité strictement nulle, car ils contiennent 2 voyelles et 7 ou 8 consonnes, patterns inexistants dans notre base 2012-2014, et donc réputés de probabilité nulle : BLANCSBECS, CULSBLANCS, FRANCFORTS, GRANDSDUCS, HALFTRACKS, NIGHTCLUBS, PLATSBORDS, SCHILLINGS, SCHTROUMPF, SCRATCHING,
SCRIPTGIRL, SHRAPNELLS, SPHINCTERS, SPRINKLERS, STRESSANTS, STRETCHING, TCHATCHANT, TRANCHANTS, TRANSCRITS, TRANSFERTS, TRANSPORTS, TREMBLANTS.

Le mot le moins probable à probabilité strictement positive est BOWWINDOWS. Sa probabilité vaut de l’ordre de 10^-15, ce qui signifie que l’on peut espérer le voir apparaître environ une fois toutes les… 26 000 milliards d’émissions !

Comme on peut s’y attendre, la probabilité est fonction décroissante du nombre de lettres, comme le montre ce graphe en échelle logarithmique :

proba_nombre_lettres.png

Je termine en donnant les mots les plus probables par nombre de lettres (de 6 à 10) :

Mots de 6 lettres

mot probabilite (%) anagrammes
AINEES 0.22 ANISEE
OISEAU 0.21
AEREES 0.20
AERIEN 0.19 ANERIE
SOIREE 0.19

Mots de 7 lettres

mot probabilite (%) anagrammes
OSERAIE 0.11
AERIENS 0.10 ANERIES
ATRESIE 0.09
SATINEE 0.08
REALISE 0.07 SALIERE

Mots de 8 lettres

mot probabilite (%) anagrammes
ARTESIEN 0.04 TANIERES,TRAINEES
ARLESIEN 0.03 LANIERES
ALTIERES 0.03 ATELIERS,REALISTE,REALITES
RATISSEE 0.03
EURASIEN 0.03 RAINEUSE

Mots de 9 lettres

mot probabilite (%) anagrammes
INALTERES 0.02 RALENTIES
NOTARIEES 0.02
INSTAUREE 0.01
CARTESIEN 0.01 CERTAINES,SENATRICE
ENLIASSER 0.01

Mots de 10 lettres

mot probabilite (%) anagrammes
ORIENTALES 0.01 ORLEANISTE
INALTEREES 0.01
NEUTRALISE 0.01 RELUISANTE
INSTAUREES 0.01 TRAINEUSES
RATIONNEES 0.01

Certains de ces mots m’étaient inconnus ou ne me seraient jamais venus à l’esprit dans le cadre du jeu (OSERAIE, ATRESIE, ENLIASSER, RAINEUSE, etc.) !

To be continued…

J’arrête ici ce post déjà bien long… mais je ne manque pas d’idée pour exploiter ma nouvelle base de données (que je vous livre en format texte csv) ! Je vous donne donc rendez-vous d’ici quelques jours pour la deuxième partie de cette analyse des mots les plus probables à “Des chiffres et des lettres”.

Télécharger la base des probabilités des mots

Image titre : logo “Des chiffres et des lettres”, © France 3

[Science] Looking for a new holiday or When does day length increase the fastest?

We’re now in the midst of the winter (at least in the northern hemisphere). Our Parisian winter, although not very cold, is quite rainy and cloudy, and it always depresses me a little no to see the sun for several days… Fortunately, at this the time of the year, every next day we can experience more daylight, reminding us we’ve never been closer to the next summer! This made me wonder: when does the length of day  increase the most? My plan was to find out which day of winter it happened, and make it my own holiday, which I would celebrate by opening a bottle of Côte Rotie 2003 I was sparing for a special occasion.

Before I went on with the math, I told a few friends about my question, and they all told me they thought the day when day length increases the most would likely be the vernal equinox. But it didn’t seem that obvious to me! Sure, on March the 21st, the length of daytime will be equal everywhere on our planet, but that does not mean it has to be the same day that day length increases the most, does it? Let’s see who’s right!

Let’s do math!

As I’m not a physicist, I need an equation describing the length of daytime  with respect to some parameters (including, hopefully, the latitude of the city where I live!). Fortunately, Robert Ferreol has already done a great job finding such an equation. For latitudes below the arctic circle (which is approximately 66.566667°), day length writes, in hours:
\(D = \frac{24}{\pi} \arccos \left( – \tan \lambda \cdot \tan \left( \arcsin \left( \sin \delta \cdot \cos \left( \frac{2\pi}{365} (r-172) \right) \right) \right) \right)\)

with:

\(\begin{align*}
\lambda &= \text{latitude} \\
\delta &= \text{Earth axial’s tilt} \\
r &= \text{day of the year, } \in [0,365]
\end{align*}\)

Sure, this equation is a tad complicated, by it doesn’t really matter if your trigonometry lessons seem a bit far, you just need to know how to read simple graphs to understand what follows.

We can plot the day length function, for example for Paris (latitude: 48.8567°).day_length_paris

I’m quite satisfied with this graph: maximum happens around day number 172, which happens to be the summer solstice, and minimum happens around day number 355, which is December 21st, aka the winter solstice. We can also read from our graph that day length on day number 1 (January 1st) is a little more than 8 hours, which matches this more precise time table.
We can also plot our equation for different latitudes. For example, at the equator (latitude 0°), we have:

day_length_equator

Again, this is what we expected, as day length when you’re located on the equator is the same across the year.

Studying the derivative

Great, it seems that we now have a nice equation describing day length with good accuracy. But, what we really want is the variations of day length. In math terms, it means we want to study the derivative of the day length function. Sure, I could compute the derivative by hand, but I am a little lazy (and with a function this ugly, it’s very likely that I would have to start over a few times before getting it right), so I’m going to ask my computer to do it instead. This is called symbolic computation. People in science usually use software like Maple or Mupad, but they’re expensive non-free software, and anyway I have a huge crush on Python, so I will use the excellent Sympy:

import math
from sympy import *
from sympy.plotting import *

r = symbols('r') # This is the variable of my function
lat = math.radians(48.8567) # Paris
delta = math.radians(23.4373408135) # Earth's axial tilt

lengthDay = 24/(math.pi) * acos(-tan(lat) * tan( asin( sin(delta) * cos(2*(math.pi)/365*(r-172)) ) ) )

# Compute the derivative
dLengthDay = diff(lengthDay, r)

print dLengthDay # We print the analytical form of the derivative (...it is as complicated as expected!)
plot(dLengthDay, (r, 0, 365))

I’m not really interested in the closed form of the derivative, but the plot should be very informative. In fact, for Paris, it gives:

deriv_paris

So, what is happening here? If you remember your calculus lessons correctly, when the derivative is greater than 0, the function increases, and when it’s lesser than 0, the function decreases. So it’s not a big surprise to see that the derivative is above the x’s axis before r = 172 (the summer solstice), and below the x’s axis between r = 172 and r = 355 (ie between the two solstices). In general, the value of the derivative gives the rate at which the function increases, which means the day we’re looking for has to be the one where the maximum of the derivative happens. In our graph, it seems to be around r = 80, but it’s not very easy to identify a clear maximum graphically. So, what is the best method to compute the maximum of a smooth function like this one? Yep, we derive it (yet again you might add. To be more specific, the day when day length increases the most is called an inflection point, which is when the second derivative equals 0).

# Second derivative
d2LengthDay = diff(dLengthDay, r)
print d2LengthDay
plot(d2LengthDay, (r, 0, 365))

The plot of the second derivative for Paris is:

deriv2_paris

To identify the day when the derivative is zero, we can do it numerically with a Newton method (seeded with r=80, which visually seems to be a good approximation of the zero):

print nsolve(d2LengthDay, 80)

The nsolve returns:

80.7499999999996

Now it’s clear that our maximum is unique, and happens around day number 80(*), which is… March 21st, the vernal equinox! Darn, it seems my friends were right!

(*) If you wonder why I turn 80.749999 into “day number 80” and not “day number 81”, check out the last paragraph of this post.

 

Wait a minute, what about other latitudes?

I was not entirely satisfied with this answer, so I started making tests for other latitudes. The shape of the derivative was very similar to the one I got with Paris, until I started tests with northern latitudes. For example, at 63.5°N, we get:

deriv_fairbanks

What does this graph mean? It appears that: for r = 80 (spring equinox), we now have a local minimum of the derivative instead of the maximum we had in Paris. The are now two local maxima, one happening between winter and spring, and the other one happening between spring and summer.
We can check this by plotting the second derivative: if we’re right, it will show three zeros between r=1 and r=172:

deriv2_fairbanks

So it seems we were right! Zeros seem to be around 40, 80 and 120. But we can compute the dates of these zeros more precisely. For example, for Fairbanks, Alaska (latitude = 64.843611°), we write:

print nsolve(d2LengthDay, 50)
print nsolve(d2LengthDay, 80)
print nsolve(d2LengthDay, 120)

and get:

39.8279585264938
80.7499999999996
121.672041473505

The first one (r approximately equals 40), means that February, 9th is (one of) the inflection points we were looking for. To ensure that my equation was accurate enough, I compared my result with an astronomical time table. If you look at the “daylength” column, you will see that March 21st is indeed a local minimum, and that there are two local maxima for day length.

 

Computing the limit latitude

After a few tests, I started realizing that for very latitude above approximately 62°N, there were always two maxima, whereas for every latitude below, the only maximum that existed was around spring. This probably means that there exists a limit latitude, above which there are always two maxima, and below which always only one.
To compute (numerically) the limit latitude, we can write a quick and simple binary search:

def searchLimitLatitude(lat, count=0, lastBelow = 0, lastAbove = math.radians(66.56278), latBackup = lat):

	print "Test for latitude = ", math.degrees(lat)
	print "Iteration number: ", count

	lengthDay = 24/(math.pi) * acos(-tan(lat) * tan( asin( sin(delta) * cos(2*(math.pi)/365*(r-172)) ) ) )
	dLengthDay = diff(lengthDay, r)
	d2LengthDay = diff(dLengthDay, r)

	# How many distinct maxima are there?
	max1 = nsolve(d2LengthDay, 50)
	max2 = nsolve(d2LengthDay, 80)

	if(count >=1 and math.fabs(lat - latBackup) <= 1e-7 ): # we ask for 7-digit precision
		print "Finished, limit latitude is: ", math.degrees(lat)
		return lat

	if( math.fabs(max1 - max2) <= 1e-16 ):
		searchLimitLatitude(0.5 * (lastAbove + lat), count+1, lat, lastAbove, lat)
	else:
		searchLimitLatitude(0.5 * (lastBelow + lat), count+1, lastBelow, lat, lat)

print "Search limit latitude...."
lat = math.radians(50) # Initialization
searchLimitLatitude(lat)

Running this program gives:

Iteration number:  22
Finished, limit latitude is:  61.2497155042

As expected, the binary search converged very quickly (just 22 iterations to give a latitude precise within 7 digits!). Finally, we have a value for our limit latitude:

\lambda_{lim} \approx 61.2497 = 61 $^{\circ}$ 14 ‘ 59 ”

which, according to wikipedia, is just above Anchorage, Alaska and slightly below Tampere, Finland

Conclusion

My friends were partly right: I won’t be able to act snob and celebrate my own holiday. In fact, where I live (Paris, France), inflection point happens on the spring equinox, which many cultures have already been celebrating for centuries.
But, if you live in Fairbanks, Alaska, I suggest you call your boss to let him know you won’t be coming to work on February, 9th (which happens to be the day this is posted!);)

A few more precisions

You might wonder why when we found an inflection point on date 80.74999 for Paris, I said the corresponding day was #80 (March 21st), and when we found 39.8279 for Fairbanks, I said the corresponding day was #40. In fact, we reach here the limits of our simplified equation. One of the assumptions we made was to consider that the Sun was located at infinite distance from Earth (ie that the Sun was just a tiny point in the sky). Of course, this is not true, and it can lead to small errors if we stick to our equation. For example, on the vernal equinox at the North Pole, day length is really about 6 minutes longer than our equation predicts. In order to decide which day was really the day of the maximum increase in day length, I had to “cheat” a little (because I needed a more precise equation) by checking an astronomical table for these locations.

Something also worth noting is that there is also a slight modification of Earth’s tilt axis every year, which means the second derivative (and thus the inflection point) also changes every year. I don’t really have time for more than a few tests, but location of the inflection point seems to be very sensitive the parameters of the equation. This means that inflection point in Fairbanks could happen on very different dates each year.

EDIT : As Boris pointed out in the comments, I could have dealt with my problem more simply, by writing day length as a discrete series of 365 days, instead of deriving a continuuous function. With the continuuous method, day number “39.9” doesn’t make much sense, and I couldn’t decide whether it was February, 8th (#39) or February, 9th (#40). Using finite sequences, I needn’t have used the astronomical table for Fairbanks, for the day I was looking for would just have been the max of the delta series.