Problemas de singularidad en el modelo de mezcla gaussiana

En el capítulo 9 del libro Reconocimiento de patrones y aprendizaje automático, hay esta parte sobre el modelo de mezcla gaussiana:

Para ser honesto, realmente no entiendo por qué esto crearía una singularidad. ¿Puede alguien explicarme esto? Lo siento, pero solo soy un estudiante universitario y un novato en el aprendizaje automático, por lo que mi pregunta puede sonar un poco tonta, pero por favor ayúdenme. Muchas gracias

Si queremos ajustar un Gaussiano a un único punto de datos usando la máxima probabilidad, obtendremos un Gaussiano muy puntiagudo que "colapsa" a ese punto. La varianza es cero cuando solo hay un punto, que en el caso gaussiano multivariado, conduce a una matriz de covarianza singular, por lo que se llama problema de singularidad.

Cuando la varianza llega a cero, la probabilidad del componente gaussiano (fórmula 9.15) llega al infinito y el modelo se sobreajusta. Esto no ocurre cuando ajustamos solo un Gaussiano a varios puntos ya que la varianza no puede ser cero. Pero puede suceder cuando tenemos una mezcla de gaussianos, como se ilustra en la misma página de PRML.

Actualización :
el libro sugiere dos métodos para abordar el problema de la singularidad, que son

1) restablecer la media y la varianza cuando se produce la singularidad

2) usando MAP en lugar de MLE agregando un previo.

Recordemos que este problema no surgió en el caso de una única distribución gaussiana. Para comprender la diferencia, tenga en cuenta que si un solo Gaussiano colapsa en un punto de datos, contribuirá con factores multiplicativos a la función de probabilidad que surge de los otros puntos de datos y estos factores irán a cero exponencialmente rápido, dando una probabilidad general que va a cero en lugar de que el infinito

También estoy un poco confundido por esta parte, y aquí está mi interpretación. Tome el caso 1D por simplicidad.

Cuando un único "colapso" gaussiano en un punto de datos , es decir, μ = x i , la probabilidad general se convierte en:$x_i$$\mu=x_i$

$$p(\mathbf{x}) = p(x_i) p(\mathbf{x}\setminus{i}) = (\frac{1}{\sqrt{2\pi}\sigma}) (\prod_{n \neq i}^N \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(x_n-\mu)^2}{2\sigma^2}} )$$

Usted ve como , el término a la izquierda p ( x i ) → ∞ , que es como el caso patológico en GMM, pero el término a la derecha, que es la probabilidad de otros puntos de datos p ( x ∖ i ) , todavía contiene términos como e - ( x n - μ ) $\sigma \to 0$$p(x_i) \to \infty$$p(\mathbf{x}\setminus{i})$ que→0exponencialmente rápido comoσ→0, por lo que el efecto general sobre la probabilidad es que vaya al cero.$e^{-\frac{(x_n-\mu)^2}{2\sigma^2}}$$\to 0$$\sigma \to 0$

El punto principal aquí es que cuando se ajusta un único gaussiano, todos los puntos de datos tienen que compartir un conjunto de parámetros , a diferencia del caso de la mezcla donde un componente puede "enfocarse" en un punto de datos sin penalizar la probabilidad general de datos .$\mu, \sigma$

Esta respuesta dará una idea de lo que está sucediendo que conduce a una matriz de covarianza singular durante el ajuste de un GMM a un conjunto de datos, por qué sucede esto y qué podemos hacer para evitarlo.

Por lo tanto, es mejor comenzar recapitulando los pasos durante el ajuste de un modelo de mezcla gaussiana a un conjunto de datos.

0. Decida cuántas fuentes / grupos (c) desea ajustar a sus datos
1. Inicialice la media de los parámetros , covarianza Σ c , y fracción_por_clase π c por grupo c $\mu_c$$\Sigma_c$$\pi_c$

$\underline{E-Step}$

