隐马尔可夫模型中的Viterbi算法的C++实现

[cpp] view plain copy print ?

1. #include <iostream>  
2. #include <vector>  
3. #include <map>  
4. #include <string>  
5. using namespace std;  
6. typedef vector< string > VecStr;  
7. typedef map< string, double > MapStrDou;  
8. typedef map< string, map< string, double > > MapStrMap;  
9. typedef vector< string > :: iterator VecStrI;  
10. typedef map< string, double > :: iterator MapStrDouI;  
11. typedef map< string, map< string, double > > :: iterator MapStrMapI;  
12. typedef struct DSD  
13. {  
14.     double d1;  
15.     string s;  
16.     double d2;  
17. }DSD;  
18. typedef map< string, DSD > MapStrDSD;  
19. #define RAINY "Rainy"  
20. #define SUNNY "Sunny"  
21. #define WALK "walk"  
22. #define SHOP "shop"  
23. #define CLEAN "clean"  
24. DSD forward_viterbi(VecStr& obs,  
25.         VecStr& states,  
26.         MapStrDou& start_p,  
27.         MapStrMap& trans_p,  
28.         MapStrMap& emit_p);  
29. int main(void)  
30. {  
31.     VecStr states;  
32.     VecStr observations;  
33.     MapStrDou start_probability;  
34.     MapStrMap transition_probability;  
35.     MapStrMap emission_probability;  
36.     states.push_back(RAINY);  
37.     states.push_back(SUNNY);  
38.     observations.push_back(WALK);  
39.     observations.push_back(SHOP);  
40.     observations.push_back(CLEAN);  
41.     start_probability[RAINY] = 0.6;  
42.     start_probability[SUNNY] = 0.4;  
43.     MapStrDou t1;  
44.     t1[RAINY] = 0.7;  
45.     t1[SUNNY] = 0.3;  
46.     MapStrDou t2;  
47.     t2[RAINY] = 0.4;  
48.     t2[SUNNY] = 0.6;  
49.     transition_probability[RAINY] = t1;  
50.     transition_probability[SUNNY] = t2;  
51.     MapStrDou e1;  
52.     e1[WALK] = 0.1;  
53.     e1[SHOP] = 0.4;  
54.     e1[CLEAN] = 0.5;  
55.     MapStrDou e2;  
56.     e2[WALK] = 0.6;  
57.     e2[SHOP] = 0.3;  
58.     e2[CLEAN] = 0.1;  
59.     emission_probability[RAINY] = e1;  
60.     emission_probability[SUNNY] = e2;  
61.       
62.     DSD ret = forward_viterbi(observations,  
63.         states,  
64.         start_probability,  
65.         transition_probability,  
66.         emission_probability);  
67.     cout << ret.d1 <<endl;  
68.     cout << ret.s << endl;  
69.     cout << ret.d2 << endl;  
70.       
71.     return 0;  
72. }  
73. DSD forward_viterbi(VecStr& obs,  
74.         VecStr& states,  
75.         MapStrDou& start_p,  
76.         MapStrMap& trans_p,  
77.         MapStrMap& emit_p)  
78. {  
79.     MapStrDSD T;  
80.     VecStrI VS_i, VS_j, VS_k;  
81.     DSD DSD_tmp;  
82.     for(VS_i = states.begin(); VS_i != states.end(); VS_i++)  
83.     {  
84.         DSD_tmp.d1 = start_p[*VS_i];  
85.         DSD_tmp.s = *VS_i;  
86.         DSD_tmp.d2 = start_p[*VS_i];  
87.         T[*VS_i] = DSD_tmp;  
88.     }  
89.     for(VS_i = obs.begin(); VS_i != obs.end(); VS_i++)  
90.     {  
91.         MapStrDSD U;  
92.         for(VS_j = states.begin(); VS_j != states.end(); VS_j++)  
93.         {  
94.             double total = 0.0;  
95.             string argmax = "";  
96.             double valmax = 0.0;  
97.             double prob = 1.0;  
98.             string v_path = "";  
99.             double v_prob = 1.0;  
100.             for(VS_k = states.begin(); VS_k != states.end(); VS_k++)  
101.             {  
102.                 DSD objs = T[*VS_k];  
103.                 prob = objs.d1;  
104.                 v_path = objs.s;  
105.                 v_prob = objs.d2;  
106.                 double p = (emit_p[*VS_k])[*VS_i] * (trans_p[*VS_k])[*VS_j];  
107.                 prob *= p;  
108.                 v_prob *= p;  
109.                 total += prob;  
110.                 if(v_prob > valmax)  
111.                 {  
112.                     argmax = v_path + "," + *VS_j;  
113.                     valmax = v_prob;  
114.                 }  
115.                 /** test */  
116.                 cout << "obs = " << *VS_i << endl;  
117.                 cout << "/tnext_state = " << *VS_j << endl;  
118.                 cout << "/t/tstate = " << *VS_k << endl;  
119.                 cout << "/t/t/tp = " << p << " = " << (emit_p[*VS_k])[*VS_i]   
120.                 << " * " << (trans_p[*VS_k])[*VS_j] << endl;  
121.                 cout << "/t/t/t{prob, v_path, v_prob} = {"   
122.                 << prob << ", " << v_path << ", " << v_prob << "}" << endl;  
123.                 cout << "/t/t/t{total, argmax, valmax} = {"   
124.                 << total << ", " << argmax << ", " << valmax << "}" << endl;  
125.                 /** end of test */  
126.             }  
127.             DSD_tmp.d1 = total;  
128.             DSD_tmp.s = argmax;  
129.             DSD_tmp.d2 = valmax;  
130.             U[*VS_j] = DSD_tmp;  
131.         }  
132.         T = U;  
133.     }  
134.       
135.     double total = 0.0;  
136.     string argmax = "";  
137.     double valmax = 0.0;  
138.     double prob = 1.0;  
139.     string v_path = "";  
140.     double v_prob = 1.0;  
141.     for(VS_i = states.begin(); VS_i != states.end(); VS_i++)  
142.     {  
143.         DSD objs = T[*VS_i];  
144.         prob = objs.d1;  
145.         v_path = objs.s;  
146.         v_prob = objs.d2;  
147.         total += prob;  
148.         if(v_prob > valmax)  
149.         {  
150.             argmax = v_path;  
151.             valmax = v_prob;  
152.         }  
153.     }  
154.     DSD_tmp.d1 = total;  
155.     DSD_tmp.s = argmax;  
156.     DSD_tmp.d2 = valmax;  
157.     return DSD_tmp;  
158. }  

