from pymc import DiscreteUniform, Exponential, Deterministic, Poisson, Uniform
import numpy as np
import disaster_model
Couldn't import dot_parser, loading of dot files will not be possible.
disasters_array = \
np.array([ 4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,
1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
disaster_model.py에 다 있으니 추가해서 사용하시오
switchpoint = DiscreteUniform('switchpoint', lower=0, upper=110, doc='Switchpoint[year]')
early_mean = Exponential('early_mean', beta=1.)
late_mean = Exponential('late_mean', beta=1.)
다음으로
@deterministic(plot=False)
def rate(s=switchpoint, e=early_mean, l=late_mean):
''' Concatenate Poisson means '''
out = np.empty(len(disasters_array))
out[:s] = e
out[s:] = l
return out
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-5-4530e4e394bc> in <module>() ----> 1 @deterministic(plot=False) 2 def rate(s=switchpoint, e=early_mean, l=late_mean): 3 ''' Concatenate Poisson means ''' 4 out = np.empty(len(disasters_array)) 5 out[:s] = e NameError: name 'deterministic' is not defined
disasters = Poisson('disasters', mu=rate, value=disasters_array, observed=True)
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-6-cfa065042c3b> in <module>() ----> 1 disasters = Poisson('disasters', mu=rate, value=disasters_array, observed=True) NameError: name 'rate' is not defined
우리가 disasters를 관측한 이후에나 이 관측 정보가 변수와 어떻게 관련이 있는지에 확인할 수 있는 정보이기 때문에 early_mean, switchpoint, late_mean에 대해 추론하기 위해서 이런 프로세스 모델을 사용할 필요가 생긴다.
from pymc.examples import disaster_model
disaster_model.switchpoint.parents
{'lower': 0, 'upper': 110}
disaster_model.disasters.parents
{'mu': <pymc.PyMCObjects.Deterministic 'rate' at 0x7fd2d1403610>}
disaster_model.rate.children
{<pymc.distributions.Poisson 'disasters' at 0x7fd2d1403750>}
disasters는 rate를 부모로 섬기기 때문에 rate는 disaster를 새끼로 키운다.
parents와 달리 children은 무순위의 객체 집단의 집합이다.
변수는 특별한 분포적 역할을 담당하는 그들의 자식과 연관되지 않는다.
모델에 다른 파라미티의 parents와 children 속성을 평가해 보라.
아래의 방향성 비순환성 그래프는 모델에서 parent-child 관계의 가시화이다.
관측되지 않은 stochastic 변수인 stochastic switchpoint, early_mean, late_mean들은 타원으로, 관측된 stochastic 변수인 disasters는 채워진 타원으로, deterministic 변수인 rate는 삼각형으로 표현된다.
부모와 자식의 관계는 화살표로 나타내고 부모에 할당된 자식은 레이블로 표시한다.
Graphing models을 상세히 보자(4장에 있음)
석탄광 재난 모델 예에서 방향을 갖는 비순환 관계 그래프이다.
위 예에서 보듯이 pymc 객체는 switchpoint, early_mean, late_mean처럼 이름이 정해져야 한다. 이 이름은 저장과 사후 처리를 위해 사용된다.
동일한 이름을 갖는 변수로 인스턴스된 모델은 이름 충돌을 피하기 위해 db저장 추적시에 err을 낸다. 그러나 일반적으로 pymc는 변수를 식별하기 위해서 객체 자체의 참조를 사용하지 이름을 사용하진 않는다.
disaster_model.disasters.value
array([4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6, 3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5, 2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0, 1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2, 3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
switchpoint, early_mean, late_mean의 값을 확인하면 PyMC가 생성한 랜덤 초기값을 볼 수 있다.
disaster_model.switchpoint.value
array(24)
disaster_model.early_mean.value
array(0.36637928740655185)
disaster_model.late_mean.value
array(1.3932648279262168)
disaster_model.rate.value
array([ 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 0.36637929, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483, 1.39326483])
이 값을 계산하기 위해서 rate는 부모의 값에 넘기는 우리가 만든 함수를 호출한다.
Stochastic 객체는 부모의 값에 의한 현재 값에서 그들의 확률 질량 함수나 확률 밀도 함수를 평가할 수 있다.
Stochastic 객체의 확률 질량/밀도 함수의 로그치는 logp 속성으로 접근할 수 있다.
disasters 같은 벡터 벨류 변수를 위해서는 logp 속성은 모든 값 인자의 결합 확률이나 밀도의 로그 합을 반환한다.
switchpoint의, disasters의 로그 확률과 early_mean의, late_mean의 로그 밀도를 평가해 봐라.
disaster_model.switchpoint.logp
-4.709530201312334
disaster_model.disasters.logp
-284.3224944872588
disaster_model.early_mean.logp
-0.36637928740655185
disaster_model.late_mean.logp
-1.3932648279262168
rate에 대한 우리의 정의를 밀착해서 보자
@deterministic(plot=False)
def rate(s=switchpoint, e=early_mean, l=late_mean):
''' Concatenate Poisson means '''
out = np.empty(len(disasters_array))
out[:s] = e
out[s:] = l
return out
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-20-4530e4e394bc> in <module>() ----> 1 @deterministic(plot=False) 2 def rate(s=switchpoint, e=early_mean, l=late_mean): 3 ''' Concatenate Poisson means ''' 4 out = np.empty(len(disasters_array)) 5 out[:s] = e NameError: name 'deterministic' is not defined
switchpoint, early_mean, late_mean은 stochastic 객체이지 수치가 아니다.
만약 수치라면 우리가 배열 out을 stochastic 객체에 썰어 넣을 때 오류가 안나겠는가?
변수는 자식변수를 위한 부모 변수로 쓰일 때마다 자식 값이나 log 확률이 계산되는 순간 PyMC는 변수를 값의 속성을 교체한다.
rate 값이 재계산되면 s.value는 switchpoint 형태로 함수에 전달된다.
rate의 부모의 값을 한꺼번에 보려면 rate.parents.value를 봐라.
from pymc.examples import disaster_model
from pymc import MCMC
M = MCMC(disaster_model)
이 경우 M은 switchpoint, early_mean, late_mean, disasters 변수를 속성으로 노출한다. M.switchpoint는 disaster_model.switchpoint와 같다.
샘플러를 돌리기 위해 iterations 수, burn-in 길이, thining interval을 아규먼트로 해서 MCMC객체의 sample()를 호출해라.
M.sample(iter=10000, burn=10000, thin=10)
[-----------------100%-----------------] 10000 of 10000 complete in 1.6 sec
수초 후에 표본이 정상적으로 추출됨을 볼 수 있다 모델이 적합된 거다.
M.trace('switchpoint')[:]
array([0])
%matplotlib inline
from pylab import hist, show
hist(M.trace('late_mean')[:])
show()
late_mean 파라미터의 마지널 사후 확률의 히스토그램.
from pymc.Matplot import plot
plot(M)
Plotting switchpoint Cannot plot autocorrelation for switchpoint Plotting early_mean Cannot plot autocorrelation for early_mean Plotting late_mean Cannot plot autocorrelation for late_mean
/usr/lib/pymodules/python2.7/matplotlib/axes.py:4511: RuntimeWarning: invalid value encountered in true_divide c /= np.sqrt(np.dot(x, x) * np.dot(y, y))
switchpoint를 위한 임시 시리즈, 자기상관 그림과 히스토 그램을 그린다.
사후 확률의 비 그래프적인 요약은 단지 M.stats()를 호출해라.
M.stats()
Could not generate output statistics for switchpoint Could not generate output statistics for early_mean Could not generate output statistics for late_mean
{'early_mean': None, 'late_mean': None, 'rate': {'95% HPD interval': array([ 0., 0.]), 'mc error': 0.0, 'mean': 0.0, 'n': 111, 'quantiles': {2.5: 0.0, 25: 0.0, 50: 0.0, 75: 0.0, 97.5: 0.0}, 'standard deviation': 0.0}, 'switchpoint': None}
Count | Site | Observer | Temperature |
---|---|---|---|
15 | 1 | 1 | 15 |
10 | 1 | 2 | NA |
0 | 1 | 1 | 15 |
x = np.array([ 4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
2, 2, 3, 4, 2, 1, 3, None, 2, 1, 1, 1, 1, 3, 0, 0,
1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
3, 3, 1, None, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
masked_values = np.ma.masked_equal(x, value=None)
masked_values
masked_array(data = [4 5 4 0 1 4 3 4 0 6 3 3 4 0 2 6 3 3 5 4 5 3 1 4 4 1 5 5 3 4 2 5 2 2 3 4 2 1 3 -- 2 1 1 1 1 3 0 0 1 0 1 1 0 0 3 1 0 3 2 2 0 1 1 1 0 1 0 1 0 0 0 2 1 0 0 0 1 1 0 2 3 3 1 -- 2 1 1 1 1 2 4 2 0 0 1 4 0 0 0 1 0 0 0 0 0 1 0 0 1 0 1], mask = [False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False True False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False False True False False False False False False False False False False False False False False False False False False False False False False False False False False False], fill_value = ?)
from pymc import Poisson
disasters = Poisson('disasters', mu=rate, value=masked_values, observed=True)
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-30-8360b20668d1> in <module>() 1 from pymc import Poisson ----> 2 disasters = Poisson('disasters', mu=rate, value=masked_values, observed=True) NameError: name 'rate' is not defined
전체 모델은 원 모델에 매우 유사해 보인다.
# Switchpoint
switch = DiscreteUniform('switch', lower=0, upper=110)
# Early mean
early_mean = Exponential('early_mean', beta=1)
# Late mean
late_mean = Exponential('late_mean', beta=1)
@deterministic(plot=False)
def rate(s=switch, e=early_mean, l=late_mean):
"""Allocate appropriate mean to time series"""
out = np.empty(len(disasters_array))
# Early mean prior to switchpoint
out[:s] = e
# Late mean following switchpoint
out[s:] = l
return out
# The inefficient way, using the Impute function:
# D = Impute('D', Poisson, disasters_array, mu=r)
#
# The efficient way, using masked arrays:
# Generate masked array. Where the mask is true,
# the value is taken as missing.
masked_values = masked_array(disasters_array, mask=disasters_array==-999)
# Pass masked array to data stochastic, and it does the right thing
disasters = Poisson('disasters', mu=rate, value=masked_values, observed=True)
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-32-a3b28b872d00> in <module>() 6 late_mean = Exponential('late_mean', beta=1) 7 ----> 8 @deterministic(plot=False) 9 def rate(s=switch, e=early_mean, l=late_mean): 10 """Allocate appropriate mean to time series""" NameError: name 'deterministic' is not defined
누락 자료의 추적기록, 자기상관도와 사후 분포의 예이다.
M.step_method_dict[disaster_model.switchpoint]
[<pymc.StepMethods.DiscreteMetropolis at 0x7fd2d13a4910>]
M.step_method_dict[disaster_model.early_mean]
[<pymc.StepMethods.Metropolis at 0x7fd2d13a4b10>]
M.step_method_dict[disaster_model.late_mean]
[<pymc.StepMethods.Metropolis at 0x7fd2d13a4b90>]
from pymc import Metropolis
M.use_step_method(Metropolis, disaster_model.late_mean, proposal_sd=2.)