베이지안 PYMC3

베이지안 라이브러리가 필요한 이유 

 

베이지안이 어려운 이유는 실제로 conjugate distribution 을 찾기 어렵기 때문이다. 

 

이로 인해 MCMC 샘플링을 사용하는 방식을 많이 사용하고 대표적인 library 로 PYMC3 가 있다. 복잡한 과정을 모두 대신해준다. 

생각보다 듀토리얼이 없어서 아주 간단한 이론적 배경과 사용법을 정리해본다.

 

해당 글을 읽었을 때의 유익 

저자는 본 글을 읽고나서 베이지안 사고방식이 어떤 의미인지, 그리고 해당 사고방식을 어떻게 접목해야하는지 깨달을 수 있는 기회가 되었다. 

 

해당 내용은 https://nbviewer.org/github/markdregan/Bayesian-Modelling-in-Python/blob/master/Section%201.%20Estimating%20model%20parameters.ipynb 

Jupyter Notebook Viewer

을 대부분 참고하였다. 본 글이 이해가 잘 가지 않는다면 위 글을 읽는 것을 추천한다. 

 

베이지안 통계학자들은 어떤 사고방식을 가지고 있을까?

아래의 시나리오가 있다고 생각해보자. 

 

한 소년이 집 앞을 지나가는 자동차를 매일 관찰하였다. 그는 매일 스쳐간 자동차의 수를 노트에 기록하였고 한주가 지나자 소년의 노트에는 2, 33, 20, 29, 20, 30, 18 의 수가 기록되었다. 

 

데이터가 관측이 완료되었다. 즉, 이로 인해 생성되는 parameter는 고정되었고 그 값이 바뀌지 않는다. 

하지만, 베이지안 통계학자들은 확률 분포를 사용하여 데이터를 통해 고정된 파라미터의 uncertainty 를 나타내고자 한다. 

 

소년의 행위는 포아산 분포를 통해 데이터를 모델링 할 수 있는 흔한 상황이다. (포아산 분포에 대해서는 따로 언급하지 않는다.)

포아산 분포는 $\mu$ 라는 parameter 하나로 데이터의 평균과 분산 모두를 설명한다. 

$$p(x \ | \ \mu) = \frac{e^{-\mu}\mu^{x}} {x!} \mbox{    for     } x = 0, 1, 2, \cdots$$

 

아래는 각각 다른 $\mu$ 를 가진 세 가지의 포아산 분포 예시이다. 

fig = plt.figure(figsize=(11,3))
ax = fig.add_subplot(111)
x_lim = 60
mu = [5, 20, 40]
for i in np.arange(x_lim):
    plt.bar(i, stats.poisson.pmf(mu[0], i), color=colors[3])
    plt.bar(i, stats.poisson.pmf(mu[1], i), color=colors[4])
    plt.bar(i, stats.poisson.pmf(mu[2], i), color=colors[5])
    
_ = ax.set_xlim(0, x_lim)
_ = ax.set_ylim(0, 0.2)
_ = ax.set_ylabel('Probability mass')
_ = ax.set_title('Poisson distribution')
_ = plt.legend(['$\mu$ = %s' % mu[0], '$\mu$ = %s' % mu[1], '$\mu$ = %s' % mu[2]])

각 $\mu$ 값에 의해 결정된 평균과 분산을 확인할 수 있다.

 

이제 실전으로 들어가보자. 

 

주어진 데이터: 메세지를 응답하는데 걸리는 시간과 걸린 시간에 따른 답장 메세지 개수

우리는 이 데이터로 포아산 분포를 모델링 하고, parameter $\mu$를 찾을 수 있다. 해당 parameter를 추정할 수 있는 두가지 방식이 있다: 빈도주의 방식, 베이지안 방식. 두 방식으로 parameter를 추정해보자.

fig = plt.figure(figsize=(11,3))
_ = plt.title('Frequency of messages by response time')
_ = plt.xlabel('Response time (seconds)')
_ = plt.ylabel('Number of messages')
_ = plt.hist(messages['time_delay_seconds'].values, 
             range=[0, 60], bins=60, histtype='stepfilled')

 

