向前-向后算法（forward-backward algorithm）

本文承接上篇博客《隐马尔可夫模型及的评估和解码问题》，用到的概念和例子都是那里面的。

学习问题

在HMM模型中，已知隐藏状态的集合S，观察值的集合O，以及一个观察序列（o₁,o₂,...,o_n），求使得该观察序列出现的可能性最大的模型参数（包括初始状态概率矩阵π，状态转移矩阵A，发射矩阵B）。这正好就是EM算法要求解的问题：已知一系列的观察值X，在隐含变量Y未知的情况下求最佳参数θ*，使得：

在中文词性标注里，根据为训练语料，我们观察到了一系列的词（对应EM中的X），如果每个词的词性（即隐藏状态）也是知道的，那它就不需要用EM来求模型参数θ了，因为Y是已知的，不存在隐含变量了。当没有隐含变量时，直接用maximum likelihood就可以把模型参数求出来。

预备知识

首先你得对下面的公式表示认同。

以下都是针对相互独立的事件，

P(A,B)=P(B|A)*P(A)

P(A,B,C)=P(C)*P(A,B|C)=P(A,C|B)*P(B)=P(B,C|A)*P(A)

P(A,B,C,D)=P(D)*P(A,B|D)*P(C|A)=P(D)*P(A,B|D)*P(C|B)

P(A,B|C)=P(D1,A,B|C)+P(D2,A,B|C)     D1,D2是事件D的一个全划分

理解了上面几个式子，你也就能理解本文中出现的公式是怎么推导出来的了。

EM算法求解

我们已经知道如果隐含变量Y是已知的，那么求解模型参数直接利用Maximum Likelihood就可以了。EM算法的基本思路是：随机初始化一组参数θ⁽⁰⁾，根据后验概率Pr(Y|X;θ)来更新Y的期望E(Y)，然后用E(Y)代替Y求出新的模型参数θ⁽¹⁾。如此迭代直到θ趋于稳定。

在HMM问题中，隐含变量自然就是状态变量，要求状态变量的期望值，其实就是求时刻t_i观察到x_i时处于状态s_i的概率，为了求此概率，需要用到向前变量和向后变量。

向前变量

向前变量 是假定的参数

它表示t时刻满足状态$S_i$，且t时刻之前（包括t时刻）满足给定的观测序列$(O_{1}O_{2}{\cdots}O_t)$的概率。

1. 令初始值$$\alpha_1(i)=\pi_ib_i(O_1)$$
2. 归纳法计算$$\alpha_{t+1}(j)=\left[\sum_{i=1}^N{\alpha_t(i)a_{ij}}\right]b_j(O_{t+1})$$
3. 最后计算$$P(O|\lambda)=\sum_{i=1}^N{\alpha_T(i)}$$

复杂度$O(N^2T)$

向后变量

向后变量

$$\beta_t(i)=p(O_{t+1}O_{t+2}{\ldots}O_T|q_t=S_i\;\lambda),1{\le}t{\le}T\;\;1{\le}i{\le}N$$

它表示在时刻t出现状态$S_i$，且t时刻以后的观察序列满足$(O_{t+1}O_{t+2}{\cdots}O_T)$的概率。

1. 初始值$\beta_T(i)=1$
2. 归纳计算$$\beta_t(i)=\sum_{j=1}^N{a_{ij}b_j(O_{t+1})\beta_{t+1}(j)}$$

E-Step

定义变量$\xi_t(i,j)$为t时刻处于状态i，t+1时刻处于状态j的概率。

$$\xi_t(i,j)=p(q_t=S_i,q_{t+1}=S_j|O;\lambda),1{\le}t{\le}T-1$$

定义变量$\gamma_t(i)$表示t时刻呈现状态i的概率。

$\gamma_t(i)=p(q_t=S_i|O;\lambda)$

实际上$$\gamma_t(i)=\sum_{j=1}^N\xi_t(i,j),1{\le}t{\le}T-1$$       

$$\gamma_t(i)=\frac{\alpha_t(i)\beta_t(i)}{\sum_{i=1}^N{\alpha_t(i)\beta_t(i)}},1{\le}t{\le}T$$ 

$\sum_{t=1}^{T}\gamma_t(i)$是从其他所有状态转移到状态i的次数的期望值。

$\sum_{t=1}^{T-1}\gamma_t(i)$是从状态i转移出去的次数的期望值。

