隱馬爾可夫模型中的Viterbi算法

1、隱馬爾可夫模型中的Viterbi算法先用一句話來(lái)簡(jiǎn)單描述一下：給出一個(gè)觀測(cè)序列o1,o2,o3 .，我們希望找到觀測(cè)序列背后的隱藏狀態(tài)序列s1, s2, s3, .；Viterbi以它的發(fā)明者名字命名，正是這樣一種由動(dòng)態(tài)規(guī)劃的方法來(lái)尋找出現(xiàn)概率最大的隱藏狀態(tài)序列（被稱為Viterbi路徑）的算法。這里需要抄一點(diǎn)有關(guān)隱馬可夫序列（HMM，Hidden Markov Model）的書(shū)頁(yè)來(lái)解釋一下觀測(cè)序列和隱藏狀態(tài)序列。首先從最簡(jiǎn)單的離散Markov過(guò)程入手，我們知道，Markov隨機(jī)過(guò)程具有如下的性質(zhì)：在任意時(shí)刻，從當(dāng)前狀態(tài)轉(zhuǎn)移到下一個(gè)狀態(tài)的概率與當(dāng)前狀態(tài)之前的那些狀態(tài)沒(méi)有關(guān)系。所以，我們可以用一

2、個(gè)狀態(tài)轉(zhuǎn)移概率矩陣來(lái)描述它。假設(shè)我們有n個(gè)離散狀態(tài)S1, S2,Sn，我們可以構(gòu)造一個(gè)矩陣A，矩陣中的元素aij表示從當(dāng)前狀態(tài)Si下一時(shí)刻遷移到Sj狀態(tài)的概率。但是在很多情況下，Markov模型中的狀態(tài)是我們觀察不到的。例如，容器與彩球的模型：有若干個(gè)容器，每個(gè)容器中按已知比例放入各色的彩球（這樣，選擇了容器后，我們可以用概率來(lái)預(yù)測(cè)取出各種彩球的可能性）；我們做這樣的實(shí)驗(yàn)，實(shí)驗(yàn)者從容器中取彩球先選擇一個(gè)容器，再?gòu)闹凶コ瞿骋粋€(gè)球，只給觀察者看球的顏色；這樣，每次取取出的球的顏色是可以觀測(cè)到的，即o1, o2,，但是每次選擇哪個(gè)容器是不暴露給觀察者的，容器的序列就組成了隱藏狀態(tài)序列S1, S2,S

3、n。這是一個(gè)典型的可以用HMM描述的實(shí)驗(yàn)。HMM有幾個(gè)重要的任務(wù)，其中之一就是期望通過(guò)觀察序列來(lái)猜測(cè)背后最有可能的隱藏序列。在上面的例子中，就是找到我們?cè)趯?shí)驗(yàn)中最有可能選擇到的容器序列。Viterbi正是用來(lái)解決這個(gè)問(wèn)題的算法。HMM另外兩個(gè)任務(wù)是：a) 給定一個(gè)HMM，計(jì)算一個(gè)觀測(cè)序列出現(xiàn)的可能性；b)已知一個(gè)觀測(cè)序列，HMM參數(shù)不定，如何優(yōu)化這些參數(shù)使得觀測(cè)序列的出現(xiàn)概率最大。解決前一個(gè)問(wèn)題可以用與Viberbi結(jié)構(gòu)非常類似的Forward算法來(lái)解決（實(shí)際上在下面合二為一），而后者可以用Baum-Welch/EM算法來(lái)迭代逼近。從Wiki上抄一個(gè)例子來(lái)說(shuō)明Viterbi算法。假設(shè)你有一個(gè)朋

4、友在外地，每天你可以通過(guò)電話來(lái)了解他每天的活動(dòng)。他每天只會(huì)做三種活動(dòng)之一Walk, Shop, Clean。你的朋友從事哪一種活動(dòng)的概率與當(dāng)?shù)氐臍夂蛴嘘P(guān)，這里，我們只考慮兩種天氣Rainy, Sunny。我們知道，天氣與運(yùn)動(dòng)的關(guān)系如下：例如，在下雨天出去散步的可能性是0.1。而天氣之前互相轉(zhuǎn)換的關(guān)系如下，（從行到列）例如，從今天是晴天而明天就開(kāi)始下雨的可能性是0.4 。同時(shí)為了求解問(wèn)題我們假設(shè)初始情況：通話開(kāi)始的第一天的天氣有0.6的概率是Rainy，有0.4概率是Sunny。OK，現(xiàn)在的問(wèn)題是，如果連續(xù)三天，你發(fā)現(xiàn)你的朋友的活動(dòng)是：Walk->Shop->Clean；那么，如何判

5、斷你朋友那里這幾天的天氣是怎樣的？解決這個(gè)問(wèn)題的python偽代碼如下： def forward_viterbi(y, X, sp, tp, ep):T = for state in X:# prob. V. path V. prob.Tstate = (spstate, state, spstate)for output in y:U = for next_state in X:total = 0argmax = Nonevalmax = 0for source_state in X:(prob, v_path, v_prob) = Tsource_statep = epsource_sta

6、teoutput * tpsource_statenext_stateprob *= pv_prob *= ptotal += probif v_prob > valmax:argmax = v_path + next_statevalmax = v_probUnext_state = (total, argmax, valmax)T = U# apply sum/max to the final states:total = 0argmax = Nonevalmax = 0for state in X:(prob, v_path, v_prob) = Tstatetotal += pr

