Monte Carlo Optimization

Monte Carlo 적분을 이전 포스트에서는 얘기를 했는데 이제는 이 적분값을 구하는 것뿐만 아니라 최적화해야하는 상황이라고 생각해보자. 

지금까지 Monte Carlo integration을 통해 직접 계산하기 어려운 적분을 표본평균으로 근사하는 방법을 살펴보았다.

그런데 통계에서는 적분값 자체를 구하는 데서 끝나지 않고, 그 적분을 포함하는 함수를 최대화하거나 최소화해야 하는 경우가 많다.

이때 Monte Carlo 방법은 최적화에 필요한 계산을 돕는 도구가 된다.

 

결측자료가 있는 상황이 대표적인 예다.

일부 자료가 관측되지 않았다면 (Censored data), 관측된 자료와 모수 아래에서 가능한 결측값들을 고려해야 한다.

특정 값 하나로 빈칸을 채우기보다는, 가능한 값들이 얼마나 그럴듯한지에 따라 가중하여 평균을 내는 것이다.

이러한 과정은 수학적으로 적분이나 기댓값의 계산으로 나타난다.

 

문제는 모형이 복잡해질수록 이 계산을 직접 수행하기 어렵다는 점이다.

여기서 Monte Carlo integration을 활용하면, 가능한 결측자료를 반복적으로 생성하고 그 결과를 평균하여 필요한 기댓값을 근사할 수 있다.

그리고 이 근사를 이용해 모수를 갱신하는 최적화를 수행한다.

대표적인 방법인 Monte Carlo EM은 결측자료에 대한 조건부 기댓값을 시뮬레이션으로 근사하는 단계와, 이를 최대화하여 모수를 갱신하는 단계를 번갈아 진행한다.

 

이처럼 Monte Carlo integration이 계산하기 어려운 값을 근사하는 방법이라면, 이를 활용한 최적화는 그 근사를 바탕으로 더 나은 모수를 찾아가는 과정으로 이해할 수 있다.

 

그럼 구체적으로 어떻게 풀어나갈 것인지 살펴보자.

Censored data에 대한 문제를 보기 전에 하나 정의를 하고 들어가자.

Missing data model
$$g(x \mid \theta) = \int_z f(x, z \mid theta)dz$$

 

그럼 이제 likelihood 부분부터 확인해보겠다.

Ex. Censored Data Likelihood
Observe $Y_1, \ldots, Y_n \overset{iid}{\sim} f(y-\theta)$
Let $\mathbf{y} = (y_1, \ldots, y_m)$: Uncensored, $\mathbf{z} = (y_{m+1}, \ldots, y_n)$: Censored

 

위와 같이 uncersored된 것은 이미 우리가 관측했다고 생각하고 뒤에 censored된 데이터는 아직 우리가 관측을 하지 않은 값들이라고 생각하자.

감이 안 잡힐 수도 있는데 쉽게 생각하면

우리가 1시간 단위로 10시간 후에 전등 10개가 꺼질지를 확인하고 있었는데, 10시간 후에 관찰을 완료한 결과 8개의 전등은 꺼져있는데

2개의 전등이 꺼져있지 않은 것이다. 자 그렇다면 이 전등 2개를 배제하고 데이터를 분석하거나 추정값을 구할 것인가 하면 조금은 위험한 선택일 수 있다는 것이다.이러한 상황에서 EM algorithm을 통해 이 2개의 전등을 censored된 데이터로 보고 parameter값에 대해서 접근을 해보겠다는 것이다.

 

이를 위해 Censored된 데이터와 Uncensored된 데이터의 Likelihood를 정의하고 시작하겠다.

The observed likelihood function is 
\begin{align*}
g(\mathbf{x}\mid\theta) = L(\theta|\mathbf{y}) &= \int \int \cdots \int \prod^n_{i=1}f(y_i \mid \theta) dy_{m+1} \cdots dy_{n} \\
&= \prod^m_{i=1} f(y_i - \theta)[1-F(a-\theta)]^{n-m}
\end{align*}
The complete likelihood function is 
$$f(\mathbf{x}, \mathbf{z} \mid \theta)L^c(\theta \mid \mathbf{y},\mathbf{z}) = \prod^m_{i=1}f(y_i -\theta) \prod^n_{i=m+1}f(z_i - \theta)$$
Define 
$$k(\mathbf{z}|\mathbf{x}, \theta) = \frac{f(\mathbf{x}, \mathbf{z} \mid \theta)}{g(\mathbf{x}\mid\theta)}$$

 

자 그럼 우리가 현재 원하는 것이 무엇이냐

$$\hat{\theta} = \operatorname{argmax}_{\theta}L(\theta \mid \mathbf{y})$$

 

저 Likelihood를 가장 크게 만드는 $\theta$를 찾고 싶은 것인데 그러나 아시다시피 우리한테는 observed likelihood만을 가지고 있을 뿐, complete likelihood를 가지고 있지 않다. 어떤 방법으로 접근을 하는게 좋을까? 당연히 다 알고있는 complete likelihood를 사용해서 접근하면 좋아보인다. 

그렇기 때문에 이전에 $k(\mathbf{z} \mid \mathbf{x}, \theta)$를 정의했는데 이제 EM algorithm을 어떤 식으로 접근을 해야하는지 하나씩 derivation 해보겠다.

 