#include <iostream> #include <vector> #include <map> #include <string> using namespace std; typedef vector< string > VecStr; typedef map< string, double > MapStrDou; typedef map< string, map< string, double > > MapStrMap; typedef vector< string > :: iterator VecStrI; typedef map< string, double > :: iterator MapStrDouI; typedef map< string, map< string, double > > :: iterator MapStrMapI; typedef struct DSD { double d1; string s; double d2; }DSD; typedef map< string, DSD > MapStrDSD; #define RAINY "Rainy" #define SUNNY "Sunny" #define WALK "walk" #define SHOP "shop" #define CLEAN "clean" DSD forward_viterbi(VecStr& obs, VecStr& states, MapStrDou& start_p, MapStrMap& trans_p, MapStrMap& emit_p); int main(void) { VecStr states; VecStr observations; MapStrDou start_probability; MapStrMap transition_probability; MapStrMap emission_probability; states.push_back(RAINY); states.push_back(SUNNY); observations.push_back(WALK); observations.push_back(SHOP); observations.push_back(CLEAN); start_probability[RAINY] = 0.6; start_probability[SUNNY] = 0.4; MapStrDou t1; t1[RAINY] = 0.7; t1[SUNNY] = 0.3; MapStrDou t2; t2[RAINY] = 0.4; t2[SUNNY] = 0.6; transition_probability[RAINY] = t1; transition_probability[SUNNY] = t2; MapStrDou e1; e1[WALK] = 0.1; e1[SHOP] = 0.4; e1[CLEAN] = 0.5; MapStrDou e2; e2[WALK] = 0.6; e2[SHOP] = 0.3; e2[CLEAN] = 0.1; emission_probability[RAINY] = e1; emission_probability[SUNNY] = e2; DSD ret = forward_viterbi(observations, states, start_probability, transition_probability, emission_probability); cout << ret.d1 <<endl; cout << ret.s << endl; cout << ret.d2 << endl; return 0; } DSD forward_viterbi(VecStr& obs, VecStr& states, MapStrDou& start_p, MapStrMap& trans_p, MapStrMap& emit_p) { MapStrDSD T; VecStrI VS_i, VS_j, VS_k; DSD DSD_tmp; for(VS_i = states.begin(); VS_i != states.end(); VS_i++) { DSD_tmp.d1 = start_p[*VS_i]; DSD_tmp.s = *VS_i; DSD_tmp.d2 = start_p[*VS_i]; T[*VS_i] = DSD_tmp; } for(VS_i = obs.begin(); VS_i != obs.end(); VS_i++) { MapStrDSD U; for(VS_j = states.begin(); VS_j != states.end(); VS_j++) { double total = 0.0; string argmax = ""; double valmax = 0.0; double prob = 1.0; string v_path = ""; double v_prob = 1.0; for(VS_k = states.begin(); VS_k != states.end(); VS_k++) { DSD objs = T[*VS_k]; prob = objs.d1; v_path = objs.s; v_prob = objs.d2; double p = (emit_p[*VS_k])[*VS_i] * (trans_p[*VS_k])[*VS_j]; prob *= p; v_prob *= p; total += prob; if(v_prob > valmax) { argmax = v_path + "," + *VS_j; valmax = v_prob; } /** test */ cout << "obs = " << *VS_i << endl; cout << "/tnext_state = " << *VS_j << endl; cout << "/t/tstate = " << *VS_k << endl; cout << "/t/t/tp = " << p << " = " << (emit_p[*VS_k])[*VS_i] << " * " << (trans_p[*VS_k])[*VS_j] << endl; cout << "/t/t/t{prob, v_path, v_prob} = {" << prob << ", " << v_path << ", " << v_prob << "}" << endl; cout << "/t/t/t{total, argmax, valmax} = {" << total << ", " << argmax << ", " << valmax << "}" << endl; /** end of test */ } DSD_tmp.d1 = total; DSD_tmp.s = argmax; DSD_tmp.d2 = valmax; U[*VS_j] = DSD_tmp; } T = U; } double total = 0.0; string argmax = ""; double valmax = 0.0; double prob = 1.0; string v_path = ""; double v_prob = 1.0; for(VS_i = states.begin(); VS_i != states.end(); VS_i++) { DSD objs = T[*VS_i]; prob = objs.d1; v_path = objs.s; v_prob = objs.d2; total += prob; if(v_prob > valmax) { argmax = v_path; valmax = v_prob; } } DSD_tmp.d1 = total; DSD_tmp.s = argmax; DSD_tmp.d2 = valmax; return DSD_tmp; }