$\sum_{t=1}^{T-1}\xi_t(i,j)$是从状态i转移到状态j的次数的期望值。

M-Step

$\overline{\pi_i}$是在初始时刻出现状态i的频率的期望值，$\overline{\pi_i}=\gamma_1(i)$

$\overline{\alpha}_{ij}$是从状态i转移到状态j的次数的期望值  除以  从状态i转移出去的次数的期望值，$$\overline{a}_{ij}=\frac{\sum_{t=1}^{T-1}\xi_t(i,j)}{\sum_{t=1}^{T-1}\gamma_t(i)}$$

$\overline{b}_j(k)$是在状态j下观察到活动为k的次数的期望值  除以  从其他所有状态转移到状态j的次数的期望值,$$\overline{b}_j(k)=\frac{\sum_{t=1,s.t.O_t=k}^T\gamma_t(j)}{\sum_{t=1}^T\gamma_t(j)}$$

然后用新的参数$(\Pi,A,B)$再来计算向前变量、向后变量、$\xi_t(i,j)$和$\gamma_t(i)$。如此循环迭代，直到前后两次参数的变化量小于某个值为止。

下面给出我的java代码：

import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.io.IOException;
import java.util.Arrays;
import java.util.HashMap;
import java.util.LinkedList;
import java.util.List;
import java.util.Map;
import java.util.Map.Entry;

/**
 * 隐马尔可夫模型参数学习。
 *
 * @Author:zhangchaoyang
 * @Since:2015年4月4日
 * @Version:1.0
 */
public class HmmLearn {

    private int stateCount;// 状态的个数
    private Map<String, Integer> observeIndexMap = new HashMap<String, Integer>();// 观察值及其索引编号
    /**
     * 通过学习得到的模型参数
     */
    private double[] stateProb;// 初始状态概率矩阵
    private double[][] stateTrans;// 状态转移矩阵
    private double[][] emission;// 混淆矩阵

    private List<String> observeSeqs = new LinkedList<String>();// 训练集中所有的观察序列

    /**
     * 迭代终止条件
     */
    private final int ITERATION_MAX = 100;
    private final double DELTA_PI = 1E-3;
    private final double DELTA_A = 1E-2;
    private final double DELTA_B = 1E-2;

    /**
     * 
     * @param stateCount
     *            指定状态取值有多少种
     * @param observeFile
     *            存储观察序列的文件，各个观察序列用空白符或换行符隔开即可
     * @throws IOException
     */
    public void initParam(int stateCount, String observeFile)
            throws IOException {
        this.stateCount = stateCount;
        int observeCount = 0;
        BufferedReader br = new BufferedReader(new FileReader(new File(
                observeFile)));
        String line = null;
        while ((line = br.readLine()) != null) {
            String[] arr = line.split("\\s+");
            for (String seq : arr) {
                if (seq.length() > 1) {// 长度为1的观察序列必须过滤掉，不然在更新stateTrans时会出现NaN的情况
                    observeSeqs.add(seq);
                    for (int i = 0; i < seq.length(); i++) {
                        String observe = seq.substring(i, i + 1);
                        if (!observeIndexMap.containsKey(observe)) {
                            observeIndexMap.put(observe, observeCount++);
                        }
                    }
                }
            }
        }
        br.close();

        stateProb = new double[stateCount];
        initWeightRandomly(stateProb, 1E5);
        // initWeightEqually(stateProb);
        stateTrans = new double[stateCount][];
        for (int i = 0; i < stateCount; i++) {
            stateTrans[i] = new double[stateCount];
            initWeightRandomly(stateTrans[i], 1E5);
            // initWeightEqually(stateTrans[i]);
        }
        emission = new double[stateCount][];
        for (int i = 0; i < stateCount; i++) {
            emission[i] = new double[observeCount];
            initWeightRandomly(emission[i], 1E9);
            // initWeightEqually(emission[i]);
        }
    }

    /**
     * 随机地初始化权重，使得各权重非负，且和为1.
     * 
     * @param arr
     * @param precision
     */
    public void initWeightRandomly(double[] arr, double precision) {
        int len = arr.length - 1;
        int[] position = new int[len];
        for (int i = 0; i < len; i++) {
            position[i] = (int) (Math.random() * precision);
        }
        Arrays.sort(position);
        int pre = 0;
        for (int i = 0; i < len; i++) {
            arr[i] = 1.0 * (position[i] - pre) / precision;
            pre = position[i];
        }
        arr[len] = 1.0 * (precision - pre) / precision;
    }

