MCMC sampling for dummies 

나는 확률 프로그래밍과 베이지안 통계를 이야기할 때, 보통 추론에 대한 자세한 설명은 얼머무리고 블랙박스로 남겨놓았다. 확률 프로그래밍의 아름다움은 모델을 만들기 위해, 추론이 어떻게 일어나는지 이해할 필요는 없다는 것이다. 물론 이해하는 것이 도움이 된다. 베이지안 통계를 모르지만 알고싶어 하는 Quantopian의 CEO, Fawce 에게 새로운 베이지안 모델을 소개할 때, 그는 내가 얼머무리던 부분에 대해서 묻기 시작했다. “토마스, 추론은 어떻게 동작하는거예요? 사후확률에서 어떻게 이런 마법과 같은 표본을 얻을 수 있죠? “

나는 이렇게 대답했다. “그건 쉬워요, MCMC 알고리즘이 평형 분포가 목표 사후확률 분포와 동일한 가역 마르코프 체인을 만듦으로써 표본을 생성해요. 질문 있나요?" 저 설명은 맞지만, 정말 이해가 되는가? 수학이나 통계를 배울 때 가장 짜증나는 점은 아무도 개념 뒤에 감춰진 직관(보통 매우 쉽다)을 말해주는 사람이 없다는 것이다. 그냥 단지 겁나는 수식만 나열한다. 이 방식이 내가 배웠던 방식이고 나는 유레카를 외치는 순간이 올 때까지, 수없이 많은 시간을 머리를 벽에 박는데 낭비해야 했다. 보통, 수학이 의미하는 것을 일단 복호화 하면, 그것은 두렵운 것도 복잡한 것도 아니었다.

이 블로그 포스트는 MCMC 샘플링 뒤에 숨겨진 직관을 설명할 것이다(더 명확하게는 Metropolis 알고리즘) 수식과 수학적 말보다는 코드 예제를 사용할 것이다. 수학으로 이해하기 전에, 예제로 시작해서 직관을 생각해보는게 더 개인적으로는 더 낫다고 생각한다.

The problem and its unintuitive solution

베이즈 공식을 보자: 
s 
s 는 데이터 s가 주어졌을 때, 모델 파라메터 s 의 확률을 말하고, 그러므로 우리가 원하는 것의 양적 수치이다.이것을 계산하기 위해, 사전확률 s(데이터를 보기 전에 s에 대해서 생각하는 것)을 우도 s(어떻게 데이터가 분포하는지 나타냄)에 곱한다.

분자는 풀기 쉽다. 그러나, 분모를 자세히 보자. s 는 소위 증거라고 불린다. 즉 데이터 s 가 모델에 의해 생성될 증거 모든 가능한 파라메터 값에 대해서 적분함으로써, 이 값을 구할 수 있다. 
s 
이는 베이즈 공식의 가장 큰 난관이다. 가장 간단한 모델에 대해서도, closed-form으로 사후 확률을 구할 수 없다. 다음과 같이 말할지 모르겠다 “좋아, 어떤 것을 정확하게 풀수 없으면, 근사해를 구하면 되지 않겠어? 예를 들어, 사후 확률에서 표본을 수집하면, Monte Carlo 근사를 할 수 있어". 불행히도, 저 분포로부터 직접적으로 샘플링을 하기 위해선, 베이즈 공식을 풀어야할 뿐만 아니라, 그 역도 구해야 하고, 그건 더 어렵다. 그럼 또 다음과 같이 말할 수 있다 “흠, 그럼 대신 사후확률 분포와 일치하는 균형 분포를 가진 마르코프 체인을 만들자".난 단지 농담을 하고 있다. 대부분의 사람들은 저런 미쳐 보이는 말을 하진 않는다.

어떤 것을 계산을 할 수 없다면, 그것에서 표본을 추출할 수 없다. 그러면, 이러한 모든 성질을 가진 마르코프 체인을 만드는 것은 휠씬 더 어렵다. 하지만, 놀라운 통찰은 이 일이 실제로는 매우 쉽고, (몬테카를로 근사를 하는 마르코프 체인을 만드는) Markov chain Monte Carlo이라 불리는 일반적인 종류의 알고리즘들이 존재한다는 것이다.

Setting up the problem

필요한 모듈들을 import하자

%matplotlib inline

import numpy as np
import scipy as sp
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from scipy.stats import norm

sns.set_style('white')
sns.set_context('talk')

np.random.seed(123)

0 중심을 가진 정규분포에서 100개의 점을 샘플링하자. 목표는 평균 mu (표준 편차는 1임을 알고 있다고 가정하자) 의 사후분포를 추청하는 것이다.

