Category Archives: Inference

Statistiques, de la théorie à la pratique (partie 2)

Il y a une quinzaine de jours (ici) on avait vu comment construire des estimateurs pour l’indice de queue de la loi de Pareto. On avait 4 valeurs (différentes). Maintenant, on peut se demander lequel devrait être le “meilleur“.

  • propriétés de l’estimateur du maximum de vraisemblance

Maintenant que l’on a un peu avancé dans le cours, on peut continuer l’étude de la distribution des couts de sinistres (étude commencée ici).

http://freakonometrics.blog.free.fr/public/perso2/pareto-100.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso2/pareto-101.gif

de telle sorte que

http://freakonometrics.blog.free.fr/public/perso2/pareto-102.gif

et

http://freakonometrics.blog.free.fr/public/perso2/pareto-204.gif

Aussi l’information de Fisher est http://freakonometrics.blog.free.fr/public/perso2/pareto-203.gif.
La variance aymptotique de http://freakonometrics.blog.free.fr/public/perso2/pareto-205.gif est alors

http://freakonometrics.blog.free.fr/public/perso2/pareto-206.gif
  • propriétés de l’estimateur de la méthode des moments

Nous avions noté que

http://freakonometrics.blog.free.fr/public/perso2/pareto-207.gif

i.e. http://freakonometrics.blog.free.fr/public/perso2/pareto-208.gif, alors http://freakonometrics.blog.free.fr/public/perso2/pareto-209.gif. Aussi, comme le théorème central limite limite

http://freakonometrics.blog.free.fr/public/perso2/pareto-210.gif

http://freakonometrics.blog.free.fr/public/perso2/pareto-211.gif

la delta method permet d’écrire, en posant http://freakonometrics.blog.free.fr/public/perso2/pareto-311.gif,

http://freakonometrics.blog.free.fr/public/perso2/correc-latex-pareto.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso2/correc-latex-pareto2.gif

La variance asymptotique est alors

http://freakonometrics.blog.free.fr/public/perso2/pareto-310.gif
  • propriétés de l’estimateur de la méthode de la médiane

Pour le dernier estimateur, le théorème de Glivenko Cantelli sur les fonctions quantiles permet d’écrire

http://freakonometrics.blog.free.fr/public/perso2/pareto-309.gif

i.e. pour la médiane,

http://freakonometrics.blog.free.fr/public/perso2/pareto-308.gif

où la variance asymptotique (pour la médiane) peut se simplifier sous la forme

http://freakonometrics.blog.free.fr/public/perso2/pareto-307.gif

Là encore, en utilisant la delta method (je peux renvoyer ici pour une application numérique de cette méthode), avec

http://freakonometrics.blog.free.fr/public/perso2/pareto-306.gif

i.e. http://freakonometrics.blog.free.fr/public/perso2/pareto-305.gif, on en déduit que

http://freakonometrics.blog.free.fr/public/perso2/pareto-304.gif

La variance asymptotique de notre estimateur est alors (en enlevant le facteur http://freakonometrics.blog.free.fr/public/perso2/IC-std-6.gifau dénominateur – car il est commun à tous),

http://freakonometrics.blog.free.fr/public/perso2/pareto-303.gif

Pour trois de nos estimateurs (je n’ai pas fait le quatrième, ça serait hors programme), on a pu obtenir facilement des variances asymptotiques. Graphiquement, on a les évolutions suivantes pour les variances asymptotiques de nos trois estimateurs, en fonction de la valeur de http://freakonometrics.free.fr/blog/latex/information09.gif, avec le maximum de vraisemblance (en rouge), celui de la méthode des moments (en bleu) et celui basé sur la médiane (en mauve)

Manifestement le premier estimateur (le maximum de vraisemblance) a une variance asymptotique plus faible que les deux autres (mais on l’avait vu en cours, c’est ce que garantit l’inégalité de Cramér-Rao). Autrement dit, l’estimateur du maximum de vraisemblance est probablement le plus fiable….
En utilisant ces résultats mentionnés, on peut aussi construire des intervalles de confiance asymptotiques pour la valeur du paramètre http://freakonometrics.free.fr/blog/latex/information09.gif, à l’aide de nos trois estimateurs, i.e. pour un intervalle de confiance à 90%,

http://freakonometrics.blog.free.fr/public/perso2/pareto-302.gif

où les estimations des variances asymptotiques sont obtenues en substituant http://freakonometrics.blog.free.fr/public/perso2/pareto-205.gifà  http://freakonometrics.free.fr/blog/latex/information09.gif dans les formules au dessus.
Si on considère que l’estimateur du maximum de vraisemblance est le meilleur que l’on puisse utiliser, on obtient comme intervalle de confiance pour http://freakonometrics.free.fr/blog/latex/information09.gif,

2.845093-1.64/sqrt(100)*2.845093
[1] 2.378498
 
2.845093+1.64/sqrt(100)*2.845093
[1] 3.311688

L’étape d’après sera d’introduire des tests. On pourra ainsi se demander si l’on peut accepter l’hypothèse http://freakonometrics.blog.free.fr/public/perso2/pareto-301.gif. Mais c’est une autre histoire… à suivre donc….

Estimateurs des moments

Dans le cours, nous avions vu que si http://freakonometrics.hypotheses.org/files/2015/12/moment-est.gif, un estimateur naturel pour le paramètre était http://freakonometrics.hypotheses.org/files/2015/12/moment-est2.gif. On appelle cette méthode la méthode des moments. Mais il est possible, que, plus généralement http://freakonometrics.hypotheses.org/files/2015/12/mmomot3.gif, et on pourrait considérer comme estimateur naturel

http://freakonometrics.hypotheses.org/files/2015/12/kvmoment4.gif

On avait noté en cours que le premier soucis était un soucis d’inversibilité. Pour utiliser cette méthode, il faut une relation bijective entre le moment et le paramètre. On avait vu avec la loi binomiale – que ce n’est pas le cas pour la variance, par exemple.

Le second soucis est un problème numérique: il faut inverser cette fonction (ce qui peut ne pas être trivial) mais avec une méthode de type sécante (programméeici), ça pourrait marcher. Par exemple pour la loi de Poisson, on sait que

http://freakonometrics.hypotheses.org/files/2015/12/poisson-mu-1.gif
http://freakonometrics.hypotheses.org/files/2015/12/poisson-m-2.gif
http://freakonometrics.hypotheses.org/files/2017/10/poisson-m-3.gif
http://freakonometrics.hypotheses.org/files/2015/12/poisson-m4.gif

(…etc). On peut numériquement construire des estimateurs de la manière suivante

secant=function(fun, x0, x1, tolerence=1e-07, niter=500){
for ( i in 1:niter ) {
	x2 <- x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0))
	if (abs(fun(x2)) < tolerence)
		return(x2)
	x0 <- x1
	x1 <- x2
}}
X=rpois(50,10)
m1=mean(X)
f=function(x){x*(1+x)-mean(X^2)}
secant(f,0,40)
m2=secant(f,0,50)
f=function(x){x*(1+3*x+x^2)-mean(X^3)}
secant(f,0,50)
m3=secant(f,0,50)
f=function(x){x*(1+7*x+6*x^2+x^3)-mean(X^4)}
secant(f,0,100)
m4=secant(f,0,20)

Le dernier soucis est un problème plus statistique: il semble que l’estimateur basé sur le premier moment soit meilleur que les autres. Sur le dessin ci-dessous, on a simulé 10,000 fois des lois de Poisson de paramètre 10 (qui est la vraie valeur du paramètre, ici représentée verticalement en pointillés), et on représente en rouge l’estimateur construit sur le premier moment, avec la croix au centre qui représente la valeur moyenne, et le segment qui est un intervalle de confiance à 80% de notre estimateur (i.e. les quantiles à 10% et 90% obtenus sur les 10,000 simulations),

En bleu, on a l’estimateur construit à l’aide du moment d’ordre deux, en vert du moment d’ordre trois, et en mauve du moment d’ordre quatre. Bref, plus on prend des moments d’ordre élevé, plus le biais semble important, et plus la variance de notre estimateur augmente.

Moralité, l’estimateur de la méthode des moments (et l’unique) sera construit sur les moments les plus faibles possibles…

Statistiques, de la théorie à la pratique

Afin d’illustrer un peu le cours, je mets en ligne ici. Il s’agit de coûts de sinistres ayant dépassé 1,000$ (qui ont été ordonnés de manière croissante). Pour importer la base en R, le code est simplement,

> X=read.table(
"http://freakonometrics.blog.free.fr/public/data/sinistres.txt")$x
> X
 [1] 1.003 1.016 1.023 1.027 1.037 1.039 1.0
