ステルスマルコフモデルにおけるViterbiアルゴリズム

この文章はViterbiアルゴリズムを簡単に説明します.1年前に名前を聞いたことがありますが、2週間前に毛皮を研究するのに少し時間がかかりました.ここで簡単な検討をします.まず一言で簡単に説明します.観測シーケンスo 1,o 2,o 3...、観測シーケンスの背後にある隠れた状態シーケンスs 1,s 2,s 3,...を見つけたいと思っています.Viterbiはその発明者の名前で命名され、動的に計画された方法によって出現確率が最も大きい隠蔽状態シーケンス(Viterbiパスと呼ばれる)を探すアルゴリズムである.
ここでは,観測シーケンスと非表示状態シーケンスを説明するために,ステルスマルコフシーケンス(HMM,Hidden Markov Model)に関するページを写す必要がある.
まず最も単純な離散Markov過程から出発し,Markovランダム過程は任意の時点で現在の状態から次の状態に移行する確率が現在の状態以前の状態とは関係ないという性質を持つことを知った.従って,状態遷移確率行列を用いてそれを記述することができる.n個の離散状態S 1,S 2,…Snがあると仮定すると,現在の状態Siから次の時点でSj状態に遷移する確率を表す行列Aを構築することができる.
しかし,多くの場合,Markovモデルにおける状態は我々には観察できない.例えば、容器とカラーボールのモデル:いくつかの容器があり、各容器に既知の割合で各色のカラーボールを入れる(このように、容器を選択した後、確率で各種カラーボールを取り出す可能性を予測することができる).私たちはこのような実験をして、実験者は容器の中からカラーボールを取ります--まず1つの容器を選んで、それから中からあるボールをつかんで、観察者にボールの色だけを見せます;このように、取り出される球の色は毎回観測できる、すなわちo 1,o 2,…,であるが、どの容器が観察者に露出しないかを選択するたびに、容器の配列が隠蔽状態配列S 1,S 2,…Snとなる.これはHMMで記述できる典型的な実験である.
HMMにはいくつかの重要なタスクがあり、その1つは、シーケンスを観察することによって、背後で最も可能性のある隠蔽シーケンスを推測することを期待することである.上記の例では,実験で最も選択可能な容器配列を見つけた.Viterbiはこの問題を解決するためのアルゴリズムです.HMMの他の2つの任務は:a)1つのHMMを与えて、1つの観測シーケンスが現れる可能性を計算します;b)観測シーケンスが知られており、HMMパラメータが不定であり、これらのパラメータをどのように最適化して観測シーケンスの出現確率を最大にするか.前の問題を解決するにはViberbi構造と非常に類似したForwardアルゴリズムを用いて解決することができ(実際には以下で2つを1つにする)、後者はBaum−Welch/EMアルゴリズムを用いて反復近似することができる.
Wikiから例を写してViterbiアルゴリズムを説明します.
もしあなたが地方に友达がいるとしたら、毎日電話で彼の毎日の活動を知ることができます.彼は毎日3つの活動の1つであるWalk,Shop,Cleanしかできない.あなたの友達がどの活動に従事する確率は現地の気候と関係があります.ここでは、Rainy、Sunnyの2つの天気しか考えていません.天気と運動の関係は以下の通りです.
 
Rainy
Sunny
Walk
0.1
0.6
Shop
0.4
0.3
Clean
0.5
0.1
例えば、雨の日に散歩に出かける可能性は0.1です.
天気の前に互いに転換する関係は以下の通りです(行から列まで)
 
Rainy
Sunny
Rainy
0.7
0.3
Sunny
0.4
0.6
例えば、今日が晴れて明日から雨が降る可能性は0.4です.
同時に問題を解くために,通話開始初日の天気が0.6の確率でRainy,0.4の確率でSunnyであると仮定した.OK、今の問題は、もし3日連続で、あなたの友达の活動を発見したら:Walk->Shop->Clean;では、あなたの友达のところのここ数日の天気はどうなっていますか.
この問題を解決するpython擬似コードは以下の通りです.