    /**
     * 均等地初始化权重，使得各权重非负，且和为1.
     * 
     * @param arr
     */
    public void initWeightEqually(double[] arr) {
        int len = arr.length;
        for (int i = 0; i < len; i++) {
            arr[i] = 1.0 / len;
        }
    }

    /**
     * BaumWelch算法学习HMM的模型参数
     */
    public void baumWelch() {
        long begin = System.currentTimeMillis();
        int iter = 0;
        while (iter++ < ITERATION_MAX) {
            double[] stateProb_new = new double[stateCount];
            double[][] stateTrans_new = new double[stateCount][];
            double[][] emission_new = new double[stateCount][];
            for (int i = 0; i < stateCount; i++) {
                stateTrans_new[i] = new double[stateCount];
            }
            for (int i = 0; i < stateCount; i++) {
                emission_new[i] = new double[observeIndexMap.size()];
            }
            for (String seq : observeSeqs) {
                int T = seq.length();
                double[][] alpha = new double[T][];
                double[][] beta = new double[T][];
                double[][] gamma = new double[T][];
                for (int i = 0; i < T; i++) {
                    alpha[i] = new double[stateCount];
                    beta[i] = new double[stateCount];
                    gamma[i] = new double[stateCount];
                }
                double[][][] xi = new double[T - 1][][];
                for (int i = 0; i < T - 1; i++) {
                    xi[i] = new double[stateCount][];
                    for (int j = 0; j < stateCount; j++) {
                        xi[i][j] = new double[stateCount];
                    }
                }
                String observe = seq.substring(0, 1);
                int observeIndex = observeIndexMap.get(observe);
                // 计算alpha
                for (int i = 0; i < stateCount; i++) {
                    alpha[0][i] = stateProb[i] * emission[i][observeIndex];
                }
                for (int t = 1; t < T; t++) {
                    observe = seq.substring(t, t + 1);// 当前时刻的观察值
                    observeIndex = observeIndexMap.get(observe);
                    for (int j = 0; j < stateCount; j++) {
                        double sum = 0;
                        for (int i = 0; i < stateCount; i++) {
                            sum += alpha[t - 1][i] * stateTrans[i][j];
                        }
                        alpha[t][j] = sum * emission[j][observeIndex];
                    }
                }
                // 计算beta
                for (int i = 0; i < stateCount; i++) {
                    beta[T - 1][i] = 1;
                }
                for (int t = T - 2; t >= 0; t--) {
                    observe = seq.substring(t + 1, t + 2);// 当前时刻的后一个观察值
                    observeIndex = observeIndexMap.get(observe);
                    for (int i = 0; i < stateCount; i++) {
                        double sum = 0;
                        for (int j = 0; j < stateCount; j++) {
                            sum += beta[t + 1][j] * stateTrans[i][j]
                                    * emission[j][observeIndex];
                        }
                        beta[t][i] = sum;
                    }
                }
                double[] denominator = new double[T];// 在计算gamma和xi时都要计算的一个变量
                // 计算gamma
                for (int t = 0; t < T; t++) {
                    double sum = 0;
                    for (int j = 0; j < stateCount; j++) {
                        sum += alpha[t][j] * beta[t][j];
                    }
                    denominator[t] = sum;
                    for (int i = 0; i < stateCount; i++) {
                        gamma[t][i] = alpha[t][i] * beta[t][i] / sum;
                    }
                }
                // 计算xi
                for (int t = 0; t < T - 1; t++) {
                    observe = seq.substring(t + 1, t + 2);// 当前时刻的后一个观察值
                    observeIndex = observeIndexMap.get(observe);
                    for (int i = 0; i < stateCount; i++) {
                        for (int j = 0; j < stateCount; j++) {
                            xi[t][i][j] = alpha[t][i] * stateTrans[i][j]
                                    * beta[t + 1][j]
                                    * emission[j][observeIndex]
                                    / denominator[t];
                        }
                    }
                }
                // 计算stateProb
                double[] curr_stateProb = new double[stateCount];
                for (int i = 0; i < stateCount; i++) {
                    curr_stateProb[i] = gamma[0][i];
                    stateProb_new[i] += curr_stateProb[i];
                }
                // 计算stateTrans
                double[][] curr_stateTrans = new double[stateCount][];
                for (int i = 0; i < stateCount; i++) {
                    curr_stateTrans[i] = new double[stateCount];
                    for (int j = 0; j < stateCount; j++) {
                        double up = 0;
                        double down = 0;
                        for (int t = 0; t < T - 1; t++) {
                            up += xi[t][i][j];
                            down += gamma[t][i];
                        }
                        if (down > 0) {
                            curr_stateTrans[i][j] = up / down;// 任何浮点操作，只要它的一个或多个操作数为NaN，其结果就是NaN。NaN不等于任何浮点数，包括它自身在内。
                            stateTrans_new[i][j] += curr_stateTrans[i][j];
                        } else {
                            stateTrans_new[i][j] += stateTrans[i][j];// 如果分母为0无法相除，就找个值来替代
                            System.err.println("up=" + up + ",down=" + down);
                        }
                    }
                }
                // 计算emission
                double[][] curr_emission = new double[stateCount][];
                for (int i = 0; i < stateCount; i++) {
                    curr_emission[i] = new double[observeIndexMap.size()];
                    for (Entry<String, Integer> entry : observeIndexMap
                            .entrySet()) {
                        String obs = entry.getKey();
                        int index = entry.getValue();
                        double up = 0;
                        double down = 0;
                        for (int t = 0; t < T; t++) {
                            observe = seq.substring(t, t + 1);
                            if (obs.equals(observe)) {
                                up += gamma[t][i];
                            }
                            down += gamma[t][i];
                        }
                        curr_emission[i][index] = up / down;
                        emission_new[i][index] += curr_emission[i][index];
                    }
                }
            }
            // 批量更新模型参数
            double delta_pi = 0;
            double delta_a = 0;
            double delta_b = 0;
            int seqCount = observeSeqs.size();
            for (int i = 0; i < stateCount; i++) {
                stateProb_new[i] /= seqCount;
                delta_pi += Math.abs(stateProb_new[i] - stateProb[i]);
            }
            for (int i = 0; i < stateCount; i++) {
                for (int j = 0; j < stateCount; j++) {
                    stateTrans_new[i][j] /= seqCount;
                    delta_a += Math
                            .abs(stateTrans_new[i][j] - stateTrans[i][j]);
                }
            }
            for (int i = 0; i < stateCount; i++) {
                for (int j = 0; j < observeIndexMap.size(); j++) {
                    emission_new[i][j] /= seqCount;
                    delta_b += Math.abs(emission_new[i][j] - emission[i][j]);
                }
            }
            System.out.println("iteration " + iter + ", delta_pi=" + delta_pi
                    + ", delta_a=" + delta_a + ", delta_b=" + delta_b);
            if (delta_pi <= DELTA_PI && delta_a <= DELTA_A
                    && delta_b <= DELTA_B) {
                break;
            } else {
                stateProb = stateProb_new;
                stateTrans = stateTrans_new;
                emission = emission_new;
            }
        }
        long end = System.currentTimeMillis();
        System.out.println("time elapse " + (end - begin) / 1000 + " seconds.");
    }

