for textmining

베이지안 추론(1)

|

이번 글에서는 베이지안 추론에 대해 살펴보도록 하겠습니다. 이 글은 ‘밑바탁부터 시작하는 데이터과학(조엘 그루스 지음, 인사이트 펴냄)’과 ‘Think Bayes(앨런 B. 다우니 지음, 권정민 옮김, 한빛미디어 펴냄)’을 정리했음을 먼저 밝힙니다. 그럼 시작하겠습니다.

문제 정의와 베타분포

동전이 하나 있습니다. 우리는 이 동전이 공평한 동전인지 아닌지 알아보고 싶습니다. 이 동전을 던졌을 때 앞면이 나올 확률이 $p$라고 한다면 이 $p$가 어떤 값인지 추정해보고 싶은 겁니다.

동전던지기 실험은 이항분포를 따릅니다. 이항분포란 성공확률이 $p$이고, 그 결과가 성공 혹은 실패뿐인 실험을 $n$번 반복시행할 때 성공횟수의 분포를 가리킵니다.

이항분포 파라메터 $p$의 사전확률베타분포를 따릅니다. 여기에서 사전확률이란 데이터를 관측하기 전 가설로 세운 확률을 뜻합니다.

이 때문에 우리가 알고 싶은 앞면이 나올 확률 $p$는 베타분포를 사전분포로 사용하게 됩니다. 베타분포도 정규분포처럼 중심이 높은 포물선 형태를 갖는데요. 베타분포의 파라메터는 $α$와 $β$이며 그 중심은 다음과 같습니다.

$α$ / ($α + β$)

$p$의 사전확률과 관련해 $α$는 성공(앞면), $β$는 실패(뒷면)와 연관이 있고 $α, β$의 크기는 믿음의 크기와 관련이 있습니다. 만약 동전에 대해 어떤 선입견도 취하고 싶지 않다면 $α, β$를 모두 1로 정하면 됩니다. 앞면이 55%의 경우로 나타난다고 ‘굳게’ 믿고 있다면 $α$를 55로, $β$를 45로 가정할 수도 있습니다.

이항분포 파라메터 $p$의 사후확률 역시 베타분포를 따릅니다. 여기에서 사후확률이란 데이터를 확인한 이후의 가설 확률을 가리킵니다. 바로 우리가 알고 싶은 값이죠. 사후확률은 사전확률의 업데이트 버전 정도로 이해하면 좋을 것 같습니다.

추론 과정

동전을 여러번 던져서 앞면이 $h$번, 뒷면이 $t$번 나왔다고 칩시다. 베이즈 규칙을 활용해 유도하면 $p$의 사후확률 분포는 파라메터가 $α+h, β+t$인 베타분포를 따른다고 합니다.

이제 동전을 10번 던져서 앞면을 3번 관측했다고 가정해 보겠습니다. 실험 대상 동전이 공평한 동전이라고 생각해 $α, β$를 모두 1로 두었다고 하면, 앞면이 나올 확률 $p$의 사후분포는 파라메터가 1+3, 1+7인 베타분포가 될 겁니다.

공평한 동전이라는 생각이 다소 강하게 들어 $α, β$를 각각 20으로 정했었다면 $p$의 사후분포는 Beta(20+3, 20+7)이 됩니다. 즉, 동전이 뒷면으로 살짝 편향되었다는 것으로 믿음이 바뀔 것을 의미합니다.

앞면이 잘 나오는 동전이라는 확신 때문에 $p$의 사전분포를 Beta(30, 10)로 두었다면 사후분포는 Beta(30+3, 10+7)이 됩니다. 이 경우 동전이 아직도 앞면이 편향되어 있다고 믿지만, 처음 생각했던 것보다는 덜 강하게 믿는다는 걸 뜻합니다.

데이터 관측 후 $α, β$ 변화에 따른 $p$의 사후분포 변화는 다음 그림과 같습니다. 예컨대 Beta(4,8)의 경우 중심이 0.33인데요, $p$의 사후확률이 0.33일 가능성이 제일 높다는 뜻으로 받아들일 수 있습니다.

동전을 많이 던져볼 수록 관측 전 가정한 사전분포는 점점 의미가 없어집니다. 예를 들어 편향된 동전이라는 사전 믿음이 아무리 강했을지라도 동전을 2000번 던져서 앞면에 1000번 나왔다면 생각을 바꿀 수밖에 없기 때문입니다. 데이터가 충분하다면 서로 다른 사전확률을 가지고 시작한다고 해도 동일한 사후확률로 수렴하는 경향이 있습니다.

