Oui, je suis fâché avec les chiffres

28/05/2011 Arthur Charpentier 7 Comments

On a tous des secrets plus ou moins difficiles à révéler… Enfin, quand je dis secret c’est un bien grand mot (c’est un peu au sens que pourrait utiliser ma fille quand elle me dit “papa j’ai un secret: ce midi, j’ai mangé deux fois du dessert“). Par exemple, quand un collègue me dit “tu pourrais m’envoyer les transparents de ta présentation“, et que je dois lui répondre “tu sais, ils sont sur mon blog“: je dévoile un faux secret… Oui, tous les profs qui tiennent un blog en parlent un jour ou l’autre: avouer à des collègues qu’on tient un blog, c’est un secret que l’on révèle encore péniblement… Péniblement n’est peut être pas le mot car si une fois sur deux cela provoque une espèce de réaction condescendante, une fois sur deux cela permet au contraire d’entamer des discussions intéressantes car les interlocuteurs peuvent être curieux…

Non, en l’occurrence, je pensais plutôt aux conversations que j’ai au moins une fois par semaine, quand il faut récupérer un enfant chez un(e) ami(e) à lui (ou à elle). On discute alors avec les parents, et invariablement arrive la question “et vous faites quoi dans la vie ?”. J’aimerais bien pouvoir dire “je suis astronaute” ou “je suis constructeur de robots” ou “je suis dompteur de lions”, enfin un truc qui fait rêver. Mais non, je dois me contenter d’un “je suis prof de maths“. “Pardon?”. “Je suis prof de maths“. “Ahhhhhhh“. Oui, il y a toujours une espèce de déception que je peux lire dans le regard qui commence soudain à se dérober, comme gêné (genre il aurait espéré “je suis chirurgien cardiaque, et j’ai encore sauvé deux enfants de la mort aujourd’hui“, qui aurait donné lieu à un autre type de “ahhhhhhh”). Et là, arrive inéluctablement une phrase qui me place dans l’embarrât: au choix “moi j’ai jamais rien compris aux chiffres” ou “ben vous devez aimer les chiffres alors“…

Et là j’avoue que je ne sais jamais quoi dire… Je pourrais le prendre à la rigolade genre “oui, mais pas tous. Genre 37 je ne l’aime pas du tout“, mais le fait est que je n’aime pas le chiffres ! C’est d’ailleurs troublant que mathématique soit associé dans l’imaginaire populaire à calculs ou arithmétique. Au XVIIIème siècle, les gens qui faisaient des mathématiques étaient des géomètres. Et, en fait je préférerais qu’on me dise “ben vous devez aimer les dessins alors”… Ou à la rigueur des lettres, car si on regarde deux minutes le blog, mes papiers ou mes notes de cours, on voit des nombres d’accidents qui s’appellent , des coûts qui s’appellent , des tailles des gens qui s’appellent , des choses inconnues qui s’appellent ou suivant le problème… bref, j’aime les lettres mais surtout pas les chiffres (sauf peut être quelques uns, je ne tiens pas à me brouiller avec tout le monde) !

Je me faisais cette réflexion ce matin, en lisant l’histoire d’une autre monsieur qui était accusé d’être fâché avec les chiffres, ici

Car il était annoncé que « les deux tiers des échecs scolaires, c’est l’échec d’enfants d’immigrés ». Loin de moi l’idée de commenter l’utilisation de ces chiffres, surtout que je pense que ce monsieur et moi n’avons en commun que le fait d’être tous deux fâchés avec les chiffres. A la lecture de l’article, nous explique ensuite que “la France compte 13% d’élèves de familles immigrées, soit 87% d’élèves de familles autochtones” selon des sources OCDE. On nous montre alors un graphique expliquant que “18% (des 87%) des enfants d’autochtones et 38% des enfants d’immigrés sont sous le niveau 2 qui correspond, selon le ministre, à une situation d’échec scolaire.”

Vient alors une explication “la France compte donc 15,7% (18% de 87%) d’enfants de 15 ans issus de familles autochtones sous le niveau 2 et 4,9% d’enfants de 15 ans issus de familles immigrés (38% de 13%) sous le niveau 2. Donc sur les 20,6% de « sous-niveau 2 » nous n’avons pas deux tiers mais plutôt un quart d’enfants d’immigrés.”

Je sais pas vous, mais c’était au réveil, j’avais encore mon café dans la main, et je relisais cette phrase sans rien comprendre… et j’avais beau relire, impossible de comprendre qui disait n’importe quoi… Car personnellement, tant qu’on n’utilise pas un peu de formalisme, je suis perdu…
On peut par exemple faire un dessin… j’ai toute ma population schématisée par le grand rectangle. Dedans j’ai mes “enfant d’immigré” en rouge et “d’autochtones” en vert, et les enfants en échec en bleu, et ceux qui ne le sont pas en violet. Et la question c’est: au sein de la population bleue, quelle est la proportion qui sontbleu et rouge (qui est en jaune en bas à droite) ?

Si on note le fait d’être en échec scolaire (population bleue) et le fait d’être “enfant d’immigré” (population rouge, je ne discuterais pas la définition, mais les calculs), et que l’on suppose que ce sont les enfants “d’autochtones” (population verte), alors on cherche , et le monsieur prétend que ça doit être 2/3.
Bon, on sait que et , et que et .
Lançons nous dans les calculs: d’après la formule de Bayes,

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

L’astuce est de noter que

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

d’où la formule

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

et ici tout est connu… on a alors

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

Aussi on retrouve les 24% évoqués dans l’article…
Ce qui est bien quand on formalise, c’est que l’on est certain des chiffres qu’on avance… mais aussi, cela permettrait d’aller plus loin (conditionnellement aux chiffres à notre disposition…). Par exemple, on peut faire un tableau de contingence. On connaît les probabilités marginales

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

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

(on a fait le calcul au dessus pour écrire le dénominateur)
contenant les probabilités jointes,

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

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

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

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

Bref, on pourrait alors faire des tests du chi-deux si on savait sur combien d’élèves les tests ont été fait pour avoir de tels chiffres…

> M=matrix(c(.38*.13,.18*.87,
+ (1-.38)*.13,(1-.18)*.87),2,2)
> n=1000
> f=function(n){chisq.test(n*M)$p.value}
> F=Vectorize(f)
> P=seq(10,1000,by=10)
> plot(P,F(P),type="l",ylim=c(0,.4))
> abline(h=.05,lty=2)

i.e. on a graphiquement en fonction de la taille de l’échantillon (en abscisse) la p-value (en ordonnée) du test d’indépendance du chi-deux.

Autrement dit, sauf si ces statistiques sont tirés d’un échantillon de 200 élèves, on rejette l’hypothèse d’indépendance entre les résultats scolaires et l’origine des élèves. Mais j’avais dit que je ne commenterais pas…

Bref, mon point est qu’il me semble qu’en formalisant un peu, au lieu de s’embrouiller dans un discours rempli de chiffres, on fait moins d’erreurs.

Actuarial science, Bayesian, Credibility, Statistics

