Algorytm EM ręcznie wdrożony

20

Chcę, aby zaimplementować algorytm EM ręcznie, a następnie porównać je do wyników działań normalmixEMz mixtoolsopakowania. Oczywiście byłbym szczęśliwy, gdyby oba doprowadziły do ​​tych samych rezultatów. Głównym odniesieniem jest Geoffrey McLachlan (2000), Finite Mixture Models .

Mam gęstość mieszanki dwóch Gaussów, w ogólnej formie, logarytmiczne prawdopodobieństwo podaje (McLachlan strona 48):

logLc(Ψ)=i=1gj=1nzij{logπi+logfi(yi;θi)}.
zij1 , jeśli obserwacja była zi XXgęstości składnika, inaczej0 . fi gęstość rozkładu normalnego. Theπ jest proporcją mieszaniny, więcπ1 to prawdopodobieństwo, że obserwacja pochodzi z pierwszego rozkładu Gaussa, aπ2 to prawdopodobieństwo, że obserwacja pochodzi z drugiego rozkładu Gaussa.

E krokiem jest teraz obliczenie warunkowa wartość oczekiwana:

Q(Ψ;Ψ(0))=EΨ(0){logLc(|Ψ)|y}.
co po kilku pochodnych prowadzi do wyniku (strona 49):

τi(yj;Ψ(k))=πi(k)fi(yj;θi(k)f(yj;Ψ(k)=πi(k)fi(yj;θi(k)h=1gπh(k)fh(yj;θh(k))
w przypadku dwóch Gaussów (strona 82):

τja(yjot;Ψ)=πjaϕ(yjot;μja,Σja)h=1solπhϕ(yjot;μh,Σh)
Mkrokiem jest maksymalizacja Q (strona 49):

Q(Ψ;Ψ(k))=ja=1soljot=1nτja(yjot;Ψ(k)){logπja+logfaja(yjot;θja)}.
Prowadzi to do (w przypadku dwóch Gaussów) (strona 82):

μja(k+1)=jot=1nτjajot(k)yjotjot=1nτjajot(k)Σja(k+1)=jot=1nτjajot(k)(yjot-μja(k+1))(yjot-μja(k+1))T.jot=1nτjajot(k)
i wiemy, że (s. 50)

πja(k+1)=jot=1nτja(yjot;Ψ(k))n(ja=1,,sol).
Powtarzamy kroki E, M, ażL.(Ψ(k+1))-L.(Ψ(k)) będzie małe.

Próbowałem napisać kod R (dane można znaleźć tutaj ).

# 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

Algorytm nie działa, ponieważ niektóre obserwacje mają prawdopodobieństwo zerowe, a logarytm z tego jest -Inf. Gdzie jest mój błąd?

Stat Tistician
źródło
Problem nie jest statystyczny, ale liczbowy. W kodzie należy dodać ewentualności mniejsze niż precyzja maszyny.
JohnRos
dlaczego nie spróbujesz zweryfikować funkcji mixtools za pomocą bardzo prostego przykładu, który można zweryfikować ręcznie, powiedz najpierw pięć lub dziesięć wartości i dwa szeregi czasowe, najpierw. następnie, jeśli okaże się, że tam działa, uogólnij swój kod i weryfikuj na każdym kroku.

Odpowiedzi:

17

Masz kilka problemów z kodem źródłowym:

  1. Jak wskazał @Pat, nie powinieneś używać log (dnorm ()), ponieważ ta wartość może łatwo przejść do nieskończoności. Powinieneś użyć logmvdnorm

  2. Kiedy używasz sumy , pamiętaj, aby usunąć nieskończone lub brakujące wartości

  3. Pętla zmiennej k jest niepoprawna, powinieneś zaktualizować loglik [k + 1], ale zaktualizujesz loglik [k]

  4. Σσ

  5. τ1τ2)

Sugeruję również, abyś umieścił pełne kody (np. Jak zainicjalizujesz loglik []) w kodzie źródłowym i wciąć kod, aby ułatwić czytanie.

W końcu dziękuję za wprowadzenie pakietu mixtools i planuję wykorzystać je w moich przyszłych badaniach.

Podaję również mój kod roboczy w celach informacyjnych:

# 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

Historgram Histogram

zhanxw
źródło
@zahnxw dzięki za odpowiedź, czy to oznacza, że ​​mój kod jest nieprawidłowy? Więc pomysł basi nie działa?
Stat Tistician
„Sugeruję również, abyś umieścił w kodzie źródłowym kompletne kody (np. Jak zainicjalizujesz loglik []) i wciąć kod, aby ułatwić czytanie.” To jest mój kod? loglik [] jest zdefiniowany tak, jak zadeklarowałem go w opublikowanym przeze mnie kodzie?
Stat Tistician
1
@StatTistician pomysł jest poprawny, ale implementacja ma wady. Na przykład nie uwzględniono niedomiaru. Również zapętlenie zmiennej k jest mylące, najpierw ustawiasz loglik [1] i loglik [2], po przejściu do pętli while ponownie ustawiasz loglik [1]. To nie jest naturalny sposób. Moja sugestia dotycząca inicjowania loglik [] oznacza kod:, loklik <- rep(NA, 100)który wstępnie przydzieli loglik [1], loglik [2] ... loglik [100]. Podnoszę to pytanie, ponieważ w twoim oryginalnym kodzie nie znalazłem delcaration loglik, może kod jest obcinany podczas wklejania?
zhanxw
Jak napisałem poniżej: Dziękuję za pomoc, ale upuszczam ten temat, ponieważ jest dla mnie zbyt zaawansowany.
Stat Tistician
Czy istnieje sposób na określenie, która część danych należy do której mieszaniny?
Kardynał
2

Ciągle pojawia się błąd podczas próby otwarcia pliku .rar, ale to może być po prostu to, że robię coś głupiego.

fa(y;θ)exp(-0,5(y-μ)2)/σ2))μyτ

Jeśli to jest problem, istnieje kilka możliwych rozwiązań:

Jednym z nich jest przeniesienie twojego τ

τlog(fa(y|θ))

oceniać

log(fa(y|θ)τ)

fa(y|θ)τ0 . Obecnie otrzymujesz:

  • 0log(0)=0(-janfa)=N.zaN.

ale z tau poruszasz się

  • log(00)=log(1)=0

00=1

Innym rozwiązaniem jest rozszerzenie elementów wewnątrz logarytmu. Zakładając, że używasz logarytmów naturalnych:

τlog(fa(y|θ))

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

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

Matematycznie to samo, ale powinien być bardziej odporny na błędy zmiennoprzecinkowe, ponieważ uniknąłeś obliczenia dużej mocy ujemnej. Oznacza to, że nie możesz już korzystać z wbudowanej funkcji oceny norm, ale jeśli nie jest to problem, prawdopodobnie jest to lepsza odpowiedź. Załóżmy na przykład, że mamy sytuację, w której

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

log(exp(-800))=log(0)=-janfa

Poklepać
źródło
mh, szczerze mówiąc: nie jestem wystarczająco dobry, aby uruchomić tę rzecz. Interesowało mnie to: czy mogę uzyskać taki sam wynik dzięki mojemu algorytmowi, jak zaimplementowana wersja pakietu mixtools. Ale z mojego punktu widzenia wydaje się, że to prosi o księżyc. Ale myślę, że wkładasz wysiłek w swoją odpowiedź, więc zaakceptuję to! Dzięki!
Stat Tistician