$\mu$ 를 추정하는 빈도주의적 방법 

빈도주의자는 포아산 분포의 $\mu$ 를 어떻게 추정할까? 

 

여기서는 함수의 likelihood를 최대로 하기 위하여 최적화 테크닉을 사용한다.

아래의 poisson_logprob() 함수는 포아산 모델과 parameter 값이 주어질 때 관측된 데이터의 총 likelihood를 return 한다. 

 

opt.minimize_scalar 함수는 주어진 관측 데이터로부터 가장 신뢰할 만한 $\mu$ (log likelihood 를 최대화하는) 를 찾는다. 

해당 최적화 테크닉은 가능한 $\mu$를 interate하며 가장 높은 likelihood를 갖는 최적의 $\mu$ 를 찾아낸다. 

 

y_obs = messages['time_delay_seconds'].values

def poisson_logprob(mu, sign=-1):
    return np.sum(sign*stats.poisson.logpmf(y_obs, mu=mu))

freq_results = opt.minimize_scalar(poisson_logprob)
%time print("The estimated value of mu is: %s" % freq_results['x'])

출력

The estimated value of mu is: 18.2307692324
CPU times: user 90 µs, sys: 45 µs, total: 135 µs
Wall time: 101 µs

 

측정된 $\mu$ 는 18.0413533867 이다. 이 최적화 테크닉은 매우 효율적으로 특정 값을 찾아 제공한다. 다만 uncertainty 에 대한 지표는 제공하지 못한다. 

 

(stats.poisson.logpmf 는 포아산 질량 함수에 log 를 취한 값을 return 한다. 각 값을 joint probability 한 결과를 얻는 것과 동일한데, log 를 취했다보니 곱대신 합이다. 바로 아래의 포아산 질량 함수를 참고해보자.)

$$p(x \ | \ \mu) = \frac{e^{-\mu}\mu^{x}} {x!} \mbox{    for     } x = 0, 1, 2, \cdots$$

 

아래 차트는 우리가 최적화 하는 함수를 시각화한 결과이다. 주어진 데이터와 모델에서 $\mu$ 의 값의 변화에 따른 log probability 를 확인할 수 있다. 

optimizer 는 마치 언덕을 올라가듯 작동한다 - 커브의 랜덤한 점에서 시작하여 더 올라갈 수 없는 지점까지 점진적으로 등산한다. 

 

x = np.linspace(1, 60)
y_min = np.min([poisson_logprob(i, sign=1) for i in x])
y_max = np.max([poisson_logprob(i, sign=1) for i in x])
fig = plt.figure(figsize=(6,4))
_ = plt.plot(x, [poisson_logprob(i, sign=1) for i in x])
_ = plt.fill_between(x, [poisson_logprob(i, sign=1) for i in x], 
                     y_min, color=colors[0], alpha=0.3)
_ = plt.title('Optimization of $\mu$')
_ = plt.xlabel('$\mu$')
_ = plt.ylabel('Log probability of $\mu$ given data')
_ = plt.vlines(freq_results['x'], y_max, y_min, colors='red', linestyles='dashed')
_ = plt.scatter(freq_results['x'], y_max, s=110, c='red', zorder=3)
_ = plt.ylim(ymin=y_min, ymax=0)
_ = plt.xlim(xmin=1, xmax=60)

위 과정은 parameter $\mu$ 를 18로 추정하였다. 포아산 분포의 parameter $\mu$ 는 평균과 분산을 동시에 나타낸다. 아래 차트는 해당 분포를 시각화한 결과이다. 

 

fig = plt.figure(figsize=(11,3))
ax = fig.add_subplot(111)
x_lim = 60
mu = np.int(freq_results['x'])
for i in np.arange(x_lim):
    plt.bar(i, stats.poisson.pmf(mu, i), color=colors[3])
    
_ = ax.set_xlim(0, x_lim)
_ = ax.set_ylim(0, 0.1)
_ = ax.set_xlabel('Response time in seconds')
_ = ax.set_ylabel('Probability mass')
_ = ax.set_title('Estimated Poisson distribution for Hangout chat response time')
_ = plt.legend(['$\lambda$ = %s' % mu])

