# -*- coding: utf-8 -*-
# 이번챕터에 사용될 유틸 클래스
from thinkbayes import Suite
class Dice(Suite):
"""Represents hypotheses about which die was rolled."""
# 우도 계산 함수. 데이터로 나온 숫자가 주사위보다 크면 확률=0, 아니면 1/눈의수
def Likelihood(self, data, hypo):
"""Computes the likelihood of the data under the hypothesis.
hypo: integer number of sides on the die
data: integer die roll
"""
if hypo < data:
return 0
else:
return 1.0/hypo
4면체, 6면체, 8면체, 12면체, 20면체 총 5개 주사위가 든 상자가 있다.
랜덤하게 하나를 집어서 던졌더니 주사위 눈이 6 이였다. 각 주사위일 확률은??
접근 방법
가설을 나타내자.
데이터를 나타내자.
우도 함수를 작성하자.
여기서 가설의 종류는 주사위 갯수인 5가지이다.
쉽게 4, 6, 8, 12, 20 가설로 표현하자.
그리고 데이터는 1~20가지의 숫자가 된다.
# 가설 등록
suite = Dice([4, 6, 8, 12, 20])
# 주사위눈 6이 나왔을때, 갱신 코드
suite.Update(6)
0.08500000000000002
sorted(suite.Items(), key=lambda (theta, prob): theta)
[(4, 0.0), (6, 0.3921568627450979), (8, 0.2941176470588235), (12, 0.19607843137254896), (20, 0.11764705882352941)]
해석해보자면,
데이터로 나온 주사위눈이 6임으로
4면체는 나올 수가 없어서 확률이 0 이되고,
나머지는 자신이 가진 눈의 갯수대로 n 등분한 값이 되어 나왔다.
#몇번 더 던져서 데이터로 6, 8, 7, 7, 5, 4 가 나왔다면.
for roll in [6, 8, 7, 7, 5, 4]:
suite.Update(roll)
sorted(suite.Items(), key=lambda (theta, prob): theta)
[(4, 0.0), (6, 0.0), (8, 0.9432484536722127), (12, 0.055206128061290875), (20, 0.001545418266496554)]
결과적으로 위와 같은 데이터가 나오게 하는 가장 높은 확률은 8면체 주사위이다.
프레드릭 모스텔러의 확률 고급문제 50선에 나오는 기관차 문제
"각 철도에는 이를 지나가는 기관차에 1부터 N까지의 순서로 번호를 붙인다. 어느 날 60호 기관차를 보았다. 이 때 이 철도에는 몇 개의 기관차가 지나가는지 추측해보자."
베이지안 추론법에 의하면 이런 문제는 2단계를 거친다.
데이터를 보기 전에 N 에 대해서 알고 있는 것은 무엇인가?
N 에 어떤 값이 주어졌을 때 관측한 데이터(60호 기관차) 의 우도는 어떻게 되는가?
첫번째 질문은 사전확률, 두번째 질문은 우도에 관한 것이다.
이 문제는 사전 확률을 선택하기에 근거가 부족하지만, 간단한 가정으로 설정해보자.
N 은 1 ~ 1000 이 동일한 확률로 선택될 수 있다고 하자.
hypos = xrange(1, 1001)
우도 함수를 만들어 보자.
N 개 기관차 중, 60호 기관차를 볼 확률은 얼마인가?
class Train(Suite):
def Likelihood(self, data, hypo):
if hypo < data:
return 0
else:
return 1.0/hypo
# 가설 등록
suite = Train(hypos)
# 데이터 적용 (60호 기관차)
suite.Update(60)
0.0028222671142652746
# 가설이 1000개 임으로 그래프로 살펴보자.
import matplotlib.pyplot as plt
xy = suite.Items()
xs, ys = zip(*xy)
plt.plot(xs, ys)
plt.show()
해석해보자면,
당연히 60개 미만인 경우에서 60호 기관차를 볼 수 없으므로 확률은 0.
나머지는 1/N 로 확률이 적용되었다.
그러면 60호까지 있는 N 의 확률이 가장 높은데... 원하는 것은 이것이 아닐것이다.
사후 확률의 평균을 계산해보자.
def Mean(suite):
total = 0
for hypo, prob in suite.Items():
total += hypo * prob
return total
print Mean(suite)
# suite 클래스 내에 이미 구현되어 있다.
print suite.Mean()
333.419893264 333.419893264
사후 평균은 333이므로 오차를 최소화 하자면 잘 추측한 것이다.
스무고개 계속하면 사후 평균을 사용해 장기적 관점에서의 평균제곱 오차를 최소화 할 수 있을 것이다.
기관차 문제 아직 안끝났다.
앞에서 우리는 기관차의 개수 N 을 1 ~ 1000 개로 정해놓고 풀었었다.
하지만, 사람마다 철도 회사에 기차를 1000대보다 많게, 또는 적게 추측할 것이다.
앞서 1~1000 까지 범위에서 사후확률 평균은 333 이였는데,
상한값이 500 이라면 207 이고,
상한값이 2000 이라면 552 가 될 것이다.
결과적으로 임의로 정한 기관차의 상한값에 민감하게 반응한다. (당연히 관찰 횟수가 1번이니)
suite500 = Train(xrange(1, 501))
suite1000 = Train(xrange(1, 1001))
suite2000 = Train(xrange(1, 2001))
suite500.Update(60)
suite1000.Update(60)
suite2000.Update(60)
print suite500.Mean()
print suite1000.Mean()
print suite2000.Mean()
207.079227983 333.419893264 552.179017165
그렇다면!!!
이후 두가지 방식으로 진행할 수 있다.
1. 데이터를 더 확보할 것
2. 배경지식을 더 확보할 것
데이터를 더 확보하면 전체적으로 수렴하는 경향을 보인다.
60호 기관차에 이어, 30호, 90호 기관차를 관측했다고 해보자.
for suite in [suite500, suite1000, suite2000]:
# 60호에 이어서 30호, 90호 관측!
for data in [30, 90]:
suite.Update(data)
print suite500.Mean()
print suite1000.Mean()
print suite2000.Mean()
151.849587959 164.305586423 171.338181092
하지만, 데이터를 더 확보할 수 없다면,
'배경지식'을 더 많이 수집해서 '사전 확률을 개선'하는 방법이 있다.
회사의 일반적인 규모에 대한 정보를 기관차 관련 전문가에게 얻을 수도 있다.
또, 이런 특화된 정보 없이 학습된 추측도 할 수 있다.
회사 규모의 분포는 '로버트 액스텔'의 논문에 따르면 '멱법칙' 을 따른다.
이 법칙에 따르면,
10대 미만의 회사가 1000 개이면,
100대 회사는 100개,
1000대 회사는 10개,
10000대 회사는 1개이다.
x 는 확률 누적 함수이고, $\alpha$ 는 보통 1에 가깝게 설정되는 매개변수이다.
from thinkbayes import Pmf
class Train2(Dice):
def __init__(self, hypos, alpha=1.0):
Pmf.__init__(self)
# 멱함수를 적용한다.
# 갯수가 많아 질수록 역수로 인해 확률은 낮아진다.
for hypo in hypos:
self.Set(hypo, hypo**(-alpha))
self.Normalize()
이제 새로운 클래스를 가지고 다시 계산해보자.
suite_list = [Train2(xrange(1, 501)),
Train2(xrange(1, 1001)),
Train2(xrange(1, 2001))]
for suite in suite_list:
for data in [30, 60, 90]:
suite.Update(data)
for suite in suite_list:
print suite.Mean()
130.708469863 133.275231375 133.997463081
차이가 많이 비슷해졌다. 사실 임의의 큰 상한값을 만들어두면 평균은 134 에 수렴한다.
따라서 멱법칙 사전확률이 보다 현실적이고, 실제로 더 잘 적용된다.
그러면 이전에 작성한 그래프와 같이 비교해보자.
suite1 = Train(xrange(1, 1001))
suite2 = Train2(xrange(1, 1001))
suite1.Update(60)
suite2.Update(60)
import matplotlib.patches as mpatches
xy = suite1.Items()
xs, ys = zip(*xy)
graph1, = plt.plot(xs, ys)
xy = suite2.Items()
xs, ys = zip(*xy)
graph2, = plt.plot(xs, ys)
plt.legend([graph1, graph2], ['uniform', 'power law'])
plt.show()
결과를 구간으로 요약하는 경우, 점으로 추정하는 경우(평균, 중간값, 최대 우도값, ...) 있다.
구간의 경우,
미확인 값이 두 값 사이에 올 확률이 90% 인 값을 나타내며, 이 범위를 '신뢰구간' 이라한다.
신뢰구간 계산하는 간단한 방법은 사후분포확률을 더한 후 5%, 95% 에 해당하는 5분위, 95분위 값을 기록하는 것이다.
# thinkbayes 에서 분위값 계산하는 함수를 제공해준다.
def Percentile(pmf, percentage):
p = percentage / 100.0
total = 0
for val, prob in pmf.Items():
total += prob
if total >= p:
return val
# 멱함수도 적용하고, 60, 30, 90 데이터까지 적용한 suite 로 살펴보자.
interval = Percentile(suite, 5), Percentile(suite, 95)
print interval
(91, 243)
가장 데이터가 많이 적용된 마지막 버전으로 확인했는데,
범위가 91 ~ 243 으로 너무 넓어서 그다지 확신할 수 없다는 것을 보여준다.
몇분위수 이상을 계산하는 경우에는, 누적분포함수나 Cdf 를 사용하는 것이 효율적이다.
Cdf 는 분포에 대해 동일한 정보를 가진다는 면에서 Pmf 와 동일하고, 언제나 서로바꿀 수 있다.
Cdf 의 장점은 분위수를 조금 더 효율적으로 계산한다는 점.
누적분포함수 계산을 위해 클래스 Cdf 를 사용하자.
from thinkbayes import Cdf
# pmf 를 cdf 로 변환해서
cdf = suite.MakeCdf()
# 분위수 계산을 더 쉽게!
interval = cdf.Percentile(5), cdf.Percentile(95)
참고로 pmf => cdf 변환은 갯수에 비례하는 len(pmf) 을 사용함으로 O(n) 시간이 든다.
cdf 는 확률을 찾는데 'O(log n) 로그 시간' 이 소요된다. 그래서 효율적!
2차 세계 대전 동안, 독일 탱크 장비 생산량을 추정하는데 통계분석기법을 사용했다.
서방 연합국에서는 각 탱크의 섀시와 엔진의 시리얼 번호를 포함한 기록, 수리 기록등을 확보했다.
분석해보니,
시리얼 번호는 제조자에 따라 할당되었고,
탱크의 유형은 100개 숫자 블록으로 이루어져 있는데 숫자가 순차적으로 사용되었으나 모든 숫자가 사용된 것은 아니었다.
따라서 독일 탱크 생산량 추정은 100개 숫자를 포함하는 블록 범위 내에 기관차 문제 형태로 단순화 된다.
이후 분석가들은 다른 정보원으로부터 이보다 낮은 추정치를 도출했고, 결과는 대체적으로 맞았다.
타이어, 트럭, 로켓, 다른 장비에 대해서도 유사한 분석을 해서 경제 기밀 정보를 도출했다.
베이지안을 통해 이 실전 문제로 향하는 간단한 길을 제공해준다.
베이지안이 선행분포를 선택하는 두가지 방법
배경지식을 잘 활용
단점은 사람마다 서로 다른 배경지식을 사용할 수 있다는 것이다.
최대한 제약을 두지 않는 것
최소한의 선행 정보를 사용하거나 일부 필요한 성질을 나타내는 고유 선행 분포를 정의할때 사용
좀 더 객관적.
저자의 경우, 정보적 선행 분포를 선호!
왜냐면, 첫째, 베이지안 분석은 항상 모델 판단에 기반해서 이루어진다.
따라서 비정보적 선행분포가 더 객관적이라해도 어차피 전체 분석은 주관적으로 이루어진다.
둘째, 실질적 문제는 보통 데이터가 굉장이 많거나, 별로 없거나 하는 2가지 유형이다.
만약 데이터가 굉장히 많다면 선행분포를 어떤 쪽으로 선택하든지 별로 상관 없고,
정보성이든 비정보성이든 거의 유사한 결과가 도출된다.
하지만, 기관차 문제처럼 데이터가 거의 없을때, 배경지식의 사용여부는 매우 큰 차이를 나타낸다.
또 독일 탱크 문제처럼 결과에 따라 생사가 달라진다면, 객관성보다는 알고있는 최대한의 지식을 쏟아 넣어야 할 것이다.
기관차 문제에서 우도함수를 만들기 위해서는 아래 질문에 답해야한다.
> "기찻길에 N 개의 기관차가 있을때, 60호 기관차를 볼 확률은 얼마인가?"