pymc3 몬테카를로 스테퍼를 사용자 정의 카테고리 배포에 사용할 수 있습니까?

Question

나는 pymc3에서 숨겨진 Markov-Chains을 구현하려고합니다. 숨겨진 상태를 구현하는 데 꽤 오래되었습니다. 아래에서는 간단한 2- 상태 마르코프 체인을 보여줍니다.

import numpy as np 
import pymc3 as pm 
import theano.tensor as tt 

# Markov chain sample with 2 states that was created 
# to have prob 0->1 = 0.1 and prob 1->0 = 0.3 
sample = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 
    0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 
    1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 
    0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0], 
dtype=np.uint8) 

이제 상태를 설명하는 클래스를 정의하고 있습니다. 입력으로 상태 0에서 상태 1로 이동하는 확률 P1과 1-> 0에서 이동하는 P2를 알아야합니다. 나는 또한 내가 theano 구글 그룹에서 배운 고급 인덱싱 닌자 트릭의 매우 자랑스럽게 생각합니다 0

class HMMStates(pm.Discrete): 
    """ 
    Hidden Markov Model States 
    Parameters 
    ---------- 
    P1 : tensor 
     probability to remain in state 1 
    P2 : tensor 
     probability to move from state 2 to state 1 

    """ 

    def __init__(self, PA=None, P1=None, P2=None, 
       *args, **kwargs): 
     super(HMMStates, self).__init__(*args, **kwargs) 
     self.PA = PA 
     self.P1 = P1 
     self.P2 = P2 
     self.mean = 0. 
     self.mode = tt.cast(0,dtype='int64') 

    def logp(self, x): 
     PA = self.PA 
     P1 = self.P1 
     P2 = self.P2 

     # now we need to create an array with probabilities 
     # so that for x=A: PA=P1, PB=(1-P1) 
     # and for x=B: PA=P2, PB=(1-P2) 

     choice = tt.stack((P1,P2)) 
     P = choice[x[:-1]] 

     x_i = x[1:] 

     ou_like = pm.Categorical.dist(P).logp(x_i) 
     return pm.Categorical.dist(PA).logp(x[0]) + tt.sum(ou_like) 

되는 최초의 상태에 대한 확률 PA를 알아야합니다. 또한 tt.switch를 사용하여 구현할 수도 있습니다. 내가 너무 확신하지 못했던 것은 자기 모드 다. 방금 testvalue 오류를 피하기 위해 0을주었습니다. 다음은 작동 여부를 테스트하는 모델에서 클래스를 사용하는 방법입니다. 이 경우 상태는 숨겨지지 않고 관찰됩니다.

with pm.Model() as model: 
    # 2 state model 
    # P1 is probablility to stay in state 1 
    # P2 is probability to move from state 2 to state 1 
    P1 = pm.Dirichlet('P1', a=np.ones(2)) 
    P2 = pm.Dirichlet('P2', a=np.ones(2)) 

    PA = pm.Deterministic('PA',P2/(P2+1-P1)) 

    states = HMMStates('states',PA,P1,P2, observed=sample) 

    start = pm.find_MAP() 

    trace = pm.sample(5000, start=start) 

출력이 데이터를 멋지게 재생합니다. 다음 모델에서는 문제를 보여 드리겠습니다. 여기에서는 상태를 직접 관찰하지는 않지만 일부 가우시안 노이즈가 추가 된 상태 (따라서 숨겨진 상태)를 관찰합니다. Metropolis stepper를 사용하여 모델을 실행하면 인덱스 오류로 인해 충돌이 발생합니다.이 오류는 대도시 배포판에서 Metropolis stepper를 사용하는 것과 관련된 문제로 거슬러 올라갑니다. 불행히도, 내 수업에 적용 할 수있는 단 Stepper는 CategoricalGibbsMetropolis 스테퍼입니다. 그러나 CategoricalGibbsMetropolis 스테퍼는 명시 적으로 Categorial Distribution이 아니기 때문에 내 수업에 참여하는 것을 거부합니다.

gauss_sample = sample*1.0 + 0.1*np.random.randn(len(sample)) 
from scipy import optimize 
with pm.Model() as model2: 
    # 2 state model 
    # P1 is probablility to stay in state 1 
    # P2 is probability to move from state 2 to state 1 
    P1 = pm.Dirichlet('P1', a=np.ones(2)) 
    P2 = pm.Dirichlet('P2', a=np.ones(2)) 

    S = pm.InverseGamma('S',alpha=2.1, beta=1.1) 

    PA = pm.Deterministic('PA',P2/(P2+1-P1)) 

    states = HMMStates('states',PA,P1,P2, shape=len(gauss_sample)) 

    emission = pm.Normal('emission', 
         mu=tt.cast(states,dtype='float64'), 
         sd=S, 
         observed = gauss_sample) 

    start2 = pm.find_MAP(fmin=optimize.fmin_powell) 
    step1 = pm.Metropolis(vars=[P1, P2, S, PA, emission]) 
    step2 = pm.ElemwiseCategorical(vars=[states], values=[0,1]) 
    trace2 = pm.sample(10000, start=start, step=[step1,step2]) 

ElemwiseCategorical은 실행되지만 내 상태에 대해 올바른 값을 할당하지 않습니다. 상태는 모두 0 또는 모두 1입니다.

어떻게 ElemwiseCategorial에 1과 0의 상태 벡터를 할당할지 또는 CategorialGibbsMetropolis에서 내 배포를 범주형으로 인식하도록하려면 어떻게해야합니까? 이는 커스텀 배포판에 공통적 인 문제 여야합니다.

Answer 1

내 질문에 아무 에게서도 소식을 듣지 못했기 때문에 직접 답변했습니다. 내가 사용한 트릭은 Chris Fonnesbeck이 pymc3 github에서 제안한 것으로 제안되었습니다. 그는 pm.Categorical을 하위 클래스로 제안했습니다.

class HMMStates(pm.Categorical): 
    """ 
    Hidden Markov Model States 
    Parameters 
    ---------- 
    P1 : tensor 
     probability to remain in state 1 
    P2 : tensor 
     probability to move from state 2 to state 1 

    """ 

    def __init__(self, PA=None, P1=None, P2=None, 
       *args, **kwargs): 
     super(pm.Categorical, self).__init__(*args, **kwargs) 
     self.PA = PA 
     self.P1 = P1 
     self.P2 = P2 
     self.k = 2 # two state model 
     self.mean = 0. 
     self.mode = tt.cast(0,dtype='int64') 

    def logp(self, x): 
     PA = self.PA 
     P1 = self.P1 
     P2 = self.P2 

     # now we need to create an array with probabilities 
     # so that for x=A: PA=P1, PB=(1-P1) 
     # and for x=B: PA=P2, PB=(1-P2) 
     PT = tt.stack((P1,P2)) 

     P = PT[x[:-1]] 

     x_i = x[1:] 

     ou_like = pm.Categorical.dist(P, shape=(N_chain-1,2)).logp(x_i) 
     return pm.Categorical.dist(PA).logp(x[0]) + tt.sum(ou_like) 

내 HMMStates 정말 pm.Categorical 슈퍼 초기화를 호출 할 수 없습니다, 그래서 pm.Discrete입니다 pm.Categorical의 슈퍼 클래스를 호출하고 있습니다. 그 트릭은 BinaryGibbsMetropolis와 CategoricalGibbsMetropolis의 테스트를 통과시킵니다.

2 상태 및 다중 상태 HMM을 구현하는 데 관심이 있으시면 here을 모두 구현하겠습니다.