Think Bayes (3) 추정 1
파이썬을 활용한 베이지안 통계, 앨런 B. 다우니 지음, 권정민 옮김/한빛미디어 |
3.1 주사위 문제
상자 안에 4, 6, 8, 12, 20면체 주사위가 들어있다고 해보자. 상자에서 임의로 주사위 하나를 집어서 던졌더니 6이 나왔다면, 이 경우 각 주사위를 선택했을 확률은 어떻게 될까?
문제를 이렇게 접근하여 보자.
- 가설을 나타내자.
- 데이터를 나타내자.
- 우도 함수를 작성하자.
가설은 다음과 같이 나타내어보자.
1 2 |
suite = Dice([4, 6, 8, 12, 20]) |
그리고 정수 1부터 20까지를 사용하여 데이터를 나타내자.
1 2 3 4 5 6 7 |
class Dice(Suite): def Likelihood(self, data, hypo): if hypo < data: return 0 else: return 1.0 / hypo |
주사위를 굴려 나온 값이 주사위가 가지는 눈보다 크다면 발생할 수 없는 일이므로 우도는 0이다. 그렇지 않다면 각 주사위가 가지는 눈의 수에 대한 확률을 가지도록 하자.
주사위를 굴려 6이 나온 경우를 적용하는 코드는 다음과 같다.
1 2 |
suite.Update(6) |
각 주사위의 사후 확률 분포는 다음과 같다.
1 2 3 4 5 6 |
4 0.0 6 0.3921568627450979 8 0.2941176470588235 12 0.19607843137254896 20 0.11764705882352941 |
만약 몇 번 더 던져서 6, 8, 7, 7, 5, 4가 나왔다면 다음과 같이 적용함
1 2 3 |
for roll in [6, 8, 7, 7, 5, 4]: suite.Update(roll) |
그 후 사후 확률 분포는 다음과 같다.
1 2 3 4 5 6 |
4 0.0 6 0.0 8 0.9432484536722127 12 0.055206128061290875 20 0.001545418266496554 |
8면제 주사위를 굴렸을 확률이 94%로 가장 높다.
코드는 http://thinkbayes.com/dice.py 에서 받을 수 있다.
3.2 기관차 문제
프레드릭 모스텔러(Frederick Mosteller)이 제시한 기관차 문제라는 것이 있다.
각 철도에는 이를 지나가는 기관차에 1부터 N까지의 순서로 번호를 붙인다. 어느 날 60호 기관차를 보았다. 이 때 이 철도에는 몇 개의 기관차가 지나가는지 추측해보자.
관측 결과에 따르면, 이 철도에는 60개 이상의 기관차가 있다는 것을 알 수 있다. 하지만 얼마나 더 많을까? 베이지안 추론법에 따르면 이런 문제는 2단계로 해결한다.
- 데이터를 보기 전에 N에 대해서 알고 있는 것이 무엇인가?
- N에 어떤 값이 주어졌을 때 관측한 데이터(60호 기관차)의 우도는 어떻게 되는가?
첫 번째 질문은 사전 확률을, 두 번째는 우도를 묻고 있다.
사전 확률을 추정하기는 어려운 일이라, 간단한 가정으로부터 시작해보도록 하자. 일단 기관차의 수가 1부터 1000까지일 상황이이 있고 이들은 동일한 확률이라고 가정하자.
1 2 |
hypos = xrange(1, 1001) |
이제 우도를 작성해보도록 하자. 주사위의 경우와 비슷하다.
1 2 3 4 5 6 7 |
class Train(Suite): def Likelihood(self, data, hypo): if hypo < data: return 0 else: return 1.0 / hypo |
60호 기관차가 지나갔으므로 이제 갱신을 하여보자.
1 2 3 |
suite = Train(hypos) suite.Update(60) |
출력하기에는 값이 너무 많긴 하지만 가장 높은 가능성은 N이 60일 경우이다. 썩 좋아보이는 답은 아니지만 확률 상 최고치를 고르자고 한다면 60이 최선일 것이다.
만약 사후 확률 분포의 평균을 계산해보면 어떻게 될까?
1 2 3 4 5 6 |
def Mean(suite): total = 0 for hypo, prob in suite.Items(): total += hypo * prob return total |
사실은 Pmf
에서도 비슷한 메서드를 제공한다.
1 2 3 4 |
print suite.Mean() 333.41989326371095 |
평균은 333 정도로 나왔다. N의 가능성을 1000까지로 생각하였을 때 오차(위험)를 최소화하고자 한다면 333을 선택해야 할 것이다. Mean Squared Error (MSE)와 연결되는 개념이다.
코드는 http://thinkbayes.com/train.py 에서 받을 수 있다.
3.3 사전 확률로 할 수 있는 것
기관차 문제를 해결하기 위해서는 임의적인 가설을 여러개 만들어 보아야 한다. 우리는 1000개까지만을 고려하였지만 어떤 사람은 그보다 적게, 혹은 많게 생각할 수도 있을 것이다. 이러한 가정에 민감하게 반응하는지 살펴보아야 한다.
1부터 1000까지의 범위라고 하였을 때는 그 사후 확률의 평균은 333이었지만, 상한이 500이라면 평균은 207이고 상한이 2000이라면 그 평균은 552일 것이다. 따라서 우리가 사용한 추론은 썩 좋지 않은 상황이다.
이러한 경우 우리는 두 가지 방법으로 더 진행할 수 있다.
- 데이터를 더 확보할 것
- 배경 지식을 더 확보할 것
데이터가 더 많다면 사후 확률의 평균 또한 수렴하는 경향을 보일 것이다.
예를 들어 60호 기관차 이후에 30, 90번 기관차도 보았다고 해보자.
1 2 3 |
for data in [60, 30, 90]: suite.Update(data) |
그 후에 사후 확률의 평균은 다음과 같이 변화한다.
상한 | 사후 평균 |
---|---|
500 | 152 |
1000 | 164 |
2000 | 171 |
3.4 사전 확률의 대안
배경지식을 더 많이 수집하는 방법도 있다. 1000대의 기관차를 운영한다는 가정은 별로 합리적이지 않다거나, 철도 회사의 목록을 찾을 수 있을지도 모른다. 기관차 운송 분야의 전문가와 인터뷰를 할 수도 있다. 혹은 로버트 액스텔(Robert Axtell)의 멱법칙에 따라 회사 규모에 따르는 분포를 생각해볼 수도 있다.
멱법칙은 주어진 규모의 회사 크기는 규모의 분포의 역수라는 뜻이다.
$$PMF(x) \propto (\frac{1}{x})^a$$
a는 보통 1에 가깝게 설정되는 매개변수이다.
멱법칙을 적용하여 사전 확률을 다음과 같이 만들 수 있다.
1 2 3 4 5 6 7 |
class Train(Dice): def __init__(self, hpos, alpha=1.0): Pmf.__init__(self) for hypo in hypos: self.Set(hypo, hypo**(-alpha)) self.Normalize() |
동일한 방식으로 사전 확률을 계산해보면, 사후 확률은 그 선택 값에 대하여 민감도가 낮아지게 된다. 또한 N이 700보다 큰 경우는 거의 대부분이 제거될 수 있다. 앞의 30, 60, 90의 관측값을 사용한 경우 사후 확률의 평균 값 변화는 다음과 같이 차이가 거의 없어지는 것을 볼 수 있다.
상한 | 사후 평균 |
---|---|
500 | 131 |
1000 | 133 |
2000 | 134 |
코드는 http://thinkbayes.com/train3.py 에서 받을 수 있다.
3.5 신뢰구간
지금까지는 결과를 단일 점으로 추정하였다. 이러한 점 추정에는 보통 평균, 중간값, 최대 우도값을 사용한다. 하지만 구간을 사용할 수도 있다. 구간은 보통 미 확인 값이 어느 구간내에 올 확률이 90%(혹은 다른 확률)인 구간의 상한값과 하한값을 의미한다. 이것을 신뢰 구간이라고 한다. 신뢰 구간을 계산하는 간단한 방법은 사후 분포 확률을 더한 후 그 값이 5%, 95%에 해당하는 값을 기록하는 것이다. 즉 5분위, 95분위값이다.
thinkbayes
에서 분위값을 계산하는 함수를 제공한다.
1 2 3 4 5 6 7 8 |
def Percentile(pmf, percentage): p = percentage / 100.0 total = 0 for val, prob in pmf.Items(): total += prob if(total >= p): return val |
다음은 이를 사용하는 예제이다.
1 2 3 |
interval = Percentile(suite, 5), Percentile(suite, 95) print interval |
앞의 멱법칙과 3개 관측 값을 사용한 예제의 경우 90% 신뢰구간은 (91, 243)이다. 이는 여전히 얼마나 많은 기관차가 있는지 확신할 수 없다는 것을 보여준다.
3.6 누적 분포 함수
앞의 코드는 분위수를 계산하는데 반복적으로 값과 확률을 더하여 계산하였다. 만약 여러번 반복해서 계산해야 한다면 누적 분포 함수나 Cdf를 사용하는 것이 훨씬 효율적이다.
Cdf는 분포에 대한 정보를 가진다는 면에서 Pmf와 동일하고, 서로 변환하는 것이 가능하다. 하지만 Cdf는 분위수를 조금 더 효율적으로 계산할 수 있다.
thinkbayes
에서는 누적 분포 함수를 나타내는 Cdf
를 제공한다. Cdf를 만들려면 다음과 같이 하자.
1 2 |
cdf = suite.MakeCdf() |
또한 Cdf
는 Percentile
이라는 함수를 제공한다.
1 2 |
interval = cdf.Percentile(5), cdf.Percentile(95) |
3.7 독일 탱크 문제
2차 세계 대전에서 영국은 독일 탱크의 장비 생산량을 추정하는데 통계 분석 기법을 사용햇다. 각 탱크의 섀시와 엔진의 시리얼 번호를 포함한 장부 및 수리 기록을 확보한 것이다. 기록을 분석해보니 시리얼 번호는 제조자에 따라 할당되어 있고, 탱크의 유형은 100개의 숫자 블록으로 이루어져 있는데 각 블록 숫자는 순서대로 증가하긴 하지만 각 블록의 모든 자리가 사용된 것은 아니다. 따라서 독일 탱크 생산량 추정은 100개의 숫자를 포함하는 블록 범위 내에서 기관차 문제로 생각할 수 있다. 결과는 도출한 내용이 대부분 맞았다고 한다.
3.8 토의
베이지안(Bayesian, 베이즈 주의자)가 선행 분포를 택하는 방법은 두 가지가 있는데, 하나는 문제의 배경 지식을 가장 잘 표현은 선행 분포를 고르는 방법이다. 이는 선행 분포가 정보적이라고 한다. 이 경우에는 사람들이 서로 다른 배경 지식을 사용하거나 해석할 수 있다는 점이 문제로 지적된다. 즉 주관적인 것이다.
다른 하나는 비정보적 선행 분포로 관측된 데이터만을 신뢰하도록 다른 제약을 최소한으로 줄이는 것이다. 때로는 추정량에 대한 최소한의 선행 정보를 나타내도록 고유한 선행 분포를 정의할 수도 있다.
비정보적 선행 분포는 조금 더 객관적이긴 하지만, 저자는 선행 분포를 선호한다. 첫번째 이유로, 베이지안 분석은 모델 판단에 기반하여 이루어지느데, 선행 분포를 선택하지 않더라도 전체 분석 자체가 거의 주관적으로 이루어지기 때문이다. 또한 데이터가 아주 많거나 아니면 아주 적거나 극단적으로 속하게 되는 경향이 있는데, 데이터가 많은 경우 문제가 되지 않지만 데이터가 적은 경우 배경 지식을 사용하는 것은 매우 큰 차이를 나타낸다. 따라서 객관성을 유지하는 것보다는 알고 있는 모든 지식을 쏟아넣어야 할 것이다.