Il y a des jours avec, et des jours sans

27/05/2011 Arthur Charpentier 1 Comment

Je ne sais pas si ça vous est déjà arrivé, mais quand on fait des crêpes, la probabilité d’en réussir une est proportionnelle à la probabilité que l’on vient de réussir les dernières. Autrement dit si j’ai réussi les premières, je suis dans un bon jour, et je réussi davantage celles qui suivent. Alors que si je commence par rater les premières, je suis dans un mauvais jour, et je les rate toutes (pour les sportifs, il est aussi possible de noter que, par exemple, au tennis, la probabilité de réussir son service est d’autant plus forte que j’ai réussi les derniers: c’est pareil).

Une première formalisation du problème

Supposons que sur les premiers essais (qui servent d’étalonnage), fois j’ai réussi ma crêpe, et fois, je l’ai ratée (trop de trous, trop irrégulière, trop cuite, que sais-je encore). Mon a priori sur la probabilité de ne pas avoir de chance (et de rater ma crêpe) est . Supposons que la probabilité de réussir à l’instant (on commence à compter après les premiers essais) soit la fréquence empirique observée sur le passé. Autrement dit

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

où est le nombre de succès. En particulier . On retrouve ici le fait que plus j’en ai raté, plus j’ai de chances de rater les prochaines… c’est la mauvaise main en quelque sorte.
Afin de mieux comprendre, supposons . Autrement dit, sur mes deux premiers essais (qui ne servent qu’à étalonner), j’ai raté et réussi une fois. Juste avant l’instant 1, j’ai deux possibilité: réussir ou rater, avec la probabilité . I.e.

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

Si on se place une crêpe plus loin, et que l’on regarde les probabilités conditionnelles

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

avec quelque-chose de similaire si on conditionne par l’autre événement,

http://freakonometrics.blog.free.fr/public/perso3/crepes-14.gif

Et donc la probabilité non conditionnelle s’écrit

http://freakonometrics.blog.free.fr/public/perso3/crepes-15.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso3/crepes-16.gif

http://freakonometrics.blog.free.fr/public/perso3/crepes-18.gif

http://freakonometrics.blog.free.fr/public/perso3/crepes-19.gif

Moralité

http://freakonometrics.blog.free.fr/public/perso3/crepes-20.gif

Autrement dit, on poursuivant la récurrence, on a que est uniformément distribué sur
Étonnant, non ? Autrement dit, si j’ai commencé par un échec et une réussite, la probabilité d’avoir 50 réussites sur 100 essais sera très exactement . Qui est la même probabilité que d’avoir 1 réussite, ou 99 !
On a aussi que

http://freakonometrics.blog.free.fr/public/perso3/crepes-22.gif

Mais que se passe-t-il si et ne sont plus égaux à 1 ?

Lien avec l’approche bayésienne et la théorie de la crédibilité

Ce qui est amusant, c’est que la probabilité

est bien connue par les personnes qui connaissent l’assurance et la théorie de la crédibilité. Supposons que je note la variable indicatrice correspondant à une réussite ou un succès au ème essai (avec probabilité ). D’un point de vu bayésien, la loi conjuguée de la loi binomiale (on va compter les réussites sur n essais) est une loi Beta. Rappelons que la loi Beta a pour densité

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

sur , et que l’espérance de cette variable est , i.e. . Donc cette modélisation devrait parfaitement convenir à notre problème. Ce qui est intéressant, c’est que comme on travaille avec la loi conjuguée, la loi ex post est encore une loi Beta. Autrement dit, ici, si je suppose que suit a priori une loi Beta , et que la loi de sachant suit une loi binomiale , alors la loi a posteriori de sachant est une loi (c’est ce que l’on avait vu ici). Et en particulier

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

Autrement dit on obtient une combinaison convexe entre la moyenne de la loi a priori (à gauche) et la fréquence empirique (à droite).
On retrouve ici ce que nous avions noté auparavant, en notant que , et .
On peut toutefois se demander si le résultat que l’on avait obtenu sur la loi de S_n reste valide ici…
Notons que

http://freakonometrics.blog.free.fr/public/perso3/crepes-38.gif

i.e. (moyennant quelques calculs)

http://freakonometrics.blog.free.fr/public/perso3/crepes-40.gif

Dans le cas où , en réécrivant aussi le premier terme, on a

http://freakonometrics.blog.free.fr/public/perso3/crepes-41.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso3/crepes-42.gif

On retrouve bien la formule obtenue auparavant… Ce qui est intéressant aussi, c’est que si on note que

http://freakonometrics.blog.free.fr/public/perso3/crepes-43.gif

http://freakonometrics.blog.free.fr/public/perso3/crepes-44.gif

et donc

http://freakonometrics.blog.free.fr/public/perso3/crepes-45.gif

et par récurrence, on obtient que

http://freakonometrics.blog.free.fr/public/perso3/crepes-46.gif

et que de plus

http://freakonometrics.blog.free.fr/public/perso3/crepes-47.gif

Autrement dit, effectivement, si j’en ai raté beaucoup au début, ma probabilité d’en rater beaucoup est grande (et elle est d’ailleurs constante). Certains parleraient de no-free lunch.

Computer, Maths

Découpage équitable de la pizza avec des olives

25/05/2011 Arthur Charpentier 3 Comments

“mais non papa, tu comprends jamais rien ! ok, c’est bien si on a autant de pizza tous les deux… mais on veut aussi avoir autant d’olives tous les deux ! et autant de morceaux d’ananas !”. Ouais, bon, on n’est pas couché là !

On va quand même essayer de la couper cette pizza (avant qu’elle ne refroidisse) !

Partage de la pizza et des olives

Commençons par supposer que seules les olives soient importantes. La première idée que j’aurais serait de faire une régression médiane, qui va découper le nuage de points (les olives) en deux…

Sous R (oui, j’ai installé R sur la table de la cuisine), c’est assez facile à mettre en œuvre… Pour des raisons de symétrie, je vais faire les régressions de l’abscisse sur l’ordonnée, et réciproquement…

> library(quantreg)
> lq=rq(y ~ x, tau = 0.5, data=D)
> abline(lq,lwd=5,col="red")
> sum(y0>predict(lq))
[1] 22
> sum(y0<=predict(lq))
[1] 22
>
> lq=rq(x ~ y, tau = 0.5, data=D)
> abline(b=1/lq$coefficients[2],
+ a=-lq$coefficients[1]/lq$coefficients[2],lwd=5,col="red")
> sum(x0>predict(lq))
[1] 22
> sum(x0<=predict(lq))
[1] 22

Le soucis est que si la régression n’est pas contrainte, il y a une répartition 50-50 des olives, mais les parts ne sont pas équitables,

> library(quantreg)
> lq=rq(y ~ 0+x, tau = 0.5, data=D)
> abline(lq,lwd=5,col="red")
> sum(y0>predict(lq))
[1] 20
> sum(y0<=predict(lq))
[1] 24
>
> lq=rq(x ~ 0+y, tau = 0.5, data=D)
> abline(b=1/lq$coefficients,a=0,lwd=5,col="red")
> sum(x0>predict(lq))
[1] 25
> sum(x0<=predict(lq))
[1] 19