def forward_viterbi(y, X, sp, tp, ep):

T = {}

for state in X:

##          prob.      V. path  V. prob.

T[state] = (sp[state], [state], sp[state])

for output in y:

U = {}

for next_state in X:

total = 0

argmax = None

valmax = 0

for source_state in X:

(prob, v_path, v_prob) = T[source_state]

p = ep[source_state][output] * tp[source_state][next_state]

prob *= p

v_prob *= p

total += prob

if v_prob > valmax:

argmax = v_path + [next_state]

valmax = v_prob

U[next_state] = (total, argmax, valmax)

T = U

## apply sum/max to the final states:

total = 0

argmax = None

valmax = 0

for state in X:

(prob, v_path, v_prob) = T[state]

total += prob

if v_prob > valmax:

argmax = v_path

valmax = v_prob

return (total, argmax, valmax)

 
  
几点说明：
 
  
 
   
 算法对于每一个状态要记录一个三元组：(prob, v_path, v_prob)，其中，prob是从开始状态到当前状态所有路径（不仅仅 是最有可能的viterbi路径）的概率加在一起的结果（作为算法附产品，它可以输出一个观察序列在给定HMM下总的出现概率，即forward算法的输出），v_path是从开始状态一直到当前状态的viterbi路径，v_prob则是该路径的概率。
 
   
 算法开始，初始化T （T是一个Map，将每一种可能状态映射到上面所说的三元组上）
 
   
 三重循环，对每个一活动y，考虑下一步每一个可能的状态next_state，并重新计算若从T中的当前状态state跃迁到next_state概率会有怎样的变化。跃迁主要考虑天气转移（tp[source_state][next_state]）与该天气下从事某种活动（ep[source_state][output]）的联合概率。所有下一步状态考虑完后，要从T中找出最优的选择viterbi路径——即概率最大的viterbi路径，即上面更新Map U的代码U[next_state] = (total, argmax, valmax)。
 
   
 算法最后还要对T中的各种情况总结，对total求和，选择其中一条作为最优的viterbi路径。
 
   
 算法输出四个天气状态，这是因为，计算第三天的概率时，要考虑天气转变到下一天的情况。
 
   
 通过程序的输出可以帮助理解这一过程：
 
  
 
  
observation=Walk
  
next_state=Sunny
  
state=Sunny
  
p=0.36
  
triple=(0.144,Sunny->,0.144)
  
state=Rainy
  
p=0.03
  
triple=(0.018,Rainy->,0.018)
  
Update U[Sunny]=(0.162,Sunny->Sunny->,0.144)
  
next_state=Rainy
  
state=Sunny
  
p=0.24
  
triple=(0.096,Sunny->,0.096)
  
state=Rainy
  
p=0.07
  
 
  
triple=(0.042,Rainy->,0.042)
  
Update U[Rainy]=(0.138,Sunny->Rainy->,0.096)
  
observation=Shop
  
next_state=Sunny
  
state=Sunny
  
p=0.18
  
triple=(0.02916,Sunny->Sunny->,0.02592)
  
state=Rainy
  
p=0.12
  
triple=(0.01656,Sunny->Rainy->,0.01152)
  
Update U[Sunny]=(0.04572,Sunny->Sunny->Sunny->,0.02592)
  
next_state=Rainy
  
state=Sunny
  
p=0.12
  
triple=(0.01944,Sunny->Sunny->,0.01728)
  
state=Rainy
  
p=0.28
  
triple=(0.03864,Sunny->Rainy->,0.02688)
  
Update U[Rainy]=(0.05808,Sunny->Rainy->Rainy->,0.02688)
  
observation=Clean
  
next_state=Sunny
  
state=Sunny
  
p=0.06
  
triple=(0.0027432,Sunny->Sunny->Sunny->,0.0015552)
  
state=Rainy
  
