在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”的概率。这样看来似乎原文的计算过程有误,读者不妨多试几个例子看看,前向算法这一章就到此为止了。

未完待续:维特比算法1

本翻译系列原文请参考:http://www.comp.leeds.ac.uk/roger/HiddenMarkovModels/html_dev/main.html

注:原创文章,转载请注明出处“我爱自然语言处理”:www.52nlp.cn

本文链接地址:https://www.52nlp.cn/hmm-learn-best-practices-five-forward-algorithm-5

作者 52nlp

《HMM学习最佳范例五:前向算法5》有57条评论
  1. 今天尝试按照关键字搜索了一下你的网站,确实很有效果,做自己想做并且喜欢做的事情,你的坚持让我很感动。

    [回复]

    admin 回复:

    恩,你也一样!

    [回复]

  2. 前辈,多谢你的辛苦翻译,让我对隐马有了更深层次的了解,最近我用HMM来对时间序列建模,有两个问题想向你请教一下,不知道能否留个联系方式?
    我的QQ:627431239,多谢~

    [回复]

  3. 我先看到英文原文,实现了算法,对比了结果,也觉得那个转移矩阵不对,正在怀疑之时,看到了你这篇文章,谢谢!

    [回复]

    52nlp 回复:

    不客气!尽信书不如无书吧!

    [回复]

  4. 请问大侠:HMM用于时间序列建立模型,特别是考虑到输入对输出结果的影响,如何解决连续性\小样本和控制输入问题?谢谢

    [回复]

    52nlp 回复:

    不好意思,这个问题我不太清楚,我对HMM的理解还仅限于离散这一块儿。

    [回复]

  5. 请问
    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
    是不是错了?是不是应:
    a[1][1] = pi[1] * b[1][1] = 0.630000 * 0.600000 = 0.378000
    a[1][2] = pi[2] * b[2][1] = 0.170000 * 0.250000 = 0.042500
    a[1][3] = pi[3] * b[3][1] = 0.200000 * 0.050000 = 0.010000

    [回复]

    52nlp 回复:

    好像没错吧,1,3,4代表的是观察序列。不过好久了,我都有点忘了,呵呵。

    [回复]

    xiao 回复:

    这个问题我认为博主有点小笔误, yxs 更正的是正确的,即,正确结果如下:
    a[1][1] = pi[1] * b[1][1] = 0.630000 * 0.600000 = 0.378000
    a[1][2] = pi[2] * b[2][1] = 0.170000 * 0.250000 = 0.042500
    a[1][3] = pi[3] * b[3][1] = 0.200000 * 0.050000 = 0.010000
    b的第二个角标都是1。因为观察序列1,3,4表示在t=1时的观察值为1,在t=2时的观察值为3,在t=3时的观察值为4,所以在求解t=1时的局部概率时,应该分别取混淆矩阵中的b[1][1], b[2][1],b[3][1]。另外,在翻译混淆矩阵的时候,个人认为混合矩阵稍显恰当。以上仅是个人拙见,请博主指点。再次感谢博主辛勤的努力,您的博客对我很有帮助。

    [回复]

    52nlp 回复:

    谢谢质疑,主要这个结果是程序给出的,我没有手写,所以不太确定,都有点糊涂了。

    [回复]

    xiao 回复:

    博主好,谢谢您的认真回复,我的理解是这样的:您程序中的initialization部分,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] );
    中的O[i]写的有误,应该是O[1],因为现在是计算t=1时的局部概率,您在这句中后面用的就是phmm->B[i][O[1]],也就是表明现在要乘的是O[1]对应的元素。
    同时,您这里的小笔误不影响计算的正确性,因为在输出语句的上一句,也就是真正局部概率的那一句alpha[1][i] = phmm->pi[i]* phmm->B[i][O[1]];您正确使用了B[i][O[1]]。以上仅是个人拙见,请大家讨论。
    再次感谢楼主认真细致的学术品格和细心翻译,我学到了很多。

    52nlp 回复:

    最近实在太忙了,我元旦的时候再看一下这个问题,谢谢你这么细致的检查!

  6. 前辈,我想问一下,UMDHMM前向算法中 初始向量都是1 怎么解释呢?代码如下:
    /* 1. Initialization */

    for (i = 1; i N; i++)
    beta[T][i] = 1.0;

    还有关于带有加权参量的向前算法,这个加权变量scale[t] 是一个时间变量,是不是可以解决部分状态转移和时间无关的缺点?

    [回复]

    52nlp 回复:

    beta[T][i] = 1.0 这个应该是后向算法里的吧?

    关于scale[t],抱歉我不太清楚是否能解决这个问题!

    [回复]

    simon 回复:

    哦 是的 猪头了一回

    [回复]

  7. 英文原文确实有误,状态转移矩阵的选择不对
    我用穷举法验证了一下
    P(obervation|HMM)= 0.026901406

    [回复]

  8. 我想用来做手势识别,但是训练以后在计算前向概率的时候概率值很低很低,一般一个序列15的长度计算得到的概率值都是-30多,不知道这样是否正常?谢谢了

    [回复]

    52nlp 回复:

    不清楚,只要有可比性,应该是正常的吧。

    [回复]

    袁军 回复:

    我也碰到了相同的问题,不知道您已经解决了吗?

    [回复]

  9. 看到umdhmm的hmm.c里面对于M的定义是 离散的观测符号的个数,请问:

    如果观测的概率分布不是离散的,而是高斯分布,还能不能使用umdhmm?

    Thanks.

    [回复]

    52nlp 回复:

    不清楚,不过即使是连续的分布,对于M也只是其中的一些观察值,代表的是一些抽样的点吧。

    [回复]

  10. 博主您好,请教您两个问题,1.隐马尔科夫模型多观察序列的概率应该怎么求?2.给定一个隐马尔科夫,用前向算法求一个观察序列的概率,但这个观察序列的T=100,算到后面alpha的值都为零了,是怎么回事?望博主赐教!

    [回复]

    52nlp 回复:

    可以考虑用scaling 的方法,建议看这个文章的第6节:HMM Scaling

    http://www.cs.sjsu.edu/~stamp/RUA/HMM.pdf

    [回复]

    袁军 回复:

    谢谢博主,这是一种数据下溢问题,需增加比例因子,程序中我调了好久也没弄好,请问博主有程序这方面资料吗?

    [回复]

    52nlp 回复:

    抱歉,不太清楚。

  11. 老师,我按照你说的执行了下,为什么在linux终端会出现:
    yuemeng@yuemeng-desktop:~/下载/umdhmm/umdhmm-v1.02$ ./testfor weather.hmm test1.seq
    bash: ./testfor: cannot execute binary file
    yuemeng@yuemeng-desktop:~/下载/umdhmm/umdhmm-v1.02$
    这样的提示呢,请问那地方出错了呢,上面有我程序的路径,请指导下啊
    非常感谢

    [回复]

    sunmoon 回复:

    我遇到和你同样的问题,请问你是怎么解决的

    [回复]

  12. 感谢博主的文章. 我刚看到这里,感觉文章中的:
    Alpha = (((0.37800002*0.5) + (0.0425*0.375) + (0.010000001*0.125)) * 0.15) = 0.03092813 应该是对的,而程序里的(0.5,0.25,0.25)是错误的. 因为是从其他隐藏状态向第一个隐藏状态转移的概率,所以使用第一列来计算,而不是第一行. 不知道我理解的对不对.

    [回复]

    mike98 回复:

    Alpha = (((0.37800002*0.5) + (0.0425*0.375) + (0.010000001*0.125)) * 0.15) = 0.03092813
    肯定是错的。0.375在后面给出的数据中应该为0.325。另外,从转移矩阵上看,行为昨天,列为今天,昨天到今天显然应该取列向量。

    [回复]

  13. 我用楼主的程序和数据实验,结果居然全是0。不知道怎么回事?

    [回复]

    William 回复:

    在Mac OS上似乎ReadHMM()无法正确的得到hmm中的数据。验证了一个hmm.pi[1] = 0.000000

    [回复]

    52nlp 回复:

    是在mac上编译的吗?

    [回复]

    William 回复:

    是的,我简单的跟了一下,应该是读取hmm数据时,数据类型(double)的问题,可能原来的实现只是在Linux平台。

    William 回复:

    好奇怪。。ReadHMM() 中的fscanf(fp, "%lf", &(phmm->A[i][j])); 居然不work。

  14. 感谢对原文中错误的指正,当时看到原文就觉得很疑惑感觉他是不是写错了。但是限于自己水平不能确定。

    [回复]

  15. 原文确实是有些问题,有可能是老外看transition matrix的顺序和我们不一样

    [回复]

  16. 我在linue下进入目录,make all ,显示结果‘all’ is up to data

    [回复]

    豹豹 回复:

    make clean
    make & make all

    [回复]

  17. 你好,楼主,,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] );,,我在linux下运行这个时,为什么报错说的a[][]和b[][]没有声明

    [回复]

  18. 老师,你好,
    forward.c: In function ‘Forward’:
    forward.c:28: error: stray ‘\342’ in program
    forward.c:28: error: stray ‘\200’ in program
    forward.c:28: error: stray ‘\234’ in program
    forward.c:28: error: ‘a’ undeclared (first use in this function)
    forward.c:28: error: (Each undeclared identifier is reported only once
    forward.c:28: error: for each function it appears in.)
    forward.c:28: error: expected expression before ‘%’ token
    forward.c:28: error: ‘pi’ undeclared (first use in this function)
    forward.c:28: error: expected expression before ‘%’ token
    forward.c:28: error: ‘b’ undeclared (first use in this function)
    forward.c:28: error: expected expression before ‘%’ token
    forward.c:28: error: expected expression before ‘%’ token
    forward.c:28: error: expected expression before ‘%’ token
    forward.c:28: error: stray ‘\’ in program
    forward.c:28: error: stray ‘\’ in program
    forward.c:28: error: stray ‘\342’ in program
    forward.c:28: error: stray ‘\200’ in program
    forward.c:28: error: stray ‘\235’ in program
    make: *** [forward.o] Error 1
    这些报错,怎么办

    [回复]

    52nlp 回复:

    你用的什么编译器?

    [回复]

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注