No post anterior, vimos o modelo de mistura de distribuições normais e fizemos um exemplo considerando duas distribuições Gaussianas. Vimos também que não há forma fechada para os EMVs dos parâmetros do modelo e recorremos ao modelo aumentado, incorporando uma variável latente para indicar de qual distribuição uma determinada observação é originada. Nesse contexto, utilizamos o algoritmo EM para estimar os parâmetros do modelo.
No post de hoje, vamos implementar o algoritmo de estimação para nosso exemplo utilizando o R. Como motivação, vamos simular dados de duas distribuições Gaussianas, a primeira com média -2 e a segunda com média 2, ambas com variância 1. A figura abaixo apresenta o histograma de 100 observações, sendo metade dos dados de cada distribuição.
set.seed(7)
x = c(rnorm(50,mean=-2,sd=1),rnorm(50,mean=2,sd=1))
par(mar=c(2,4,.5,.5))
hist(x,main='',xlab='')
n = length(x)
gamaic = matrix(nrow=n,ncol=2);
#theta = (mu1,mu2,sigma2.1,sigma2.2,p1)
theta.curr = c(1,1,10,1,.2)
theta.trace = matrix(nrow=1,ncol=length(theta.curr));
theta.trace[1,] = theta.curr
max_ite = 100;
O que fizemos:
- Computamos o tamanho da amostra,
- Criamos uma matriz para conter todos os valores de $\gamma_{ic} = P(Y_i = c|x_i)$.
- Computamos valores iniciais para o vetor de parâmetros $(\mu_1,\mu_2,\sigma^2_1,\sigma^2_2,p_1)$.
- Criamos uma matriz para armazenar a trajetória dos parâmetros.
- Colocamos os valores iniciais na primeira linha da matriz de trajetórias.
- Definimos o número máximo de iterações do algoritmo como 100.
Nosso algoritmo é iterativo, logo deve estar dentro de uma estrutura de repetição:
for(j in 1:max_ite){
(...)
}
Vamos começar a programar nosso método de estimação. Primeiro vamos armazenar os valores de theta.curr em variáveis com os nomes dos parâmetros para o código ficar mais legível (mas não é necessário).
mu1 = theta.curr[1];
mu2 = theta.curr[2];
sigma2.1 = theta.curr[3];
sigma2.2 = theta.curr[4];
p1 = theta.curr[5];
p2 = 1-p1;
Agora, vamos implementar o passo E. Nesse passo, temos que calcular os valores de $\gamma_{ic}$ dados os valores de $x$ e $\theta^{(j-1)}$.
#PASSO E - CALCULO DE GAMMAIC
for(i in 1:n){
d1 = p1*dnorm(x[i],mean=mu1,sd=sqrt(sigma2.1));
d2 = p2*dnorm(x[i],mean=mu2,sd=sqrt(sigma2.2));
gamaic[i,1] = d1/(d1+d2);
gamaic[i,2] = d2/(d1+d2);
}
Agora, dados os valores de $\gamma_{ic}$, vamos maximizar a função $Q$ e estimar os parâmetros (passo M), calculando as estimativas dos parâmetros.
#PASSO M - ESTIMACAO
sg1 = sum(gamaic[,1]);
sg2 = sum(gamaic[,2]);
mu1.hat = sum(gamaic[,1]*x)/sg1;
mu2.hat = sum(gamaic[,2]*x)/sg2;
delta = (x-mu1.hat)^2;
sigma2.1.hat = sum(gamaic[,1]*delta)/sg1;
delta = (x-mu2.hat)^2;
sigma2.2.hat = sum(gamaic[,2]*delta)/sg2;
p1.hat = sg1/(sg1+sg2);
Com as novas estimativas dos parâmetros, precisamos atualizar o vetor theta.curr e armazenar essas estimativas na matriz de trajetórias.
theta.curr = c(mu1.hat,mu2.hat,sigma2.1.hat,sigma2.2.hat,p1.hat)
theta.trace = rbind(theta.trace,theta.curr);
O algoritmo está pronto. Mas temos que executar todo o bloco de código cada vez que quisermos realizar a estimação. Dessa forma, o código abaixo apresenta o algoritmo dentro de uma função.
GMM_EM = function(x,theta_ini,max_ite){
n = length(x)
gamaic = matrix(nrow=n,ncol=2);
theta.curr = theta_ini;
theta.trace = matrix(nrow=1,ncol=length(theta.curr));
theta.trace[1,] = theta.curr
for(j in 1:max_ite){
mu1 = theta.curr[1];
mu2 = theta.curr[2];
sigma2.1 = theta.curr[3];
sigma2.2 = theta.curr[4];
p1 = theta.curr[5];
p2 = 1-p1;
#PASSO E - CALCULO DE GAMMAIC
for(i in 1:n){
d1 = p1*dnorm(x[i],mean=mu1,sd=sqrt(sigma2.1));
d2 = p2*dnorm(x[i],mean=mu2,sd=sqrt(sigma2.2));
gamaic[i,1] = d1/(d1+d2);
gamaic[i,2] = d2/(d1+d2);
}
#PASSO M - ESTIMACAO
sg1 = sum(gamaic[,1]);
sg2 = sum(gamaic[,2]);
mu1.hat = sum(gamaic[,1]*x)/sg1;
mu2.hat = sum(gamaic[,2]*x)/sg2;
delta = (x-mu1.hat)^2;
sigma2.1.hat = sum(gamaic[,1]*delta)/sg1;
delta = (x-mu2.hat)^2;
sigma2.2.hat = sum(gamaic[,2]*delta)/sg2;
p1.hat = sg1/(sg1+sg2);
theta.curr = c(mu1.hat,mu2.hat,sigma2.1.hat,sigma2.2.hat,p1.hat)
theta.trace = rbind(theta.trace,theta.curr);
}
return(theta.trace);
}
Vamos testar nosso código com os dados que geramos inicialmente.
set.seed(7)
x = c(rnorm(50,mean=-2,sd=1),rnorm(50,mean=2,sd=1))
theta_ini = c(1,1,10,1,.2);
max_ite = 100;
res = GMM_EM(x,theta_ini,max_ite);
>head(res);
[,1] [,2] [,3] [,4] [,5]
1.000000 1.000000 10.000000 1.000000 0.2000000
theta.curr -1.139293 1.070248 4.817979 2.227314 0.4216040
theta.curr -1.194942 1.261527 3.356893 2.722751 0.4570911
theta.curr -1.313885 1.418818 2.698869 2.608024 0.4684452
theta.curr -1.477107 1.584125 2.149425 2.196990 0.4721721
theta.curr -1.655599 1.751241 1.567991 1.660437 0.4733255
Agora vamos plotar os gráficos com as trajetórias das estimativas dos parâmetros.
par(mfrow=c(2,3))
plot(res[,1],type='l',main=expression(mu[1])); abline(h=-2,col='blue');
plot(res[,2],type='l',main=expression(mu[2])); abline(h=2,col='blue');
plot(res[,3],type='l',main=expression(sigma[1]^2 )); abline(h=1,col='blue');
plot(res[,4],type='l',main=expression(sigma[2]^2 )); abline(h=1,col='blue');
plot(res[,5],type='l',main=expression(p[1])); abline(h=0.5,col='blue');
>res[(max_ite+1),] #mu1 mu2 sigma2.1 sigma2.2 p1
[1] -1.9704849 1.8669399 0.6421497 1.0473874 0.4503654
>1 - res[(max_ite+1),5] #p2
0.5496346
Agora temos um algoritmo para estimar os parâmetros do modelo de misturas de duas distribuições Gaussianas! Incrível, não acha?
Porém, há um problema que precisamos discutir. O que aconteceria se assumíssemos valores iniciais iguais para os parâmetros das duas distribuições? Vamos testar.
theta_ini = c(1,1,1,1,.2);
max_ite = 100;
res = GMM_EM(x,theta_ini,max_ite);
> head(res);
[,1] [,2] [,3] [,4] [,5]
1.0000000 1.0000000 1.000000 1.000000 0.2
theta.curr 0.1386966 0.1386966 4.510061 4.510061 0.2
theta.curr 0.1386966 0.1386966 4.510061 4.510061 0.2
theta.curr 0.1386966 0.1386966 4.510061 4.510061 0.2
theta.curr 0.1386966 0.1386966 4.510061 4.510061 0.2
theta.curr 0.1386966 0.1386966 4.510061 4.510061 0.2
Observe que a estimativa de $p_1$ não se altera e, após a primeira iteração, as estimativas dos demais parâmetros permanecem constantes. Por que isso ocorre? Observe que assumimos $\mu_1^{(0)} = \mu_2^{(0)}$ e $(\sigma^2_1)^{(0)} = (\sigma^2_2)^{(0)}$ e ,como consequência, $f(x_i|\theta_1) = f(x_i|\theta_2)$. Dessa forma, $\gamma_{i1} = p_1^{(0)}$ e $\gamma_{i2} = p_2^{(0)}$ e, como consequência,
\begin{eqnarray} \hat{\mu}_1 = \hat{\mu}_2 &=& \dfrac{\sum_{i-1}^{n} x_i}{n}, \nonumber\\ \hat{\sigma}^2_1 = \hat{\sigma}^2_2 &=& \dfrac{\sum_{i-1}^{n}(x_i-\hat{\mu}_1)^2}{n}, \nonumber\\ \hat{p}_1 &=& p_1^{(0)}. \nonumber \end{eqnarray}
Logo, para que o algoritmo funcione, precisamos assumir valores iniciais distintos para as distribuições. Se consideramos mesma média, devemos assumir variâncias diferentes e vice-versa.
Na próxima aula vamos expandir a ideia do modelo para uma mistura de $n$ distribuições Gaussianas.
Espero que tenha gostado da aula.
Até a próxima aula!
Nenhum comentário:
Postar um comentário