簡介

隱馬爾可夫模型（Hidden Markov Model，HMM）創(chuàng)立于20世紀70年代。主要用于行為識別，語音識別，文字識別等。

原理簡述

隱馬爾可夫模型由五個部分組成：狀態(tài)空間S，觀測空間O，初始狀態(tài)概率空間PI，狀態(tài)概率轉(zhuǎn)移矩陣P以及觀測值生成概率矩陣Q。另外，隱馬爾可夫模型還包括一條觀測鏈，一條隱藏鏈。（后面將詳述）下面是隱馬爾可夫模型示意圖：

因此整個過程就是觀測值隨狀態(tài)的轉(zhuǎn)移而生成，而我們所關(guān)心的是通過已有的觀測值來判斷其隱藏的狀態(tài)，即通過一長串的觀測序列推算導(dǎo)致這一結(jié)果的可能的狀態(tài)序列。例如：有兩枚不同的硬幣（一枚正面拋擲后正面朝上的概率比較大，另一個反面朝上的概率比較大），現(xiàn)在一個人按照其習(xí)慣每次選擇其中的一枚硬幣拋擲，共拋擲N次，將結(jié)果記錄下來（設(shè)正面為1，反面為0），之后你就可以利用隱馬爾可夫模型，通過已有結(jié)果反推這個人每次使用哪枚硬幣進行投擲的。

實現(xiàn)方法

要實現(xiàn)上面所述原理就必須解決三個問題：評估問題（evaluation），解碼問題（decoding）和學(xué)習(xí)問題（learning）

1、評估問題，即評估當(dāng)前狀態(tài)為真實狀態(tài)的可能性。

最簡單的方法有前向算法和后向算法（當(dāng)然也可以聯(lián)合使用這兩種算法）。

前向算法：從前遞歸，一層一層計算概率，最后再求總和。

1）t=0（事實上 t 的首項應(yīng)該為1，但是考慮到編程的方便這里就設(shè)首項為0）

，即?alpha(i,t)=PI(i)*Q(i,t)（偽代碼，這只是為了表示方便易懂，與之后的代碼可能會有出入）

PS：alpha(i,t)指t時刻狀態(tài)為Si的概率（下面同義），PI(i)為狀態(tài)Si的初始概率，Q(i,t)指的是 t 時刻觀測值Vt由狀態(tài)Si生成的概率。

2）t>0 && t<=n

，即alpha(t,i)=Sum[ alpha(j,t-1)*P(j,i)*Q(i,t) ]

PS：P(j,i)指由狀態(tài) Sj 轉(zhuǎn)移到 Si 的概率

3）

，即將3）所算的所有狀態(tài)Si的結(jié)果再求和。

PS：對應(yīng)后面的Java類為?AlgorithmFront.java

后向算法：與前向算法相反

1）t>=0 && t<n

，即beta(i,t)=Sum[ beta(j,t+1)*P(i,j)*Q(j,t+1) ]

PS：beta(i,t)表示t時刻狀態(tài)為Si的條件下，從t+1時刻到n生成相應(yīng)觀測序列的概率。

2）t=n

PS：因為下一個時刻就已結(jié)束，所以無論是什么狀態(tài)都是確定的，所以概率都為1。

3）

，與前向算法相似，最后也是將所有結(jié)果進行求和。

PS：對應(yīng)后面的Java類為?AlgorithmBack.java

2、解碼問題，即如何根據(jù)觀測值，狀態(tài)轉(zhuǎn)移概率矩陣，生成概率矩陣得到真正的狀態(tài)序列。（有時候你完全可以根據(jù)先驗知識給參數(shù)設(shè)值，這樣就無需 學(xué)習(xí)步驟（Learning） 便可以解碼了）

Viterbi算法：基本原理就是計算概率每一步最高時對應(yīng)的狀態(tài)序列

1）初始化

2）遞歸

3）終止

，

PS：表示 n 時刻沿著X1,X2,...Xn 且在 n 時刻狀態(tài)Xn=Si 產(chǎn)生相應(yīng)觀測序列的最大概率

保存著狀態(tài)序列信息。

4）回溯

根據(jù)??的結(jié)果便可知道相應(yīng)的狀態(tài)序列了。

PS：對應(yīng)的后面的Java類為?HMMDecisionVbImp.java

3、學(xué)習(xí)問題，即如何通過觀測值來獲取初始狀態(tài)概率，狀態(tài)轉(zhuǎn)移概率矩陣以及生成概率矩陣。

