데이터과학 유망주의 매일 글쓰기 — 쉬어가는 한주 — 5
베이지안과 PyMC3
# 빈도주의, #베이지안, #PyMC3
오늘 한일:
쉬어가는 한주의 마지막 날인 오늘은 첫 번째 섹션에서 추가적인 과제였지만 내가 시간을 내서 미쳐 학습하지 못한 PyMC3에 대해 다루기로 했다. 이 도구는 반복적이면서도 소모적인 확률계산을 자동화 할 수 있는 유용한 도구로써, 특히 베이지안 정리를 활용하는데 상당히 유용한 도구다.
이 도구를 제대로 활용하기 위해서는 베이지안 정리에 대해 다시 한번 짚고 넘어가는 것이 필요하다고 판단되었다.
PyMC3
빈도주의와 베이지안의 차이
빈도주의(frequentist)의 관점에서 특정한 이벤트가 발생할지 안 할지를 무한대로 시험하는 것은 불가능하다. 예를 들어, 내일 비가올 확률이 80%라고 한다면, 이 것이 어떠한 의미인지 설명하기가 쉽지 않다.
바로 여기서 베이지안 정리(Bayesian Theorem)가 좋은 해결책이 될 수 있다. 베이지안 정리는 이벤트가 신빙성 있는지에 대한 측정을 하게 해 줌으로써, 이벤트가 얼마나 높은 확률로 나타날 지를 예측할 수 있게 해주기 때문이다.
즉, 베이지안 정리는 “그 어떤 것도” 절대적으로 확실하지 않다는 것을 전제로 한다. 그러나, 이벤트가 얼마나 높은 확률로 발생할 지는 말해 줄 수 있다. 데이터가 더 많을 수록, 확신에 대한 신빙성은 더 커질 수 있다.
데이터 사이언티스트로써, 모든 것을 의심해야하는 과학자의 시선에서 본다면, 베이지안 정리가 훨씬 더 직관적이라 할 수 있다.
그러나, 베이지안 정리를 사용한다는 것은 상당히 계산적으로도, 개념적으로도 복잡할 수 있다. 상당히 많은 경우 수학적으로 매우 어렵고 긴 계산이 필요할 수 있기 때문이다. 수학자로써도 이러한 계산은 상당히 반복적이고 소모적이라고 느껴진다.
PyMC3란?
이러한 소모적인 계산을 편리하게 할 수 있게 해주는 도구가 PyMC3라 할 수 있다. 그렇다면 이 라이브러리는 도대체 어떤 도구이며, 어떤 경우에 활용할 수 있을까?
예제 — 동전
동전 하나를 3번 던졌다고 가정하자. 뒷면(tail)이 0, 앞면(head)이 1이라고 가정한다. 그리고 결과는 [0,1,1]이 나왔다고 하자. 이렇게 나왔다면, 과연 동전이 공정하다고 할 수 있을까? 동전의 앞면이 나올 확률을 “t”라고 가정한다면, “t=12”라는 것이 사실이라고 증명할 수 있을까?
사실 위 한 번의 실험 결과만 가지고는, 동전의 그 어떤 정보도 알 수 없다. 그러므로, 아직은 어떤 것도 확실하다고 할 수 없다. 빈도주의의 관점에서 t는 아래와 같이 계산 될 수 있다.
이러한 빈도주의 방식의 계산은, 각 앞면과 뒷면이 어느 확률로 나올 지에 대해서는 설명하지 못한다는 단점이 있다. 바로 여기서 베이지안 정리가 큰 도움이 될 수 있다. 베이지안 관점에서는 아직 t에 대한 정보를 모르기 때문에 0과 1 사이의 그 어떤 값도 될 수 있다고 가정한다. 즉, t는 0과 1 사이의 균등분포(uniform distribution)를 갖게 된다. 이 균등분포의 PDF(Probability Density Function)을 계산한다면 아래와 같다.
여기서 부터는 증거와 관측값을 통해 t의 분포를 업데이트 하면된다.
D를 증거(이 예에서는 동전을 던진 행위에 대한 결과를 의미한다.)라고 할 때, 베이지안 정리에 의거하여 사후확률(posterior probability)는 아래와 같이 계산될 수 있다.
p(D/t)는 확률 함수를 뜻하며, p(t)는 사전 확률(prior probability, 이 예에서는 0과 1 사이의 균등분포를 말한다. 여기서는 2가지의 방법을 사용할 수 있다.
방법 1. 손으로 계산하기
t가 주어졌을 때, 3번을 던져 2번 앞면이 나올 확률은 아래와 같이 계산할 수 있다.
p(t)를 1(앞면으로 나올 확률이 100%)이라고 가정한다면, 아래와 같이 분모를 표현할 수 있다
위 적분을 계산한다면 0.25라는 값을 얻을 수 있다. 그러므로 p(t/D)값은 아래와 같다.
만약 t의 사전확률이 beta 분포일 때(alpha와 beta 값을 가질 때), p(t) = B(alpha, beta)로 정의된다. 샘플 사이즈는 N이며, 이 중 k가 앞면이다. t의 사후확률은 B(alpha+k, beta+N-k)이다. 이 예제에서는, alpha=beta=1, N=3, k=2로 정의된다.
방법 2. PyMC3 사용
위 방법은 conjugate prior라는 원칙을 사용하는데, 현실을 충분히 반영하지 못할 수 있다는 단점이 있다. 현실에서는 이러한 conjugate prior를 찾기 어려운 경우가 있기 때문이다.
이 것은 Markov Chain Monte Carlo(MCMC)라는 방법을 사용하여 극복할 수 있는데, 사후확률을 예측하는 방법이다. 이 방법은 PyMC3를 이용하여 구현할 수 있다.
일단 설명을 위해 아래와 같은 코드가 필요하다.
import pymc3 as pmimport scipy.stats as statsimport pandas as pdimport matplotlib.pyplot as pltimport numpy as np%matplotlib inlinefrom IPython.core.pylabtools import figsize
다음으로 t의 사전확률을 초기화해야한다. PyMC3를 이용해 아래의 코드로 구현할 수 있다.
with pm.Model() as model:theta=pm.Uniform('theta', lower=0, upper=1
이후 관측된 데이터를 이용하여 모델을 학습시킨다. 아래의 코드를 활용한다.
occurrences=np.array([1,1,0]) #our observationwith model:obs=pm.Bernoulli("obs", theta, observed=occurrences) #input the observationsstep=pm.Metropolis()trace=pm.sample(18000, step=step)burned_trace=trace[1000:]
PyMC3는 내부적으로 Metropolis-Hastings 알고리즘이란 것을 사용하여 사후확률을 예측한다. Trace 함수는 사후확률에서 추출한 샘플의 수를 나타낸다. 마지막으로, 알고리즘은 초기에 불안정할 수 있으므로, 몇번 반복을 한 이후 샘플을 추출하는 것이 좋다. 시각화를 위해 사후확률에서 얻은 샘플을 histogram으로 그려내어, TDF(True Density Function)과 비교한다.
from IPython.core.pylabtools import figsizep_true=0.5figsize(12.5, 4)plt.title(r"Posterior distribution of $\theta$")plt.vlines(p_true,0, 2, linestyle='--', label=r"true $\theta$ (unknown)", color='red')plt.hist(burned_trace["theta"], bins=25, histtype='stepfilled', density=True, color='#348ABD')x=np.arange(0,1.04,0.04)plt.plot(x, 12*x*x*(1-x), color='black')plt.legend()plt.show()
그래프에서 볼 수 있는 것처럼, PyMC3를 통해 계산한 방법이 사후 확률과 매우 비슷한 것을 알 수 있다.
그렇다면 샘플 사이즈를 높이면 어떨까?
이전에 언급한 것처럼, 데이터가 더 많을수록, t의 확률에 대한 확신이 더 높아질 수 있다. 이번에는 동전을 1000번 던졌다고 가정하자. PyMC3를 이용해 t의 사후확률을 구하고, 이 확률에서 얻은 샘플을 histogram으로 그려본다.
N=1000 #the number of samplesoccurences=np.random.binomial(1, p=0.5, size=N)k=occurences.sum() #the number of head#fit the observed datawith pm.Model() as model1:theta=pm.Uniform('theta', lower=0, upper=1)with model1:obs=pm.Bernoulli("obs", theta, observed=occurrences)step=pm.Metropolis()trace=pm.sample(18000, step=step)burned_trace1=trace[1000:]#plot the posterior distribution of theta.p_true=0.5figsize(12.5, 4)plt.title(r"Posterior distribution of $\theta for sample size N=1000$")plt.vlines(p_true,0, 25, linestyle='--', label="true $\theta$ (unknown)", color='red')plt.hist(burned_trace1["theta"], bins=25, histtype='stepfilled', density=True, color='#348ABD')plt.legend()plt.show()
사후확률이 t값을 중심으로 분포된 모습을 볼 수 있다. 이 상태에서 샘플의 평균을 수하면 t값을 예측할 수 있다.
burned_trace1[‘theta’].mean()위 코드를 실행시키면 약 0.5정도의 값을 얻을 수 있는데, t값과 상당히 비슷한 것을 알 수 있다.
보이는 것처럼 PyMC3는 상당히 통계적인 계산을 잘 하는 모습을 볼 수 있다. 이 도구를 이용하면 소모적인 확률 계산 프로그래밍의 수고를 덜 수 있다.
앞으로 할일:
오늘은 확률 계산을 자동화하여 수학자들과 데이터 사이언티스트에게 불편한 수고를 덜어주는 PyMC3에 대해 배워보았다. 수십년 전만해도 이렇게 복잡한 계산을 손으로 다 했을텐데, 컴퓨터의 발달로 이제 더 이상 그러한 수고를 하지 않아도 된다는 것이 참 감사하다.
결국, 갈수록 개념과 원리 자체를 이용하는 것이 더 중요하다는 생각이든다. 프로그램코드로 복잡한 계산은 얼마든지 할 수 있는 만큼, 인간은 창의력을 이용하여 컴퓨터에게 잘 설계된 알고리즘에 따라 계산을 하게 하는 것이 중요하다고 느낀다. 그런 점에서, 베이지안 같은 수학적 개념을 이해하는 것이 정말 중요함을 다시 느낀다.
다음주 부터 시작되는 딥러닝 과정에서도, 물론 힘들겠지만 수학적인 원리를 조금이라도 더 잘 이해할 수 있게되기를 바라며, 또 한번의 스프린트를 달려보자고 다짐한다.
이 글은 대부분 Introduction to PyMC3: A Python package for probabilistic programming 이라는 블로그를 참조하였습니다.
참조:
