PyMC3で離散HMMを構築する際の問題

Question

PyMC3でおもちゃのディスクリートHMMモデルを使用していますが、いくつかの問題があります。

このような私の最初のコードを見て：

# Known transition and emission 
with pymc3.Model() as hmm1: 
    T = tt.as_tensor_variable(Tr) 
    E = tt.as_tensor_variable(Er) 
    # State models 
    p0 = np.ones(num_states)/num_states 
    # No shape, so each state is a scalar tensor 
    states = [pymc3.Categorical('s0', p=p0)] 
    emissions = [pymc3.Categorical('z0', p=E[:,states[0]], observed=Zr[:,0])] 
    for i in range(1, num_times): 
     states.append(pymc3.Categorical('s{0}'.format(i), p=T[:,states[i-1]])) 
     emissions.append(pymc3.Categorical('z{0}'.format(i), p=E[:,states[i]], observed=Zr[:,i])) 

ここ

TrとErは本当の推移と排出量行列です。

私の問題は、次のとおりです。

1. このモデルは、状態の完全な値を探索するようではありません、それは各状態のための単一の値にとどまる（ノートを参照してください）。

2. statesとemissionsをより厳密に定義する方法が見つかりませんでした。 shape=...を使用してください。

3. さらに私は、遷移または放出マトリックス不明を説明するためにモデルを拡張するとき、私は、インデックスの問題に走った、と私は次のコードのように、theano.tensor.clipを使用するように強制しています：

   # Unknown transitions and emissions 
   with pymc3.Model() as hmm3: 
       # Transition "matrix" 
       a_t = np.ones((num_states,num_states)) 
       T = pymc3.Dirichlet('T', a=a_t, shape=(num_states,num_states)) 
       # Emission "matrix" 
       a_e = np.ones((num_emissions, num_states)) 
       E = pymc3.Dirichlet('E', a=a_e, shape=(num_emissions, num_states)) 
       # State models 
       p0 = np.ones(num_states)/num_states 
       # No shape, so each state is a scalar tensor 
       states = [pymc3.Categorical('s0', p=p0)] 
       clp_state = tt.clip(states[0], 0, num_states-1) 
       emissions = [pymc3.Categorical('z0', p=E[:,clp_state], observed=Zr[0])] 
       for i in range(1, num_times): 
        clp_prev_state = tt.clip(states[i-1], 0, num_states-1) 
        states.append(pymc3.Categorical('s{0}'.format(i), p=T[:,clp_prev_state])) 
        clp_state = tt.clip(states[i], 0, num_states-1) 
        emissions.append(pymc3.Categorical('z{0}'.format(i), p=E[:,clp_state], observed=Zr[i])) 
   

notebook with a complete codeを参照してください。 （あるいは、それ以上のtheanoウェイ）

Answer 1

よりpymc3の方法以下のように隠れマルコフ状態モデルが定義する：

class HMMStatesN(pm.Categorical): 
    """ 
    Hidden Markov Model States 
    Parameters 
    ---------- 
    P : tensor 
     transition probability 
     shape = (N_states,N_states) 

    PA : tensor 
     equilibrium probabilities 
     shape = (N_states) 

    """ 

    def __init__(self, PA=None, P=None, 
       *args, **kwargs): 
     super(pm.Categorical, self).__init__(*args, **kwargs) 
     self.P = P 
     self.PA = PA 
     self.k = N_states 
     self.mode = tt.cast(0,dtype='int64') 

    def logp(self, x): 
     P = self.P 
     PA = self.PA 

     # calculate equilibrium 

     # 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) 
     #P = tt.switch(x[:-1],P1,P2) 

     PS = P[x[:-1]] 

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

Iも同様のガウス発光クラスを作成した：

class HMMGaussianEmmisionsN(pm.Continuous): 
    """ 
    Hidden Markov Model Gaussian Emmisions 
    Parameters 
    ---------- 
    A : tensor 
     prior for Gaussian emmision mu 
     shape = (2,N_states) 

    S : tensor 
     prior for Gaussian emmision width 
     shape = (2,N_states) 

    states : tensor 
     equilibrium probabilities 
     shape = (N_states) 

    """ 

    def __init__(self, A=None, S=None, states=None, 
       *args, **kwargs): 
     super(HMMGaussianEmmisionsN, self).__init__(*args, **kwargs) 
     self.A = A 
     self.S = S 
     self.states = states 
     self.mean = 0. 

    def logp(self, x): 
     A = self.A 
     S = self.S 
     states = self.states 

     AS = A[states]   
     SS = S[states] 

     ou_like = pm.Normal.dist(mu=AS,sd=SS).logp(x) 
     return tt.sum(ou_like) 

以下の方法で呼び出すことができ

：

from scipy import optimize 
with pm.Model() as model: 
    # N_states state model 

    P = pm.Dirichlet('P', a=np.ones((N_states,N_states)), shape=(N_states,N_states)) 
    A = pm.InverseGamma('A',alpha=alphaA, beta=betaA, shape=(N_states)) 
    S = pm.InverseGamma('S', alpha=alphaS, beta=betaS, shape=(N_states)) 

    AA = tt.dmatrix('AA') 
    AA = tt.eye(N_states) - P + tt.ones(shape=(N_states,N_states)) 
    PA = pm.Deterministic('PA',sla.solve(AA.T,tt.ones(shape=(N_states)))) 

    states = HMMStatesN('states',PA=PA, P=P, shape=(len(measurement))) 

    emmision = HMMGaussianEmmisionsN('emmision', A=A, S=S, states=states, observed=measurement) 

    start = pm.find_MAP(fmin=optimize.fmin_powell) 
    step1 = pm.Metropolis(vars=[P,A,S,PA,emmision]) 
    step2 = pm.CategoricalGibbsMetropolis(vars=[states]) 
    trace = pm.sample(10000, start=start, step=[step1,step2]) 

は、私はあなたのモデルができると思うだろう同様の方法でコード化する。 HMMstatesのCategoricalGibbsMetropolisステッパー（詳しくはquestionを参照）を使用することが重要です。より詳細なコード例hereを参照してください。