---

执行结果：

obs = walk
    next_state = Rainy
        state = Rainy
            p = 0.07 = 0.1 * 0.7
            {prob, v_path, v_prob} = {0.042, Rainy, 0.042}
            {total, argmax, valmax} = {0.042, Rainy,Rainy, 0.042}
obs = walk
    next_state = Rainy
        state = Sunny
            p = 0.24 = 0.6 * 0.4
            {prob, v_path, v_prob} = {0.096, Sunny, 0.096}
            {total, argmax, valmax} = {0.138, Sunny,Rainy, 0.096}
obs = walk
    next_state = Sunny
        state = Rainy
            p = 0.03 = 0.1 * 0.3
            {prob, v_path, v_prob} = {0.018, Rainy, 0.018}
            {total, argmax, valmax} = {0.018, Rainy,Sunny, 0.018}
obs = walk
    next_state = Sunny
        state = Sunny
            p = 0.36 = 0.6 * 0.6
            {prob, v_path, v_prob} = {0.144, Sunny, 0.144}
            {total, argmax, valmax} = {0.162, Sunny,Sunny, 0.144}
obs = shop
    next_state = Rainy
        state = Rainy
            p = 0.28 = 0.4 * 0.7
            {prob, v_path, v_prob} = {0.03864, Sunny,Rainy, 0.02688}
            {total, argmax, valmax} = {0.03864, Sunny,Rainy,Rainy, 0.02688}
obs = shop
    next_state = Rainy
        state = Sunny
            p = 0.12 = 0.3 * 0.4
            {prob, v_path, v_prob} = {0.01944, Sunny,Sunny, 0.01728}
            {total, argmax, valmax} = {0.05808, Sunny,Rainy,Rainy, 0.02688}
obs = shop
    next_state = Sunny
        state = Rainy
            p = 0.12 = 0.4 * 0.3
            {prob, v_path, v_prob} = {0.01656, Sunny,Rainy, 0.01152}
            {total, argmax, valmax} = {0.01656, Sunny,Rainy,Sunny, 0.01152}
obs = shop
    next_state = Sunny
        state = Sunny
            p = 0.18 = 0.3 * 0.6
            {prob, v_path, v_prob} = {0.02916, Sunny,Sunny, 0.02592}
            {total, argmax, valmax} = {0.04572, Sunny,Sunny,Sunny, 0.02592}
obs = clean
    next_state = Rainy
        state = Rainy
            p = 0.35 = 0.5 * 0.7
            {prob, v_path, v_prob} = {0.020328, Sunny,Rainy,Rainy, 0.009408}
            {total, argmax, valmax} = {0.020328, Sunny,Rainy,Rainy,Rainy, 0.009408}
obs = clean
    next_state = Rainy
        state = Sunny
            p = 0.04 = 0.1 * 0.4
            {prob, v_path, v_prob} = {0.0018288, Sunny,Sunny,Sunny, 0.0010368}
            {total, argmax, valmax} = {0.0221568, Sunny,Rainy,Rainy,Rainy, 0.009408}
obs = clean
    next_state = Sunny
        state = Rainy
            p = 0.15 = 0.5 * 0.3
            {prob, v_path, v_prob} = {0.008712, Sunny,Rainy,Rainy, 0.004032}
            {total, argmax, valmax} = {0.008712, Sunny,Rainy,Rainy,Sunny, 0.004032}
obs = clean
    next_state = Sunny
        state = Sunny
            p = 0.06 = 0.1 * 0.6
            {prob, v_path, v_prob} = {0.0027432, Sunny,Sunny,Sunny, 0.0015552}
            {total, argmax, valmax} = {0.0114552, Sunny,Rainy,Rainy,Sunny, 0.004032}
0.033612
Sunny,Rainy,Rainy,Rainy
0.009408

 

---

参考文献：

http://en.wikipedia.org/wiki/Viterbi_algorithm
http://hi.baidu.com/liqinghuisi/blog/item/a11cd4f54394f22fbc3109cf.html