2. Calcule para cada punto de datos la probabilidad r i c de que el punto de datos x i pertenece al grupo c con: r i c = π c N ( x i | μ c , Σ c$x_i$$r_{ic}$$x_i$
   
   dondeN(x|μ,Σ)describe la Gaussiana multivariada con: N(xi,μc,Σc)= $$r_{ic} = \frac{\pi_c N(\boldsymbol{x_i} \ | \ \boldsymbol{\mu_c},\boldsymbol{\Sigma_c})}{\Sigma_{k=1}^K \pi_k N(\boldsymbol{x_i \ | \ \boldsymbol{\mu_k},\boldsymbol{\Sigma_k}})}$$ $N(\boldsymbol{x} \ | \ \boldsymbol{\mu},\boldsymbol{\Sigma})$
   
   ricnos da para cada punto de datosxila medida de:Probabilitythatxibelongstocla $$N(\boldsymbol{x_i},\boldsymbol{\mu_c},\boldsymbol{\Sigma_c}) \ = \ \frac{1}{(2\pi)^{\frac{n}{2}}|\Sigma_c|^{\frac{1}{2}}}exp(-\frac{1}{2}(\boldsymbol{x_i}-\boldsymbol{\mu_c})^T\boldsymbol{\Sigma_c^{-1}}(\boldsymbol{x_i}-\boldsymbol{\mu_c}))$$
   
   $r_{ic}$$x_i$ por lo tanto, sixiestá muy cerca de una c gaussiana, obtendrá un altovalor dericpara este gaussiano y valores relativamente bajos de lo contrario. M-Step_ Para cada grupo c: Calcule el peso totalm$\frac{Probability \ that \ x_i \ belongs \ to \ class \ c}{Probability \ of \ x_i \ over \ all \ classes}$$x_i$$r_{ic}$
   
   $\underline{M-Step}$
   
   $m_c$(en términos generales, la fracción de puntos asignados al grupo c) y actualice , μ c y Σ c utilizando r i c con: m c = Σ i r i c π c = m $\pi_c$$\mu_c$$\Sigma_c$$r_{ic}$
   
   $$m_c \ = \ \Sigma_i r_ic$$
   μc= $$\pi_c \ = \ \frac{m_c}{m}$$
   Σc= $$\boldsymbol{\mu_c} \ = \ \frac{1}{m_c}\Sigma_i r_{ic} \boldsymbol{x_i}$$
   Tenga en cuenta que debe usar los medios actualizados en esta última fórmula. Repita iterativamente los pasos E y M hasta que la función de probabilidad logarítmica de nuestro modelo converja donde la probabilidad logarítmica se calcula con: lnp(X|π,μ,Σ)=Σ N i = 1 ln(Σ $$\boldsymbol{\Sigma_c} \ = \ \frac{1}{m_c}\Sigma_i r_{ic}(\boldsymbol{x_i}-\boldsymbol{\mu_c})^T(\boldsymbol{x_i}-\boldsymbol{\mu_c})$$
   
   
   
   
   $$ln \ p(\boldsymbol{X} \ | \ \boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) \ = \ \Sigma_{i=1}^N \ ln(\Sigma_{k=1}^K \pi_k N(\boldsymbol{x_i} \ | \ \boldsymbol{\mu_k},\boldsymbol{\Sigma_k}))$$

$X$$AX = XA = I$

$$\begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}$$

$A$$X$$I$$\boldsymbol{\Sigma_c^{-1}}$
$\boldsymbol{0}$matriz de covarianza anterior si el gaussiano multivariado cae en un punto durante la iteración entre el paso E y M. Esto podría suceder si tenemos, por ejemplo, un conjunto de datos en el que queremos ajustar 3 gaussianos, pero que en realidad consiste solo en dos clases (grupos) de tal manera que, en términos generales, dos de estos tres gaussianos atrapan su propio grupo mientras que el último gaussiano solo lo maneja para atrapar un solo punto en el que se asienta. Veremos cómo se ve esto a continuación. Pero paso a paso: suponga que tiene un conjunto de datos bidimensional que consta de dos grupos, pero no lo sabe y desea ajustar tres modelos gaussianos, es decir c = 3. Inicializa sus parámetros en el paso E y traza los gaussianos en la parte superior de sus datos que se ve algo. como (tal vez pueda ver los dos grupos relativamente dispersos en la parte inferior izquierda y superior derecha): $\mu_c$$\pi_c$

$r_{ic}$$cov$$r_{ic}$
$$r_{ic} = \frac{\pi_c N(\boldsymbol{x_i} \ | \ \boldsymbol{\mu_c},\boldsymbol{\Sigma_c})}{\Sigma_{k=1}^K \pi_k N(\boldsymbol{x_i \ | \ \boldsymbol{\mu_k},\boldsymbol{\Sigma_k}})}$$ $r_{ic}$$r_{ic}$$x_i$$x_i$$x_i$$r_{ic}$$x_i$$r_{ic}$$r_{ic}$ $$\boldsymbol{\Sigma_c} \ = \ \Sigma_i r_{ic}(\boldsymbol{x_i}-\boldsymbol{\mu_c})^T(\boldsymbol{x_i}-\boldsymbol{\mu_c})$$ $r_{ic}$$x_i$$(\boldsymbol{x_i}-\boldsymbol{\mu_c})$$\boldsymbol{\mu_c}$$x_i$$j$$\boldsymbol{\mu_j}$$\boldsymbol{\mu_j} = \boldsymbol{x_n}$$r_{ic}$

$$\begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}$$

