개시하다

이것은 지바대학 노스페이어유한공사천구보입니다. 이번에 저는 EM 알고리즘.EM 알고리즘이 데이터의 가능성을 최대화하는 알고리즘이라고 설명합니다.

수법의 설명

문제 설정

관측 데이터 벡터는 $y$입니다. $z$의 확률 변수와 파라미터 벡터가 없는 $\theta$의 외곽 밀도 함수가 $f(y\mid\theta)$매개 변수 $\theta의 추정 방법으로 $\logf(y\mid\theta)달러의 대수 유사 함수를 최대화했다고 가정합니다.즉, 가장 비슷한 방법을 고려한다. $f(y\mid\theta), $y와 $z의 동시 밀도 $f(y,z\mid\theta)=f(y\mid,\theta)f(z\mid\theta)$, $y의 주변 밀도를 계산할 때의 포인트 평가
$$
f(y\mid\theta) =\int f(y,z\mid\theta) dz
$$
이 경우 사용하기 편한 것은 지금부터 해설되는 EM 알고리즘이다.

EM 알고리즘의 단계

EM 알고리즘은 expectation-maximization 알고리즘의 약자로, 기대치를 계산하는 E-step과 최대화하는 M-step을 번갈아 가며 수렴될 때까지 반복 계산하는 알고리즘이다.

E-step: $k회 중복된 매개 변수 평가치를 $\theta^{(k)}달러로 설정하여 다음 기대치 계산

\begin{align}
Q(\theta, \theta^{(k)}) &= E_{\theta^{(k)}} \left[ \log f(y,z \mid \theta) \mid y \right] \\
&= \int \log \{ f(y,z \mid \theta) \} \ f(z \mid y, \theta^{(k)}) dz.
\end{align}

하지만 $E{theeta^ {(k)} [\cdot\midy] 달러는 주어진 참값이 $\theta^ {(k)} 달러라는 조건부 분포의 기대치 연산을 나타낸다.

M-step: E-step에서 평가한 $Q(\theeta,\theeta ^((k)})는 $\theta 달러를 최대화하고, 최대 $\theta 달러의 값을 $\theeta ^ {k+1)} 달러로 설정합니다.

E-step과 M-step의 교체를 수렴조건($\theeta^{(k)}달러의 값이 거의 움직이지 않을 때까지)을 충족시킬 때까지 반복하고, 최종적으로 얻은 $\theeta^{(k)}는 $\theta$를 추정값으로 하는 것이 EM 알고리즘이다. 여기서 $f(y\mid\theta)가 단봉 유사함수라면 EM 산법으로 얻은 $\the^{(k)달러는 $f(y\mid\theta)}이다.수렴을 최대화하는 최대 유사 추정치를 표시합니다.
EM 알고리즘의 매번 교체에서 가장 유사한 외곽 대수를 최대화하는 것은 $\logf(y\mid\theta)$(전체 데이터의 대수는 비슷하고complete data likelihood)이지만 $때문에이것은 $z\midy]의 값을 조건부 기대값으로 대체하는 그림입니다. 엄밀히 말하면 $E[z\midy]의 전체 데이터의 대수 유사성을 대입하는 것이 아니라 전체 데이터의 대수 유사도를 $z\midy달러의 분포로 기대치를 얻는 것입니다.

EM 알고리즘 적용

EM 알고리즘은 데이터 분석을 측정하는 데 사용될 뿐만 아니라 잠재적 변수를 도입하여 완전한 데이터 가능성을 쓰기 쉬운 모델에서 가장 유사한 방법도 사용할 수 있다. 전형적인 예는 혼합 분포의 추정, 잠재적 변수의 벡터 $z달러를 도입하여 모든 관측이 혼합 분포의 어느 component 표본에서 왔는지 표시하는 것이다.이 $z$를 측정 데이터로 EM 알고리즘을 적용합니다. 혼합 분포 중의 EM 알고리즘은 많은 설명서에 소개되어 있으며, 본고는 선형 혼합 모델이 미지의 매개 변수에 대한 응용 실례를 소개하고자 합니다.

선형 블렌드 모델에 적용

이전 기사에서 소개한 선형 혼합 모델의 일반적인 형식을 고려해 보자.
$$
y = X\beta + Zb +\varepsilon,\quad b\sim\mathrm{N}_m(0,G),\quad\varepsilon\sim\mathrm{N}_n(0,R),\tag{1}
$$
이 모델의 알 수 없는 매개 변수는 EM 알고리즘으로 추정되며, 이 EM 알고리즘은 변수 효과 $b$$에 적용됩니다. 여기서 간단하게 보기 위해 $G와 $R달러에 포함된 분산 성분을 모두 알고 있으며, 알 수 없는 매개 변수는 $\beta$에 불과합니다. 물론 분산 성분이 알 수 없더라도EM 알고리즘으로 회귀계수와 방차분량을 추정할 수 있다. 최대화된 목표 함수를 바탕으로 한 완전한 데이터의 대수 유사도는