data = np.random.randn(20)
ax = plt.subplot()
sns.distplot(data, kde=False, ax=ax)
_ = ax.set(title='Histogram of observed data', xlabel='x', ylabel='# observations');

다음으로, 모델을 정의해야 한다. 이 간단한 경우에서, 이 데이터는 정규 분포라 가정한다. (즉 모델의 우도는 정규분포이다) 아다시피, 정규분포는 두 개의 (하이퍼)파라메터를 가지고 있다. 평균 s 와 표준 편차 s

간소함을 위해, s 임을 가정하고 s에 대한 사후분포를 추론한다고 하자. 추론하고 하는 각 파라메터에 대해서, 사전 확률을 골라야 한다. 간소함을 위해, s에 대한 사전 확률로서 정규 분포를 가정하자. 그러므로, 통계학적으로 말해서 우리의 모델은 다음과 같다. 
s 
What is convenient, is that for this model, we actually can compute the posterior analytically. 
편리한 것은, 이 모델에 대해서, 사후분포를 실재로는 해석적으로 구할 수 있다는 것이다. 왜냐하면, 알려진 표준편차를 가진 정규분포의 우도에 대해서, 정규분포를 따르는 사전확률 s는 conjugate다. (conjugate 는 사후 확률이 사전확률과 같은 분포를 가진다는 의미이다) 그래서, s에 대한 사후 확률이 정규분포임을 안다.

사후확률의 파라메터를 어떻게 계산하는 지는 위키피디아를 쉽게 참고할 수 있다. 
수학적인 유도는 여기를 보라

def calc_posterior_analytical(data, x, mu_0, sigma_0):
    sigma = 1.
    n = len(data)
    mu_post = (mu_0 / sigma_0**2 + data.sum() / sigma**2) / (1. / sigma_0**2 + n / sigma**2)
    sigma_post = (1. / sigma_0**2 + n / sigma**2)**-1
    return norm(mu_post, np.sqrt(sigma_post)).pdf(x)

ax = plt.subplot()
x = np.linspace(-1, 1, 500)
posterior_analytical = calc_posterior_analytical(data, x, 0., 1.)
ax.plot(x, posterior_analytical)
ax.set(xlabel='mu', ylabel='belief', title='Analytical posterior');
sns.despine()

이것은 우리가 알고자하는, 사전 확률 정보를 고려하고, 데이터를 보고난 후의 s의 값들의 확률의 양적 수치를 보여준다. 그러나, 사전확률이 conjugate 가 아니고, 손으로 사후 확률을 구할 수 없다고 가정하자.

Explaining MCMC sampling with code

샘플링 로직을 보자. 맨처음에는, (랜덤으로 선택할 수 있는) 시작 파라메터 지점을 임의로 지정한다.

mu_current = 1.

그 후, 그 지점을 어느곳으로 옮길지 알아야한다. 이것이 마르코프가 하는 일이다. 독자들은 어떻게 옮길 지점을 알지에 대해서 잘 알수도, 모를 수도 있다.

Metropolis sampler는 매우 멍청하고, 단지 현재 mu 값(즉, mu_current)이 중심이고 특정 표준편차(proposal_width)를 가진 정규 분포에서 표본을 추출한다(모델에 대해서 가정한 정규분포와는 관계없음). proposal_width는 얼마나 멀리 새자리로 옮길지를 결정할 것이다.

proposal = norm(mu_current, proposal_width).rvs()

다음에는, 옮긴 장소가 좋은 곳인지 평가한다. 제안된 mu를 가진 정규 분포가 예전 mu보다 데이터를 좀 더 잘 설명한다면, 당연히 그곳으로 가고 싶을 것이다. “데이터를 좀 더 잘 설명"한다는 것은 무슨 의미일까? 제안된 파라메터 값(제안된 mu와 고정된 sigma)을 가진 우도가 주어졌을때 (정규분포) 데이터의 확률을 계산함으로써 적합함을 수치화한다.

scipy.stats.normal(mu, sigma).pdf(data)을 사용해서 각 데이터 포인트에 대한 확률을 계산하고 각 확률을 곱함으로써 이를 구할 수 있다. 즉 우도를 계산한다(보통 log 확률을 쓰지만 여기서는 생략)

likelihood_current = norm(mu_current, 1).pdf(data).prod()
likelihood_proposal = norm(mu_proposal, 1).pdf(data).prod()

# Compute prior probability of current and proposed mu        
prior_current = norm(mu_prior_mu, mu_prior_sd).pdf(mu_current)
prior_proposal = norm(mu_prior_mu, mu_prior_sd).pdf(mu_proposal)

