Notebook

A Hierarchical Model¶

계층적 모형의 특징: Hyperprior를 도입함으로써 개인차가 기본적으로 집단 수준 분포에서 옴을 반영하면서도, 개인 수준에서 차이가 있을 수 있음을 인정하는 것

14.1 가이거 계수기 문제¶

방사성 물질이 가이거 계수기에서 평균 초당 r개의 입자를 발산한다고 가정할 때 계수기는 반응하는 입자의 분수 f만을 등록한다. 만약 f가 10%이고 계수기가 1초 간격으로 15개의 입자를 등록했다면 실제 계수기에 충돌한 입자의 실제 개수 n의 사후 분포와 평균 방출 입자 비율 r은 어떻게 될까?
이 문제를 풀기 위해 시스템 변수로부터 시작해서 관측 데이터로 끝나기까지 인과관계를 생각해보자.

가이거 계수기는 손으로 들고 다닐 수 있어 널리 사용되는 방사능 측정장비이다. 불활성 기체를 담은 가이거-뮐러 계수관을 이용하여 알파 입자, 베타 입자, 감마선과 같은 방사능에 의해 불활성 기체가 이온화되는 정도를 표시하여 방사능을 측정한다. (위키피티아: https://ko.wikipedia.org/wiki/%EA%B0%80%EC%9D%B4%EA%B1%B0_%EA%B3%84%EC%88%98%EA%B8%B0)

방사성 물질은 평균 r의 비율로 방출한다.
주어진 초 동안 방사성 물질은 계수기에 n개의 입자를 방출한다.
n개의 입자 중 일부인 k개가 기록된다.

원자의 붕괴 확률은 어떤 시간 지점에서든 동일하므로 방사능 동위 원소 붕괴도 포아송 프로세스로 모델링 할 수 있다. R에 대하여 n의 분포는 변수 r에 대한 포아송 분포다.

그리고 각 입자에 대한 감지 확률은 다른 입자의 확률과 독립적이라고 가정하면 k의 분포는 n과 f의 변수를 갖는 이항분포다.

시스템 변수가 주어지면, 데이터의 분표를 알 수 있다. 이를 정방향 문제(forward problem)라 한다.

반대로 데이터가 주어지고 이에 대한 변수의 분포를 알고자 하는 것이 있는데 이를 역방향 문제(inverse problem)이라고 한다. 그리고 정문제를 풀 수 있다면 베이지안 방식을 사용해서 역문제를 풀 수 있다.

14.2 단순하게 시작하기¶

변수 r의 값을 알고 있다면 문제의 단순한 형태를 시작할 수 있다. 주어진 f의 값을 알고 있으므로 n을 추정하기만 하면 된다.

이를 위해 검출기의 행동을 모델링하고 n을 추정하는 Detector스윗을 정의한다.

In [1]:

import thinkbayes
import thinkplot

from math import exp

FORMATS = ['pdf', 'eps', 'png']

In [2]:

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의 상한값과 가설 값 사이의 단계 크기를 결정한다.

우리는 이제 우도함수(likelihood function)가 필요하다.

In [3]:

class Detector(thinkbayes.Suite):

    def Likelihood(self, data, hypo):
        k = data
        n = hypo
        p = self.f

        return thinkbayes.EvalBinomialPmf(k, n, p)

데이터는 검출된 입자의 수이고 hypo 는 방출된 입자의 가설 수 n이다.

실제로 n개의 입자가 있을 때 이 중 하나가 검출될 확률은 f로 k개의 입자가 검출될 활귤은 이항분포로 주어진다.

그러면 Detector를 위해 필요한 것은 다 준비되었다. 이제 r의 값 범위 내에서 실행해볼 수 있다.

In [4]:

    f = 0.1
    k = 15  
   

    for r in [100, 250, 400]:
        suite = Detector(r, f, step=1)
        suite.Update(k)
        print suite.MaximumLikelihood()

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-4-da803ef36b00> in <module>()
      4 
      5 for r in [100, 250, 400]:
----> 6     suite = Detector(r, f, step=1)
      7     suite.Update(k)
      8     print suite.MaximumLikelihood()

TypeError: __init__() got an unexpected keyword argument 'step'

그럼 14-1 에는 주어진 몇 가지 r에 대한 n의 사후 분포가 나타나 있다. r과 n의 사후 분포는 유사하다. 유일한 차이점은 상대적으로 n에 대해서는 확신도가 아주 조금 떨어진다는 것이다. 특정 초에 방출될 입자의 수 n보다 일반적으로 장기간 방출 비율 r에 대해 더 확신할 수 있다.

그림 14-1 세 가지 r값에 대한 n의 사후 분포

14.3 계층적으로 만들기¶

앞에서 r이 알려져 있다고 가정했다. 그러면 이제 이 가정을 접어두자. 여기서는 방출체의 동작을 모델링하고 r을 추정하는 Emitter스윗을 정의한다.

In [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 는 Detector 객체(objects)의 수열(sequence)로 각각 r의 값 안에 들어 있다. 스윗의 값은 Detector로 Emitter는 meta-Suite이다. 메타 스윗은 다른 스윗을 값으로 가지는 것이다.

Emitter 갱신(update)을 위해서는 각 r에 대한 가설 값하의 데이터의 우도를 계산해야 한다. 하지만 각 r의 값은 n의 범위를 가지는 Detector로 나타난다. 주어진 Detector의 데이터의 우도 계산을 위해서는 n의 값에 대해 반복문을 돌려서 k의 건체 확률을 더해야 한다. 이 기능은 SuiteLikelihood가 수행한다.

In [6]:

class Detector(thinkbayes.Suite):

    def SuiteLikelihood(self, data):
        total = 0
        for hypo, prob in self.Items():
            like = self.Likelihood(data, hypo)
            total += prob * like
        return total

이제 Emitter의 Likelihood 함수를 작성할 수 있다.

In [7]:

class Detector(thinkbayes.Suite):

    def Likelihood(self, data, hypo):
        detector = hypo
        like = detector.SuiteLikelihood(data)
        return like

각 hypo는 Detector로 가설하의 데이터의 우도를 얻기 위해서 SuiteLikelihood를 호출할 수 있다.

Emitter갱신 후에는 각 Detector도 갱신해 주어야 한다.

In [8]:

class Detector(thinkbayes.Suite):

    def Update(self, data):
        thinkbayes.Suite.Update(self, data)
        
        for detector in self.Values():
            detector.Update()

이 모델처럼 Suites이 multiple levels로 사용되는 것을 계층적(hierarchical)이라고 한다.

14.4 약간 최적화하기¶

SuiteLikelihood가 계산하는 전체 확률은 완전히 정규화된 상수(normalizing constant)가 되어서 Update된 것이므로 실제로는 필요하지 않다. .

Emitter를 updating한 후에 Detector를 updating하는 방법 대신 Detector.Update의 결과를 Emitter의 likelihood로 사용해서 한 번에 계산할 수 있다.

다음은 Emitter.Likelihood이 새 버전이다.

In [9]:

class Emitter(thinkbayes.Suite):

    def Likelihood(self, data, hypo):
        return hypo.Update(data)

Likelihood의 새 버전을 사용해서 Update의 기본 버전을 사용할 수 있다. 따라서 이 버전은 코드 줄 수가 더 적고 상수 정규화(normalizing constant)를 두 번 하지 않으므로 더 빨리 돌아간다.

14.5 사후 분포( posteriors) 추출하기¶

Emitter를 갱신한 후, Detector와 각각의 확률을 반복하면서 r의 사후 확률 분포를 구할 수 있다.

In [10]:

class Emitter(thinkbayes.Suite):

    def DistOfR(self):
        items = [(detector.r, prob) for detector, prob in self.Items()]
        return thinkbayes.MakePmfFromItems(items)

items는 r의 값과 이에 해당하는 확률의 리스트이다. 결과는 r의 Pmf이다.

n의 사후 분포를 구하려면 Detector의 혼합을 계산해야 한다. 각 분포를 이에 해당하는 확률에 연결하는 meta-Pmf that maps from each distribution to its probability. And thatâs exactly what the Emitter is:

In [11]:

class Emitter(thinkbayes.Suite):

    def DistOfN(self):
        return thinkbayes.MakeMixture(self)

그림은 아래에 나와 있다. 예상했겠지만 가장 가능성 높은 n의 값은 150이다. 주어진 f와 n에 대해서 예상 수치는 k=fn으로 구하므로 주어진 f와 k에 대해 n의 계상 수치는 k/f로 150이 된다.

그리고 만약 150개의 입자가 1초에 방출돈다면 r의 가장 가능성 높은 값은 초당 150개의 입자다.

r과 n의 사후 분포는 유사하다. 유일한 차이점은 상대적으로 n에 대해서는 확신도가 아주 조금 떨어진다는 것이다. 특정 초에 방출될 입자의 수 n보다 일반적으로 장기간 방출 비율 r에 대해 더 확신할 수 있다.

토의¶

최상위 단계에서 r의 가설 값의 범위로부터 시작한다.
각 r 값에 대해 n의 값의 범위를 구하고 r에 따른 n의 사전 분포를 구한다.
모델을 갱신할 때는 아래에서 위로 올라간다. 각 r의 값에 대한 n의 사후 분포를 구한 후 r의 사후 분포를 계산한다.

즉 인과 정보는 계층을 따라서 내려오게 되고 추론은 위로 향하게 된다.

""" http://www.greenteapress.com/thinkbayes/jaynes.py%EC%97%90%EC%84%9C 가져온 코드 This file contains code used in "Think Stats", by Allen B. Downey, available from greenteapress.com

import thinkbayes import thinkplot

from math import exp

FORMATS = ['pdf', 'eps', 'png']

class Emitter(thinkbayes.Suite): """Represents hypotheses about r."""

def __init__(self, rs, f=0.1):
    """Initializes the Suite.

    rs: sequence of hypothetical emission rates
    f: fraction of particles registered
    """
    detectors = [Detector(r, f) for r in rs]
    thinkbayes.Suite.__init__(self, detectors)

def Update(self, data):
    """Updates the Suite based on data.

    data: number of particles counted
    """
    thinkbayes.Suite.Update(self, data)

    for detector in self.Values():
        detector.Update()

def Likelihood(self, data, hypo):
    """Likelihood of the data given the hypothesis.

    Args:
        data: number of particles counted
        hypo: emission rate, r

    Returns:
        probability density of the data under the hypothesis
    """
    detector = hypo
    like = detector.SuiteLikelihood(data)
    return like

def DistOfR(self, name=''):
    """Returns the PMF of r."""
    items = [(detector.r, prob) for detector, prob in self.Items()]
    return thinkbayes.MakePmfFromItems(items, name=name)

def DistOfN(self, name=''):
    """Returns the PMF of n."""        
    return thinkbayes.MakeMixture(self, name=name)

class Emitter2(thinkbayes.Suite): """Represents hypotheses about r."""

def __init__(self, rs, f=0.1):
    """Initializes the Suite.

    rs: sequence of hypothetical emission rates
    f: fraction of particles registered
    """
    detectors = [Detector(r, f) for r in rs]
    thinkbayes.Suite.__init__(self, detectors)

def Likelihood(self, data, hypo):
    """Likelihood of the data given the hypothesis.

    Args:
        data: number of counted per unit time
        hypo: emission rate, r

    Returns:
        probability density of the data under the hypothesis
    """
    return hypo.Update(data)

def DistOfR(self, name=''):
    """Returns the PMF of r."""
    items = [(detector.r, prob) for detector, prob in self.Items()]
    return thinkbayes.MakePmfFromItems(items, name=name)

def DistOfN(self, name=''):
    """Returns the PMF of n."""        
    return thinkbayes.MakeMixture(self, name=name)

class Detector(thinkbayes.Suite): """Represents hypotheses about n."""

def __init__(self, r, f, high=500, step=5):
    """Initializes the suite.

    r: known emission rate, r
    f: fraction of particles registered
    high: maximum number of particles, n
    step: step size between hypothetical values of n
    """
    pmf = thinkbayes.MakePoissonPmf(r, high, step=step)
    thinkbayes.Suite.__init__(self, pmf, name=r)
    self.r = r
    self.f = f

def Likelihood(self, data, hypo):
    """Likelihood of the data given the hypothesis.

    data: number of particles counted
    hypo: number of particles hitting the counter, n
    """
    k = data
    n = hypo
    p = self.f

    return thinkbayes.EvalBinomialPmf(k, n, p)

def SuiteLikelihood(self, data):
    """Adds up the total probability of the data under the suite.

    data: number of particles counted
    """
    total = 0
    for hypo, prob in self.Items():
        like = self.Likelihood(data, hypo)
        total += prob * like
    return total

def main(): k = 15 f = 0.1

# plot Detector suites for a range of hypothetical r
thinkplot.PrePlot(num=3)
for r in [100, 250, 400]:
    suite = Detector(r, f, step=1)
    suite.Update(k)
    thinkplot.Pmf(suite)
    print suite.MaximumLikelihood()

thinkplot.Save(root='jaynes1',
               xlabel='Number of particles (n)',
               ylabel='PMF',
               formats=FORMATS)

# plot the posterior distributions of r and n
hypos = range(1, 501, 5)
suite = Emitter2(hypos, f=f)
suite.Update(k)

thinkplot.PrePlot(num=2)
post_r = suite.DistOfR(name='posterior r')
post_n = suite.DistOfN(name='posterior n')

thinkplot.Pmf(post_r)
thinkplot.Pmf(post_n)

thinkplot.Save(root='jaynes2',
               xlabel='Emission rate',
               ylabel='PMF',
               formats=FORMATS)

if name == 'main': main()

In [ ]: