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

Les doctorants, le sexe et le directeur de thèse

Pourquoi certains chercheurs publient beaucoup, et d’autres moins ? ou de l’hétérogénéité dans une loi de Poisson (oui, le but principal du billet est de parler de la régression de Poisson).

Il y a une vingtaine d’années, Scott Long publiait un article intitulé ‘The Origins of Sex Differences in Science‘. Mais surtout, les données sont téléchargeables ici.

> library(foreign)
> base=read.dta("http://www.ats.ucla.edu/stat/stata/examples/long/couart2.dta")
> head(base)
  art   fem     mar kid5  phd ment
1   0   Men Married    0 2.52    7
2   0 Women  Single    0 2.05    6
3   0 Women  Single    0 3.75    6
4   0   Men Married    1 1.18    3
5   0 Women  Single    0 3.75   26
6   0 Women Married    2 3.59    2

On dispose de plusieurs informations, recueillies auprès de 950 docteurs en biochimie. On a ainsi

  • le nombre de publications au cours des 3 dernières années de leur doctorat,
  • le sexe du docteur,
  • son statut marital (marié(e) ou célibataire), a priori à la fin du doctorat,
  • le nombre d’enfant(s) du docteur agé(s) de moins de 5 ans, là encore à la fin du doctorat
  • le nombre de publications du directeur de thèse (appelé ici mentor) au cours des 3 dernières années

(on a aussi le prestige du département, mais je trouve cette variable trop subjective). Comme on fait un comptage de publications, on pourrait espérer voir une loi de Poisson.

> n=nrow(base)
> lambda=mean(base$art)
> compte=table(base$art)
> N=as.numeric(names(compte))
> empirique=as.numeric(compte[as.character(0:19)])/n
> theorique=dpois(0:19,lambda)
> cbind(0:19,empirique,theorique)
        empirique theorique
 [1,]    0 0.3005  0.1839
 [2,]    1 0.2688    0.3114
 [3,]    2 0.1945    0.2636
 [4,]    3 0.0918    0.1487
 [5,]    4 0.0732  0.0629
 [6,]    5 0.0295  0.0213
 [7,]    6 0.0185  0.0060
 [8,]    7 0.0131  0.0014
 [9,]    8 0.0010  0.0003
[10,]    9 0.0021  0.0000
[11,]   10 0.0010  0.0000
[12,]   11 0.0010  0.0000
[13,]   12 0.0021  0.0000
[14,]   13     NA  0.0000
[15,]   14     NA  0.0000
[16,]   15     NA  0.0000
[17,]   16 0.0010  0.0000
[18,]   17     NA  0.0000
[19,]   18     NA  0.0000
[20,]   19 0.0010  0.0000
> plot(0:19,theorique,col="red",type="p")
> lines(0:19,empirique,type="p",col="blue ")

Mais non… graphiquement, on est relativement éloignés de la loi de Poisson

avec en rouge la fréquence que l’on aurait eu avec une loi de Poisson (de paramètre la moyenne obtenue sur l’ensemble de la population), et en bleu, la fréquence empirique.
Un test rapide permettra de convaincre les sceptiques,

> library(vcd)
> summary(goodfit(base$art,type="poisson"))

	 Goodness-of-fit test for poisson distribution

                      X^2 df     P(> X^2)
Likelihood Ratio 296.3715 13 1.381259e-55

Manifestement, il y a de l’hétérogénéité dans la population des doctorants. Il y a ceux qui publient beaucoup, ceux qui publient peu et ceux qui ne publient pas.

  • Le sexe et le nombre de publications

Peut-être qu’en discriminant par le sexe, on aurait, par sous-populations, des lois de Poisson. Si on étudie la distribution du nombre de publications pour les hommes

versus la distribution du nombre de publications pour les femmes

ça ne semble pas suffire pour expliquer toute l’hétérogénéité. Même s’il semble y avoir une différence ne serait-ce que sur le nombre moyen de publications par sexe,

> tapply(base$art,base$fem,mean)
     Men    Women 
1.882591 1.470309

Là encore, on peut faire un test rapide de comparaison de moyennes

> t.test(base$art~base$fem)

	Welch Two Sample t-test

data:  base$art by base$fem 
t = 3.3298, df = 885.945, p-value = 0.0009049
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 0.1692780 0.6552866 
sample estimates:
  mean in group Men mean in group Women 
           1.882591            1.470309

Visiblement les deux sous-populations sont différentes, mais le critère du sexe ne suffit pas à expliquer l’hétérogénéité.
Maintenant, notons que le critère retenu est le nombre de publications. La stratégie des personnes de sexe féminin peut être de viser moins de publications, mais dans de meilleures revues. En
tous les cas, le sexe pourrait expliquer en partie (et en partie seulement) l’hétérogénéité observée.

  • Le statut marital et le nombre de publications

On peut aussi comparer le statut marital, en particulier marié(e)

versus célibataire

Là encore, on peut comparer les moyennes, et faire un test d’égalité. Mais cette fois, l’influence est beaucoup plus faible que le sexe.

> t.test(base$art~base$mar)

	Welch Two Sample t-test

data:  base$art by base$mar 
t = -1.1869, df = 710.307, p-value = 0.2356
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.4033993  0.0994165 
sample estimates:
 mean in group Single mean in group Married 
             1.592233              1.744224

la différence entre les moyennes n’étant pas, ici, significative.

  • Les enfants et le nombre de publications

Une variable probablement plus intéressante que le statut marital pourrait être le nombre d’enfants (moi qui finissait la rédaction de ma thèse alors que ma deuxième apprenait à marcher On a la distribution suivante pour le nombre de publications pour un(e) docteur(e) ayant au moins un enfant,

versus un(e) docteur(e) sans enfant,

Mais la encore, le test de comparaison de moyenne n’est pas vraiment concluant,

> base$kid=base$kid5>0
> t.test(base$art~base$kid)

	Welch Two Sample t-test

data:  base$art by base$kid 
t = 0.6136, df = 646.359, p-value = 0.5397
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.1803154  0.3442384 
sample estimates:
mean in group FALSE  mean in group TRUE 
           1.721202            1.639241

Certes, les docteur(e)s sans enfant publient plus, mais ce n’est pas non plus flagrant. On pourrait se dire que ce sont surtout les femmes qui ont eu des enfants qui sont pénalisées (la encore, car l’objectif est d’expliquer une quantité de publications), mais tout d’abord autant les hommes que les femmes répercutent le fait d’avoir eu, ou non, un enfant; mais surtout la variable de sexe est beaucoup plus discriminante que le fait d’avoir eu des enfants,

> baseH=base[base$fem=="Men",]
> baseF=base[base$fem=="Women",]
> t.test(baseH$art~baseH$kid)

	Welch Two Sample t-test

data:  baseH$art by baseH$kid 
t = 1.3831, df = 491.989, p-value = 0.1673
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.1136025  0.6538074 
sample estimates:
mean in group FALSE  mean in group TRUE 
           2.011628            1.741525 

> t.test(baseF$art~baseF$kid)

	Welch Two Sample t-test

data:  baseF$art by baseF$kid 
t = 0.9525, df = 137.906, p-value = 0.3425
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.1764049  0.5043375 
sample estimates:
mean in group FALSE  mean in group TRUE 
           1.501466            1.337500
  • Impact du directeur de thèse sur le nombre de publications

Une fois analysées les variables personnelles, si on regardait l’impact du directeur de thèse sur le travail de l’étudiant ? Qui sait, les variables pourraient être corrélées (je ne parle pas ici de relation causale). Par exemple, on peut regarder la distribution du nombre de publications pour un(e) docteur(e) ayant eu un directeur de thèse qui publiait peu (moins que la moyenne ici)

versus un directeur de thèse qui publiait beaucoup (plus que la moyenne),

Cette fois, on a une nette différence

> base$publis=base$ment>=9
> t.test(base$art~base$publis)

	Welch Two Sample t-test

data:  base$art by base$publis 
t = -6.7568, df = 497.017, p-value = 3.96e-11
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -1.2506429 -0.6871632 
sample estimates:
mean in group FALSE  mean in group TRUE 
           1.341338            2.310241

Un directeur de thèse qui publie semble avoir des étudiant(e)s qui publient aussi (ou réciproquement).

  • Régression de Poisson et nombre de publications

La population était initialement très hétérogène. Trop hétérogène pour une loi de Poisson, avec une variance valant plus du double de la moyenne empirique,

> var(base$art)/mean(base$art)
[1] 2.191358

> summary(glm(art~1,data=base,family=quasipoisson))

Call:
glm(formula = art ~ 1, family = quasipoisson, data = base)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8401  -1.8401  -0.5770   0.2294   7.5677  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.52644    0.03761   14.00   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for quasipoisson family taken to be 2.191389)

    Null deviance: 1817.4  on 914  degrees of freedom
Residual deviance: 1817.4  on 914  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

Si on régresse sur les différentes variables on retrouve les interprétations que nous avions eues,

> summary(glm(art~mar+fem+I(kid5>0)+ment,data=base,family=quasipoisson))

Call:
glm(formula = art ~ mar + fem + I(kid5 > 0) + ment, family = quasipoisson, 
    data = base)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.5080  -1.5615  -0.3626   0.5614   5.4494  

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      0.338732   0.081482   4.157 3.53e-05 ***
marMarried       0.149527   0.085110   1.757  0.07928 .  
femWomen        -0.218643   0.074151  -2.949  0.00327 ** 
I(kid5 > 0)TRUE -0.250286   0.085834  -2.916  0.00363 ** 
ment             0.026000   0.002656   9.788  < 2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for quasipoisson family taken to be 1.837376)

    Null deviance: 1817.4  on 914  degrees of freedom
Residual deviance: 1641.1  on 910  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

On peut d’ailleurs faire une régression GAM avec une estimation par splines de l’effet du nombre de publications du directeur de thèse. Ce qui permet de légitimer le fait que la variable intervient de manière linéaire

On pourrait noter la concavité de la courbe qui semble sous-entendre qu’un directeur qui publie énormément peut toujours s’occuper de ses étudiants – et les faire publier plus qu’un directeur qui publie moins, mais en première approximation, la courbe peut être supposée linéaire.
On peut aussi regarder rapidement qui ne publie pas pendant sa thèse,

> summary(glm((art==0)~mar+fem+I(kid5>0)+ment,data=base,family=binomial))

Call:
glm(formula = (art == 0) ~ mar + fem + I(kid5 > 0) + ment, family = binomial, 
    data = base)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.1609  -0.9043  -0.7019   1.2697   2.5221  

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -0.28843    0.18591  -1.551   0.1208    
marMarried      -0.33565    0.18800  -1.785   0.0742 .  
femWomen         0.23697    0.15829   1.497   0.1344    
I(kid5 > 0)TRUE  0.42957    0.19032   2.257   0.0240 *  
ment            -0.08141    0.01261  -6.455 1.08e-10 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1118.7  on 914  degrees of freedom
Residual deviance: 1052.0  on 910  degrees of freedom
AIC: 1062.0

Number of Fisher Scoring iterations: 5

Certain(e)s ont eu un enfant, et semblent avoir été pénalisé(e)s, mais si un doctorant n’a pas publié, c’est presque uniquement à cause de son directeur (qui publie peu).
Pour reprendre l’interprétation sur l’analyse du nombre de publications, par exemple, le fait d’être marié, ou pas, a peu d’influence sur le nombre de publications. Par contre les femmes – toutes choses étant égales par ailleurs – publient 20% de moins que les hommes,

> exp(glm(art~mar+fem+I(kid5>0)+ment,data=base,
+ family=quasipoisson)$coefficients[3])
 femWomen 
0.8036087

Par contre le fait d’avoir – ou pas – des enfants semble ici beaucoup plus fort, car les docteur(e)s qui ont eu au moins un enfant publient 25% de moins que ceux qui n’en ont pas eu.

> exp(glm(art~mar+fem+I(kid5>0)+ment,data=base,
+ family=quasipoisson)$coefficients[4])
I(kid5 > 0)TRUE 
      0.7785778

Aussi, une femme ayant eu un enfant aura publié 37% de moins qu’un homme qui n’a pas eu d’enfants. Toutes choses étant égales par ailleurs.

  • Toutes choses étant égales par ailleurs ?

C’est ici qu’est l’arnaque, car rien ne garantie que toutes choses puissent être égales par ailleurs. Par exemple, une femme avec enfant (ou pas) peut elle avoir un directeur de thèse qui publie beaucoup ? N’y a-t-il pas un peu d’endogénéité cachée ?
Si on regarde le nombre de publications du directeur de thèse versus le sexe de l’étudiant,

> t.test(base$ment~base$fem)

	Welch Two Sample t-test

data:  base$ment by base$fem 
t = 2.7101, df = 901.2, p-value = 0.006854
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 0.4587075 2.8673525 
sample estimates:
  mean in group Men mean in group Women 
           9.532389            7.869359

on note que les directeurs de thèse des hommes publient plus que les directeurs de thèse des femmes. Autrement dit, si les femmes publient moins, c’est peut être parce que les directeurs de thèses qui pourraient les aider à publier ne les prennent pas en thèse, non ? Car il n’y a pas de lien entre le fait d’avoir ou pas des enfants à la fin de sa thèse, et le nombre de publications de son directeur,

> t.test(base$ment~base$kid)

	Welch Two Sample t-test

data:  base$ment by base$kid 
t = -1.5626, df = 499.06, p-value = 0.1188
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -2.5485856  0.2905184 
sample estimates:
mean in group FALSE  mean in group TRUE 
           8.377295            9.506329

Bref, c’est manifestement le directeur de these qui influence le plus le nombre de publication d’un doctorant. Mais il serait intéressant d’avoir les memes statistiques 5 ans apres la soutenance…

ACT2040, exercices de tarification

