EM алгоритм реализован вручную

Я хочу реализовать алгоритм EM вручную , а затем сравнить его с результатами normalmixEMиз mixtoolsпакета. Конечно, я был бы счастлив, если бы они оба привели к одинаковым результатам. Основное упоминание - Джеффри МакЛахлан (2000), Модели конечных смесей .

У меня плотность смеси двух гауссианов, в общем виде, логарифмическая вероятность определяется (McLachlan стр. 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}$ являются

1

$1$ , если наблюдение было от

i

$i$ ^{- го}компонента плотности,противном случае

0

$0$ .

f_{i}

$f_i$ плотность нормального распределения.

π

$\pi$ - пропорция смеси, поэтому

π_{1}

$\pi_1$ - вероятность того, что наблюдение происходит из первого гауссовского распределения, а

π_{2}

$\pi_2$ - вероятность того, что наблюдение происходит из второго гауссовского распределения.

Е шаг теперь, вычисление условного ожидания:

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

$Q(\Psi;\Psi^{(0)}) = E_{\Psi(0)}\{\log L_c(|\Psi)|y\}.$ что приводит после нескольких выводов к результату (стр. 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}$ в случае двух гауссовых (стр 82):

τ_{i} (y_{j}; Ψ) = \frac{π_{i} ϕ (y_{j}; μ_{i}, Σ_{i})}{\sum_{h = 1}^{g} π_{h} ϕ (y_{j}; μ_{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)}$ Мшаг теперь максимизация Q (страница 49):

Q (Ψ; Ψ^{(k)}) = \sum_{i = 1}^{g} \sum_{j = 1}^{n} τ_{i} (y_{j}; Ψ^{(k)}) {\log π_{i} + \log f_{i} (y_{j}; θ_{i})} .

$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)\}.$ Это приводит к (в случае двух гауссиан) (стр. 82):

\begin{aligned} μ_{i}^{(k + 1)} & = \frac{\sum_{j = 1}^{n} τ_{i j}^{(k)} y_{j}}{\sum_{j = 1}^{n} τ_{i j}^{(k)}} \\ Σ_{i}^{(k + 1)} & = \frac{\sum_{j = 1}^{n} τ_{i j}^{(k)} (y_{j} - μ_{i}^{(k + 1)}) (y_{j} - μ_{i}^{(k + 1)})^{T}}{\sum_{j = 1}^{n} τ_{i j}^{(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}$ , и мы знаемчто (стр. 50)

π_{i}^{(k + 1)} = \frac{\sum_{j = 1}^{n} τ_{i} (y_{j}; Ψ^{(k)})}{n} (i = 1, \dots, g) .

$\pi_i^{(k+1)} = \frac{\sum_{j=1}^n \tau_i(y_j;\Psi^{(k)})}{n}\qquad (i = 1, \ldots, g).$ Мы повторяем шаги E, M до тех пор, пока

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

$L(\Psi^{(k+1)})-L(\Psi^{(k)})$ станет маленьким.

Я пытался написать код R (данные можно найти здесь ).

# 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

Алгоритм не работает, так как некоторые наблюдения имеют вероятность нуля, и лог этого -Inf. Где моя ошибка?

r expectation-maximization gaussian-mixture Стат Тистициан
источник

Проблема не статистическая, а количественная. Вы должны добавить непредвиденные обстоятельства для вероятностей, меньших, чем точность машины, в ваш код.

JohnRos

почему бы вам не попробовать выполнить функцию mixtools на очень простом примере, который можно проверить вручную, скажем, сначала пять или десять значений и два временных ряда, во-первых. затем, если вы обнаружите, что он работает там, обобщите ваш код и проверяйте его на каждом этапе.

Ответы:

У вас есть несколько проблем в исходном коде:

Как указал @Pat, вы не должны использовать log (dnorm ()), так как это значение может легко переходить в бесконечность. Вы должны использовать logmvdnorm
Когда вы используете сумму , будьте внимательны, чтобы удалить бесконечные или отсутствующие значения
Вы зацикливаете переменную k неправильно, вы должны обновить loglik [k + 1], но вы обновляете loglik [k]
$\Sigma$ $\sigma$
$\tau_1$ $\tau_2$

Я также предлагаю вам вставить полные коды (например, как вы инициализируете loglik []) в свой исходный код и сделать отступ для кода, чтобы его было легко читать.

В конце концов, спасибо за введение пакета mixtools , и я планирую использовать их в своих будущих исследованиях.

Я также поставил свой рабочий код для вашей справки:

# 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 Гистограмма

zhanxw
источник

@zahnxw спасибо за ваш ответ, значит ли это, что мой код неверен? Так что основная идея не работает?

Стат Тистициан

«Я также предлагаю вам вставить полные коды (например, как вы инициализируете loglik []) в свой исходный код и сделать отступ для кода, чтобы его было легко читать». Ну это мой код? loglik [] определяется так, как я объявил в коде, который я разместил?

Стат Тистициан

@StatTistician идея верна, но реализация имеет недостатки. Например, вы не учли недостаток. Кроме того, зацикливание переменной k сбивает с толку, вы сначала устанавливаете loglik [1] и loglik [2], после входа в цикл while вы снова устанавливаете loglik [1]. Это не естественный способ сделать. Мое предложение об инициализации loglik [] означает код:, loklik <- rep(NA, 100)который будет предварительно выделять loglik [1], loglik [2] ... loglik [100]. Я поднимаю этот вопрос, потому что в вашем исходном коде я не нашел делкарации loglik, может быть, код обрезается при вставке?

zhanxw

Как я писал ниже: Спасибо за вашу помощь, но я опускаю эту тему, так как она слишком сложна для меня.

Стат Тистициан

Есть ли теперь способ определить, какая часть данных принадлежит какой смеси?

кардинал

Я продолжаю получать сообщение об ошибке при попытке открыть ваш .rar файл, но это может быть просто из-за того, что я делаю что-то глупое.

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

Если это проблема, есть несколько возможных решений:

$\tau$

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

оценивать

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

$f(y|\theta)$ $\tau$ $\approx 0$

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

но с тау переехал ты получаешь

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

$0^0 = 1$

Другое решение состоит в том, чтобы расширить содержимое внутри логарифма. Предполагая, что вы используете натуральные логарифмы:

$\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}$

Математически то же самое, но должно быть более устойчивым к ошибкам с плавающей запятой, поскольку вы избегаете вычисления большой отрицательной мощности. Это означает, что вы больше не можете использовать встроенную функцию оценки норм, но если это не проблема, возможно, это лучший ответ. Например, скажем, у нас есть ситуация, когда

$-0.5\frac{(y-\mu)^2}{\sigma^2} = -0.5*40^2 = -800$

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

похлопывание
источник

хм, если честно: я недостаточно хорош, чтобы заставить эту вещь работать. Меня заинтересовало следующее: могу ли я получить тот же результат с моим алгоритмом, что и реализованная версия пакета mixtools. Но с моей точки зрения это похоже на луну. Но я думаю, что вы приложили усилия к своему ответу, поэтому я приму его! Благодарность!

Стат Тициан