Análisis cluster y clasificación basados en modelos

Facultad de Ciencias M´ aster en matem´ aticas TRABAJO DE FIN DE M´ ASTER: AN´ ALISIS CLUSTER Y CLASIFICACI´ ON BASADOS EN MODELOS AUTORA: Carla Perucha Jurjo TUTORES: Carlos Matr´an Bea, Luis ´ Angel Garc´ıa Escudero Julio 2023 ´ Indice 1. Introducci´on 2 2. Presentaci´on del an´alisis cl´uster basado en modelos 4 2.1. Planteamiento intuitivo del modelo finito de mixturas . . . . . . . . . . . . . . . . . 4 2.2. Formulaci´on del modelo finito de mixturas . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.1. Estimaci´on de los par´ametros del modelo . . . . . . . . . . . . . . . . . . . . 9 2.2.2. El modelo de mezclas en el An´alisis Cluster . . . . . . . . . . . . . . . . . . . 11 2.3. Aplicaciones del modelo finito de mixturas . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3.1. Conjuntos de datos apropiados para el modelo de mezclas . . . . . . . . . . . 13 2.3.2. Tratamiento de datos desequilibrados . . . . . . . . . . . . . . . . . . . . . . 15 2.3.3. Detecci´on de valores at´ıpicos . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.4. Superposici´on entre los grupos . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3. Estudio te´orico del modelo de mixturas 22 3.1. Planteamiento del modelo general de mixturas . . . . . . . . . . . . . . . . . . . . . 22 3.2. Propiedades matem´aticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2.1. Existencia de soluci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2.2. Identificabilidad .................................. 30 3.3. Consistencia......................................... 34 4. Aspectos Metodol´ogicos 40 4.1. Restricciones geom´etricas en el modelo de mezclas . . . . . . . . . . . . . . . . . . . 40 4.2. Elecci´on del N´umero de Componentes y el modelo de Clustering . . . . . . . . . . . 42 4.2.1. Componentes no gaussianas y fusi´on de clusters . . . . . . . . . . . . . . . . . 46 4.3. ElalgoritmoEM ...................................... 49 4.3.1. Formulaci´on del algoritmo EM . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.3.2. El algoritmo EM en el modelo de mezclas gaussiano . . . . . . . . . . . . . . 53 4.3.3. Convergencia del Algoritmo EM . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.3.4. Inicializaci´on del Algoritmo EM . . . . . . . . . . . . . . . . . . . . . . . . . . 56 5. Ap´endice 60 5.1. Clases de Glivenko Cantelli y Vapnik-Chervonenkis . . . . . . . . . . . . . . . . . . . 60 5.2. Algunos resultados auxiliares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 1 1. Introducci´on En el contexto actual, donde la disponibilidad de datos masivos y la complejidad de las aplicaciones en diversos campos han aumentado exponencialmente, el An´alisis Cluster y la Clasificaci´on se han convertido en herramientas fundamentales para descubrir estructuras subyacentes y organizar grandes conjuntos de datos de manera significativa. Estas t´ecnicas multivariantes se encargan del estudio de m´etodos y algoritmos ideados para agrupar objetos en conjuntos de acuerdo con sus caracter´ısticas internas, a la par que asignar nuevos individuos a estos grupos en base a caracter´ısticas compartidas. Una gran cantidad de los enfoques del Aprendizaje Autom´atico se sustentan en t´ecnicas heur´ısticas para particionar los datos, pudiendo dar la sensaci´on de que las similitudes intr´ınsecas entre individuos no se basan en un modelo probabil´ıstico subyacente. Sin embargo, las agrupaciones que surgen en un conjunto de datos dependen fuertemente del m´etodo de agrupamiento elegido, siendo su actuaci´on muy dependiente de la morfolog´ıa de las observaciones. El An´alisis Cluster y la Clasificaci´on basados en modelos incorporan modelos estad´ısticos de los datos haciendo posible modelar distribuciones subyacentes complejas y capturar diferentes formas de agrupamiento. En este documento presentaremos los modelos finitos de mezclas, incidiendo especialmente en la familia de modelos en los que las componentes de mixtura consideradas son gaussianas. Por su naturaleza, estos modelos admiten un estudio matem´atico que permite demostrar propiedades como la existencia o la consistencia de soluciones, realizar comparaciones entre ellos, resaltar sus bondades y prever su comportamiento y sus inconvenientes o sus limitaciones. En las ´ultimas dos d´ecadas, se han producido progresos significativos en el problema de ajustar modelos de mezclas a conjuntos de datos. Muchos de los avances m´as destacados en este campo se han producido hace apenas 20 a˜nos, en particular la aplicaci´on del algoritmo de Expectation- Maximization (EM) para el c´alculo del estimador m´aximo veros´ımil revolucion´o el modelado a trav´es de modelos de mezclas y catapult´o su uso en diversas aplicaciones. En este trabajo hemos procurado introducir el modelo de mezclas junto con sus m´as recientes avances y retos metodol´ogicos respecto a la selecci´on de modelos, estimaci´on del n´umero de componentes e inicializaci´on del algoritmo EM, planteados en [2]. Adem´as, hemos completado esta exposici´on con el estudio de algunos de los resultados matem´aticos m´as interesantes ya que a menudo no encontramos resultados te´oricos en los manuales sobre el modelo de mezclas. Por ´ultimo, destacar que en las matem´aticas, en especial en el An´alisis Cluster, y en general en Clasificaci´on, una visualizaci´on por medio de ejemplos gr´aficos siempre resulta de gran ayuda para conseguir comprender mejor qu´e est´a sucediendo. Para ello, hemos empleando el lenguaje de programaci´on R, ampliamente utilizado en Estad´ıstica, donde existen m´ultiples paquetes que implementan algoritmos de clustering y funciones para visualizar sus resultados. En este documento se emplean los siguientes: rgl: Proporciona funciones para visualizaciones interactivas en 3D, ya sea a la hora de proporcionar gr´aficos (plot3d()) como funciones para contruir representaciones de objetos geom´etricos 2 tridimensionales (ellipse3d()). Los resultados se visualizan a trav´es de una ventana denominada RGL device 5 que permite efectuar rotaciones para poder explorar las diferentes perspectivas de los gr´aficos generados. mclust: Permite el ajuste de modelos de mezclas gaussianos por medio del algoritmo EM para el An´alisis Cluster basado en modelos, clasificaci´on y estimaci´on de funciones de densidad, incluyendo la regularizaci´on Bayesiana, inferencia basada en resampling y reducci´on de dimensiones para una mejor visualizaci´on. tclust: Incluye funciones de recorte propias del An´alisis Cluster robusto. La funci´on tclust() busca k(o menos) clusters con diferentes matrices de covarianza en un conjunto de datos dado. Se pueden especificar dos par´ametros muy interesantes: restr.fact, que restringe el ratio entre los autovalores de las matrices de covarianza y alpha, que permite un porcentaje de observaciones a recortar o eliminar a la hora de ajustar el modelo con el fin de hacer la estimaci´on m´as robusta. mixtools: Permite el an´alisis de los modelos finitos de mixtura para varias familias de densidades param´etricas, entre ellas la gaussiana. Por medio de funciones como normalmixEM() o mvnormalmixEM(), se puede estudiar el comportamiento del algoritmo EM al ajustar el MLE (como por ejemplo el valor de la funci´on de log-verosimilitud en cada iteraci´on, n´umero de iteraciones...) MASS: Recoge numerosos conjuntos de datos y las distribuciones de probabilidad m´as utilizadas en Estad´ıstica. En concreto, con la funci´on mvrnorm() se pueden generar muestras aleatorias de una normal multivariante especificando el vector de medias y la matriz de covarianzas. Rmixmod: Interfaz del software “MIXMOD” para clasificaci´on supervisada, semi-supervisada y no supervisada por medio del modelo de mixturas. Contiene la funci´on mixmodCluster() que ofrece la posibilidad de hallar el iterante inicial del EM en el con el algoritmo smallEM. Tambi´en est´a presente mixmodGaussianModel(), donde podemos introducir una lista de modelos seg´un su clasificaci´on VSO para el ajuste de los datos con el fin de compararlos entre s´ı. 3 2. Presentaci´on del an´alisis cl´uster basado en modelos Entre las numerosas t´ecnicas que existen para abordar el problema de discriminar individuos a partir de la medici´on de algunas de sus caracter´ısticas, el modelo finito de mixturas resulta una herramienta muy poderosa para identificar y modelar estas subpoblaciones. En esta secci´on presentaremos unos primeros ejemplos donde se ver´a patente el inter´es del m´etodo y daremos unas primeras definiciones relativas al problema de ajuste de un modelo de mixturas. Por otro lado, expondremos algunas de las situaciones en las que podr´ıa ser interesante utilizar esta metodolog´ıa. Como ya hemos comentado, los ejemplos expuestos a lo largo de este trabajo ser´an procesados con el lenguaje de programaci´on R, de uso extendido en el ´ambito de la Estad´ıstica. 2.1. Planteamiento intuitivo del modelo finito de mixturas Supongamos que contamos con un conjunto de datos constituido por 600 observaciones de individuos a los que hemos medido valores en las variables V1 y V2. Mostramos una representaci´on bivariante de los datos a continuaci´on: Figura 1: Representaci´on gr´afica de las observaciones en el plano para el primer conjunto de datos Las observaciones de este conjunto de datos est´an distribuidas en el plano formando una figura elipsoidal. Adem´as, existe una mayor densidad de puntos a medida que nos acercamos al centro de la nube de puntos. Construimos dos histogramas (uno para cada variable), donde dividimos el rango de valores que toman las observaciones en intervalos y mostramos la frecuencia con la que los datos caen en cada intervalo. 4 Figura 2: Histogramas para ambas variables Observamos que ambos histogramas tienen una forma bastante sim´etrica y acampanada que, junto con la informaci´on que recog´ıamos en el primer gr´afico, nos hace pensar que tal vez una distribuci´on normal bivariante con los par´ametros adecuados sea capaz de describir la variabilidad de las observaciones. Esto puede ser interesante ya que la distribuci´on gaussiana es ampliamente utilizada en Estad´ıstica, por lo que el modelo resultar´ıa simple y familiar, a la par que nos permitir´ıa utilizar propiedades conocidas de esta distribuci´on relativas a la inferencia y a la predicci´on. Un posible ajuste de los datos mediante una distribuci´on normal bivariante podr´ıa ser mediante una funci´on de densidad normal con vector de medias µ= (5,03 5,15) y matriz de covarianzas Σ = 6,93 −1,87 −1,87 6,78 estimadas por el m´etodo de m´axima verosimilitud. En el gr´afico siguiente, representamos el conjunto de datos se˜nalando con una “X” el centro de la distribuci´on y con una l´ınea verde la elipse de confianza del 95 % asociada a Σ : 5 Figura 3: Ajuste por medio de una distribuci´on Gaussiana Consideremos ahora otro conjunto de datos constituido por otras 600 observaciones bivariantes. Si nos fijamos en su representaci´on en el plano como puntos de R2, observamos que los individuos toman valores bastante diferentes en las variables V1 y V2. Adem´as, al dibujar la funci´on de densidad de los datos para cada variable vemos que claramente son bimodales (observamos dos cumbres, dos modas). Figura 4: Representaci´on bivariante de las observaciones del segundo conjunto de datos 6 En el gr´afico, se han dibujado por separado las dos densidades normales junto con la densidad de probabilidad total, que es una combinaci´on lineal de dos densidades normales ponderadas. Los puntos que aparecen en la parte inferior del gr´afico son una muestra de esta distribuci´on y est´an coloreados en funci´on de qu´e componente los ha generado. Observamos adem´as que la funci´on de densidad es m´as baja entre las dos jorobas, por lo que existe algo de separaci´on entre las dos componentes en cada una de las dos variables. Sin embargo, esta separaci´on no es total y existe incertidumbre a la hora de saber qu´e puntos provienen de cada una de las componentes. Parece razonable tratar de extender la idea anterior de resumir la distribuci´on y variabilidad de los datos al caso en el que contamos con Gdistribuciones en lugar de solo una. Este procedimiento se denomina modelo de mixturas y trata de representar cada subgrupo o componente del conjunto de datos mediante una distribuci´on de probabilidad, en el ejemplo que estamos planteando esta ser´ıa la distribuci´on Gaussiana. Si consideramos el siguiente conjunto de datos, por la colocaci´on que tienen los puntos en el plano, parece coherente pensar que existen tres subploblaciones, por lo que tratamos de describir a los individuos por medio de tres distribuciones Gaussianas. El objetivo ser´a elegir los par´ametros de las funciones de densidad y las proporciones en las que aparecen con el fin de lograr la mejor representaci´on de los datos. Figura 5: Tres componentes para ajustar los datos En las situaciones que hemos expuesto, las agrupaciones podr´ıan parecer triviales y factibles de detectar con casi cualquier procedimiento de An´alisis cluster ya que existen diferencias bastante evidentes entre los individuos de distintas poblaciones. Sin embargo, podemos encontrarnos con configuraciones mucho m´as complejas donde no es evidente una separaci´on y donde las observaciones se presentan como una mezcla indistinguible sin el conocimiento de los grupos. Un ejemplo de esto 7 podr´ıa ser el siguiente conjunto con 600 observaciones donde, por sus caracter´ısticas, el modelo de mixturas puede desentra˜nar su estructura latente. Figura 6: Un conjunto de datos donde los grupos no son tan triviales 2.2. Formulaci´on del modelo finito de mixturas El objetivo del modelo de mezclas es estimar los par´ametros de cada componente y las proporciones con las que se presentan para explicar de la mejor manera posible la variablidad de los datos. Concretamos un poco m´as a qu´e nos referimos con todo esto en la siguiente definici´on: Definici´on 2.1. Sea A={x1, .., xn}un conjunto con nobservaciones multivariantes pertenecientes aRd, de tal manera que xi= (xi,1, .., xi,d). Sea G un entero positivo y πjpesos tales que πj≥0para j= 1, .., G con PG j=1 πj= 1. Sea F={f(·;θ)/θ ∈Θ}una familia de funciones de densidad con vector de par´ametros θen el correspondiente espacio de par´ametros Θ. Un modelo finito de mezclas es un modelo estad´ıstico que representa la funci´on de densidad de una observaci´on xi,f(xi), como una combinaci´on lineal ponderada de Gfunciones de densidad fj(·;θj)∈F: f(xi) = G X j=1 πjf(xi;θj) (1) Las funciones de densidad f(·;θj) = fjse denominan componentes de la mixtura y los valores πj proporciones de mezcla o pesos. Observaci´on 2.2. Dado que las funciones f1(·;θ1), .., fG(·;θG)son densidades, es claro que la expresi´on (1) define una funci´on de densidad. El vector que recoge los par´ametros desconocidos a estimar es Ψ = (π1, .., πG−1, θ1, .., θG), donde hemos omitido πGpor ser redundante ya que las proporciones πjsuman uno. Denotamos por Θc 8 Figura 9: Diferentes ajustes de densidades multimodales En la figura, observamos varios ejemplos de conjuntos de datos unidimensionales con distribuciones multimodales. En el primer caso, se ven claramente dos protuberancias en la funci´on de densidad, con un valle entre medias. Adem´as, parece que el tama˜no de las mismas es id´entico, por lo que si utiliz´asemos un modelo de mezclas gaussiano, ser´ıa adecuado ajustar un modelo de dos mezclas con la misma varianza. En el segundo caso, uno de los picos es m´as alto que el otro, debido a varianzas diferentes en las subpoblaciones. A continuaci´on, en el tercer caso, tenemos un conjunto de puntos que ha sido generado por tres distribuciones normales, a pesar de que en la funci´on de densidad resultante al combinarse pudiera parecer que solo existen dos modos. Para finalizar, observamos una funci´on de densidad con tres modos de diferentes tama˜nos, por lo que parecer´ıa razonable utilizar una mixtura con tres componentes de distinta varianza. 2.3.2. Tratamiento de datos desequilibrados Dado que para ajustar el modelo de mezclas tenemos que estimar las proporciones de las componentes, este procedimiento puede resultar interesante en conjuntos de datos con clases desequilibradas. Puede ser ´util a la hora de compensar la falta de observaciones de un grupo en el contexto de la clasificaci´on o bien para capturar la asimetr´ıa de una funci´on de densidad, model´andola a trav´es de varias componentes con diferentes pesos. Consideremos el siguiente conjunto de datos desequilibrados cuya separaci´on es relativamente buena en el que contamos con 8 grupos de observaciones bidimensionales, tres de los cuales tienen 2000 individuos mientras que los cinco restantes tienen 100. La flexibilidad a la hora de ajustar los 15 pesos π1, .., π8consigue evitar la situaci´on desfavorable en la que el m´etodo de mezclas gaussiano trate de fragmentar uno de los grupos m´as densos. Observamos como la partici´on que se genera es bastante buena. Figura 10: Conjunto de datos desequilibrado Por otro lado, si tenemos inter´es en modelar un conjunto de datos cuya distribuci´on tiene asimetr´ıa, es posible capturar esta concentraci´on de la masa de probabilidad hacia alguna direcci´on por medio del modelo de mixturas. Figura 11: Modelado de funciones de densidad variando su sesgo En la imagen podemos observar c´omo utilizando dos componentes en la mezcla y otorgando mayor peso a una de ellas conseguimos capturar la asimetr´ıa de las funciones de densidad consideradas, algo que no hubiera sido posible utilizando una sola densidad gaussiana debido a su car´acter sim´etrico. 2.3.3. Detecci´on de valores at´ıpicos De manera general en el An´alisis Estad´ıstico, un outlier se define como un elemento cuyas caracter´ısticas no siguen el mismo patr´on que la mayor´ıa del resto de observaciones. Estos puntos, en el 16 contexto del modelo de mixturas, pueden ser observaciones at´ıpicas que est´an alejadas de los grupos definidos por las componentes o que est´an dispuestas entremedias de ellos, en particular cuando existe una buena separaci´on entre clusters. Suponen un problema ya que al distorsionar la distribuci´on real de los datos, influyen en la estimaci´on de par´ametros. Adem´as, los modelos estad´ısticos pueden ser sensibles a los valores extremos, por lo que incluso unos pocos outliers pueden tener un impacto desproporcionadamente grande en los resultados del modelo de mezclas. Es por lo tanto importante identificar y tratar adecuadamente las observaciones at´ıpicas para obtener resultados m´as precisos, robustos y confiables. El modelo de mezclas, debido a la flexibilidad que tiene a la hora de elegir las proporciones de las componentes, es una herramienta ´util en algunas ocasiones para detectar estas observaciones at´ıpicas. Consideramos el siguiente conjunto de datos donde se tienen dos grupos, de 60 y 120 observaciones cada uno, elipsoidales y f´acilmente distinguibles. Por otro lado, generamos 20 observaciones con distribuci´on normal pero con una matriz de covarianza de determinante mayor. Algunos de estos outliers son f´acilmente identificables, mientras que aquellos que se encuentran en las vecindades de los grupos resultan m´as dif´ıciles de detectar. Si ajustamos tres componentes gaussianas, el modelo de mixturas pr´acticamente es capaz de recuperar las dos agrupaciones iniciales con bastante ´exito, asignando a la tercera componente aquellos puntos que son considerados outlier. Figura 12: Detecci´on de observaciones at´ıpicas utilizando el modelo de mezclas. En ocasiones, el m´etodo no tiene capacidad suficiente para tratar con ´exito este tipo de observaciones. Existen multitud de enfoques para abordar esta situaci´on, pero nosotros presentaremos dos de ellos en este apartado. Componente uniforme para los valores at´ıpicos Una forma de enfrentarnos a los outliers o al ruido en el modelo de mixturas consiste en a˜nadir una componente adicional al modelo para representar la aportaci´on de estas observaciones at´ıpicas. Por ejemplo, se puede elegir una componente uniforme para tratar de reflejar esa dispersi´on equitativa de los valores at´ıpicos en la regi´on del espacio donde se sit´uan los grupos, 17 de tal manera que se localicen tanto en las afueras de los clusters como entremedias de ellos. De este modo, el modelo que hab´ıamos definido antes tendr´ıa ahora la siguiente forma: f(xi) = π0 V+ G X j=1 πjf(xi;θj),1≤i≤n donde Ves el volumen aproximado de la regi´on que ocupan los datos y π0es la proporci´on esperada de valores at´ıpicos. Podemos encontrar m´as informaci´on sobre este planteamiento en [2] y sobre la sensibilidad del algoritmo EM dependiendo de la inicializaci´on en el caso de existir outliers. Procedimientos de recorte Podemos encontrar estas ideas en art´ıculos como [4], donde se presenta una t´ecnica que consiste en eliminar parte de los datos de manera no arbitraria con el fin de hacer m´as robusto el procedimiento de mixturas, en concreto el caso en el que las funciones de densidad son gaussianas. Para ello, se consideran ciertas restricciones para asegurar la robustez de los par´ametros estimados y se construye la funci´on de verosimilitud recortada n X i=1 z(xi) logG X j=1 πjϕ(xi;µj,Σj) Donde z(·) es una funci´on de recorte que vale 0 si la observaci´on ha sido eliminada y 1 en el caso contrario. Previamente, se elige un n´umero α∈(0,1) denominado nivel de recorte que ser´a la proporci´on de datos que se desea recortar de la muestra {x1, .., xn}. Esa proporci´on de datos ser´a la que asumimos como observaciones espurias que no queremos incluir cuando estimemos los par´ametros del modelo de mixturas. Por lo tanto, se tiene que Pn i=1 z(xi)=[n(1 −α)], donde [·] denota el redondeo hacia arriba. Este art´ıculo proporciona adem´as una versi´on robusta del algoritmo EM que incluye un paso de recorte adicional. El criterio que se utiliza para decidir qu´e observaciones parecen ser espurias es calcular di= m´ax{π1ϕ1(yi), .., πGϕG(yi)}y y seleccionar la proporci´on αde puntos que toman los valores m´as peque˜nos de di. 18 Figura 13: Recorte de nivel α= 0,2 2.3.4. Superposici´on entre los grupos La existencia de superposici´on entre los grupos en un conjunto de datos supone un desaf´ıo para el modelo de mezclas ya que puede tener dificultades para distinguir las subpoblaciones. En la siguiente figura observamos la actuaci´on del modelo de mezclas gaussiano al tratar de distinguir los 15 grupos de la figura, cada uno con 300 observaciones. Dado que los clusters son elipsoidales y no existe una superposici´on muy notoria, con pocos puntos entre los clusters, la partici´on resulta bastante satisfactoria. Figura 14: Grupos bien separados Sin embargo, al considerar un conjunto de datos con mayor grado de superposici´on, la estimaci´on de par´ametros se vuelve m´as sensible a peque˜nas variaciones y adem´as el m´etodo fracasa a la hora de detectar los grupos. Observamos el siguiente ejemplo donde se ha realizado un ajuste por medio de un modelo de mezclas gaussiano a un conjunto de datos con tres subpoblaciones. En la figura de la derecha vemos una comparaci´on gr´afica entre los par´ametros reales utilizados para generar los datos y los par´ametros estimados por el modelo de mezclas gaussiano. 19 Figura 15: Elecci´on de par´ametros con superposici´on entre tres grupos A la hora de elegir los vectores de medias, la actuaci´on del modelo es bastante buena, mientras que el ajuste de las matrices de covarianza no resulta del todo exacto. Aun as´ı, el m´etodo consigue capturar relativamente bien la disposici´on y orientaci´on de las componentes. Figura 16: Diferencia entre la partici´on generada por el modelo y la agrupaci´on real 20 Sin embargo, la presencia de puntos con probabilidades a posteriori parecidas de pertenecer a varios grupos produce una gran incertidumbre a la hora de asignar las observaciones a un cluster o a otro. Esto genera una mayor variabilidad en los resultados al variar la muestra y dificulta la interpretaci´on de los grupos obtenidos, ya que en la agrupaci´on real los individuos est´an entremezclados en ciertas zonas de su representaci´on en el plano y el modelo en cambio realiza una separci´on total de los grupos de manera que no existe solapamiento. Lo mismo ocurre al considerar conjuntos incluso con mayor grado de superposici´on. En el gr´afico siguiente vemos se˜nalados con una equis los vectores de medias estimados por el modelo de mezclas gaussiano y con c´ırculos rojos los verdaderos centros. De nuevo, la elecci´on de los centros parece razonable, mientras que la actuaci´on a la hora de crear los clusters de nuevo es bastante deficiente, ya que las fronteras son variables y poligonales. Figura 17: Variabilidad y rigidez en las fronteras en el caso de superposici´on entre grupos Es importante tener en cuenta que la creaci´on de estas particiones tan fluctuantes en presencia de superposici´on entre grupos no siempre implica un fallo del modelo de mezclas. En algunos casos, puede ser un reflejo de la naturaleza de las observaciones: si existe un solapamiento real entre los grupos es posible que simplemente sus individuos sean dif´ıciles de distinguir. 21 3. Estudio te´orico del modelo de mixturas Hemos planteado la idea de configurar una descripci´on de un conjunto de npuntos de la mejor manera posible utilizando Gfunciones de densidad pertenecientes a una familia de distribuciones {f(·;θ), θ ∈Θ}. En esta secci´on, buscamos enfocar este problema desde una perspectiva m´as general, a trav´es de un estudio te´orico que incluya tanto conjuntos de datos como distribuciones de probabilidad. Para ello, generalizaremos de manera natural la expresi´on de la funci´on de verosimilitud utilizando la esperanza matem´atica. A continuaci´on, particularizando en la familia de distribuciones Gaussianas, hablaremos de algunas propiedades matem´aticas de este problema como la existencia, la identificabilidad o la consitencia. Esta familia de distribuciones es especialmente interesante ya que encontramos de manera frecuente el problema en el que los individuos provienen de diferentes poblaciones normales. Tanto para hablar de la existencia de un vector de par´ametros que maximice la expresi´on de la verosimilitud como para hablar de la consistencia del vector de estimadores muestrales es necesario imponer ciertas restricciones en el problema, algunas de ellas especialmente confeccionadas para las funciones de densidad normales, siendo esta otra raz´on por la que nos enfocamos en el caso de la familia de distribuciones Gaussianas. A lo largo de esta secci´on trataremos de ilustrar gr´aficamente estos resultados te´oricos a fin de facilitar la visualizaci´on y la comprensi´on de los resultados planteados. 3.1. Planteamiento del modelo general de mixturas Inicialmente hab´ıamos introducido el ajuste un modelo de Gfunciones de densidad a un conjunto de nobservaciones {x1, .., xn}. Supongamos que en lugar de tener un conjunto con ndatos contamos con una distribuci´on de probabilidad: Xser´a una variable aleatoria definida en un espacio probabil´ıstico (Ω, σ, P) y Pla ley de probabilidad que induce en R. Podemos reescribir la funci´on de verosimilitud para esta probabilidad P en lugar de para nobservaciones de Xpor medio de la esperanza matem´atica, que no es sino la generalizaci´on de la media a trav´es del concepto de integral respecto de una probabilidad: E(X) = ZR xdP(x) (3) Para variables aleatorias multidimensionales, su esperanza o valor esperado se define componente a componente, esto es, dado un vector aleatorio X= (X1, .., Xp):Ω−→ Rp, definimos su esperanza como E[(X1, .., Xp)] = [E(X1), .., E(Xp)] De este modo, para dicho vector aleatorio X, la funci´on de log-verosimilitud que trat´abamos de maximizar en la secci´on anterior se reescribir´ıa como 22 L(Ψ, P) = EP log   G X j=1 πjf(x;θj)  , f ∈F={f(·;θ) : θ∈Θ} Este planeamiento del problema mediante distribuciones de probabilidad generaliza el caso de tener un conjunto con npuntos. En el caso de tener x1, .., xnelementos de Rd, podemos considerar Pn, la probabilidad muestral que asigna probabilidad 1 na cada xiy utilizar que la esperanza no es sino el promedio de los npuntos. Estamos en disposici´on de plantear el problema de ajuste de un modelo de mezclas a la distribuci´on P: Definici´on 3.1. Sea X: (Ω, σ, P)−→ Rdun vector aleatorio y sea Pla probabilidad inducida por Xen Rd. Sea Gun entero positivo y π1, .., πG−1≥0pesos tales que PG j=1 πj= 1. Sea F= {f(·;θ)/θ ∈Θ}una familia de funciones de densidad donde Θdenota el espacio param´etrico y Ψ = (π1, .., πG, θ1, .., θG)el vector que recoge los par´ametros a estimar. El problema de ajustar a la variable Xun modelo de Gdistribuciones respecto a la familia Fconsiste en encontrar Ψ∗que maximice la expresi´on L(Ψ, P) = EP log   G X j=1 πjf(x;θj)  (4) 3.2. Propiedades matem´aticas Abordaremos en esta secci´on el estudio te´orico del modelo de mixturas para el caso particular en el que la familia de distribuciones sea la gaussiana: fjser´a la funci´on de densidad normal multivariante ϕj, parametrizada por su vector de medias µjy su matriz de covarianzas Σj, con la forma ϕ(xi;µj,Σj) = 1 (2π)n 2|Σj|1 2 exp {−1 2(xi−µj)TΣ−1 j(xi−µj)} Es decir, consideraremos el problema en el que F={ϕ(·;µ, σ):(µ, σ)∈Θ}, donde Θ = Rp×Mp×p. En primer lugar, plantearemos las hip´otesis necesarias para poder garantizar la existencia de una soluci´on al problema de maximizaci´on de la log-verosimilitud. Despu´es, definiremos el concepto de identificabilidad de un modelo y demostraremos que los modelos de mezclas gaussianos poseen esta propiedad. Por ´ultimo, estudiaremos la consistencia del m´etodo, que resulta interesante ya que a menudo el problema aparece en un ´ambito estad´ıstico en el que contamos con conjuntos de datos que podemos interpretar como una muestra. 3.2.1. Existencia de soluci´on A continuaci´on, abordaremos el problema de existencia de un vector de par´ametros que minimicen la expresi´on (4). La existencia de una soluci´on no est´a garantizada sin imponer ciertas restricciones sobre la medida de probabilidad Py la relaci´on entre los autovalores de las matrices de covarianzas. La funci´on de verosimilitud no est´a acotada en el caso en el que nos propusiesemos ajustar un modelo de Gmixturas a un conjunto de Gpuntos o a un conjunto de datos en el que su disposici´on se aproxime demasiado a un subespacio de dimensi´on menor, ya que en estos casos el determinante de la 23 matriz de covarianzas tender´ıa a cero y, en consecuencia, ϕjtender´ıa a infinito para esa componente en concreto. Si excluimos el caso en el que manejemos una probabilidad concentrada en Gpuntos e impedimos que las proporciones entre los tama˜nos de los autovalores de las matrices de convarianza sean arbitrarias, seremos capaces de asegurar la existencia y consistencia de soluciones. Los resultados que encontramos en esta secci´on son una adaptaci´on de los que aparecen en los art´ıculos de Garc´ıa Escudero et al. (2008, 2015) [6], [5], donde se trabaja con un problema m´as general denominado ajuste de mixturas Gaussianas recortado. En ´el, se presenta una t´ecnica que consiste en eliminar una proporci´on αde los datos de manera no arbitraria con el fin de hacer m´as robusto el procedimiento. En esta versi´on del problema, una de las restricciones que expondremos a continuaci´on para garantizar la existencia y consistencia de un vector de par´ametros que maximice la verosimilitud dejar´ıa de ser necesaria, por lo que cuando lleguemos a ese punto resaltaremos este hecho. Dada la matriz de covarianzas Σj, correspondiente a la j−´esima componente de la mixtura, dado que es sim´etrica y definida positiva, podemos encontrar matrices EjyDjde tal manera que Σj=EjDjET j Donde Djes la matriz diagonal que recoge los autovalores λj,1, .., λj,G de ΣjyEjes la matriz de cambio de base cuyas columnas son los autovectores unitarios vj,1, .., vjGasociados a los autovalores anteriores (al ser los autovectores ortonormales, ET jcoincide con E−1 j). Esta descomposici´on nos resultar´a ´util a la hora de probar ciertos resultados posteriores y nos ayudar´a a entender una de las condiciones que imponemos sobre los autovalores de las matrices de covarianzas. Dicho esto, presentamos a continuaci´on las restricciones necesarias para poder garantizar la existencia y la consistencia: 1. (PR) Restricciones sobre la distribuci´on P: La distribuci´on P no debe est´ar concentrada en Gpuntos. P debe de tener momento de segundo orden finito, esto es: EP(∥·∥2)<∞. En el caso de estar trabajando con recortes, esta hip´otesis no ser´ıa necesaria. 2. (ER) Restricciones sobre el ratio entre autovalores: Debe cumplirse, para una constante c≥1, que Mn/mn≤c Donde Mn= m´ax 1≤j≤Gm´ax 1≤l≤pλjl mn= m´ın 1≤j≤Gm´ın 1≤l≤pλjl Es decir, queremos tener controlado el ratio entre el autovalor m´as grande y el m´as peque˜no entre todas las matrices de covarianza. En el ´ambito de la b´usqueda de agrupaciones o de la clasificaci´on de individuos, esta restricci´on controla simultaneamente las diferencias entre grupos y cu´anto distan estos de ser esf´ericos. 24 F(x) = G X j=1 πjF(x;θj) Trabajaremos con esta expresi´on an´aloga para aprovechar las propiedades de las funciones de distribuci´on a la hora de dar definiciones o planear resultados. Consideremos F={F(·;θ)/θ ∈Θ}una familia de funciones de distribuci´on parm´etricas y definamos Mla familia de todas las mixturas finitas de la clase F, que coincide con su envoltura convexa: M={F(·; Ψ) = G X j=1 πjF(·;θj) : πj≥0, G X j=1 πj= 1, θj∈Θ, G ∈N,1≤j≤G} Esta familia es la que consideramos cuando tratamos de ajustar un modelo de mezclas a un conjunto de datos o a una distribuci´on de probabilidad. En este caso, decimos que F(·; Ψ) es identificable si es invariante bajo las G! permutaciones de las etiquetas de Ψ: Definici´on 3.6. Sea Mla familia param´etrica de distribuciones propia del modelo de mezclas. Sean F(·; Ψ) = PG j=1 πjF(·;θj)yF′(·; Ψ′) = PH j=1 π′ jF(·;θ′ j)dos elementos de M. Decimos que Mes identificable para Ψsi F(·; Ψ) ≡F′(·; Ψ′)si y solo si G=Hy podemos permutar las etiquetas de las componentes de mixtura de tal manera que πj=π′ jyF(·;θ′ j) = F(·;θj)para 1≤j≤G. Figura 18: Diferentes permutaciones de un modelo identificable Para evitar la falta de identificabilidad de Ψ debido al intercambio de etiquetas entre las componentes, se suelen imponer restricciones sobre Ψ como π1≤π2≤... ≤πG. Se tiene la siguiente caracterizaci´on sobre las familias de funciones de distribuci´on param´etricas identificables. Teorema 3.7. Una condici´on necesaria y suficiente para que la clase Msea identificable es que F sea un conjunto linealmente independiente sobre el cuerpo de los n´umeros reales. 31 Demostraci´on Denotamos por Fj=F(·;θj) y por ⟨A⟩el conjunto de las combinaciones lineales finitas de elementos de Asobre el cuerpo de los reales ⟨A⟩=nk X i=1 αiai:k∈N, ai∈A, αi∈Rpara 1≤i≤ko Supongamos que Fno es linealmente independiente, esto es, existen α1, .., αk∈Rno todos nulos tales que Pk j=1 αjFj= 0. Estos αjpueden ser estrictamente negativos, cero o estrictamente positivos. Si hay M∈ {1, .., k}coeficientes estrictamente negativos, reordenemos los elementos de tal manera que coloquemos al principio: aj<0 si 1 ≤j≤M. As´ı, se tendr´ıa M X j=1 |αj|Fj= k X j=M+1 |αj|Fj Dado que los elementos de Fson funciones de distribuci´on, sabemos que l´ımx→∞ F(x) = 1, por lo que haciendo tender xa infinito en ambos lados de la igualdad, llegamos a que M X j=1 |αj|= k X j=M+1 |αj|=β > 0 ya que al menos uno de los αjes no nulo. Si tomamos γj=|αj| β, entonces PM j=1|γ|jFj=Pk j=M+1|γ|jFj ser´ıan dos representaciones distintas de la misma mixtura finita, por lo que Mno ser´ıa identificable. Por otro lado, supongamos que Fes linealmente independiente. Eso quiere decir que es una base de ⟨F⟩, por lo que dos representaciones distintas de la misma mixtura supondr´ıan una contradicci´on de las propiedades de representaci´on ´unica que tienen las bases, teniendo en cuenta que Mser´ıa un subespacio de ⟨F⟩. □ Teorema 3.8. La familia Fde funciones de distribuci´on gaussianas n−variantes generan mixturas finitas identificables. Supongamos que Fno es identificable y lleguemos a una contradicci´on. Por el teorema visto anteriormente, es equivalente suponer que la clase de funciones de distribuci´on normales multivariantes de dimensi´on nes un conjunto linealmente dependiente sobre R. Esto implica que existen α1, .., αkno todos nulos tal que Pk j=1 αjFj= 0. Reescribimos como hac´ıamos en el resultado anterior, colocando en los primeros ´ındices los Melementos cuyo coeficiente es estrictamente negativo para tener M X j=1 αjFj= k X j=M+1 αjFj(5) Es decir, en los dos sumatorios considerados, los coeficientes son no negativos. Utilizaremos dos resultados relativos a la funci´on caracter´ıstica de una variable aletoria y la relaci´on entre la funci´on de distribuci´on de una variable y su funci´on generadora de momentos. El 32 Teorema de Unicidad de la funci´on caracter´ıstica asegura que distribuciones de probabilidad distintas tienen diferentes funciones caracter´ısticas (ver Teorema 26.2 en Billingsley [1]). Como consecuencia, teniendo en cuenta que la integral de una funci´on respecto de una mixtura de probabilidades es igual a la mixtura de las integrales, dadas F1, F2, ... distribuciones de probabilidad con funciones caracter´ısticas φ0, φ1, ... yβi≥0 con Pr i=1 βi= 1, la mixtura U=Pr i=1 βiFies una distribuci´on de probabilidad con funci´on caracter´ıstica φ=Pr i=1 βiφi. Volviendo a nuestra expresi´on, consideremos PM j=1 αj=γ > 0. Si dividimos entre γa ambos lados de la expresi´on (5), obtendr´ıamos M X j=1 αj γFj= k X j=M+1 αj γFj Es claro que los coeficientes de la expresi´on de la izquierda son pesos (mayores o iguales que cero cuya suma es uno), por lo que U=PM j=1 αj γFjser´ıa una distribuci´on de probabilidad. Adem´as, como las Fjson funciones de distribuci´on, se tiene que l´ım x→∞ k X j=M+1 αj γFj(x) = k X j=M+1 αj γ= M X j=1 αj γ= l´ım x→∞ M X j=1 αj γFj(x) Pero por c´omo hemos definido γ,PM j=1 αj γ= 1, por lo que los coeficientes αM+1 γ, .., αk γson pesos y se tiene que V=Pk j=M+1 αj γFjes una distribuci´on de probabilidad tambi´en. Dado que UyVson distribuciones de probabilidad iguales, necesariamente sus funciones caracter´ısticas tambi´en lo son. Si φjdenota la funci´on caracter´ıstica de Fj, se tendr´ıa entonces M X j=1 αj γφj= k X j=M+1 αj γφj⇒ M X j=1 αjφj= k X j=M+1 αjφj⇒ k X j=1 αjφj= 0 Denotando por µjy Σjel vector de medias y la matriz de coviarianzas de Fjrespectivamente, dado que Fjes la funci´on de distribuci´on de una normal N(µj,Σj), su funci´on caracter´ıstica es φj= exp{1 2xTΣjx+xTµj}. Por lo tanto k X j=1 αjexp {1 2tTΣjt+tTµj}= 0 Si tomamos t=cu donde ces un escalar y ues un vector, obtenemos k X j=1 αjexp{c2 2uTΣju+uTµj}= 0 (6) Si todas las matrices de covarianza son iguales, los vectores de medias deben ser todos diferentes entre s´ı para que exista diferencia de par´ametros entre los ´ındices. Por lo tanto, salvo para un n´umero finito de hiperplanos, los pares de n´umeros reales (uTΣju, uTµj) son distintos para cada j∈ {1, .., k}. De otro modo, supongamos que Σ1, .., ΣNson las ´unicas matrices de covarianza distintas entre 33 Σ1, .., Σk. Entonces, para un vector uque no pertenezca a un n´umero finito de c´onicas, los n´umeros reales utΣjuson distintos para 1 ≤j≤N. Para los ´ındices N+ 1, .., k, los µjson diferentes entre s´ı, por lo que para los vectores que no pertenezcan a un n´umero finito de hiperplanos y c´onicas, los pares (uTΣju, uTµj) son distintos dos a dos, 1 ≤j≤k. Si elegimos un ucon estas caracter´ısticas, lo que estar´ıamos diciendo en (6) es que la clase de las mixturas finitas de normales univariantes no es identificable, lo cual contradice los resultados demostrados por Teicher [11]. 3.3. Consistencia Una vez planteado el problema de existencia, buscamos evaluar si el modelo finito de mixturas produce resultados estables y consistentes a medida que se aumenta el tama˜no de la muestra. El marco de trabajo caracter´ıstico de la Estad´ıstica Matem´atica es el escenario en el que tenemos un conjunto de nindividuos {x1, .., xn}que toman medidas en pvariables y se considera que estas observaciones se obtienen como muestra en un modelo probabil´ıstico donde el mecanismo generador de los datos Pes desconocido. A la hora de ajustar un modelo de mixturas a una variable aleatoria Xcon distribuci´on P, los vectores de medias y las matrices de covarianza que hallamos son concretos y no dependen de la aleatoriedad. Sin embargo, estimar estos mismos par´ametros para una muestra dada depende del ωescogido, por lo que obtenemos valores distintos al generar muestras nuevas. Buscamos por lo tanto estudiar la consistencia estad´ıstica del modelo, es decir, analizar la capacidad del m´etodo para obtener resultados consistentes y estables a medida que tama˜no de la muestra aumenta indefinidamente. En t´erminos pr´acticos, nos preguntamos si los vectores de medias y matrices de covarianza que hallamos para una muestra se acercan a aquellos de la distribuci´on de probabilidad te´orica. Figura 19: Vectores de medias y matrices de covarianza seg´un aumenta el tama˜no de la muestra En esta figura podemos ver por ejemplo el resultado de estimar los vectores de medias y las ma- 34 trices de covarianza en el modelo de mixturas cuando generamos una muestra de tama˜no 60, 150 y 1500 obsevaciones respectivamente de una combinaci´on balanceada de tres normales multivariantes X1∼N((0,0),1 0 0 1 ), X2∼N((6,5),4−2 −2 4 ) y X3∼N((12,5),1 0 0 10 ). Podemos observar a grandes rasgos que, a medida que aumenta el tama˜no de la muestra considerada, se estabilizan los valores de µy Σ para cada distribuci´on normal. Este es el resultado que probaremos de manera te´orica en esta subsecci´on, cuya demostraci´on ser´a muy similar al que llev´abamos a cabo en el apartado de existencia. Por ´ultimo, decir que las propiedades de consistencia permiten extraer informaci´on sobre la capacidad del m´etodo para identificar patrones subyacentes en los datos y generar agrupamientos coherentes y significativos. Definici´on 3.9. Sea X0un vector aleatorio definido en el espacio probabil´ıstico (Ω, σ, P)con llegada en Rd. Sean X1, ..., Xn, ... vectores aleatorios i.i.d. definidas en el mismo espacio con L(X) = L(Xi) = P. Se define Pw nla distribuci´on de probabilidad emp´ırica como aquella funci´on de masa de probabilidad discreta que otorga masa 1 na cada Xi(ω),i= 1, .., n. Es claro que, dado que Pω ndepende del ω∈Ω elegido, tambi´en lo har´an los par´ametros que estamos estimando. Denotamos por {Ψω n}∞ n=1 ={(πn,ω 1, .., πn,ω G, µn,ω 1, .., µn,ω G,Σn,ω 1, ..Σn,ω G)}∞ n=1 la sucesi´on de estimadores emp´ıricos donde cada Ψω nes el vector de par´ametros que maximiza L(Ψ, P, ω) = EPn log   G X j=1 πn,w jϕn,w j  = n X i=1  log   G X j=1 πn,w jϕ(Xi(w); µn,w j,Σn,w j)  (7) donde el tama˜no de la muestra es n. La propiedad de identificabilidad de la familia de mixturas considerada anteriormente es necesaria si queremos asegurar la consistencia del problema de m´axima verosimilitud, como veremos en el siguinte lema. Lema 3.10. Consideramos M′una familia identificable de mixturas finitas definida como M′={f(·; Ψ) = G X j=1 πjf(·;θj) : πj≥0, G X j=1 πj= 1, θj∈Θ, G ∈N,1≤j≤G} Sea Xuna variable aleatoria definida en un espacio probabil´ıstico con densidad f(·; Ψ0)∈M′. Denotemos por EΨ0(Y)el valor esperado de una variable Ycuando el valor verdadero del par´ametro es Ψ0. Entonces, se tiene que EΨ0log f(X, Ψ) <EΨ0log f(X, Ψ0) Para cualquier Ψ= Ψ0en el espacio param´etrico. 35 Demostraci´on Dado que Mes identificable, para vectores de par´ametros distintos Ψ0,Ψ1,f(x; Ψ) =f(x; Ψ0) para al menos un valor de x. Como las funciones de Fson continuas, podemos asegurar que existe un entorno de xdonde se tiene bien f(x; Ψ1)> f(x; Ψ0) o bien f(x; Ψ1)< f(x; Ψ). Esto quiere decir que el conjunto {x∈Rp:f(x; Ψ) =f(x; Ψ0)}tiene medida de Lebesgue estrictamente positiva. Si escribimos log(f(X; Ψ))−log(f(X; Ψ0)) = logf(X;Ψ) f(X;Ψ0), como la funci´on logaritmo es estrictamente c´oncava, la desigualdad de Jensen estricta establece que EΨ0logf(X; Ψ) f(X; Ψ0)<log EΨ0f(X; Ψ) f(X; Ψ0)=Zf(X; Ψ) f(X; Ψ0)dP(Ψ0) = =Zf(X; Ψ) f(X; Ψ0)f(X; Ψ0)dx =Zf(X; Ψ)dx = log(1) = 0 Donde la ´ultima igualdad es consecuencia de que f(X; Ψ) sea una funci´on de densidad. □ Teorema 3.11. En el marco descrito anteriormente, sea {Ψω n}∞ n=1 una sucesi´on de estimadores emp´ıricos donde ndenota el tama˜no de la muestra. Suponemos que X0es un vector aleatorio definido en el espacio probabil´ıstico (Ω, σ, P)con llegada a Rdsiendo Ψ0el ´unico vector de par´ametros que maximiza la verosimilitud (4). Entonces, se tiene que Ψω n−→ Ψ0para ω∈Ω0⊂Ωcon P(Ω0) = 1. Observaci´on 3.12. Se puede calcular la sucesi´on de estimadores emp´ıricos {Ψw n}∞ n=1 (cuya existencia est´a garantizada) resolviendo el problema de maximizaci´on de la funci´on de verosimilitud correspondiente para cada n, a menudo mediante procedimientos iterativos. Cabe destacar que cada uno de los Ψw npuede no ser ´unico. Para demostrar este teorema, probaremos en primer lugar que existe un conjunto compacto K en el espacio param´etrico tal que para un nsuficientemente avanzado, Ψw n∈Kcasi seguro. Sean X1, ..., Xn, ... vectores aleatorios i.i.d. definidos en el mismo espacio con L(X) = L(Xi) = P, y sea Pw nla distribuci´on de probabilidad emp´ırica para X1(ω), ..., Xn(ω). De nuevo, impondremos las restricciones (PR) y(ER) sobre P que enunci´abamos en el apartado anterior. Se verificar´a una versi´on de los dos lemas que utiliz´abamos previamente: Lema 3.13. Dada la sucesi´on {Ψw n}∞ n=1 y dadas las restricciones (ER) y(PR), no puede darse ni mn→0ni Mn→ ∞, donde mnyMndenotaban el autovalor m´as peque˜no y el autovalor m´as grande de las matrices de covarianza respectivamente. Este lema se prueba con un razonamiento similar al que llev´abamos a cabo en el lema (3.3). En este caso, tambi´en necesitamos una cota del tipo EPnm´ax 1≤j≤G∥x−µn,w j∥2≥h′ Para un h′>0. Este resultado auxiliar est´a disponible en el Ap´endice en (5.6) y utiliza el argumento de que la clase de las bolas en Rpes una clase de Glivenko-Cantelli. Por otro lado, tambi´en se verifica el siguiente resultado: Lema 3.14. Dada la sucesi´on {Ψw n}∞ n=1 y dadas las restricciones (ER) y(PR), es posible elegir los vectores de medias emp´ıricos µn,w 1, .., µn,w Gde tal manera que sus normas est´en uniformemente acotadas con probabilidad 1. 36 De nuevo, la demostraci´on de este lema es muy similar a la que expon´ıamos anteriormente en (3.5). Es decir, dado que para los pesos πn,w jera sencillo encontrar un compacto que los contuviese, podemos afirmar que existen n0yKtal que Ψw n∈Ksi n≥n0. A continuaci´on, ahora que sabemos que la sucesi´on {Ψw n}∞ n=1 es uniformemente ajustada, buscamos probar la convergencia casi seguro de Ψw na Ψ0. Para ello, recurriremos a varios resultados presentes en el libro Vaart and Wellner [12], donde se utiliza una extensi´on del Teorema de Glivenko- Cantelli para clases de conjuntos y funciones m´as generales. En el Ap´endice se recoge m´as informaci´on sobre los conceptos necesarios para plantear estos resultados, como las clases de Glivenko-Cantelli o las de Vapnik-Chervonenkis. Consideremos, dado un compacto Ken el espacio param´etrico, la clase de funciones H=nlogG X j=1 πjϕj:θ∈Ko El Teorema 3.2.3 presente en Vaart and Wellner [12] nos permite asegurar que se da la convergencia de Ψw na Ψ0con probabilidad 1 si demostramos que la clase Hes de Glivenko Cantelli. Es decir, buscamos demostrar la siguiente proposici´on. Proposici´on 3.15. La clase Hdefinida anteriormente satisface ∥Pw n−P∥H= sup h∈H |Pw n(h)−P(h)|= sup h∈H |ZhdPw n−ZhdP| −→ 0c.s. Demostraci´on Para probar esto, definiremos la clase G=nIB(·) logG X j=1 πjϕj:θ∈Ko Donde Bes un conjunto compacto fijo. Demostraremos que esta clase de funciones es de Glivenko- Cantelli, lo que nos ayudar´a a concluir que Htambi´en lo es. Se tiene que Φ = {ϕ:µ∈Rp,Σ∈Mp×p} es un espacio vectorial de funciones medibles de dimensi´on finita. Por el lema 2.6.15 de [12] , sabemos que su dimensi´on de Vapnik-Chervonenkis es finita y por lo tanto es una clase V C. Dado que π1, .., πGson pesos, PG j=1 πjϕj:θ∈Kes la envoltura convexa de Φ, por lo que es una clase envolvente (o clase H) de Vapnik Chervonenkis. Si aplicamos el Teorema 2.10.20 presente en [12] con ϕ(x) = IB(x) log(x), se tiene que Gsatisface la condici´on de entrop´ıa uniforme. Como est´a acotada uniformemente, Ges una clase de Glivenko-Cantelli. A continuaci´on, veremos que existen aybpara acotar el tama˜no de cada elemento hde H: |h(x)| ≤ a∥x∥2+b Cota inferior: 37 Dado que Kes un compacto y las matrices son definidas positivas, existen constantes positivas myMque acotan superior e inferiormente todos los autovalores de las Gmatrices de covarianza. Resulta entonces: G X j=1 πjϕj≥ G X j=1 πj 1 (2πM)−p/2exp {−1 2m∥x−µj∥2} Por la identidad del paralelogramo, se tiene que ∥x−µ∥2≤2(∥x∥2+∥µ∥2). Por lo tanto, resulta que G X j=1 πj 1 (2πM)−p/2exp {−1 2m∥x−µj∥2} ≥ exp {−1 m∥x∥2} G X j=1 πj 1 (2πM)−p/2exp {−1 m∥µj∥2} Hab´ıamos demostrado anteriormente que m´ax1≤j≤G∥µj∥<∞, as´ı que existen a′yb′tales que log   G X j=1 πjϕj ≥a′∥x∥2+b′ Cota superior: Sabiendo que existen m, M > 0, escribimos G X j=1 πjϕj≤ G X j=1 πj 1 (2πm)−p/2exp {−1 2M∥x−µj∥2} ≤ G X j=1 πj 1 (2πm)−p/2 Donde la ´ultima desigualdad se tiene ya que el exponente de la funci´on exponencial es negativo y por lo tanto, toda ella es un factor menor que uno. Llegamos a que existen a′′ yb′′ con log   G X j=1 πjϕj ≤a′′∥x∥2+b′′ De esta manera, tomando adecuadamente aybteniendo en cuenta las dos desigualdades, llegamos a que |h(x)| ≤ a∥x∥2+b∀h∈H. Sabiendo esto, para cada h∈Hy para Bun compacto prefijado de Rp, se tiene |EPw n[h(·)] −EP[h(·)]|=|EPw n[h(·)IB(·)] + EPw n[h(·)IBc(·)] −EP[h(·)IB(·)] −EP[h(·)IBc(·)]| ≤ ≤ |EPw n[h(·)IB(·)] −EP[h(·)IB(·)]|+ 38 +|EPw n[(a∥x∥2+b∥x∥+c)IBc(·)] −EP[(a∥x∥2+b∥x∥+c)IB(·)]| La expresi´on anterior tiende a 0 cuando ntiende a infinito teniendo en cuenta que h(·)IB(·)∈G, que es una clase de Glivenko-Cantelli. Esto provoca que ∥Pw n−P∥H= sup h∈H |ZhdPw n−ZhdP|= sup h∈H |EPw n(h(·)) −EP[h(·)]| −→ 0c.s. Y que por lo tanto Hsea una clase de Glivenko-Cantelli. □ Una vez probados estos lemas, procedemos a demostrar el Teorema 3.11 de Consistencia que enunci´abamos anteriormente. Demostraci´on Teniendo en cuenta que Ψ0es ´unico y que los lemas anteriores que aseguran que {Ψw n}∞ n=1 es uniformemente ajustada y Hes una clase de Glivenko-Cantelli, podemos afirmar que Ψw n→Ψ0con probabilidad uno utilizando el Corolario 3.2.3 presente en [12]. □ 39 4. Aspectos Metodol´ogicos Tras presentar el modelo de mixturas y haber estudiado sus principales propiedades matem´aticas, abordamos en esta secci´on el c´alculo pr´actico de estimaciones de los par´ametros del modelo, limit´andonos al modelo de mezclas gaussiano. En primer lugar, hablaremos sobre las restricciones geom´etricas que se pueden introducir en el modelo para aliviar el n´umero de par´ametros a estimar, tras lo cual introduciremos algunos criterios ´utiles para inferir las caracter´ısticas geom´etricas de los grupos. Otra cuesti´on importante es determinar el n´umero de componentes del modelo cuando esta informaci´on no es conocida a priori y debemos inferirla de los datos, al igual que considerar el modelado de datos con agrupaciones no gaussianas fusionando clusters creados a partir del modelo de mezclas gaussiano. Recientemente, se han propuesto multitud de criterios para todas estas cuestiones, muchos de los cuales tienen resultados alentadores en estudios emp´ıricos para estudiar su actuaci´on. Por otro lado, presentaremos el algoritmo EM incidiendo en el impacto que tuvo en el problema de ajuste del modelo de mezclas y detallaremos los pasos del mismo. Comentaremos someramente sus propiedades de convergencia bajo ciertas restricciones e introduciremos algunos procedimientos relativos a la incializaci´on del algoritmo. 4.1. Restricciones geom´etricas en el modelo de mezclas Volvamos a la situaci´on en la que consideramos un conjunto de nobservaciones p-variantes {x1, .., xn}y buscamos expresar su funci´on de densidad a trav´es de una combinaci´on convexa de Gfunciones de densidad gaussianas: f(xi) = G X j=1 πjϕ(xi;µj,Σj) = G X j=1 πjϕj(xi),1≤i≤n(8) Analicemos el n´umero de par´ametros que debemos estimar en este modelo: 1. En primer lugar, ser´a necesario estimar (G−1) de los pesos πjya que πGse puede obtener como 1 −PG−1 j=1 πj. 2. En el caso de los vectores de medias, ser´a necesario estimar Gvectores p−dimensionales, por lo que tendr´ıamos que calcular Gp par´ametros. 3. Para finalizar, dado que contamos con Gmatrices de covarianza (sim´etricas, solo es necesario obtener el tri´angulo superior de la matriz para reconstruirla) que son de tama˜no p×p, en este caso acumulamos Gp(p+1) 2par´ametros a estimar. En total, el n´umero de par´ametros para el modelo planteado en (8) ser´ıa (G−1) + Gp +Gp(p+1) 2. En un caso muy sencillo en el que cont´asemos con datos bidimensionales y busc´asemos ajustar un modelo con dos componentes, tendr´ıamos que calcular 11 par´ametros, algo que sin algoritmos y computadoras que realicen los c´alculos resultar´ıa bastante tedioso. Pero no solo eso: este n´umero puede crecer r´apidamente y ser bastante elevado en el momento en el que aumentemos el n´umero de variables p, considerando conjuntos de datos en dimensiones mayores. Esto hace patente la necesidad 40 de pertencer a cualquiera de los Ggrupos). Resulta que Entr(G)≥Entr(G−1), ya que al a˜nadir grupos se tienen m´as posibilidades a la hora de asignar un individuo a un cluster y por lo tanto mayor incertidumbre para cada observaci´on. La idea fundamental de este procedimiento es escoger el par de grupos que minimizan el incremen- to de entrop´ıa al ser fusionados. Sean CkyChdos clusters de la partici´on en Ggrupos {C1, .., CG} yh, k sus respectivos ´ındices. Si unimos CkyCh, dado que son disjuntos, el nuevo cluster Ch∪Ck con ´ındice {h, k}tiene la siguiente probabilidad a posteriori: ˆzG i,{h,k}= ˆzG ik + ˆzG ih Mientras que el resto de valores ˆzG ij se mantienen iguales, con j=k, h. La entrop´ıa resultante ser´a −n X i=1 (ˆzG ih + ˆzG ik) log(ˆzG ih + ˆzG ik) + X j=h,k ˆzG ij log(ˆzG ij) Si restamos ambas dos entrop´ıas, buscamos dos clusters ChyCkque maximicen la diferencia − n X i=1 (ˆzG ih + ˆzG ik) log(ˆzG ih + ˆzG ik) + n X i=1 ˆzG ik log(ˆzG ik) + ˆzG ih log(ˆzG ih) Una vez seleccionado el par de clusters, se renombra su uni´on y se acualizan las probabilidades a posteriori ˆzG−1 ij (Ahora tenemos G−1 grupos). El objetivo es continuar uniendo grupos hasta que se provoca un incremento grande en la entrop´ıa. Usualmente, se puede utilizar el criterio del codo en un gr´afico que compare el n´umero de componentes con la entrop´ıa para decidir en qu´e momento parar. Expondremos el funcionamiento del m´etodo de Entrop´ıa de Clasificaci´on con el siguiente conjunto de datos. En la siguiente figura, observamos su representaci´on bidimensional y su partici´on en clusters a medida que variamos el n´umero de grupos. 47 Figura 23: Agrupaciones variando G Observamos c´omo este conjunto de datos tiene 5 grupos m´as o menos separados, alguno de apariencia elipsoidal mientras que otros son claramente no gaussianos. Si construimos un gr´afico que muestre la entrop´ıa de clasificaci´on seg´un el n´umero de componentes consideradas, observamos como su variaci´on es peque˜na hasta que consideramos 5 clusters. A partir de ese momento, se dispara, por lo que el m´etodo de la Entrop´ıa de Clasificaci´on sugerir´ıa que existen 5 grupos en este conjunto de datos. Como coment´abamos, el procedimiento parte del modelo escogido por el BIC, que a menudo contiene una partici´on m´as exhaustiva en componentes el´ıpticas de los clusters no gaussianos con el fin de representarlos mejor. Es decir, en este caso el criterio BIC determin´o que el n´umero adecuado de grupos es 8. Por otro lado, es interesante estudiar cu´antas componentes preve´ıamos esperar en el modelo con mejor ICL, por lo que producimos dos gr´aficos que muestran el BIC e ICL respectivamente para los 14 modelos seg´un la clasificaci´on VSO variando el n´umero de componentes. 48 Figura 24: Selecci´on de modelos utilizando BIC e ICL Observamos c´omo el criterio ICL consigue recuperar el n´umero de clusters que visualmente se aprecian en el conjunto de datos, es decir, cinco. 4.3. El algoritmo EM Hace unas d´ecadas, a pesar de contar con ordenadores que permit´ıan realizar c´alculos con mucha mayor rapidez, exist´ıa cierta reticencia a utilizar los modelos de mezclas en el caso multidimensional, posiblemente por la falta de informaci´on sobre algunas cuestiones abiertas que a´un no hab´ıan sido estudiadas en profundidad, como la presencia de m´ultiples m´aximos en la funci´on de verosimilitud y la no acotaci´on de esta en el caso de utilizar componentes gaussianas con distintas matrices de covarianza. En el momento en el que estas dificultades se comprendieron satisfactoriamente y se consensuaron algunas soluciones, se increment´o el uso de estos modelos en la pr´actica. En los sesenta, el ajuste de modelos de mezclas por m´axima verosimilitud hab´ıa sido estudiado en numerosos art´ıculos, pero fue la publicaci´on de Dempster, Laird y Rubin en 1977 sobre el algoritmo EM lo que puso el foco de inter´es en este m´etodo para modelar datos heterog´eneos. El algoritmo EM es una herramienta de uso extendido para el c´alculo del estimador m´aximo verosimil especialmente provechosa en problemas con datos incompletos donde m´etodos como el de Newton-Raphson pueden resultar m´as complicados. Este procedimiento iterativo puede utilizarse no solo en situaciones en las que hay datos faltantes, distribuciones truncadas o agrupaciones en las observaciones, sino en circunstancias donde el car´acter incompleto de los datos no es tan evidente, como puede ser el modelo de mixtura, los modelos log-lineales, convoluciones... Abordaremos su estudio en el problema de la estimaci´on de par´ametros en el modelo de mezclas y observaremos c´omo se ve simplificado por el EM, adem´as de proporcionar en algunos casos expresiones cerradas de los estimadores. 49 4.3.1. Formulaci´on del algoritmo EM Consideramos el caso general del modelo de mixturas en el que la familia de densidades de probabilidad es arbitraria. M´as tarde, particularizaremos estos resultados para el caso del modelo de mezclas gaussiano. Es decir, nos encontramos en el caso en el que buscamos estimar el vector de par´ametros Ψ = (π1, .., πG, θ1, .., θG) dada una muestra aleatoria A={x1, .., xn}a la que hemos ajustado un modelo de mixturas f(xi; Ψ) = G X j=1 πjf(xi;θj), f ∈F={f(·;θ) : θ∈Θ},1≤i≤n Su funci´on de log-verosimilitud correspondiente se escribe como L(Ψ, A) = n X i=1 logG X j=1 πjf(xi;θj) El estimador m´aximo versol´ımil de Ψ se calcula resolviendo la ecuaci´on de verosimilitud ∂L(Ψ, A) ∂Ψ= 0 A menudo en la pr´actica, la funci´on de log-verosimilitud no puede ser maximizada de forma anal´ıtica. En muchos casos, es posible calcular de forma iterativa el MLE de Ψ utilizando el procedimiento de Newton-Raphson o alguno similar cuando el n´umero de par´ametros en el modelo no es demasiado grande. Consideremos las etiquetas zij no observadas que toman el valor 1 si la observaci´on xiproviene de la componente j-´esima y 0 en otro caso (1 ≤i≤n, 1≤j≤G), y denotamos x= (x1, .., xn). Los vectores zi= (zi1, .., ziG) se consideraban una realizaci´on de Z1, .., Znvectores aleatorios independientes igualmente distribuidos definidos en un mismo espacio probabil´ıstico que siguen una distribuci´on multinomial, donde se elige una de Gcategor´ıas con probabilidades π1, .., πG. Asumimos que la densidad de una observaci´on xidado zies G Y j=1 f(xi;θj)zij Es decir, si xiproviene de la componente j, todos los factores distintos de jser´an 1 y solo quedar´a f(xi;θj). La log-verosimilitud de los datos completos es por lo tanto Lc(Ψ, A) = G X j=1 n X i=1 zij(log πj+ log f(xi;θj)) El algoritmo EM cuenta con dos pasos que describimos a continuaci´on: E Step (Expectation Step): Sea Ψ(0) un valor inicial para Ψ. En este paso, se procede al c´alculo de Q(Ψ; Ψ(0)) = EΨ(0) Lc(A, Ψ)|x 50 Es decir, calculamos el valor esperado de la verosimilitud completa condicionado a x, utilizando el ajuste Ψ(0) para Ψ. Como hemos visto, la log-verosimilitud de los datos completos es una funci´on lineal en zij , por lo que basta con calcular la esperanza condicional de Zij dados los datos observados X. Por la definici´on de esperanza condicional, podemos reescribir EΨ(0) Zij|X=x= 1 ·P(Zij = 1|X=x)+0·P(Zij = 0|X=x) = P(Zij = 1, X =x) P(X=x) Por el Teorema de Bayes, la ´ultima expresi´on puede reescribirse como P(Zij = 1, X =x) P(X=x)=P(X=x|Zij = 1)P(Zij = 1) P(X=x)=πjf(xi;θj) f(xi; Ψ(0))=τj(xi; Ψ(0)) Para 1 ≤i≤n, 1≤j≤G. La cantidad τj(xi; Ψ(0)) es la probabilidad a posteriori de que la observaci´on i−´esima pertenezca a la j−´esima componente. Sustituyendo esto en la esperanza condicional de la log-verosimilitud, se tiene Q(Ψ; Ψ(0)) = G X j=1 n X i=1 τj(xi; Ψ(0))(log πj+ log f(xi;θj)) M Step (Maximization Step): El primer paso consiste en elegir Ψ(1) el valor que maximiza Q(Ψ; Ψ(0)) respecto de Ψ. Si las zij fueran datos observados, el MLE de los pesos se calcular´ıa como ˆπj= n X i=1 zij n Por lo tanto, en este primer paso π(1) j=Pn i=1 τj(xi;Ψ(0)) n. Es decir, al formar la estimaci´on del peso de la componente j−´esima en la primera iteraci´on, π(1) j, se tiene en cuenta la contribuci´on de cada observaci´on xipor medio de su probabilidad posterior (evaluada en ese momento) de pertenecer a la j-´esima componente del modelo de mezclas. Esta estimaci´on por medio del EM tiene por lo tanto una interpretaci´on bastante intuitiva. A la hora de estimar los par´ametros θ(1) jcorrespondientes, debemos buscar una ra´ız de la ecuaci´on G X j=1 n X i=1 τj(xi; Ψ(0))∂(log f(xi;θj)) ∂θ = 0 Donde θdenota (θ1, .., θG). La raz´on por la que el algoritmo EM resulta interesante en este contexto se debe al hecho de que a menudo la soluci´on de esta ecuaci´on tiene una forma cerrada, como veremos m´as adelante en el caso del modelo de mezclas gaussiano. 51 Hemos descrito los dos pasos del algoritmo EM en la primera iteraci´on. De manera general, se tendr´ıa: 1. Inicializamos el vector de par´ametros Ψ(0). 2. En la iteraci´on k≥1, se efect´uan E-Step: Se calcula Q(Ψ; Ψ(k−1)) = EΨ(k−1) Lc(A, Ψ)|x= G X j=1 n X i=1 τj(xi; Ψ(k−1))(log πj+ log f(xi;θ(k−1) j)) M-Step: Se busca el vector de par´ametros Ψ(k)tal que Ψ(k)= arg max Ψ Q(Ψ; Ψ(k−1)) Se calcula L(Ψ(k))−L(Ψ(k−1)) 3. Cuando la diferencia calculada anteriormente es suficientemente peque˜na o no ha habido un cambio significativo en los par´ametros, se detiene el proceso. En caso contrario, se actualiza el valor del vector de par´ametros Ψ(k−1) := Ψ(k)y se repite el segundo paso. Figura 25: Ejemplo de actuaci´on del algoritmo EM En esta imagen observamos la actuaci´on del algoritmo EM a la hora de ajustar tres componentes de mixtura al conjunto de datos de la izquierda. En el gr´afico central, observamos los valores que toma la funci´on de log-verosimilitud en cada una de las iteraciones. Observamos c´omo a partir de la iteraci´on 15, los cambios en la log-verosimilitud son muy peque˜nos. Una vez superada la tolerancia, el algoritmo se detiene y ofrece la agrupaci´on final en tres componentes que vemos a la derecha. 52 Por otro lado, en aquellas situaciones en las que existe un nivel alto de superposici´on entre los grupos, la convergencia del algoritmo es mucho m´as lenta y necesita un n´umero mayor de iteraciones. En este caso, vemos como la actuaci´on no resulta del todo satisfactoria a la hora de hallar las tres componentes del conjunto de datos y observamos que fueron necesarias m´as de 300 iteraciones para la convergencia del algoritmo. Figura 26: Ejemplo de actuaci´on del algoritmo EM 4.3.2. El algoritmo EM en el modelo de mezclas gaussiano Si consideramos el caso particular en el que las componentes de la distribuci´on son gaussianas, existe una expresi´on cerrada para el c´alculo de los par´ametros por medio del algoritmo EM. Recordemos que en este caso trabajamos con Ψ = (π1, .., πG, µ1, .., µG,Σ1, .., ΣG) el vector de par´ametros yϕ(·;µj,Σj) la j−´esima componente de la mixtura. Consideramos en primer lugar el caso heteroced´astico (no hay restricciones sobre las matrices de covarianza): Si nos encontramos en la iteraci´on k−´esima en la etapa de maximizaci´on, los vectores de medias y matrices de covarianzas se actualizan del siguiente modo: µ(k) j= n X i=1 z(k−1) ij xi Pn i=1 z(k) ij ,1≤j≤G Σ(k) j= n X i=1 z(k−1) ij (xi−µ(k) j)(xi−µ(k) j)T Pn i=1 z(k−1) ij ,1≤j≤G En la pr´actica, cuando la naturaleza de los datos lo permite, se trabaja bajo la hip´otesis de que las componentes tienen matriz de covarianza com´un Σ desconocida. En este caso de homocedasticidad, se sustituye la actualizaci´on correspondiente de las Σ(k) jpor Σ(k)= G X j=1 n X i=1 z(k−1) ij (xi−µ(k) j)(xi−µ(k) j)T n 53 En este caso, el MLE existe como el maximizador global de la verosimilitud y adem´as es fuertemente consistente. En el caso heterosced´astico, la maximizaci´on de la verosimilitud sin ninguna restricci´on podr´ıa llevarnos a una soluci´on degenerada. Una vez seleccionamos un modelo, podemos plantear condiciones sobre el ratio entre los autovalores de las matrices de covarianza en cada iteraci´on M(k) m(k)≥c, c ≥1 M(k)mayor autovalor de las matrices de covarianza en la iteraci´on k-´esima m(k)menor autovalor de las matrices de covarianza en la iteraci´on k-´esima De la misma forma, cuando buscamos elegir un modelo por medio del BIC, podr´ıamos descartar aquellos modelos cuyo ratio entre autovalores de las matrices finales de covarianza est´an por debajo del umbral c. En este caso, la cuesti´on se reduce a c´omo elegir este par´ametro para asegurar que la restricci´on que hemos hecho del espacio param´etrico contienen el verdadero valor del vector de par´ametros Ψ, o para garantizar que no descartamos modelos que podr´ıan ajustar mejor los datos pero que no superan esta restricci´on. La elecci´on de cdepender´a de las particularidades del conjunto de datos en cada caso. 4.3.3. Convergencia del Algoritmo EM Nuestro objetivo es recoger br´evemente los resultados m´as importantes disponibles sobre la teor´ıa b´asica relativa a la convergencia del algoritmo EM. Para ello, utilizaremos como gu´ıa el libro de McLachlan [10] donde se encuentran las pruebas de los teoremas que enunciaremos. En primer lugar, demostraremos que la funci´on de log-verosimilitud L(Ψ) que buscamos maximizar es no decreciente al actualizar los iterantes. Proposici´on 4.1. (Motonon´ıa del EM) Consideremos x= (x1, .., xn)una muestra aleatoria, Ψ = (π1, .., πG, θ1, .., θG)el vector de par´ametros a estimar y L(Ψ, A)la funci´on de verosimilitud L(Ψ, x) = n Y i=1G X j=1 πjf(xi;θj)= n Y i=1 f(xi,Ψ) Consideramos z= (z1, .., zn)los datos no observados donde zi= (zi1, .., ziG)es el vector de etiquetas de la observaci´on i−´esima. Sea {Ψ(k)}k≥0una sucesi´on de iterantes obtenidos por medio del algoritmo EM. Entonces, para todo k≥0, se tiene que L(Ψ(k+1))≥L(Ψ(k)) Demostraci´on Recordemos que la funci´on de verosimilitud completa ten´ıa la expresi´on Lc(Ψ, x, z) = n Y i=1 n Y i=1πjf(xi;θj)zij = n Y i=1 fc(Ψ, xi, zi) 54 Sea k(z|x, Ψ) la densidad condicional de Zdado X=x, es decir k(z|x, Ψ) = Qn i=1 fc(Ψ, xi, zi) Qn i=1 f(xi,Ψ) La log-verosimilitud se escribe entonces L(Ψ, x) = log(L(Ψ, x)) = log n Y i=1 f(xi,Ψ)= log n Y i=1 fc(Ψ, xi, zi)−log k(z|x, Ψ) Si tomamos esperanzas condicionales respecto de Zdado X=xa ambos lados y utilizando el ajuste Ψ(k)para Ψ, se tiene que L(Ψ, x) = EΨ(k)log n Y i=1 f(xi,Ψ)−EΨ(k)log k(z|x, Ψ)=Q(Ψ; Ψ(k))−H(Ψ; Ψ(k)) De esta manera, se tiene que L(Ψ(k+1), x)−L(Ψ(k), x) = Q(Ψ(k+1); Ψ(k))−Q(Ψ(k); Ψ(k))−H(Ψ(k+1); Ψ(k))−H(Ψ(k); Ψ(k)) (10) El iterante Ψ(k+1) se elige de tal manera que maximice Q(Ψ; Ψ(k)), por lo que Q(Ψ(k+1); Ψ(k))≥Q(Ψ(k); Ψ(k)) Y con esto, la resta del primer par´entesis es no negativa. Si logramos ver que H(Ψ(k+1); Ψ(k))−H(Ψ(k); Ψ(k))≤0 se tendr´ıa que la expresi´on (10) es no negativa y habr´ıamos acabado. Para un Ψ cualquiera H(Ψ; Ψ(k))−H(Ψ(k); Ψ(k)) = EΨ(k)log k(z|x, Ψ) −log k(z|x, Ψ(k))|x= =EΨ(k)logk(z|x, Ψ) k(z|x, Ψ(k))x≤logEΨ(k)k(z|x, Ψ) k(z|x, Ψ(k))x Donde la ´ultima desigualdad es consecuencia de la desigualdad de Jensen ya que la funci´on logaritmo es c´oncava. Esta ´ultima expresi´on, recordando la definici´on de esperanza condicional respecto de Z dado xy utilizando el ajuste Ψ(k)para Ψ, puede reescribirse como logEΨ(k)k(z|x, Ψ) k(z|x, Ψ(k))x=Zk(z|x, Ψ) k(z|x, Ψ(k))k(z|x, Ψ(k))dx = log(1) = 0 Por lo tanto, hemos probado que la verosimilitud no decrece al iterar el algoritmo EM. □ Una consecuencia de esto es que si ˆ Ψ es el MLE que maximiza L(Ψ, x), necesariamente Q(ˆ Ψ; ˆ Ψ) ≥Q(Ψ; ˆ Ψ) Para cualquier Ψ del espacio param´etrico. Si no fuese as´ı, llegar´ıamos a una contradicci´on ya que ˆ Ψ no maximizar´ıa L(Ψ, x). 55 Bajo condiciones no demasiado restrictivas, podemos asegurar la convergencia de los iterantes hacia puntos estacionarios, es decir, donde el gradiente de la funci´on de verosimilitud se anula. Teorema 4.2. Supongamos que la funci´on Q(Ψ; Φ) es continua en ambas variables. Entonces todos los puntos l´ımites de cualquier sucesi´on de iterantes {Ψ(k)}k≥0del algoritmo EM son puntos estacionarios de L(Ψ, x)yL(Ψ(k))converge mon´otonamente hacia un valor L∗=L(Ψ∗, x)para un punto estacionario Ψ∗ Esto provoca que en numerosas situaciones pr´acticas podamos encontrar un m´aximo local L∗= L(Ψ∗, x) por medio del algoritmo EM. En general, si L(Ψ, x) tiene varios puntos estacionarios, la convergencia de los iterantes hacia cualquier punto que anula el gradiente de la funci´on de verosimilitud (m´aximos locales o globales, puntos de silla) depende de la elecci´on del iterante inicial Ψ(0). En el caso de considerar el problema de mezclas gaussiano, este teorema nos garantiza que imponiendo una condici´on sobre el ratio de los autovalores de las matrices de covarianza como la ya planteada, podemos asegurar que el algoritmo converger´a hacia un punto estacionario. Por ´ultimo, si la funci´on de verosimilitud es unimodal con un solo punto estacionario en el interior del espacio param´etrico y se cumplen ciertas condiciones de regularidad, tenemos el siguiente resultado. Corolario 4.3. Supongamos que L(Ψ, x)es unimodal con Ψ∗es su ´unico punto estacionario y ∂Q(Ψ;Φ) ∂Ψes continua en ambas variables. Entonces cualquier sucesi´on {Ψ(k)}k≥0de iterantes del algoritmo EM converge al ´unico maximizador Ψ∗de L(Ψ; x), es decir, su MLE. 4.3.4. Inicializaci´on del Algoritmo EM Dado que usualmente existen m´aximos locales en la funci´on de verosimilitud, a menudo el vector de par´ametros obtenido por medio del algoritmo EM depende del valor inicial Ψ(0) elegido. Es por ello que existen diferentes t´ecnicas de inicializaci´on m´as o menos sofisticadas con el fin de seleccionar un vector de par´ametros inciales que permitan un buen punto de partida a la hora de comenzar el algoritmo. En la figura siguiente, mostramos la actuaci´on del algoritmo EM cuando la inicializaci´on se realiza por medio de una partici´on aleatoria en dos grupos equilibrados. En este caso, observamos que el algoritmo tiene una convergencia r´apida ya que a partir de la octava iteraci´on se producen cambios muy peque˜nos en la log-verosimilitud. En la figura de la derecha podemos observar la partici´on resultante efectuada por el EM, donde las l´ıneas muestran los valores que han tomado los dos vectores de medias en cada iteraci´on. Dado que la partici´on inicial en grupos es aleatoria, los iterantes iniciales est´an muy pr´oximos entre s´ı, pr´acticamente en medio de la nube de puntos. Los tri´angulos azul y verde son las estimaciones finales de los centros de las distribuciones. 56 que pueda separarlos de todas las formas posibles. Por lo tanto, la dimensi´on VC de Les 2. 2. La dimensi´on VC de la clase Fes infinita, ya que dado cualquier conjunto finito de puntos A en R, sus propios subconjuntos ya pertenecen a F. Dado que podemos aumentar el tama˜no de Aarbitrariamente, la dimensi´on VC de Fes infinita. En este documento utilizaremos el siguiente resultado que permite relacionar las clases VC con las clases de Glivenko-Cantelli. Teorema 5.4. Una clase Ces de Glivenko-Cantelli si y solo si es una clase Vapnik-Chervonenkis La idea intuitiva que subyace bajo este resultado es que algunas clases de conjuntos son “demasiado grandes” para que se pueda verificar el Teorema de Glivenko-Cantelli. Por ejemplo, si consideramos la clase de las bolas en Rp, que denotamos por B, es claro que su dimensi´on VC es finita ya que a partir de un cardinal del conjunto A, una bola no es capaz de separar todos los subconjuntos de A. Esto provoca que Bsea una clase de Glivenko Cantelli, y que por lo tanto sup B∈B |Pn(B)−P(B)| −→ 0c.s. 5.2. Algunos resultados auxiliares Recordemos que las restricciones que establec´ıamos al plantear la existencia y consistencia de soluciones en el modelo de mezclas eran las siguientes: 1. (PR) Restricciones sobre la distribuci´on P: La distribuci´on P no debe est´ar concentrada en Gpuntos. P debe de tener momento de segundo orden finito, esto es: EP(∥·∥2)<∞. En el caso de estar trabajando con recortes, esta hip´otesis no ser´ıa necesaria. 2. (ER) Restricciones sobre el ratio entre autovalores: Debe cumplirse, para una constante c≥1, que Mn/mn≤c Donde Mn= m´ax 1≤j≤Gm´ax 1≤l≤pλjl mn= m´ın 1≤j≤Gm´ın 1≤l≤pλjl Lema 5.5. Sea {Ψn}∞ n=1 una sucesi´on de estimadores tal que l´ımn→∞ L(Ψn, P)>−∞. Dadas las restricciones (PR) y(ER), existe h > 0tal que EPm´ax1≤j≤G∥x−µn j∥2≥h > 0. 63 Demostraci´on Dado que Pno est´a concentrada en Gpuntos, existen y1, .., yG+1 ∈Rptales que P(B(yj, ϵ)) > δϵ>0, donde B(y, ϵ) denota la bola de centro ycon radio ϵ. Consideramos ahora ϵ0<m´ın1≤j<k≤G+1 ∥yj−yk∥ 2y las bolas B(yj, ϵ0) con j= 1, .., G + 1. Observamos que, dado que son bolas abiertas, son disjuntas dos a dos. Dado que tenemos Gvectores µn j yG+ 1 puntos yj, sabemos que existe un ´ındice, vamos a decir que es el G+ 1 tras reordenarlos, de tal manera que la bola B(yG+1, ϵ0) no contiene a ninguno de los µn j. Denotamos como BG+1 esta bola. De esta manera, podemos escribir EP( m´ax 1≤j≤G∥x−µn j∥2)≥EP( m´ın 1≤j≤G∥x−µn j∥2) = =EP( m´ın 1≤j≤G∥x−µn j∥2·IBG+1 ) + EP( m´ın 1≤j≤G∥x−µn j∥2·IBc G+1 )≥EP( m´ın 1≤j≤G∥x−µn j∥2·IBG+1 ) La ´ultima desigualdad es consecuencia de que m´ın1≤j≤G∥X−µn j∥2sea una variable positiva. Dado que x∈BG+1 yµn jno, se tiene que EP( m´ın 1≤j≤G∥x−µn j∥2·IBG+1 )≥EP m´ın 1≤j<k≤G+1 ∥yj−yk∥ 22 ·IBG+1 !≥ϵ2 0·P(BG+1)>0 □ Lema 5.6. Sea {Pn}∞ n=1 una sucesi´on de distribuciones de probabilidad emp´ıricas que satisfacen (ER). Sea {Ψn}∞ n=1 una sucesi´on de estimadores muestrales. Si P satisface (PR), entonces existe h′>0tal que EPnm´ax1≤j≤G∥x−µn j∥2≥h′. Demostraci´on Hemos comentado que supB∈B|Pn(B)−P(B)| −→ 0c.s. por ser Bde Glivenko-Cantelli. La expresi´on anterior se puede reescribir como sup B∈B |Pn(B)−P(B)|= sup B∈B |ZIBdPn−ZIBdP|= sup B∈B |EPn(IB)−EP(IB)| −→ 0c.s. Esto nos permite afirmar que, de manera uniforme en B, para un nelevado EPn(IB)dista muy poco de EP(IB). Por lo tanto, tomando un ´ındice suficientemente avanzado podr´ıamos sustituir una esperanza por la otra y conseguir´ıamos demostrar EPnm´ax1≤j≤G∥x−µn j∥2≥h′siguiendo unos pasos an´alogos a los que utiliz´abamos en el lema anterior para probar EPm´ax1≤j≤G∥x−µn j∥2≥h. □ 64 Referencias [1] Billingsley, P. (2013). Convergence of probability measures. John Wiley and Sons. [2] Bouveyron, C., Celeux, G., Murphy, T. B., and Raftery, A. E. (2019). Modelbased clustering and classification for data science: with applications in R (Vol. 50). Cambridge University Press. [3] Fritz, H., Garc´ ıa-Escudero, L. A., Mayo-Iscar, A. (2013). Computational Statistics and Data Analysis, 61, 124-136. Computational Statistics and Data Analysis, 61, 124-136. [4] Garc´ ıa-Escudero, L. A., Gordaliza, A., and Mayo-´ Iscar, A. (2014). A constrained robust proposal for mixture modeling avoiding spurious solutions.. Advances in Data Analysis and Classification, 8(1), 27-43. [5] Garc´ ıa-Escudero, L. A., Gordaliza, A., Matr´ an, C., and Mayo-Iscar, A. (2008). A general trimming approach to robust cluster analysis. 1324-1345. [6] Garc´ ıa-Escudero, L. A., Gordaliza, A., Matr´ an, C., and Mayo-Iscar, A. (2015). Avoiding spurious local maximizers in mixture modeling. Statistics and Computing, 25, 619-633. [7] Hathaway, R. J. (1985). A constrained formulation of maximum-likelihood estimation for normal mixture distributions. . The Annals of Statistics, 13(2), 795-800. [8] McLachlan, G.J and Peel, D.(1999) Computing Issues for the EM Algorithm in Mixture Models. Interface’99, Schaumbury Illinois, Virginia: Interface Foundation of North America. [9] McLachlan, G. J., Lee, S. X., and Rathnayake, S. I. (2019). Finite mixture models. Annual review of statistics and its application, 6, 355-378. [10] McLachlan, G. J., and Krishnan, T. (2007). The EM algorithm and extensions. John Wiley and Sons. [11] Teicher, H. (1963): Identifiability of Finite Mixtures. Mathematical Statistics, 34,4,1265–1269 [12] Van der Vaart, A., and Wellner, J. A. 1996. Weak convergence and empirical processes. Springer New York. [13] Yakowitz, S. J., and Spragins, J. D. (1968). On the identifiability of finite mixtures. The Annals of Mathematical Statistics, 39(1), 209-214. 65