Les exercices commenceront ce vendredi, et les énoncés (ainsi que pour quelques uns des pistes de correction) sont en ligne ici. Les exercices (et les corrections) ont été tapés pour ce cours par Jean-Philippe Boucher. La premiere série d’exercices sont des exercices de remise en route sur les calculs d’espérances, et d’espérance conditionnelle. On reviendra dessus mardi prochain.

C’est pas ma faute: c’est l’autre là, devant !

Les vacances sont l’occasion de se poser des questions que l’on n’a pas l’occasion de se poser d’ordinaire. Par exemple quand je m’étais demandé pourquoi, dans un petit village perdu, je continuais de fermer la porte de la voiture a clé (ici), il y a maintenant plusieurs mois.

Mais surtout, les longs trajets en voiture sont une opportunité pour discuter avec le copilote (qui en l’occurrence était la mère de mes enfants), et d’essayer de la convaincre du bien-fondé de certaines de nos décisions. Après une petite heure de route, la semaine dernière, je l’ai entendue m’interpeller “ralentis,  il y un radar et tu est à 115 !”. Bon, je dois avouer que l’argument “non non, je respecte la vitesse, mais de ton point de vue tu dois avoir une vision déformée et biaisée du tableau de bord car si tu étais pile en face, tu verrais que je suis un peu en dessous de 105 km/h” me marche plus…. car elle avait raison, j’étais bien à 115 km/h. Mais instinctivement, ma réponse a été la même que des millions d’automobilistes qui doivent faire face aux foudres de leur copilote: “oui, mais c’est pas ma faute”… Sans y réfléchir, cette excuse plus-que-centenaire, que je sortais a mes parents après avoir cassé la porte d’entrée après avoir donné un coup de pied dedans, ou que me ressortent mes enfants après lancé leur avion Playmobil par la fenêtre, est sortie toute seule de ma bouche. Et au moment de la sortir, j’ai compris immédiatement dans quel pétrin je me mettais… surtout que dans les yeux de ma copilote, je sentais pétiller le goût de la victoire (ou plutôt du je-sens-qu’on-va-rigoler-un-peu-avec-une-histoire-a-dormir-debout). “ah oui ? mais c’est la faute de qui alors… ?”.

Bon, alors la, on commence à visualiser en une fraction de seconde la liste des boucs émissaires possibles (le voisin, les enfants, mes parents, la belle-mère, la petite vieille du coin de la rue, la police, la pluie , le soleil, la fatigue, internet, la société, Nicolas Sarkozy, les normes IFRS, que sais-je encore….), et je ne sais pourquoi, j’ai retenu “ben c’est de la faute de l’autre, là, devant” (encore un 75 tiens !).
C’est vrai quoi, si ça se trouve, avec un effet d’entraînement, il m’aspire et sans que je fasse rien, il me fait gagner 10 km/h. Bon, mes compétences en mécanique des fluides sont limitées. Donc j’aurai du mal a justifier cette théorie (un peu fumeuse, j’en conviens). J’ai personnellement opté pour “quand on conduit. on ne peut pas passer son temps le nez sur le tableau de bord… donc on regarde très occasionnellement… Et la plupart du temps, on se positionne par rapport aux autres, et en particulier la voiture juste devant… on essaye de maintenir une distance à peu près constante, et donc on roule a la même vitesse. Si j’allais trop vite, c’est que lui allait trop vite !”. Et toc… pas mal ma théorie, non ?

Bon, maintenant va falloir le vérifier.

Mais ça tombe bien, car il y a quelque temps, j’avais travaillé sur des données de passage de voitures, dans la banlieue parisienne (ici).
Dans un premier temps, on peut faire un petit test d’indépendance, dichotomique (je dépasse la vitesse autorisée, versus la voiture devant dépasse la vitesse autorisée). Plaçons nous un moment ou ça ne bouchonne pas (disons entre 15 et 16 heures, un jour de semaine) et regardons la fréquence de véhicules en excès de vitesse, ainsi que la même fréquence pour les véhicules immédiatement derrières,

> M
     [,1] [,2]
[1,]  879  104
[2,]  103   36

avec par lignes le véhicule qui suit, et par colonne celui qui précède. Ici, 103 véhicules dépassaient la vitesse autorisée sans que le véhicule qui précédait ne dépasse les 110 km/h. Un test d’indépendance du chi-deux donne ici,

> chisq.test(M)
 
	Pearson's Chi-squared test with Yates' continuity correction
 
data:  M 
X-squared = 24.7862, df = 1, p-value = 6.405e-07
 
 
> chisq.test(M)$expected
         [,1]      [,2]
[1,] 860.3440 122.65597
[2,] 121.6560  17.34403

Bref, il semble que l’on rejette l’hypothèse d’indépendance.

On peut d’ailleurs le faire à plusieurs heures de la journée, et étudier l’évolution de la p-value du test d’indépendance,

Moralité, on rejette indépendance entre les vitesses de véhicules qui se suivent aux heures normales (i.e. en dehors des heures de pointe, et la nuit).

On peut aussi regarder en fonction de la vitesse du premier véhicule, la vitesse du véhicule juste derrière, avec une régression lissée,

Enfin, on peut aussi tenter une régression logistique, afin de quantifier la probabilité de faire un excès de vitesse, conditionnellement a la vitesse du véhicule qui le précède,

Si au lieu de le faire sur une journée, on regarde sur 38 jours (c’est ce dont je dispose dans ma base), on a

Autrement dit, la probabilité que je dépasse la vitesse autorisée (ici 110 km/h) dépend très fortement – et de manière sensiblement croissante – de la vitesse du véhicule qui me précède…

Moralité, mon argument pourrait être valide: si on va vite, c’est en partie parce que la voiture devant va vite, Mais si mon argumentaire me convainc (et pourrait – presque – convaincre ma femme), on verra s’il convainc la police. Car s’ils me contactent pour m’avoir flashé a une vitesse excessive, j’aurais envie de leur sortir ma théorie…. En plus ma théorie est basée sur une étude de corrélation (je fais des régressions toutes simples), et pas de causalité… ma conclusion pourrait alors parfaitement être inversée, ou complétée: je pourrais tout aussi légitimement expliquer à ma copilote “ben si j’allais trop vite, c’est de la faute de l’autre, la, derrière… il me collait, alors j’allais vite”.
Bref, il semble donc que l’excuse du “c’est pas ma faute, c’est l’autre” a de beaux jours devant elle… surtout que les élections approchent en France, donc je pense qu’on n’a pas fini de l’entendre !

Les traders sont des animaux curieux…

Et je dis ça d’autant plus sereinement que j’ai travaillé un temps en salle de marche, que j’ai pas mal d’anciens élèves qui y sont, ainsi que la majorité de mes copains de promo…

