Tag Archives: mélange

Maximum de Vraisemblance, approche numérique

Histoire de voir des choses plus amusantes dans le cours de statistique, on va aborder le cas où le maximum de la vraisemblance n’a pas d’expression analytique simple, pour insister sur les aspects numériques. Par exemple, si on a un mélange de lois, comme un mélange de lois normales,

Cette loi peut être utile, comme loi paramétrique, pour décrire la loi de la taille de gens. En guise d’illustration, on pourra utiliser les données déjà utilisée au début du cours

> Davis=read.table(
+ "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
> Davis[12,c(2,3)]=Davis[12,c(3,2)]
> X=Davis$height

La première étape est d’écrire la fonction de log-vraisemblance du mélange (ou plutôt l’opposé de la log-vraisemblance, les fonctions d’optimisation de R cherchant des minimums de fonctions)

> logL=function(theta){
+ p=theta[1]
+ m1=theta[2]
+ s1=theta[3]
+ m2=theta[4]
+ s2=theta[5]
+ logL=-sum(log(p*dnorm(X,m1,s1)+(1-p)*dnorm(X,m2,s2)))
+ return(logL)
+ }

La fonction la plus simple, pour maximiser la log-vraisemblance est la suivante (à partir d’un jeu de paramètre initiaux, histoire d’avoir un point de départ pour une descente de gradient)

> opt<-optim(c(.5,160,1,180,1),logL) > (theta=opt$par)
[1] 0.5987635 165.2547700 5.9410993 178.4856961 6.3547038

On reviendra dessus en cours car on peut faire ‘mieux’ en utilisant un algorithme d’optimisation sous contrainte, pour tenir compte du fait, par exemple, que la probabilité est forcément comprise entre 0 et 1 (c’est le cas ici, mais on a eu de la chance….).

Pour visualiser la densité estimée, on utilise

> hist(X,col="light green",probability=TRUE)
> lines(density(X))
> p=theta[1]
> m1=theta[2]
> s1=theta[3]
> m2=theta[4]
> s2=theta[5]
> v=p*dnorm(u,m1,s1)+(1-p)*dnorm(u,m2,s2)
> lines(u,v,col="blue",lwd=2)

Une autre solution est d’utiliser l’algorithme EM (sur lequel on reviendra en cours). Le point de départ part du constat que le problème de notre mélange est qu’on ne sait pas à quelle classe appartiennent nos gens. On va partir de valeurs initiales pour nos paramètres, et comparer les chances d’appartenir à chacunes des classes

> theta=c(.5,150,1,200,1)
> p1=dnorm(X,theta[2],theta[3])
> p2=dnorm(X,theta[4],theta[5])
> p=p1/(p1+p2)

A partir de ces probabilités d’appartenir à chacune des classes, on va estimer les paramètres de nos deux lois normales. Soit à l’aide de la formule fermée, soit en faisant du maximum de vraisemblance

> m1=sum(p*X)/sum(p)
> m2=sum((1-p)*X)/sum(1-p)
> s1=sqrt(sum(p*(X-m1)^2)/sum(p))
> s2=sqrt(sum((1-p)*(X-m2)^2)/sum(1-p))
> MlogL=function(a){
+ m1=a[1]
+ s1=a[2]
+ m2=a[3]
+ s2=a[4]
+ logL=-sum(log(p*dnorm(X,m1,s1)+(1-p)*dnorm(X,m2,s2)))
+ return(logL)
+ }
> opt<-optim(theta[2:5],MlogL)

L’idée est en fait d’avoir une boucle: on estime les probabilités d’appartenir dans les classes (compte tenu des paramètres des lois normales), et une fois qu’on a ces probabilités, on peut réestimer les paramètres des lois. Puis on recommence,

> theta=c(.5,150,1,200,1)
> for(s in 1:100){
+ p1=dnorm(X,theta[2],theta[3])
+ p2=dnorm(X,theta[4],theta[5])
+ p=p1/(p1+p2)
+ m1=sum(p*X)/sum(p)
+ m2=sum((1-p)*X)/sum(1-p)
+ s1=sqrt(sum(p*(X-m1)^2)/sum(p))
+ s2=sqrt(sum((1-p)*(X-m2)^2)/sum(1-p))
+ theta=c(mean(p),m1,s1,m2,s2)
+ }

On récupère alors les paramètres ‘optimaux’ de notre loi mélange

> hist(X,col="light green",probability=TRUE)
> lines(density(X))
> u=seq(130,210,by=.25)
> p=theta[1]
> m1=theta[2]
> s1=theta[3]
> m2=theta[4]
> s2=theta[5]
> v=p*dnorm(u,m1,s1)+(1-p)*dnorm(u,m2,s2)
> lines(u,v,col="red",lwd=2)

qui est relativement proche de notre estimation globale.

On reviendra au prochain cours sur ces méthodes ‘pratiques’ de calculs du maximum de vraisemblance.

Quelques heures d’effervescence sur Twitter

Pour ceux qui dormaient entre 11 heures et 15 heures (à Montréal, soit 17 heures – 21 heures en France), hier, Bernard Clavel est mort. Auparavant, j’allais sur le site de Libé, ou du Monde, pour me tenir informé de ce qui se passe en dehors de mon bureau (ou de ma cuisine),

Maintenant je suis sur Twitter… et autant l’avouer, j’ai essayé des outils comme TweetDeck mais au bout de 2 heures, force est de constater que j’ai arrêté l’application car j’en avais marre d’être harcelé par ces fenêtres…
Sur Twitter, tout le monde a quelque chose à dire sur tout, et ne peut s’empêcher de vouloir le dire à tout le monde…. Par exemple “Bernard Clavel” ne faisait pas franchement parti des mots tapé régulièrement sur Twitter…. Mais à l’annonce de sa mort, tout le monde veut diffuser l’information à son réseau…. Ça allait tellement vite que j’étais obligé de faire des backups des pages internets tous les quarts d’heure (sinon l’historique disparaît, ce qui limite l’intérêt des outils statistiques). Pour être précis, ma recherchait était

http://twitter.com/#search?q=bernard%20clavel

et pendant près de 4 heures, ça allait du tweet original, au retweet, au message recopié sans citation, etc…. J’ai récupéré ainsi plus de 400 tweets, mais sans l’horaire…. Finalement, j’ai découvert

http://search.tweetreports.com/search_results.php?search=bernard+clavel&sear=Y

Pendant l’heure qui suit les premiers tweets, c’est l’effervescence, puis ça retombe… avec un pic à 20 heures, visiblement certains l’apprenant aux informations du soir…

Manifestement la loi gamma (la date initiale étant le premier tweet annonçant le décès) n’est pas adaptée…

Moralité, pour poursuivre le billet de la semaine dernière (ici), on peut tenter un mélange, en l’occurrence un mélange gamma/exponentielle…. la première loi gamma permet d’intégrer le délai de transmission de l’information au réseau, alors que pour le second, il semble que ce soit une annonce au 20 heures qui ait déclenché la soudaine envie de tweeter..

Allez, le jour où je maîtrise mieux les outils linguistiques, je regarde en plus quand les tweets anglophone et hispaniques sont arrivés…..

Faut-il se fâcher avec tout le monde ?

(s’il le faut allons-y) Samedi matin, j’ai aperçu avec un grand plaisir le dernier numéro de Fakir chez mon buraliste préféré (i.e. le “journal fâché avec tout le monde“). Je n’ai pas vocation à me fâcher avec tout le monde, loin de là, mais je dois avoir un côté pitbull que j’ai du mal à cacher…

Et si quelques pages ont rapidement retenu mon attention, je dois avouer que j’en suis sorti très déçu (par ces quelques pages). Bon, autant le dire tout de suite, les pages en questions ne sont pas écrite par les journalistes de Fakir, mais correspondent à un “courrier des lecteur” (les “meilleurs” paraît-il… je me sens flaté, je suis un lecteur régulier). Le “papier” de deux pages part d’une attention louable (me semble-t-il), à savoir essayer de repérer des classes (au sens marxiste) en étudiant la distribution des revenus. On nous explique que l’auteur est “directeur de recherche au CNRS“, ce qui vise à donner a priori une légitimité scientifique à ses propos. N’ayant pas encore l’accord pour mettre en ligne les deux pages en question (je vais quand même demander l’accord à la rédaction), je me contente de ces deux copies de mauvaise qualité.

Pour résumer, l’auteur nous explique que l’INSEE ment, ou triche, dans sa façon de produire ses graphiques (comme je l’avais dit ici, je rejoins l’auteur sur ce point: j’avais déjà dit que la façon de présenter des informations chiffrées n’était pas neutre sur le message que l’on souhaite diffuser).

  • Densité ou fonction quantile ?

L’auteur commence par nous présenter deux graphiques. Le premier tiré de Max Weber montre la répartition de la taille de fourmis (géantes !), et observe que la distribution bimodale s’interprète comme la présence de deux classes de fourmis.

Ensuite, il nous présente la fonction quantile du revenu des français, et s’étonne de ne pas retrouver de distribution bimodale.

Il nous explique qu’avec ce genre de graphique, on ne trouverait pas les modes observés sur les tailles des fourmis,

plus précisément, l’auteur explique

Pour formaliser un peu (et revenir sur ce qu’est la “moulinette de l’insee“), rappelons le premier graphique est une densité https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-01.png (je passe sur le fait que l’on travaille sur des échantillons, je vais faire une explication sur la modélisation sous-jacente, i.e. en travaillant sur la densité sous-jacente). Pour interpréter le second, rappelons que l’on construit la fonction de répartition en intégrant la densité,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-02.png

La fonction quantile est alors l’inverse de cette fonction

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-03.png

autrement dit, d’un point de vue géométrique, le symétrique par rapport à la première bissectrice de la fonction de répartition.  Avec des graphiques un peu plus clairs, je pense, voici une densité de mélange, en noir

Les deux populations sont visualisées avec les deux couleurs, rouge et bleu. Si je trace la fonction cumulative, j’obtiens

Malheureusement, on ne dispose pas de données précises sur les revenus (et je suis le premier à m’en désoler), et on observe une version discrète de cette fonction, à savoir les niveaux des déciles,

En inversant cette focntion, on obtient alors la fonction quantile, correspondant au graphique publié par l’insee,

Ma première observation est que si https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-01.png est effectivement bimodale, cela signifie que la dérivée de https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-04.png possède deux modes, ou plutôt deux maximum locaux. Si on se souvient des trucs que l’on apprend au lycée (et non pas à un niveau “bac+5” comme l’indique le haut du document), si on fait le tableau de signe (la première ligne étant le signe de la dérivée de la densité) on obtient

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-tableau-signe.png

Autrement dit, on ne voit rien. Formellement, cela se traduira uniquement sur la pente de https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-04.png qui diminuera après le premier mode, qui augmentera juste avant d’atteindre le deuxième, avant de diminuer à nouveau après le second. Sur les courbes en escaliers, c’est vrai qu’il est difficile de visualiser une distribution bimodale, i.e. des accroissements sur la “pente” (ce n’était déjà pas forcément simple avec l’ensemble des valeurs possibles, mais sur une fonction en escalier….). Mais quoi qu’il en soit, je ne suis pas sûr que cette recherche de la bimodalité soit fondamental dans cette étude….

  • Bimodalité et mélange

En fait, ma seconde critique est un peu plus technique. Pour faire simple, la présence de “classes” de traduit par une présence de “mélange” (c’est la base de la tarification où l’on essaye d’identifier les classes tarifa
ires, en assurance, comme j’en ai parlé à de nombreuses reprises ici ou , par exemple). Mais ce “mélange” n’a aucune raison de se traduire par une distribution bimodale. On retrouve ici un argument fallacieux que l’on retrouve couramment, y compris dans des revues scientifiques, par exemple dans ce papier sur l’âge des tentatives de suicides, publié dans Journal of Psychiatric Research (ici). On y retrouve une tentative de retrouver des mélanges,

(je reviendrais un jour sur cette étude). Pour essayer de mieux comprendre en quoi l’argument est fallacieux, regardons quelques cas simples.
Le premier point est que la densité d’un mélange n’a aucune raison d’être bimodale. Même avec des lois sous-jacentes unimodales, comme la loi normale par exemple,

La courbe en noir est – comme auparavant – la distribution sur l’ensemble de la population. Elle est bimodale si les groupes sont bien discriminés, mais pas si les groupes sont proches. Comme je l’explique dans toutes mes interventions sur  les modèles linéaires généralisés, il est difficile de se faire une idée sur la loi par sous-groupe en regardant la distribution globale. Si la loi n’est pas unimodale (à en croire les travaux de Vilfredo Pareto – mentionnés ici ou – j’aurais tendance à croire que c’est le cas quand on parle des revenus), cette recherche des modes a encore moins de sens, me semble-t-il…. Le graphique ci-dessous correspond à mélanges de lois exponentielles (les traits en pointillés correspondant, comme auparavant, aux moyennes par sous groupe).

De plus, pour revenir aux exemples mentionnés dans le papier, à on avis, il y a une différence fondamentale entre le revenu et la taille: si la taille peut être supposé suivre une loi normale (par sous-population) je doute que ce soit le cas pour le revenu…Et encore, j’ai fait simple, car on suppose qu’il existe uniquement deux classes, à peu près équiréparties (ici 40-60 dans les deux animations précédentes).

Sur l’exemple ci-dessus, on suppose qu’il existe un petit groupe  hétérogène (i.e. avec une plus forte variance). Cette hypothèse ne me semble pas absurde si on raisonne sur les revenus. Bref, on ne trouve presque jamais de distribution bimodale… Et ce d’autant plus que le petit groupe est hétérogène. Si en plus on suppose qu’il existe non pas deux groupes, mais trois (voire plus), c’est encore pire: il devient difficile de visualiser les sous-groupes…

  • Faute de logique ou de statistique ?

J’ai l’impression que l’auteur part d’une affirmation qui me semble juste, à savoir “une distribution bimodale peut être modélisée par un mélange de deux distributions unimodales“. Sauf qu’il fait un raisonnement étrange en réécrivant cette affirmation sous la forme “une distribution qui n’est pas bimodale ne peut pas être un mélange” (je passe sur le premier point que j’avais mentionné comme quoi chercher un mélange sur une discrétisation de la fonction quantile, ou de la fonction de répartition, me semble étrange). Bref, je pense que l’ensemble de l’étude statistique repose sur des arguments fallacieux, ce qui a malheureusement tendance à décrébiliser l’ensemble de l’analyse. C’est dommage… Pourquoi ne pas avoir fait appel à un chercheur connaissant statistiques pour faire l’étude ? Mais beaucoup sont passés par l’INSEE, ce qui semblerait ôter tout légitimité à être critique (à en croire le petit encadré en bas)… Je crois qu’il serait temps d’arrêter de stigmatiser ainsi l’insee, beaucoup de gens y travaillant ayant gardé un esprit ouvert et critique ! Il serait dommage de ne pas s’appuyer sur leur expertise. L’autre regret est que la caution intellectuelle viennent d’un ” directeur de recherche du CNRS” (éventuellement en linguistique ou en anthropologie, on n’a malheureusement pas plus d’information).  Je trouve regrettable qu’un titre académique permette de se dédouaner de toute relecture critique… par des statisticiens….