\(\color{lightseagreen}{\text{EM алгоритам}}\)
EM (Expectation-Maximization) алгоритам је још један алгоритам за одређивање оцена методом максималне веродостојности, који се често користи у Бајесовској статистици и машинском учењу. Примењује се у случајевима са недостајућим подацима.
\(\color{lightseagreen}{\text{Пример 1}}\)
Нека су \(X_1, X_2, \ldots, X_n\) независне случајне величине из eкспоненцијалне расподеле са непознатим параметром \(\lambda\). Тада је оцена параметра \(\lambda\) добијена методом максималне веродостојности, за реализован узорак \(x_1, x_2, \ldots, x_n\), једнака \(\widehat \lambda=\frac{1}{\overline x_n}\).
Међутим, претпоставимо да уместо датих конкретних постигнутих вредности \(x_1, x_2, \ldots, x_n\), располажемо само са вредностима доњих и горњих граница интервала у којима се те вредности налазе. На пример, претпоставимо да је дат хистограм из експоненцијалне расподеле са непознатим параметром:
set.seed(2)
h <- hist(rexp(1000, rate = rexp(1, 0.1)), freq=FALSE)
h## $breaks
## [1] 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
##
## $counts
## [1] 592 243 91 45 20 5 4
##
## $density
## [1] 11.84 4.86 1.82 0.90 0.40 0.10 0.08
##
## $mids
## [1] 0.025 0.075 0.125 0.175 0.225 0.275 0.325
##
## $xname
## [1] "rexp(1000, rate = rexp(1, 0.1))"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"Сада, једноставним рачуном, користећи границе интервала којем припада \(x_i\), можемо наћи оцену вредности \(x_i\) као условно очекивање: \[\widehat x_i=E(X_i|X_i \in (u_i, v_i))=\frac{1}{\lambda}+\frac{u_ie^{-\lambda u_i}-v_ie^{-\lambda v_i}}{e^{-\lambda u_i}-e^{-\lambda v_i}}.\]
Примећујемо да се у оцени вредности \(x_i\) појављује параметар \(\lambda\), који је, такође, непознат, па ћемо на почетку узети неку произвољну оцену \(\widehat\lambda\) и итеративно понављати следећа два корака до конвергенције некој оцени непознатог параметра \(\lambda\).
- (Е корак) Оцену \(\widehat\lambda\) користимо да одредимо очекивање
\[\widehat x_i=E(X_i|X_i \in (u_i, v_i))=\frac{1}{\widehat\lambda}+\frac{u_ie^{-\widehat\lambda u_i}-v_ie^{-\widehat\lambda v_i}}{e^{-\widehat\lambda u_i}-e^{-\widehat\lambda v_i}}, \ i= 1, \ldots, n;\]
- (М корак) На основу добијених оцена вредности \(x_i, \ i= 1, \ldots, n\), одређујемо нову оцену параметра \(\lambda\)
\[\widehat\lambda=\arg\max_{\lambda}\log L(\lambda;\ {\bf \widehat x})=\frac{n}{\sum_{i=1}^n\widehat x_i}.\]
Имплементираћемо наведени алгоритам, а потом га применити на податке које смо добили на основу хистограма.
Е корак:
expected_x <- function(lambda, u, v) {
1/lambda +
(u*exp(-lambda*u) - v*exp(-lambda*v))/
(exp(-lambda*u) - exp(-lambda*v))
}М корак:
max_logL <- function(x) {
1/mean(x)
}Функција која итеративно спроводи кораке Е и М и тражи оцену ЕМ алгоритмом, на основу прослеђене почетне оцене параметра \(\lambda\) и вектора доњих (\(u\)) и горњих (\(v\)) граница интервала:
EM_estimate <- function(lambda_0, u, v, tol = 1e-8, maxiter = 1000) {
x <- expected_x(lambda_0, u, v)
lambda <- max_logL(x)
print(lambda)
iter <- 1
while((abs(lambda - lambda_0) > tol) &&
iter < maxiter) {
iter <- iter + 1
lambda_0 <- lambda
x <- expected_x(lambda_0, u, v)
lambda <- max_logL(x)
print(lambda)
}
lambda
}Одредићемо векторе граница интервала у нашем примеру:
u <- rep(head(h$breaks, -1), h$counts)
v <- rep(tail(h$breaks, -1), h$counts)Сада је оцена параметра \(\lambda\) добијена ЕМ алгоритмом једнака:
est=EM_estimate(1, u, v)## [1] 16.88001
## [1] 17.86531
## [1] 17.92865
## [1] 17.93273
## [1] 17.93299
## [1] 17.93301
## [1] 17.93301
## [1] 17.93301
## [1] 17.93301est## [1] 17.93301Проверићемо да ли се добијена оцена уклапа са подацима које имамо додвањем густине експоненцијалне расподеле са оцењеним параметром на хистограм.
set.seed(2)
h <- hist(rexp(1000, rate = rexp(1, 0.1)), freq=FALSE)
curve(dexp(x, rate = est), add=TRUE, col="red")
Видимо да крива добро прати хистограм, па можемо закључити да оцена приближно одговара правој вредности.
\(\color{lightseagreen}{\text{Пример 2}}\)
Експеримент се састоји од бацања неисправног новчића и уколико падне глава узима се узорак из нормалне \(\mathscr{N}(\mu_1, 1)\) расподеле, а уколико падне писмо (\(p=0.7\)) узима се узорак из нормалне \(\mathscr{N}(\mu_2, 1)\) расподеле.
Коришћењем ЕМ алгоритма треба одредити вредности параметара \(\mu_1\) и \(\mu_2\).
Посматрамо узорак из мешавине две нормалне расподеле са густином \[f(x)=(1-p)f_{\mu_1}(x)+pf_{\mu_2}.\]
Најпре генеришимо такав узорак, за вредности параметра \(\mu_1=-1\) и \(\mu_2=3\) (то су стварне вредности, које у наставку, на основу узорка, желимо да оценимо).
rmix <- function(n) {
u <- runif(n)
(u < 0.3) * rnorm(n, -1) +
(u >= 0.3) * rnorm(n, 3)
}
uzorak <- rmix(1000)Скривена променљива која се овде јавља је индикаторска \(W \sim Ber(p)\), која се односи на исход бацања новчића. Вредност тог индикатора можемо оценити са \[\widehat w_i=E(W_i|X_i=x_i).\] Даље се добија \[\widehat w_i=E(W_i=1|X_i=x_i)=\frac{pf_{\mu_2}(x_i)}{(1-p)f_{\mu_1}(x_i)+pf_{\mu_2}(x_i)}.\]
Овако добијена оцена представља Е корак алгоритма и може се једноставно имплементирати.
expectation <- function(mi1, mi2, xs) {
0.7*dnorm(xs, mi2) /
(0.3*dnorm(xs, mi1) +
0.7*dnorm(xs, mi2))
}Оцену параметра \(\mu=(\mu_1, \mu_2)^T\) ћемо добити методом максималне веродостојности и то је М корак алгоритма.
\[\begin{split} \log L(\mu; {\bf x}, {\bf w})& =\log(\prod_{w_i=0}f_{\mu_1}(x_i)\prod_{w_i=1}f_{\mu_2}(x_i)) \\ &=\sum_{w_i=0}\log f_{\mu_1}(x_i)+\sum_{w_i=1}\log f_{\mu_2}(x_i) \\ &= \sum_{i=1}^n((1-w_i)\log f_{\mu_1}(x_i)+w_i\log f_{\mu_2}(x_i)), \end{split}\]
\[\widehat \mu = \arg \max_{\mu} \log L(\mu; {\bf x}, {\bf w}).\]
Сада ћемо имплементирати цео поступак.
logLxw <- function(mi1, mi2, xs, ws) {
sum((1-ws)*log(dnorm(xs, mi1)) +
ws*log(dnorm(xs, mi2)))
}
maximization <- function(Ew, xs, mi_0) {
nlm(function(mi) -logLxw(mi[1], mi[2], xs, Ew), mi_0)$estimate
}EM_estimate_mix <- function(mi_0, xs, tol=1e-6, maxiter=100) {
iter <- 0
Ew <- expectation(mi_0[1], mi_0[2], xs)
mi <- maximization(Ew, xs, mi_0)
print(mi)
while(any(abs(mi - mi_0) > tol)){
if(iter == maxiter)
break
mi_0 <- mi
Ew <- expectation(mi_0[1], mi_0[2], xs)
mi <- maximization(Ew, xs, mi_0)
print(mi)
iter <- iter + 1
}
mi
}Применом алгоритма, за почетну вредност параметра \(\mu=(0,1)^T\), добија се
ocene=EM_estimate_mix(c(0,1), uzorak)## Warning in nlm(function(mi) -logLxw(mi[1], mi[2], xs, Ew), mi_0): NA/Inf
## replaced by maximum positive value
## Warning in nlm(function(mi) -logLxw(mi[1], mi[2], xs, Ew), mi_0): NA/Inf
## replaced by maximum positive value## [1] -0.7072833 2.5826884## Warning in nlm(function(mi) -logLxw(mi[1], mi[2], xs, Ew), mi_0): NA/Inf
## replaced by maximum positive value## [1] -0.9533337 3.0170478
## [1] -0.9193623 3.0513780
## [1] -0.9056396 3.0583254
## [1] -0.9013947 3.0603282
## [1] -0.900104 3.060931
## [1] -0.8997121 3.0611137
## [1] -0.8995932 3.0611691
## [1] -0.8995571 3.0611859
## [1] -0.899547 3.061191
## [1] -0.8995424 3.0611940
## [1] -0.8995415 3.0611950Две последње вредности представљају оцене \(\widehat \mu_1\) и \(\widehat\mu_2\)