    public int getStateCount() {
        return stateCount;
    }

    public double[] getStateProb() {
        return stateProb;
    }

    public double[][] getStateTrans() {
        return stateTrans;
    }

    public double[][] getEmission() {
        return emission;
    }

    public Map<String, Integer> getObserveIndexMap() {
        return observeIndexMap;
    }

}

测试代码：

import java.io.BufferedWriter;
import java.io.File;
import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayDeque;
import java.util.ArrayList;
import java.util.Deque;
import java.util.LinkedList;
import java.util.List;
import java.util.Map;
import java.util.Map.Entry;

import org.junit.Assert;
import org.junit.BeforeClass;
import org.junit.Test;

import xxzl.dm.core.sequence.HmmInference;
import xxzl.dm.core.sequence.HmmLearn;
import xxzl.dm.utility.FileUtil;
import xxzl.dm.utility.Pair;

public class TestHmmLearn {

    private static String corpusPath;
    private static String modelPath;

    @BeforeClass
    public static void setup() throws IOException {
        corpusPath = TestHmmLearn.class.getResource("/").getPath()
                + "/resources/corpus/";
        modelPath = TestHmmLearn.class.getResource("/").getPath()
                + "/resources/model/hmm/";
    }

    /**
     * 先通过BaumWelch算法去学习模型参数,并把学习到的参数写入文件
     */
    @Test
    public void printParam() throws IOException {
        HmmLearn hmmLearn = new HmmLearn();
        hmmLearn.initParam(4, corpusPath + "msr_test");
        hmmLearn.baumWelch();
        double[] PI = hmmLearn.getStateProb();
        double[][] A = hmmLearn.getStateTrans();
        double[][] B = hmmLearn.getEmission();
        int stateCount = PI.length;
        List<String> stateSet = new ArrayList<String>();
        for (int i = 0; i < stateCount; i++) {
            stateSet.add(getLabelName(i + 1));
        }
        FileUtil.writeLines(modelPath + "state", stateSet);
        List<String> observeSet = new LinkedList<String>();
        Map<String, Integer> observeIndex = hmmLearn.getObserveIndexMap();
        for (Entry<String, Integer> entry : observeIndex.entrySet()) {
            observeSet.add(entry.getKey());
        }
        FileUtil.writeLines(modelPath + "observe", observeSet);
        List<String> stateProb = new ArrayList<String>();
        for (double pi : PI) {
            stateProb.add(String.valueOf(pi));
        }
        FileUtil.writeLines(modelPath + "PI", stateProb);
        BufferedWriter bw = new BufferedWriter(new FileWriter(new File(
                modelPath + "A")));
        for (int i = 0; i < A.length; i++) {
            for (int j = 0; j < A[i].length; j++) {
                bw.write(A[i][j] + "\t");
            }
            bw.newLine();
        }
        bw.close();
        bw = new BufferedWriter(new FileWriter(new File(modelPath + "B")));
        for (int i = 0; i < B.length; i++) {
            for (int j = 0; j < B[i].length; j++) {
                bw.write(B[i][j] + "\t");
            }
            bw.newLine();
        }
        bw.close();
    }

