Hoje vou falar um pouco sobre um algoritmo de estimação bastante conhecido e utilizado quando temos dados faltantes ou variáveis latentes no modelo que estamos trabalhando: o algoritmo EM.
O algoritmo EM foi proposto por Dempster et al. (1977) e é um algoritmo iterativo de estimação de máxima verossimilhança. Este método envolve duas etapas: uma de Esperança (Expectation) e outra de Maximização (Maximization), daí seu nome. A ideia é, no primeiro passo, calcular os valores esperados das variáveis latentes e, no segundo passo, maximizar a função de interesse em relação aos parâmetros do modelo. Essa função é chamada de função Q, como veremos mais adiante.
Como dito, estamos em uma situação em que há dados ausentes, seja por haver missing data ou por modelarmos variáveis latentes para simplificar nosso modelo (o exemplo que darei reflete o segundo caso). Sem perda de generalidade, vamos considerar o caso em que há apenas uma variável latente. Seja $y$ o vetor de valores observados da variável resposta e $u$ a variável latente (ou variável que representa os dados faltantes). Trabalhamos com a função de log-verossimilhança completa dos dados, dada por
$l_c(\theta) = \log f(y,u) = \log f(y|u) + \log(u).$
A função acima seria realmente completa se tivéssemos observado os valores de $u$, mas não temos tais valores. Dessa forma, o primeiro passo do algoritmo EM - passo E - é calcular o valor esperado da função de log-verossimilhança completa dados os valores de $y$ e dos parâmetros de interesse. Para isso, definimos a função $Q$ como
$Q(\theta | \theta^{(k)},y) = E[\log f(y,u)|\theta^{(k)},y].$
Como dito, o algoritmo EM é iterativo, e a função acima calcula a esperança da função de log-verossimilhança completa dados os valores de $y$ e os valores de $\theta$ na $k$-ésima iteração. Ou seja, precisamos dos valores dos parâmetros de interesse para realizar esse passo. Para ficar claro: o algoritmo é iniciado com valores iniciais para os parâmetros.
O segundo passo do algoritmo EM - passo M - é realizado maximizando a função $Q$ em relação aos parâmetros de interesse, isto é,
$\hat{\theta}^{(k+1)} = \text{argmax} Q(\theta | \theta^{(k)},y) .$
Os dois passos são repetidos até que haja convergência na estimação dos parâmetros ou um número máximo de iterações seja atingida.
É importante dizer que os estimadores EM convergem para os estimadores de máxima verossimilhança (EMV). De fato, utilizamos o algoritmo EM quando não há forma fechada para os EMVs do modelo ou quando facilitamos o problema trabalhando com o modelo aumentado (verossimilhança completa).
Isso tudo é bem teórico e acredite quando falo que estou resumindo bastante as coisas aqui. Para tudo que fazemos há fundamentos mas esse blog é voltado para a parte prática, então vou contar que você pesquisará mais sobre o assunto se não é familiarizado com o algoritmo EM. Aliás, não queremos ser apertadores de botão, não é mesmo?
Para entender o método, vamos ver um exemplo simples e didático.
Exemplo: estimação dos parâmetros da distribuição t-Student (com graus de liberdade fixos)
Suponha que temos uma amostra aleatória de tamanho $n$ de uma população com distribuição t-Student com $\nu$ graus de liberdade e suponha que nosso interesse é estimar os parâmetros de média $\mu$ e de escala $\sigma^2$.
Uma opção seria escrever a função de log-verossimilhança da distribuição t e encontrar os estimadores de máxima verossimilhança (EMV). Mas queremos utilizar o algoritmo EM, então vamos escrever a variável $Y$, que segue uma distribuição t, na seguinte forma hierárquica:
Note que, marginalmente (em relação a $U$), $Y$ possui uma distribuição t-Student com $\nu$ graus de liberdade.
Bom, agora vamos construir o algoritmo para nosso exemplo. A função de log-verossimilhança completa é dada por
$\begin{eqnarray} l_c(\theta) &\propto& \frac{n}{2} \log(u) - \frac{n}{2} \log(\sigma^2) - \frac{u}{2\sigma^2} \sum_{i=1}^{n} (y_i-\mu)^2 - \frac{\nu}{2} u + \left(\frac{\nu}{2}-1\right) \log(u). \nonumber \end{eqnarray}$
Escrevi "proporcional" por não precisarmos das constantes (em relação aos parâmetros). A função $Q$ será dada por
$\begin{eqnarray} Q(\theta | \theta^{(k)},y) &\propto& \frac{n}{2} E[\log(u)| \theta^{(k)},y] - \frac{n}{2} \log(\sigma^2) - \frac{E[u| \theta^{(k)},y]}{2\sigma^2} \sum_{i=1}^{n} (y_i-\mu)^2 \nonumber\\ && - \frac{\nu}{2} E[u| \theta^{(k)},y] + \left(\frac{\nu}{2}-1\right) E[\log(u)| \theta^{(k)},y]. \nonumber \end{eqnarray}$
O segundo passo é maximizar a função $Q$ em relação aos parâmetros de interesse, no caso, $\mu$ e $\sigma^2$. Para isso, derivamos a função $Q$ em relação ao parâmetro e igualamos a zero para encontrar seu estimador. Note, na função acima, que há termos que são constantes em relação aos parâmetros e, após a derivação, serão anulados. Dessa forma, podemos simplificar ainda mais a função $Q$ mantendo apenas os termos que apresentam $\mu$ e $\sigma^2$, ou seja,
$\begin{eqnarray} Q(\theta | \theta^{(k)},y) &\propto& - \frac{n}{2} \log(\sigma^2) - \frac{E[u| \theta^{(k)},y]}{2\sigma^2} \sum_{i=1}^{n} (y_i-\mu)^2. \nonumber \end{eqnarray}$
Isso nos ajuda (e muito), pois não precisamos calcular a distribuição de $\log(u)| \theta^{(k)},y$ e diminuímos nosso trabalho. Para simplificar a notação, vamos escrever $E[u| \theta^{(k)},y] = E_{U|Y}$.
Agora sim vamos para o segundo passo do algoritmo. Primeiro, vamos encontrar o estimador para a média:
$\begin{eqnarray} \dfrac{d}{d\mu}Q() = \dfrac{E_{U|Y}}{\sigma^2} (\sum_{i=1}^{n} y_i- n\mu). \nonumber \end{eqnarray}$
Igualando a derivada acima a zero, obtemos
$\begin{eqnarray} \hat{\mu} = \dfrac{\sum_{i=1}^{n} y_i}{n}. \nonumber \end{eqnarray}$
Agora vamos encontrar o estimador para $\sigma^2$:
$\begin{eqnarray} \dfrac{d}{d\sigma^2}Q() = -\dfrac{n}{2\sigma^2} + \dfrac{E_{U|Y}}{2\sigma^4} \sum_{i=1}^{n} (y_i- \mu)^2. \nonumber \end{eqnarray}$
o que nos leva ao seguinte estimador:
$\begin{eqnarray} \hat{\sigma}^2 = E_{U|Y} \dfrac{\sum_{i=1}^{n} (y_i- \mu)^2}{n}. \nonumber \end{eqnarray}$
Você deve ter notado falta de alguma coisa: como calculamos $E_{U|Y}$? Para isso, vamos descobrir a distribuição de $U|Y$ e, em seguida, calcular seu valor esperado. Note que
$\begin{eqnarray} f(u|y) &\propto& f(y|u)f(u) \nonumber \\ &\propto& (u)^{n/2} \exp \left\{ -\dfrac{u}{2\sigma^2} \sum_{i=1}^{n} (y_i- \mu)^2 \right\} \times \exp\left\{-\dfrac{\nu}{2}u\right\} u^{\nu/2-1} \nonumber \\ &\propto& (u)^{(n+\nu)/2 -1} \exp \left\{ -u\left(\dfrac{\nu}{2} + \dfrac{\sum_{i=1}^{n} (y_i- \mu)^2}{2\sigma^2}\right) \right\}. \nonumber \end{eqnarray}$
A última expressão corresponde ao núcleo de uma distribuição Gamma. Logo,
$$U|Y \sim Gamma \left( \dfrac{1}{2}(n+\nu), \dfrac{1}{2} (\nu + \sum_{i=1}^{n} (y_i- \mu)^2/\sigma^2) \right). $$
O valor esperado de uma distribuição Gamma($\alpha,\beta$) é $\alpha/\beta$. Então,
$$E_{U|Y} = \dfrac{(n-\nu)\sigma^2}{\nu\sigma^2 + \sum_{i=1}^{n} (y_i- \mu)^2}.$$
Agora temos o valor de $E_{U|Y}$ e os estimadores dos parâmetros de interesse. Estamos prontos para executar o algoritmo EM:
- Assuma valores iniciais $\hat{\mu}^{(0)}$ e $(\hat{\sigma}^2)^{(0)}$ e um número de iterações máxima.
- E: Dados os valores de $y$, $\hat{\mu}^{(k)}$ e $(\hat{\sigma}^2)^{(k)}$, calcule $E_{U|Y}$.
- M: Dados os valores de $y$ e $E_{U|Y}$, calcule $\hat{\mu}^{(k+1)}$ e $(\hat{\sigma}^2)^{(k+1)}$.
- Se houve convergência ou o número máximo de iterações foi atingida, retorne as estimativas dos parâmetros. Caso contrário, repita os passos 2-3.
Com esse exemplo é fácil ver como o algoritmo EM funciona. Gastamos um tempo para explicar e exemplificar o método para que fique claro seu funcionamento, pois iremos implementá-lo no próximo post. E para programar precisamos entender o que estamos fazendo, não é verdade?
Então espero você na próxima aula!
Referência citada
Dempster, A.P., Laird, N.M., and Rubin, D.B. (1977). Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society. Series B, 39(1):1-38.