[13] 1.061 1.072 1.078 1.082 1.087 1.094 1.0
[25] 1.110 1.112 1.117 1.132 1.138 1.141 1.1
[37] 1.180 1.186 1.187 1.190 1.193 1.203 1.2
[49] 1.321 1.326 1.338 1.342 1.343 1.344 1.3
[61] 1.428 1.432 1.442 1.457 1.463 1.466 1.4
[73] 1.551 1.553 1.566 1.584 1.632 1.695 1.6
[85] 1.881 1.893 1.897 1.958 2.000 2.175 2.2
[97] 3.045 3.103 4.495 5.614

Pour information, l’histogramme ressemble à ça

On va supposer que l’on peut ajuster une loi de Pareto, dont la fonction de survie est donnée par

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-1.gif

et de densité

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-2.gif

L’espérance est

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-3.gif

et la variance

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-5.gif
  • Estimateur de la méthode des moments
http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-6.gif

Numériquement, on obtient

> mean(X)/(mean(X)-1)
[1] 2.946827
  • Estimateur du maximum de vraisemblance

On peut écrire la vraisemblance, et on recherche – numériquement – la maximum (ou le minimum de l’opposé de la log-vraisemblance car la fonction recherche ici les minimum),

> f=function(x,a){
+ a/x^(a+1)
+ }
> LogV = function(a,echantillon){
+ -sum(log(f(echantillon,a)))
+ }
> optim(fn=LogV,par=2,echantillon=X,method="BFGS")$par
[1] 2.845093

Mais on aurait pu faire mieux… car ici le maximum de vraisemblance peut s’écrire analytiquement

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-7.gif
> length(X)/(sum(log(X)))
[1] 2.845093

(tout va bien, on trouve pareil).

  • Estimateur par régression linéaire

On peut noter que

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-8.gif

ce qui veut dire qu’on peut tracer le nuage de points du log des couts, et s’ils sont alignés suivant une droite (passant par l’origine), la pente nous donnera un estimateur de notre paramètre. Le plus simple est de prendre la pente de la droite qui passe par l’origine, et par le centre de gravité du nuage,

> Z=X
> Y=log((length(Z):1)/length(Z))
> X=log(Z)
> plot(X,Y)
> -mean(Y)/mean(X)
[1] 2.753414
  • Estimation à l’aide de la médiane

Enfin, notons que la médiane est

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-9.gif

Un estimateur naturel est alors

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-11.gif

Numériquement, c’est alors tout simple,

> log(2)/log(quantile(X,.5))
     50% 
2.413719

Moralité, on a quatre estimateurs, qui nous donne quatre valeurs numériques différentes

2.9468 2.8451 2.7534 2.4137

Pour ceux qui veulent trouver d’autres estimateurs (car on peut en construire des centaines), une référence intéressante est la suivante,
Quandt, R. E. (1966). Old and new methods of estimation and the Pareto distribution, Metrika, 10, 55-82 (ici)
Le but de la statistique mathématique est de mieux connaître les propriétés (fondamentales) de ces estimateurs afin de choisir le “meilleur“. Et c’est un problème important, par exemple si on souhaite mettre en place un contrat de réassurance (et transférer les très grands sinistres). Par exemple si on veut transférer les sinistres dépassant 100,000$, la probabilité de dépasser ce seuil est 10 fois plus grande avec le dernier estimateur par rapport au premier…

> 1/100^2.4137
[1] 1.48799e-05
 
> 1/100^2.9468
[1] 1.277615e-06

Moralité, les résultats théoriques que l’on aborde en ce moment vont nous aider à résoudre des questions très pratiques…. à suivre.

Information de Fisher

Nous avons défini dans le cours l’information de Fisher par

http://freakonometrics.free.fr/blog/latex/information01.gif

où http://freakonometrics.free.fr/blog/latex/information02.gif. On verra la semaine prochaine que cette information peut s’écrire (moyennant quelques conditions de régularité)

http://freakonometrics.free.fr/blog/latex/information03.gif

Dans un cadre multivarié, on a que http://freakonometrics.free.fr/blog/latex/information04.gif avec

http://freakonometrics.free.fr/blog/latex/information05.gif

ou encore

http://freakonometrics.free.fr/blog/latex/information06.gif

(avec une écriture analogue au cas univarié en utilisant la matrice Hessienne).
Étant donné un échantillon, l’information de Fisher empirique est la pente de la log-vraisemblance en http://freakonometrics.hypotheses.org/files/2015/12/information09.gif, i.e. dans le cas d’échantillons de 20 tirages de pile/face (avec une pièce non truquée)

