2014-03-27 3 views
3

나는 포아송 과정의 비율을 추산하려고한다. 최대 비율의 사후 견적을 사용하여 시간에 따라 비율이 변한다. 다음은 선형으로 변화하는 비율 (λ = ax + b)의 간단한 예입니다.변하기 쉬운 비율로 덮인 포아송 과정을 붙이는 것

import numpy as np 
import pymc 

# Observation 
a_actual = 1.3 
b_actual = 2.0 
t = np.arange(10) 
obs = np.random.poisson(a_actual * t + b_actual) 

# Model 
a = pymc.Uniform(name='a', value=1., lower=0, upper=10) 
b = pymc.Uniform(name='b', value=1., lower=0, upper=10) 


@pymc.deterministic 
def linear(a=a, b=b): 
    return a * t + b 

r = pymc.Poisson(mu=linear, name='r', value=obs, observed=True) 

model = pymc.Model([a, b, r]) 
map = pymc.MAP(model) 
map.fit() 
map.revert_to_max() 

print "a :", a._value 
print "b :", b._value 

이것은 잘 동작합니다. 그러나 나의 실제 포아송 과정은 결정 론적 가치에 의해 제한됩니다. 나는 결정적 함수 내 관찰 된 값을 연결할 수 없기 때문에, 나는 나의 관찰을위한 작은 변화와 정상 확률 기능을 추가 해요 :

import numpy as np 
import pymc 

# Observation 
a_actual = 1.3 
b_actual = 2.0 
t = np.arange(10) 
obs = np.random.poisson(a_actual * t + b_actual).clip(0, 10) 

# Model 
a = pymc.Uniform(name='a', value=1., lower=0, upper=10) 
b = pymc.Uniform(name='b', value=1., lower=0, upper=10) 


@pymc.deterministic 
def linear(a=a, b=b): 
    return a * t + b 

r = pymc.Poisson(mu=linear, name='r') 


@pymc.deterministic 
def clip(r=r): 
    return r.clip(0, 10) 

rc = pymc.Normal(mu=r, tau=0.001, name='rc', value=obs, observed=True) 

model = pymc.Model([a, b, r, rc]) 
map = pymc.MAP(model) 
map.fit() 
map.revert_to_max() 

print "a :", a._value 
print "b :", b._value 

다음과 같은 오류 생산이 코드 :

Traceback (most recent call last): 
    File "pymc-bug-2.py", line 59, in <module> 
    map.revert_to_max() 
    File "pymc/NormalApproximation.py", line 486, in revert_to_max 
    self._set_stochastics([self.mu[s] for s in self.stochastics]) 
    File "pymc/NormalApproximation.py", line 58, in __getitem__ 
    tot_len += self.owner.stochastic_len[p] 
KeyError: 0 

내가 뭘 잘못하고 있는지에 대한 아이디어가 있습니까?

답변

3

"뚜껑 달린 소리"는 잘린 포아송임을 의미합니까? 당신이 말하는 그게 나타납니다. 왼쪽 잘림 (더 일반적인) 인 경우 TruncatedPoisson 배포를 사용할 수 있지만 오른쪽 잘림을 수행하고 있기 때문에 더 일반적으로 만들었어야합니다. 노력하고있는 것은 작동하지 않습니다. Poisson 객체에는 clip() 메소드가 없습니다. 당신이 할 수있는 것은 잠재적 인 요소를 사용하는 것입니다. 그 결과는 다음과 같습니다이 10 미만이 Potential 클래스에 대한 정보 pymc docs를 참조 할 r의 값을 제한합니다

@pymc.potential 
def clip(r=r): 
    if np.any(r>10): 
     return -np.inf 
    return 0 

.