Think Bayes (14) 계층 모델
 파이썬을 활용한 베이지안 통계, 앨런 B. 다우니 지음, 권정민 옮김/한빛미디어 |
14.1 가이거 계수기 문제
다음 문제를 생각해보자.
방사성 물질이 가이거 계수기에서 평균 초당 r개의 입자를 발산한다고 가정할 때 계수기는 반응하는 입자의 분수 f만을 등록한다. 만약 f가 10%이고 계수기가 1초 간격으로 15개의 입자를 등록햇다면 실제 계수기에 충돌한 입자의 실제 개수 n의 사후 분포와 평균 방출 입자 비율 r은 어떻게 될까?
이런 문제를 풀기 위해, 시스템 변수로부터 시작해서 관측 데이터로 끝나기까지의 과정을 살펴보자.
- 방사성 물질은 평균 r의 비율로 입자를 방출한다.
- 주어진 초 동안, 방사성 물질은 계수기에 n개의 입자를 방출한다.
- n개의 입자 중 일부인 k개가 기록된다.
원자의 붕괴 확률은 어떤 시간 지점에서든 동일하므로 포아송 프로세스를 모델링할 수 있다. r에 대해서 n의 분포는 변수 r에 대한 포아송 분포이다.
각 입자의 대한 감지 확률은 다른 입자의 확률과 독립적이라고 가정하면, k의 분포는 n과 f의 변수를 갖는 이항 분포이다.
시스템의 변수가 주어지면, 데이터의 분포를 알 수 있다. 이를 정문제(forward problem)이라고 한다. 반대로 데이터가 주어지고 이에 대한 변수의 분포를 알고자 하는 것을 역문제(backward problem)이라고 한다.
정문제를 풀 수 있다면, 베이지안 방식을 사용해서 역문제를 풀 수 있다.
14.2 단순하게 시작하기
r의 값을 알고 있다면 문제를 단순하게 생각하여 시작할 수 있다. 주어진 f의 값을 알고 있으므로 n을 추정하기만 하면 된다. n을 추정하는 Detector 스윗을 정의하자.
|
1 2 3 4 5 6 7 |
class Detector(thinkbayes.Suite): def __init__(self, r, f, high=500, step=1): pmf = thinkbayes.MakePoissonPmf(r, high, step=step) thinkbayes.Suite.__init__(self, pmf, name=r) self.r = r self.f = f |
평균 방출 비율이 초 당 입자 r개라면, n의 분포는 r을 변수로 갖는 포아송 분포다. high와 step은 n의 상한 값과 가설 값 사이의 단계 크기를 결정한다. 우도 함수를 정의하자.
|
1 2 3 4 5 6 7 8 |
# class Detector def Likelihood(self, data, hypo): k = data n = hypo p = self.f return thinkbayes.EvalBinomialPmf(k, n, p) |
data는 검출된 입자의 수이고, hypo는 방출된 입자의 가설 수 n이다. n개의 입차가 있을 때 이 중 하나가 검출된 확률은 f로, k개의 입자가 검출될 확률은 이항 분포라고 가정하였다.
이제 실행해보자.
|
1 2 3 4 5 6 7 8 |
f = 0.1 k = 15 for r in [100, 250, 400]: suite = Detector(r, f, step=1) suite.Update(k) print suite.MaximumLikelihood() |
14.3 계층적으로 만들기
앞에서 미리 정해두었던 r을 이제는 모른다고 가정하고 방출체의 동작을 모델링하고 r을 추정하는 Emitter 스윗을 정의하자.
|
1 2 3 4 5 |
class Emitter(thinkbayes.Suite): def __init__(self, rs, f=0.1) detectors = [Detector(r, f) for r in rs] thinkbayes.Suite.__init__(self, detectors) |
rs는 r의 가설 값으로 사용하려는 값이다. detectors는 각 r에 대한 Detector 객체들이다.
이 Emitter를 업데이트 하기 위해서 모든 가설 r에 대하여 우도를 계산하여 더하여야 한다.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
# class Detector def SuiteLikelihood(self, data): total = 0 for hypo, prob in self.Items() like = self.Likelihood(data, hypo) total += prob * like return total # class Emitter def Likelihood(self, data, hypo): detector = hypo like = detector.SuiteLikelihood(data) return like def Update(self, data): thinkbayes.Suite.Update(self, data) for detector in self.Values(): detector.Update() |
14.4 약간 최적화하기
사실 SuiteLikelihood가 계산하는 확률은 Update에서 정규화되어 업데이트되었으므로 실제로는 필요하지 않다. 따라서 Emitter를 업데이트한 후 Detector를 업데이트하는 대신 Detector.Update의 결과를 Emitter의 우도로 사용해서 한번에 계산할 수 있다. 다음은 수정된 함수이다.
|
1 2 3 4 |
# class Emitter def Likelihood(self, data, hypo): return hypo.Update(data) |
14.5 사후 분포 추출하기
Emitter 업데이트 후 Detector와 각각의 확률을 반복하면서 r의 사후 분포를 구할 수 있다.
|
1 2 3 4 5 6 7 8 9 |
# class Emitter def DistOfR(self): items = [(detector.r, prob) for detector, prob in self.Items()] return thinkbayes.MakePmfFromItems(item) def DistOfN(self): return thinkbayes.MakeMixture(self) |
items는 r의 값과 이에 해당하는 확률의 리스트다. 만들어진 값은 r에 대한 PMF이다.
n의 사후 확률 분포는 Detector의 혼합 확률을 계산하여 구한다.
결과로 가장 가능성 높은 n의 값은 150이다. 그리고 150개의 입자가 1초에 방출된다면 r의 가장 가능선 높은 값 또한 초당 150개이다.
r과 n의 사후 분포는 유사하지만 상대적으로 n에 대해서 확신도가 조금 떨어진다.
사용한 코드는 http://thinkbayes.com/jaynes.py
14.6 토의
방출 비율 r은 입자의 수 n에 대한 인과 효과를 가지고, n은 입자 수치 k에 대한 인과 효과를 가진다.
계층 모델은 상위에서 원인을 제공하면 하위에서 영향을 받는 시스템의 구조를 반영한다.
- 최상위 단계에서 r의 가설 값의 범위로부터 시작한다.
- 각 r 값에 대해 n의 값의 범위를 구하고 r에 따른 n의 사전 분포를 구한다.
- 모델을 갱신할 때는 아래에서 위로 올라온다. 각 r의 값에 대한 n의 사후 분포를 구한 후 r의 사후 분포를 계산한다.
인과 정보는 계층을 따라 내려오게 되고 추론은 위로 향하게 된다.