這篇文章簡單描述一下Viterbi算法——一年之前我聽過它的名字,直到兩周之前才花了一點時間研究了個皮毛,在這里做個簡單檢討。先用一句話來簡單描述一下:給出一個觀測序列o1,o2,o3 …,我們希望找到觀測序列背后的隱藏狀態序列s1, s2, s3, …;Viterbi以它的發明者名字命名,正是這樣一種由動態規劃的方法來尋找出現概率最大的隱藏狀態序列(被稱為Viterbi路徑)的算法。
這里需要抄一點有關隱馬可夫序列(HMM,Hidden Markov Model)的書頁來解釋一下觀測序列和隱藏狀態序列。
首先從最簡單的離散Markov過程入手,我們知道,Markov隨機過程具有如下的性質:在任意時刻,從當前狀態轉移到下一個狀態的概率與當前狀態之前的那些狀態沒有關系。所以,我們可以用一個狀態轉移概率矩陣來描述它。假設我們有n個離散狀態S1, S2,…Sn,我們可以構造一個矩陣A,矩陣中的元素aij表示從當前狀態Si下一時刻遷移到Sj狀態的概率。
但是在很多情況下,Markov模型中的狀態是我們觀察不到的。例如,容器與彩球的模型:有若干個容器,每個容器中按已知比例放入各色的彩球(這樣,選擇了容器后,我們可以用概率來預測取出各種彩球的可能性);我們做這樣的實驗,實驗者從容器中取彩球——先選擇一個容器,再從中抓出某一個球,只給觀察者看球的顏色;這樣,每次取取出的球的顏色是可以觀測到的,即o1, o2,…,但是每次選擇哪個容器是不暴露給觀察者的,容器的序列就組成了隱藏狀態序列S1, S2,…Sn。這是一個典型的可以用HMM描述的實驗。
HMM有幾個重要的任務,其中之一就是期望通過觀察序列來猜測背后最有可能的隱藏序列。在上面的例子中,就是找到我們在實驗中最有可能選擇到的容器序列。Viterbi正是用來解決這個問題的算法。HMM另外兩個任務是:a) 給定一個HMM,計算一個觀測序列出現的可能性;b)已知一個觀測序列,HMM參數不定,如何優化這些參數使得觀測序列的出現概率最大。解決前一個問題可以用與Viberbi結構非常類似的Forward算法來解決(實際上在下面合二為一),而后者可以用Baum-Welch/EM算法來迭代逼近。
從Wiki上抄一個例子來說明Viterbi算法。
假設你有一個朋友在外地,每天你可以通過電話來了解他每天的活動。他每天只會做三種活動之一——Walk, Shop, Clean。你的朋友從事哪一種活動的概率與當地的氣候有關,這里,我們只考慮兩種天氣——Rainy, Sunny。我們知道,天氣與運動的關系如下:
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,現在的問題是,如果連續三天,你發現你的朋友的活動是: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 = 0argmax = Nonevalmax = 0for source_state in X:(prob, v_path, v_prob) = T[source_state]p = ep[source_state][output] * tp[source_state][next_state]prob *= pv_prob *= ptotal += probif v_prob > valmax:argmax = v_path + [next_state]valmax = v_probU[next_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) = T[state]total += probif v_prob > valmax:argmax = v_pathvalmax = v_probreturn (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=Walknext_state=Sunnystate=Sunnyp=0.36triple=(0.144,Sunny->,0.144)state=Rainyp=0.03triple=(0.018,Rainy->,0.018)Update U[Sunny]=(0.162,Sunny->Sunny->,0.144)next_state=Rainystate=Sunnyp=0.24triple=(0.096,Sunny->,0.096)state=Rainyp=0.07triple=(0.042,Rainy->,0.042)Update U[Rainy]=(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 U[Sunny]=(0.04572,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 U[Rainy]=(0.05808,Sunny->Rainy->Rainy->,0.02688)observation=Cleannext_state=Sunnystate=Sunnyp=0.06triple=(0.0027432,Sunny->Sunny->Sunny->,0.0015552)state=Rainyp=0.15triple=(0.008712,Sunny->Rainy->Rainy->,0.004032)Update U[Sunny]=(0.0114552,Sunny->Rainy->Rainy->Sunny->,0.004032)next_state=Rainystate=Sunnyp=0.04triple=(0.0018288,Sunny->Sunny->Sunny->,0.0010368)state=Rainyp=0.35triple=(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<string> & ob, viterbi_triple_t & vtriple)
{
//alias
map<string, double>& sp = start_prob;
map<string, map<string, double> > & tp = transition_prob;
map<string, map<string, double> > & ep = emission_prob;
// initialization
InitParameters();
map<string, viterbi_triple_t> T;
for (vector<string>::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<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 (vector<string>::iterator itNextState = states.begin();
itNextState != states.end();
++itNextState)
{
cout << “\tnext_state=” << *itNextState << endl;
total = 0;
argmax.clear();
valmax = 0;
for (vector<string>::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<string>::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;
}