p=0.15
  
triple=(0.008712,Sunny->Rainy->Rainy->,0.004032)
  
Update U[Sunny]=(0.0114552,Sunny->Rainy->Rainy->Sunny->,0.004032)
  
next_state=Rainy
  
state=Sunny
  
p=0.04
  
triple=(0.0018288,Sunny->Sunny->Sunny->,0.0010368)
  
state=Rainy
  
p=0.35
  
triple=(0.020328,Sunny->Rainy->Rainy->,0.009408)
  
Update U[Rainy]=(0.0221568,Sunny->Rainy->Rainy->Rainy->,0.009408)
  
final triple=(0.033612,Sunny->Rainy->Rainy->Rainy->,0.009408)
  
だから、最終的な結果は、友達のところでここ数日最も可能性のある天気はSunny->Rainy->Rainy->Rainyで、0.009408の確率で現れます.我々のアルゴリズムのもう一つの付随的な結論は,我々が観察した友人のここ数日の活動シーケンス:Walk->Shop->Cleanが我々のステルスマルコフモデルの下で現れる総確率は0.033612である.
参考文献
http://www.ece.ucsb.edu/Faculty/Rabiner/ece259/Reprints/tutorial%20on%20hmm%20and%20applications.pdf   
http://en.wikipedia.org/wiki/Viterbi_algorithm   
http://googlechinablog.com/2006/04/blog-post_17.html   
添付:c++プライマリコード断片
void forward_viterbi(const vector & ob, viterbi_triple_t & vtriple) 
{ 
//alias 
map& sp = start_prob; 
map> & tp = transition_prob; 
map> & ep = emission_prob; 
//initialization 
InitParameters(); 
map T; 
for (vector::iterator it = states.begin(); 
it != states.end(); 
++it) 
{ 
viterbi_triple_t foo; 
foo.prob = sp[*it]; 
foo.vpath.push_back(*it); 
foo.vprob = sp[*it]; 
T[*it] = foo; 
} 
map U; 
double total = 0; 
vector argmax; 
double valmax = 0; 
double p = 0; 
for (vector::const_iterator itob = ob.begin(); itob != ob.end();++itob) 
{ 
cout << “observation=” << *itob << endl; 
U.clear(); 
for (vector::iterator itNextState = states.begin(); 
itNextState != states.end(); 
++itNextState) 
{ 
cout << “\tnext_state=” << *itNextState << endl; 
total = 0; 
argmax.clear(); 
valmax = 0; 
for (vector::iterator itSrcState = states.begin(); 
itSrcState != states.end(); 
++itSrcState) 
{ 
cout << “\t\tstate=” << *itSrcState << endl; 
viterbi_triple_t foo = T[*itSrcState]; 
p = ep[*itSrcState][*itob] * tp[*itSrcState][*itNextState]; 
cout << “\t\t\tp=” << p << endl; 
foo.prob *= p; 
foo.vprob *= p; 
cout << “\t\t\ttriple=” << foo << endl; 
total += foo.prob; 
if (foo.vprob > valmax) 
{ 
foo.vpath.push_back(*itNextState); 
argmax = foo.vpath; 
valmax = foo.vprob; 
} 
} 
U[*itNextState] = viterbi_triple_t(total, argmax, valmax); 
cout << “\tUpdate U["<< *itNextState << "]=” << U[*itNextState] << “” << endl; 
} 
T.swap(U); 
} 
total = 0; 
argmax.clear(); 
valmax = 0; 
for (vector::iterator itState = states.begin(); 
itState != states.end(); 
++itState) 
{ 
viterbi_triple_t foo = T[*itState]; 
total += foo.prob; 
if (foo.vprob > valmax) 
{ 
argmax.swap(foo.vpath); 
valmax = foo.vprob; 
} 
} 
vtriple.prob = total; 
vtriple.vpath = argmax; 
vtriple.vprob = valmax; 
cout << “final triple=” << vtriple << endl; 
}