포아산 모델과 측정된 $\mu$ 는 10보다 작거나 30보다 큰 관측치가 낮은 확률로 존재할 수도 있다는 것을 보여준다. 대부분의 확률 밀도는 10과 30 사이에 위치한다. 하지만, 우리는 이 결과가 우리가 관측한 결과와 다르다는 것을 알 수 있다. 

$\mu$ 를 추정하는 베이지안 방법 

$$\overbrace{p(\mu \ |\ Data)}^{\text{posterior}} = \dfrac{\overbrace{p(Data \ | \ \mu)}^{\text{likelihood}} \cdot \overbrace{p(\mu)}^{\text{prior}}}{\underbrace{p(Data)}_{\text{marginal likelihood}}}$$a

 

간단하게 위 식에 대해 이해해보자. 

 

(marginal likelihood 는 발생할 수 있는 모든 경우의 확률을 의미한다. 즉 발생할 수 있는 모든 경우를 분포로 그려본다면 이 때 해당하는 너비이다. )

 

위 그림은 다음과 같이 해석될 수 있다. (위 그림의 아래부터 대응된다.)

 

• 관측을 통해 데이터를 얻는다. ($y_{i}$)
• 관측된 데이터는 알 수 없는 어떠한 과정을 통해 생성 되었다. 하지만 우리는 이 과정이 포아산 분포를 통해 나타낼 수 있다고 믿는다. (Likelihood)
• 포아산 분포는 parameter를 한개 갖는다. ($\mu$ ) $\mu$는 0 ~ 60 사이 값이다. (Prior)
• 우리는 Prior $\mu$ 를 uniform 분포로 모델링한다. 이는 $\mu$ 가 해당 범위내 어디쯤에 어느정도 위치할지 예상할 수 없기 때문이다. 

놀라운 MCMC 의 메커니즘 

The process of Markov Chain Monte Carlo (MCMC) is nicely illustrated in the below animation. The MCMC sampler draws parameter values from the prior distribution and computes the likelihood that the observed data came from a distribution with these parameter values.

Markov Chain Monte Carlo (MCMC)의 과정은 아래 에니메이션 시각화를 통해 살펴볼 수 있다. MCMC sampler는 prior 분포로부터 parameter 값들을 뽑고 이 파라미터 값들을 갖는 분포로부터 관측한 데이터가 나올 수 있는 likelihood를 계산한다. 

 

$$\overbrace{p(\mu \ |\ Data)}^{posterior} \varpropto \overbrace{p(Data \ | \ \mu)}^{likelihood} \cdot \overbrace{p(\mu)}^{prior}$$

 

위 공식은 MCMC sampler 을 이끄는 한줄기 빛 역할을 한다. MCMC sampler 가 prior 에서 parameters 를 뽑아낸 이후,  주어진 데이터에 대하여 prior에서 뽑아낸 parameters의 likelihood를 계산한다. (그리고 sampler를 더 높은 확률의 영역으로 향하도록 이끈다.)

 

사전에 언급했던 빈도주의적 최적화 테크닉과 개념적으로 유사한 면이 있다. MCMC sampler 도 더 높은 확률의 likelihood 영역을 향해 해맨다. 그러나 베이지안 방법은 절대적 최대치 값을 찾는 것보다 확률이 가장 높은 영역의 sample들을 탐색하고 모으는 것에 관심이 있다. 이렇게 하여 얻은 모든 sample들은 믿을 만한 parameter로 고려된다. 

 

with pm.Model() as model:
    mu = pm.Uniform('mu', lower=0, upper=60)
    likelihood = pm.Poisson('likelihood', mu=mu, observed=messages['time_delay_seconds'].values)
    
    start = pm.find_MAP()
    step = pm.Metropolis()
    trace = pm.sample(200000, step, start=start, progressbar=True)

Applied interval-transform to mu and added transformed mu_interval_ to model.
  2%|▏         | 3769/200000 [01:13<1:05:21, 50.04it/s]

 