Baum-Welch算法：

Step1: 隨機產(chǎn)生一組參數(shù)，并代入評估函數(shù)（evaluation，例如前向算法），計算結(jié)果。

Step2: 利用參數(shù)估算初始狀態(tài)概率，狀態(tài)轉(zhuǎn)移概率矩陣以及生成概率矩陣

由于：

，即kis(i,j,t)=alpha(i,t)*P(i,j)*Q(j,t+1)*beta(j,t+1)

，即gamma(i,t)=alpha(i,t)*beta(i,t)

PS：a）偽代碼中并沒有除以，這主要是為了減少運算量，因為之后計算狀態(tài)概率矩陣、生成矩陣這項都會被約掉。

b）kis(i,j,t)即，表示t時刻為狀態(tài)Si，t+1時刻為狀態(tài)Sj的概率

c）gamma(i,t)即，表示t時刻狀態(tài)為Si的概率

d）相應(yīng)的Java類為?Gammas.java，Ksis.java

所以：

1）估計概率轉(zhuǎn)移概率矩陣

2）估計初始狀態(tài)概率

，（注意，實際編程實現(xiàn)時這里還需除以之前漏除的）

3）估計概率生成矩陣

Step3: 將剛估計的參數(shù)代入 評估函數(shù) 進行計算，并與上一次評估的結(jié)果做比較，若差異小于某個閾值（thresh，例如 0.05）則接受。否則繼續(xù)迭代計算。

PS：相應(yīng)的Java類為?HMMLearnBwImp.java

Java具體實現(xiàn)

基礎(chǔ)類 Package lxwo.utils

1、AlgorithmFront

package lxwo.utils;
public class AlgorithmFront{
private double[] Api;
private double[][] AP;
private double[][] AQ;
private int[] V;
public AlgorithmFront(double[] Api, double[][] AP, double[][] AQ, int[] V) {
this.Api = Api;
this.AP = AP;
this.AQ = AQ;
this.V = V;
}
public double calculate(int step) {
double Result = 0.0;
for (int Pindex = 0; Pindex < this.AP.length; Pindex++)
Result += this.alpha(Pindex, step);
return Result;
}
public double alpha(int toI, int step) {
double tempValue = 0.0;
if (step > 0) {
for (int pindex = 0; pindex < this.AP.length; pindex++)
tempValue += this.alpha(pindex, step - 1)* this.AP[pindex][toI] * this.AQ[toI][this.V[step]];
} else?
tempValue = Api[toI] * this.AQ[toI][this.V[step]];

return tempValue;
}
public double[] getApi() {
return Api;
}
public void setApi(double[] api) {
Api = api;
}
public double[][] getAP() {
return AP;
}
public void setAP(double[][] aP) {
AP = aP;
}
public double[][] getAQ() {
return AQ;
}
public void setAQ(double[][] aQ) {
AQ = aQ;
}
public int[] getV() {
return V;
}
public void setV(int[] v) {
V = v;
}
}

2、AlgorithmBack

package lxwo.utils;
public class AlgorithmBack{

private double[] Api;
private double[][] AP;
private double[][] AQ;
private int[] V;

public AlgorithmBack(double[] Api, double[][] AP, double[][] AQ, int[] V){
this.Api = Api;
this.AP = AP;
this.AQ = AQ;
this.V = V;
}

public double calculate(int step){
double Result = 0.0;

for (int Pindex = 0; Pindex < this.AP.length; Pindex++)
Result += this.belta(Pindex, step);

return Result;
}

public double belta(int fromI,int step){
double tempValue = 0.0;

if(step<this.V.length-1){
for(int pindex=0;pindex<this.AP.length;pindex++)
tempValue += this.belta(pindex, step+1)*this.AP[fromI][pindex]*this.AQ[pindex][this.V[step+1]];
}else{
tempValue = 1.0;
}

return tempValue;
}


public double[] getApi() {
return Api;
}


public void setApi(double[] api) {
Api = api;
}


public double[][] getAP() {
return AP;
}


public void setAP(double[][] aP) {
AP = aP;
}


public double[][] getAQ() {
return AQ;
}


public void setAQ(double[][] aQ) {
AQ = aQ;
}


public int[] getV() {
return V;
}


public void setV(int[] v) {
V = v;
}


}

3、?Ksis?

