2. Markov Chain Monte Carlo

 1. Bayesian 추론에서 MCMC를 왜 쓰는가?

   Markov Chain Monte Carlo는 각각의 시도가 이전 시도에 영향을 받지 않는(memoryless, Markov Chain) 랜덤서치(Monte Carlo)를 통하여 답을 구하는 방식이다.
   1번 포스트에서 썼던 예제는 prior belief와 posterior가 같은 형태의 분포를 갖는 conjugate prior형태이므로 닫힌 해를 이용하여 posterior를 구할 수 있었다. 하지만 대부분의 문제는 $\theta$로 표현하는 인수의 수가 많고, conjugate prior가 아닌 경우가 많아 닫힌 해를 이용하기가 힘들다. 다시 한 번 Bayes Rule의 식을 보자.
$$
P(\theta|D)=\frac{P(D|\theta)P(\theta)}{P(D)}
$$
   위의 식에서 보면 posterior의 $\theta$를 계산하기 위해서는 $P(D)$를 계산해야 하고, 이를 위해 다음의 적분을 해야 한다.
$$
P(D) = \int_{\Theta} P(D, \theta)d\theta
$$
   하지만 인수를 구성하는 $\theta$의 수가 많아질 경우 앞서 설명한대로 닫힌 해를 이용한 적분이 매우 어렵고, 그로인해 닫힌해를 가지고 posterior의 $\theta$를 추정하기가 힘들다.
 그러므로 이런 경우 posterior의 $\theta$를 구하기 위해 MCMC를 적용한다. 물론 여기서 위 적분을 직접 구하는 것은 아니다.

2. MCMC 알고리즘

   MCMC에 속하는 알고리즘은 다음과 같다. Metropolis, Metropolis-Hastings, Gibbs Sampler, Hamiltonian MCMC, No-U-Turn Sampler(NUTS), 그 외.
   이 챕터에서는 가장 오래되고 널리 쓰이는 Metropolis를 우선 설명하고, 이후 챕터에서 다른 알고리즘을 소개한다.

3. Metropolis 알고리즘

   일반적으로 MCMC 알고리즘은 다음의 순서를 따른다.

1. 현재 인수에서 시작($\theta_{c}$)
2. 새로운 인수로의 점프를 제안($\theta_{n}$)
3. 이전 정보와 데이터를 이용하여 확률적으로 새로운 인수의  허용/거부를 결정
4. 허용을 결정하면 새로운 인수로 업데이트 후 1번부터 반복
5. 거부를 결정하면 인수를 업데이트하지 않고 1번부터 반복
6. 지정된 횟수만큼의 시도를 한 후 허용을 결정한 인수를 모두 모음

   위의 과정에서 '어떻게 점프를 하는가'와 '어떻게 점프 여부(허용/거부)를 결정하는가'가 알고리즘마다 다르다.

   Metropolis 알고리즘은 새로운 인수를 찾아 제안할 때, 인수에 대해 가정한 분포와 상관없이 '정규 분포'를 따라 다음 인수를 제안한다. 그래서 다음 제안을 위한 $\mu$와 $\sigma$가 있는데, $\mu$는 기존의 $\theta_{c}$로 정한다. 여기서 $\sigma$를 어찌 지정하는가에 따라 수렴에 영향을 미친다. $\sigma$가 큰 값일 경우 더 넓은 영역을 살펴볼 수 있으나 허용 확률이 높은 부분을 놓치고 지나갈 수 있고, 값이 작을 경우 큰 영역을 살피는데 시간이 오래 걸리므로 수렴 시간이 더 길어질 가능성이 있다.

   이렇게 새로운 인수에 대한 제안이 만들어지면 이제 이 제안을 받아들일지 여부를 결정해야 한다. 이를 위해 가정한 $\theta$의 분포에 따라 $P(\theta_{n})$와 Likelihood $P(D|\theta)$를 계산한다. 그리고 다음의 비율을 구한다.

$$
p = \frac{ P(D|\theta_{n})P(\theta_{n}) }{ P(D|\theta_{c})P(\theta_{c}) }
$$

   그리고 $[0, 1]$에 속하는 균일 분포의 난수를 만들고, 만들어진 수가 $[0, p]$안에 들어오면 이 인수를 허용하고, 아니면 거부한다. 

   결국 Bayesian 추론에서는 새로운 인수로 만들어진 posterior $P(\theta_{n}|D)$를 시뮬레이션을 통하여 구해가는 과정인 것이다. 

   실제로 앞서 제시한 Coin toss 문제를 Metropolis를 이용하여 푼 예제를 보자. 여기서 50번의 동전 던지기에서 10번의 앞면이 나왔을 경우 posterior를 구하는 예제이다.

   위 그림에서 녹색선으로 그려진 부분이 닫힌 해를 이용한 답이고, 붉은 박스로 그려진 히스토그램이 Metropolis를 이용한 결과이다. 나름 비슷한 모양으로 답을 찾은 것을 볼 수 있다. 하지만 이 시뮬레이션의 경우 $\theta$의 수가 하나인데도, 이 결과를 얻기 위해 십만 번의 시뮬레이션을 수행해야 했다. 그러니 인수의 수가 많아질수록 더 큰 컴퓨팅 파워를 필요로 하게될 것이다.

   위 그래프는 PyMC3라는 라이브러리를 이용하여 python으로 구현한 내용이다. 코드를 다음에 첨부한다.

 
# mcmc_coin.py
# Markov Chain Monte Carlo Test for coin flipping

import matplotlib.pyplot as plt
import numpy as np
import pymc3
import scipy.stats as stats

# gg plot style의 plotting
plt.style.use("ggplot")

# parmeters
n = 50 # num of coin flip
z = 10 # num of heads
alpha = 12 # alpha for beta-distribution
beta = 12 # beta for beta-distribution
alpha_post = 22 # known updated parameter
beta_post = 52 # known updated parameter

# Metropolis algorithm에서 수행할 시뮬레이션 횟수
iterations = 100000

# PyMC3를 사용하여 model context 만들기
basic_model = pymc3.Model()
with basic_model:
    # Beta Distribution을 이용하여 prior belief 만들기
    theta = pymc3.Beta("theta", alpha=alpha, beta=beta)

    # Bernoulli likelihood 함수 지정
    y = pymc3.Binomial("y", n=n, p=theta, observed=z)

    #Metrobpolis Algorithm을 이용하여 MCMC수행
    # Maximum A Posteriori (MAP) 최적화를 initial value를 찾는데 사용
    start = pymc3.find_MAP()
    step = pymc3.Metropolis()

    # trace 계산
    trace = pymc3.sample(iterations, step, start, random_seed=1, progressbar=True)

# posterio histogram 그리기
bins = 50
plt.hist(trace["theta"], bins, histtype="step", normed=True, label="Posterior (MCMC)", color="red")

# analytic 찍기
x = np.linspace(0, 1, 100)
plt.plot(x, stats.beta.pdf(x, alpha, beta), "--", label="Prior", color="blue")
plt.plot(x, stats.beta.pdf(x, alpha_post, beta_post), label="Posterior (Analytic)", color="green")

#graph label
plt.legend(title="Parameters", loc="best")
plt.xlabel("$\\theta$, Fairness")
plt.ylabel("Density")
plt.show()

MCMC에 관해 조금 더 쉽게 설명된 내용은 http://twiecki.github.io/blog/2015/11/10/mcmc-sampling/   를 참고하기 바란다.