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.

Métadonnées et vie privée, l’équation insoluble?

Dans le cadre de la série d’articles (sur le thème des données) que j’édite sur variances.eu (le blog d’ENSAE Alumni), Yves-Alexandre de Montjoye –  @yvesalexandre – revient sur le thème métadonnées et vie privée, l’équation insoluble?

À l’image du texto, les comptes-rendus d’appels ou CRA, autrefois uniquement destinés à des fins de facturation, ont de loin dépassé leur but premier. Il est devenu clair au cours des derniers cinq ans, que ces métadonnées, littéralement “données à propos des données”, étaient utiles pour bien plus que juste la facturation. Connaître, même de manière anonyme, qui téléphone ou qui envoie un texto, à quel moment, et de quel endroit est une source d’information exceptionnelle sur nous-même et la société dans laquelle nous vivons.

Les comptes-rendus d’appels (CRA ou CDR en anglais) ont déjà été utilisés par des chercheurs pour étudier les vecteurs de propagation d’épidémies comme la malaria [1] ou de réaliser des recensements de population en temps réel [2]. Au-delà du monde de la recherche, l’utilisation commerciale de ces données est déjà en plein essor aux États-Unis et ailleurs, comme en témoigne le programme “precision insight” de Verizon ou FluxVision d’Orange.

Le potentiel de ces données ne doit cependant pas nous faire oublier que chacune d’entre elles est produite par une personne passant un appel, recevant un texto d’un ami ou effectuant un déplacement. Les CRA contiennent des informations détaillées et potentiellement sensibles sur le comportement, les habitudes de déplacement, ou encore le style de vie d’une personne.

Permettre la collecte et l’utilisation des données en respectant la vie privée des utilisateurs passe souvent, d’un point de vue pratique mais aussi légal, par leur “anonymisation”. L’idée est en effet que si les données ne sont pas associées à un individu, les informations qu’elles contiennent ne peuvent pas lui nuire.

[à suivre…]