Analyse de mélanges gaussiens pour de petits échantillons : applications à la cinématique stellaire

The estimation of the parameters of a mixture distribution is particularly delicate, especially when the sample size is small. We present in this paper two competing methods to deal with this problem, SAEM and Bayesian sampling. SAEM appears as a simulated annealing version of EM, while Bayesian sampling proposes a Monte-Carlo approximation of the Bayes estimator. We consider a joint use of the two approaches : SAEM provides some initial values for Bayesian sampling and Bayesian sampling somehow selects the best solution among these values. The model we consider in the paper is suggested by an astrophysical problem. We are comparing the two algorithms on three samples of bidimensional star velocities, in order to exhibit a mixture of two normal distributions. Key-words : mixture, small samples, stochastic algorithm, EM algorithm, Bayesian sampling, Bayes estimator, stellar kinematics. 18 C. SOUBIRAN, G. CELEUX, J. DIEBOLT, C. ROBERT 1. Motivations de l’étude La distribution des étoiles dans la Galaxie est fonction de plusieurs paramètres : position par rapport au centre de la Galaxie, vitesse (résiduelle et systématique), composition chimique, âge, masse. Dans les modèles de représentation de notre galaxie, les étoiles sont séparées en 2 populations. La Population I, ou population du disque, comprend des étoiles de même métallicité que le Soleil, concentrées dans le plan galactique, avec des orbites peu excentriques autour du centre de la Galaxie. On trouve dans la Population II des étoiles vieilles, de faible métallicité, ayant une distribution spatiale pratiquement sphérique et des orbites très excentriques. Certains modèles prennent en compte l’existence d’une population ayant des propriétés intermédiaires, et que l’on nomme disque épais. La cinématique stellaire joue un rôle très important dans l’étude des populations de notre Galaxie et permet de développer des théories sur sa formation et son évolution. On considère le mouvement des étoiles dans un système de référence centré sur le Soleil. Les directions des 3 axes sont notés : U dans la direction du centre galactique, V dans la direction de la rotation galactique, W dans la direction perpendiculaire au plan galactique. L’étude de la vitesse spatiale des étoiles du disque a montré que chacune de ces 3 vitesses a une distribution très proche d’une loi gaussienne, dont la variance diffère selon le type d’étoiles. La cinématique des étoiles naines de type spectral A normal a été étudiée récemment par Gômez et al. (1990). Les étoiles A sont des étoiles jeunes, d’âge inférieur à 109 ans. Du fait de leurs caractéristiques physiques (composition chimique, température effective,...), de leur distribution spatiale et de leurs propriétés cinématiques, elles appartiennent à la Population I. L’histogramme de la vitesse U de ces étoiles suggérait la possibilité d’un mélange, et a incité les auteurs à utiliser l’algorithme SEM (voir Celeux et Diebolt (1986)), qui est une version stochastique de la méthode itérative EM de résolution des équations du maximum de vraisemblance, introduite par Dempster, Laird et Rubin (1977). Effectivement, SEM leur a permis de mettre en évidence l’existence de 2 ou 3 sous-populations, leur proportion étant fonction de l’âge moyen des étoiles considérées. La conclusion de cette étude était que chacune de ces 2 sous-populations est représentative d’une même génération d’étoiles qui seraient nées en même temps, dans une même "bouffée" de formation. Il existe, parmi les étoiles A, des étoiles qui présentent des anomalies spectrales. Ces étoiles, bien que de composition chimique de surface différente, admettent une distribution spatiale identique à celle des étoiles de type spectral A normal, et devraient présenter le même comportement cinématique car elles ont le même âge. Nous disposons de 2 échantillons de taille N de ces étoiles particulières : AMI (N=51) défini par Gômez et al. (1981), 19 ANALYSE DE MÉLANGES GA USSIENS Ap4-Ap5 (N=38) défini par Grenier et al. (1981). A caase du petit nombre de données, les algorithmes EM et SEM n’ont pas permis de séparer de manière évidente d’éventuelles sous-populations gaussiennes dans ces 2 échantillons (voir Soubiran (1988)). Par contre, les 2 algorithmes (EB et SAEM) que nous présentons ici sont particulièrement adaptés à l’étude des petits échantillons. Dans un premier temps, nous avons testé la fiabilité de ces 2 algorithmes sur un échantillon d’étoiles A normales A2V (N=97). Cet échantillon a déjà été analysé par des techniques graphiques, ainsi que par les algorithmes EM et SEM (voir Bougeard et al. (1989), Soubiran et al. (1989), Gômez et al. (1990)), et a toujours conduit à des résultats très stables car il a 2 composantes bien séparées. Nous avons ensuite analysé la cinéinatique de ces étoiles A particulières à l’aide de l’algorithme SAEM. La méthode bayésienne nous a alors permis, en soumettant les données à un point de vue totalement différent, de confirmer ou d’infirmer les solutions proposées par SAEM. Dans cette étude, nous avons travaillé sur les vitesses U et V, qui sont les plus discriminantes (voir Gómez et al. (1990)). Mais il est possible d’adapter ces méthodes au cas tri-dimensionnel et à des distributions autres que les lois gaussiennes. 2. L’algorithme SAEM L’estimation par maximum de vraisemblance constitue l’approche la plus répandue pour l’identification de paramètres d’un mélange. Cette approche nécessite l’emploi d’algorithmes itératifs dont le plus utilisé est l’algorithme EM (voir Dempster et al. (1977) et Redner et Walker (1984)). Cet algorithme conduit dans de nombreux cas à des estimations satisfaisantes des paramètres. Néanmoins, il peut s’avérer décevant (lenteur excessive, dépendance très forte de la position initiale, convergence vers un col de la vraisemblance, etc.), particulièrement dans les situations délicates (composantes du mélange imbriquées, proportions des composantes déséquilibrées, etc.). Une modification stochastique de l’algorithme EM, l’algorithme SEM (voir, par exemple, Celeux et Diebolt (1986)), permet en général de dépasser les limitations précédentes. Cet algorithme comporte une étape stochastique (étape ’S’) qui vient s’intercaler entre l’étape ’E’ de calcul des probabilités a posteriori d’appartenance des points aux composantes du mélange et l’étape ’M’ de maximisation de la vraisemblance conditionnelle. L’étape ’S’ est fondée sur un principe d’affectation aléatoire qui consiste à remplacer les provenances (inconnues) des points de l’échantillon de l’une des composantes du mélange par des affectations aléatoires des observations à ces composantes, suivant les probabilités a posteriori calculées à l’étape ’E’. L’étape ’M’ se construit alors sur la base de ces affectations aléatoires. Les perturbations ainsi introduites par SEM lui permettent, en général, d’éviter les écueils de EM. 20 C. SOUBIRAN, G. CELEUX, J. DIEBOLT, C. ROBERT Les défauts de l’algorithme EM sont particulièrement génants pour les petits échantillons. En effet, pour N petit, la fonction de vraisemblance est souvent entachée de nombreux maxima locaux, instables et sans intérêt statistique, mais qui accentuent la dépendance de EM vis-à-vis de sa position initiale et rendent difficile l’utilisation de cet algorithme. De plus, l’algorithme SEM ne pallie pas les limitations de EM dans ces situations car les perturbations aléatoires de l’étape ’S’ prennent alors trop d’importance, et l’algorithme SEM risque alors de faire disparaître à tort une ou plusieurs des composantes du mélange analysé. L’algorithme SAEM (Celeux et Diebolt, 1990) a été conçu précisément pour conserver les qualités de SEM tout en évitant les défauts de EM, même pour de petits échantillons. Cet algorithme occupe en fait une place intermédiaire entre les deux méthodes précédentes. Comme SEM, il utilise un principe d’affectation aléatoire pour produire des perturbations dans la suite des estimations des paramètres, mais, grâce à une suite de nombres réels positifs qui décroît vers 0 quand n tend vers l’infini, il réduit l’intensité des perturbations aléatoires au fil des itérations et est ainsi plus fiable pour traiter de petits ensembles de données. De plus, comme EM, il propose à la convergence un estimateur ponctuel des paramètres. Il est par conséquent d’une plus grande simplicité d’utilisation que SEM dont l’estimateur converge seulement en distribution. La suite (03B3n) joue en fait le même rôle que la température dans les algorithmes de type recuit simulé (voir, par exemple, van Laarhoven (1988)), d’où le nom SAEM : Simulated Annealing EM. Nous décrivons maintenant l’algorithme SAEM en nous limitant au cadre qui nous intéresse ici, à savoir l’identification d’un mélange gaussien bidimensionnel. On dispose d’un échantillon Xl, ..., xN d’une variable aléatoire à valeurs dans R2, de densité avec 0 pi 1 1,p2 = 1 pl, et ~(.|03B8,03A3) désignant la densité d’une loi normale de moyenne 0 et de matrice de variance E. Il s’agit d’estimer le paramètre 0 (pl, 01, El, 03B82,03A32). Une itération de l’algorithme SAEM, qui à 0(’) associe 0(’+’), se décompose ainsi : Etape E : pour j = 1, 2, calcul des probabilités a posteriori d’appartenance des points xi à chaque composante du mélange (i == 1, ..., N). Etape S : tirage pour chaque xi de la variable aléatoire zn1(xi) suivant une loi de Bernoulli de paramètre tn1(xi). On pose zn2(xi) = 1 zn1(xi). Si zn1(xi) = 1, Xi est affecté à la première composante, sinon il est affecté à la seconde composante. Si l’une des deux classes ainsi obtenues a un cardinal plus petit que 2, l’algorithme s’arrête. Sinon, Etape A : calcul des quantités 21 ANALYSE DE MÉLANGES GA USSIENS pour i = 1,..., N et j = 1, 2. Etape M : calcul de ~(n+1), estimateur du maximum de vraisemblance de ~, les rj (Xi) étant considérés artificiellement comme les probabilités a posteriori d’appartenance des points xi aux composantes du mélange. Pour 3 1, 2, on en déduit les actualisations où at désigne le transposé