博文
HMM模型的初步学习
|||
这几天学习了HMM模型，得到了下面的一些粗浅认识，以及对于一个简单的HMM模型的实现（Perl语言实现的）。code写得很冗余，要是把计算过程写成函数会更好。
本文只是自己的一个简单的总结，需要了解的更详细可以参考Lawrence.R. Rabiner的经典论文"A tutorial on Hidden Markov Models and selected applications in Speech recognition"。
1 什么是HMM模型？隐马尔科夫模型(Hidden Markov Model),它是一种统计学模型，用来描述随机事件的发生。该模型与马尔科夫模型的差别在于它的状态是隐含的，我们看到的是事件的发生。隐状态之间有转换概率，从隐状态到事件之间有映射概率。
2 HMM模型能够干什么？
三个问题：
1）给出HMM模型和观测序列，计算在该模型下得到该观测序列的概率
2）给出HMM模型和观测序列，计算在该模型下得到该观测序列最有可能的隐状态序列
3）给定观测序列，估计HMM模型
3 对于每个问题的解答以及例子
对于问题1，一般采用的是forward algorithm(也可以用backward algorithm)。该算法的实现如下（Perl语言）：
#指定HMM模型的参数
my $stateNum=3;
my $obsNum=10;
my @initialP=qw(0.1 0.5 0.4);
my @state=qw(A B C);
my @obsKind=qw(1 2 3 4);
my $obsKind2num;
for(my $i=0;$i<scalar @obsKind;$i++){
$obsKind2num->{$obsKind[$i]}=$i;
}
my @Atransition=(0.8,0.1,0.1);
my @Btransition=(0.2,0.6,0.2);
my @Ctransition=(0.3,0.1,0.6);
my $transitionRef->{'A'}=\@Atransition;
$transitionRef->{'B'}=\@Btransition;
$transitionRef->{'C'}=\@Ctransition;
my @Atransmit=(0.1,0.2,0.3,0.4);
my @Btransmit=(0.2,0.1,0.2,0.5);
my @Ctransmit=(0.3,0.4,0.1,0.2);
my $transmitRef->{'A'}=\@Atransmit;
$transmitRef->{'B'}=\@Btransmit;
$transmitRef->{'C'}=\@Ctransmit;
my @obs=qw(1 4 2 3 1 4 1 1 4 3);
my $proRef; #前向向量
for(my $i=0;$i<scalar @obs;$i++){
my $obsKindNum=$obsKind2num->{$obs[$i]};
if($i ==0){
for(my $i=0;$i<scalar @state;$i++){
$proRef->{0}->{$state[$i]}=$initialP[$i]*$transmitRef->{$state[$i]}->[$obsKindNum];
}
next;
}
for(my $k=0;$k<scalar @state;$k++){
my $pro=0;
for(my $h=0;$h<scalar @state;$h++){
$pro+=$proRef->{$i-1}->{$state[$h]}*$transitionRef->{$state[$h]}->[$k];
}
$proRef->{$i}->{$state[$k]}=$pro*$transmitRef->{$state[$k]}->[$obsKindNum];
}
}
#把T时刻各个状态的概率相加得到观测序列的概率
my $proFinal=0;
foreach my$state(@state){
$proFinal+=$proRef->{(scalar @obs-1)}->{$state};
}
print "The final probability is: ",$proFinal,"\n";
对于问题2,同样是forward algorithm和backword algorithm都可以用。forward algorithm的实现如下：
#指定HMM模型的参数
my $stateNum=3;
my $obsNum=10;
my @initialP=qw(0.1 0.5 0.4);
my @state=qw(A B C);
my @obsKind=qw(1 2 3 4);
my $obsKind2num;
for(my $i=0;$i<scalar @obsKind;$i++){
$obsKind2num->{$obsKind[$i]}=$i;
}
my @Atransition=(0.8,0.1,0.1);
my @Btransition=(0.2,0.6,0.2);
my @Ctransition=(0.3,0.1,0.6);
my $transitionRef->{'A'}=\@Atransition;
$transitionRef->{'B'}=\@Btransition;
$transitionRef->{'C'}=\@Ctransition;
my @Atransmit=(0.1,0.2,0.3,0.4);
my @Btransmit=(0.2,0.1,0.2,0.5);
my @Ctransmit=(0.3,0.4,0.1,0.2);
my $transmitRef->{'A'}=\@Atransmit;
$transmitRef->{'B'}=\@Btransmit;
$transmitRef->{'C'}=\@Ctransmit;
my @obs=qw(1 1 4 2 2 3 3);
my $current2preMax; #用来回溯状态序列
my $proRef; #$proRef->{$index}->{$state} 前向向量
for(my $index=0;$index <scalar @obs;$index++){
my $obsNum=$obsKind2num->{$obs[$index]};
my $currentObs=$obs[$index];
if($index==0){
for(my $i=0;$i<scalar @state;$i++){
$proRef->{$index}->{$state[$i]}=$initialP[$i]*$transmitRef->{$state[$i]}->[$obsNum];
}
next;
}
for(my $i=0;$i<scalar @state;$i++){
my $subRef;
for(my $j=0;$j<scalar @state;$j++){
$subRef->{$state[$j]}=$proRef->{$index-1}->{$state[$j]}*$transitionRef->{$state[$j]}->[$i];
}
my @previousState=sort {$subRef->{$b}<=>$subRef->{$a}}keys %{$subRef};
$proRef->{$index}->{$state[$i]}=$subRef->{$previousState[0]}*$transmitRef->{$state[$i]}->[$obsNum];
$current2preMax->{$index}->{$state[$i]}=$previousState[0];
}
}
my @index=sort {$b<=>$a}keys %{$current2preMax};
my @finalStateSort=sort {$proRef->{$index[0]}->{$b} <=> $proRef->{$index[0]}->{$a}}keys %{$proRef->{$index[0]}};
print "The final step: ",$index[0],"\tthe state with the largest probiligy: ",$finalStateSort[0],"\n";
#回溯过程
my @stateSeq;
my $finalState=$finalStateSort[0];
push(@stateSeq,$finalState);
for(my $i=$index[0];$i>0;$i--){
my $cur=$stateSeq[-1];
my $next=$current2preMax->{$i}->{$cur};
push(@stateSeq,$next);
}
@stateSeq=reverse @stateSeq;
print join("-",@stateSeq),"\n";
print join("-",@obs),"\n";
对于问题3，采用forward-backward algorithm，也称为 Baum-Welch算法。它是一种贪婪算法，往往只能够达到局部最优。它的实现如下：
##############################################################################
#采用Forward-backward algorithm （Baum-Welch算法）估计HMM模型的参数；
##############################################################################
#模型参数的初始化
my $stateNum=3;
my $obsNum=10;
my @initialP=qw(0.1 0.5 0.4);
my @state=qw(A B C);
my @obsKind=qw(1 2 3 4);
my $obsKind2num;
for(my $i=0;$i<scalar @obsKind;$i++){
$obsKind2num->{$obsKind[$i]}=$i;
}
my @Atransition=(0.8,0.1,0.1);
my @Btransition=(0.2,0.6,0.2);
my @Ctransition=(0.3,0.1,0.6);
my $transitionRef->{'A'}=\@Atransition;
$transitionRef->{'B'}=\@Btransition;
$transitionRef->{'C'}=\@Ctransition;
my @Atransmit=(0.1,0.2,0.3,0.4);
my @Btransmit=(0.2,0.1,0.2,0.5);
my @Ctransmit=(0.3,0.4,0.1,0.2);
my $transmitRef->{'A'}=\@Atransmit;
$transmitRef->{'B'}=\@Btransmit;
$transmitRef->{'C'}=\@Ctransmit;
my @obs=qw(1 1 4 1 2 1 1 2 3 4 2 2 3 3 4 2); #观测序列
#初始模型下得到观测序列的概率;
my $proBefore=0;
my $proRef;
for(my $i=0;$i<scalar @obs;$i++){
my $obsKindNum=$obsKind2num->{$obs[$i]};
if($i ==0){
for(my $i=0;$i<scalar @state;$i++){
$proRef->{0}->{$state[$i]}=$initialP[$i]*$transmitRef->{$state[$i]}->[$obsKindNum];
}
next;
}
for(my $k=0;$k<scalar @state;$k++){
my $pro=0;
for(my $h=0;$h<scalar @state;$h++){
$pro+=$proRef->{$i-1}->{$state[$h]}*$transitionRef->{$state[$h]}->[$k];
}
$proRef->{$i}->{$state[$k]}=$pro*$transmitRef->{$state[$k]}->[$obsKindNum];
}
}
foreach my$state(keys %{$proRef->{(scalar @obs)-1}}){
$proBefore+=$proRef->{(scalar @obs)-1}->{$state};
}
$proBefore=log($proBefore);
open(OUT,">test")or die"$!";
my $proAfter=$proBefore;
my $threshold=1e-10; #判断迭代过程终止的cutoff
for(my $time=1;$time < 1000;$time++){
my $epsonRef; #$epsonRef->{$t}->{$i}->{$j}
my $alphaRef; #$alphaRef->{$t}->{$i}
my $betaRef; #$betaRef->{$t}->{$i}
#calculate $alphaRef;
for(my $t=0;$t <scalar @obs;$t++){
my $obsKindNum=$obsKind2num->{$obs[$t]};
if($t==0){
for(my $i=0;$i<scalar @state;$i++){
$alphaRef->{0}->{$state[$i]}=$initialP[$i]*$transmitRef->{$state[$i]}->[$obsKindNum];
}
next;
}
for(my $i=0;$i<scalar @state;$i++){
my $pro=0;
for(my $j=0;$j<scalar @state;$j++){
$pro+=$alphaRef->{$t-1}->{$state[$j]}*$transitionRef->{$state[$j]}->[$i];
}
$alphaRef->{$t}->{$state[$i]}=$pro*$transmitRef->{$state[$i]}->[$obsKindNum];
}
}
#calculate $betaRef;
for(my $t=(scalar @obs-1);$t >=0;$t--){
if($t==((scalar @obs)-1)){
for(my $i=0;$i<scalar @state;$i++){
$betaRef->{$t}->{$state[$i]}=1;
}
next;
}
my $obsKindNum=$obsKind2num->{$obs[$t+1]};
for(my $i=0;$i<scalar @state;$i++){
my $pro=0;
for(my $j=0;$j<scalar @state;$j++){
$pro+=$betaRef->{$t+1}->{$state[$j]}*$transitionRef->{$state[$i]}->[$j]*$transmitRef->{$state[$j]}->[$obsKindNum];
}
$betaRef->{$t}->{$state[$i]}=$pro;
}
}
#calculate $epsonRef;
for(my $t=0;$t <scalar @obs-1;$t++){
my $obsKindNum=$obsKind2num->{$obs[$t+1]};
my $denominator=0;
for(my $i=0;$i<scalar @state;$i++){
for(my $j=0;$j<scalar @state;$j++){
$denominator+=$alphaRef->{$t}->{$state[$i]}*$transitionRef->{$state[$i]}->[$j]*$transmitRef->{$state[$j]}->[$obsKindNum]*$betaRef->{$t+1}->{$state[$j]};
}
}
for(my $i=0;$i<scalar @state;$i++){
for(my $j=0;$j<scalar @state;$j++){
my $numerator=0;
$numerator+=$alphaRef->{$t}->{$state[$i]}*$transitionRef->{$state[$i]}->[$j]*$transmitRef->{$state[$j]}->[$obsKindNum]*$betaRef->{$t+1}->{$state[$j]};
$epsonRef->{$t}->{$state[$i]}->{$state[$j]}=$numerator/$denominator;
}
}
}
#calculate $rRef;
my $rRef;
for(my $t=0;$t<scalar @obs;$t++){
for(my $i=0;$i<scalar @state;$i++){
my $r=0;
foreach my$state(keys %{$epsonRef->{$t}->{$state[$i]}}){
$r+=$epsonRef->{$t}->{$state[$i]}->{$state};
}
$rRef->{$t}->{$state[$i]}=$r;
}
}
#Re-estimation of the parameters;
for(my $i=0;$i<scalar @state;$i++){
$initialP[$i]=$rRef->{0}->{$state[$i]};
}
for(my $i=0;$i<scalar @state;$i++){
my $denominator=0;
for(my $t=0;$t<(scalar @obs)-1;$t++){
$denominator+=$rRef->{$t}->{$state[$i]};
}
for(my $j=0;$j<scalar @state;$j++){
my $numerator=0;
for(my $t=0;$t < (scalar @obs)-1;$t++){
$numerator+=$epsonRef->{$t}->{$state[$i]} ->{$state[$j]};
}
$transitionRef->{$state[$i]}->[$j]=$numerator/$denominator;
}
}
for(my $i=0;$i<scalar @state;$i++){
my $denominator=0;
for(my $t=0;$t<(scalar @obs);$t++){
$denominator+=$rRef->{$t}->{$state[$i]};
}
for(my $j=0;$j<scalar @obsKind;$j++){
my $numerator=0;
for(my $t=0;$t<scalar @obs;$t++){
if($obs[$t]==$obsKind[$j]){
$numerator+=$rRef->{$t}->{$state[$i]};
}
}
$transmitRef->{$state[$i]}->[$j]=$numerator/$denominator;
}
}
#计算更新参数后的HMM模型得到观测序列的概率；
my $proRef;
$proAfter=0;
for(my $i=0;$i<scalar @obs;$i++){
my $obsKindNum=$obsKind2num->{$obs[$i]};
if($i ==0){
for(my $i=0;$i<scalar @state;$i++){
$proRef->{0}->{$state[$i]}=$initialP[$i]*$transmitRef->{$state[$i]}->[$obsKindNum];
}
next;
}
for(my $k=0;$k<scalar @state;$k++){
my $pro=0;
for(my $h=0;$h<scalar @state;$h++){
$pro+=$proRef->{$i-1}->{$state[$h]}*$transitionRef->{$state[$h]}->[$k];
}
$proRef->{$i}->{$state[$k]}=$pro*$transmitRef->{$state[$k]}->[$obsKindNum];
}
}
foreach my$state(keys %{$proRef->{(scalar @obs)-1}}){
$proAfter+=$proRef->{(scalar @obs)-1}->{$state};
}
$proAfter=log($proAfter);
my $delta=$proAfter-$proBefore;
$proBefore=$proAfter;
#输出更新后的HMM参数
print OUT $time,"\t",$delta,"\n";
print OUT "initial: ",join("\t",@initialP),"\n";
print OUT "transition matrix: \n";
print OUT "\t",join("\t",@state),"\n";
for(my $i=0;$i<scalar @state;$i++){
print OUT $state[$i];
for(my $j=0;$j<scalar @state;$j++){
print OUT "\t",$transitionRef->{$state[$i]}->[$j];
}
print OUT "\n";
}
print OUT "transmit matrix: \n";
print OUT "\t",join("\t",@obsKind),"\n";
for(my $i=0;$i<scalar @state;$i++){
print OUT $state[$i];
for(my $j=0;$j<scalar @obsKind;$j++){
print OUT "\t",$transmitRef->{$state[$i]}->[$j];
}
print OUT "\n";
}
print OUT "\n";
if($delta < $threshold){
print OUT $time,"\t",$delta,"\tOver\n";
last;
}
}
close OUT or die"$!";
更多语言版本的实现见网友的总结：http://hi.baidu.com/shinchen2/item/86fdadcd67d63c7189ad9ee3。此外，我在学习的过程中参考了下面这位网友的学习经验(http://www.52nlp.cn/hmm-learn-best-practices-seven-forward-backward-algorithm-5).他写得很通俗易懂，推荐大家看看。
另外，在R语言中专门有一个包HMM用来实现该模型。
https://blog.sciencenet.cn/blog-54276-648486.html
上一篇：R中用于排序的几个函数
下一篇：perl中的pos()函数处理重复匹配的问题
