Tuesday morning, Antoine Ly will defend his PhD thesis.
And earlier in the morning (breakfast time) I will give a talk with Molly James at the Chaire PARI Seminar. Slides are available online.
Tuesday morning, Antoine Ly will defend his PhD thesis.
And earlier in the morning (breakfast time) I will give a talk with Molly James at the Chaire PARI Seminar. Slides are available online.
Entendre « il y a 10% de chances de pluie aujourd’hui » ou « le test médical a une valeur prédictive positive de 75% » montre que les probabilités sont aujourd’hui partout[i]. Une probabilité est une grandeur difficile à appréhender, mais incontournable quand on cherche à théoriser et mesurer le hasard. Et si la théorie mathématique est finalement arrivée très tard, comme le rappelle Hacking (2006), cela n’a pas empêché l’assurance de se développer assez tôt, et d’avoir les premières tables (actuarielles) de mortalité avant même que la « probabilité de décès » ou « l’espérance de vie » n’ait de fondement mathématique. Et de la même manière, de nombreuses techniques ont été inventées pour « générer du hasard », avant l’explosion des méthodes dites de Monte Carlo, en parallèle avec le développement de l’informatique (et du fait qu’une machine pouvait générer du hasard). Continue reading (Petite) histoire du hasard et de la simulation
Il y a quelques mois, Baptiste Coulmont m’avait contacté avec une question passionnante (comme chaque fois qu’il me contacte). Baptiste voulait travailler sur la proportion de personnes qui ont un homonymes dans une population de taille donnée, ou sur la probabilité de ne pas avoir d’homonymes dans un bureau de vote, par exemple. Or ce dernier problème n’est pas sans rappeler le “paradoxe des anniversaires“.
Dans un groupe de 23 personnes, il y a plus d’une chance sur deux pour que deux personnes de ce groupe aient leur anniversaire le même jour. Ou dit autrement, il est plus surprenant de ne pas avoir deux personnes qui ont leur anniversaire le même jour que d’avoir deux personnes qui ont leur anniversaire le même jour.
On parle de “paradoxe” parce que le nombre de jours possibles est relativement grand (k=365), mais que la taille du groupe est elle relativement petite (n=23). Ici c’est un peu la même chose, sauf qu’on est sur des nombres beaucoup plus grands: on a des milliers de noms possibles (en fait, on s’intéresse au couple prénom+nom, ce qui augmente très largement le nombre de possibilités, que l’on notera k ici), on considère des groupes relativement grands (la taille du groupe sera elle notée n), et surtout, l’hypothèse (implicite) d’uniformité dans la distribution des dates ne marche pas sur les noms de famille, et les prénoms. On peut malgré tout tenter de formaliser le problème.
Considérons une collection de k éléments, notée \mathcal{X}=\lbrace x_1,\cdots,x_k\rbrace, et faisons n tirages avec remise. Les tirages obtenus sont notés ici X_1,\cdots,X_n et correspondent à des variables aléatoires iid, suivant des lois multinomiales sur \mathcal{X} avec comme vecteur de probabilités \mathbf{p}=(p_1,\cdots,p_k). Dans le cas des anniversaires, \mathcal{X} correspond aux dates, avec k=365, et \mathbf{p} est uniforme sur \mathcal{X} – ou encore \mathbf{p}=\mathbf{1}/k . Le problème classique est de calculer la probabilité que dans un groupe de taille donnée tout le monde ait des dates d’anniversaire différentes
\mathbb{P}[\forall j\neq i, X_i\neq X_j]mais on peut aussi calculer (comme le voulait Baptiste) la proportion de personne ayant un alter-ego (une personne dans le groupe qui a son anniversaire le même jour)
\frac{1}{n}\sum_{i=1}^n \mathbf{1}(\exists j\neq i,X_i=X_j)
On a donc tenté d’approcher le second, quand on a un groupe de taille n. Notons Z_i le nombre de personnes (parmi les n-1 restantes) qui partagent la même date d’anniversaire que i,
Z_i=\sum_{j\in\{1,\cdots,n\}\backslash\{i\}} \mathbf{1}(X_j \neq X_i)
Aussi, la personne i a un alter-ego si Z_i\geq 1. Et donc, la proportion de personnes avec un alter-ego dans le groupe s’écrit
P_{n,k}=\frac{1}{n}\sum_{i=1}^n \mathbf{1}(Z_i\geq 1)
Si \mathbf{p} correspond à la loi uniforme sur \mathcal{X}, alors
1-P_{n,k}=\binom{k}{n}\frac{1}{k^n}\sim \exp\left[\binom{n+1}{2}\frac{1}{k}\right]\sim e^{-n^2/2k}
comme le montre Chatterjee et al. (2004). Je pourrais aussi renvoyer aux articles de Inoue & Aki (2008) et plus récemment de Cortina Borja (2013) pour des revues de littérature sur les aspects computationnels.
Mais la plupart des résultats sont obtenu dans le cas où \mathbf{p} correspond à une distribution uniforme. Les propriétés plus générales sont beaucoup plus rares(on pourra lire Munford (1977), Das Gupta (2004), Inoue & Aki (2008) ou Nunnikhoven (1992) pour des tentatives. Utiliser des méthodes de Monte Carlo est probablement le plus naturel (si on n’est pas trop pressé) en se donnant un ensemble \mathcal{X} (qui sont les dates dans le problème des anniversaires et les noms dans le cas des homonymes) de taille k, et une distribution \mathbf{p} sur \mathcal{X}. On a ainsi pu approcher P_{n,k}, la proportion de personnes ayant un alter-ego dans une population de taille n: sur les figures ci-dessous, on a cette proportion en fonction de n (en abscisse), avec des courbes différentes pour différents k. Histoire d’être un peu réaliste, on va chercher le taux de personnes ayant un homonyme dans une popoluation de plusieurs centaine de milliers de personnes, et quand aussi des milliers, des dizaines de milliers, de possibilités.
Dans le cas d’une distribution uniforme, la proportion évolue comme sur la figure ci-dessous
On a tenté de quitter l’indépendance en utilisant une loi paramétrique, plus précisément une loi de type Pareto, ou Zipf, pour \mathbf{p}. En fait, c’est le type de distribution qu’avait observé Li (2012) sur la distribution des prénoms, avec p_i \propto i^{-\alpha} (ou plus récemment, Höhle (2017)). On a alors l’évolution suivante
On est ici en échelle logarithmique pour la taille n parce que c’est naturel. Mais pour l’instant, on était resté sur P_{n,k}. Baptiste insistait pour passer aussi à une échelle logarithmique sur l’axe des ordonnées, mais ca n’etait pas concluant. On a alors tenté deux autres représentations, plus classiques quand on travaille sur des probabilités (ou des proportions) : représenter \text{logit}(P_{n,k}) en fonction du logarithme de n (où \text{logit}(u)=\log[u/(1-u)]) et une version probit, à savoir \Phi^{-1}(P_{n,k}) (où ici \Phi est la fonctionc cumulative de la loi normale centrée et réduite).
En grande dimension, on a l’impression d’avoir une relation linéaire entre la transformée probit de la proportion de personnes ayant un homonyme et le logarithme de la taille de la population considérée
\Phi^{-1}(P_{n,k}) \sim a_{k}+b_{k} \log[n]C’est joli, mais on se demandait comment utiliser ca sur nos données, pour travailler vraiment sur les homonymes. Le soucis, c’est qu’on n’a pas de fichier national, avec les noms et les prénoms. Alors Baptiste a sorti de son sac les listes électorales parisiennes. Le soucis est que la base est (relativement) petite. Le nombre de couple est très limité. Au lieu de faire des tirages avec remise dans les couples (“prénom”, “nom”) – encore une fois, la base est trop petite – on a voulu tenter de tirer indépendemment le “prénom” et le “nom”. On y gagne en ayant à disposition une base plus grande, mais on introduit un biais colossal, les “prénom” et “nom” étant assez corrélés. Dans notre cas, on a k_1 “prénom” et k_2 “prénom”. Il y a quelques années Chalabi & Flowers (2014)} (repris plus en détails sur ce blog) était revenu sur la non indépendance des noms et des prénoms, en s’interrogeant sur le nom le plus populaire aux Etats-Unis. Selon les auteurs “Michael Smith” n’était probablement pas le nom le plus populaire, même si “Michael” et “Smith” étaient respectivement les “prénom” et “nom” les plus populaires. En France “Jean” et “Marie” sont les prénoms les plus populaires pour les garcons et les filles, respectivement, et “Martin” est le nom le plus populaire (au moins sur les listes électorales parisiennes), mais en fait, le nom le plus populaire est “Thi Nguyen”, avec 350 personnes à Paris, alors que seulement 84 “Jean Martin” sont observés, le troisième plus haut score. Ah oui, je dois préciser qu’on se limite ici au premier prénom, en cas de prénoms composés. Tout comme Chalabi & Flowers (2014), on a fait un test du chi-deux pour tester l’indépendance, sur les 20 “noms” et “prénoms” les plus populaires. Sur la figure ci-dessous, on a les résidus normalisés du test. Par exemple “Michel Michel”est en rouge, car c’est beaucoup moins donné que ce qu’on aurait si les deux composantes étaient indépendantes
Bref, on n’a pas trop le droit de tirer le “prénom” et le “nom” indépendament. Mais on va quand même le faire, pour voir. Ici \mathcal{X} est l’ensemble \mathcal{X}_1\times\mathcal{X}_2 avec k_1 choix pour le “prénom” et k_2 pour le “nom”.
Soit \mathbf{p}=(p_{i_1,i_2}) la probabilité empirique sur \mathcal{X}. On imagine que \mathbf{p} contient énormément de 0, car beaucoup de couples ne sont pas observés. Pour les simulations, on va noter \mathbf{p}^{\perp} la pobabilité jointe sous hypothèse d’indépendance
p^{\perp}_{i_1,i_2}=\left(\frac{1}{k_2}\sum_{j=1}^{k_2} p_{i_1,j}\right)\cdot\left(\frac{1}{k_1}\sum_{j=1}^{k_1} p_{j,i_2}\right)=p_{i_1,\cdot}p_{\cdot,i_2}
Lors de nos simulations, on va considérer deux cas : on dire des \mathbf{x}=(x_1,x_2) soit suivant la distribution \mathbf{p}, soit suivant \mathbf{p}^{\perp}. Bon, j’avais formalisé pour clarifier un peu dans ma tête et aussi pour impressionner Baptiste, mais techniquement, dans un cas on tire une paire observée – genre “Baptiste Coulmont” – et dans le second cas, on tire d’un côté un prénom – genre “Baptiste” – et de l’autre un nom de famille – genre “Charpentier” (même s’il n’existe aucun “Baptiste Charpentier”, on dit juste qu’il serait possible qu’il y en ait un).
Bref, on tirant des échantillons (plusieurs milliers de taille donnée, puis en moyennant), on obtient la figure ci-dessous (avec une échelle logarithme pour l’axe des ordonnées)
La courbe noir sur-estime sûrement à cause du faible nombre de choix possibles (on va tirer jusqu’à 50 millions de noms). En revanche, l’indépendance n’étant pas réaliste, la courbe rouge sous-estime probablement la proportion. Mais le plus intéressant est le graphique ci-dessous, obtenu en représentant la version probit de la probabilité
La courbe rouge est incroyablement linéaire ! Bien plus linéaire que ce qu’on espérait avoir suite à nos simulations dans le cas de la loi de Pareto.
En fait, on note que si on se limite au cas où n n’est pas trop grand, la courbe noire semble aussi linéaire. Ensuite on perd la linéarité parce qu’on tire avec remise dans une base trop petite, mais au début (lorsque n n’est pas trop grand), on pourrait accepter l’hypothèse d’une relation linéaire entre la transformée probit, et le logarithme de n. Bref, on estime une relation linéaire, et on extrapole pour les valeurs de n plus grandes
Si on revient sur la représentation plus classique, de la proportion de personnes ayant un homonyme en fonction de n, on obtient
Autrement dit, en France, on peut estimer que 80% des francais ont un homonyme.
Deux points pour poursuivre, sur la robustesse de l’approche. Tout d’abord, si on change un peu la taille de l’intervalle des n pour lesquels la courbe noire semble linéaire, la pente change peu
et on a un ordre de grandeur assez proche. Bref, ce 80% semble assez robuste,
Ensuite, comme cette relation linéaire me surprenait, on a aussi tenté sur des données américaines… et on a exactement le même genre de résultats…
Pour aller plus loin, on a tout mis sur github (avec des graphs en ggplot2 pour faire plus joli…) et une version (en anglais) est en ligne sur arxiv.
Pour les rappels, j’utiliserais a priori des transparents car il y a beaucoup de matière à couvrir.
Je repasserai au tableau une fois les rappels finis.
J’ai découvert avec un peu de retard le joli billet Les anniversaires de vos amis sur Facebook, qui tentait de répondre à la question
Si je possède amis, quelle est la probabilité qu’il y ait au moins un jour dans l’année où je n’ai pas d’anniversaire à souhaiter ?
Ce problème, on peut aussi l’analyser en posant plutôt la question suivante,
Il y a en tout 365 jours dans une année. Combien d’amis doit-on espérer avoir pour qu’il n’y ait pas un jour sans fêter un anniversaire d’un de ses amis ?
Le problème est différent. Mais la résolution est assez jolie, si on regarde attentivement.
Si je note le nombre d’amis qu’il me faut pour couvrir tous les anniversaires, et que je commence à regarder, ami après ami le nombre de jours couverts. Supposons que j’ai jours de l’année couverts, et notons le nombre de nouveaux amis qu’il va me falloir avant de couvrir jours de l’année. On peut voir que
Le point intéressant avec cette notation, c’est qu’on connaît la loi de . En fait, on reconnaît une loi géométrique
Aussi, on peut assez facilement avoir le nombre moyen d’amis qu’il va falloir pour couvrir tous les jours de l’année,
Je laisse les amateurs de calculs essayer de simplifier, mais ça me semble compliqué. Par contre, numériquement on obtient facilement
> sum(365/(1:365)) [1] 2364.646
En fait, la somme, on peut quand même la simplifier un peu, car on reconnaît
Et on sait que
où est la constante d’Euler. Et effectivement, avec cette approximation, on obtient
> (log(365)+ 0.57721)*365 [1] 2364.144
Ce nombre d’amis est à rapprocher du calcul fait sur https://blogdemaths.wordpress.com/, qui obtenait
on a une chance sur deux pour qu’il y ait au moins un jour sans anniversaire à souhaiter, il faut avoir environ 2285 amis
Calculer un quantile est un peu plus compliqué pour moi (je ne peux calculer que des moyennes), mais on voit que les ordres de grandeurs sont assez proches.
Vendredi, suite du cours ACT2121, de préparation pour l’examen P de la SOA (probability). Un nouveaux jeu d’exercices, sur les thèmes 10, 11, 12 et 16 (tels que classifiés dans le livre de Jacques Labelle, qui sert de référence pour ce cours)
L’examen de lundi est en ligne ici, avec des éléments de correction là (avec comme toujours des statistiques sur les réponses). Les notes seront publiées bientôt. Toute personne qui trouve des erreurs peut me contacter avant que je ne valide les notes. L’examen correspondait à l’examen de pratique 15 du livre de Jacques Labelle (je n’ai rien inventé). Pour une des questions, la réponse n’était pas parmi les réponses proposées. J’ai donc finalement noté l’épreuve sur 29 (et appliqué un coefficient multiplicatif pour ramener à une note sur 30). Comme toujours, ceux qui ont prédit correctement leur nombre de bonnes réponses ont eu un point bonus.
La saison 1 de Mad Men se termine avec l’élection de 1960, Nixon contre Kennedy, qui s’est jouée dans un mouchoir de poche.
L’autre jour, Benoit Rittaud évoquait sur http://images.math.cnrs.fr/ la possibilité que des candidats soient ex-æquo dans une élection, et les aspects légaux (car visiblement la constitution française n’a pas vraiment prévu ce cas). Le plus intéressant aurait été – à mon avis (mais mon avis est fortement biaisé) – de calculer la probabilité qu’un tel évènement survienne. Par exemple, qu’à l’issu du premier tour, le deuxième et le troisième candidat soient ex-æquo, de telle sorte qu’on ne puisse déterminer qui ira au second tour. Justement, un des commentaire affirme “la probabilité qu’un tel événement arrive est environ l’inverse du nombre typique de voix séparant le deuxième homme (ou femme) du troisième”. N’ayant réussi à démontrer ce résultat, je me suis dit qu’on pourrait se lancer dans des simulations pour quantifier cette probabilité.
Commençons par supposer qu’il y a quatre candidats (et la possibilité de s’abstenir). On suppose que le taux d’abstention est de 20%, et que les intentions de vote pour chacun des candidats soit de 20%.
EXE=function(N=100,ns=1000000){ ExE=rep(NA,ns) for(s in 1:ns){ M=sample(c("A","B","C","D","E"),size=N,prob=rep(.2,5), replace=TRUE) tb=table(M) Ms=rev(sort(tb[-(names(tb)=="A")])) ExE[s]=(Ms[2]==Ms[3]) } return(mean(ExE))}
Le soucis est que faire plusieurs centaines de millions de tirages de plus de 40 millions de lois multinomiales pourrait prendre du temps (genre vraiment beaucoup). Une stratégie intuitive peut être de faire des simulations pour une population plus faible, puis d’extrapoler.
E=function(N) EXE(N) taille=c(20,50,100,200,500,1000,2000) proba=Vectorize(E)(taille) plot(taille,proba,type="b",log="xy") base=data.frame(y=proba,x=taille) reg=lm(log(y)~I(log(x)),data=base) s2=summary(reg)$sigma m=exp(predict(reg,newdata= data.frame(x=41194689)))*exp(s2^2/2)
L’extrapolation linéaire n’étant peut-être pas parfaite, on peut tenter un ajustement quadratique (comme cela a été fait sur la figure ci-dessus)
reg2=lm(log(y)~I(log(x))+I(log(x)^2),data=base) s2=summary(reg2)$sigma m2=exp(predict(reg2,newdata=data.frame(x=41194689)))*exp(s2^2/2)
La différence de prédiction avec les deux modèles est significative comme on le voit sur le graphique. Toutefois, avoir des ex-æquo avec une très très grande population est tellement rare qu’il faut réellement beaucoup de simulations pour avoir un ordre de grandeur valide. L’extrapolation semble être une alternative valide ici (mais je suis ouvert à toute critique, ou discussion).
Maintenant que l’on a une méthodologie (ou presque), utilisons des ordres de grandeurs réalistes. Par exemple en se basant sur l’élection présidentielle de 2002. 41.19 millions d’électeurs, 28.4% d’abstentions, et pour les votes exprimés, un candidat en tête avec 19.9% des votes, un second à 16.8% et le troisième à 16.2% (et les autres plus loin). Les intentions de votes pour les principaux candidats figurent dans la loi multinomiale ci-dessous
EXE=function(N=100,ns=1000){ ExE=rep(NA,ns) for(s in 1:ns){ M=sample(c("abs","JC","JMLP","LJ","FB","AL","JPC" ,"NM","autres"),size=N,prob=c(.284,.1988*.716, .1686*.716,.1618*.716,.0684*.716,.0572*.716, .0533*.716,.0525*.716,.1714104),replace=TRUE) tb=table(M); Ms=rev(sort(tb[-(names(tb)%in%c("abs","autres"))])); ExE[s]=(Ms[2]==Ms[3]) } return(mean(ExE))}
Si on en croit l’affirmation du commentaire, la probabilité d’avoir eu des ex-æquo (entre le second et le troisième candidat) serait de l’ordre de
> 1/(41194689*(.1686*.716-.1618*.716)) [1] 4.985823e-06
Je dois avouer que les calculs n’ont rien donné (on verra tout à l’heure ce qui se passe si on change la définition d’ex-æquo). Numériquement, les calculs tournent encore (après 4 jours) sur le serveur… Par contre, si la population ne contenait que 10000 électeurs, la probabilité d’avoir des ex-æquo (entre le 2ème et le 3ème candidat) serait de l’ordre de 3 sur 10000 (avec un intervalle de confiance que je n’ai toujours pas calculé)
> EXE(10000,100000000) [1] 3.245562e-05
A titre d’illustration, on peut regarder ce qui s’est passé, commune par commune en France (les données sont en ligne sur le site du Ministère de l’Intérieur). L’exercice est biaisé parce qu’il y a de l’hétérogénéité spatiale… mais ça permet toujours d’illustrer… Les nombres de voix obtenues, pour le second et le troisième candidat, par commune est le suivant
> election=read.table( + "/Users/UQAM/Documents/Pres2002Tour1.csv",sep=";") + election=election[-1,] > VOIX=election[,seq(14,74,by=4)] > DEUXIEME=TROISIEME=TOTAL=rep(NA,nrow(VOIX)) > for(i in 1:nrow(VOIX)){ + u=rev(sort(as.numeric(as.character( + election[i,seq(14,74,by=4)])))) + TOTAL[i]=sum(u) + DEUXIEME[i]=u[2] + TROISIEME[i]=u[3] + }
Sur 37000 communes,
> nrow(election) [1] 37513
ou plutôt sur les 24000 communes qui ont eu plus de 5000 votes exprimés
> sum((TOTAL>5000)) [1] 24230
il y a eu 56 cas avec de deuxième et troisième candidats ex-aequo,
> base=data.frame(TOTAL,DEUXIEME,TROISIEME,COMMMUNE=election[,3]) > selection=(TOTAL>5000) > sousbase=base[selection,] > sum(sousbase$DEUXIEME==sousbase$TROISIEME) [1] 56
On a même le détail,
> sousbase[indice,] TOTAL DEUXIEME TROISIEME COMMMUNE 1134 5819 699 699 Crouy 1507 5386 717 717 Louroux-Bourbonnais 2603 6235 1151 1151 Saint-Peray 2916 6346 935 935 Coucy 2970 7529 1413 1413 Haraucourt 5089 5567 717 717 Vignats
(etc). On peut alors regarder, en fonction de la taille des communes (ici le nombre de votes exprimés) la probabilité d’avoir un ex-æquo, dans une commune donnée. On va subdivisionner en 50 classes (pour commencer), et compter à chaque fois le nombre d’ex-æquo,
> n=50 > niveau=seq(1/(2*n),1-1/(2*n),by=1/n) > base$CLASSE=cut(TOTAL,quantile(TOTAL,seq(0,1,by=1/n)), + labels=niveau) > base$EXAEQUO=(base$DEUXIEME==base$TROISIEME)*1 > m=by(base$EXAEQUO, base$CLASSE, sum)
Par exemple pour les 2% des plus grosses villes, on a eu 5 cas de communes avec des ex-aequos
------------------------------------------------------------ base$CLASSE: 0.99 [1] 5
Si on regarde plus en détails, ce sont les communes avec plus de 7500 votes exprimés,
> quantile(TOTAL,.98) 98% 7427.76 > selection=(TOTAL>quantile(TOTAL,.98)) > sousbase=base[selection,] > sum(sousbase$DEUXIEME==sousbase$TROISIEME) [1] 5 > indice=which(sousbase$DEUXIEME==sousbase$TROISIEME) > sousbase[indice,] TOTAL DEUXIEME TROISIEME COMMMUNE CLASSE EXAEQUO 2970 7529 1413 1413 Haraucourt 0.99 1 14359 7496 1392 1392 Renac 0.99 1 25028 8104 1381 1381 Heloup 0.99 1 30014 8087 1340 1340 Barnay 0.99 1 36770 7435 764 764 Pont-sur-Yonne 0.99 1
On peut représenter ces nombres de cas d’ex-æquo et faire une régression de Poisson (car on fait du comptage ici),
> taille=length(TOTAL)/n > plot(quantile(TOTAL,niveau),m/taille) > sb=data.frame(niveau,exaequo=as.vector(m)) > library(splines) > reg=glm(exaequo~bs(niveau),family=poisson,data=sb) > u=seq(0,1,by=.001) > lines(quantile(TOTAL,u),predict(reg, + newdata=data.frame(niveau=u), + type="response")/taille,lwd=2,col="red")
pour des facilités de lecture, on a en ordonnées des probabilités et en abscisse les tailles moyennes des villes
On peut faire la même chose avec 200 classes,
Bref, pour une commune de taille moyenne (disons moyenne supérieure, avec 5000 votes exprimés), on a 2 chances sur 1000 d’avoir un second et un troisième candidat ex-æquo. Ce qui laisserait entendre que pour 40 millions, la probabilité devrait être faible, voire très faible…
Bon, au delà du problème conceptuel des ex-æquo (ou légal, car je doute qu’il soit possible d’avoir une égalité parfaite, le Ministère de l’Intérieur ferait un recompte et trouverait forcément un résultat différent), il semble que la probabilité semble beaucoup plus faible que celle obtenue par la règle ad-hoc….
Mais peut-être serait-il réaliste de dire que l’élection n’a pas su dégager une majorité (c’est le but d’une élection il me semble) si le nombre de voix séparant deux candidats est trop faible. Par exemple, si moins de 100 voix séparent les candidats, on refait un vote ! Pour revenir sur mon exemple précédant, avec 10000 votants
> EXE100(10000,1000000) [1] 0.012641
soit une chance sur 100 (je me suis autorisé à faire moins de simulations car l’évènement me semble moins improbable… donc il y a besoin de moins de simulations pour l’observer à l’occasion). Si on va plus loin, et que l’on fait des simulations que l’on projette, comme auparavant, pour estimer la probabilité d’avoir des ex-æquo non pas pour une ville, mais au niveau de la France entière, on est plutôt aux alentours de
> (m=exp(predict(reg,newdata= + data.frame(x=41194689)))*exp(s2^2/2)) 1 1.505422e-12
pour un modèle linéaire, et
> (m=exp(predict(reg2,newdata= + data.frame(x=41194689)))*exp(s2^2/2)) 1 4.074588e-98
pour un modèle quadratique (avec la confiance – relativement faible – que l’on peut accorder à une telle projection). Ce qui donne une différence énorme… Si on visualise tout ça,
lines(c(1,100000000),exp(predict(reg,newdata= data.frame(x=c(1,100000000))))*exp (s2^2/2),lty=2,col="light blue")
que l’on représente en bleu, en vert un ajustement quadratique
lines(10^(seq(1,8,by=.1)),exp(predict(reg,newdata= data.frame(x=10^(seq(1,8,by=.1)))))*exp (s2^2/2),lty=2,col="dark green")
alors que le calcul suggéré dans le commentaire donnerait la courbe rouge,
lines(c(1,100000000),1/(c(1,100000000)* (.1686*.716-.1618*.716)),col="red",lty=2)
Peut-on parler d’évènement improbable ? A titre de comparaison, quand on parle de quelque chose d’improbable, ou qui arrive avec une probabilité infinitésimale, on pense au loto. A l’ancien loto (français), il fallait choisir 6 numéros parmi 49 (sans remise). La probabilité de gagner le gros lot était de
Pour l’Euromillions, il faut choisir 5 numéros parmi 50 numéros auxquels il faut ajouter les combinaisons possibles avec les Étoiles, soit 2 numéros à choisir parmi 9. La probabilité de gagner le gros lot est alors de
Bref, tomber sur des ex-æquo lors d’une élection nationale, en France, c’est a priori plus rare que gagner le gros lot en jouant une fois au loto… On peut donc dire qu’il ne devrait pas y avoir d’ex-æquo à l’issu du premier tour !
Je suis tombé un peu par hasard sur un vieux papier de Peter Enis, Enis (1973), qui revient sur une relation que l’on prend toujours pour acquis (moi le premier) à savoir . Mais il faut des hypothèses complémentaires pour que cette relation soit vraie. Le fait que existe par exemple ne suffit pas ! (je laisse les amateurs de résultats contre-intuitifs le temps de réfléchir à cette affirmation)
Supposons par exemple que , conditionnellement à , suive une loi normale centrée, d’écart-type ,
où est une variable (strictement) positive. Presque sûrement quelle que soit la loi de . Et donc . Mais si la loi de présente des queues “trop épaisses” (sans pour autant être extrême pour autant d’ailleurs), il est fort possible que ne soit plus intégrable ! Par exemple (c’est ce que propose Enis (1973)), on peut supposer que suit une loi Gamma de paramètres 1/2 et 2, i.e. de densité
Dans ce cas, suit une loi de Cauchy, qui n’est pas d’espérance finie. C’est pénible, hein ?
Bon, d’un autre côté, c’est très lié à la définition de l’espérance conditionnelle retenue. En l’occurrence, j’utilisais l’espérance de la loi conditionnelle. Mais je ferais peut être un billet un jour sur les diverses constructions de l’espérance conditionnelle.
Petite précision liminaire: je suis nul en sondage. Je n’ai jamais rien compris aux sondages. C’est la matière où j’ai eu ma pire note quand j’étais étudiant. Mais si on voit un sondage comme l’estimation d’une probabilité de succès, dans des tirages de lois binomiales (ou multinomiales), alors je me sens plus à l’aise pour commenter des choses lues ici ou là. Car en ces temps d’élections (qui approchent), on lit beaucoup de choses sur les sondages ou sur les intervalles de confiance des sondages, voire sur des comparaison de sondages… mais rares sont les justifications. Donc on va creuser un peu (histoire de voir ce qui peut être dit, dans un cadre largement simplifié).
comme le disait @adelaigue (sur Twitter). Faire un sondage politique (si on résume simplement) c’est estimer la probabilité de voter pour A, sachant que l’on a alors une probabilité pour qu’une personne prise au hasard ne vote par pour A. Sur personnes interrogées, le nombre de personnes qui déclarera voter pour A est , qui suit une loi binomiale de paramètres et . Or par approximation (ou convergence) de la loi binomiale par une loi Gaussienne, on peut dire que
Or l’estimateur naturel de est , i.e.
Aussi, on en déduit l’intervalle de confiance classique, tel que appartient à cet intervalle avec une probabilité ,
où désigne la fonction quantile de la loi normale centrée réduite; ou encore, en remplaçant les quantités inconnues par des estimateurs,
Les bornes de l’intervalle de confiance de niveau sont alors
C’est ce qu’on peut visualiser sur la figure ci-dessous,
Comme le notait @adelaigue, si on souhaite avoir
(car on parle d’un point de base 100), il faudrait que
Classiquement, est aux alentours de 1000. Et ce qu’on cherche c’est une probabilité d’appartenir à un intervalle, i.e. ,
Graphiquement, en fonction de on a alors
p=seq(0,1,by=.0025) plot(p,2*(1-pnorm(.01/sqrt(p*(1-p)/1000))))
(ca sera notre premier graphique). Effectivement, entre 35% et 65%, il y a environ 50% de chances d’être autour de la fréquence observée, à un point près. Ou 50% de chances d’avoir plus d’un point de différence. Autrement dit, si on a observé avec 1000 personnes une fréquence de 40%, il y a une chance sur 4 que la probabilité soit entre 39% et 40%, une chance sur 4 qu’elle soit entre 40% et 41%, et surtout une chance sur 4 qu’elle soit inférieure à 39%, et 1 chance sur 4 qu’elle soit supérieure à 41%. Comme le notait @adelaigue, sur les bords, il semble que la formule ne soit pas valide. Mais si on prend justement les cas limites (probabilité très proche de 0, voire franchement nulle), est-ce qu’on peut dire quelque-chose ?
ou plus précisément, si en interrogeant personnes et que personne n’a déclaré vouloir voter pour A, quelle serait l’intervalle de confiance de la probabilité de voter pour A ? J’avais évoqué ce point dans un très vieux billet. La première piste est la réponse apporté par Pearson-Clopper (Pearson & Clopper (1934)), basée sur l’utilisation de la loi binomiale, et non plus d’une approximation.
library(Hmisc) bornesup=function(N){ binom.test(x=0, n=N, p = 0, alternative="less")$conf.int[2]*100 } BS=Vectorize(bornesup) vn=seq(10,500,by=10) plot(vn,BS(vn))
Une alternative est de faire du bayésien. On utilise le fait que la loi conjuguée de la loi binomiale est la loi beta (comme évoqué dans des précédents billets). Si on suppose que la loi a priori de est une loi alors la loi a posteriori, sachant que sur lancés, on n’a observé que des 0, devient
bornesup=function(N){ qbeta(.95,1,1+N)*100 } BS=Vectorize(bornesup) lines(vn,BS(vn),lwd=2,col="blue")
La courbe ci-dessus correspond au calcul dans le cas du modèle de Pearson-Clopper, en rouge, et avec le quantile de la loi Beta en bleu (qui sont confondues sur ce graphique). Autrement dit, avec 300 personnes interrogées, si aucune ne dit vouloir voter pour A, l’intervalle de confiance pour est [0%;1%] (environ). Ca sera notre second graphique.
c’était déjà sur cette question que je concluais mon précédant billet, mais il est possible d’exploiter davantage cette idée de modèle bayésien avec un a priori suivant une loi beta (comme au dessus) ou de Dirichlet (si on a plus de deux alternatives). La première chose que l’on peut noter, c’est que s’il y a trois alternatives (A, B et C), et que la probabilité de voter pour A est 40%, la probabilité que A gagne… dépend des autres probabilités. On ne peut pas faire comme si A était seul, contre tous les autres:
– si les probabilités de voter pour B et C sont de l’ordre de 30%, on se dit que A a de fortes chances de gagner (oui, pour l’instant ce n’est pas très quantitatif).
– si les probabilités de voter pour B et C sont respectivement 55% et 5%, alors on se dit que A a de faibles chances de gagner.
Plus précisément, si on souhaite savoir si A va l’emporter ou pas, et que 200 personnes ont été interrogées, on peut essayer de faire une abaque, avec des courbes de niveau… dans le simplexe. Pour rappel, les scores obtenus par A, B et C se lisent sur le triangle ci-dessous (le point représenté correspond au triplet (40%,35%,25%))
Les probabilités sont les suivantes, avec n=200 (personnes qui se seraient exprimés)
où la région rouge correspond à une (quasi) certitude de gagner, pour A, et la région blanche, au contraire, une probabilité infinitésimale de gagner.
vx=v[i] vz=v[j] vy=1-vz-vx alpha=c(vx,vy,vz) alpha=c(alpha,1-sum(alpha))*200 RD=rdirichlet(ns, alpha) PROBA[i,j]=mean(RD[,1]>=apply( cbind(RD[,3],RD[,2]),1,max))
Pour ceux que ne sont pas très à l’aise avec les excursions dans le simplexe, on peut se ramener dans un repère orthonormé. Ici, les scores obtenus par A et B se lisent sur le triangle ci-dessous (le point représenté correspond au triplet (40%,35%,25%), le même qu’auparavant
On a alors les probabilités de gagner suivantes, avec à droite la représentation sur le simplexe, en zoomant autour de notre triplet,
|
(ça sera le troisième et dernier graph). Si on regarde plus en détails, l’intervalle de confiance sur la probabilité de gagner dépend de la région où on se trouve dans le simplexe (un peu comme l’intervalle de confiance dépendait de la probabilité de l’emporter, dans la première partie de ce billet).
Maintenant, comme le rappelait @Christian dans un commentaire, tout cela est la version très simpliste de la pratique des sondages. Dans les vrais sondages (et pas des exercices de probabilités comme on vient d’en faire ici), on utilise des quotas, on fait du redessement, et c’est beaucoup plus complexe (enfin… c’est ce que j’essaye de me dire, pour me convaincre que la mauvaise note que j’avais obtenue lorsque j’étais étudiant était motivée).
Ce qui est pénible avec les paradoxes, en tant que scientifique, c’est qu’ils nous obligent à remettre en cause notre intuition, sans cesse. Car mine de rien, on fonctionne tous à l’intuition, c’est un peu comme ça qu’on monte un programme de recherche: on a l’intuition que l’on pourra obtenir des résultats, mais tout reste à prouver.
Prenons la question suivante “les garçons on-t-ils plus de sœurs que les filles ?” en moyenne. N’importe qui avec un peu de bon sens répondrait “oui“, l’idée naturelle est que la fille ne peut pas se compter elle même comme sœur, donc les filles devraient en compter un peu moins… Et comme toujours la formulation est importante. On ne demande pas, par exemple “les garçons ont-t-ils plus chance d’avoir au moins une sœur que les filles ?” (auquel cas la réponse serait différente).
Maintenant réfléchissons un peu, et faisons un petit modèle où on suppose que la taille de la fratrie ne dépend pas des sexes de la fratrie. Et comptons, pour mieux voir ce qui se passe. Et commençons par un cas simple, par exemple avec trois enfants,
Si on récapitule pour les garçons, on a 3 cas avec 2 garçons ayant 1 sœur et 3 cas avec 1 garçon ayant 2 sœurs, soit (à un facteur de normalisation près) un nombre espéré de sœurs de 12. Et pour les filles, on a 3 cas avec 2 filles ayant 1 sœur et 1 cas avec 3 filles ayant 2 sœurs, soit (au même facteur de normalisation près) un nombre espéré de sœur de 12.
Posons la taille de la fratrie, que l’on va supposer fixé. Essayons de formaliser les calculs qu’on vient de faire. S’il y a filles, et garçons,
Bref, dans le calcul de l’espérance pour les garçons, on doit avoir des et pour les filles, . Et en faisant intervenir le nombre de combinaisons possibles de filles, et garçons, et par qui était notre facteur de proportionnalité de tout à l’heure (i.e. le nombre de possibilités d’agencement des sexes dans les familles avec enfants, quand on se focalise sur le sexe de tous, sauf un, qui est choisi au hasard, et qui est la personne qui compte ses soeurs), et en notant l’espérance du nombre de sœurs chez les garçons, et chez les filles, on a (merci à Pierre Antoine qui m’a fait remarqué une erreur dans mes calculs initiaux)
alors que
Un peu de calculs permet de montrer que ces deux quantités sont en fait égales… et valent
ce qui a du sens car il reste -1 personnes dans la fratrie (je peux renvoyer au livre de Ruma Falk, Understanding Probability and Statistics: A Book of Problems, page 46, qui passe par les fonctions génératrices pour une démonstration propre). Amusant, non ? En tous les cas suffisamment pour attirer mon attention alors que j’essaye de récupérer d’un soucis de santé, en dormant dans le lit de ma (grande) fille au milieu des peluches et des poupées…
Hier soir, nous sommes allé mangé au restaurant, en famille, tous les 5. Le serveur nous a placé à une table ronde. Tout le monde s’est installé. Alors que nous commencions l’apéro, ma fille s’est exclamé “ça alors, c’est dingue, on est rangé par age…“. Vraiment ? c’est une coïncidence exceptionnelle ?… allons, quelle est la probabilité que si 5 personnes se mettre autour d’une table ronde, les personnes soient ordonnées par age ? Je vais laisser un peu de temps pour trouver la solution…
Personnellement je dirais 1 chance sur 12….
oui, car 12, c’est 2 fois 3 fois 2..
Non ? bon, la réponse en images alors…
Bon, et je remets en jeu le titre de Sheldon d’or au premier qui reconnaît le morceau…
En mars 1996, Parade Magazine publiait une question intéressante relative à un calcul de probabilité (la question est reproduite ici, ainsi qu’une discussion autour de l’existence de ce professeur… mais comme le note le site là, il semble qu’il s’agisse d’une légende urbaine). De manière un peu surprenante, la réaction de la personne est que – si les deux élèves ont menti – la probabilité que deux personnes désigne le même pneu d’une voiture serait de 1/16. Pourtant, si on note A, B, C et D les quatre pneus de la voiture, un dénombrement rapide donne 16 situations possibles,
mais 4 paires concordantes. Soit un chance sur quatre que les étudiants indiquent le même pneu. Intuitivement, on demande à un étudiant quel pneu a crevé (jusque là pas de calculs de probabilités) puis on veut voir si l’autre indique le même pneu (s’il ment, il a alors une chance sur quatre de donner le même que son voisin). Cette intuition que la probabilité serait 1/16 est alors un peu surprenante. Cela correspondrait au cas où le professeur sait quel pneu avait éclaté, mais pas les étudiants (on retrouve le mythe que l’enseignant sait toujours tout). Ou au cas où trois étudiants étaient dans la voiture et que l’on demande que les trois répondent la même chose.
Considérons un simple jeu de pile ou face. On suppose que l’on tombe sur pile avec une probabilité (histoire de formaliser un peu). Le but du jeu est dedeviner sur quelle face la pièce va tomber: si on devine bien, on gagne, sinon on perd. Et le but de mon billet est de comprendre quelle est la stratégie qu’il convient de mettre en œuvrer si on connaît, ou si on estime .
Supposons que l’on connaisse parfaitement . A priori deux stratégies sont possibles:
Notons la variable correspondant à la valeur prise au lancer, disons 1 si on tomber sur pile et 0 sur face. Et notons la valeur prédite. On gagne si , et la probabilité de gagner est alors
i.e. pour la première stratégie
En revanche pour la seconde, la probabilité de gagner est
Graphiquement, on peut représenter en ligne la prédiction que l’on fait, et en colonne la face obtenue. En bleu, on gagne, et en rouge on perd,
ou encore
i.e.
(avec une inégalité stricte si est strictement inférieur à 1). Bref, il est toujoursplus intéressant de parier toujours la même chose (et de ne pas essayer de jouer au plus malin).
Supposons maintenant que la probabilité soit inconnue, mais estimable. Autrement dit, on a assisté à lancers, on a un échantillon et on estime la probabilité de tomber sur pile par
Autrement dit, si on estime que la probabilité de tomber sur pile dépasse 50%, on prédit pile, et sinon on prédit face.
On gagne si , et la probabilité de gagner est alors
Dans le premier cas, la probabilité de gagner est
Graphiquement, on obtient, cette fois,
Si on se place au bout de 20 lancers, on obtient la courbe rouge ci-dessous pour le cas d’information imparfaite, et verte dans le cas de l’information parfaite, pour le ratio de la probabilité de gagner avec le second cas divisé par le premier cas (nous avions vu auparavant que le ratio été inférieur à 1)
à partir du code suivant
Supposons que l’on répète le jeu un certain nombre de fois, avec un apprentissage au fur et à mesure que les tirages avancent, i.e. l’algorithme s’écrit
> n=10 > p=.64 > X=sample(c("a","b"),size=100,prob=c(p,1-p),replace=TRUE) > S1=S2=rep(NA,100) > for(s in (n+1):100){ + P=c(mean(X[1:(s-1)]=="a"),mean(X[1:(s-1)]=="b")) + S1[s]=(sample(c("a","b"),size=1,prob=P)==X[s]) + S2[s]=(names(sort(table(X[1:s])))[2]==X[s]) + } > sum(S1[(n+1):100])/sum(S2[(n+1):100])
Graphiquement, si on suppose que l’on apprend sur 10 tirages (et que l’on calcule le nombre moyen de victoires sur les 90 lancers suivants), on obtient un graphique assez proche de ce que l’on avait auparavant,
I.e. comme le ratio est inférieur à 1, la seconde stratégie l’emporte toujours.
Moralité, j’ai un peu l’impression qu’il est souvent inutile d’élaborer une stratégie trop complexe, car les stratégies les plus simples permettent de gagner le plus souvent…