et si je supprime la constante les parts sont équitables (la régression passe par le centre la pizza), mais le partage des olives n’est pas équilibré…

On peut aussi calculer pour tous les angles possibles ce que donnerait des partages de la pizza (qui passent par le centre). C’est ce que représente la courbe rouge ci-dessous: l’extérieur du disque le ratio des olives dépasse 1, à l’intérieur il est inférieur à 1. Les découpages optimaux correspondent au cas où le ratio vaut précisément 1… ce qui est possible si on coupe comme indiqué ci-dessous, à droite,

Je pense qu’on doit pouvoir montrer qu’il est toujours possible de couper la pizza en deux, de manière équitable, avec autant d’olives de part et d’autre (s’il y a un nombre pair d’olives).
“oui mais papa… tu nous avais dit que tu savais pas couper en passant par le centre…”. Ah oui, c’est vrai ! Damned…. bon, utilisons la méthode décrite dans le billet précédant (ici), a priori ça doit être possible de trouver un sommet pour lequel le partage des olives est équitable. En fait, on pourrait aussi tirer au hasard des points dans la pizza, et compter les olives en faisant la répartition de manière alternée,

> sommet=matrix(NA,500000,3)
> for(s in 1:nrow(sommet)){
+ a0=runif(10)*2-1; b0=runif(10)*2-1
+ if(a0^2+b0^2)<=1){
+ slope1=(y0-b0)/(x0-a0)
+ I1=((slope1>0)&(slope1<1))|((slope1<(-1)))
+ sommet[s,]=c(sum(I1),a0,b0)
}}

On a alors en rouge l’ensemble des points à partir desquels on peut faire des découpages en 8, avec des pointes à 45° (en commençant par couper horizontalement)

Et ça semble marcher !

“oui mais papa… tu as oublié les ananas…”.

Partage de la pizza et des olives et des ananas

En utilisant la méthode précédente, rien de plus simple ! En plus de compter la répartition des olives, on va compter la répartition des ananas…! Sur le dessin à gauche on a les olives et les ananas. A droite, les points rouges correspondent à des répartitions équitables des olives, et les bleus des répartitions équitables des ananas… Reste à trouver des points à la fois rouges et bleus,

ce sont finalement les points violets ci-dessous ! On a alors une méthode pour couper la pizza en deux parts de même taille, avec autant d’olives et autant d’ananas de part et d’autres !

Maths

Découpage équitable de la pizza

25/05/2011 Arthur Charpentier 5 Comments

Bon, maintenant que l’année scolaire est passée, je vais arrêter un peu les billets sérieux pour revenir pour mes mésaventures à la maison. Par exemple, hier soir, c’était pizza. Et mes enfants aiment partager à deux une pizza (sur laquelle on a rajouté plein de choses qu’ils aiment, genre des olives, ou des morceaux d’ananas). Mais il faut un découpage équitable. Les enfants ne veulent pas un partage juste, mais un partage moitié-moitié.

Classiquement, je coupe la pizza en deux (cf figure ci-contre, avec un côté en bleu pour mon fils, et un en rouge pour ma fille, histoire de simplifier un peu). Mais comme me l’a fait noté mon fils, si mon couteau ne passe pas par le centre de la pizza (on va dire le centre du cercle parce qu’il va falloir faire des maths… pour les recettes de cuisine, il faut changer de site, genre ici).
Donc pour faire simple, la proba d’avoir deux parts rigoureusement égale est très très faible… car il faut couper suivant une droite qui passe exactement au centre (sinon, il y a quelqu’un qui est défavorisé, alors on tire au sort… a prioric’est juste juste, mais ex-post, il y en a toujours un ou une qui trouve le système injuste et qui veut qu’on retire au hasard…)
Fort heureusement, il existe une technique simple, et rigoureusement juste, permettant d’avoir des parts équitables…
L’idée est la suivante: on commence par couper horizontalement, n’importe où. Ensuite, on coupe verticalement (à angle droit, en prenant son équerre) de manière à récupérer 4 parts. Au point d’intersection, on coupe à 45° direction nord-est, et 45° direction nord-ouest…. Une fois le découpage effectué, on réparti les parts, à tour de rôle, les parts bleues pour mon fils, et celles en rouges pour ma fille,

Car le truc génial est que le découpage – en plus d’être fort élégant – est parfaitement équitable ! Comme les enfants sont toujours sceptiques, il faut bien sûr leur prouver que oui, ils en auront exactement pareil. Pour ça on reprend le couteau…

Si on recoupe les parts suivant les traits-ci-dessus, chaque part bleue pourra être superposée à une part rouge (la preuve mathématique repose sur des simple notions de symétrie…). Génial non ? On a même le théorème du mardi soir, à savoir que “si on coupe une pizza en 8 parts, avec des pointes à 45° partir d’un point choisi arbitrairement dans la pizza, alors la somme alternée des parts sera égale” (cf ici). En 1967, Upton a posé ce problème,

Larry Carter et Stan Wagon affirmaient que ce résultat resterait vrai si on coupait en 12 ou 16 parts (avec des angles égaux pour les pointes), mais pas pour 6, 10, 14 ou 18 par exemple. Le blog du chou romanesco (ici et là) prétend que c’est vrai pour les découpages en un nombres pairs. En fait, Rick Mabry et Paul Deiermann ont un théorème assez général (là pour la version complète) et qui généralise aussi en dimension 3 (pour une calzone en quelque sorte). Mais je sens que ça n’intéresse plus trop les enfants, qui ont déjà fini la pizza….

Les plus curieux peuvent aussi consulter le court article ci-dessous,

Statistics

Comment parle un président français le 31 décembre ?

13/05/2011 Arthur Charpentier 1 Comment

Il y a pas mal d’agitation présidentielle sur les blogs (français) cette semaine, que ce soit au sujet des 30 ans de l’élection de François Mitterrand, ou l’annonce de la date de la prochaine élection présidentielle (et probablement d’autres raisons). Un peu naïvement je me suis demandé si les présidents étaient vraiment si différents que ça ? Ou plus précisément, est-ce qu’une fois investis par la fonction, tous ces hommes ne se comporteraient-ils pas – finalement – tous pareils ?

Pour cela je me suis simplement demandé si on pouvait distinguer les discours des différents Présidents. J’ai donc pris quelques discours entendus les 31 décembre, lors de vœux retransmis à la télévision (en ligne ici, en excluant les discours aux forces armées, qui me semblaient biaiser un peu mon échantillon). On commence par récupérer la liste de tous les mots prononcés (les fichiers csv sont en ligne ici)

setwd("/home/arthur/Documents/r-data/voeux/")
FILES=list.files()
library(tau)
BASE=scan(FILES[1],what="character")
WORD=fixEncoding(BASE,latin1=TRUE)
 for(i in 2:length(FILES)){
BASE=scan(FILES[i],what="character")
 WORD=c(WORD,fixEncoding(BASE,latin1=TRUE))
 }