On a beaucoup lu sur les cracks en maths qui partaient faire de la finance. Je ne pense pas qu’il y ait besoin d’être un crack en maths pour aller jouer en salle de marche (comme je l’ai dit à plusieurs reprises ici ou , pas mal de problème de finance sont des problèmes de buveurs d’alcool qui ont du mal à rentrer chez eux), mais c’est vrai que lors des entretiens d’embauche, il n’est pas rare de mettre des petits exos de maths. Le but est surtout de se faire une idée sur la manière dont les candidats résolvent un problème, et de vérifier que les candidats peuvent faire preuve d’un peu de logique (cela dit, j’ai fait la même chose lorsque j’ai recruté un collègue lorsque je travaillais à Hong Kong dans une compagnie d’assurance).
J’y pensais l’autre jour en lisant chasseur de tête (Hodejegerne en vo) de Jo Nesbø (sans l’inspecteur Harry Hole), où un des personnages demandait

En fait, la réponse proposée est fausse, ou partiellement fausse… Et si je le sais c’est parce qu’il y a quelques mois, je discutais avec un ancien élève qui me racontait qu’il posait la même question lors des entretiens d’embauche (ou de stages). Plus précisément, ça commence calmement…

« donnez nous deux nombres premiers qui soient deux nombres impairs consécutifs »

Effectivement, pour l’instant, ce n’est pas trop dur…  on peut penser à 3 et 5, ou 5 et 7, 11 et 13, ou encore 17 et 19. Les plus tordus penseront à 31721 et 31723, voire 49919 et 49921. Bref, cette propriété est assez commune chez les nombres premiers, donc c’est facile…

« donnez nous trois nombres premiers qui soient trois nombres impairs consécutifs »

Si on a pensé à 3 et 5 avant, on pensera naturellement à 3, 5 et 7… Et ça marche… Mais ensuite vient la question plus vicieuse…

« pourriez-vous en trouver trois autres »

plus dur déjà… Un indice ? Comme c’est un entretien d’embauche dans une banque (et pas un exercice d’agrégation de mathématiques), soit il y a une réponse simple… soit la réponse est que ça n’est pas possible….

Et effectivement, ça n’est pas possible. Car si on prend trois nombres impairs consécutifs, il y en a toujours un qui est un multiple de 3. Formellement, on a retenu n, n+2 et n+4

Si les deux premiers (n et n+2) sont premiers, alors n=3k-1 et n+2=3k+1. Et donc n+4=3k+3=3(k+1) qui est un multiple de trois…

Si les deux derniers (n+2 et n+4) sont premiers, alors n+2=3k-1 et n+4=3k+1. Et donc n=3k-3=3(k-1) qui est un multiple de trois…

Enfin, si le premier et le dernier sont premiers, alors n=3k-2 et n+4=3k+2. Et donc n+2=3k qui est un multiple de trois… Donc la seule solution est que le multiple de 3 qui doit apparaître soit précisément 3. Donc la solution unique est 3, 5 et 7.

Je sais pas vous… mais moi ça me laisse un petit goût de pas-fini… On aurait envie d’aller plus loin, genre

« donnez nous n nombres premiers qui soient nombres consécutifs, en comptant de p en p »

en notant que, pour l’instant, on a regardé (uniquement) le cas où n vaut 3, et p vaut 2. Et il y a peut-être un résultat plus général derrière, non ?

Comme me l’a fait note Jürgen, la question un peu mieux formulée (en termes mathématiques) serait

« donnez nous une progression arithmétique de longueur n et de raison p formée de nombres premiers »

Intuitivement, on se doute qu’il faut regarder si n et p sont premiers entre eux (car c’est dans ce cas que la réponse risque de ne pas exister) et c’est ce que m’a confirme Agnès. En fait

  • dans le cas où n et p sont premiers entre eux

si n n’est pas premier, alors il n’y a aucune solution. Par exemple (n,p)=(4,2).
En revanche, si n est premier, alors il n’y a qu’un nombre fini de solution (qui peut être nul, ou 1 comme dans la question initiale). L’idée est que le premier terme est à chercher dans l’ensemble des entiers

{n-(n-1)p, n-(n-2)p, …, np, n}

En particulier, n doit être un des termes de la progression.

Par exemple, si (n,p)=(3,2) il y a une seule solution (3,5,7) alors que pour (n,p)=(5,2) il n’y a pas de solution.Et si (n,p)=(5,6) il y a une seule solution avec (5, 11, 17, 23, 29).

  • dans le cas où n et p ne sont pas premiers entre eux

Dans ce cas, la question doit être difficile (toujours selon Agnès), et on conjecture qu’il y a un nombre infini de solutions. Je ne sais pas si ça aide, car je ne sais pas si cet infini est grand ou pas… au sens où il faudrait peut être attendre un temps long infiniment long pour avant de tomber sur une telle suite…
En fait, Jürgen me suggérait de regarder un article de Green et Tao (ici) qui ont montré qu’il existe une infinité de progressions arithmétiques (formées de premiers) de longueur arbitraire. Plus précisément, on conjecture que le nombre de progressions de longueur n et de raison p (formées de nombres premiers) dans l’ensemble des nombres entiers de 1 à m serait de l’ordre de m/log(m)^n. Pour rappel, le théorème des nombres premiers dit qu’il y a m/log(m) nombres premiers entre 1 et m (Hadamard, La Vallée Poussin, 1896, ici).

Intéressant problème, non ? En tous les cas, au moins pour repérer des nerds… pour embaucher des gens qui pourront éviter d’acheter ou de vendre n’importe quoi, je ne sais pas trop le rapport…

Cela dit, pour ceux que ces problèmes amusent, un conteur dans le Lot nous a posé une question que je ne serais pas surpris de voir réapparaître dans les entretiens d’embauche, « Dans un lointain pays, un roi souhaite commander plusieurs bagues en or parfaitement identiques a 10 joailliers (et il peut en demander autant qu’il le souhaite à chacun) de 10 grammes chacune . Si 9 sont parfaitement honnêtes, 1 d’entre eux ne peut s’empêcher d’enlever 1 gramme d’or à chacune des bagues qu’il fabrique… Mais le roi ne sait pas lequel. Comment peut-il le démasquer en faisant une unique pesée..? ».

Avis aux amateurs…

En fait, comme il peut demander autant qu’il le souhaite, il demande 1 bague au premier joaillier, 2 au second, 3 au troisième, …., 10 au dernier. Il pèse alors toutes les bagues ensemble. Si tous étaient honnêtes, il y aurait 550 grammes. Le nombre de grammes qu’il manque dit quel joaillier a triché. S’il manque 3 grammes, c’est que 3 bagues sont trop légères, et c’est donc le troisième joaillier qui a triché.

We keep breaking records ? so what ?… Get statistical perspective….

This summer, we have been told that some financial series broke some records (here, in French)

For instance, the French CAC40 had negative return for 11 consecutive days (which has never been seen, so far).

> library(tseries)
> x<-get.hist.quote("^FCHI")
> Y=x$Close
> Z=diff(log(Y))
> RUN=rle(as.character(Z>=0))$lengths
> n=length(RUN)
> LOSS=RUN[seq(2,n,by=2)]
> GAIN=RUN[seq(1,n,by=2)]
> TG=sort(table(GIN))
> TG[as.character(1:13)]
GAIN
   1    2    3    4    5    6    7    8    9 <NA> <NA> <NA>   13
 645  336  170   72   63   21    7    3    4   NA   NA   NA    1
 