http://freakonometrics.blog.free.fr/public/perso2/info-fisher.gif
X=sample(c("PILE","FACE"),size=20,prob=c(.5,.5),replace=TRUE)
X=(X=="PILE")*1
vecP=seq(0,1,by=.01)
logL=function(p){sum(log(dbinom(X,size=1,prob=p)))}
L=lapply(vecP,logL)
plot(vecP,L,type="l",col="blue",lwd=2,xlab="",ylab="")
abline(v=.5,col="red",lty=2)

En répétant plusieurs milliers de fois, on peut étudier la distribution de cette pente

En note que le score est centré, et ici la variance est

> var(SCORE)
[1] 80.00168
> 20/(.5*(1-.5))
[1] 80
20/(.5*(1-.5))

(ce qui correspond effectivement à l’information de Fisher). Mais cette informationest un concept plus général qu’une variance…
Bon, maintenant essayons de creuser un peu cette histoire d’information (c’est Ronald Fisher qui a – le premier – utilisé ce nom en 1925, bien avant les travaux de Shannon et Wiener sur des mesures logarithme d’information dans un contexte de théorie de la communication).

On cherche a quantifier le mesure d’information d’un échantillon http://freakonometrics.free.fr/blog/latex/information07.gif correspondant à une collection de variables aléatoires http://freakonometrics.free.fr/blog/latex/information08.gif, i.i.d. de loi http://freakonometrics.free.fr/blog/latex/information10.gif, et de densité/loi de probabilité http://freakonometrics.free.fr/blog/latex/information11.gif.

Une hypothèse est qu’avoir http://freakonometrics.free.fr/blog/latex/information12.gif observations apporte deux fois plus d’information qu’en avoir seulement http://freakonometrics.free.fr/blog/latex/information13.gif. On a alors quelque chose qui pourrait s’écrire

http://freakonometrics.free.fr/blog/latex/information14.gif

Aussi on suppose que pas d’observation signifie pas d’information… J’ai l’impression que les bayésiens ne seraient pas d’accord, mais c’est une autre histoire…

On suppose aussi que l’information dépend de la loi  http://freakonometrics.free.fr/blog/latex/information11.gif, et pas de l’échantillon http://freakonometrics.free.fr/blog/latex/information07.gif.

On peut aussi noter (c’est une idée peut être plus simple à comprendre dans le contexte de la communication) que la compression ne peut que détériorer l’information. Autrement dit, si on utiliser http://freakonometrics.blog.free.fr/public/perso2/infofi-01.gif et non pas http://freakonometrics.blog.free.fr/public/perso2/infori-02.gif, l’information est moins bonne, i.e.

http://freakonometrics.blog.free.fr/public/perso2/infofi-03.gif

Et le seul cas pour lequel on n’a pas de perte d’information est le cas d’une statistique exhaustive.

On notera aussi qu’une information doit être positive (au pire on n’a pas d’information).

Mais pour relier cette information a une propriété sur la variance de notre estimateur, il faut creuser un peu. Considérons un estimateur sans biais http://freakonometrics.free.fr/blog/latex/information15.gif de http://freakonometrics.free.fr/blog/latex/information09.gif, i.e.

http://freakonometrics.free.fr/blog/latex/infromation16.gif

ou encore (en supposant que l’échantillon soit réduit à une unique observation)

http://freakonometrics.free.fr/blog/latex/information17.gif

soit, en notant que

http://freakonometrics.free.fr/blog/latex/information18.gif

(car on a une densité)

http://freakonometrics.free.fr/blog/latex/information19.gif

En différenciant par rapport à http://freakonometrics.free.fr/blog/latex/information09.gif, on peut en déduire que

http://freakonometrics.free.fr/blog/latex/information20.gif

soit (en admettant que l’on puisse intervertir intégrale et dérivée)

http://freakonometrics.free.fr/blog/latex/information21.gif

Or on peut réécrire

http://freakonometrics.free.fr/blog/latex/information22.gif

ce qui donne

http://freakonometrics.free.fr/blog/latex/information23.gif

ou encore

http://freakonometrics.free.fr/blog/latex/information24.gif

Et là, on invoque l’inégalité de Cauchy-Schwarz afin de garantir que

http://freakonometrics.free.fr/blog/latex/information25.gif

On reconnait alors deux espérances, i.e.

http://freakonometrics.free.fr/blog/latex/information26.gif