Ensuite on regarde le top 500 (j’exclue les mots communs genre le ou la, mais je garde les pronoms, comme mon ou notre),

> WORD=toupper(WORD)
> LW=rev(sort(table(WORD)))
> LISTE=names(LW[1:500])
> LISTE=LISTE[-c(1,2,3,4,5,6,7,8,9,10,11,14,15,18,19,20,21,22,23,25,
               27,43,40,36,41,60)]
> LISTE[1:25]
 [1] "NOUS"   "JE"     "VOUS"   "NOTRE"  "PLUS"   "FRANCE" "PAS"    "NOS"   
 [9] "EST"    "AVEC"   "SE"     "MES"    "CEUX"   "MAIS"   "SONT"   "OU"    
[17] "CETTE"  "AUSSI"  "TOUS"   "SUR"    "SI"     "PAYS"   "LEUR"   "FAIRE" 
[25] "VOTRE"

et créer une base de décompte, par discours,

ANNEE=as.numeric(substr(FILES,6,7))
AN=1900+ANNEE*(ANNEE>15)+(100+ANNEE)*(ANNEE<15)
PRES=c("CHIRAC","CHIRAC","CHIRAC","SARKOZY","SARKOZY",
"SARKOZY","VGE","VGE","VGE","MITTERRAND","MITTERRAND",
"MITTERRAND","MITTERRAND","MITTERRAND")
NOM=paste(PRES,AN,sep="-")
COMPTE=matrix(NA,length(FILES),length(LISTE))
for(i in 1:length(FILES)){
BASE=scan(FILES[i],what="character")
WORD=toupper(fixEncoding(BASE,latin1=TRUE))
for(j in 1:length(LISTE)){
COMPTE[i,j]=sum(WORD==LISTE[j])
}}
colnames(COMPTE)=LISTE
rownames(COMPTE)=NOM

On peut alors faire tourner des procédures d’analyse de données, comme une ACP

library("FactoMineR")
res.pca = PCA(COMPTE, scale.unit=TRUE, ncp=5, graph=FALSE)
plot.PCA(res.pca, axes=c(1, 2), choix="ind")

Les projections des individus (les discours) sont les suivantes

alors que celles pour les variables (les mots) sont, toujours sur les deux premiers axes factoriels,

On peut jouer des heures à se demander qui dit quoi, mais ce qui m’intéressait était de voir si les présidents étaient proches, en terme de discours. On peut alors faire des classifications hiérarchiques,

DISTANCE = dist(COMPTE)
h = hclust(DISTANCE) 
plot(as.dendrogram(h), horiz = TRUE)

ou ici avec un autre algorithme de classification

On notera les regroupements dépendent peu des distances et des méthodes de constitution de classes… et que les regroupements se font par président: tous les discours de François Mitterrand semblent proches, et tous ceux de Jacques Chirac aussi (en tous les cas ceux que j’ai téléchargé, mais l’analyse n’est pas exhaustive). On notera aussi la proximité entre Valéry Giscard d’Estaing l’année de son élection, et Nicolas Sarkozy l’année de son élection….

C’est finalement un peu rassurant de voir que les Présidents ne sont pas tous les mêmes, et qu’on peut les reconnaitre en fonction de leur discours… c’est tout du moins la prochaine étape, ou on va voir si on peut affecter des nouveaux discours (que je n’ai pas utilisés pour l’analyse) dans des classes, et repérer qui l’a écrit…

Inference, Regression

Multivariate probit regression using (direct) maximum likelihood estimators

11/05/2011 Arthur Charpentier 4 Comments

Consider a random pair of binary responses, i.e. with taking values 1 or 2. Assume that probability can be function of some covariates.

The Gaussian vector latent structure

A standard model is based a latent Gaussian structure, i.e. there exists some random vector such that if is lower than a given threshold, and 1 otherwise.
As in standard probit models, assume that

http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-09.gif

where we can assume that is a Gaussian random vector. This assumption can be used to derive the likelihood of a sample .

> logV=function(parameter){
+ CORRELATION=parameter[1]
+ BETA=matrix(parameter[2:length(parameter)],ncol(Y),ncol(X))
+ z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]))
+ sigma=matrix(c(1,CORRELATION,CORRELATION,1),2,2)
+     a11=pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
+ for(i in 2:nrow(z)){a11=c(a11,pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
+     a10=pnorm(z[1,1],sd=sqrt(sigma[1,1]))-pmnorm(z[1,],varcov=sigma)
+ for(i in
+ 2:nrow(z)){a10=c(a10,pnorm(z[i,1],sd=sqrt(sigma[1,1]))-pmnorm(z[i,],varcov=sigma))}
+     a01=pnorm(z[1,2],sd=sqrt(sigma[2,2]))-pmnorm(z[1,],varcov=sigma)
+ for(i in
+ 2:nrow(z)){a01=c(a01,pnorm(z[i,2],sd=sqrt(sigma[2,2]))-pmnorm(z[i,],varcov=sigma))}
+     a00=1-a10-a01-a11
+ -sum(((Y[,1]==1)&(Y[,2]==1))*log(a11) +
+     ¹*log(a01) +
+     ²*log(a10) +
+     ³*log(a00) )
+ }
> OPT=optim(fn=logV,par=c(0,1,1,1,1,1,1),method="BFGS")$par

(the code is a bit long since I had trouble working properly with matrices – or more precisely to vectorize my functions – so I used loops… I am sure it is possible to write a better code).
It is possible to generate samples (based on that specific model) to check that we can actually derive proper maximum likelihood estimators,

> library(mnormt)
> set.seed(1)
> n=1000
> r=0.5
> X1=runif(n)
> X2=rnorm(n)
> Y1S=1+5*X1
> Y2S=8-5*X1
> RES=rmnorm(n,mean=c(0,0),varcov=matrix(c(1,r,r,1),2,2))
> YS=cbind(Y1S,Y2S)+RES
> Y1=(YS[,1]>quantile(YS[,1],.5))*1
> Y2=(YS[,2]>quantile(YS[,2],.5))*1
> base=data.frame(i,Y1,Y2,X1,X2,YS)
> head(base)
  i Y1 Y2        X1          X2      Y1S      Y2S
1 1  0  0 0.2655087  0.07730312 3.177587 5.533884
2 2  0  0 0.3721239 -0.29686864 1.935307 5.089524
3 3  1  0 0.5728534 -1.18324224 4.757848 5.172584
4 4  1  0 0.9082078  0.01129269 4.600029 3.878225
5 5  0  1 0.2016819  0.99160104 2.547362 6.743714
6 6  1  0 0.8983897  1.59396745 5.309974 4.421523

(the two columns on the right are latent observations, that cannot be used since theoretically they are unobservable). Note that it is a simple regression, one of the component is here only to bring some noise. First of all, let us look at marginal probit regression

>  reg1=glm(Y1~X1+X2,data=base,family=binomial)
>  reg2=glm(Y2~X1+X2,data=base,family=binomial)
> summary(reg1)
 
Call:
glm(formula = Y1 ~ X1 + X2, family = binomial, data = base)
 
Deviance Residuals:
Min        1Q    Median        3Q       Max
-2.90570  -0.50126  -0.00266   0.49162   2.78256
 
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.291725   0.267149 -16.065   <2e-16 
X1           8.656836   0.510153  16.969   <2e-16 ***
X2           0.007375   0.090530   0.081    0.935
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1386.29  on 999  degrees of freedom
Residual deviance:  726.48  on 997  degrees of freedom
AIC: 732.48

Number of Fisher Scoring iterations: 5
> summary(reg2)
Call:
glm(formula = Y2 ~ X1 + X2, family = binomial, data = base)
Deviance Residuals:
Min        1Q    Median        3Q       Max
-2.74682  -0.51814  -0.00001   0.57969   2.58565
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  3.91709    0.24399  16.054   <2e-16 ***
X1          -7.89703    0.46277 -17.065   <2e-16 ***
X2           0.18360    0.08758   2.096    0.036 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)

