Algorytm EM ręcznie wdrożony

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):

\log L_{c} (Ψ) = \sum_{i = 1}^{g} \sum_{j = 1}^{n} z_{i j} {\log π_{i} + \log f_{i} (y_{i}; θ_{i})} .

$\log L_c(\Psi) = \sum_{i=1}^g \sum_{j=1}^n z_{ij}\{\log \pi_i + \log f_i(y_i;\theta_i)\}.$

z_{i j}

$z_{ij}$ są

1

$1$ , jeśli obserwacja była z

i

$i$ ^XXgęstości składnika, inaczej

0

$0$ .

f_{i}

$f_i$ gęstość rozkładu normalnego. The

π

$\pi$ jest proporcją mieszaniny, więc

π_{1}

$\pi_1$ to prawdopodobieństwo, że obserwacja pochodzi z pierwszego rozkładu Gaussa, a

π_{2}

$\pi_2$ to prawdopodobieństwo, że obserwacja pochodzi z drugiego rozkładu Gaussa.

E krokiem jest teraz obliczenie warunkowa wartość oczekiwana:

Q (Ψ; Ψ^{(0)}) = E_{Ψ (0)} {\log L_{c} (| Ψ) | y} .

$Q(\Psi;\Psi^{(0)}) = E_{\Psi(0)}\{\log L_c(|\Psi)|y\}.$ co po kilku pochodnych prowadzi do wyniku (strona 49):

\begin{aligned} τ_{i} (y_{j}; Ψ^{(k)}) & = \frac{π_{i}^{(k)} f_{i} (y_{j}; θ_{i}^{(k)}}{f (y_{j}; Ψ^{(k)}} \\ = \frac{π_{i}^{(k)} f_{i} (y_{j}; θ_{i}^{(k)}}{\sum_{h = 1}^{g} π_{h}^{(k)} f_{h} (y_{j}; θ_{h}^{(k)})} \end{aligned}

$\begin{align} \tau_i(y_j;\Psi^{(k)}) &= \frac{\pi_i^{(k)}f_i(y_j;\theta_i^{(k)}}{f(y_j;\Psi^{(k)}} \\[8pt] &= \frac{\pi_i^{(k)}f_i(y_j;\theta_i^{(k)}}{\sum_{h=1}^g \pi_h^{(k)}f_h(y_j;\theta_h^{(k)})} \end{align}$ w przypadku dwóch Gaussów (strona 82):

τ_{ja} (y_{jot}; Ψ) = \frac{π_{ja} ϕ (y_{jot}; μ_{ja}, Σ_{ja})}{\sum_{h = 1}^{sol} π_{h} ϕ (y_{jot}; μ_{h}, Σ_{h})}

$\tau_i(y_j;\Psi) = \frac{\pi_i \phi(y_j;\mu_i,\Sigma_i)}{\sum_{h=1}^g \pi_h\phi(y_j; \mu_h,\Sigma_h)}$ Mkrokiem jest maksymalizacja Q (strona 49):

Q (Ψ; Ψ^{(k)}) = \sum_{ja = 1}^{sol} \sum_{jot = 1}^{n} τ_{ja} (y_{jot}; Ψ^{(k)}) {\log π_{ja} + \log {fa}_{ja} (y_{jot}; θ_{ja})} .

$Q(\Psi;\Psi^{(k)}) = \sum_{i=1}^g\sum_{j=1}^n\tau_i(y_j;\Psi^{(k)})\{\log \pi_i + \log f_i(y_j;\theta_i)\}.$ Prowadzi to do (w przypadku dwóch Gaussów) (strona 82):

\begin{aligned} μ_{ja}^{(k + 1)} & = \frac{\sum_{jot = 1}^{n} τ_{ja jot}^{(k)} y_{jot}}{\sum_{jot = 1}^{n} τ_{ja jot}^{(k)}} \\ Σ_{ja}^{(k + 1)} & = \frac{\sum_{jot = 1}^{n} τ_{ja jot}^{(k)} (y_{jot} - μ_{ja}^{(k + 1)}) (y_{jot} - μ_{ja}^{(k + 1)})^{T.}}{\sum_{jot = 1}^{n} τ_{ja jot}^{(k)}} \end{aligned}

$\begin{align} \mu_i^{(k+1)} &= \frac{\sum_{j=1}^n \tau_{ij}^{(k)}y_j}{\sum_{j=1}^n \tau_{ij}^{(k)}} \\[8pt] \Sigma_i^{(k+1)} &= \frac{\sum_{j=1}^n \tau_{ij}^{(k)}(y_j - \mu_i^{(k+1)})(y_j - \mu_i^{(k+1)})^T}{\sum_{j=1}^n \tau_{ij}^{(k)}} \end{align}$ i wiemy, że (s. 50)

π_{ja}^{(k + 1)} = \frac{\sum_{jot = 1}^{n} τ_{ja} (y_{jot}; Ψ^{(k)})}{n} (ja = 1, \dots, sol) .

$\pi_i^{(k+1)} = \frac{\sum_{j=1}^n \tau_i(y_j;\Psi^{(k)})}{n}\qquad (i = 1, \ldots, g).$ Powtarzamy kroki E, M, aż

L (Ψ^{(k + 1)}) - L (Ψ^{(k)})

$L(\Psi^{(k+1)})-L(\Psi^{(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?

r expectation-maximization gaussian-mixture 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:

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

Jak wskazał @Pat, nie powinieneś używać log (dnorm ()), ponieważ ta wartość może łatwo przejść do nieskończoności. Powinieneś użyć logmvdnorm
Kiedy używasz sumy , pamiętaj, aby usunąć nieskończone lub brakujące wartości
Pętla zmiennej k jest niepoprawna, powinieneś zaktualizować loglik [k + 1], ale zaktualizujesz loglik [k]
$\Sigma$ $\sigma$
$\tau_1$ $\tau_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

@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ł

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.

$f(y;\theta)$ $\exp(-0.5(y-\mu)^2/\sigma^2)$ $\mu$ $y$ $\tau$

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

Jednym z nich jest przeniesienie twojego $\tau$

$\tau \log(f(y|\theta))$

oceniać

$\log \left( f(y|\theta)^\tau \right)$

$f(y|\theta)$ $\tau$ $\approx 0$ . Obecnie otrzymujesz:

$0 \log (0) = 0 (-Inf) = NaN$

ale z tau poruszasz się

$\log \left( 0^0\right) = \log(1) = 0$

$0^0 = 1$

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

$\tau \log(f(y|\theta))$

$= \tau \log(\exp(-0.5(y-\mu)^2/\sigma^2)/\sqrt{2\pi\sigma^2})$

$= -0.5\tau \log(2 \pi\sigma^2) - 0.5 \tau \frac{(y-\mu)^2}{\sigma^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\frac{(y-\mu)^2}{\sigma^2} = -0.5*40^2 = -800$

$\log(\exp(-800)) = \log(0) = -Inf$

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