이러한 과정을 통해 $p$가 가지는 값을 확률적으로 추정할 수 있게 됩니다. 이를 베이지안 추론이라고 합니다. 이는 p-value을 이용한 가설검정 결과와는 접근 방식이 근본적으로 다릅니다. 예컨대 다음과 같습니다.

베이지안 추론 : 관측한 데이터와 사전분포를 고려해볼 때 동전의 앞면이 나올 확률이 49~51%일 경우는 5%밖에 되지 않는다

가설검정 : 동전이 공평하다면 이렇게 편향된 데이터를 관측할 경우는 5%밖에 되지 않는다

베이지안 추론 과정은 사전에 가설을 세운 뒤 실험(데이터 관측) 결과로 가설을 조금씩 업데이트하면서 완성됩니다. 평소 우리가 경험적으로 체감하고 있는 확률들(예컨대 전철에서 자리가 날 확률, 휴강을 할 확률 등등)은 사실 베이지안 추론과 본질상 다르지 않은 것 같습니다.

시각화

위 그림을 만드는 데 사용한 파이썬 코드는 다음과 같습니다.

import math
from matplotlib import pyplot as plt
def B(alpha, beta):
    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta)
def beta_pdf(x, alpha, beta):
    # [0, 1] 구간 밖에서는 밀도가 없음
    if x < 0 or x > 1:
        return 0
    return x ** (alpha - 1) * (1 - x) ** (beta - 1) / B(alpha, beta)

xs = [x / 100.0 for x in range(0,100)]
plt.plot(xs,[beta_pdf(x,alpha=4,beta=8) for x in xs],'-',label='Beta(4,8)')
plt.plot(xs,[beta_pdf(x,alpha=23,beta=27) for x in xs],'--',label='Beta(23,27)')
plt.plot(xs,[beta_pdf(x,alpha=33,beta=17) for x in xs],':',label='Beta(33,17)')
plt.legend()
plt.title('Beta pdfs')
plt.show()

구현

동전을 250번 던졌는데 앞면은 140번, 뒷면은 110번 나왔다고 칩시다. 앞면이 나온 확률 x 100을 $x$라고 할 때 코드는 다음과 같습니다. thinkbayes.py는 이곳에서 내려받을 수 있습니다.

import thinkbayes as tb

class Euro(tb.Suite):
    """Represents hypotheses 
    about the probability of heads."""

    def Likelihood(self, data, hypo):
        """Computes the likelihood of the data 
        under the hypothesis.
        hypo: integer value of x, 
        the probability of heads (0-100)
        data: tuple of (number of heads, number of tails)
        """
        x = hypo / 100.0
        heads, tails = data
        # 이항분포의 우도함수
        like = x**heads * (1-x)**tails
        return like

Euro 클래스는 앞면과 뒷면이 나올 확률이 동일하다는 가정에서 출발해 관측 결과를 업데이트한 값을 반환합니다. 다음과 같이 관측 결과를 업데이트한 뒤 suite 객체에 포함된 MaximumLikelihood 함수를 호출하면 56(=140/250*100)이 반환됩니다.

한편 suite.Mean() 함수는 가설(x)에 사후확률을 곱한 값을 모두 더한 평균값을 반환하는데요. 호출 결과 55.95238095가 나왔습니다.

# 가설 설정(x=1~100) 및 객체 선언
suite = Euro(xrange(0,101))
# 관측치 업데이트
data = (140, 110)
suite.Update(data)
# 최대 우도에 해당하는 x 반환
suite.MaximumLikelihood()
# 평균 반환
suite.Mean()

$x$의 사전분포는 베타분포입니다. 베타분포의 모양은 $α, β$ 에 따라 달라지는데요. 만약 사전분포가 $α, β$에 대한 베타분포이고, 앞면 $h$와 뒷면 $t$의 데이터를 가지고 있다면 사후확률은 $α+h, β+t$에 대한 베타분포가 됩니다. 이를 구현한 코드는 다음과 같습니다.

# 베타분포 선언, alpha=beta=1
beta = tb.Beta()
# 관측치 업데이트
# alpha = 1 + 140
# beta = 1 + 110
beta.Update(data)
# 베타분포의 중심(평균) 반환
print beta.Mean()

beta.Mean()은 베타분포의 중심을 출력해줍니다. 그 결과는 0.5595238로 suite.Mean()의 결과와 동일합니다.

Comments