package lxwo.utils;
public class Ksis {

private double[] Api;
private double[][] AP;
private double[][] AQ;
private int[] V;

public Ksis(double[] Api, double[][] AP, double[][] AQ, int[] V){
this.Api = Api;
this.AP = AP;
this.AQ = AQ;
this.V = V;
}

public double calculate(int i,int j,int step){
AlgorithmFront f1 = new AlgorithmFront(this.Api,this.AP,this.AQ,this.V);
AlgorithmBack f2 = new AlgorithmBack(this.Api,this.AP,this.AQ,this.V);
// Considering the amount of calculation, we don't divide the result by p(V|lambda)
return f1.alpha(i, step)*this.AP[i][j]*this.AQ[j][this.V[step+1]]*f2.belta(j, step+1);
}

public double sumKsi(int i,int j, int T){
double tempValue = 0.0;
for(int pindex=0;pindex<T;pindex++)
tempValue += this.calculate(i, j, pindex);
return tempValue;
}


}

4、Gammas?

package lxwo.utils;
public class Gammas {


private double[] Api;
private double[][] AP;
private double[][] AQ;
private int[] V;

public Gammas(double[] Api, double[][] AP, double[][] AQ, int[] V){
this.Api = Api;
this.AP = AP;
this.AQ = AQ;
this.V = V;
}

public double calculate(int i,int step){
AlgorithmFront f1 = new AlgorithmFront(this.Api,this.AP,this.AQ,this.V);
AlgorithmBack f2 = new AlgorithmBack(this.Api,this.AP,this.AQ,this.V);
// Considering the amount of calculation, we don't divide the result by p(V|lambda)
return f1.alpha(i, step)*f2.belta(i, step);//step+1
}

public double sumGamma(int i, int T){
double tempValue = 0.0;
for(int pindex=0;pindex<T;pindex++)
tempValue += this.calculate(i, pindex);
return tempValue;
}

}

核心類 Package lxwo.core

1、HMMDecision &??HMMDecisionVbImp

package lxwo.core;
public interface HMMDecision {
public int[] recognize(int step);
}

package lxwo.core;
public class HMMDecisionVbImp implements HMMDecision {
private double[] pi;
private double[][] P;
private double[][] Q;
private int[] V;
private int[]phi;
public HMMDecisionVbImp(double[] pi, double[][] P, double[][] Q, int[] V) {
this.pi = pi;
this.P = P;
this.Q = Q;
this.V = V;
this.phi = new int[this.V.length];
for(int i=0;i<this.phi.length;i++)
this.phi[i]=-1;
}
@Override
public int[] recognize(int step) {
int[] tempFlag = new int[this.phi.length];
double sumTempMax = 0.0;
for(int dindex=0;dindex<this.P.length;dindex++){
double tempVal = this.delta(dindex, step);
if(tempVal>sumTempMax){
sumTempMax = tempVal;
tempFlag = this.phi.clone();
tempFlag[step]=dindex;
}
}
return tempFlag;
}
private double delta(int toI,int step) {
double tempValue = 1.0;
if (step == 0) {
tempValue = this.pi[toI]*this.Q[toI][step];
} else {
double tempMax = 0.0;
for(int jindex=0;jindex<this.P.length;jindex++){
double tempV = delta(jindex,step-1)*this.P[jindex][toI];
if(tempV>tempMax){
tempMax = tempV;
this.phi[step-1]=jindex;
}
}
tempValue = tempMax*this.Q[toI][this.V[step]];
}
return tempValue;
}
}

2、HMMLearn&?HMMLearnBwImp

package lxwo.core;
public interface HMMLearn {
public boolean learn();
}

package lxwo.core;
import lxwo.utils.AlgorithmFront;
import lxwo.utils.Gammas;
import lxwo.utils.Ksis;