# Nominator of Bayes formula
p_current = likelihood_current * prior_current
p_proposal = likelihood_proposal * prior_proposal

지금까지, 본질적으로 랜덤한 방향으로 움직임을 제안하고 mu_proposal가 mu_current보다 높은 우도를 가지는 경우에만 그 움직임를 받아드리는 언덕 오르기 알고리즘을 시행했다. 결국, 더이상 이동할 수 없는 mu=0 (또는 0에 가까운) 를 얻을 것이다.

그러나, 사후 확률을 얻고 가끔은 다른 방향으로의 움직임을 승인해야 한다. 가장 중요한 비결은 두 확률을 나누는 것이다

p_accept = p_proposal / p_current

승인 확률을 얻었다. p_proposal가 크면, 저 확률은 >1 일것이며, 승인 된다. 그러나, p_current가 크면, 예를 들어 두배 크면, 움직일 확률은 50%가 될 것이다.

accept = np.random.rand() < p_accept

if accept:
    # Update position
    cur_pos = proposal

이 간단한 프로시져가 사후확률에서 표본을 얻을수 있게 해준다.

Why does this make sense?

한 걸음 뒤로 물러나서, 위의 승인 비율이 이 것이 동작하고, 적분을 피할수 있는 이유임을 주목하자. 정규화된 사후확률에 대한 승인률을 구하고 정규화되지 않은 사후확률의 승인률과 같음을 보임으로써 이를 보일 수 있다. s를 현재 지점이고 s를 제안지점이라고 하자. 
s 
말로 하면, 제안된 파라메터 setting의 사후확률을 현재 파라메터 setting의 사후확률로 나누면, 계산할 수 없는 양 s는 소거된다. 그래서, 한 지점에서 full posterior를 다른 지점에서의 full posterior와 실제로 나눈 것과 같은 것임을 직감할 수 있다.

저 방법으로, 상대적으로 높은 사후확률을 가진 지점을 좀 더 자주 방문 할 수 있다.

Putting it all together

def sampler(data, samples=4, mu_init=.5, proposal_width=.5, plot=False, mu_prior_mu=0, mu_prior_sd=1.):
    mu_current = mu_init
    posterior = [mu_current]
    for i in range(samples):
        # suggest new position
        mu_proposal = norm(mu_current, proposal_width).rvs()

        # Compute likelihood by multiplying probabilities of each data point
        likelihood_current = norm(mu_current, 1).pdf(data).prod()
        likelihood_proposal = norm(mu_proposal, 1).pdf(data).prod()

        # Compute prior probability of current and proposed mu        
        prior_current = norm(mu_prior_mu, mu_prior_sd).pdf(mu_current)
        prior_proposal = norm(mu_prior_mu, mu_prior_sd).pdf(mu_proposal)

        p_current = likelihood_current * prior_current
        p_proposal = likelihood_proposal * prior_proposal

        # Accept proposal?
        p_accept = p_proposal / p_current

        # Usually would include prior probability, which we neglect here for simplicity
        accept = np.random.rand() < p_accept

        if plot:
            plot_proposal(mu_current, mu_proposal, mu_prior_mu, mu_prior_sd, data, accept, posterior, i)

        if accept:
            # Update position
            mu_current = mu_proposal

        posterior.append(mu_current)

    return posterior

