Tag Archives: probabilité

A la recherche des homonymes…

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.

Les anniversaires de vos amis sur Facebook

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 n 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’analyses 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.

Moments et variables aléatoires

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)

Examen final ACT2121

L’examen de lundi est en ligne ici, avec des éléments de correction  (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.

Ex-æquo aux élections ?

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 !

Probabilité et contre-intuition

Je suis tombé un peu par hasard sur un vieux papier de Peter Enis, Enis (1973), qui revient sur une relation que l’on prend toujours pour acquis (moi le premier) à savoir http://freakonometrics.blog.free.fr/public/perso5/proba07.gif. Mais il faut des hypothèses complémentaires pour que cette relation soit vraie. Le fait que http://freakonometrics.blog.free.fr/public/perso5/proba09.gif existe par exemple ne suffit pas ! (je laisse les amateurs de résultats contre-intuitifs le temps de réfléchir à cette affirmation)

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

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

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

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

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

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

Sondages et probabilités en trois graphiques

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é).

  • “Il y a 1 chance sur 4 qu’un sondage donne une variation d’un point en + ou en – sur une proportion inchangée”

comme le disait @adelaigue (sur Twitter). Faire un sondage politique (si on résume simplement) c’est estimer la probabilité https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic01.gif de voter pour A, sachant que l’on a alors une probabilité https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic02.gif pour qu’une personne prise au hasard ne vote par pour A. Sur https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic03.gif personnes interrogées, le nombre de personnes qui déclarera voter pour A est https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic04.gif, qui suit une loi binomiale de paramètres https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic03.gif et https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic01.gif. Or par approximation (ou convergence) de la loi binomiale par une loi Gaussienne, on peut dire que

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic05.gif

Or l’estimateur naturel de https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic01.gif est https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic06.gif, i.e.

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic07.gif

Aussi, on en déduit l’intervalle de confiance classique, tel que https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic01.gif appartient à cet intervalle avec une probabilité https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic11.gif,

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic08.gif

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic12.gif 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,

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic09.gif

Les bornes de l’intervalle de confiance de niveau https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic11.gif sont alors

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic10.gif

C’est ce qu’on peut visualiser sur la figure ci-dessous,

Comme le notait @adelaigue, si on souhaite avoir

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic13.gif

(car on parle d’un point de base 100), il faudrait que

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic14.gif

Classiquement, https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic03.gif est aux alentours de 1000. Et ce qu’on cherche c’est une probabilité d’appartenir à un intervalle, i.e. https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic11.gif,

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adi15.gif

Graphiquement, en fonction de https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic01.gif 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 ?

  • peut-on dire quelque chose (malgré tout) quand on ne peut rien dire ?

ou plus précisément, si en interrogeant https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic03.gif 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 https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic01.gifest une loi alors la loi a posteriori, sachant que sur https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic03.gif 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 https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/adic01.gif est [0%;1%] (environ). Ca sera notre second graphique.

  • si on a trois alternatives, quelles sont les vraies probabilités de gagner ?

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.

Ces résultats ont été obtenus par simulation,
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).

Les garçons ont-ils plus de sœurs que les filles ?

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,

  • (G,G,G) – personne n’a de sœur
  • (G,G,F), (G,F,G) ou (F,G,G) – la fille n’a pas de sœur, par contre, il y a deux garçons avec une sœur
  • (G,F,F), (F,G,F) ou (F,F,G) – dans chaque cas, il y a deux filles avec une sœur, et un garçon avec deux sœurs,
  • (F,F,F) – ce qui fait trois filles avec chacune deux sœurs.

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 http://freakonometrics.blog.free.fr/public/perso4/soeur07.gif 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 http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif filles, et http://freakonometrics.blog.free.fr/public/perso4/soeur11.gif garçons,

  • http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif filles et http://freakonometrics.blog.free.fr/public/perso4/soeur11.gif garçons – chacun des http://freakonometrics.blog.free.fr/public/perso4/soeur11.gif garçons a http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif sœurs, alors que chacune des http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif filles a http://freakonometrics.blog.free.fr/public/perso4/soeur13.gif sœurs,

Bref, dans le calcul de l’espérance pour les garçons, on doit avoir des http://freakonometrics.blog.free.fr/public/perso4/soeur12.gif et pour les filles, http://freakonometrics.blog.free.fr/public/perso4/soeur14.gif. Et en faisant intervenir le nombre de combinaisons possibles de http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif filles, et http://freakonometrics.blog.free.fr/public/perso4/soeur11.gif 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 http://freakonometrics.blog.free.fr/public/perso4/soeur07.gif 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 http://freakonometrics.blog.free.fr/public/perso4/soeur06.gif l’espérance du nombre de sœurs chez les garçons, et http://freakonometrics.blog.free.fr/public/perso4/soeur05.gif 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 http://freakonometrics.blog.free.fr/public/perso4/soeur07.gif-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…