(pour le terme de gauche, on a utilisé le fait qu’on avait un estimateur sans biais). Bref, il existe plein d’application de cette information, cette dernière inégalité étant probablement la plus importante.
Pour aller un peu plus loin, je peux renvoyer vers des notes de cours de Tran Viet Chi, en ligne ici. Et pour ceux qui veulent aller vraiment plus loin, John Pierce Robert Ash et surtout Solomon Kullback ont écrit des livres entiers sur statistiques et mesure d’information….

Estimer une variance… ou un écart-type ?

Lors du cours de demain, nous allons parler (longuement) de l’estimation de la variance (comme cela avait commencé en TD vendredi dernier). Mais en pratique, c’est souvent l’écart-type qui nous intéresse. On va passer par mal de temps à chercher un estimateur sans biais de la variance, et montrer que

devrait convenir, mais cela ne nous donne pas vraiment d’information sur le biais de l’estimateur naturel de l’écart-type, à savoir

car http://freakonometrics.blog.free.fr/public/maths/variance-est03.png.
Bref, il y a sûrement un biais, et ça pourrait être intéressant de le quantifier. Classiquement, on suppose disposer d’un échantillon http://freakonometrics.free.fr/blog/sigma01.png i.i.d. On va faire une hypothèse forte, à savoir que leur loi est une loi normale http://freakonometrics.free.fr/blog/sigma03.png, où http://freakonometrics.free.fr/blog/sigma02.png est inconnu.
Naturellement, on voudrait considérer

http://freakonometrics.free.fr/blog/sigma04.png

On va poser, afin de simplifier un peu,

http://freakonometrics.free.fr/blog/sigma07.png
et on va chercher à calculer

http://freakonometrics.free.fr/blog/sigma06.png

où http://freakonometrics.free.fr/blog/sigma08.png est la densité de cette variable. Or on sait que http://freakonometrics.free.fr/blog/sigma20.png suit une loi du chi-deux à http://freakonometrics.free.fr/blog/sigma11.png degrés de liberté, dont la densité est la loi gamma suivante

http://freakonometrics.free.fr/blog/sigma13.png

En notant que l’on peut réécrire

http://freakonometrics.free.fr/blog/sigma22.png

on a que cette espérance s’écrit

http://freakonometrics.free.fr/blog/sigma24.png

soit en faisant un peu de calculs,

http://freakonometrics.free.fr/blog/sigma26.png

ou encore

http://freakonometrics.free.fr/blog/sigma27.png

Ce qui va se simplifier en utilisant la fonction gamma

http://freakonometrics.free.fr/blog/sigma28.png

On a ainsi

http://freakonometrics.blog.free.fr/public/perso/stdvariance.png

http://freakonometrics.free.fr/blog/sigma35.png

EM and mixture estimation

Following my previous post on optimization and mixtures (here), Nicolas told me that my idea was probably not the most clever one (there).
So, we get back to our simple mixture model,

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

In order to describe how EM algorithm works, assume first that both https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png and  https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.pngare perfectly known, and the mixture parameter is the only one we care about.

  • The simple model, with only one parameter that is unknown

Here, the likelihood is

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

so that we write the log likelihood as

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

which might not be simple to maximize. Recall that the mixture model can interpreted through a latent variate https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png (that cannot be observed), taking value when https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM07.png is drawn from https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png, and 0 if it is drawn from https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png. More generally (especially in the case we want to extend our model to 3, 4, … mixtures), https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM08.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM09.png.
With that notation, the likelihood becomes

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

and the log likelihood

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

the term on the right is useless since we only care about p, here. From here, consider the following iterative procedure,
Assume that the mixture probability https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png is known, denoted https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png. Then I can predict the value of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png (i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM08.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM09.png) for all observations,

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

So I can inject those values into my log likelihood, i.e. in

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

having maximum (no need to run numerical tools here)

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

that will be denoted https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM17.png. And I can iterate from here.
Formally, the first step is where we calculate an expected (E) value, where https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.pngis the best predictor of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM19.png given my observations (as well as my belief in https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png). Then comes a maximization (M) step, where using https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png, I can estimate probability https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png.

  • A more general framework, all parameters are now unkown

So far, it was simple, since we assumed that https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png and  https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png were perfectly known. Which is not reallistic. An there is not much to change to get a complete algorithm, to estimate https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM30.png. Recall that we had https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png which was the expected value of Z_{1,i}, i.e. it is a probability that observation i has been drawn from https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png.
If https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png, instead of being in the segment https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM31.png was in https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM32.png, then we could have considered mean and standard deviations of observations such that https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png=0, and similarly on the subset of observations such that https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png=1.
But we can’t. So what can be done is to consider https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png as the weight we should give to observation i when estimating parameters of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png, and similarly, 1-https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.pngwould be weights given to observation i when estimating parameters of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png.
So we set, as before

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