# Function to display
def plot_proposal(mu_current, mu_proposal, mu_prior_mu, mu_prior_sd, data, accepted, trace, i):
    from copy import copy
    trace = copy(trace)
    fig, (ax1, ax2, ax3, ax4) = plt.subplots(ncols=4, figsize=(16, 4))
    fig.suptitle('Iteration %i' % (i + 1))
    x = np.linspace(-3, 3, 5000)
    color = 'g' if accepted else 'r'

    # Plot prior
    prior_current = norm(mu_prior_mu, mu_prior_sd).pdf(mu_current)
    prior_proposal = norm(mu_prior_mu, mu_prior_sd).pdf(mu_proposal)
    prior = norm(mu_prior_mu, mu_prior_sd).pdf(x)
    ax1.plot(x, prior)
    ax1.plot([mu_current] * 2, [0, prior_current], marker='o', color='b')
    ax1.plot([mu_proposal] * 2, [0, prior_proposal], marker='o', color=color)
    ax1.annotate("", xy=(mu_proposal, 0.2), xytext=(mu_current, 0.2),
                 arrowprops=dict(arrowstyle="->", lw=2.))
    ax1.set(ylabel='Probability Density', title='current: prior(mu=%.2f) = %.2f\nproposal: prior(mu=%.2f) = %.2f' % (mu_current, prior_current, mu_proposal, prior_proposal))

    # Likelihood
    likelihood_current = norm(mu_current, 1).pdf(data).prod()
    likelihood_proposal = norm(mu_proposal, 1).pdf(data).prod()
    y = norm(loc=mu_proposal, scale=1).pdf(x)
    sns.distplot(data, kde=False, norm_hist=True, ax=ax2)
    ax2.plot(x, y, color=color)
    ax2.axvline(mu_current, color='b', linestyle='--', label='mu_current')
    ax2.axvline(mu_proposal, color=color, linestyle='--', label='mu_proposal')
    #ax2.title('Proposal {}'.format('accepted' if accepted else 'rejected'))
    ax2.annotate("", xy=(mu_proposal, 0.2), xytext=(mu_current, 0.2),
                 arrowprops=dict(arrowstyle="->", lw=2.))
    ax2.set(title='likelihood(mu=%.2f) = %.2f\nlikelihood(mu=%.2f) = %.2f' % (mu_current, 1e14*likelihood_current, mu_proposal, 1e14*likelihood_proposal))

    # Posterior
    posterior_analytical = calc_posterior_analytical(data, x, mu_prior_mu, mu_prior_sd)
    ax3.plot(x, posterior_analytical)
    posterior_current = calc_posterior_analytical(data, mu_current, mu_prior_mu, mu_prior_sd)
    posterior_proposal = calc_posterior_analytical(data, mu_proposal, mu_prior_mu, mu_prior_sd)
    ax3.plot([mu_current] * 2, [0, posterior_current], marker='o', color='b')
    ax3.plot([mu_proposal] * 2, [0, posterior_proposal], marker='o', color=color)
    ax3.annotate("", xy=(mu_proposal, 0.2), xytext=(mu_current, 0.2),
                 arrowprops=dict(arrowstyle="->", lw=2.))
    #x3.set(title=r'prior x likelihood $\propto$ posterior')
    ax3.set(title='posterior(mu=%.2f) = %.5f\nposterior(mu=%.2f) = %.5f' % (mu_current, posterior_current, mu_proposal, posterior_proposal))

    if accepted:
        trace.append(mu_proposal)
    else:
        trace.append(mu_current)
    ax4.plot(trace)
    ax4.set(xlabel='iteration', ylabel='mu', title='trace')
    plt.tight_layout()
    #plt.legend()

Visualizing MCMC

샘플링을 시각화하기 위해, 계산된 수치들을 그려볼 것이다. 밑의 각 행은 Metropolis sampler를 1 iteration를 의미한다. 첫번째 열은 사전 확률 분포이다 - 데이터를 보기전에 s 에 대해서 믿는 것을 의미한다.

어떻게 분포가 고정됐는지 볼 수 있다. 제안된 s 를 끼워넣었다. 가로 선은 현재 s를 파란색으로, 제안된 s를 빨간색 또는 초록색(거절 또는 승인)으로 표현했다.

두번째 열은 우도이고 얼마나 모델이 데이터를 잘 설명하는가를 평가하기 위해 사용한다. 우도함수는 제안 s에 반응해서 바뀐다는 것을 알 수 있다. 파란색 히스토그람은 우리의 데이터이다. 초록색 또는 빨간색의 두꺼운 선은 현재 제안된mu로 계산한 우도이다. 직관적으로, 우도와 데이터가 많이 겹칠수록, 모델이 데이터를 잘 설명하고 있는 것이다. 같은 색의 점선은 제안된 mu이고 점선 파란선은 현재 mu이다.

3번째 열은 우리의 사후확률이다. 여기에, 정규화된 사후확률을 표시했다. 그러나 위에서 말했듯, 정규화되지 않은 사후확률 값을 얻기 위해 현재와 제안된 s에 각각에 대한 사전확률을 두 s에 대한 우도로 곱하고 서로 나우어서 승인 확률을 구할 수 있다.

4번째 열은 trace (생성된 s의 사후확률 표본들)이다. 승인여부가 상관없이 저장한다. 사후확률 밀도 높은 쪽으로 많이 이동하고, 때로는 낮은 쪽으로도 이동함을 볼 수 있다.

MCMC의 마법은 오랜 기간 샘플링을 하면, 이 방법으로 얻어진 표본은 모델의 사후확률 분포에서 나온 것이라는 것이다. 엄격한 수학적 증명도 있지만, 자세히 살펴보지는 않는다.

이 과정을 엿보기 위해서 많은 표본을 추출한 후 그려보자.

