No post anterior apresentei o algoritmo EM e expliquei como o algoritmo funciona. Também apresentei um exemplo em que calculamos as quantidades necessárias para sua execução. Vamos relembrar esse exemplo.
Suponha que temos uma amostra aleatória de tamanho $n$ de uma população com distribuição t-Student com $\nu$ graus de liberdade. Nosso interesse é estimar os parâmetros de média $\mu$ e de escala $\sigma^2$. Seja $Y$ nossa variável de interesse. Escrevemos a distribuição de $Y$ na forma hierárquica
$$\begin{eqnarray}Y|U &\sim& N(\mu,u^{-1}\sigma^2), \nonumber \\ U &\sim& Gamma(\nu/2,\nu/2). \nonumber \end{eqnarray}$$
Vimos no post anterior que precisamos da quantidade $E[U|Y,\theta]$ e a calculamos:
$$E_{U|Y} = \dfrac{(n-\nu)\sigma^2}{\nu\sigma^2 + \sum_{i=1}^{n} (y_i- \mu)^2}.$$
Também encontramos os estimadores dos parâmetros ao maximizar a função $Q$, obtendo
$$\begin{eqnarray} \hat{\mu} &=& \dfrac{\sum_{i=1}^{n} y_i}{n}. \nonumber \\ \hat{\sigma}^2 &=& E_{U|Y} \dfrac{\sum_{i=1}^{n} (y_i- \mu)^2}{n}. \nonumber \end{eqnarray}$$
Lembrando que o primeiro passo - E - é calcular o valor esperado da variável latente dados os valores observados $y$ e os valores dos parâmetros na $k$-ésima iteração, e o segundo passo - M - é maximizar a função $Q$ dados os valores observados e o valor esperado da variável latente para obter os valores dos parâmetros da iteração $k+1$. Os passos são repetidos até que se observe convergência ou um número máximo de iterações seja atingido.
IMPORTANTE RELEMBRAR: os estimadores EM convergem para os estimadores de máxima verossimilhança (EMV).
Agora que relembramos as quantidades que serão necessárias, vamos implementar o algoritmo EM para nosso exemplo. Vamos começar gerando nossos dados:
n = 100;
nu = 5;
mu = 0;
sigma2 = 1;
set.seed(7);
u = rgamma(1,shape=n/2,rate=n/2);
y = rnorm(1000,mean = 0, sd = sqrt(sigma2)*(u^-1));
Estamos gerando $n=100$ observações de uma distribuição t-Student com 5 graus de liberdade, com média igual a 0 e escala igual a 1. Para visualizar os dados, computamos seu histograma, por exemplo.
hist(y,main='')
Qual a estimativa dos parâmetros pelos EMVs?
emv_mu = sum(y)/n
emv_sigma2 = sum((y-emv_mu)^2)/n
> emv_mu
[1] 0.09044438
> emv_sigma2
[1] 0.4832229
Observe que a estimativa da média é um pouco maior que zero e que, apesar de assumirmos que $\sigma^2 =1$, temos que $\hat{\sigma}^2 = 0,48$, menos que a metade do valor original. Isso é possível pois estamos trabalhando com uma amostra e essa foi a informação que a amostra nos forneceu.
Vamos partir para o algoritmo EM. O começo do código:
theta.curr = c(1,2);
max_ite = 100;
theta.trace = matrix(ncol=2,nrow=1);
theta.trace[1,] = theta.curr;
O que estamos fazendo?
- theta.curr armazenará os valores correntes do vetor de parâmetros. Iniciamos com $\mu=3$ e $\sigma^2=5$.
- Definimos o número máximo de iterações como 100.
- theta.trace é onde armazenaremos os valores estimados dos parâmetros em cada iteração do algoritmo.
- O primeiro valor armazenado em theta.trace são os valores iniciais.
Ao final do algoritmo, faremos um plot para visualizar as estimativas de cada parâmetro, observando visualmente se o algoritmo convergiu ou não. Essa é uma maneira de observar convergência, mas não a melhor. Discutiremos isso mais adiante.
A seguir, construímos o algoritmo propriamente dito. Os passos E e M são realizados até o número de iterações atingir o número máximo. Portanto
for(i in 1:max_ite){
#Passo 1 do EM - calcular E(U|Y)
alfa = 0.5*(n + nu)
beta = 0.5*(nu + sum((y-theta.curr[1])^2)/sigma2);
u.y = alfa/beta;
#Passo 2 do EM - Maximizar a função Q
# em relação aos parâmetros
mu.hat = sum(y)/n;
sigma2.hat = u.y*(sum((y-mu.hat)^2)/n);
theta.curr = c(mu.hat,sigma2.hat);
theta.trace = rbind(theta.trace,theta.curr);
}
O que estamos fazendo?
- Calculamos o valor de $E_{U|Y}$ - passo E.
- Calculamos os valores de $\hat{\mu}^{(k+1)}$ e $(\hat{\sigma}^2)^{(k+1)}$ - passo M.
- Atualizamos o vetor theta.curr com as novas estimativas dos parâmetros.
- Armazenamos essas estimativas em theta.trace.
- Esse bloco é repetido até que i seja igual a max_ite.
Se executarmos o código acima, obtemos:
> head(theta.trace)
[,1] [,2]
1.00000000 2.0000000
theta.curr 0.09044438 0.7194313
theta.curr 0.09044438 0.7030660
theta.curr 0.09044438 0.6881575
theta.curr 0.09044438 0.6745351
theta.curr 0.09044438 0.6620535
Observe que a estimativa da média não mudou após a primeira iteração. Isso ocorre porque seu estimador não depende do valor de $u$ nem da estimativa de $\sigma^2$. Em relação ao parâmetro de escala, perceba que temos uma diferença grande nas estimativas nas iterações iniciais e essa diferença diminui. É melhor visualizar isso em um gráfico de trajetória, construído com estimativa versus iteração.
par(mar=c(4,2,2,2),mfrow=c(1,2))
plot(theta.trace[,1],main=expression(mu),ylab='',
xlab='Iteração',type='l',ylim=c(-0.5,1),lwd=2)
abline(h=mu,col='blue',lwd=2)
plot(theta.trace[,2],main=expression(sigma^2),ylab='',
xlab='Iteração',type='l',lwd=2)
abline(h=sigma2,col='blue',lwd=2)
Nos gráficos acima, a linha azul representa o valor verdadeiro do parâmetro e a linha vermelha, o valor do EMV. Então perceba que as estimativas pelo algoritmo EM estão convergindo para os valores do EMV.
Simples, não? É claro que a programação do método não é tão simples para modelos mais complexos.
Você pode estar se perguntando: por que os estimadores do EM convergem para os valores do EMV e não para os valores verdadeiros? A resposta é bem teórica, mas vou resumir a ideia. O EM é uma ferramenta utilizada para estimar os parâmetros do modelo em um modelo aumentado onde há dados faltantes. A verossimilhança incompleta (ou observada) é obtida integrando a verossimilhança completa em relação às variáveis latentes. Dessa forma, os estimadores do EM convergem para os EMVs. Como os EMVs são não viciados para os valores dos parâmetros, então os estimadores do EM convergem para os parâmetros da população. Essa é só a ideia, se quiser aprofundar sobre o assunto aconselho estudar as provas matemáticas do método.
É importante executar o algoritmo EM com diferentes valores iniciais para os parâmetros para verificar se a convergência se dá na mesma região. Caso haja diferença da região de convergência a partir de valores iniciais distintos, pode significar que o algoritmo está convergindo para máximos locais e não para o máximo global. Ou seja, sempre faça essa verificação.
Para terminar, vamos organizar nosso código e colocá-lo em uma função. Dessa forma, não precisaremos executar o bloco de código do for toda vez que quisermos rodar o algoritmo. Além disso, vamos
- implementar um critério de parada, para que não seja necessário chegar ao máximo de iterações se o algoritmo convergir antes,
- plotar as trajetórias das estimativas se o usuário quiser.
O código da função fica como:
algoritmoEM_tStudent = function(par_ini,max_ite,y,nu,plot.trace=0){
theta.curr = c(par_ini[1],par_ini[2]);
n = length(y);
theta.trace = matrix(ncol=2,nrow=1);
theta.trace[1,] = theta.curr;
dif_parada = 10^-4;
for(i in 1:max_ite){
#Passo 1 do EM - calcular E(U|Y)
alfa = 0.5*(n + nu)
beta = 0.5*(nu + sum((y-theta.curr[1])^2)/theta.curr[2]);
u.y = alfa/beta;
#Passo 2 do EM - Maximizar a função Q
# em relação aos parâmetros
mu.hat = sum(y)/n;
sigma2.hat = u.y*(sum((y-mu.hat)^2)/n);
theta.curr = c(mu.hat,sigma2.hat);
theta.trace = rbind(theta.trace,theta.curr);
dif = abs(theta.trace[i+1,]-theta.trace[i,])
if(dif[1]<dif_parada & dif[2]<dif_parada) {
cat("Algoritmo convergiu com",i,"iterações.\n")
break;
}
}
if(plot.trace==1){
par(mfrow=c(1,2))
plot(theta.trace[,1],main=expression(mu),ylab='',
xlab='Iteração',type='l',lwd=2)
plot(theta.trace[,2],main=expression(sigma^2),ylab='',
xlab='Iteração',type='l',lwd=2)
}
return(theta.trace)
}
1: Para dizer que o algoritmo convergiu, verificamos se a diferença das estimativas de cada parâmetro da iteração atual para anterior é menor que 0,0001. Caso isso ocorra, paramos o algoritmo e evitamos iterações desnecessárias. Obviamente, se alterarmos o valor do critério de parada, alteramos os números de iterações para o mesmo problema.
2: Caso o usuário passe 1 no parâmetro plot.trace, após a execução do EM um gráfico de trajetória para cada parâmetro é construído.
3: Ao final do algoritmo, retornamos theta.trace, que contém todas as estimativas dos parâmetros.
Vamos testar nossa função.
> res = algoritmoEM_tStudent(c(1,2),100,y,nu,1)
Algoritmo convergiu com 91 iterações.
Armazenamos theta.trace (objeto retornado pela nossa função) em res. Os atributos passados são os mesmos do começo do post, então devemos ter os mesmos gráficos de trajetórias. Observe que o algoritmo parou com 91 iterações, antes do número máximo definido. Pense que, se tivéssemos atribuído um número máximo de 1000 iterações, esse critério de parada nos pouparia muito trabalho (e tempo e memória). Por fim, os gráficos plotados são
Perceba que são os mesmos gráficos que construímos anteriormente (sem as personalizações). Agora temos nossa função executando perfeitamente para o nosso problema.
Espero que tenha gostado da aula!
Até a próxima!
Nenhum comentário:
Postar um comentário