HMM Viterbi算法SAS代码_Andy的原创帖

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

起源:https://en.wikipedia.org/wiki/Viterbi_algorithm
在维基上看了关于维特比算法的解释和代码，发现没有SAS的实现代码，在百度上也搜索不到，我用SAS IML写了一段，可以对比原文阅读。帮助使用SAS的朋友们理解维特比。
*PS：代码可能还不够漂亮，但我觉得毛坯比较容易理解和检查(对于初学者)
/*Trial Complete*/
proc iml;
so={normal cold dizzy}; /*观察序列*/
/*-->so转换成列代码*/
ivec=i(2);
so_t={1 2 3 };
/*1 2 3 代表 normal cold dizzy
  1 2 代表 Healthy Fever*/
s={Healthy,Fever};/* 隐含状态定义 不加引号就会变成纯大写*/
o={normal,cold,dizzy};/*观察状态定义 observation*/
pi={0.6 0.4};/*start_probability 初始矩阵*/
tp={0.7 0.3,0.4 0.6};  /*transition_probability 隐含状态转移矩阵~健康状态*/
ep={0.5 0.4 0.1,0.1 0.3 0.6};/*emission_probability 发射状态概率分布矩阵~马尔科夫状态转移矩阵~病人的感觉:正常、眩晕、发烧*/


path={0};/*路径*/


pr={0,0};/*对应的概率*/
/*1 - 第一次 */
h=pi;
obs=so_t[1];
com=h#t(ep[,obs]);
mpos=com[<:>];
path=path||mpos;
pr=pr||t(com);




/*2 - 第二次 */
h1=t(pr[,2]);
h=h1*tp;
obs=so_t[2];
com=h#t(ep[,obs]);
mpos=com[<:>];
path=path||mpos;


wpos=path[2];
wvec=ivec[wpos,];


h2=pr[wpos,2]#wvec;
h=h2*tp;
com=h#t(ep[,obs]);
pr=pr||t(com);


/*3 第三次 */
h1=t(pr[,3]);
h=h1*tp;
obs=so_t[3];
com=h#t(ep[,obs]);
mpos=com[<:>];
path=path||mpos;


wpos=path[3];
wvec=ivec[wpos,];


h2=pr[wpos,3]#wvec;
h=h2*tp;
com=h#t(ep[,obs]);
pr=pr||t(com);


print path;
print pr;
quit;

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：sas代码 ITER Andy HMM vit 维特比 HMM SAS 隐马尔科夫