7、obif v_prob > valmax:argmax = v_pathvalmax = v_probreturn (total, argmax, valmax)幾點(diǎn)說(shuō)明：1、算法對(duì)于每一個(gè)狀態(tài)要記錄一個(gè)三元組：(prob, v_path, v_prob)，其中，prob是從開(kāi)始狀態(tài)到當(dāng)前狀態(tài)所有路徑（不僅僅 是最有可能的viterbi路徑）的概率加在一起的結(jié)果（作為算法附產(chǎn)品，它可以輸出一個(gè)觀察序列在給定HMM下總的出現(xiàn)概率，即forward算法的輸出），v_path是從開(kāi)始狀態(tài)一直到當(dāng)前狀態(tài)的viterbi路徑，v_prob則是該路徑的概率。2、算法開(kāi)始，初始化T （T是一個(gè)Map，

8、將每一種可能狀態(tài)映射到上面所說(shuō)的三元組上）3、三重循環(huán)，對(duì)每個(gè)一活動(dòng)y，考慮下一步每一個(gè)可能的狀態(tài)next_state，并重新計(jì)算若從T中的當(dāng)前狀態(tài)state躍遷到next_state概率會(huì)有怎樣的變化。躍遷主要考慮天氣轉(zhuǎn)移（tpsource_statenext_state）與該天氣下從事某種活動(dòng)（epsource_stateoutput）的聯(lián)合概率。所有下一步狀態(tài)考慮完后，要從T中找出最優(yōu)的選擇viterbi路徑即概率最大的viterbi路徑，即上面更新Map U的代碼Unext_state = (total, argmax, valmax)。4、算法最后還要對(duì)T中的各種情況總結(jié)，對(duì)tota

9、l求和，選擇其中一條作為最優(yōu)的viterbi路徑。5、算法輸出四個(gè)天氣狀態(tài)，這是因?yàn)?，?jì)算第三天的概率時(shí)，要考慮天氣轉(zhuǎn)變到下一天的情況。6、通過(guò)程序的輸出可以幫助理解這一過(guò)程：observation=Walknext_state=Sunnystate=Sunnyp=0.36triple=(0.144,Sunny->,0.144)state=Rainyp=0.03triple=(0.018,Rainy->,0.018)Update USunny=(0.162,Sunny->Sunny->,0.144)next_state=Rainystate=Sunnyp=0.24tri

10、ple=(0.096,Sunny->,0.096)state=Rainyp=0.07triple=(0.042,Rainy->,0.042)Update URainy=(0.138,Sunny->Rainy->,0.096)observation=Shopnext_state=Sunnystate=Sunnyp=0.18triple=(0.02916,Sunny->Sunny->,0.02592)state=Rainyp=0.12triple=(0.01656,Sunny->Rainy->,0.01152)Update USunny=(0.045

11、72,Sunny->Sunny->Sunny->,0.02592)next_state=Rainystate=Sunnyp=0.12triple=(0.01944,Sunny->Sunny->,0.01728)state=Rainyp=0.28triple=(0.03864,Sunny->Rainy->,0.02688)Update URainy=(0.05808,Sunny->Rainy->Rainy->,0.02688)observation=Cleannext_state=Sunnystate=Sunnyp=0.06triple

12、=(0.0027432,Sunny->Sunny->Sunny->,0.0015552)state=Rainyp=0.15triple=(0.008712,Sunny->Rainy->Rainy->,0.004032)Update USunny=(0.0114552,Sunny->Rainy->Rainy->Sunny->,0.004032)next_state=Rainystate=Sunnyp=0.04triple=(0.0018288,Sunny->Sunny->Sunny->,0.0010368)state=

13、Rainyp=0.35triple=(0.020328,Sunny->Rainy->Rainy->,0.009408)Update URainy=(0.0221568,Sunny->Rainy->Rainy->Rainy->,0.009408)final triple=(0.033612,Sunny->Rainy->Rainy->Rainy->,0.009408)所以，最終的結(jié)果是，朋友那邊這幾天最可能的天氣情況是Sunny->Rainy->Rainy->Rainy，它有0.009408的概率出現(xiàn)。而我們算法的

14、另一個(gè)附帶的結(jié)論是，我們所觀察到的朋友這幾天的活動(dòng)序列：Walk->Shop->Clean在我們的隱馬可夫模型之下出現(xiàn)的總概率是0.033612。附：c+主要代碼片斷void forward_viterbi(const vector<string> & ob, viterbi_triple_t & vtriple)./aliasmap<string, double>& sp = start_prob;map<string, map<string, double> > & tp = transition_

15、prob;map<string, map<string, double> > & ep = emission_prob;/ initializationInitParameters();map<string, viterbi_triple_t> T;for (vector<string>:iterator it = states.begin();it != states.end();+it).viterbi_triple_t foo;b = sp*it;foo.vpath.push_back(*it);foo.vprob =

16、 sp*it;T*it = foo;map<string, viterbi_triple_t> U;double total = 0;vector<string> argmax;double valmax = 0;double p = 0;for (vector<string>:const_iterator itob = ob.begin(); itob != ob.end(); +itob).cout << "observation=" << *itob << endl;U.clear();for (

17、vector<string>:iterator itNextState = states.begin();itNextState != states.end();+itNextState).cout << " next_state=" << *itNextState << endl;total = 0;argmax.clear();valmax = 0;for (vector<string>:iterator itSrcState = states.begin();itSrcState != states.end(

18、);+itSrcState).cout << " state=" << *itSrcState << endl;viterbi_triple_t foo = T*itSrcState;p = ep*itSrcState*itob * tp*itSrcState*itNextState;cout << " p=" << p << endl;b *= p;foo.vprob *= p;cout << " triple=" << foo << endl;total += b;if (foo.vprob > valmax).foo.vpath.push_back(*itNextState);argmax = foo.vpath;valma