> TL=sort(table(LOSS))
> TL[as.character(1:15)]
LOSS
   1    2    3    4    5    6    7    8    9 <NA>   11 <NA> <NA> 
 664  337  186   68   42   14    5    3    1   NA    1   NA   NA 
 
> TR=sort(table(RUN))
> TR[as.character(1:15)]
RUN
   1    2    3    4    5    6    7    8    9 <NA>   11 <NA>   13 
1309  673  356  140  105   35   12    6    5   NA    1   NA    1

Indeed 11 consecutive days of negative returns is a record. But one should keep in mind the fact that the real records for runs is 13 consecutive days with positivereturns…
But what does that mean ? Can we still assume time independence of log-returns (since today, a lot of financial models are still based on that assumption) ?
Actually. if financial series were time-independence, such a probability, indeed, should be rather small. At least on 11 or 10 runs. Something like

http://freakonometrics.blog.free.fr/public/perso3/cacindp03.gif

(assuming that each day, the probability to observe a negative return is 50%). But maybe not over 25 years (6250 trading days): the probability to observe a sub-sequence of 10 consecutive negative value (with daily probability of one half) over 6250 observations will be much larger. My guess is that is would be

http://freakonometrics.blog.free.fr/public/perso3/cacindp02.gif

where at the numerator we have the number of favourable cases over the total number of cases. At the numerator, the first number the number of cases where the first 10 (at least) are negative, then for the second one, we count the number of cases where the first is positive, then the next 10 (at least) are negative (and then the second is positive and then the next 10 are negative, the third is positive etc). For those interested by more details (and a more general formula on runs), an answer can be found here.
But note that the probability is quite large… So it is not that unlikely to observe such a sequence over 25 years.
A classical idea when looking at time series is to look at the autocorrelation function of the returns,

which might suggest that there is no correlation with past returns. But it should be possible to do more advanced tests.
On the CAC40 series, we can run an independence run test on the latest 100 consecutive days, and look at the p-value,

> library(lawstat)
> u=as.vector(Z[(n-100):n])
> runs.test(u,plot=TRUE)
 
	Runs Test - Two sided
 
data:  u 
Standardized Runs Statistic = -0.4991, p-value = 0.6177

The B’s here are returns lower than the median (almost null, so they might be considered as negative returns). With such a high p-value, we accept the null hypothesis, i.e. time independence.
If we consider a moving-time window

we can see that we accept the assumption of independence most the the time.
Actually, here, the time window is 100 days (+/- 50 days). But it is possible to consider 200 days,

or even 400 days,

So, except if we focus on 2006, it looks like we should reject the idea of time dependence in financial markets.
It is also possible to look more carefully at the distribution of runs, and to compare it with the case of independent samples (here we consider monte carlo generation of sequences having the same size),

> m=length(Z)
> ns=100000
> HIST=matrix(NA,ns,15)
> for(j in 1:ns){
+ XX=sample(c("A","B"),size=m,replace=TRUE)
+ RUNX=rle(as.character(XX))$lengths
+ S=sort(table(RUNX))
+ HIST[j,]=S[as.character(1:15)]
+ }
> meana=function(x){sum(x[is.na(x)==FALSE])/length(x)}
> cbind(TR[as.character(1:15)],apply(HIST,2,meana),
+       round(m/(2^(1+1:15))))
 
     [,1]       [,2] [,3]
1    1309 1305.12304 1305
2     673  652.46513  652
3     356  326.21119  326
4     140  163.05101  163
5     105   81.52366   82
6      35   40.74539   41
7      12   20.38198   20
8       6   10.16383   10
9       5    5.09871    5
10     NA    2.56239    3
11      1    1.26939    1
12     NA    0.63731    1
13      1    0.31815    0
14     NA    0.15812    0
15     NA    0.08013    0

The first column above is the empirical frequency of runs of length 1,2,3, etc. The second one is the average frequencies obtained on random simulation of independent sample. The third one is the theoretical frequency based on a (geometric distribution with mean 1).

Here again, it looks like our time series behave like an independent sample. Here is also a nice paper by Mark Schilling on the longest run of heads.
So it is not that odd to observe such a series of losses on financial markets….

Pauvres américains

Pendant nos vacances dans l’Aubrac avec des amis, Christian avait acheté Libé, et machinalement, j’ai entrepris de le survoler le lendemain matin (en sirotant mon café). Je suis tombé sur le paragraphe suivant qui a retenu mon attention pendant plusieurs jours…

L’auteur n’est pas n’importe qui, puisqu’il s’agit de Kenneth Rogoff (ici), grand spécialiste de l’économie américaine. Relisons la phrase afin de mieux comprendre ce qu’il dit: pour “25% des propriétaires immobiliers aux États-Unis” […] “la valeur de leur maison serait inférieure à leur crédit immobilier” 1. Je me permettrais de réécrire la phrase sous la forme suivante “pour un quart des propriétaires immobiliers américains n’ayant pas fini de rembourser leur crédit, la vente de leur maison ne leur permettrait pas de rembourser leur crédit” (c’est en tous les cas comme ça que je la comprends).
Cette petite phrase pourrait être intéressante, En tous les cas, elle semble importante dans l’argumentation visant à expliquer que les américains sont beaucoup trop endettés2. Mais 25%, en quoi est-ce vraiment exceptionnel ou incroyable voire inquiétant ? C’est quoi le pourcentage acceptable ou normal que l’on s’attendrait à avoir ?
N’ayant pas de statistiques sur le sujet, faisons des calculs.

  • un peu de calcul d’actualisation de crédits

Car intuitivement, si un acheteur emprunte avec un faible apport, et sur une longue durée, son crédit va lui coûter cher, éventuellement plus cher que la maison. Au moins au début. Car avec le temps, la valeur du crédit diminue, alors que le prix de la maison, habituellement, augmente.
Considérons une maison de valeur 1 (a l’achat, histoire de simplifier, et de raisonner en pourcentage pour l’apport initial, par exemple). On dispose d’un capital initial http://freakonometrics.blog.free.fr/public/perso3/credit-03.gif (correspondant à l’apport), on contracte un crédit pour une durée http://freakonometrics.blog.free.fr/public/perso3/credit-02.gif, et on suppose que le taux pour le crédit est http://freakonometrics.blog.free.fr/public/perso3/credit-04.gif, et que le taux d’inflation est http://freakonometrics.blog.free.fr/public/perso3/credit-06.gif (la valeur de la maison peut augmenter dans le temps, mais aussi éventuellement baisser si la valeur de http://freakonometrics.blog.free.fr/public/perso3/credit-06.gif est négative). A la date http://freakonometrics.blog.free.fr/public/perso3/credit-01.gif, à son actif, le propriétaire possède la maison, d’une valeur http://freakonometrics.blog.free.fr/public/perso3/credit-09.gif (ce qu’il touche s’il revend la maison, si l’on oublie les frais associés); et au passif, il doit rembourser à la banque un montant http://freakonometrics.blog.free.fr/public/perso3/credit-08.gif, où http://freakonometrics.blog.free.fr/public/perso3/credit-07.gif est le montant des remboursement effectues tous les ans, i.e. solution de