and then

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM34.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM35.png

and for the variance, well, it is a weighted mean again,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM36.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM37.png

and this is it.

  • Let us run the code on the same data as before

Here, the code is rather simple: let us start generating a sample
> X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z  = sample(c(1,2,2),size=n,replace=TRUE)
> X2=4+X20
> X = c(X1[Z==1],X2[Z==2])
then, given a vector of initial values (that I called https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png and then https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM99.png before),
>  s = c(0.5, mean(X)-1, var(X), mean(X)+1, var(X))
I define my function as,
>  em = function(X0,s) {
+  Ep = s[1]*dnorm(X0, s[2], sqrt(s[4]))/(s[1]*dnorm(X0, s[2], sqrt(s[4])) +
+  (1-s[1])*dnorm(X0, s[3], sqrt(s[5])))
+  s[1] = mean(Ep)
+  s[2] = sum(Ep*X0) / sum(Ep)
+  s[3] = sum((1-Ep)*X0) / sum(1-Ep)
+  s[4] = sum(Ep*(X0-s[2])^2) / sum(Ep)
+  s[5] = sum((1-Ep)*(X0-s[3])^2) / sum(1-Ep)
+  return(s)
+  }
Then I get https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png, or https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM99.png. So this is it ! We just need to iterate (here I stop after 200 iterations) since we can see that, actually, our algorithm converges quite fast,
> for(i in 2:200){
+ s=em(X,s)
+ }

Let us run the same procedure as before, i.e. I generate samples of size 200, where difference between means can be small (0) or large (4),

Ok, Nicolas, you were right, we’re doing much better ! Maybe we should also go for a Gibbs sampling procedure ?… next time, maybe….

Optimization and mixture estimation

Recently, one of my students asked me about optimization routines in R. He told me he that R performed well on the estimation of a time series model with different regimes, while he had trouble with a (simple) GARCH process, and he was wondering if R was good in optimization routines. Actually, I always thought that mixtures (and regimes) was something difficult to estimate, so I was a bit surprised…

Indeed, it reminded me some trouble I experienced once, while I was talking about maximum likelihooh estimation, for non standard distribution, i.e. when optimization had to be done on the log likelihood function. And even when generating nice samples, giving appropriate initial values (actually the true value used in random generation), each time I tried to optimize my log likelihood, it failed. So I decided to play a little bit with standard optimization functions, to see which one performed better when trying to estimate mixture parameter (from a mixture based sample). Here, I generate a mixture of two gaussian distributions, and I would like to see how different the mean should be to have a high probability to estimate properly the parameters of the mixture.

The density is here https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-01.png proportional to

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

The true model is https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-03.png, and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-04.png being a parameter that will change, from 0 to 4.
The log likelihood (actually, I add a minus since most of the optimization functions actually minimize functions) is
> logvraineg <- function(param, obs) {
+ p <- param[1]
+ m1 <- param[2]
+ sd1 <- param[3]
+ m2 <- param[4]
+  sd2 <- param[5]
+  -sum(log(p * dnorm(x = obs, mean = m1, sd = sd1) + (1 – p) *
+ dnorm(x = obs, mean = m2, sd = sd2)))
+  }
The code to generate my samples is the following,
>X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z  = sample(c(1,2,2),size=n,replace=TRUE)
> X2=m+X20
> X = c(X1[Z==1],X2[Z==2])
Then I use two functions to optimize my log likelihood, with identical intial values,
> O1=nlm(f = logvraineg, p = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)), obs = X)
> logvrainegX <- function(param) {logvraineg(param,X)}
> O2=optim( par = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)),
+   fn = logvrainegX)
Actually, since I might have identification problems, I take either https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-05.png or https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-06.png, depending whether https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-07.png or https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-08.png is the smallest parameter.

On the graph above, the x-axis is the difference between means of the mixture (as on the animated grap above). Then, the red point is the median of estimated parameter I have (here https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-05.png), and I have included something that can be interpreted as a confidence interval, i.e. where I have been in 90% of my scenarios: theblack vertical segments. Obviously, when the sample is not enough heterogeneous (i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-09.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-04.png rather different), I cannot estimate properly my parameters, I might even have a probability that exceed 1 (I did not add any constraint). The blue plain horizontal line is the true value of the parameter, while the blue dotted horizontal line is the initial value of the parameter in the optimization algorithm (I started assuming that the mixture probability was around 0.2).
The graph below is based on the second optimization routine (with identical  starting values, and of course on the same generated samples),

