Questões de singularidade no modelo de mistura gaussiano

No capítulo 9 do livro Reconhecimento de padrões e aprendizado de máquina, há uma parte sobre o modelo de mistura gaussiano:

Para ser sincero, não entendo por que isso criaria uma singularidade. Alguém pode me explicar isso? Sinto muito, mas sou apenas um graduado e um novato em aprendizado de máquina, então minha pergunta pode parecer um pouco boba, mas por favor me ajude. Muito obrigado

Se quisermos ajustar um gaussiano a um único ponto de dados usando a máxima probabilidade, obteremos um gaussiano muito espetado que "entra em colapso" nesse ponto. A variância é zero quando há apenas um ponto, que no caso gaussiano multi-variável, leva a uma matriz de covariância singular, por isso é chamado de problema de singularidade.

Quando a variância chega a zero, a probabilidade do componente Gaussiano (fórmula 9.15) vai para o infinito e o modelo fica sobreajustado. Isso não ocorre quando ajustamos apenas um Gaussiano a um número de pontos, pois a variação não pode ser zero. Mas isso pode acontecer quando temos uma mistura de gaussianos, conforme ilustrado na mesma página do PRML.

Atualização :
O livro sugere dois métodos para abordar o problema da singularidade, que são

1) redefinir a média e a variância quando a singularidade ocorre

2) usando MAP em vez de MLE adicionando um prior.

Lembre-se de que esse problema não surgiu no caso de uma única distribuição gaussiana. Para entender a diferença, observe que, se um único Gaussiano entrar em colapso em um ponto de dados, ele contribuirá com fatores multiplicativos para a função de probabilidade decorrente dos outros pontos de dados e esses fatores chegarão a zero exponencialmente rápido, fornecendo uma probabilidade geral de zero em vez de zero. do que infinito.

Também estou meio confuso com esta parte, e aqui está a minha interpretação. Tome 1D case para simplificar.

Quando um único "colapsos" Gaussian sobre um ponto de dados , isto é, μ = x i , a probabilidade global torna-se:$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}} )$$

Você vê como , o termo à esquerda p ( x i ) → ∞ , que é como o caso patológico no GMM, mas o termo à direita, que é a probabilidade de outros pontos de dados p ( x ∖ i ) , ainda contém termos como e - ( x n - μ ) $\sigma \to 0$$p(x_i) \to \infty$$p(\mathbf{x}\setminus{i})$ que→0 éexponencialmente rápido comoσ→0, então o efeito geral sobre a probabilidade é que ele vá ao zero.$e^{-\frac{(x_n-\mu)^2}{2\sigma^2}}$$\to 0$$\sigma \to 0$

O ponto principal aqui é que, ao ajustar um único Gaussiano, todos os pontos de dados precisam compartilhar um conjunto de parâmetros , diferente do caso de mistura em que um componente pode "focar" em um ponto de dados sem penalizar a probabilidade geral de dados .$\mu, \sigma$

Esta resposta fornecerá uma visão do que está acontecendo que leva a uma matriz de covariância singular durante a adaptação de um GMM a um conjunto de dados, por que isso está acontecendo e o que podemos fazer para evitar isso.

Portanto, é melhor começar recapitulando as etapas durante a adaptação de um Modelo de Mistura Gaussiana a um conjunto de dados.

0. Decida quantas fontes / clusters (c) você deseja ajustar aos seus dados
1. Inicialize os parâmetros como , covariância Σ c e fração_per_classe π c por cluster c $\mu_c$$\Sigma_c$$\pi_c$

$\underline{E-Step}$

2. Calcule para cada ponto de dados a probabilidade r i c que o ponto de dados x i pertence ao cluster c com: r i c = π c N ( x i | μ c , Σ c$x_i$$r_{ic}$$x_i$
   
   queN(x|μ,Σ)descreve o gaussiano multivariado com: 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 fornece, para cada ponto de dadosxi,a 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 conseguinte, sexié muito próximo de um c de Gauss, que terá uma elevadaricvalor para esses valores gaussianos e relativamente baixos, caso contrário. M-Step_ Para cada cluster c: Calcule o 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$(falando livremente a fracção de pontos atribuídos ao aglomerado c) e atualizar , μ c , e Σ c utilizando r i c com: 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}$$
   Lembre-se de que você deve usar os meios atualizados nesta última fórmula. Repita iterativamente os passos E e M até que a função de probabilidade logarítmica do nosso modelo converja para a qual a probabilidade logarítmica é calculada com: 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 covariância acima se o Gaussiano multivariado cair em um ponto durante a iteração entre as etapas E e M. Isso pode acontecer se, por exemplo, tivermos um conjunto de dados no qual queremos ajustar 3 gaussianos, mas que na verdade consiste apenas de duas classes (clusters), de modo que, falando livremente, dois desses três gaussianos capturem seu próprio cluster enquanto o último gaussiano o gerencia apenas para pegar um único ponto em que está assentado. Vamos ver como isso se parece abaixo. Mas passo a passo: suponha que você tenha um conjunto de dados bidimensional que consiste em dois clusters, mas você não sabe disso e deseja ajustar três modelos gaussianos, ou seja, c = 3. Você inicializa seus parâmetros na etapa E e na plotagem os gaussianos em cima dos seus dados, que parecem smth. como (talvez você possa ver os dois grupos relativamente dispersos no canto inferior esquerdo e no canto superior direito): $\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. Isso é feito adicionando um valor muito pequeno (no GaussianMixture do sklearn, esse valor é definido como 1e-6) ao digonal da matriz de covariância. Também existem outras maneiras de impedir a singularidade, como perceber quando um gaussiano entra em colapso e definir sua matriz de média e / ou covariância para um novo valor arbitrariamente alto. Essa regularização de covariância também é implementada no código abaixo com o qual você obtém os resultados descritos. Talvez você precise executar o código várias vezes para obter uma matriz de covariância singular, pois, como dito. isso não deve acontecer a cada vez, mas também depende da configuração inicial dos 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()

Imho, todas as respostas perdem um fato fundamental. Se observarmos o espaço de parâmetro para um modelo de mistura gaussiana, esse espaço é singular ao longo do subespaço, onde há menos do que o número total de componentes na mistura. Isso significa que as derivadas são automaticamente zero e, geralmente, todo o subespaço será exibido como um mle. Mais filosoficamente, o subespaço de covariâncias inferiores à classificação completa é o limite do espaço de parâmetro e deve-se sempre desconfiar quando a mle ocorre no limite - geralmente indica que há um espaço de parâmetro maior à espreita no qual se pode encontrar o 'real' mle. Há um livro chamado "Algebraic Statistics" de Drton, Sturmfeld e Sullivant. Essa questão é discutida nesse livro com mais detalhes. Se você está realmente curioso, você deveria olhar para isso.

$x_n$

$$\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}$$ $\sigma_j\rightarrow 0$

$x_m$$\sigma_j$

$$\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}$$ e agora o argumento do exponencial diverge (e é negativo) no limite $\sigma_j\rightarrow 0$. Como resultado, o produto desses dois termos na função de probabilidade desaparecerá.