http://freakonometrics.blog.free.fr/public/perso3/credit-11.gif

Si on veut faire les choses proprement, il faudrait intégrer les frais de notaire (disons 7% de la valeur de la maison), ici notés http://freakonometrics.blog.free.fr/public/perso3/credit-10.gif,

http://freakonometrics.blog.free.fr/public/perso3/credit-12.gif

La valeur de la maison est inférieure à la valeur du crédit si

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

(les frais de notaires étant payés à l’achat comme on l’évoquait auparavant, mais aussi en cas de revente3). On peut faire le calcul facilement, sous R,

valeur = function(t,T,a,r=.05,i=0,delta=.07){
k=(1-a+delta)/sum(1/(1+r)^(0:(T-1)))
s=(1+i)^t
v=(T-t)*k
return(c(s*(1-delta),v))}

Par exemple, si http://freakonometrics.blog.free.fr/public/perso3/credit-01.gif est nul, on compare la valeur du crédit à la valeur de la maison au moment de l’achat. Pour quelqu’un ayant un apport de 25%, prenant un crédit avec 20 échéances (sur 20 ans) que l’on commence à rembourser le jour de la signature, la valeur de son crédit (sur une somme empruntée de 0.75) est de 1.2533, environ, si le taux de crédit est de l’ordre de 5%. C’est plus que la valeur (brute) de la maison (ici 1), voire beaucoup plus que ce que rapporterait la revente la maison, qui rapporterait 0.93, ce qui ne lui permet pas de rembourser le crédit….

> valeur(0,20,.25,.05,0)
[1] 0.9300 1.2533

Sur la même durée, vu après la 5ieme échéance (i.e. au bout de 25% des échéances) avec toujours un apport initial de 25%, la valeur du crédit restant à payer à la banque est de l’ordre de 0.94, c’est à dire à peu de choses près, la valeur de revente de la maison s’il n’y a pas d’inflation (ou de perte de valeur du bien immobilier).

> valeur(5,20,.25,.05,0)
[1] 0.9300000 0.9399846

Autrement dit, dans un monde avec une inflation nulle, avec des cohortes d’acheteurs constantes dans le temps, qui prendraient des crédits de 20 ans avec un apport de 25% de la valeur de la maison, 25% des emprunteurs ont, en moyenne, un crédit à rembourser supérieur à la valeur de revente de leur maison, comme le dit l’article. Cette proportion augment

  • quand les taux d’emprunt augmentent
  • quand la durée des emprunts augmente
  • quand l’apport initial diminue

Mais on peut essayer de visualiser tout ça,

  • visualisation des valeurs du crédit, et de la maison
dessin=fonction(T=20,a=.333,r=.05,i=.02,delta=.07){
S=V=rep(NA,T)
for(j in 1:T){
S[j]=valeur(j-1,T,a,r,i)[1]
V[j]=valeur(j-1,T,a,r,i)[2]}
YL=range(S,V)
plot(1:T-.5,V,type="b",col="red",ylim=YL)
lines(1:T-.5,S,col="blue",type="b")
}

Comme on le voit sur le dessin ci-dessous, la proportion des acheteurs dont la valeur du crédit excède la valeur de la maison est d’environ 20% (même avec un apport non négligeable, ici un tiers, et une inflation non nulle, ici 2%). On le voit sur le graphique ci-dessous, avec en bleu la valeur de la maison, et en rouge, la valeur du crédit,

On peut d’ailleurs faire varier les différents paramètres, comme le taux d’emprunt, avec une baisse (passant de 5% à 3%),

ou avec une hausse (passant de 5% à 7%),

On peut aussi faire varier l’apport initial (passant à 50%),

On peut enfin supprimer l’inflation, et supposer que le prix de la maison n’augmente pas vraiment…

Moralité? 25% semble effectivement important, trop important (pour une économie en bonne santé). Mais il ne faut pas se leurrer, car un pourcentage raisonnable (ou viable) semblerait être davantage aux alentours de 15% que de 0%.

  • du crédit immobilier au crédit automobile

Et cela dit, 25% serait un pourcentage relativement faible si on regardait non pas les crédit immobiliers, mais les crédit automobiles. Car par rapport à la situation précédente, on est dans un cas où les taux sont élevés, et où la valeur du bien ne cesse de se dégrader. Par contre la durée est souvent plus courte. Une déflation de 10% n’est peut-être pas la meilleure modélisation qui soit de la perte de valeur du véhicule, mais en première approximation, ça devrait convenir,,,
Graphiquement, on a

Bref, dans le cas du crédit auto (où l’acheteur achèterait intégralement à crédit), dans une situation normale entre 70% et 80% des acheteurs de voiture à crédit sont dans une situation où la revente de leur voiture ne permettrait pas de rembourser leur crédit… Ne faudrait-il pas s’en inquiéter également ? Acheter à crédit un bien dont la valeur ne cesse de baisser, n’est-ce pas dangereux ?

1 au début de l’été, en discutant avec des couples d’amis, dont deux venaient d’avoir des postes de profs à l’autre bout de la France (et qui devaient revendre leur maison), j’ai été surpris de voir que quand ils parlaient de “ne pas perdre d’argent lors de la revente“, ils valorisaient la maison à partir du prix initial, auquel ils ajoutaient les frais de notaires, mais oubliaient le coût du crédit.
2 je me contenterais ici de discuter ce chiffre de 25%, et non pas de savoir si c’est grave que la revente de la maison ne permette pas de rembourser le crédit.
3 je préfère prendre en compte ces frais, car sinon, comme je l’avais déjà évoqué ici, l’achat d’une maison semble toujours une opération gagnante.

Erratum to the paper on heatwave modelling

An erratum to the paper published on heatwave modelling should appear soon, on http://link.springer.com/article/…

In my article, On the return period of the 2003 heat wave, the sentence “summer of2003 was by far the hottest summer since 1500” is wrongly attributed to ‘Luterbacher et al. (2003)’. My probable intent was to cite Casty et al. (2005) which actually states “the years 1994, 2000, 2002 and particularly 2003 were the warmest since 1500”. Note that a similar statement can also be found in Stott et al. (2004), “The summer of 2003 was probably the hottest in Europe since at latest ad 1500”.