(just to be honest, in many cases, it did not converge, so the loop stopped, and I had to run it again… so finally, my study is based on a bit less than 500 samples (times 15 since I considered several values for the mean of my second underlying distribution), with 200 generated observations from a mixture).
The graph below compares the two (empty circles are the first algorithm, while plain circles the second one),

On average, it is not so bad…. but the probability to be far away from the tru value is not small at all… except when the difference between the two means exceeds 3…
If I change starting values for the optimization algorithm (previously, I assumed that the mixture probability was 1/5, here I start from 1/2), we have the following graph

which look like the previous one, except for small differences between the two underlying distributions (just as if initial values had not impact on the optimization, but it might come from the fact that the surface is nice, and we are not trapped in regions of local minimum).
Thus, I am far from being an expert in optimization routines in R (see here for further information), but so far, it looks like R is not doing so bad… and the two algorithm perform similarly (maybe the first one being a bit closer to the trueparameter).

Hotline, modèles paramétriques et angoisse insoutenable

Après le billet sexe de l’été (ou presque, ici), je vais faire un  rapide un billet angoissant (pas angoisse insécuritaire comme dans les journaux, mais pire). En tout cas, moi ça m’angoisse, et il n’y a pas de raison que je sois le seul à angoisser quand j’appelle une hotline…
Afin d’illustrer mon billet, je me suis permis de faire une petite expérience: j’ai appelé 7 fois de suite une hotline (gratuite1) et j’ai noté le temps avant qu’un opérateur ne décroche2. J’ai obtenu les durées suivantes

1:21,10:07,2:21,8:58,0:11,0:34,1:13

Bref, à partir de ces données, on peut ajuster une loi pour modéliser le temps d’attente avant que l’opérateur ne décroche,

> T
[1]  1.3500000 10.1166667  2.3500000  8.9666667  0.1833333  0.5666667  1.2166667
> mean(T)
[1] 3.535714
> var(T)
[1] 17.40319
> library(MASS)
> (Exp <- fitdistr(T,"exponential"))
     rate   
  0.2828283 
 (0.1068990)
> (Gam <- fitdistr(T,"gamma"))
     shape       rate   
  0.7928017   0.2242283 
 (0.3653143) (0.1404783)

où j’ai retenu deux lois classiques sur les modèles de durées: la loi exponentielle et la loi gamma. Ces lois donnent des modèles assez comparable en regardant les fonctions de répartition.

> u=seq(0,12,.1)
> plot(ecdf(T))
> lines(u,pexp(u,rate=Exp$estimate),col="blue",lwd=3)
> lines(u,pgamma(u,shape=Gam$estimate[1],rate=Gam$estimate[2]),col="red",lwd=3)

Toutefois,ces modèles donnes des choses assez différentes si l’on regarde les espérances de vie résiduelles, i.e. les mean remaining lifetime en anglais.

> t2=as.POSIXlt(naissance$date)> F=function(x){1-pgamma(x,shape=Gam$estimate[1],rate=Gam$estimate[2])}
> e=function(x){
+ integrate(F,x,Inf)$value/F(x)
+ }
> v=rep(NA,length(u))
> for(i in 1:length(u)){
+ v[i]=e(u[i])
+ }

avec la loi exponentielle en bleu, et la loi gamma en rouge. La loi exponentielle est bien connue pour son absence de mémoire: peu importe le temps déjà attendu, l’espérance d’attente reste constante: ici, même si j’ai déjà attendu 3 minutes, ou 8 minutes, mon espérance d’attente est constante, à savoir 3 minutes et demi. On parle de constant mean remaining lifetime. En revanche, pour une loi gamma dont le shape paramètre est inférieur à 1, on a une loi dite increasing mean remaining lifetime (IMRL, dont la théorie peut se retrouver ici ou pour le chapitre 2 du livre de Maxim Finkelstein, en ligne )
En fait, si on regarde les évolutions des espérances de vie résiduelles, on a la loi exponentielle en bleu, constante, et la loi gamma, avec une fonction monotone: si le paramètre shape est inférieur à un (en rouge), on a une loi IMRL alors que s’il excède 1 (en vert), on a une loi DMRL (decreasing mean remaining lifetime).

