出处:http://www.52nlp.cn/hmm-learn-best-practices-five-forward-algorithm-5
在HMM这个翻译系列的原文中,作者举了一个前向算法的交互例子,这也是这个系列中比较出彩的地方,但是,在具体运行这个例子的时候,却发现其似乎有点问题。
先说一下如何使用这个交互例子,运行时需要浏览器支持java,我用的是firefox。首先在Set按钮前面的对话框里上观察序列,如“Dry,Damp, Soggy” 或“Dry Damp Soggy”,观察符号间用逗号或空格隔开;然后再点击Set按钮,这样就初始化了观察矩阵;如果想得到一个总的结果,即Pr(观察序列|隐马尔科夫模型),就点旁边的Run按钮;如果想一步一步观察计算过程,即每个节点的局部概率,就单击旁边的Step按钮。
原文交互例子(即天气这个例子)中所定义的已知隐马尔科夫模型如下:
1、隐藏状态 (天气):Sunny,Cloudy,Rainy;
2、观察状态(海藻湿度):Dry,Dryish,Damp,Soggy;
3、初始状态概率: Sunny(0.63), Cloudy(0.17), Rainy(0.20);
4、状态转移矩阵:
weather today
Sunny Cloudy Rainy
weather Sunny 0.500 0.375 0.125
yesterday Cloudy 0.250 0.125 0.625
Rainy 0.250 0.375 0.375
5、混淆矩阵:
observed states
Dry Dryish Damp Soggy
Sunny 0.60 0.20 0.15 0.05
hidden Cloudy 0.25 0.25 0.25 0.25
states Rainy 0.05 0.10 0.35 0.50
为了UMDHMM也能运行这个例子,我们将上述天气例子中的隐马尔科夫模型转化为如下的UMDHMM可读的HMM文件weather.hmm:
——————————————————————–
M= 4
N= 3
A:
0.500 0.375 0.125
0.250 0.125 0.625
0.250 0.375 0.375
B:
0.60 0.20 0.15 0.05
0.25 0.25 0.25 0.25
0.05 0.10 0.35 0.50
pi:
0.63 0.17 0.20
——————————————————————–
在运行例子之前,如果读者也想观察每一步的运算结果,可以将umdhmm-v1.02目录下forward.c中的void Forward(…)函数替换如下:
——————————————————————–
void Forward(HMM *phmm, int T, int *O, double **alpha, double *pprob)
{
int i, j; /* state indices */
int t; /* time index */
double sum; /* partial sum */
/* 1. Initialization */
for (i = 1; i <= phmm->N; i++)
{
alpha[1][i] = phmm->pi[i]* phmm->B[i][O[1]];
printf( “a[1][%d] = pi[%d] * b[%d][%d] = %f * %f = %f\\n”,i, i, i, O[i], phmm->pi[i], phmm->B[i][O[1]], alpha[1][i] );
}
/* 2. Induction */
for (t = 1; t < T; t++)
{
for (j = 1; j <= phmm->N; j++)
{
sum = 0.0;
for (i = 1; i <= phmm->N; i++)
{
sum += alpha[t][i]* (phmm->A[i][j]);
printf( “a[%d][%d] * A[%d][%d] = %f * %f = %f\\n”, t, i, i, j, alpha[t][i], phmm->A[i][j], alpha[t][i]* (phmm->A[i][j]));
printf( “sum = %f\\n”, sum );
}
alpha[t+1][j] = sum*(phmm->B[j][O[t+1]]);
printf( “a[%d][%d] = sum * b[%d][%d]] = %f * %f = %f\\n”,t+1, j, j, O[t+1], sum, phmm->B[j][O[t+1]], alpha[t+1][j] );
}
}
/* 3. Termination */
*pprob = 0.0;
for (i = 1; i <= phmm->N; i++)
{
*pprob += alpha[T][i];
printf( “alpha[%d][%d] = %f\\n”, T, i, alpha[T][i] );
printf( “pprob = %f\\n”, *pprob );
}
}
——————————————————————–
替换完毕之后,重新“make clean”,“make all”,这样新的testfor可执行程序就可以输出前向算法每一步的计算结果。
现在我们就用testfor来运行原文中默认给出的观察序列“Dry,Damp,Soggy”,其所对应的UMDHMM可读的观察序列文件test1.seq:
——————————————————————–
T=3
1 3 4
——————————————————————–
好了,一切准备工作就绪,现在就输入如下命令:
testfor weather.hmm test1.seq > result1
result1就包含了所有的结果细节:
——————————————————————–
Forward without scaling
a[1][1] = pi[1] * b[1][1] = 0.630000 * 0.600000 = 0.378000
a[1][2] = pi[2] * b[2][3] = 0.170000 * 0.250000 = 0.042500
a[1][3] = pi[3] * b[3][4] = 0.200000 * 0.050000 = 0.010000
…
pprob = 0.026901
log prob(O| model) = -3.615577E+00
prob(O| model) = 0.026901
…
——————————————————————–
黑体部分是最终的观察序列的概率结果,即本例中的Pr(观察序列|HMM) = 0.026901。
但是,在原文中点Run按钮后,结果却是:Probability of this model = 0.027386915。
这其中的差别到底在哪里?我们来仔细观察一下中间运行过程:
在初始化亦t=1时刻的局部概率计算两个是一致的,没有问题。但是,t=2时,在隐藏状态“Sunny”的局部概率是不一致的。英文原文给出的例子的运行结果是:
Alpha = (((0.37800002*0.5) + (0.0425*0.375) + (0.010000001*0.125)) * 0.15) = 0.03092813
而UMDHMM给出的结果是:
——————————————————————–
a[1][1] * A[1][1] = 0.378000 * 0.500000 = 0.189000
sum = 0.189000
a[1][2] * A[2][1] = 0.042500 * 0.250000 = 0.010625
sum = 0.199625
a[1][3] * A[3][1] = 0.010000 * 0.250000 = 0.002500
sum = 0.202125
a[2][1] = sum * b[1][3]] = 0.202125 * 0.150000 = 0.030319
——————————————————————–
区别就在于状态转移概率的选择上,原文选择的是状态转移矩阵中的第一行,而UMDHMM选择的则是状态转移矩阵中的第一列。如果从原文给出的状态转移矩阵来看,第一行代表的是从前一时刻的状态“Sunny”分别到当前时刻的状态“Sunny”,“Cloudy”,“Rainy”的概率;而第一列代表的是从前一时刻的状态“Sunny”,“Cloudy”,“Rainy”分别到当前时刻状态“Sunny”的概率。这样看来似乎原文的计算过程有误,读者不妨多试几个例子看看,前向算法这一章就到此为止了。