Regarding the statement (in reference to Luterbacher et al., which was actually published in 2004) that “their estimate of the return period of that event is 250 years”, the attribution again is incorrect. Such an estimate can be found in Brown et al. (2005), “best estimate of a 1 in 250 year event” or Lowe et al. (2004) “the same event is likely to occur much more frequently; most likely, once every 250 years”, among others.
I regret these inappropriate attributions and inadequate referencing.
  • Brown S, Stott P, Clark R, Temperature Extremes (2011) the Past and the Future. Hadley Centre for Climate Prediction and Research Working Paper. [online http://stabilisation.metoffice.com/…]
  • Casty C, Wanner H, Luterbacher J, Esper J, Böhm R (2005) Temperature and precipitation variability in the European Alps since 1500. Int J Climatol 25(14):1855–1880, 30 November 2005 CrossRef
  • Lowe J, Smith F, Jenkins G, Pope V (2004) Uncertainty, risk and dangerous climate change. Hadley Centre [online http://cedadocs.badc.rl.ac.uk/247/
  • Luterbacher J, Dietrich D, Xoplaki E, Grosjean M, Wanner H (2004) European seasonal and annual temperature variability, trends, and extremes since 1500. Science 303:1499–1503 CrossRef
  • Stott PA, Stone DA, Allen MR (2004) Human contribution to the European heatwave of 2003. Nature 432:610–614 (2 December 2004) CrossRef

Comparer des notes entre pays

A Montréal, l’an passé, le principal souvenir que je garderai de mon cours est le mal que j’ai eu à comprendre comment mettre des notes. Je ne ferais pas un billet sur la difficulté à évaluer le niveau d’un étudiant (car j’ai toujours détesté les examens, aussi bien en tant qu’étudiant que prof), mais j’ai une une réelle difficulté à comprendre ce que vaut une note… Par exemple, un C, ça correspond à quoi (avec mon référentiel de prof français) ? C’est bien ou pas ? C’est quoi la moyenne ?
La seule information dont on dispose pour mettre les notes est le tableau ci-contre. Et autant je pense voir ce qu’est un travail “excellent“, ou “faible“, mais “passable” ? Car le tableau me dit qu’un élève qui a 50 sur 100 a une note passable, ce qui est aussi la référence que l’on utilise en France. Sauf que les élèves ont commence à me dire que 50, c’est pas terrible… qu’une note moyenne ici c’était plutôt 60, voire 65. Bref, sous-entendu (voire entendu, tout court), “monsieur, vous devez remonter les notes“. Et comme je l’ai prouve ici, il ne faut jamais croire les élèves…
Et effectivement, dans un vieux rapport datant des années 90 (en maths et informatiques), je suis tombé sur le graphique suivant
qui donne la distribution “typique” des notes. Autrement dit, la médiane est à 63,ce qui signifie que plus de moitie a une note qui est “bien“,”très bien” voire “excellent“. En revanche, une distribution comme celle dessous est considérée comme “dure et discriminatoire
Pour mieux comprendre les notations au Canada et en France, j’ai voulu faire une petite étude (non, je ne reviendrai pas sur les notes). On continue aujourd’hui àtrouver ce genre de distribution à Montréal. Par exemple, si on considère les cours, analyse (en première année), algèbre linéaire et calcul, on a les histogrammes suivants

ce qui donne le boxplot suivant

i.e. pour les cours algèbre linéaire et calcul on est effectivement sur des moyennes autour de 60. Prenons maintenant, a titre de comparaison, des cours de maths, dans une université française, en l’occurrence a Rennes, comme analyse (en première année), arithmétique et algèbre linéaire, on a des histogrammes beaucoup plus plats,

A titre de comparaison, j’ai été voir du cote du département de Sciences Économiques, avec macroéconomie, microéconomie et mathématiques appliquées, on a

i.e. on retrouve de belles gaussiennes pour les cours d’économie, mais une distribution assez différente pour le cours de maths,

On a effectivement des distributions assez différentes, généralement, en première année à l’université, entre le Canada et la France. Ce qui explique que les profs, lorsque les étudiants souhaitent passer d’un pays a l’autre ne regarde plus vraiment les notes, mais plus le rang de l’étudiant au sein de sa promotion. D’ailleurs, si on fait un graphique quantile-quantile, avec les cours de maths à Montréal versus les cours de maths a Rennes, on a le graphique suivant

autrement dit, une note de 50 à Montréal correspond à 33 en France (soit 6.5 sur 20), alors que 50 en France (10 sur 20) correspond a 59 à Montréal. On retrouve la encore le fait qu’une note médiane correspond à 60 ici… c’est a dire un C. Il va falloir que je m’adapte l’année prochaine….

Comments

De la créativité des gangsters

Pendant mon séjour récent en Nouvelle Angleterre, j’ai survolé le livre de Leonard Mlodinow, the drunkard’s walk. Et au hasard de mes lectures, je suis tombé sur la petite histoire suivante

Autrement dit, en utilisant les cinq derniers chiffres d’une quantité économique comme le us treasury balance, on aurait un générateur de nombre aléatoire…
Par contre la suite est un peu plus surprenante,


http://freakonometrics.hypotheses.org/files/2015/12/sammy.gifLa loi de Benford (que j’avais pu évoquer ici ou ) parle des premiers chiffres, mais cette fois on parle deslast five digits. Donc visiblement l’évoquation n’est pas pertinente ici. Mais qui sait ? Ca reste malgré tout une histoire intéressante. Considérons – histoire de tester cette légende – les deux sources de données suivantes,http://treasurydirect.gov/ et http://economagic.com/. Sur ce dernier site, un petit travail de mise en forme des données est nécessaire.

b1=read.table("http://freakonometrics.free.fr/debtus1.txt",
   header=TRUE,sep="\t")
b2=read.table("http://freakonometrics.free.fr/debtus2.txt",
   header=TRUE)
X1=as.character(b1$Dollar.Amount)
n1=nchar(X1)
Y1=substr(X1,n1-8,n1-3)
X1=as.numeric(substr(Y1,1,2))*1000+as.numeric(substr(Y1,4,6))
x=X1/100000
X2=b2$DEBT
Y2=trunc(as.numeric(X2))
X2=as.character(Y2)
n2=nchar(X2)
Y2=substr(X2,n2-4,n2)
y=as.numeric(Y2)/100000
y=y[y<1]

Pour rappel, un générateur aléatoire (standard) vérifie deux propriétés importantes

  • les nombres doivent être tirés suivant une loi uniforme sur [0,1], i.e. ici, si on divise les nombres à 5 chiffres par 10000,
  • les tirages doivent être indépendants entre eux.

La première propriété semble assez naturelle, et correspond à l’histoire racontée dans un commentaire posté ici (expliquant comment un casino avait été au bord de la faillite car une roulette faisait sortir certains chiffres trop souvent, et j’essayais de comprendre comment utiliser l’information qu’un chiffre sort plus souvent). La seconde est probablement encore plus importante.

  • Visualisation de la distribution

La première idée est de visualiser la densité de nos séries de chiffres. Pour éviter les problèmes de bord (et comme c’est juste en introduction) on va utiliser un histogramme, et pas une estimation à noyau.

hist(x,col="red")
hist(y,col="blue")

On obtient pour la première série la courbe rouge, et pour la seconde la courbebleue,

On note qu’a priori, pour la première série, l’hypothèse d’uniformité n’est peut être pas la plus réaliste…

  • Test de Kolmogorov-Smirnov

On peut aussi mettre en œuvre un test de Kolmogorov-Smirnov afin de tester si la loi uniforme est adaptée.:

> ks.test(x,"punif")
 
	One-sample Kolmogorov-Smirnov test
 
data:  x 
D = 0.1047, p-value = 0.01645
alternative hypothesis: two-sided 
 
> ks.test(y,"punif")
 
	One-sample Kolmogorov-Smirnov test
 
data:  y 
D = 0.0456, p-value = 0.3581
alternative hypothesis: two-sided

On retrouve ici confirmée l’intuition précédante: la loi uniforme est pertinente pour la seconde série, pas la première.

  • Les autocorrélations de la série

Travaillons uniquement sur la seconde série. On peut étudier l’autocorrélation de notre série de nombres, ou peut-être un peu plus malin, sur les quantiles gaussiens associés (les autocorrélations étant intéressantes pour les séries gaussiennes),:

plot(acf(y))
plot(acf(qnorm(y)))

ie. pour la série brute

et pour la série normalisée,

Bref, on pourrait être tenté de valider l’hypothèse d’indépendance entre les tirages.
  • Run test (de Wald–Wolfowitz)

L’idée est de comparer une série de chiffres à la médiane, s’ils sont plus grands, on note + (ou A) et sinon – (ou B). On crée alors une série du genre “+++−−++−−++++++−−−” et on compte les séries de + et les séries de -, les runs,

library(lawstat)
runs.test(y,plot=TRUE)

Runs Test - Two sided

data:  y
Standardized Runs Statistic = -0.2462, p-value = 0.8055

Bref, la légende me semble à prendre avec des pincettes (car fonction de la source considérée), même si l’idée est intéressante (si l’on met de côté les aspects d’aléa moral). Et l’analyse sur la loi de Benford ne semble pas valide: les derniers chiffres sur les grands nombres ne se comportent pas du tout comme les premiers.

Loi des grands nombres et théorème central limite en assurance

Dans le dernier numéro de la revue Risques figure un petit article sur la loi des grands nombres et le théorème central limite en assurance (ici) Le but était de remettre un peu les pendules à l’heure sur la différence entre les implications de ces deux résultats, et surtout d’amorcer une discussion sur l’importance des hypothèses (en particulier celle de risques indépendants).

Valeur de l’euro de rente

Il y a quelques années, j’avais été sollicité pour rédiger les Annexes de l’ouvrageévaluation du préjudice corporel, de Max le Roy. Comme l’indique le descriptif du livre, j’ai un peu la pression car pas mal de jugements rendus sont basés sur mes calculs, “La résolution 75-7 du comité des ministres du Conseil de l’Europe est intervenue pour faciliter l’harmonisation des législations et des jurisprudences en matière de réparation des dommages en cas de lésion corporelle ou de décès. Elle a proposé un certain nombre de principes qui, bien que non obligatoires pour les États membres, sont très largement suivis par les juridictions françaises. Ce sont ces principes et leur application en droit français qui sont exposés et commentés dans cet ouvrage qui rend compte, notamment, de la très importante réforme apportée en matière de recours subrogatoire des caisses de sécurité sociale contre les tiers par la loi n° 2006-1640 du 21 décembre 2006. Le barème de capitalisation de rentes tient compte des plus récentes tables de mortalité et prend pour base un taux d’intérêt de 5 % conforme aux données économiques et financières actuelles. Accompagné d’un barème fonctionnel indicatif des incapacités en droit commun, de tableaux de jurisprudence et de formules de jugement, l’ouvrage fait le point sur toutes les questions auxquelles sont confrontés les magistrats, avocats, experts, médecins, assureurs.”

Avec la mise à jour des tables de mortalité en France, j’ai du remonter mes manches et reprogrammer tous les calculs. Sauf que cette fois, j’ai proposé un petit code R pour générer à nouveau ces tables. Le but est de calculer le prix de mille euros de rente viagère ou temporaire limitée. Une petite option permet de dire si on veut la nouvelle table, ou l’ancienne (ce qui m’a permis de vérifier que je ne disais pas des bêtises).

LF=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",
header=TRUE,sep=";")
LH=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",
header=TRUE,sep=";") 
PRIX=function(annee=2011,age,sex="HOM",taux=0.04,
              duree,C=1000,old=FALSE){
if(old==FALSE){
file  = paste("http://freakonometrics.blog.free.fr/public/data/",
              sex,"-table-SPLx.csv",sep="")
B     = read.table(file,header=TRUE,sep=",")
an    = annee-age; if(an>2005){an=2005}
nom   = paste("X",an,sep="")
L     = B[,nom]}
if(old==TRUE){
if(sex=="HOM"){L = c(LH$Lx,rep(0,20))}
if(sex=="FEM"){L = c(LF$Lx,rep(0,20))}}
Q     = L[(age+1):length(L)]/L[(age+1)]
actualisation = (1+taux)^(0:min(duree,120-age))
prixsup = sum(Q[2:(min(duree,120-age)+1)]/
          actualisation[2:(min(duree,120-age)+1)] )
prixinf = sum(Q[1:(min(duree,120-age))]/
          actualisation[1:(min(duree,120-age))] )
return(C*c(prixsup,prixinf))}