On tourne en rond

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…

Petit calcul de probabilités

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

(A,A) (A,B) (A,C) (A,D)
(B,A) (B,B) (B,C) (B,D)
(C,A) (C,B) (C,C) (C,D)
(D,A) (D,B) (D,C) (D,D)

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) http://freakonometrics.blog.free.fr/public/perso3/101_3.jpgpuis 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.

Que signifie “utiliser de l’information” ?

Considérons un simple jeu de pile ou face. On suppose que l’on tombe sur pile avec une probabilité http://freakonometrics.blog.free.fr/public/perso3/pf-20.gif (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 http://freakonometrics.blog.free.fr/public/perso3/pf-20.gif.

  • Connaissance parfaite des probabilités

Supposons que l’on connaisse parfaitement http://freakonometrics.blog.free.fr/public/perso3/pf-20.gif. A priori deux stratégies sont possibles:

  • jouer au hasard, en utilisant les différentes probabilités: on prédit pile avec probabilité http://freakonometrics.blog.free.fr/public/perso3/pf-20.gif et face avec probabilité http://freakonometrics.blog.free.fr/public/perso3/pf-21.gif.
  • ne pas jouer au hasard, et prédire la valeur la plus probable, autrement dit si http://freakonometrics.blog.free.fr/public/perso3/pf-22.gif, on prédit pile, sinon on prédit face (à chaque lancer si on fait plusieurs lancers consécutifs).

Notons http://freakonometrics.blog.free.fr/public/perso3/pf-24.gif la variable correspondant à la valeur prise au lancer, disons 1 si on tomber sur pile et 0 sur face. Et notons http://freakonometrics.blog.free.fr/public/perso3/pf-23.gif la valeur prédite. On gagne si http://freakonometrics.blog.free.fr/public/perso3/pf-25.gif, et la probabilité de gagner est alors

http://freakonometrics.blog.free.fr/public/perso3/pf-26.gif

i.e. pour la première stratégie

http://freakonometrics.blog.free.fr/public/perso3/pf-27.gif

En revanche pour la seconde, la probabilité de gagner est

http://freakonometrics.blog.free.fr/public/perso3/pf-28.gif

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,

Supposons que http://freakonometrics.blog.free.fr/public/perso3/pf-29.gif. Considérons le polynôme http://freakonometrics.blog.free.fr/public/perso3/pf-30.gif, dont les racines sont http://freakonometrics.blog.free.fr/public/perso3/pf-31.gif et http://freakonometrics.blog.free.fr/public/perso3/pf-32.gif, i.e sur http://freakonometrics.blog.free.fr/public/perso3/pf-33.gif,
http://freakonometrics.blog.free.fr/public/perso3/pf-34.gif

ou encore

http://freakonometrics.blog.free.fr/public/perso3/pf-35.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso3/pf-36.gif

