背景
承接上文,繼續 HMM 的公式推導過程,本文主要介紹後向機率演算法以及前向機率演算法的 Java 實現版本。
後向機率公式及推導過程1.後向機率定義給定 λ,定義到時刻 t ,狀態為 si 的前提下,從 t+1 至 T 時刻出現的觀測序列為 qt+1,qt+2… qT 的機率。記作: βt(i) = P(qt+1,qt+2,…qT|st=si)。
2.後向機率的初始條件根據定義,第 T 個時刻的後向機率是 βT(i)=P(qT+1,sT=si|)=1 。
βT(i) = 1 是指 T+1 時刻觀測序列是指定值的機率,既然 T 時刻已經是最後時刻了,那麼 T+1 時刻就是必然發生的事件了,必然發生的事件,其機率值是 1。
3.後向機率的遞推公式向後遞推,就是從第 T 個時刻 βT(i) 值,向後推導,計算出 β1(i) ,再將第一個時刻的所有狀態累加就得到 HMM 計算問題的解。
遞推的核心是根據 t 時刻的後向機率的定義,引入 t+1 時刻的狀態,再根據聯合機率和全機率計算公式,得到 t 時刻和 t-1 時刻的後向機率之間的遞推關係。
4.後向機率遞推過程這是我在簡書上看到的一個遞推計算公式,HMM 計算問題 原文只有黑色的公式,沒有每一步的縮減依據。這裡補充我的推導過程,有助於深刻理解演算法本身。
我們看下推導過程:
第 1 步,是公式的基本定義;第 2 步,根據全機率公式,引入 t+1 時刻的狀態事件,將全機率計算轉換為求解獨立事件 st+1 的機率和。3.第 2 步和第 3 步之前的中間步驟,是聯合條件的機率公式,P(AB|C)=P(A|BC)*P(B|C) ,經過該公式轉換後得到中間步驟為 P(qt+1...qT|it+1=Sj,it=Si)*P(it+1=Sj|it=Si)。第 3 步,是根據 HMM 的獨立觀測假設,t+1 時刻的觀測值只與 t+1 時刻的狀態無關,與其他時刻的狀態如 st=Si 無關,所以可以消掉左側機率中的 it=si 這個條件。即 P(qt+1...qT|it+1=Sj,it=Si) = P(qt+1...qT|it+1=Sj)第 3 步和第 4 步,之間省略了一個推導過程是聯合條件機率的轉換公式,中間轉換的結果為:P(qt+2...qT|qt+1,it+1=Sj)*P(qt+1|it+1=Sj)第 4 步是根據 HMM 的獨立觀測假設,消掉左側聯合條件機率中的 qt+1 這個條件得到的簡化值,即 P(qt+2...qT|qt+1,it+1=Sj) = P(qt+2...qT|it+1=Sj)第 5 步,再次根據後向機率的定義,以及 HMM 的模型引數的關係,最終得到 t 時刻和 t-1 時候的遞推關係。這樣的話,後向機率演算法就轉換為一個初始值和一個遞推公式。
前向機率演算法從上一篇介紹的前向機率的四個基本資訊已經構成了一個完整的演算法,用陣列很容易實現。演算法的基本資訊為:
初始條件,四個陣列 π、 A、B、Q;初始值,α1(i) = πi * Bi(q1);遞推公式,αt(i) = Sum(αt-1(j) * A[j][[i])* B[i][qt];目標函式 ,sum(αT(i)),i=1,2....n;這裡來實現經典的取球案例的機率,案例描述:
有三個盒子,每個盒子有紅球白球,從三個盒子中隨機取一個盒子的機率是 πi,三個盒子狀態轉換機率為 A ,從盒子中取出對應球的觀測轉移機率為 B 。目標:抽取三次,結果為:紅白紅,求其出現的機率。
Java 實現前向機率 αt(i) 是一個元素,儲存 t 時刻狀態為 i 、目標觀測值 q1,q2…qt 出現的機率,所以需要定義一個二維陣列,大小為 T * N (N 為狀態的總數)。
實現程式碼為:
import java.util.Arrays;/*** @Title:HMMForward * @Description: 馬爾可夫模型的基本實現過程<br/>* 1、初始模型引數 π,A,B,最終觀測序列 O=[q1,q2,q3....qT]* 2、初始值 α[1][i]= π[i] * B[i][q1]* 3、迴圈 t=1....T ,計算 α[t][i]= (∑α[t-1][j]*a[j][i])*b[i][qt]* 4、目標 P(Q|λ)= ∑α[T][i] ,i=1....M* @author:wang_ll* @date :2019年7月29日 */public class HMMForward { /** * 狀態個數 */ private int M = 0; /** * 狀態初始的機率 */ private double PI [] = {0.2,0.4,0.4}; /** * 狀態轉移矩陣 M X M */ private double A [][] = {{0.5,0.2,0.3},{0.3,0.5,0.2},{0.2,0.3,0.5}}; /** * 觀察指轉移矩陣 */ private double B[][] = {{0.5,0.5},{0.4,0.6},{0.7,0.3}}; /** * 最終求解的觀測序列值 */ private int [] observeSequence = {0,1,0}; /** * 儲存最終計算的前向機率的值的陣列 T X M */ private double forwardRate [][] = null; /** * 最終的計算目標:P(Q|λ)= ∑α[T][i] ,i=1....M */ private double target ; /** * 初始化模型引數 */ public void initModel() { M = 3; forwardRate = new double[observeSequence.length][M]; } /** * HMM 計算問題求解 */ public void calculate() { //初始化初值 α[1][i]= π[i] * B[i][q1] for(int i=0;i<M;i++) { forwardRate[0][i] = PI[i]*B[i][observeSequence[0]]; } //向前遞推 for(int t=1;t<observeSequence.length;t++) { for(int i=0;i<M;i++) { double sumAllState = 0.0; for(int j=0;j<M;j++) { //累加各個狀態 t-1 時刻的機率和 sumAllState+= forwardRate[t-1][j]*A[j][i]; } forwardRate[t][i] = sumAllState* B[i][observeSequence[t]]; } } //計算最終的結果 for(int i=0;i<M;i++) { target+=forwardRate[forwardRate.length-1][i]; } } /** * 列印最終的結果 */ public void printResult() { System.out.println("初始模型引數 π:"+Arrays.toString(PI)); System.out.println("狀態轉移矩陣 A 為:"); for(int i=0;i<M;i++) { System.out.println(Arrays.toString(A[i])); } System.out.println("觀測值轉移矩陣 B 為:"); for(int i=0;i<M;i++) { System.out.println(Arrays.toString(B[i])); } System.out.println("最終觀測序列 :"+Arrays.toString(observeSequence)); System.out.println(); System.out.println("HMM 計算得到的機率值為:"+target); } public static void main(String[] args) { HMMForward instance = new HMMForward(); instance.initModel(); instance.calculate(); instance.printResult(); }}
執行結果:
初始模型引數 π:[0.2, 0.4, 0.4]狀態轉移矩陣 A 為:[0.5, 0.2, 0.3][0.3, 0.5, 0.2][0.2, 0.3, 0.5]觀測值轉移矩陣 B 為:[0.5, 0.5][0.4, 0.6][0.7, 0.3]最終觀測序列值:[0, 1, 0]最終觀測序列描述:[紅,白,紅]HMM 計算得到的機率值為:0.130218
參考附錄
(一)HMM 模型計算問題完整的推導說明(二)HMM Java 實現的例子(三)HMM 解碼演算法實現
這三篇文章,對我理解 HMM 的過程有很大的幫助,《HMM Java 實現的例子》這篇文章中的 HMM 演算法實現,如果不瞭解演算法的推導過程直接看程式碼是不好懂的;理解過程後,自己實現的程式碼結果跟它是一致的,也是一種印證。
《HMM Java 實現的例子》這篇文章給我的啟示是:原來 HMM 用 Java 程式碼實現是這麼簡單。最初的教程使用 python 的 hmmlearn 庫,但是沒有安裝成功,現在看來用 Java 和 python 是殊途同歸,自己實現一遍後發現 python 的庫顯得沒那麼高深了!