posterior = sampler(data, samples=15000, mu_init=1.)
fig, ax = plt.subplots()
ax.plot(posterior)
_ = ax.set(xlabel='sample', ylabel='mu');

 
이것은 보통 trace라 불린다. 사후확률의 근사를 얻기 위해 (이 것을 하는 이유), 간단히 이 trace의 히스토그램을 취하면된다. 모델을 fit하기 위해서 위에서 추출된 데이터와 모양이 비슷하지만, 그 두 데이터는 완전히 다른 것임을 아는 것이 중요하다. 밑의 그림은 mu에 대한 우리의 믿음을 표시한다. 이 경우, 정규분포의 모양이지만, 다른 모델에 대해서, 우도 또는 사전확률과는 완전히 다른 모양이 될 수 있다.

ax = plt.subplot()

sns.distplot(posterior[500:], ax=ax, label='estimated posterior')
x = np.linspace(-.5, .5, 500)
post = calc_posterior_analytical(data, x, 0, 1)
ax.plot(x, post, 'g', label='analytic posterior')
_ = ax.set(xlabel='mu', ylabel='belief');
ax.legend();

보았듯이, 위의 과정을 따름으로써, 해석학적으로 유도된 것과 같은 분포의 표본을 얻을 수 있다.

Proposal width

위에서 proposal width를 0.5로 설정했다. 그것은 아주 좋은 수치임이 밝혀졌다.

일반적으로, 너무 좁은 너비를 원하지는 않는다 왜냐하면 샘플링이 비효율적으로 탐색에 너무 오랜 시간이 걸리기 때문이다.

posterior_small = sampler(data, samples=5000, mu_init=1., proposal_width=.01)
fig, ax = plt.subplots()
ax.plot(posterior_small);
_ = ax.set(xlabel='sample', ylabel='mu');

 
그러나 , 절대 승인되지 않을 움직임 정도로 커지길 원치 않는다.

posterior_large = sampler(data, samples=5000, mu_init=1., proposal_width=3.)
fig, ax = plt.subplots()
ax.plot(posterior_large); plt.xlabel('sample'); plt.ylabel('mu');
_ = ax.set(xlabel='sample', ylabel='mu');

 
그러나, 수학적 증명이 보장하는 것처럼 목표 사후확률 분포로부터 여전히 샘플링 되고 있음을 알수 있다, 단지 좀 비효율적일 뿐이다.

sns.distplot(posterior_small[1000:], label='Small step size')
sns.distplot(posterior_large[1000:], label='Large step size');
_ = plt.legend();

 
좀 더 많이 샘플링 하면, 진짜 사후확률에 더 가깝게 포인다. 중요한 것은 여기서 명백히 사실이 아닌 각자에 대하서 독립인 샘플을 원한다는 것이다. 그러므로, 샘플러의 효율성을 평가하는 중요 메트릭은 autocorrelation 이다. - 즉 표본 i는 표본 i-1, i-2, 등과 얼마나 연관되어 있는지

from pymc3.stats import autocorr
lags = np.arange(1, 100)
fig, ax = plt.subplots()
ax.plot(lags, [autocorr(posterior_large, l) for l in lags], label='large step size')
ax.plot(lags, [autocorr(posterior_small, l) for l in lags], label='small step size')
ax.plot(lags, [autocorr(posterior, l) for l in lags], label='medium step size')
ax.legend(loc=0)
_ = ax.set(xlabel='lag', ylabel='autocorrelation', ylim=(-.1, 1))

 
명백히, 올바른 step width를 자동으로 정하는 똑똑한 방법을 원한다. 한가지 주로 쓰는 방법은 제안을 50%정도로 기각하도록 proposal width를 계속 조절하는 것이다.

Extending to more complex models

두번째 파라메터 sigma에 대해서도 위와 같은 과정을 동일하게 적용할 수 있음을 쉽게 생각할 수 있다. 이 경우, mu와sigma에 대한 제안을 생성할 것이다. 하지만 알고리즘 로직은 거의 동일하다. 또한, Binomial 같은 매우 다른 분포에 대해서도 데이터를 얻을 수 있고, 사후 확률 분포를 얻기 위해 같은 알고리즘을 쓸 수 있다.

매우 cool 하다. 확률적 프로그래밍의 큰 이점이다: 단지 원하는 모델을 정의하고, MCMC로 추론해라. 예를 들어, 밑의 모델은 PyMC3로 매우 쉽게 작성할 수 있다. 역시 마찬가지로 (자동으로 proposal width를 튜닝하는) Metropolis sampler를 사용했고, 같은 결과를 얻었다.