算法概述
隐马尔可夫模型(Hidden Markov Model, HMM)是一种统计模型,用于描述由隐藏的马尔可夫链随机生成观测序列的过程。HMM假设系统在每个时刻处于某个隐藏状态,隐藏状态之间的转移服从马尔可夫性质,而每个隐藏状态会以一定概率产生观测值。
HMM由三个核心要素定义:
– 初始概率分布pi:系统在t=1时处于各状态的概率
– 状态转移矩阵A:隐藏状态之间的转移概率
– 观测概率矩阵B:各状态下生成不同观测值的概率
HMM三大经典问题
- 评估问题:给定模型和观测序列O,计算P(O|lambda) —— 前向算法解决
- 解码问题:给定模型和观测序列O,找出最可能的隐藏状态序列 —— 维特比算法解决
- 学习问题:给定观测序列O,估计模型参数 —— Baum-Welch算法解决
完整Java实现
import java.util.*;
/**
* 隐马尔可夫模型(HMM)完整实现
* 包含前向算法、维特比解码和Baum-Welch学习
*/
public class HiddenMarkovModel {
int numStates;
int numObservations;
double[] pi;
double[][] A;
double[][] B;
public HiddenMarkovModel(int numStates, int numObservations) {
this.numStates = numStates;
this.numObservations = numObservations;
this.pi = new double[numStates];
this.A = new double[numStates][numStates];
this.B = new double[numStates][numObservations];
initializeRandom();
}
void initializeRandom() {
Random rand = new Random(42);
double sum = 0;
for (int i = 0; i < numStates; i++) {
pi[i] = rand.nextDouble();
sum += pi[i];
}
for (int i = 0; i < numStates; i++) pi[i] /= sum;
for (int i = 0; i < numStates; i++) {
sum = 0;
for (int j = 0; j < numStates; j++) {
A[i][j] = rand.nextDouble();
sum += A[i][j];
}
for (int j = 0; j < numStates; j++) A[i][j] /= sum;
}
for (int i = 0; i < numStates; i++) {
sum = 0;
for (int j = 0; j < numObservations; j++) {
B[i][j] = rand.nextDouble();
sum += B[i][j];
}
for (int j = 0; j < numObservations; j++) B[i][j] /= sum;
}
}
public double forward(int[] observations) {
int T = observations.length;
double[][] alpha = new double[T][numStates];
for (int i = 0; i < numStates; i++) {
alpha[0][i] = pi[i] * B[i][observations[0]];
}
for (int t = 1; t < T; t++) {
for (int j = 0; j < numStates; j++) {
double sum = 0;
for (int i = 0; i < numStates; i++) {
sum += alpha[t - 1][i] * A[i][j];
}
alpha[t][j] = sum * B[j][observations[t]];
}
}
double probability = 0;
for (int i = 0; i < numStates; i++) {
probability += alpha[T - 1][i];
}
return probability;
}
public double forwardLog(int[] observations) {
int T = observations.length;
double[][] logAlpha = new double[T][numStates];
for (int i = 0; i < numStates; i++) {
logAlpha[0][i] = Math.log(pi[i]) + Math.log(B[i][observations[0]]);
}
for (int t = 1; t < T; t++) {
for (int j = 0; j < numStates; j++) {
double[] terms = new double[numStates];
for (int i = 0; i < numStates; i++) {
terms[i] = logAlpha[t - 1][i] + Math.log(A[i][j]);
}
logAlpha[t][j] = logSumExp(terms) + Math.log(B[j][observations[t]]);
}
}
return logSumExp(logAlpha[T - 1]);
}
double logSumExp(double[] values) {
double max = values[0];
for (double v : values) if (v > max) max = v;
double sum = 0;
for (double v : values) sum += Math.exp(v - max);
return max + Math.log(sum);
}
public int[] viterbi(int[] observations) {
int T = observations.length;
double[][] delta = new double[T][numStates];
int[][] psi = new int[T][numStates];
for (int i = 0; i < numStates; i++) {
delta[0][i] = Math.log(pi[i]) + Math.log(B[i][observations[0]]);
psi[0][i] = 0;
}
for (int t = 1; t < T; t++) {
for (int j = 0; j < numStates; j++) {
double maxProb = Double.NEGATIVE_INFINITY;
int maxState = 0;
for (int i = 0; i < numStates; i++) {
double prob = delta[t - 1][i] + Math.log(A[i][j]);
if (prob > maxProb) {
maxProb = prob;
maxState = i;
}
}
delta[t][j] = maxProb + Math.log(B[j][observations[t]]);
psi[t][j] = maxState;
}
}
double maxFinalProb = Double.NEGATIVE_INFINITY;
int finalState = 0;
for (int i = 0; i < numStates; i++) {
if (delta[T - 1][i] > maxFinalProb) {
maxFinalProb = delta[T - 1][i];
finalState = i;
}
}
int[] path = new int[T];
path[T - 1] = finalState;
for (int t = T - 2; t >= 0; t--) {
path[t] = psi[t + 1][path[t + 1]];
}
return path;
}
public void baumWelch(int[][] observationSequences, int maxIterations) {
for (int iter = 0; iter < maxIterations; iter++) {
double[] newPi = new double[numStates];
double[][] newA = new double[numStates][numStates];
double[][] newB = new double[numStates][numObservations];
for (int[] observations : observationSequences) {
int T = observations.length;
double[][] alpha = new double[T][numStates];
for (int i = 0; i < numStates; i++) {
alpha[0][i] = pi[i] * B[i][observations[0]];
}
for (int t = 1; t < T; t++) {
for (int j = 0; j < numStates; j++) {
double sum = 0;
for (int i = 0; i < numStates; i++) {
sum += alpha[t - 1][i] * A[i][j];
}
alpha[t][j] = sum * B[j][observations[t]];
}
}
double[][] beta = new double[T][numStates];
for (int i = 0; i < numStates; i++) {
beta[T - 1][i] = 1.0;
}
for (int t = T - 2; t >= 0; t--) {
for (int i = 0; i < numStates; i++) {
double sum = 0;
for (int j = 0; j < numStates; j++) {
sum += A[i][j] * B[j][observations[t + 1]] * beta[t + 1][j];
}
beta[t][i] = sum;
}
}
double[][] gamma = new double[T][numStates];
double[][][] xi = new double[T - 1][numStates][numStates];
for (int t = 0; t < T; t++) {
double denominator = 0;
for (int i = 0; i < numStates; i++) {
denominator += alpha[t][i] * beta[t][i];
}
for (int i = 0; i < numStates; i++) {
gamma[t][i] = alpha[t][i] * beta[t][i] / denominator;
}
}
for (int t = 0; t < T - 1; t++) {
double denominator = 0;
for (int i = 0; i < numStates; i++) {
for (int j = 0; j < numStates; j++) {
denominator += alpha[t][i] * A[i][j] *
B[j][observations[t + 1]] * beta[t + 1][j];
}
}
for (int i = 0; i < numStates; i++) {
for (int j = 0; j < numStates; j++) {
xi[t][i][j] = alpha[t][i] * A[i][j] *
B[j][observations[t + 1]] * beta[t + 1][j] / denominator;
}
}
}
for (int i = 0; i < numStates; i++) {
newPi[i] += gamma[0][i];
}
for (int i = 0; i < numStates; i++) {
for (int j = 0; j < numStates; j++) {
double sum = 0;
for (int t = 0; t < T - 1; t++) {
sum += xi[t][i][j];
}
newA[i][j] += sum;
}
}
for (int j = 0; j < numStates; j++) {
for (int k = 0; k < numObservations; k++) {
double sum = 0;
double denom = 0;
for (int t = 0; t < T; t++) {
denom += gamma[t][j];
if (observations[t] == k) {
sum += gamma[t][j];
}
}
newB[j][k] += sum;
}
}
}
double piSum = 0;
for (double v : newPi) piSum += v;
for (int i = 0; i < numStates; i++) pi[i] = newPi[i] / piSum;
for (int i = 0; i < numStates; i++) {
double rowSum = 0;
for (int j = 0; j < numStates; j++) rowSum += newA[i][j];
for (int j = 0; j < numStates; j++) A[i][j] = newA[i][j] / rowSum;
}
for (int i = 0; i < numStates; i++) {
double rowSum = 0;
for (int j = 0; j < numObservations; j++) rowSum += newB[i][j];
for (int j = 0; j < numObservations; j++) B[i][j] = newB[i][j] / rowSum;
}
}
}
public int[] generateSequence(int length, long seed) {
Random rand = new Random(seed);
int[] states = new int[length];
int[] observations = new int[length];
double r = rand.nextDouble();
double cumSum = 0;
for (int i = 0; i < numStates; i++) {
cumSum += pi[i];
if (r <= cumSum) {
states[0] = i;
break;
}
}
r = rand.nextDouble();
cumSum = 0;
for (int j = 0; j < numObservations; j++) {
cumSum += B[states[0]][j];
if (r <= cumSum) {
observations[0] = j;
break;
}
}
for (int t = 1; t < length; t++) {
r = rand.nextDouble();
cumSum = 0;
for (int j = 0; j < numStates; j++) {
cumSum += A[states[t - 1]][j];
if (r <= cumSum) {
states[t] = j;
break;
}
}
r = rand.nextDouble();
cumSum = 0;
for (int j = 0; j < numObservations; j++) {
cumSum += B[states[t]][j];
if (r <= cumSum) {
observations[t] = j;
break;
}
}
}
return observations;
}
public static void main(String[] args) {
HiddenMarkovModel hmm = new HiddenMarkovModel(3, 4);
hmm.pi = new double[]{0.6, 0.3, 0.1};
hmm.A = new double[][]{
{0.7, 0.2, 0.1},
{0.3, 0.5, 0.2},
{0.2, 0.3, 0.5}
};
hmm.B = new double[][]{
{0.5, 0.3, 0.1, 0.1},
{0.2, 0.4, 0.2, 0.2},
{0.1, 0.2, 0.3, 0.4}
};
int[] testObs = new int[]{0, 1, 2};
System.out.println("=== HMM演示 ===");
System.out.println("观测序列: 散步 -> 购物 -> 打扫");
double prob = hmm.forwardLog(testObs);
System.out.println("观测序列概率(log): " + prob);
int[] bestPath = hmm.viterbi(testObs);
String[] stateNames = {"晴天", "多云", "雨天"};
System.out.print("最可能状态序列: ");
for (int s : bestPath) System.out.print(stateNames[s] + " ");
System.out.println();
System.out.println("\n=== Baum-Welch学习 ===");
int[][] trainData = new int[50][];
for (int i = 0; i < 50; i++) {
trainData[i] = hmm.generateSequence(20, i);
}
HiddenMarkovModel learned = new HiddenMarkovModel(3, 4);
learned.baumWelch(trainData, 100);
System.out.println("学习后的转移矩阵A:");
for (int i = 0; i < 3; i++) {
System.out.printf("%.3f %.3f %.3f%n",
learned.A[i][0], learned.A[i][1], learned.A[i][2]);
}
}
}
算法复杂度分析
- 前向算法:时间复杂度O(T * N^2),空间复杂度O(T * N),其中T为序列长度,N为状态数。
- 维特比算法:时间复杂度O(T * N^2),空间复杂度O(T * N)。
- Baum-Welch算法:每轮迭代时间复杂度O(K * T * N^2),其中K为训练序列数量。
应用场景
HMM是语音识别(如HTK工具包)、自然语言处理(词性标注、命名实体识别)、生物信息学(DNA序列分析、蛋白质结构预测)和手势识别等领域的核心技术。理解HMM是掌握更高级概率图模型(如CRF、MEMM)的基础。