(avec une inégalité stricte si http://freakonometrics.blog.free.fr/public/perso3/pf-20.gif 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).

  • Connaissance imparfaite des probabilités

Supposons maintenant que la probabilité http://freakonometrics.blog.free.fr/public/perso3/pf-20.gif soit inconnue, mais estimable. Autrement dit, on a assisté à http://freakonometrics.blog.free.fr/public/perso3/pf-03.gif lancers, on a un échantillon http://freakonometrics.blog.free.fr/public/perso3/pf-01.gif et on estime la probabilité de tomber sur pile par

http://freakonometrics.blog.free.fr/public/perso3/pf-05.gif

Comme auparavant pour la première stratégie, on utilise un générateur de nombre aléatoire http://freakonometrics.blog.free.fr/public/perso3/pf-04.gif (qui tire une loi uniforme sur le segment unité), pour la prochaine valeur, on prédit

http://freakonometrics.blog.free.fr/public/perso3/pf-06.gif

Pour la seconde stratégie

http://freakonometrics.blog.free.fr/public/perso3/pf-07.gif

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 http://freakonometrics.blog.free.fr/public/perso3/pf-08.gif, et la probabilité de gagner est alors

http://freakonometrics.blog.free.fr/public/perso3/pf-09.gif

Dans le premier cas, la probabilité de gagner est

http://freakonometrics.blog.free.fr/public/perso3/pf-11.gif
dont l’espérance est http://freakonometrics.blog.free.fr/public/perso3/pf-12.gif (comme auparavant). En revanche, dans le second cas, le calcul de la probabilité de gagner est

http://freakonometrics.blog.free.fr/public/perso3/pf-13.gif

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

> p=seq(.5,1,by=.01)
> q=(p^2+(1-p)^2)/(p*(1-pbinom(10,20,p))+
+ (1-p)*pbinom(10,20,p))
> plot(p,q)
  • Répétons le jeu plusieurs fois et apprentissage des probabilités

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…

Mal calculer les probabilités tue les probabilités

Suite à un tweet d’Olivier, j’ai découvert avec stupéfaction l’article paru dans Libé sur le risque nucléaire (ici).

On apprend que la probabilité d’accident majeur est d’environ 0.0003 0,0003 par an pour chaque réacteur. Or l’Europe possède 143 réacteurs et donc, sur 30 ans, la probabilité d’un incident majeur est

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/nucle-05.gif

et les auteurs de l’article concluent que “la probabilité d’occurrence d’un accident majeur sur ces parcs serait donc […] de plus de 100% pour l’Union européenne”

Etienne Ghys a déjà – brillamment – discuté ce chiffre avant hier (ici). Le soucis est que si https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/nucle-02.gif est la proba annuelle d’incident par an et par réacteur, en supposant les centrales indépendantes, la probabilité d’avoir un incident sur https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/nucle-3.gif réacteurs sur une durée de https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/nucle-04.gif années est

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/nucle-01.gif

Comme https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/nucle-02.gif est supposé petit, un développement limité donnerait

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/nucle-06.gif

Malheureusement, dans les développements limités de ce genre, s’il faut que https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/nucle-02.gif soit petit, encore faut-il que la puissance ne soit pas trop grande….
En fait, si on retient que pour un réacteur nucléaire à eau pressurisée (REP) tels ceux exploités en Europe de l’Ouest, le risque de fusion du cœur est estimé à 0.00005 par centrale et par an (en reprenant les estimateurs mentionnés sur wikipedia, ici), la vraie probabilité d’avoir un incident majeur en Europe serait plutôt de l’ordre de

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/nuclea-07.gif

(avec une probabilité 10 fois supérieure, Etienne arrivait à 72%). C’est marrant car cette erreur commise en croyant qu’un développement limité serait légitime quand on utilise des probabilité annuelle sur des horizons plus longs était évoqué ici il y a un an jour pour jour, sur des problèmes d’actualisation….

Trop de probabilités tue les probabilités

Hier matin, au moment de quitter la maison, j’ai jeté un œil à la météo (car à Montréal le temps est très changeant, pas seulement d’un jour sur l’autre, comme je l’avis mentionné en arrivant en septembre, ici). En voyant les gros nuages gris dans le ciel, en me réveillant, je me suis demandé s’il fallait mettre un manteau et des bottes en caoutchouc aux enfants (pour aller l’école). Sur le site de météo, on a en effet la probabilité de pluie, i.e. de P.D.P.

Si on suppose que les probabilités sont indépendantes d’heure en heure, la probabilité qu’il va pleuvoir dans la matinée (dans les 5 prochaines heures pour faire simple, ou plus concret), va s’écrire comme le complémentaire de la probabilité qu’il va faire beau toute la matinée, autrement dit, avec des probabilités de pluie de l’ordre de 25%, par heure, la probabilité qu’il pleuve pendant les prochaines heures est de 75% environ

> 1-((1-.25))^5
[1] 0.7626953

Bref, on est loin de 25%. Essayons d’aller un peu plus loin pour mieux comprendre…
Notons https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie03.gif la variable indicatrice indiquant s’il pleut, ou pas (1 s’il pleut, et 0 sinon) à l’instant https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie04.gif. Supposons qu’il existe un processus latent sous-jacent, https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie05.gif tel que

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie06.gif

 

A une date https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie04.gif, on se demande ici s’il va pleuvoir dans la matinée, i.e. on veut

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie07.gifqui peut encore s’écrire

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie08.gifou, ce qui sera pleut être plus utile ensuite

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie09.gifAvec notre modèle latent, cette probabilité peut s’écrire

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie10.gifc’est à dire en notant https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie13.gif la fonction de répartition (jointe) du vecteur aléatoire https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie11.gif cette probabilité s’écrit
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie12.gifOr Wassily Hoeffding a montré qu’il était possible d’encadrer une fonction de répartition, avec les bornes suivantes: pour tout https://latex.codecogs.com/gif.latex?\boldsymbol{x}\in\mathbb{R}^d,
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie14.gif

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie15.gif
et
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie16.gif
(en notant que ces bornes sont atteignables). Autrement dit, ici, la probabilité qu’il fasse beau 5 heures de suite est encadrée par la probabilité qu’il pleuve sur une heure… et 1.
Bref, on n’a pas grand chose à dire*.

Une possibilité est de supposer un modèle Gaussien pourhttps://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie05.gif. Supposer les lois marginales n’est en rien contraignant. En revanche, supposer que la dynamique sous-jacente est Gaussienne est une hypothèse très forte. Et pour faire simple, supposons même que la dynamique soit auto-régressive (à l’ordre 1), i.e.
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie18.gifoù https://latex.codecogs.com/gif.latex?(\eta_t) est un bruit blanc, indépendant du processus latent. L’intérêt est que l’on a spécifié ainsi la structure de la matrice de variance covariance. En particulier, la matrice de corrélation du vecteur https://latex.codecogs.com/gif.latex?\boldsymbol{X}^\star est alors
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie01.gif
On peut d’ailleurs visualiser facilement la forme de cette matrice en fonction de https://latex.codecogs.com/gif.latex?\rho,
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/animationcorrelationpluie.gif
Avec ce formalisme, on peut alors en déduire https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie12.gif car on connaît très bien la fonction de répartition d’un vecteur gaussien multivarié via

On peut alors étudier l’impact de http://freakonometrics.blog.f ree.fr/public/perso3/slatex.gif (qui est associée à la probabilité d’avoir de la pluie pendant une heure en moyenne) et de la corrélation https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/correlation-pluie21.gif (qui est associée à la dynamique du processus). Pour simplifier, si on prend https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/slatex.gif tel que la probabilité de pluie soit de l’ordre de 25%, on a la probabilité d’avoir de la pluie suivante,

> s = .25
> r = .85
> i = 1:n
> I = matrix(rep(i,each=n),n,n)
> J = matrix(rep(i,n),n,n) 
> C = r^(abs(I-J))
> 1-pmnorm(rep(qnorm(1-s),n),mean=rep(0,n),varcov=C)[1]
[1] 0.458074

i.e. en faisant une boucle sur la valeur de la corrélation

Et plus généralement si on croise le seuil (i.e. la probabilité horaire de pleuvoir) et la probabilité, on a l’abaque suivante

On pourrait toutefois objecter que l’important n’est pas forcément qu’il ne pleuve pas de la matinée, mais qu’il ne pleuve pas quand ils sont dehors, i.e. (pour simplifier), à 8 heures, midi et 16 heures (c’est à dire toutes les 4 heures). On s’intéresse alors à
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/pluiecorrel-31.gif
qui est toujours un vecteur Gaussien, mais avec cette fois une structure de corrélation un peu différente,
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/ppluie-coreel.gif
Dans ce cas, la probabilité de ne pas avoir de plus à ces trois moments de la journée est alors de la forme suivante

La corrélation joue alors un rôle relativement faible car il intervient à la puissance 4 dans la matrice d’auto-corrélation.
Dans ce cas, on notera que l’hypothèse d’indépendance n’est peut être pas absurde. Et avec une P.D.P. de 30%, sur toute la journée, on peut probablement se dire que la probabilité d’avoir de la pluie soit le matin, soit le midi, soit le soir sera de l’ordre de

> 1-(1-.3)^3
[1] 0.657

A tester !
* En fait si, on pourrait dire beaucoup de choses (et écrire une dizaine de billets) en particulier sur la borne supérieure, ici 1, dès lors que la probabilité horaire d’avoir de la pluie dépasse 20% – car la borne supérieure est

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/prooooobacommentaire.gif

En fait, la borne inférieur n’est pas une fonction de répartition de manière générale. Toutefois, la borne est atteignable pour une certaine structure de dépendance. Ce qui explique que par la suite, avec un modèle Gaussien, on sera relativement loin de cette borne supérieure. La borne inférieure (i.e. la borne supérieure sur les fonction de répartition) est elle atteinte dès lors que les risques sont variables sont comonotones, c’est à dire dans le cas Gaussien si la corrélation entre deux instants (et entre tous les instants) vaut 1. La petite animation sur la matrice de corrélation permet de mieux comprendre la difficulté de comprendre le cas où la corrélation rho est négative: dans le cas où https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/rhofleche.gif1, on a une corrélation proche de 1 partout dans la matrice. En revanche, si https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/rhofleche.gif -1, on a des corrélations proches de -1 dans la moitié (environ) des cellules de la matrice, et +1 dans l’autre moitié. La borne supérieure sur la probabilité savoir de la pluie dans matinée correspond donc à une structure plus complexe que la borne inférieure, car on ne peut pas avoir une corrélation négative forte avec tous les lags.