Null deviance: 1386.29  on 999  degrees of freedom
Residual deviance:  777.61  on 997  degrees of freedom
AIC: 783.61
Number of Fisher Scoring iterations: 5

Here, the optimization yields,

> OPT=optim(fn=logV,par=c(0,1,1,1,1,1,1),method="BFGS")$par
> OPT[1]
[1] 0.5261382
> matrix(OPT[2:7],2,3)
          [,1]      [,2]       [,3]
[1,] -2.451721  4.908633 0.01600769
[2,]  2.241962 -4.539946 0.10614807

Note that the coefficients we have obtained are almost identical to the ones obtained with R standard procedure,

>  library(Zelig)
>  REG= zelig(list(mu1=Y1~X1+X2,
+             mu2=Y2~X1+X2,
+     rho=~1),
+     model="bprobit",data=base)
>  summary(REG)
 
Call:
zelig(formula = list(mu1 = Y1 ~ X1 + X2, mu2 = Y2 ~ X1 + X2,
    rho = ~1), model = "bprobit", data = base)
 
Pearson Residuals:
                 Min        1Q     Median      3Q     Max
probit(mu1) -10.5442 -0.377243  0.0041803 0.36709 8.60398
probit(mu2)  -7.8547 -0.376888  0.0083715 0.42923 5.88264
rhobit(rho) -13.8322 -0.091502 -0.0080544 0.37218 0.85101
 
Coefficients:
                  Value Std. Error   t value
(Intercept):1 -2.451699   0.135369 -18.11116
(Intercept):2  2.241964   0.125072  17.92536
(Intercept):3  1.169461   0.189771   6.16249
X1:1           4.908617   0.252683  19.42602
X1:2          -4.539951   0.233632 -19.43203
X2:1           0.015992   0.050443   0.31703
X2:2           0.106154   0.049092   2.16235
 
Number of linear predictors:  3
 
Names of linear predictors: probit(mu1), probit(mu2), rhobit(rho)
&n
bsp;
Dispersion Parameter for binom2.rho family:   1
 
Residual Deviance: 1460.355 on 2993 degrees of freedom
 
Log-likelihood: -730.1774 on 2993 degrees of freedom
 
Number of Iterations: 3

> matrix(coefficients(REG)[c(1:2,4:7)],2,3)
          [,1]      [,2]       [,3]
[1,] -2.451699  4.908617 0.01599183
[2,]  2.241964 -4.539951 0.10615443

The correlation here is also the same

> (exp(summary(REG)@coef3[3])-1)/(exp(summary(REG)@coef3[3])+1)
[1] 0.5260951

That procedure works well an can be extended to ordinal responses (not only binary ones, or to three dimensional problems,

logV=function(beta){
BETA=matrix(beta[4:(3+ncol(Y)*ncol(X))],ncol(Y),ncol(X))
z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]),X%*%(BETA[3,]))
r12=beta[1]
r23=beta[2]
r31=beta[3]
s1=s2=s3=1
sigma=matrix(c(s1^2,r12*s1*s2,r31*s1*s3,
               r12*s1*s2,s2^2,r23*s2*s3,
               r31*s1*s3,r23*s2*s3,s3^2),3,3)
sigma1=matrix(c(s2^2,r23*s2*s3,
                r23*s2*s3,s3^2),2,2)
sigma2=matrix(c(s1^2,r31*s1*s3,
                r31*s1*s3,s3^2),2,2)
sigma3=matrix(c(s1^2,r12*s1*s2,
                r12*s1*s2,s2^2),2,2)
    a111=pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a111=c(a111,pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
    a011=pmnorm(z[1,2:3],varcov=sigma1)-pmnorm(z[1,],varcov=sigma)
for(i in 2:nrow(z)){a011=c(a011,pmnorm(z[i,2:3],varcov=sigma1)-pmnorm(z[i,],varcov=sigma))}
    a101=pmnorm(z[1,c(1,3)],varcov=sigma2)-pmnorm(z[1,],varcov=sigma)
for(i in 2:nrow(z)){a101=c(a101,pmnorm(z[i,c(1,3)],varcov=sigma2)-pmnorm(z[i,],varcov=sigma))}
    a110=pmnorm(z[1,1:2],varcov=sigma3)-pmnorm(z[1,],varcov=sigma)
for(i in 2:nrow(z)){a110=c(a110,pmnorm(z[i,1:2],varcov=sigma3)-pmnorm(z[i,],varcov=sigma))}
    a100=pnorm(z[1,1],sd=s1)-pmnorm(z[1,c(1,2)],varcov=sigma3)-pmnorm(z[1,c(1,3)],varcov=sigma2)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a100=c(a100,pnorm(z[i,1],sd=s1)-pmnorm(z[i,c(1,2)],varcov=sigma3)-pmnorm(z[i,c(1,3)],varcov=sigma2)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
    a010=pnorm(z[1,2],sd=s2)-pmnorm(z[1,c(1,2)],varcov=sigma3)-pmnorm(z[1,c(2,3)],varcov=sigma1)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a010=c(a010,pnorm(z[i,2],sd=s2)-pmnorm(z[i,c(1,2)],varcov=sigma3)-pmnorm(z[i,c(2,3)],varcov=sigma1)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
    a001=pnorm(z[1,3],sd=s3)-pmnorm(z[1,c(2,3)],varcov=sigma1)-pmnorm(z[1,c(1,3)],varcov=sigma2)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a001=c(a001,pnorm(z[i,3],sd=s3)-pmnorm(z[i,c(2,3)],varcov=sigma1)-pmnorm(z[i,c(1,3)],varcov=sigma2)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
    a000=1-a111-a011-a101-a110-a001-a010-a100
 
a111[a111<=0]=1e-50
a110[a110<=0]=1e-50
a101[a101<=0]=1e-50
a011[a011<=0]=1e-50
a100[a100<=0]=1e-50
a010[a010<=0]=1e-50
a001[a001<=0]=1e-50
a000[a000<=0]=1e-50
 
-sum(((Y[,1]==0)&(Y[,2]==0)&(Y[,3]==0))*log(a111) +
    ⁴*log(a011) +
    ⁵*log(a101) +
    ⁶*log(a110) +
    ⁷*log(a001) +
    ⁸*log(a010) +
    ⁹*log(a100) +
    ¹⁰*log(a000) )
}

A strong assumption in that bivariate model is that residuals have a Gaussian structure. It is possible to change that assumption

marginally: for instance if we use a logistic cumulative distribution function, then we will have a bivariate logit regression
in terms of dependence structure: it is possible to consider another copula than the gaussian one, e.g. Gumbel’s copula (also called the bivariate logistic copula), or Clayton’s

Here, the following code can be used to extend the model to non Gaussian structures,

> F=function(x,r){pmnorm(x,rep(0,length(x)),
+                 varcov=matrix(c(1,r,r,1),2,2))}
> Fx=function(x1){F(c(x1,1e40),0)}
> Fy=function(x2){Fx(x2)}
> 
> logVgen=function(parameter){
+ CORRELATION=parameter[1]
+ BETA=matrix(parameter[2:length(parameter)],ncol(Y),ncol(X))
+ z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]))
+     a11=F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a11=c(a11,F(z[i,],r=CORRELATION))}
+     a10=Fx(z[1,1])-F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a10=c(a10,Fx(z[i,1])-F(z[i,],r=CORRELATION))}
+     a01=Fy(z[1,2])-F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a01=c(a01,Fy(z[i,2])-F(z[i,],r=CORRELATION))}
+     a00=1-a10-a01-a11
+ -sum(((Y[,1]==1)&(Y[,2]==1))*log(a11) +
+     ¹¹*log(a01) +
+     ¹²*log(a10) +
+     ¹³*log(a00) )
+ }
>
> beta0=c(0,1,1,1,1,1,1)
> (OPT=optim(fn=logVgen,par=beta0,method="BFGS")$par)
[1]  0.52613820 -2.45172059  2.24196154  4.90863292 -4.53994592  0.01600769
[7]  0.10614807
There were 23 warnings (use warnings() to see them)