\begin{align*} \hat{\theta} &\in \underset{\theta}{\operatorname{argmax}}\, \log g(\mathbf{x} \mid \theta). \end{align*} \begin{align*} f(\mathbf{x}, \mathbf{z} \mid \theta) &= g(\mathbf{x} \mid \theta) k(\mathbf{z} \mid \mathbf{x}, \theta), \\ \log g(\mathbf{x} \mid \theta) &= \log f(\mathbf{x}, \mathbf{z} \mid \theta) - \log k(\mathbf{z} \mid \mathbf{x}, \theta). \end{align*} \begin{align*} \log g(\mathbf{x} \mid \theta) &= \operatorname{E}_{\mathbf{z} \mid \mathbf{x}, \theta^t} \left[\log f(\mathbf{x}, \mathbf{z} \mid \theta)\right] \\ &\quad - \operatorname{E}_{\mathbf{z} \mid \mathbf{x}, \theta^t} \left[\log k(\mathbf{z} \mid \mathbf{x}, \theta)\right]. \end{align*} \begin{align*} \theta^{t+1} &\in \underset{\theta}{\operatorname{argmax}}\, \operatorname{E}_{\mathbf{z} \mid \mathbf{x}, \theta^t} \left[\log f(\mathbf{x}, \mathbf{z} \mid \theta)\right] \\ &= \underset{\theta}{\operatorname{argmax}}\, \int \log f(\mathbf{x}, \mathbf{z} \mid \theta) k(\mathbf{z} \mid \mathbf{x}, \theta^t) \,d\mathbf{z}. \end{align*}

 

왜 이게 성립하는지 확인해보겠다.

WTS: $\log L(\hat{\theta} \mid \mathbf{x}) \geq \log L(\theta_0 \mid \mathbf{x})$
\begin{align*} \hat{\theta} &\in \underset{\theta}{\operatorname{argmax}}\, \operatorname{E}_{\mathbf{z}\mid\mathbf{x},\theta_0} \left[\log L^c(\theta\mid\mathbf{x},\mathbf{z})\right]\end{align*} \begin{align*} \log L(\theta_0\mid\mathbf{x}) &=\operatorname{E}_{\mathbf{z}\mid\mathbf{x},\theta_0} \left[\log L^c(\theta_0\mid\mathbf{x},\mathbf{z})\right]\\ &\quad-\operatorname{E}_{\mathbf{z}\mid\mathbf{x},\theta_0} \left[\log k(\mathbf{z}\mid\mathbf{x},\theta_0)\right],\\ \log L(\hat{\theta}\mid\mathbf{x}) &=\operatorname{E}_{\mathbf{z}\mid\mathbf{x},\theta_0} \left[\log L^c(\hat{\theta}\mid\mathbf{x},\mathbf{z})\right]\\ &\quad-\operatorname{E}_{\mathbf{z}\mid\mathbf{x},\theta_0} \left[\log k(\mathbf{z}\mid\mathbf{x},\hat{\theta})\right]. \end{align*} \begin{align*} &\operatorname{E}_{\mathbf{z}\mid\mathbf{x},\theta_0} \left[\log L^c(\hat{\theta}\mid\mathbf{x},\mathbf{z})\right] -\operatorname{E}_{\mathbf{z}\mid\mathbf{x},\theta_0} \left[\log L^c(\theta_0\mid\mathbf{x},\mathbf{z})\right]\geq 0 \cdots \text{trivial}. \end{align*} \begin{align*} &\operatorname{E}_{\mathbf{z}\mid\mathbf{x},\theta_0} \left[\log k(\mathbf{z}\mid\mathbf{x},\hat{\theta})\right] -\operatorname{E}_{\mathbf{z}\mid\mathbf{x},\theta_0} \left[\log k(\mathbf{z}\mid\mathbf{x},\theta_0)\right]\\ &=\int \log\frac{k(\mathbf{z}\mid\mathbf{x},\hat{\theta})} {k(\mathbf{z}\mid\mathbf{x},\theta_0)} k(\mathbf{z}\mid\mathbf{x},\theta_0)\,d\mathbf{z}\\ &\leq\log\int \frac{k(\mathbf{z}\mid\mathbf{x},\hat{\theta})} {k(\mathbf{z}\mid\mathbf{x},\theta_0)} k(\mathbf{z}\mid\mathbf{x},\theta_0)\,d\mathbf{z} \because \text{Jensen's Inequality }\\ &=\log 1=0. \end{align*} \begin{align*} &\log L(\hat{\theta}\mid\mathbf{x})-\log L(\theta_0\mid\mathbf{x})\\ &=\operatorname{E}_{\mathbf{z}\mid\mathbf{x},\theta_0} \left[\log\frac{L^c(\hat{\theta}\mid\mathbf{x},\mathbf{z})} {L^c(\theta_0\mid\mathbf{x},\mathbf{z})}\right]\\ &\quad-\operatorname{E}_{\mathbf{z}\mid\mathbf{x},\theta_0} \left[\log\frac{k(\mathbf{z}\mid\mathbf{x},\hat{\theta})} {k(\mathbf{z}\mid\mathbf{x},\theta_0)}\right]\geq 0. \end{align*} \[ \boxed{\log L(\hat{\theta}\mid\mathbf{x}) \geq\log L(\theta_0\mid\mathbf{x}).} \]

 

그렇게 해서 EM algorithm은 다음과 같아진다. 

 

EM algorithm
1. Compute $\operatorname{E}_{\hat{\theta_{m}}}[\log L^c(\theta \mid \mathbf{x}, \mathbf{z})]$, where the expectation is w.r.t. $k(\mathbf{z}\mid \hat{\theta}_m, \mathbf{x})$.
2. Maximize $\operatorname{E}_{\hat{\theta_{m}}}[\log L^c(\theta \mid \mathbf{x}, \mathbf{z})]$ in $\theta$ and take $$\hat{\theta}_{(m+1)} = \operatorname{argmax}_{\hat{\theta}_(m)}\operatorname{E}[\log 
L^c(\theta \mid \mathbf{x}, \mathbf{z}
)]$$