Bref, à chaque fois que j’appelle une hotline et que j’attends plus de 3 minutes, je me demande toujours: quelle est la stratégie à adopter ? Pour ma part, le plus souvent, j’ai l’angoisse de la loi IMRL, autrement dit, plus j’ai attendu, plus le temps qu’il me reste à attendre est grand ! Et je raccroche assez vite…

Moralité, en statistique paramétrique, le choix du modèle n’est jamais neutre, donc il vaut mieux connaître les propriétés de différentes lois avant de se lancer dans un modèle paramétrique…. Et j’espère, ami  lecteur, que je ne serais plus le seul à me demander si je suis tombé sur une hotline IMRL la prochaine fois que j’essaye de connaître le prix d’un transfert bancaire…
1 à l’intention des jeunes chercheurs en herbe qui souhaitent faire avancer la science, l’expérience ne marche que sur des hotline gratuites….
2 au début je posais une question quelconque, sur la facturation d’un service bancaire ou autre, mais je dois avouer que sur la fin, je raccrochais sans rien dire… oui, la science autorise ce genre de fantaisie….

Intervalle de confiance pour une proportion

Je voulais poster un petit billet puisque c’est la deuxième fois que l’on me pose une question similaire. Pour déterminer un intervalle de confiance pour une proportion, on connaît tous la formule

La question qui m’avait été posée est simple: que se passe-t-il si n est petit, et si – empiriquement – on a trouvé une proportion de 100% ? Dit autrement: “on fait 5 lancers de “pile” ou “face” et que on obtient 5 fois “pile”: que peut-on dire sur p, la probabilité de tomber sur “pile” ?”

On utilisant la formule précédente, on est un peu bloqué…

  • La première réponse pourrait être fréquentiste, et consiste à proposer une lecture duale de l’intervalle de confiance, et des abaques.

Le graphique de droite montre, en trait violet, l’intervalle de confiance approché, c’est à dire avec l’intervalle de confiance indiqué ci-dessus. La courbe bleue, encadrant la zone verte correspond à l’intervalle de confiance calculé à partir des quantiles d’une loi binomiale de paramètre la fréquence empirique observée. Ici, on aurait envie de dire que p a 95% de chance d’être supérieure à 47.8%.

En fait, cette “lecture duale” des abaques revient à construire – me semble-t-il – l’intervalle de confiance de Clopper-Pearson. L’intervalle de confiance a ici la forme suivante

Notons que c’est cette technique qui est utilisée sous R, via la fonction binom.test du package de base. Ici, on obtient

> binom.test(x, n, p = 0.5,conf.level = 0.95)
Exact binomial test
data: x and n
number of successes = 5, number of trials = 5, p-value = 0.0625
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.4781762 1.0000000
sample estimates:
probability of success
1

…. ce qui correspond exactement à l’intervalle de confiance calculé précédemment. On retrouve d’ailleurs cette fonction dans binconf i.e.

> binconf(x = 5, n = 5, alpha = .05, method = "exact")
PointEst Lower Upper
1 0.4781762 1

qui donne (c’est rassurant) toujours la même chose. On peut aussi programme l’intervalle de confiance proposé par Heldge Blacker, qui donne ici un intervalle de confiance de la forme [50%; 100%]. Pour aller plus loin, je pourrais renvoyer au papier de Lawrence Brown, Tony Cai et Anirban DasGupta, dans Statistical Science sur ce sujet.

  • Une seconde réponse pourrait être une approche bayésienne, qui d’ailleurs correspond exactement à l’expérience de Bayes (1763) ou Laplace (1786).

On suppose que p est une variable aléatoire. On se donne une information a priori, et on peut se demander ce que devient la loi a posteriori de p, sachant que 5 “piles” ont été observés. On suppose que les tirages sont indépendants, et donc la probabilité d’avoir k fois “piles” sur n tirages suit une loi binomiale, conditionnellement à p, i.e.

Le plus simple est de considérer comme loi a priori la loi conjuguée de la loi binomiale, à savoir une loi beta. Dans ce cas, si p suit a priori une loi B(a,b), alors la loi a posteriori de p, sachant que x “pile” ont été observés est une loi B(a+x,b+n-x).

expérience décrite

La courbe en rouge est la densité de la loi a priori (uniforme ou beta de paramètres 2 et 2), et la courbe bleue, la loi a posteriori, avec un intervalle de confiance à 95%. Sur l’exemple de droite, par exemple, on a 95% de chance que p soit compris entre 47,3% et 96,8%.