E.g.

> library(copula)
> F=function(x,r){pcopula(pnorm(x),
               claytonCopula(2, r))}
> Fx=function(x1){F(c(x1,1e40),0)
}
> Fy=function(x2){Fx(x2)}

An application to school tests

Consider the following dataset,

hsb2=read.table("http://freakonometrics.free.fr/hsb2.csv",
        header=TRUE, sep=",")
math_male=hsb2$math[female==0]
write_male=hsb2$write[female==0]
math_female=hsb2$math[female==1]
write_female=hsb2$write[female==1]
plot(math_female, write_female, type="p",
     pch=19,col="red",xlab="maths",ylab="writing",cex=.8)
points(math_male, write_male, cex=1.2, col="blue")

with here maths versus writing, with girls in red and boys in blue, where variables here are

  female :
    0: male
    1: female
  race :
    1: hispanic
    2: asian
    3: african-amer
    4: white
  ses :
    1: low
    2: middle
    3: high
  schtyp : type of school
    1: public
    2: private
  prog : type of program
    1: general
    2: academic
    3: vocation
  read : reading score
  write : writing score
  math : math score
  science : science score
  socst : social studies score

We can try to understand correlation between math and writing skills. Covariates can be the sex of the child, and his reading skills. The question will then be: are good students in maths and writing simply students that can read well ?

Here the code is simply

> W=hsb2$write>=50
> M=hsb2$math>=50
> base=data.frame(Y1=W,Y2=M,
+             X1=hsb2$female,X2=hsb2$read)
>
> library(Zelig)
> REG= zelig(list(mu1=Y1~X1+X2,
+             mu2=Y2~X1+X2,
+     rho=~1),
+     model="bprobit",data=base)
> summary(REG)
 
Call:
zelig(formula = list(mu1 = Y1 ~ X1 + X2, mu2 = Y2 ~ X1 + X2,
    rho = ~1), model = "bprobit", data = base)
 
Pearson Residuals:
                Min        1Q  Median      3Q    Max
probit(mu1) -4.7518 -0.502594 0.15038 0.53038 1.8592
probit(mu2) -3.4243 -0.653537 0.23673 0.67011 2.6072
rhobit(rho) -4.9821  0.010481 0.13500 0.40776 2.9171
 
Coefficients:
                  Value Std. Error  t value
(Intercept):1 -5.484711   0.787101 -6.96825
(Intercept):2 -4.061384   0.633781 -6.40818
(Intercept):3  1.332187   0.322175  4.13497
X1:1           1.125924   0.233550  4.82092
X1:2           0.167258   0.202498  0.82598
X2:1           0.103997   0.014662  7.09286
X2:2           0.082739   0.012026  6.88017
 
Number of linear predictors:  3
 
Names of linear predictors: probit(mu1), probit(mu2), rhobit(rho)
 
Dispersion Parameter for binom2.rho family:   1
 
Residual Deviance: 364.51 on 593 degrees of freedom
 
Log-likelihood: -182.255 on 593 degrees of freedom
 
Number of Iterations: 3
> (exp(summary(REG)@coef3[3])-1)/(exp(
summary(REG)@coef3[3])+1)
[1] 0.5824045

with a remaining correlation among residuals of 0.58. So with only the sex of the student, and his or her reading skill, we cannot explain the correlation between maths and writing skills. With our previous code, we have here

> beta0=c((exp(summary(REG)@coef3[3])-1)/(exp(summary(REG)@coef3[3])+1),
+      summary(REG)@coef3[c(1:2,4:7),1])
> beta0
              (Intercept):1 (Intercept):2          X1:1          X1:2
0.58240446   -5.48471133   -4.06138412    1.12592427    0.16725842
X2:1          X2:2
0.10399668    0.08273879
> (OPT=optim(fn=logV,par=beta0,method="BFGS")$par)
(Intercept):1 (Intercept):2          X1:1          X1:2
0.5824045    -5.4847113    -4.0613841     1.1259243     0.1672584
X2:1          X2:2
0.1039967     0.0827388

i.e. we obtain (almost) exactly the same estimators. But here I have used as starting values for the optimization procedure the estimators given by R. If we change them, hopefully we have a robust maximum likelihood estimator,

> (OPT=optim(fn=logV,par=beta0/2,method="BFGS")$par)
              (Intercept):1 (Intercept):2          X1:1          X1:2
   0.58233360   -5.49428984   -4.06839571    1.12696594    0.16760347
         X2:1          X2:2
   0.10417767    0.08287409
There were 12 warnings (use warnings() to see them)

So once again, it is possible to optimize numerically a likelihood function, and it works.

