作者：jostree 轉載請注明出處 http://www.cnblogs.com/jostree/p/4335810.html
一個例子：
韋小寶使用骰子進行遊戲，他有兩種骰子一種正常的骰子，還有一種不均勻的骰子，來進行出千。
開始遊戲時他有2/5的機率出千。
對于正常的骰子A，每個點出現的機率都是1/6.
對于不均勻的骰子B，5，6兩種出現的機率為3/10，其餘為1/10.
出千的随機規律如下圖所示：
我們觀測到的投擲結果為：ob={1，3，4，5，5，6，6，3，2，6}
請判斷韋小寶什麼時候出千了？
我們可以這樣模組化$x_i$表示第$i$次投擲的骰子的種類，$y_i$表示第$i$次投擲出的點數，$\lambda$表示各個機率參數。
那麼第$t$次使用第$i$種骰子投擲的機率$\delta_t(i)$等于
\begin{equation} \delta_t(i)=\max_{x_1,\dots,x_{t-1}}P(x_1,\dots,x_{t-1},x_t=i,y_1,\dots,y_t|\lambda) \end{equation}
其實$\delta_{t+1}(i)$可以由$\delta_t(i)$推倒得出：
\begin{eqnarray} \delta_{t+1}(i) &=& \max_{x_1,\dots,x_{t}}P(x_1,\dots,x_{t},x_{t+1}=i,y_1,\dots,y_{t+1}|\lambda)\\ &=& \max_j \delta_t(j)\alpha_{ji}\beta_i(y_{t+1})\end{eqnarray}
其中$\alpha_{ji}$表示從第$j$個骰子轉移到第$i$個骰子的機率。
$\beta_i(y_{t+1})$表示使用第i個骰子投出點$y_{t+1}$的機率。
進而可以使用上述利用動态規劃算法進行逐次遞推計算。
得到的結果為：
| t | $y_t$ | $\delta_t(A)$ | $\Psi_t(A)$ | $\delta_t(B)$ | $\Psi_t(B)$ |
| 1 | 0.1 | A | 0.04 | ||
| 2 | 3 | 0.0133333 | 0.0036 | B | |
| 4 | 0.00177778 | 0.000324 | |||
| 5 | 0.000237037 | 0.000106667 | |||
| 3.16049e-05 | 2.88e-05 | ||||
| 6 | 4.21399e-06 | 7.776e-06 | |||
| 7 | 5.61866e-07 | 2.09952e-06 | |||
| 8 | 7.49154e-08 | 1.88957e-07 | |||
| 9 | 9.98872e-09 | 1.70061e-08 | |||
| 10 | 1.33183e-09 | 4.59165e-09 |
因為最後一步$\delta_t(B)$的值大于$\delta_t(A)$，是以一次使用B骰子的機率最大，進而一直向上回溯，得到的使用骰子的序列為：AAABBBBBBB
代碼如下所示：
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <string.h>
4 #include <string>
5 #include <iostream>
6 using namespace std;
7 double initP[2] = {0.6, 0.4};//骰子A，B的初始機率
8 double transferMatrix[2][2] = {{0.8, 0.2}, {0.1, 0.9}};//骰子之間的轉移機率
9 double EmissionP[2][6]={{1/6.0, 1/6.0, 1/6.0, 1/6.0, 1/6.0, 1/6.0},//骰子A的發射機率
10 {0.1, 0.1, 0.1, 0.1, 0.3, 0.3}};//骰子B的發射機率
11 double dp[10][2];//dp[i][j]第i步時，使用第j個骰子的最大機率
12 double dpS[10][2];//dpS[i][j]第i步時，使用第j個骰子，得到的最大機率時，使用的骰子種類, 0->A, 1->B
13 int ob[10] = {1, 3, 4, 5, 5, 6, 6, 3, 2, 6};//觀測點數
14 bool diceArray[10];//預測骰子使用序列
15 void Viterbi()
16 {
17 memset(dp,0,sizeof(dp));
18 memset(dpS,0,sizeof(dpS));
19 memset(diceArray,0,sizeof(diceArray));
20 dp[0][0] = initP[0]* EmissionP[0][ob[0]-1];
21 dp[0][1] = initP[1]* EmissionP[1][ob[0]-1];
22 for( int i = 1 ; i < 10 ; i++ )//投擲次數
23 {
24 for( int j = 0 ; j < 2 ; j++ )//目前狀态
25 {
26 for( int k = 0 ; k < 2 ; k++ )//上一個狀态
27 {
28 double tempP = dp[i-1][k] * transferMatrix[k][j] * EmissionP[j][ob[i]-1] ;
29 if( dp[i][j] < tempP )
30 {
31 dp[i][j] = tempP;
32 dpS[i][j] = k;
33 }
34 }
35 }
36 }
37 int maxState = 0;
38 if( dpS[9][0] < dpS[9][1] )
39 {
40 maxState = 1;
41 }
42 for( int i = 9 ; i >=0 ; i-- )
43 {
44 diceArray[i] = maxState;
45 maxState = dpS[i][maxState];
46 }
47 }
48 int main(int argc, char *argv[])
49 {
50 Viterbi();
51 cout<<"每步每個狀态下的機率和骰子種類:"<<endl;
52 for( int i = 0 ; i < 10 ; i++ )
53 {
54 for( int j = 0 ; j < 2 ; j++ )
55 {
56 cout<<dp[i][j]<<" "<<dpS[i][j]<<" ";
57 }
58 cout<<endl;
59 }
60 cout<<"預測骰子種類，0->A, 1->B : "<<endl;
61 for( int i = 0 ; i < 10 ; i++ )
62 {
63 cout<<diceArray[i]<<" ";
64 }
65 cout<<endl;
66 }
67 /* result:
68 每步每個狀态下的機率和骰子種類:
69 0.1 0 0.04 0
70 0.0133333 0 0.0036 1
71 0.00177778 0 0.000324 1
72 0.000237037 0 0.000106667 0
73 3.16049e-05 0 2.88e-05 1
74 4.21399e-06 0 7.776e-06 1
75 5.61866e-07 0 2.09952e-06 1
76 7.49154e-08 0 1.88957e-07 1
77 9.98872e-09 0 1.70061e-08 1
78 1.33183e-09 0 4.59165e-09 1
79 預測骰子種類，0->A, 1->B :
80 0 0 0 1 1 1 1 1 1 1
81 */ View Code