    /**
     * 再用学习好的参数去做inference
     */
    @Test
    public void testInference() {
        HmmInference hmmInference = new HmmInference();
        hmmInference.initStateSet(modelPath + "state");
        hmmInference.initObserveSet(modelPath + "observe");
        hmmInference.initStateProb(modelPath + "PI");
        hmmInference.initStateTrans(modelPath + "A");
        hmmInference.initConfusion(modelPath + "B");
        String sentence = "他们的成长与优越的家庭背景没有任何正相关的关系";
        List<String> obs_seq = str2List(sentence);
        Pair<Double, LinkedList<String>> states = hmmInference.viterbi(obs_seq);
        System.out.println("最可能的状态序列：" + states.second + "，其概率为："
                + states.first);
        double ratio = hmmInference.estimate(obs_seq);
        System.out.println("句子1出现的概率：" + ratio);
        sentence = "为报倾城随太守亲射虎看孙郎会挽雕弓如满月西北望";// 要保证两个句子长度相等，比较其出现概率才有意义
        obs_seq = str2List(sentence);
        double ratio2 = hmmInference.estimate(obs_seq);
        System.out.println("句子2出现的概率：" + ratio2);
        Assert.assertTrue(ratio > ratio2);
    }

    /**
     * 按照Excel的规则，给一个列号，给出列名称。比如1对应A，26对应Z，27对应AA，53对应BA
     * 
     * @param a
     * @return
     */
    private static String getLabelName(int a) {
        StringBuilder sb = new StringBuilder();
        char[] arr = new char[26];
        arr[0] = 'A';
        for (int i = 1; i < 26; i++) {
            arr[i] = (char) (arr[i - 1] + 1);
        }
        Deque<String> deque = new ArrayDeque<String>();
        while (a > 0) {
            char c = arr[(a - 1) % 26];
            deque.push(String.valueOf(c));
            a = (a - 1) / 26;
        }
        while (deque.size() > 0) {
            sb.append(deque.pop());
        }
        return sb.toString();
    }

    private List<String> str2List(String sentence) {
        List<String> obs_seq = new ArrayList<String>();
        int i = 0;
        while (i < sentence.length()) {
            obs_seq.add(sentence.substring(i, i + 1));
            i++;
        }
        return obs_seq;
    }
}

 其中用到的类HmmInference原代码在上篇博客中。