Y[,1]==0)&(Y[,2]==1 [↩]
Y[,1]==1)&(Y[,2]==0 [↩]
Y[,1]==0)&(Y[,2]==0 [↩]
Y[,1]==1)&(Y[,2]==0)&(Y[,3]==0 [↩]
Y[,1]==0)&(Y[,2]==1)&(Y[,3]==0 [↩]
Y[,1]==0)&(Y[,2]==0)&(Y[,3]==1 [↩]
Y[,1]==1)&(Y[,2]==1)&(Y[,3]==0 [↩]
Y[,1]==1)&(Y[,2]==0)&(Y[,3]==1 [↩]
Y[,1]==0)&(Y[,2]==1)&(Y[,3]==1 [↩]
Y[,1]==1)&(Y[,2]==1)&(Y[,3]==1 [↩]
Y[,1]==0)&(Y[,2]==1 [↩]
Y[,1]==1)&(Y[,2]==0 [↩]
Y[,1]==0)&(Y[,2]==0 [↩]

Computer, Statistics

Maudits français

11/05/2011 Arthur Charpentier 5 Comments

Tous les malouins le savent, Jacques Cartier est parti de Saint Malo pour découvrir le Canada en 1536 (si l’on considère son second voyage, c’est à dire celui où il a vraiment pénétré le royaume de Kanata en voguant sur le Saint Laurent). Il y a d’ailleurs un musée du Québec dans les remparts. Et comme s’en vente le musée, les québécois viennent de Saint-Malo.

Or, pour ceux qui ne connaissent pas Saint-Malo, c’est à deux pas du Mont Saint Michel, et du Couesnon qui sépare la Bretagne et la Normandie. Étant né en Normandie, avant de venir vivre à Rennes, j’étais convaincu que j’avais des racines communes avec les Québécois^*. Mais j’ai été surpris en voyant les noms de famille des copines de classe de ma fille: les noms ne m’étaient pas du tout familiers… Pourtant, si nous avons des racines communes, il ne serait pas surprenant que les noms de familles se retrouvent (c’est le principe de base de la généalogie, si j’ai bien tout suivi). Aussi j’ai voulu creuser davantage….

Que disent les experts en généalogie ?

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

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

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

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

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

On peut trouver des éléments de réponse ici ou là. Parmi les 3800 Français qui ont immigré (ou émigrer, ça dépend du point de vue) entre 1608 et 1700 :

La Normandie a fourni 1 pionnier sur 5.
Le Poitou, l’Aunis et la Saintonge ont, ensemble, donné 1 pionnier sur 4.
Paris et l’Ile-de-France, 1 pionnier sur 7.
La Bretagne, un peu plus de 3 pionniers sur 100.
L’Anjou, 3 pionniers sur 100.
La Champagne, presque 3 pionniers sur 100.
La Picardie, un peu plus de 2 pionniers sur 100.

Entre 1700 et 1765, près de 5000 autres hommes, femmes et enfants originaires de France viennent s’établir au Canada, mais “très peu proviennent de la Bretagne et des régions baignées par la mer“. Ce sont d’ailleurs les ordres de grandeur que j’ai pu avoir dans l’ouvrage de Michel Lambert (qui n’est pas un livre de généalogie à proprement parler mais qui est incroyablement intéressant, même s’il ne donne pas de source – ou de détails sur le sens – pour ces chiffres),

Normandie, 22.6%
Aunis, 16.4%
Perche, 11.4%
Ile de France (région parisienne) 10.5%
Poitou, 7.5%
Maine, 5.2%
Saintonge 5.2%
Bretagne 3.5%

Moralité, même si Jacques Cartier est parti de Bretagne, rares semblent être les bretons qui ont émigrés au Québec.

Calculs de probabilités conditionnelles à partir des noms de famille

Il existe des sites qui donnent les “classements” des noms de familles, au Québec par exemple (ici, avec les 1000 noms les plus portés, avec les proportions respectives), ou dans les régions françaises (là, où j’ai pris à chaque fois les 2000 noms les plus portés avec les proportions respectives, en changeant de région). Les données peuvent être récupérées avec le code suivant

quebec=read.table("http://freakonometrics.free.fr/nom-quebec2.txt",   
                   header=TRUE,dec=",")
bretagne=read.table("http://freakonometrics.free.fr/nom-bretagne.txt",
                    header=TRUE,dec=",",sep="\t")
poitou=read.table("http://freakonometrics.free.fr/nom-poitou.txt",
                  header=TRUE,dec=",",sep="\t")
normandie=read.table("http://freakonometrics.free.fr/nom-basse-normandie.txt",
                     header=TRUE,dec=",",sep="\t")
aquitaine=read.table("http://freakonometrics.free.fr/nom-aquitaine.txt",
                     header=TRUE,dec=",",sep="\t")
alsace=read.table("http://freakonometrics.free.fr/nom-alsace.txt",  
                  header=TRUE,dec=",",sep="\t")
loire=read.table("http://freakonometrics.free.fr/nom-loire.txt",
                 header=TRUE,dec=",",sep="\t")

J’ai retenu 6 régions françaises, et j’ai calculé la proportions de québécois dont le nom de famille était dans le top 2000 d’une région donnée, i.e.

> head(quebec)
  Rang Nomdefamille Pourcentage
1    1     Tremblay       1.076
2    2       Gagnon       0.790
3    3          Roy       0.753
4    4         Cote       0.692
5    5     Bouchard       0.530
6    6     Gauthier       0.522

> minquebec=tolower(quebec$Nomdefamille)
> rangquebec=quebec$Rang
> pctquebec=quebec$Pourcentage/sum(quebec$Pourcentage)
> head(bretagne)
  Rang Patronyme Naissances
1    1   LE GALL      18126
2    2   LE GOFF      16093
3    3   LE ROUX      15905
4    4    THOMAS      15705
5    5    MARTIN      12094
6    6    TANGUY      12045

> minbretagne=tolower(bretagne$Patronyme)
> rangbretagne=bretagne$Rang
> pctbretagne=bretagne$Naissances/sum(bretagne$Naissances)
> I=minquebec%in%minbretagne
> sum(pctquebec[I])/sum(pctquebec)
[1] 0.2579641

Région	Proportion
Bretagne	23,39%
Normandie	31,07%
Loire	32,83%
Aquitaine	27,65%
Poitou	32,50%
Alsace	9,71%

Si on retrouve peu de noms alsacien au Québec, on peut être un peu “surpris” que les régions qui ressortent le plus sont la Loire et le Poitou et pas la Bretagne et la Normandie…. Bon, j’en conviens, je travaille sur les noms de familles bruts, c’est à dire que je ne tiens pas compte du fait que certains noms se sont déformés. Il s’agit de la version simple… on peut bien entendu aller plus loin, par exemple en utilisant des modèles de mélange par exemple… à suivre donc…

^*en fait, mes origines sont davantage bourguignonnes, de part mes quatre grands parents, comme le savent tous ceux qui m’ont déjà payé un verre de vin… mais ayant passé toute mon enfance en Normandie, je peux me considérer un peu comme normand. Et adorant la galette saucisse, je suis définitivement breton….

Computer

Playing with robots

03/05/2011 Arthur Charpentier 2 Comments

My son would be extremely proud if I tell him I can spend hours building robots. Well, my robots are not as fancy as Dr Tenma’s, but they usually do what I ask them to do. For instance, it is extremely simple to build a robot with R, to extract data from websites. I have mentioned it here (one tennis matches), but it failed there (on NY Marathon). To illustrate the use of robots, assume that one wants to build his own dataset to study prices of airline tickets. First, we have to choose a departure city (e.g. Paris) and an arrival city (e.g. Montreal). Then, one wants to look at all possible dates from April first (I ran it last month) till the end of December (so we create a vector with all leaving dates, namely a vector for the day, one for the month, and one for the year). Then, we choose a return date (say 3 days after).

DEP="Paris"
ARR="Montreal"
DATE1D=rep(c(1:30,1:31,1:30,1:31,1:31,1:30,1:31,1:30,
1:31,1:31,1:29),3)
DATE1M=rep(c(rep(4,30),rep(5,31),rep(6,30),rep(7,31),
rep(8,31),rep(9,30),rep(10,31),rep(11,30),rep(12,31),
rep(1,31),rep(2,29)),3)
DATE1Y=rep(c(rep(2011,30+31+30+31+31+30+31+
30+31+31+28),rep(2012,31+29)),3)
k=3
DATE3D=c((1+k):30,1:31,1:30,1:31,1:31,1:30,1:31,
1:30,1:31,1:31,1:29,1:k)
DATE3M=c(rep(4,30-k),rep(5,31),rep(6,30),rep(7,31),rep(8,31),
rep(9,30),rep(10,31),rep(11,30),rep(12,31),rep(1,31),rep(2,29),
rep(3,k))
DATE3Y=c(rep(2011,30+31+30+31+31+30+31+30+31+
31+28-k),re
p(2012,31+29+k))

It is also possible (for a nice robot), to skip all prior dates

skip=max(as.numeric(Sys.Date()-as.Date("2011-04-01")),1)

Then, we need a website where requests can be written nicely (with cities and dates appearing explicitly). Here, I cannot not mention the website that I used since it is stated on the website that it is strictly forbidden to run automatic requests… Anyway, consider a loop create a url address (actually I chose the value of the date randomly, since I had been told that those websites had memory: if you ask too many times for the same thing during a short period of time, prices would go up),

URL=paste("http://www.♦♦♦♦/dest.dll?qscr=fx&flag=q&city1=",
DEP,"&citd1=",ARR,"&",
"date1=",DATE1D[s],"/",DATE1M[s],"/",DATE1Y[s],
"&date2=",DATE3D[s],"/",DATE3M[s],"/",DATE3Y[s],
"&cADULT=1",sep="")

then, we just have to scan the webpage, looking for ticket prices (just looking for some specific names)

page=as.character(scan(URL,what="character"))
I=which(page%in%c("Price0","Price1","Price2"))
if(length(I)>0){
PRIX=substr(page[I+1],2,nchar(page[I+1]))
if(PRIX[1]=="1"){PRIX=paste(PRIX,page[I+2],sep="")}
if(PRIX[1]=="2"){PRIX=paste(PRIX,page[I+2],sep="")}

Here, we have to be a bit cautious, if prices exceed 1000. Then, it is possible to start a statistical study. For instance, if we compare to destination (from Paris), e.g. Montréal and New York, we obtain the following patterns (with high prices during holidays),

It is also possible to run the code twice (here it was run last month, and a couple of days ago), for the same destination (from Paris to Montréal),

Of course, it would be great if I could run that code say every week, to build up a nice dataset, and to study the dynamic of prices…

The problem is that it is forbidden to do this. In fact, on the website, it is mentioned that if we want to extract data (for an academic purpose), it is possible to ask for an extraction. But if we do tell that we study specific prices, data might be biased. So the good idea would be to use several servers, to make several requests, randomly, and to collect them (changing dates and destination). But here, my computing skills – unfortunately – reach a limit….

Computer, Statistics

Oscar awards: good actor versus good actress

02/05/2011 Arthur Charpentier 2 Comments

I am not a big fan of those ceremonies, where some actors pretend that they are extremely happy to be there, and then some win a trophy, some don’t, and those who win start to cry, and those who did not get a trophy try to pretend that they are not affected, etc. The other reason is that, since I have several kids, I do not go to see the movies that often (I mean apart from Shrek, Toy Story… Harry Potter is probably the only movie I’ve seen with real actors – or at least human actors).

But I remember being surprised when I looked at the nominees in newspapers,

Actresses are beautiful and look young, while actors are more experienced. So I have try to see how old were those who win an Oscar, as best actor (here) or best supporting actor (there), and best actress (here) and best supporting actress (there).

OSCAR=read.table("http://freakonometrics.blog.free.fr/public/data/OSCAR.csv",
sep=",",header=TRUE,dec=".")
actor=OSCAR[,1]
suppactor=OSCAR[,2]
actress=OSCAR[,3]
suppactress=OSCAR[,4]
actor=actor[is.na(actor)==FALSE]
actor=actor[actor>0]
actress=actress[is.na(actress)==FALSE]
actress=actress[actress>0]
suppactor=suppactor[is.na(suppactor)==FALSE]
suppactor=suppactor[suppactor>0]
suppactress=suppactress[is.na(suppactress)==FALSE]
suppactress=suppactress[suppactress>0]
 
boxplot(actor,suppactor,actress,suppactress,col=c("blue","blue","red","red"),
names=c("actor","supp. actor","actress","supp. actress"))

On average, a best actress is 36 years old, while a best actor is 44 years old. Which is quite a difference… Perhaps because it takes more time to an actor to be a good one ? Assuming that they start acting at 18, it takes 18 more years for an actress to be recognized as a good one (here the best one), and 26 for an actor. Or perhaps it is simply because leading actresses have to look young…
The oldest actor who won an Oscar was Henry Fonda (at the age of 76) and the oldest actress was Jessica Tendy (nearing 81). Tatum O’Neal became the youngest person to win the best suppo
rting actress award at the age of 10 (she was 8 when she was acting). The youngest best actress was Marlee Matlin, 21. The distribution was be seen below, with actors in blue, and actresses in red, best supporting actors in dotted lines, and best actors in plain lines,

plot(density(actor),xlim=c(10,80),axes=FALSE,
col="blue",names="",ylab="",xlab="",ylim=c(0,.051))
lines(density(suppactor),col="blue",lty=2)
lines(density(actress),col="red")
lines(density(suppactress),col="red",lty=2)
axis(1)

Note that the age of supporting actors is older that leading ones. E.g. the average age for supporting actors winning an Oscar is 50, while it is 44 for actors. Similarly, it is 40 for supporting actresses, and 36 for actresses.

> mean(suppactor)
[1] 50.23762
 
> mean(actor)
[1] 44.29982
 
> mean(suppactress)
[1] 40.55766
 
> mean(actress)
[1] 36.39733

Here, I have to admit that I was surprised. I always thought that being a supporting actor was a first step before being a leading one. So winners of supporting awards should have been younger that winners of leading ones. But this is not the case.

And the dynamic here is rather stable, with actors,

and actresses,

except that the age difference between supporting roles and leading roles have increased in the 80’s for actors, while it decreased in the 80’s for actresses.

Freakonometrics

Monthly Archives: May 2011