( 정리하자면, 만약 사전분포의 분산이 커서 precision 이 낮아지게 되면 prior 에서 뽑는 parameter 의 value 의 범위도 넓어지게 되고, 이러면 관측 데이터의 likelihood를 찾을 수 있는 범위가 넓어져 관측 데이터에 영향을 더욱 많이 받게 되는 것으로 보인다. 반면 분산이 작아 precision 이 크게되면 prior 에서 뽑는 parameter 의 value 의 범위도 좁아지고 관측 데이터의 likelihood를 찾을 수 있는 범위도 좁아져 관측 데이터에 영향을 거의 받지 않고 사전 정보에 남아있게 되는 것으로 보인다. )

 

상단의 코드는 높은 likelihood의 $\mu$ 사후 분포를 탐색하여 200,000개의 믿을만한 $\mu$ sample 을 수집했다. 하단의 왼쪽 차트는 $\mu$ 를 위해 수집된 분포를 보여준다. 해당 분포의 평균은 빈도주의 방식의 예측값과 거의 동일하다. 그러나, 이 뿐 만이 아니라 uncertainty 측정값 또한 얻을 수 있고, 이를 통해 $\mu$는 17 ~ 19 사이에 충분히 믿을 만한 값들이 존재함을 알 수 있다. 추후 다루겠지만 uncertainty 값은 정말 중요하다. 

 

초기 sample 버리기 

상단 MCMC code 의 pm.find_MAP() 에 대하여 간략하게 설명해본다. MAP 은 Maximum a posteriori estimation 을 뜻한다. MCMC sampler 가 sampling을 시작하기 좋은 위치를 찾을 수 있도록 돕는다. 이상적으로는 likelihood 가 높은 영역에서 모델을 시작하지만, 때때로는 이런 일은 일어나지 않는다. 이에 따라, 초기에 수집된 samples 은 종종 버려진다. (이를 burnin samples 라 한다.)

_ = pm.traceplot(trace, varnames=['mu'], lines={'mu': freq_results['x']})

fig = plt.figure(figsize=(11,3))
plt.subplot(121)
_ = plt.title('Burnin trace')
_ = plt.ylim(ymin=16.5, ymax=19.5)
_ = plt.plot(trace.get_values('mu')[:1000])
fig = plt.subplot(122)
_ = plt.title('Full trace')
_ = plt.ylim(ymin=16.5, ymax=19.5)
_ = plt.plot(trace.get_values('mu'))

 

Model convergence

상단 모델이 $\mu$의 값을 추정했다고 해서 모델이 꼭 주어진 데이터에 대해서 좋은 결과를 찾았다는 뜻은 아니다. 몇가지 추천할만한 확인 과정이 있다. 첫번째로, output 의 trace 를 확인해본다. trace가 확확 뛰어야 하며 마치 송충이 같아 보여야한다. 만약 trace 가 위아래로 움직이는 뱀 흔적처럼 보이거나 한 곳에 갇힌것처럼 보이면 수렴 문제가 있거나, MCMC sampler 의 결과물이 믿을 만 하지 않은 것일 수 있다. 

 

Autocorrelation plot

두번째 확인 해볼 수 있는 방법은 autocorrelation test 이다. 이는 연속적인 MCMC sampling chain 의 sample들 간의 correlation을 측정하는 방법이다. sample들이 서로 낮은 correlation을 가지고 있으면 높은 correlation을 일 때보다 parameter 측정 값에 새로운 정보를 더하는 모양이다. 

시각적으로 보면 상대적으로 빠르게 0으로 감소한 다음 0 상관관계 위아래로 진동하는 autocorrelation 차트가 우리가 원하는 결과이다. 만약 autocorrelation 차트가 떨어지지 않으면 주로 혼합 불량 신호이므로 모델을 다시 선정하고 (ex. likelihood) 샘플링 방식 (ex. matropolis)도 다시 고려해야 한다. 

 

_ = pm.autocorrplot(trace[:2000], varnames=['mu'])