public class HMMLearnBwImp implements HMMLearn {


private double[] pi;
private double[][] P;
private double[][] Q;
private int[] V;
private double thresh;
private int deadline;


public HMMLearnBwImp(double[] pi, double[][] P, double[][] Q, int[] V,
double thresh, int deadline) {
this.pi = pi;
this.P = P;
this.Q = Q;
this.V = V;
this.thresh = thresh;
this.deadline = deadline;
}


@Override
public boolean learn() {
// TODO Auto-generated method stub
double flag1 = 0.0;
double flag2 = 0.0;
double flag3 = 0.0;
double[] tpi = new double[this.pi.length];
double[][] tP = new double[this.P.length][this.P[0].length];
double[][] tQ = new double[this.Q.length][this.Q[0].length];
int count = 0;
double diff = 1000.0;
flag3 = new AlgorithmFront(pi, P, Q, V).calculate(this.V.length - 1);
do {
count++;
// evaluate
flag1 = flag3;
// recalculate pi
double tempM1 = new AlgorithmFront(this.pi, this.P, this.Q, this.V)
.calculate(this.V.length - 1);
for (int i1 = 0; i1 < tpi.length; i1++)
tpi[i1] = (new Gammas(this.pi, this.P, this.Q, this.V)
.calculate(i1, 0)) / tempM1; // 'cause we don't divide?it before, so we?should make up here

// recalculate P
for (int i2 = 0; i2 < this.P.length; i2++)
for (int j2 = 0; j2 < this.P[0].length; j2++)
tP[i2][j2] = (new Ksis(this.pi, this.P, this.Q, this.V)
.sumKsi(i2, j2, this.V.length - 1))
/ (new Gammas(this.pi, this.P, this.Q, this.V)
.sumGamma(i2, this.V.length - 1));


// recalculate Q
for (int i3 = 0; i3 < this.Q.length; i3++) {
double tempM2 = new Gammas(this.pi, this.P, this.Q, this.V)
.sumGamma(i3, this.V.length);
for (int j3 = 0; j3 < this.V.length; j3++)
tQ[i3][this.V[j3]] += (new Gammas(this.pi, this.P, this.Q,
this.V).calculate(i3, j3)) / tempM2;
}

// re-evaluate
flag2 = new AlgorithmFront(tpi, tP, tQ, V)
.calculate(this.V.length - 1);
flag3 = flag2;

// reset args
this.pi = tpi.clone();
this.P = tP.clone();
this.Q = tQ.clone();
tQ = new double[this.Q.length][this.Q[0].length];

diff = Math.abs(flag1 - flag2);

} while (diff > thresh && count < this.deadline);

System.out.println("count: "+count);

if (count == this.deadline && diff > this.thresh)
return false;
else
return true;
}


public double[] getPi() {
return pi;
}


public void setPi(double[] pi) {
this.pi = pi;
}


public double[][] getP() {
return P;
}


public void setP(double[][] p) {
P = p;
}


public double[][] getQ() {
return Q;
}


public void setQ(double[][] q) {
Q = q;
}


}

測試類 Package lxwo.test

package lxwo.test;
import lxwo.core.HMMDecision;
import lxwo.core.HMMDecisionVbImp;
import lxwo.core.HMMLearnBwImp;
public class Test {
/**
* @param args

?* 實驗：用兩種骰子（0，1）投擲，其中一個骰子為正常的（0），另一個為灌鉛（1），出現(xiàn)456的可能性較大。

?* 代碼中用012345代替123456
*/
public static void main(String[] args) {
double[] api = { 0.5, 0.5 };
double[][] P = { { 0.9, 0.1 }, { 0.2, 0.8 }};
double[][] Q = { { 0.2, 0.16, 0.16, 0.16, 0.16, 0.16 }, {0, 0, 0.10, 0.30, 0.30, 0.30 } };
int[] V = {5,1,2,4,5,4,2,1,0,5};// -5,1,2,-4,-5,-4,2,1,0,5 這里標(biāo)記符號的表示用第二種骰子投擲的

HMMLearnBwImp hlbi = new HMMLearnBwImp(api, P, Q, V, 0.05, 100);
if (hlbi.learn()) {
System.out.println("result:");
HMMDecision hd = new HMMDecisionVbImp(hlbi.getPi(), hlbi.getP(), hlbi.getQ(), V);
int[] result = hd.recognize(V.length-1);
for(int r:result)
System.out.print(r+"\t");
System.out.println();
} else {
System.out.println("Fail!");
}
}
}

觀測序列：5，1，2，4，5，4，2，1，0，5

結(jié)果：0，0，0，1，1，1，0，0，0，0

除了第一項估計有誤，其余都正確。（這里只是用一條觀測值序列做的測試，如果有多條觀測值，預(yù)測結(jié)果會好很多）

PS：由于這只是簡單的實現(xiàn)HMM，因此其真正的實用性還不強（如果你把觀測值加到>20個，其運算時間將是巨大的，因此實際應(yīng)用中還需對上述代碼進行改進）

創(chuàng)作挑戰(zhàn)賽新人創(chuàng)作獎勵來咯，堅持創(chuàng)作打卡瓜分現(xiàn)金大獎

總結(jié)

以上是生活随笔為你收集整理的隐马尔可夫模型 HMM 原理及实现的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。