Les tables TGH-05 et TGF-05 sont issues d’une population de rentiers et on peut consulter les tables ici (pour les 38 pages du Journal Officiel). Ces tables ont été construites sur une population différente des tables TH et TF (construites sur l’ensemble de la population française). La méthodologie est décrite ici (merci d’ailleurs à Fred qui a approuvé mon utilisation de ces tables).

Par exemple, on peut retrouver la valeur figurant dans les anciennes Annexes (ici),

> PRIX(age=50,sex="HOM",duree=150,taux=.05,old=TRUE)
[1] 13452.15 14452.15

En bouclant, on peut générer les tables données en Annexe.

AGES=c(150,65,60,55,25,20,18,16)
tableH=matrix(0,106,length(AGES))
for(i in 1:length(AGES)){
for(a in 0:min(c(AGES[i]-1,105))){
tableH[a+1,i]=PRIX(age=a,sex="FEM",duree=
        AGES[i]-a,taux=.04,old=FALSE)[1]
}}
write(tableH,"table-femme-nouvelle-04.csv",sep=";")

(on utilise la première valeur var on veut des rentes versées en fin d’année, en cas de vie). Avec un taux de 4%, on a pour les hommes la table ici, et pour les femmes la table . Si quelqu’un voit une coquille, merci de me signaler avant que le livre ne parte à l’éditeur…

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

http://freakonometrics.hypotheses.org/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 http://freakonometrics.hypotheses.org/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 http://freakonometrics.hypotheses.org/files/2015/12/nucle-3.gif réacteurs sur une durée de http://freakonometrics.hypotheses.org/files/2015/12/nucle-04.gif années est

http://freakonometrics.hypotheses.org/files/2015/12/nucle-01.gif

Comme http://freakonometrics.hypotheses.org/files/2015/12/nucle-02.gif est supposé petit, un développement limité donnerait

http://freakonometrics.hypotheses.org/files/2015/12/nucle-06.gif

Malheureusement, dans les développements limités de ce genre, s’il faut que http://freakonometrics.hypotheses.org/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

http://freakonometrics.hypotheses.org/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….

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