$\boldsymbol{0}$$\boldsymbol{0}$matriz. Esto se hace agregando un valor muy pequeño (en la Mezcla Gaussiana de sklearn, este valor se establece en 1e-6) a la digonal de la matriz de covarianza. También hay otras formas de prevenir la singularidad, como darse cuenta cuando un gaussiano colapsa y establecer su matriz de media y / o covarianza en un nuevo valor arbitrariamente alto. Esta regularización de covarianza también se implementa en el siguiente código con el que obtiene los resultados descritos. Quizás tenga que ejecutar el código varias veces para obtener una matriz de covarianza singular desde entonces, como se dijo. Esto no debe suceder cada vez, sino que también depende de la configuración inicial de los gaussianos.

import matplotlib.pyplot as plt
from matplotlib import style
style.use('fivethirtyeight')
from sklearn.datasets.samples_generator import make_blobs
import numpy as np
from scipy.stats import multivariate_normal


# 0. Create dataset
X,Y = make_blobs(cluster_std=2.5,random_state=20,n_samples=500,centers=3)

# Stratch dataset to get ellipsoid data
X = np.dot(X,np.random.RandomState(0).randn(2,2))


class EMM:

    def __init__(self,X,number_of_sources,iterations):
        self.iterations = iterations
        self.number_of_sources = number_of_sources
        self.X = X
        self.mu = None
        self.pi = None
        self.cov = None
        self.XY = None



    # Define a function which runs for i iterations:
    def run(self):
        self.reg_cov = 1e-6*np.identity(len(self.X[0]))
        x,y = np.meshgrid(np.sort(self.X[:,0]),np.sort(self.X[:,1]))
        self.XY = np.array([x.flatten(),y.flatten()]).T


        # 1. Set the initial mu, covariance and pi values
        self.mu = np.random.randint(min(self.X[:,0]),max(self.X[:,0]),size=(self.number_of_sources,len(self.X[0]))) # This is a nxm matrix since we assume n sources (n Gaussians) where each has m dimensions
        self.cov = np.zeros((self.number_of_sources,len(X[0]),len(X[0]))) # We need a nxmxm covariance matrix for each source since we have m features --> We create symmetric covariance matrices with ones on the digonal
        for dim in range(len(self.cov)):
            np.fill_diagonal(self.cov[dim],5)


        self.pi = np.ones(self.number_of_sources)/self.number_of_sources # Are "Fractions"
        log_likelihoods = [] # In this list we store the log likehoods per iteration and plot them in the end to check if
                             # if we have converged

        # Plot the initial state    
        fig = plt.figure(figsize=(10,10))
        ax0 = fig.add_subplot(111)
        ax0.scatter(self.X[:,0],self.X[:,1])
        for m,c in zip(self.mu,self.cov):
            c += self.reg_cov
            multi_normal = multivariate_normal(mean=m,cov=c)
            ax0.contour(np.sort(self.X[:,0]),np.sort(self.X[:,1]),multi_normal.pdf(self.XY).reshape(len(self.X),len(self.X)),colors='black',alpha=0.3)
            ax0.scatter(m[0],m[1],c='grey',zorder=10,s=100)


        mu = []
        cov = []
        R = []


        for i in range(self.iterations):               

            mu.append(self.mu)
            cov.append(self.cov)


            # E Step
            r_ic = np.zeros((len(self.X),len(self.cov)))

            for m,co,p,r in zip(self.mu,self.cov,self.pi,range(len(r_ic[0]))):
                co+=self.reg_cov
                mn = multivariate_normal(mean=m,cov=co)
                r_ic[:,r] = p*mn.pdf(self.X)/np.sum([pi_c*multivariate_normal(mean=mu_c,cov=cov_c).pdf(X) for pi_c,mu_c,cov_c in zip(self.pi,self.mu,self.cov+self.reg_cov)],axis=0)
            R.append(r_ic)

            # M Step

            # Calculate the new mean vector and new covariance matrices, based on the probable membership of the single x_i to classes c --> r_ic
            self.mu = []
            self.cov = []
            self.pi = []
            log_likelihood = []

            for c in range(len(r_ic[0])):
                m_c = np.sum(r_ic[:,c],axis=0)
                mu_c = (1/m_c)*np.sum(self.X*r_ic[:,c].reshape(len(self.X),1),axis=0)
                self.mu.append(mu_c)

                # Calculate the covariance matrix per source based on the new mean
                self.cov.append(((1/m_c)*np.dot((np.array(r_ic[:,c]).reshape(len(self.X),1)*(self.X-mu_c)).T,(self.X-mu_c)))+self.reg_cov)
                # Calculate pi_new which is the "fraction of points" respectively the fraction of the probability assigned to each source 
                self.pi.append(m_c/np.sum(r_ic)) 



            # Log likelihood
            log_likelihoods.append(np.log(np.sum([k*multivariate_normal(self.mu[i],self.cov[j]).pdf(X) for k,i,j in zip(self.pi,range(len(self.mu)),range(len(self.cov)))])))



        fig2 = plt.figure(figsize=(10,10))
        ax1 = fig2.add_subplot(111) 
        ax1.plot(range(0,self.iterations,1),log_likelihoods)
        #plt.show()
        print(mu[-1])
        print(cov[-1])
        for r in np.array(R[-1]):
            print(r)
        print(X)

    def predict(self):
        # PLot the point onto the fittet gaussians
        fig3 = plt.figure(figsize=(10,10))
        ax2 = fig3.add_subplot(111)
        ax2.scatter(self.X[:,0],self.X[:,1])
        for m,c in zip(self.mu,self.cov):
            multi_normal = multivariate_normal(mean=m,cov=c)
            ax2.contour(np.sort(self.X[:,0]),np.sort(self.X[:,1]),multi_normal.pdf(self.XY).reshape(len(self.X),len(self.X)),colors='black',alpha=0.3)