\begin{align}
\log f(y,b \mid \beta) =& \ \log f(y \mid b,\beta) + \log f(b) \\
=& -{n \over 2}\log(2\pi) -{1 \over 2}\log |R| - {1 \over 2} (y - X\beta - Zb)^\top R^{-1}(y - X\beta - Zb) \\
& \ -{m \over 2}\log(2\pi) - {1 \over 2}\log|G| - {1 \over 2}b^\top G^{-1}b
\end{align}

$(k+1) 두 번째 E-step에서 $Q(\beta,\beta^{(k)})=E{beta^ {(k)} {\log f(y, b\mid\beta)\midy] 달러 확정, M-step은 $Q(\beta,\beta^ {(k)}) 달러를 $\beta 달러로 최대화하고,
$$
{d\over d\beta}Q(\beta,\beta^{(k)}) = -(X^\top X)\beta + X^\top (y - ZE_{\beta^{(k)}}[b\mid y])
$$
따라서 실제로 E-step에서는 $b\midy 달러 분포의 회전 모멘트만 알면 충분합니다.
$$
E_{\beta^{(k)}}[b\mid y] = GZ^\top\Sigma^{-1}(y - X\beta^{(k)}) =:\tilde{b}^{(k)}
$$
그러나 달러\Sigma=ZGZ^\top+R달러는 $y$의 주위에 분산되어 분포되어 있기 때문에 $(k+1)달러의 M-step에서 가장 좋은 해결 방안은
$$
\beta^{(k+1)} = (X^\top X)^{-1} X^\top (y - Z\tilde{b}^{(k)})
$$
네.

수치 시험

(1) 모델의 특수한 형식으로 다음과 같은 모델을 고려한다.
$$
y_{ij} =\beta_0 + x_{ij}\beta_1 + b_i +\varepsilon_{ij},\quad b_i\sim\mathrm{N}(0,\tau^2),\quad\varepsilon_{ij}\sim\mathrm{N}(0,\sigma^2),\quad (i=1,\dots,m,\j=1,\dots,n_i)
$$
이 모델에서 우리는 $\ta^2=0.5^2,\sigma^2=1^2달러를 이미 알고 있으며 EM 알고리즘으로 $(\beta0,\beta1)=(1,2)달러를 추산한다. 다음은 R의 코드 예시이다.

library(Matrix)

RR <- 100  # 実験の繰り返し数
m  <- 30
ni <- 5
n  <- m * ni

# 真のパラメータの値
tau2 <- 0.5^2
sigma2 <- 1^2
beta <- c(1,2)

# 行列XとZの作成
X <- cbind(1, rbinom(n, 1, 0.3))
i <- 1:n
j <- rep(1:m, each = ni)
Z <- sparseMatrix(i, j, x = 1, dims = c(n,m))

# Sigmaとその逆行列の計算
Sigma <- tau2 * Z %*% t(Z) + sigma2 * diag(n)
Sigma_inv <- solve(Sigma)

res <- data.frame(beta0 = rep(0,RR), beta1 = rep(0,RR), iter = rep(0,RR))
# 実験の繰り返し
for(rr in 1:RR){
  # 乱数の発生
  b <- rnorm(m, sd = sqrt(tau2))
  y <- X %*% beta + Z %*% b + rnorm(n, sd = sqrt(sigma2))

  d <- 1
  iter <- 0
  beta_new <- c(0,0)  # 初期値
  # EMの繰り返し
  while(d > 0.001){
    iter <- iter + 1
    beta_old <- beta_new  # 1期前のM-stepで求めたbeta^kの値

    # E-step
    b_tilde <- tau2 * t(Z) %*% Sigma_inv %*% (y - X %*% beta_old)

    # M-step
    beta_new <- solve(t(X) %*% X) %*% t(X) %*% (y - Z %*% b_tilde)

    # 収束判定
    d <- sqrt( sum( ( beta_new - beta_old )^2 ) )
  }
  # 結果の格納
  res[rr,1:2] <- beta_new
  res[rr,3] <- iter
}

실험을 100번 반복한 결과 데이터 프레임res에 저장되었다. res의 1열과 2열은 $(\beta 0,\beta 1)의 추정치를 저장하고, 3열은 EM 알고리즘을 여러 번 수렴하여

d = \| \beta^{(k+1)} - \beta^{(k)} \|

결과는 다음과 같습니다.

> mean(res$beta0)
[1] 1.00784
> mean(res$beta1)
[1] 2.025926
> var(res$beta0)
[1] 0.0196979
> var(res$beta1)
[1] 0.03556165
>
> table(res$iter)

12 13 
26 74 

실제 매개변수 값은 $(\beta 0,\beta 1)=(1,2)달러이기 때문에 합리적인 결과를 얻었다고 할 수 있습니다.

끝말

주식회사 노스페어는 통계학 각 분야에 전문적으로 종사하는 연구원으로 통계 고문과 상업 데이터에 대한 분석은 문의주식회사하십시오.