그럼 MCMC는 도대체 뭔가?

 

그 동안 MCMC 이야기하면서 그럼 도대체 MCMC가 뭐냐고 답답하게 생각하셨을 분들이 있었을겁니다. 그래서 간단히 MCMC를 소개하겠습니다.

 

목적 확률분포를 f(x)라 하고, proposal 분포를 g(x)하고 하죠.

 

x(n)는 n번 째 f(x)에서 뽑힌 표본값이고 x*는 그 다음 뽑힐 가능성이 있는 표본 값이라고 하죠. 이 x*이 어떤 확률값을 가지고 다음 x(n+1)일 될 수도 있고, 아니면 그냥 기각될 수 있습니다. 기각될 경우 다음 가능성이 있는 x*를 g(x*|x)에 뽑고 위의 방법으로 받아 들이거나 아니고 또 기각을 합니다. U는 균등분포, 즉 랜덤넘버입니다.

 

x(n+1)=x*           if U<min(1, )

          =x(n)         otherwise: 즉 x*을 기각하는 것입니다.

 

복잡하죠. 이 알고리즘이 기본이 되고 여기서 수 많은 변형들이 나옵니다. 그래서 설명을 늦춘겁니다. 이 MCMC의 역사적 기원을 보면

 

옛날에 Metropolis와 동료들이 이 알고리즘을 개발합니다. 이 알고리즘을 메트로폴리스 알고리즘이라 합니다. 이는 위에 식에서

 

 

일 경우입니다. 정규분포가 대표적으로 이런 조건을 만족시킵니다. f(x)는 Boltzman 분포이고요. 제 기억으로는요. 이 메트로폴리스 알고리즘은 통계역학이나 열역학에서 물질의 입자의 상태를 표현하는 알고리즘으로 물리학에서 이런 쪽 공부하는 사람에게만 일부 알려져 있었습니다. 이 알고리즘이 Hasting과 동료들이 위의 식으로 일반화해서 이젠 Metropolis-Hastings(MH) 알고리즘으로 알려져 있습니다.

 

이게 크게 관심을 받게 된 논문이 이미지 복원 논문인 German & Gemrman 논문입니다. 베이지안 추론시 복잡한 계산을 이 알고리즘의 변형을 이용하여 풀어냅니다. 이 German & German 알고리즘을 Gibbs Sampler라 합니다. 이 논문만 보면 German & German은 Metropolis의 알고리즘은 알고 있었지만 일반화된 위의 MH 알고리즘은 모르고 있었던 것 같습니다. 나중에 이 깁스 샘플러가 이 MCMC의 특수한 형태임이 밝혀집니다.

 

통계학에서는 MLE 이후 모수 추정치를 계산하는데는 큰 발전이 없었습니다. 그러나다 Dempster와 동료들이 EM 알고리즘을 발표합니다. 앞에서 EM 알고리즘을 소개하였죠. 그러다가 Tanner & Wong이 EM 알고리즘을 도입하는 대신 조건부 확률에서 표본을 추출해 반복해서 계산하는 알고리즘을 발표합니다. data augmentation 또는 imputed 알고리즘이라 이야기를 했습니다. 지금 제 기억으로는 그렇습니다. 나중에 다시 확인해 보겠습니다. 이게 나중에 German & German의 깁스 샘플러나는게 밝혀졌습니다. 이후 통계학이나 기타 분야에서 통계 방법론을 쓰는 사람들이 베이지안의 복잡한 계산에서 이 알고리즘을 쓰는 것이 대 유행으로 번졌습니다. 최근에는 인공지능에 응용하려고 하고요.

    

 

이 MCMC의 맛을 좀 보려는 사람은 Navarro & Perfors의 MCMC 소개 글을 보시기 바랍니다. 인터넷에 돌아 다닙니다. 이 사람들이 예로 들은 확률분포입니다.

복잡한 것 같지만 별거 아닙니다. 그냥 적분가능한 양의 함수를 만들고 아래 분포에 함수의 적분값으로 나눈 것입니다. 이렇게 적분 값으로 나누면 이 함수의 적분값은 1, 즉 확률분포가 됩니다. 이걸 정규화한다고 합니다. 그리고 아래 분모는 적분값, 그냥 숫자이기 때문에 무시해도 됩니다. 이 확률분포의 그래프는 다음과 같습니다.

    이 분포의 최빈, 중앙값, 평균, 또는 다양한 기대값을 구하고 싶을 때 이 분포 모양을 하는 표본을 뽑아서 한다는 것이거든요.

 

그리고 이 사람들이 proposal 분포로 정규분포를 사용했습니다. 왜냐면 정규분포는 g(x|y), 즉 조건부 분포도 정규분포이고 이미 널리 알려져 있고, 정규분포에서 표본을 뽑는 것은 웬만한 통계관련, S, R, matlab, Gauss, Eview 등 행렬 계산 하는 프로그램에 그냥 명령문만 치면 나옵니다.

 

그래서 이 proposal 정규분포에서 몇가지 값을 주어 MCMC를 했습니다. 이 proposal 정규분포의 초기값으로 평균은 -1로 주고 표준편차인 시그마를 1, 0.025, 50으로 세 가지 경우를 돌렸습니다. 그 결과가 아래 그림입니다.

    그림을 보면 맨 왼쪽인 시그마가 1인 경우는 표본이 매우 잘 뽑혔습니다. 목적 확률분포의 모양과 매우 흡사하죠. 그러나 두 번째인 정규분포가 매우 집중되어 있는 시그마가 0.025는 -1 근처에서 많이 몰려 있습니다. 즉 local minimum이나 maximum에 빠진 모양이라 비슷합니다. 마지막으로 시그마가 50인 매우 평평한 정규분포를 제안분포로 했을 경우는 목적 분포인 f(x) 전 지역을 카바하지만 accept 되는 비율이 너무 낮습니다.

 

 

그러나 이 MCMC을 계속 시도하면, 즉 n을 엄청나게 많이 하면 아래 그림이 됩니다. 효율성에는 좀 문제가 되지만 어떤 경우라도 추출된 표본값이 f(x)라 비슷하게 갑니다.

 

    여기서 결론은 n을 크게 하면 원래 목적 분포 모양을 하는 표본을 뽑을 수는 있다, 그러나 제안분포 g(x)을 잘못 설정하면 매우 비효율적이고 추출 회수 n을 작게 하면 처음 결과의 두 번째인 로칼 미니멈이나 맥시멈에 빠질 수가 있다는 것입니다. 그러나 우리가 최초에 어떤 평균이나 그마를 줘야 한다는 것은 모른다는 것이 문제이죠. 그래서 안전을 위해 초반에 나오는 값은 무시하고 어느 정도 진행된 후 나중에 나오는 표본을 사용해야 합니다. 이렇게 초반 표본을 무시하는 것을 burn-in 이라고 합니다.

    

 

이런 작업이 힘들가요? 아닙니다. 저도 이런 일 안하지 십오년이 넘어가지만 하루면 충분히 할 수 있는 작업입니다. 단 R, S, Matlab, Gauss 같은 행렬 조작 프로그램을 능숙하게 작업할 줄 알아야 합니다. 하여간 구글이 생겨나서 공부하는 사람들이 참 편해졌습니다. 제가 공부할 때는 이렇게 쉽게 이해하기 힘들어거든요.