EMM = EMM(X,3,100)     
EMM.run()
EMM.predict()

En mi humilde opinión, todas las respuestas pierden un hecho fundamental. Si uno mira el espacio de parámetros para un modelo de mezcla gaussiana, este espacio es singular a lo largo del subespacio donde hay menos del número total de componentes en la mezcla. Eso significa que los derivados son automáticamente cero y, por lo general, todo el subespacio se mostrará como un archivo. Más filosóficamente, el subespacio de covarianzas de rango inferior al completo es el límite del espacio de parámetros y uno siempre debe sospechar cuando el mle ocurre en el límite; generalmente indica que hay un espacio de parámetros más grande al acecho en el que uno puede encontrar el 'real' mle. Hay un libro llamado "Estadísticas algebraicas" por Drton, Sturmfeld y Sullivant. Este tema se discute en ese libro con cierto detalle. Si eres realmente curioso, deberías mirar eso.

Para un solo gaussiano, la media posiblemente sea igual a uno de los puntos de datos ($x_n$ por ejemplo) y luego está el siguiente término en la función de probabilidad:

$$\begin{equation} {\cal N}(x_n|x_n,\sigma_j 1\!\!1)\rightarrow \lim_{\sigma_j\rightarrow x_n}\frac{1}{(2\pi)^{1/2}\sigma_j} \exp \left( -\frac{1}{\sigma_j}|x_n-\sigma_j|^2 \right)= \frac{1}{(2\pi)^{1/2}\sigma_j} \end{equation}$$ El límite $\sigma_j\rightarrow 0$ ahora es claramente divergente ya que el argumento de lo exponencial se desvanece.

Sin embargo para un punto de datos $x_m$ diferente de la media $\sigma_j$, tendremos

$$\begin{equation} {\cal N}(x_m|x_m,\sigma_j 1\!\!1)= \frac{1}{(2\pi)^{1/2}\sigma_j} \exp \left( -\frac{1}{\sigma_j}|x_m-\sigma_j|^2 \right) \end{equation}$$ y ahora el argumento de la exponencial diverge (y es negativo) en el límite $\sigma_j\rightarrow 0$. Como resultado, el producto de estos dos términos en la función de probabilidad desaparecerá.