每日算法 — 使用java实现HMM:隐马尔可夫模型与序列状态推断

算法概述

隐马尔可夫模型(Hidden Markov Model, HMM)是一种统计模型,用于描述由隐藏的马尔可夫链随机生成观测序列的过程。HMM假设系统在每个时刻处于某个隐藏状态,隐藏状态之间的转移服从马尔可夫性质,而每个隐藏状态会以一定概率产生观测值。

HMM由三个核心要素定义:
– 初始概率分布pi:系统在t=1时处于各状态的概率
– 状态转移矩阵A:隐藏状态之间的转移概率
– 观测概率矩阵B:各状态下生成不同观测值的概率

HMM三大经典问题

  1. 评估问题:给定模型和观测序列O,计算P(O|lambda) —— 前向算法解决
  2. 解码问题:给定模型和观测序列O,找出最可能的隐藏状态序列 —— 维特比算法解决
  3. 学习问题:给定观测序列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)的基础。