HMM学习最佳范例七：前向-后向算法5

[…] 隐马尔科夫模型(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],”tthestatewiththelargestprobiligy:”,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 mystate(keysproRef->{(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 mystate(keysepsonRef->{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 mystate(keysproRef->{(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,”tOvern”; last; } } close OUT or die”$!”; 更多语言版本的实现见网友的总结：http://hi.baidu.com/shinchen2/item/86fdadcd67d63c7189ad9ee3。此外，我在学习的过程中参考了下面这位网友的学习经验(https://www.52nlp.cn/hmm-learn-best-practices-seven-forward-backward-algorithm-5).他写得很通俗易懂，推荐大家看看。 另外，在R语言中专门有一个包HMM用来实现该模型。 […]