Algoritmo EM implementado manualmente

20

Eu quero implementar o algoritmo EM manualmente e depois compará-lo com os resultados normalmixEMdo mixtoolspacote. Claro, eu ficaria feliz se os dois tivessem os mesmos resultados. A referência principal é Geoffrey McLachlan (2000), Modelos de Mistura Finita .

Eu tenho uma densidade de mistura de dois gaussianos, de forma geral, a probabilidade de log é dada por (McLachlan página 48):

logLc(Ψ)=i=1gj=1nzij{logπi+logfi(yi;θi)}.
Ozij são1 , se a observação foi a partir doEu thdensidade componente, caso contrário0 0 . AfEu é a densidade da distribuição normal. oπ é a proporção da mistura, entãoπ1 é a probabilidade de que uma observação seja da primeira distribuição gaussiana eπ2 é a probabilidade de que uma observação seja da segunda distribuição gaussiana.

A etapa E é agora, cálculo da expectativa condicional:

Q(Ψ;Ψ(0 0))=EΨ(0 0){registroeuc(|Ψ)|y}.
o que leva, após algumas derivações ao resultado (página 49):

τEu(yj;Ψ(k))=πEu(k)fEu(yj;θEu(k)f(yj;Ψ(k)=πEu(k)fEu(yj;θEu(k)h=1gπh(k)fh(yj;θh(k))
no caso de dois gaussianos (página 82):

τEu(yj;Ψ)=πEuϕ(yj;μEu,ΣEu)h=1gπhϕ(yj;μh,Σh)
AetapaMagora é a maximização de Q (página 49):

Q(Ψ;Ψ(k))=i=1gj=1nτi(yj;Ψ(k)){logπi+logfi(yj;θi)}.
This leads to (in the case of two Gaussians) (page 82):

μi(k+1)=j=1nτij(k)yjj=1nτij(k)Σi(k+1)=j=1nτij(k)(yjμi(k+1))(yjμi(k+1))Tj=1nτij(k)
and we know that (p. 50)

πi(k+1)=j=1nτi(yj;Ψ(k))n(i=1,,g).
We repeat the E, M steps until L(Ψ(k+1))L(Ψ(k)) is small.

I tried to write a R code (data can be found here).

# EM algorithm manually
# dat is the data

# initial values
pi1       <-  0.5
pi2       <-  0.5
mu1       <- -0.01
mu2       <-  0.01
sigma1    <-  0.01
sigma2    <-  0.02
loglik[1] <-  0
loglik[2] <- sum(pi1*(log(pi1) + log(dnorm(dat,mu1,sigma1)))) + 
             sum(pi2*(log(pi2) + log(dnorm(dat,mu2,sigma2))))

tau1 <- 0
tau2 <- 0
k    <- 1

# loop
while(abs(loglik[k+1]-loglik[k]) >= 0.00001) {

  # E step
  tau1 <- pi1*dnorm(dat,mean=mu1,sd=sigma1)/(pi1*dnorm(x,mean=mu1,sd=sigma1) + 
          pi2*dnorm(dat,mean=mu2,sd=sigma2))
  tau2 <- pi2*dnorm(dat,mean=mu2,sd=sigma2)/(pi1*dnorm(x,mean=mu1,sd=sigma1) + 
          pi2*dnorm(dat,mean=mu2,sd=sigma2))

  # M step
  pi1 <- sum(tau1)/length(dat)
  pi2 <- sum(tau2)/length(dat)

  mu1 <- sum(tau1*x)/sum(tau1)
  mu2 <- sum(tau2*x)/sum(tau2)

  sigma1 <- sum(tau1*(x-mu1)^2)/sum(tau1)
  sigma2 <- sum(tau2*(x-mu2)^2)/sum(tau2)

  loglik[k] <- sum(tau1*(log(pi1) + log(dnorm(x,mu1,sigma1)))) + 
               sum(tau2*(log(pi2) + log(dnorm(x,mu2,sigma2))))
  k         <- k+1
}


# compare
library(mixtools)
gm <- normalmixEM(x, k=2, lambda=c(0.5,0.5), mu=c(-0.01,0.01), sigma=c(0.01,0.02))
gm$lambda
gm$mu
gm$sigma

gm$loglik

The algorithm is not working, since some observations have the likelihood of zero and the log of this is -Inf. Where is my mistake?

Stat Tistician
fonte
The problem is not a statistical one, but rather a numerical one. You should add contingencies for likelihoods smaller than machine precision in your code.
JohnRos
why dont you try veryfying the mixtools function with a very simple example that can be verified by hand , say just five or ten values and two timeseries,first. then, if you find it works there, generalize your code and verify at each step.

Respostas:

17

You have several problems in the source code:

  1. As @Pat pointed out, you should not use log(dnorm()) as this value can easily go to infinity. You should use logmvdnorm

  2. When you use sum, be aware to remove infinite or missing values

  3. You looping variable k is wrong, you should update loglik[k+1] but you update loglik[k]

  4. The initial values for your method and mixtools are different. You are using Σ in your method, but using σ for mixtools(i.e. standard deviation, from mixtools manual).

  5. Your data do not look like a mixture of normal (check histogram I plotted at the end). And one component of the mixture has very small s.d., so I arbitrarily added a line to set τ1 and τ2 to be equal for some extreme samples. I add them just to make sure the code can work.

Eu também sugiro que você coloque códigos completos (por exemplo, como você inicializa o loglik []) no seu código-fonte e indente o código para facilitar a leitura.

Afinal, obrigado por apresentar o pacote mixtools e pretendo usá-los em minhas pesquisas futuras.

Eu também coloquei meu código de trabalho para sua referência:

# EM algorithm manually
# dat is the data
setwd("~/Downloads/")
load("datem.Rdata")
x <- dat

# initial values
pi1<-0.5
pi2<-0.5
mu1<--0.01
mu2<-0.01
sigma1<-sqrt(0.01)
sigma2<-sqrt(0.02)
loglik<- rep(NA, 1000)
loglik[1]<-0
loglik[2]<-mysum(pi1*(log(pi1)+log(dnorm(dat,mu1,sigma1))))+mysum(pi2*(log(pi2)+log(dnorm(dat,mu2,sigma2))))

mysum <- function(x) {
  sum(x[is.finite(x)])
}
logdnorm <- function(x, mu, sigma) {
  mysum(sapply(x, function(x) {logdmvnorm(x, mu, sigma)}))  
}
tau1<-0
tau2<-0
#k<-1
k<-2

# loop
while(abs(loglik[k]-loglik[k-1]) >= 0.00001) {
  # E step
  tau1<-pi1*dnorm(dat,mean=mu1,sd=sigma1)/(pi1*dnorm(x,mean=mu1,sd=sigma1)+pi2*dnorm(dat,mean=mu2,sd=sigma2))
  tau2<-pi2*dnorm(dat,mean=mu2,sd=sigma2)/(pi1*dnorm(x,mean=mu1,sd=sigma1)+pi2*dnorm(dat,mean=mu2,sd=sigma2))
  tau1[is.na(tau1)] <- 0.5
  tau2[is.na(tau2)] <- 0.5

  # M step
  pi1<-mysum(tau1)/length(dat)
  pi2<-mysum(tau2)/length(dat)

  mu1<-mysum(tau1*x)/mysum(tau1)
  mu2<-mysum(tau2*x)/mysum(tau2)

  sigma1<-mysum(tau1*(x-mu1)^2)/mysum(tau1)
  sigma2<-mysum(tau2*(x-mu2)^2)/mysum(tau2)

  #  loglik[k]<-sum(tau1*(log(pi1)+log(dnorm(x,mu1,sigma1))))+sum(tau2*(log(pi2)+log(dnorm(x,mu2,sigma2))))
  loglik[k+1]<-mysum(tau1*(log(pi1)+logdnorm(x,mu1,sigma1)))+mysum(tau2*(log(pi2)+logdnorm(x,mu2,sigma2)))
  k<-k+1
}

# compare
library(mixtools)
gm<-normalmixEM(x,k=2,lambda=c(0.5,0.5),mu=c(-0.01,0.01),sigma=c(0.01,0.02))
gm$lambda
	gm$mu
gm$sigma

gm$loglik

Historograma Histograma

zhanxw
fonte
@zahnxw obrigado pela sua resposta, isso significa que meu código está errado? Então a ideia básica não está funcionando?
Stat Tistician
"Eu também sugiro que você coloque códigos completos (por exemplo, como você inicializa o loglik []) no seu código-fonte e indente o código para facilitar a leitura." Bem, este é o meu código? o loglik [] é definido como eu o declarei no código que publiquei?
Stat Tistician
1
@StatTistician, a ideia está correta, mas a implementação tem falhas. Por exemplo, você não considerou o fluxo insuficiente. Além disso, o loop da variável k é confuso; você primeiro define loglik [1] e loglik [2]; depois de inserir o loop while, você define loglik [1] novamente. Esta não é a maneira natural de fazer. Minha sugestão sobre a inicialização do loglik [] significa código:, loklik <- rep(NA, 100)que pré-alocará o loglik [1], loglik [2] ... loglik [100]. Eu levanto essa pergunta porque, no seu código original, não encontrei o delcaration do loglik, talvez o código esteja truncado durante a colagem?
Zhanxw
Como eu postei abaixo: Obrigado por sua ajuda, mas estou descartando este tópico, pois ele é muito avançado para mim.
Stat Tistician
Existe agora uma maneira de determinar qual parte dos dados pertence a qual mistura?
Cardeal
2

Eu continuo recebendo um erro ao tentar abrir o arquivo .rar, mas isso pode ser apenas eu que estou fazendo algo bobo.

Não vejo erros óbvios no seu código. Um possível motivo para você obter zeros é devido à precisão do ponto flutuante. Lembre-se, quando você calculaf(y;θ), você está avaliando exp(-0,5(y-μ)2/σ2). Não é preciso uma grande diferença entreμ e ypara que isso seja arredondado para 0 quando você o faz em um computador. Isso é duplamente perceptível nos modelos de mistura, pois alguns de seus dados não serão "atribuídos" a cada componente da mistura e, portanto, podem ficar muito distantes dele. Em teoria, esses pontos também devem acabar com um baixo valor deτ quando você avalia a probabilidade do log, combatendo o problema - mas, graças ao erro de ponto flutuante, a quantidade já foi avaliada como -Inf nesse estágio, para que tudo quebre :).

Se esse for o problema, existem algumas soluções possíveis:

Um é mover o seu τdentro do logaritmo. Então, em vez de avaliar

τregistro(f(y|θ))

Avalie

registro(f(y|θ)τ).

Matematicamente o mesmo, mas pense no que acontece quando f(y|θ) e τ são 0 0. Atualmente você obtém:

  • 0 0registro(0 0)=0 0(-Eunf)=NumaN

mas com tau mudou você começa

  • registro(0 00 0)=registro(1)=0 0

assumindo que R avalia 0 00 0=1 (Não sei se funciona ou não, pois tenho tendência a usar o matlab)

Outra solução é expandir as coisas dentro do logaritmo. Supondo que você esteja usando logaritmos naturais:

τregistro(f(y|θ))

=τregistro(exp(-0,5(y-μ)2/σ2)/2πσ2)

=-0,5τregistro(2πσ2)-0,5τ(y-μ)2σ2.

Matematicamente o mesmo, mas deve ser mais resistente a erros de ponto flutuante, pois você evitou calcular uma grande potência negativa. Isso significa que você não pode mais usar a função de avaliação de norma incorporada, mas se isso não for um problema, essa provavelmente é a melhor resposta. Por exemplo, digamos que temos a situação em que

-0,5(y-μ)2σ2=-0,540.2=-800.

Avalie isso como sugeri, e você receberá -800. No entanto, no matlab, se expomos o take the log, obtemosregistro(exp(-800))=registro(0 0)=-Eunf.

Pat
fonte
Para ser sincero: não sou bom o suficiente para fazer essa coisa funcionar. O que me interessava é: Posso obter o mesmo resultado com meu algoritmo que a versão implementada do pacote mixtools. Mas, do meu ponto de vista, isso parece estar pedindo a lua. Mas acho que você se esforçou em sua resposta, então eu aceitarei! Obrigado!
Stat Tistician