PyMC3のシンプルな隠れマルコフモデルがTheanoエラーを投げる

Question

私はPyMC3、Theano、numpyを使い慣れていません。スタンのマニュアルで最初の「隠された」マルコフモデルを複製しようとしていました - 実際に状態が観察されたもの。しかし、私はTheanoと間違い、numpy、そしておそらく私にとって少し不思議そうなPyMC3ディストリビューションの後ろで何が起こっているのかを考えています。モデルのための私のコード以下である：

import pandas as pd 
dat_hmm = pd.read_csv('hmmVals.csv') 
emission=dat_hmm.emission.values 
state=dat_hmm.state.values 

from pymc3 import Model, Dirichlet, Categorical 
import numpy as np 

basic_model = Model() 

with basic_model: 
    #Model constants: 
    #num unique hidden states, num unique emissions, num instances 
    K=3; V=9; T=10 
    alpha=np.ones(K); beta=np.ones(V) 
    # Priors for unknown model parameters 
    theta = np.empty(K, dtype=object) #theta=transmission 
    phi = np.empty(K, dtype=object) #phi=emission 
    #observed emission, state: 
    w=np.empty(T, dtype=object); z=np.empty(T, dtype=object); 
    for k in range(K): 
     theta[k]=Dirichlet('theta'+str(k), alpha) 
     phi[k]=Dirichlet('phi'+str(k), beta) 
    # Likelihood (sampling distribution) of observations 
    for t in range(T): 
     w[t]=Categorical('w'+str(t),theta[state[t]], shape=1, observed=emission[t]) 
    for t in range(2, T): 
     z[t]=Categorical('z'+str(t),phi[state[t-1]], shape=1, observed=state[t]) 

線「W [T] =カテゴリ（ 'W' + STR（T）=放射観察シータ[状態[T]、形状= 1、[ t]） "はエラーを生成しますが、w0を埋め込むt = 0ではなく、t = 1でインデックス外のエラーを生成します。 state [1]、theta [state [t]]、およびemission [t]がすべて存在するため、コードライン自体にインデックスがありません。エラーメッセージは以下のとおりです。

Traceback (most recent call last): 
    File "/usr/local/lib/python3.4/dist-packages/pymc3/distributions/distribution.py", line 25, in __new__ 
    return model.Var(name, dist, data) 
    File "/usr/local/lib/python3.4/dist-packages/pymc3/model.py", line 306, in Var 
    var = ObservedRV(name=name, data=data, distribution=dist, model=self) 
    File "/usr/local/lib/python3.4/dist-packages/pymc3/model.py", line 581, in __init__ 
    self.logp_elemwiset = distribution.logp(data) 
    File "/usr/local/lib/python3.4/dist-packages/pymc3/distributions/discrete.py", line 400, in logp 
    a = tt.log(p[value]) 
    File "/usr/local/lib/python3.4/dist-packages/theano/tensor/var.py", line 532, in __getitem__ 
    lambda entry: isinstance(entry, Variable))) 
    File "/usr/local/lib/python3.4/dist-packages/theano/gof/op.py", line 668, in __call__ 
    required = thunk() 
    File "/usr/local/lib/python3.4/dist-packages/theano/gof/op.py", line 883, in rval 
    fill_storage() 
    File "/usr/local/lib/python3.4/dist-packages/theano/gof/cc.py", line 1707, in __call__ 
    reraise(exc_type, exc_value, exc_trace) 
    File "/usr/local/lib/python3.4/dist-packages/six.py", line 686, in reraise 
    raise value 
IndexError: index out of bounds 

私はPyMC3分布にnumpyのオブジェクトを貼り付けるか、別の分布をパラメータ化しようとすることの結果を使用しての知恵を知りませんが、私は、ウェブ上やや似たコードを見てきました、最後の部分をマイナスします。 PyMC3でこのような隠れマルコフモデルをコーディングする良い方法はないでしょうか？

Answer 1

私は上記のエラーを修正する方法を発見しました。次のコードは動作します。エラーはなく、少なくともメトロポリスで正しいパラメータ推定値を得ることができます。

私は2つの間違いを犯しましたが、私はテアノで何か複雑なことが起こると予想していたので、彼らはとてもシンプルだったとは思いませんでした。 1つは、私のデータがStan用に設定されていて、0ではなく1で始まるように索引付けされているということです。Pythonはすべてを0にインデックスします。もう一つの誤差は、フィー・マトリックスの個々のエミッションを計算するために伝達行列であるシータを使用し、その逆も同様でした。シータは排出量が足りなかった。

私が今思ったのは、NUTSサンプラーがMAP推定値を与えているにもかかわらず、非正定のスケーリングをしていることを伝え続けている理由でした。メトロポリスは機能しますが、これらの300回の観測と1000回のサンプルでは約11分です。もう一つの謎は、PyMC3がサンプルの計算に数秒しかかかりませんでした。

import pandas as pd 

dat_hmm = pd.read_csv('hmmVals.csv') 

emission=dat_hmm.emission.values 
state=dat_hmm.state.values 

from pymc3 import Model, Dirichlet, Categorical 
import numpy as np 

basic_model = Model() 

with basic_model: 
    #Model constants: 
    K=3; V=9; T=300 #num unique hidden states, num unique emissions, num instances 
    alpha=np.ones(K); beta=np.ones(V) 
    # Priors for unknown model parameters 
    theta = np.empty(K, dtype=object) #theta=transmission 
    phi = np.empty(K, dtype=object) #phi=emission 
    w=np.empty(T, dtype=object); z=np.empty(T, dtype=object); #observed emission, state 
    for k in range(K): 
     theta[k]=Dirichlet('theta'+str(k), alpha) 
     phi[k]=Dirichlet('phi'+str(k), beta) 
    #Likelihood (sampling distribution) of observationss 
    for t in range(2, T): 
     z[t]=Categorical('z'+str(t),theta[state[t-1]], shape=1, observed=state[t]) 
    for t in range(T): 
     w[t]=Categorical('w'+str(t),phi[state[t]], shape=1, observed=emission[t]) 

Answer 2

また、pymc3でHMMを実装しようとしましたが、同様の問題が発生しました。私はちょうどベクトル化された方法で2つのレベルのHMMを実装する方法を見つけました（私のモデルは隠されていませんが、隠れた部分は簡単に追加できます - 私は状態変数の説明をベクトル化しました）。私はこれが最も効率的な方法かどうかはわかりませんが、このコードを状態を定義するための単純なforループに対してテストしました。以下のコードは1000データポイントで1分未満で実行されますが、forループは数時間かかりました。ここで

はコードです：

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

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. 

    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) 
     length = x.shape[0] 
     P1T = tt.tile(P1,(length-1,1)).T 
     P2T = tt.tile(P2,(length-1,1)).T 

     P = tt.switch(x[:-1],P1T,P2T).T 

     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) 

このクラスは、HMMの状態を作成します。それを呼び出すには、次の操作を行うことができます

theta = np.ones(2) # prior for probabilities 
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=theta) 
    P2 = pm.Dirichlet('P2', a=theta) 

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

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

    start = pm.find_MAP() 

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

だけで古いコードがどのように見えるかを示すために：このソリューションは、多くの人々が実行されているように見える（参考になりました場合は

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)) 

    state = pm.Categorical('state0',PA, observed=data[0]) 
    for i in range(1,N_chain): 
     state = pm.Categorical('state'+str(i), tt.switch(data[i-1],P1,P2), observed=data[i]) 

    start = pm.find_MAP() 

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