HMM Forward算法SAS代码_Andy的原创帖12

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

感谢下两个资料:
外文原版:
http://www.comp.leeds.ac.uk/roger/HiddenMarkovModels/html_dev/main.html
中文翻译版:
http://blog.csdn.net/daringpig/article/details/8072794
****************************************************************
首先觉得Forward算法并不难，关键是理解起来有点绕。特别是放到矩阵运算下解释很简单，但是写成公式就真的真的很费解。
老外原版上把HMM说的很明白，不过在实例的applet小程序里有点小错误，这点我和翻译版对过的结果是一致的，以下我会梳理一遍，然后把我的理解和sas 代码附上。(代码没有写成循环或迭代的形式，后续我再整理的时候会写)
具体的理论就不说了，参照外文原版即可，以下我把例子的参数和验证过程写一下，我觉得代码看懂了就很简单。
以下三个是HMM的基本参数，定义了隐含及观察状态数，初始向量，转移矩阵，发射矩阵(混淆矩阵)


算法的结构图:
原文如下:

我理解的结构如下:

*中间的一堆线事实上就是向量和转移矩阵相乘。
*P(O|HMM)=P(O|SUNNY)+P(O|CLOUDY)+P(O|RAINY)
*根据观察到的状态进行发射(确定的隐藏状态向量乘以对应状态的发射向量)
*每个空格都创建了存储空间，以便递归调用上一步的概率
---动态规划的本质也就是用空间换时间
说说原版示例中碰到的bug: 行向量*列向量 错写为 行向量*行向量
如果用矩阵操作的话，一般是不会出这个错的，因为行向量乘以转移矩阵得到的还是行向量 1*3 X 3*3= 1*3
翻译版中的解释:(图的大小调不来，有兴趣可以到翻译版中原文看)

我的验证:

第二部的sunny也是

“0.0303188” 与翻译版中的“0.030319"应该是一致的。
最终的概率结果:



和翻译版中 “0.026901"应该也是一致的。
据说(本机暂时跑不了原版的java applet)原版的结果是0.027386915
所以原版应该出错了。
以下是sas代码,有兴趣可以把print的部分放到对应的步骤后执行(我还是觉得矩阵比数组写起来节约太多代码了)

proc iml;

/*****************定义部分开始**********************/

/*sunny->sunny cloudy rainy*/

tp1={0.5 0.375 0.125};

/*cloudy->sunny cloudy rainy*/

tp2={0.25 0.125 0.625};

/*rainy->sunny cloudy rainy*/

tp3={0.25 0.375 0.375};

/*转移矩阵TP*/

tp=tp1//tp2//tp3;

/*2 发射概率矩阵*/

/*sunny->drydryish damp soggy*/

ep1={0.6 0.2 0.15 0.05};

/*cloudy->drydryish damp soggy*/

ep2={0.25 0.25 0.25 0.25};

/*rainy->drydryish damp soggy*/

ep3={0.05 0.10 0.35 0.50};

/*发射概率EP*/

ep=ep1//ep2//ep3;

pi={0.63 0.17 0.2};

/*****************定义部分结束**********************/

/*序列 dry damp soggy*/

/*第一次观察_dry*/

i=1;

o=ep[,i]`;

/*第一次隐含状态是pi*/

h=pi#o;

/*计入空间*/

t=h`;

/*进入第二轮_damp*/

i=3;

o=ep[,i]`;

/*第二次的隐含状态从第一次取*/

tem=t[,1]`;

h=(tem*tp)#o;

t=t||h`;

/*进入第三轮_soggy*/

i=4;

o=ep[,i]`;

/*第三次的隐含状态从第二次取*/

tem=t[,2]`;

h=(tem*tp)#o;

t=t||h`;

p=t[+,];

print o;

print h;

print t;

print p;

quit;

扫码加我 拉你入群

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

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

分享0 收藏1 回帖

关键词：forward sas代码 ward Andy RWA 中文翻译 程序 外文 资料

把前向算法打包成宏
%macro forward(os);
proc iml;
/*sunny cloudy rainy*/
a1={0.5 0.375 0.125};
a2={0.25 0.125 0.625};
a3={0.25 0.375 0.375};
a=a1//a2//a3;
/*dry dryish damp soggy*/
b1={0.6 0.2 0.15 0.05};
b2={0.25 0.25 0.25 0.25};
b3={0.05 0.1 0.35 0.5};
b=b1//b2//b3;
pi={0.63 0.17 0.2};
os={&os};

/*n=1*/
/*隐含状态*/
hs=pi;
/*观察向量*/
obs=b[,os[1]]`;
alpha=(hs#obs)`;
do i=2 to ncol(os);
/*n>=2*/
/*通过前一时刻alpha和转移矩阵获得新的隐藏状态*/
hs=(alpha[,i-1])`*a;
/*获得新的观察向量*/
obs=b[,os[i]]`;
alpha=alpha||(hs#obs)`;
end;
apr=alpha[+,];
print  alpha apr;
quit;
%mend;
%forward(1 3 4 )