Normalmente utilizamos uma determinada distribuição para modelar dados. Escolhida uma distribuição para a modelagem, precisamos estimar seus parâmetros por algum método de estimação. Por exemplo, se estamos trabalhando com a distribuição normal, é de interesse estimar os parâmetros $\mu$ e $\sigma^2$, que são a média e a variância da distribuição. Note que, nesse caso, supomos que os dados são originados de uma distribuição apenas.
Mas e se os dados vierem de duas ou mais distribuições? Confuso? Imagine clusters (ou agrupamentos), como na figura abaixo. Percebemos que uma parte dos dados parecem ser originados de uma distribuição e a outra parte, de outra distribuição. Os dados observados (representados pelo histograma) contemplam os dados das duas distribuições mas, na prática, não sabemos quais são (ou quantas são). Dessa forma, a modelagem de mistura de distribuições se torna bastante útil.
Nesse post, ensino como trabalhar com mistura de distribuições Gaussianas, ou seja, quando os dados pertencem a mais de uma distribuição normal (médias e/ou variâncias distintas). Para a estimação dos parâmetros, utilizaremos o algoritmo EM (já explicado aqui no blog).
Gaussian Mixture Models
Para começarmos, vamos supor que nossos dados sejam provenientes de duas distribuições Gaussianas com parâmetros desconhecidos. Seja $X$ nossa variável de interesse. A função de densidade do modelo é dada por
$$\begin{eqnarray}\label{eq:model} f(x|\theta) = p_1f(x|\theta_1) + p_2f(x|\theta_2) \end{eqnarray}$$
em que $\theta_1 = (\mu_1,\sigma^2_1)$ e $\theta_2 = (\mu_2,\sigma^2_2)$ são os vetores de parâmetros das distribuições normais 1 e 2, respectivamente, e $p_1$ e $p_2$ são os pesos relacionados a ambas distribuições. Podemos interpretar os pesos como as probabilidades de uma observação ser proveniente de uma determinada distribuição, ou seja, $p_1$ é a probabilidade dessa observação pertencer à distribuição 1 e $p_2$, a probabilidade dela pertencer à distribuição 2. O vetor de parâmetros do modelo é $\theta = (\theta_1,\theta_2,p_1)$. Por que apenas $p_1$? Porque $p_2$ é completamente determinada por $p_1$, isto é, $p_2 = 1 - p_1$.
Seja $X_1, X_2, ..., X_n$ uma amostra aleatória da distribuição em ($\ref{eq:model}$). A função de log-verossimilhança do modelo é dada por
$$\begin{eqnarray} l(\theta) &=& \log \left[ \prod_{i=1}^n f(x_i|\theta) \right] = \sum_{i=1}^n \log f(x_i|\theta) = \sum_{i=1}^n \log \big[ p_1 f(x_i|\theta_1) + p_2 f(x_2|\theta_2) \big]\end{eqnarray}$$
Se prosseguirmos para encontrar os estimadores de máxima verossimilhança (EMV) a partir da função acima, veremos que não há forma fechada para os estimadores. Logo, outra abordagem é necessária.
Modelo aumentado
Para isso, imagine que para cada observação $x_i$ conhecemos a qual distribuição ela pertence. Isso é equivalente a dizer que conhecemos uma variável $Y_i$ que indica se a observação $x_i$ pertence a distribuição $C$, $C = 1,2$. A função de densidade de $Y_i$ é dada por
$$f(y_i) = \prod_{y_i \in C} p_{y_i}^{1 \{y_i = c\}}, \quad C = \{1,2\}$$
em que $1\{\}$ é a função indicadora. Dado o valor de $y_i$, a função de densidade de $X_i$ será
$$f(x_i|y_i) = \prod_{y_i \in C} f(x_i|\theta_{y_i})^{1\{y_i = c\}}$$
Na prática não observamos os valores de $Y_i$. Dessa forma, $Y_i$ é uma variável latente e trabalhamos com o modelo aumentado $(X,Y)$. A função de densidade conjunta de $(X_i,Y_i)$ é dada por
\begin{eqnarray}
f(x_i,y_i) &=& f(x_i|y_i) f(y_i) = \prod_{y_i \in C} \left( p_{y_i} f(x_i|\theta_{y_i}) \right)^{1\{y_i = c\}} \nonumber
\end{eqnarray}
A função de verossimilhança do modelo aumentado, dada uma amostra aleatória $(X,Y) = ((X_1,Y_1),...,(X_n,Y_n))$, será
\begin{eqnarray}
L(\theta) &=& \prod_{i=1}^n f(x_i|y_i) f(y_i) = \prod_{i=1}^n \left[ \prod_{y_i \in C} \left( p_{y_i} f(x_i|\theta_{y_i}) \right)^{1\{y_i = c\}} \right] \nonumber
\end{eqnarray}
e a função de log-verossimilhança é dada por
\begin{eqnarray}
l(\theta) &=& \sum_{i=1}^n \sum_{y_i \in C} {1\{y_i = c\}} \left[ \log (p_{y_i}) + \log (f(x_i|\theta_{y_i})) \right] \nonumber
\end{eqnarray}
Vamos refletir um pouco sobre o que nós fizemos.
- Primeiro modelamos o problema e chegamos na função de verossimilhança em (\ref{eq:model}).
- Percebemos que não há como estimar os parâmetros do modelo dessa forma e outra abordagem se tornou necessária.
- Consideramos uma variável $Y$ para indicar em que grupo (cluster) a variável $X$ pertence.
- Escrevemos o modelo completo para $(X,Y)$.
Observe que $X$ é variável observada, enquanto que $Y$ é não observada (latente). Dessa forma, precisamos de um método para estimar os parâmetros do modelo nesse contexto. Como falamos no começo do post, utilizaremos o algoritmo EM. Mas por quê? Se você acompanha o blog, já deve saber a resposta.
Temos um modelo com variáveis observadas e variáveis não observadas. O algoritmo EM é excelente no contexto de dados faltantes, em que a função de verossimilhança completa é mais simples de se trabalhar do que a função de verossimilhança dos dados observados. Dessa forma, vamos prosseguir para a estimação dos parâmetros.
Estimação dos parâmetros - Algoritmo EM
Passo E
A função $Q$ será dada por
\begin{eqnarray}
Q(\theta|\theta^{(k)}) &=& \sum_{i=1}^n \sum_{y_i \in C} E[{1\{y_i = c\}}|x,\theta_{(k)}] \left[ \log (p_{y_i}) + \log (f(x_i|\theta_{y_i})) \right] \nonumber
\end{eqnarray}
Para facilitar a notação, escrevemos $\gamma_{ic} = E[{1\{y_i = c\}}|x,\theta_{(k)}] $. Note que $\gamma_{ic}$ é a probabilidade de $Y_i$ ser igual a $c$ (ou seja, indicar que $x_i$ pertence a distribuição $c$) dado o valor de $x_i$. Ou seja,
\begin{eqnarray} \gamma_{ic} = P(Y_i = c|x_i) = f(y_i = c|x_i) \nonumber \end{eqnarray}
Para calcularmos o valor de $\gamma_{ic}$, precisamos encontrar a distribuição de $Y_i|X_i$. Temos então
\begin{eqnarray}
f(y_i|x_i) &=& \dfrac{f(x_i|y_i)f(y_i)}{f(x_i)} = \dfrac{\prod_{y_i \in C} \left[p_{y_i}f(x_i|\theta_{y_i}) \right] ^{1\{y_i = c\}}}{p_{1}f(x_i|\theta_{1})+p_{2}f(x_i|\theta_{2})} \nonumber
\end{eqnarray}
Dessa forma,
\begin{eqnarray} \gamma_{ic} &=& P(Y_i = c|x_i) = \dfrac{p_{c}f(x_i|\theta_{c}) }{p_{1}f(x_i|\theta_{1})+p_{2}f(x_i|\theta_{2})} \nonumber \end{eqnarray}
Para encontrar os estimadores, precisamos derivar a função $Q()$ em relação aos parâmetros de interesse. Vamos realizar primeiro a estimação de $\mu_1$:
\begin{eqnarray}
\dfrac{d}{d\mu_1}Q() &=&\sum_{i=1}^{n} \gamma_{i1} \dfrac{d}{d\mu_1} \log f(x_i|\theta_1) = \sum_{i=1}{n} \gamma_{i1} \dfrac{(x_i-\mu_1)}{\sigma^2_1} \nonumber
\end{eqnarray}
Igualando $(d/d\mu_1)Q()$ a zero, obtemos $\hat{\mu}_1 = \dfrac{\sum_{i=1}^n (\gamma_{i1} x_i)}{\sum_{i=1}^n \gamma_{i1}}$.
De forma análoga, $\hat{\mu}_2 = \dfrac{\sum_{i=1}^n (\gamma_{i2} x_i)}{\sum_{i=1}^n \gamma_{i2}}$
Em relação à variância da distribuição 1, $\sigma^2_1$, temos
\begin{eqnarray}
\dfrac{d}{d\sigma^2_1}Q() &=&\sum_{i=1}{n} \gamma_{i1} \left(\dfrac{-1}{2\sigma^2_1} + \dfrac{1}{2\sigma^4_1}(x_i-\mu_1)^2 \right)\nonumber \\
\end{eqnarray}
Igualando $(d/d\sigma^2_1)Q()$ a zero, obtemos $\hat{\sigma}^2_1 = \dfrac{\sum_{i=1}^n \gamma_{i1}(x_i-\mu_1)^2 }{\sum_{i=1}^n \gamma_{i1}}$.
De forma análoga, $\hat{\sigma}^2_2 = \dfrac{\sum_{i=1}^n \gamma_{i2}(x_i-\mu_2)^2 }{\sum_{i=1}^n \gamma_{i2}}$.
Em relação aos pesos $p_1$ e $p_2$, temos que lembrar que
$$p_1 + p_2 = 1 \rightarrow p_2 = 1 - p_1$$
Dessa forma, vamos encontrar o estimador de $p_1$.
\begin{eqnarray}
\dfrac{d}{d p_1}Q() &=&\sum_{i=1}^{n} \left( \gamma_{i1} \dfrac{1}{p_1} - \gamma_{i2} \dfrac{1}{1-p_1} \right) \nonumber
\end{eqnarray}
Igualando $(d/dp_1)Q()$ a zero, obtemos $\hat{p}_1 = \dfrac{\sum_{i=1}^n \gamma_{i1} }{\sum_{i=1}^n (\gamma_{i1}+\gamma_{i2})}$.
Logo, o estimador para $p_2$ é $\hat{p}_2 = 1 - \hat{p}_1$.
Dessa forma, encontramos uma metodologia para estimar os parâmetros da distribuição de mistura de normais. A implementação deste exemplo será feita no próximo post.
Espero que tenha gostado da aula.
Até a próxima aula!
