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