隐马尔可夫模型及的评估和解码问题

HMM介绍

Hidden Markov Models是一种统计信号处理方法,模型中包含2个序列和3个矩阵:状态序列S、观察序列O、初始状态矩阵P、状态转移矩阵A、混淆矩阵B。举个例子来说明。

你一个异地的朋友只做三种活动：散步、看书、做清洁。每天只做一种活动。假设天气只有两种状态：晴和兩。每天只有一种天气。你的朋友每天告诉你他做了什么，但是不告诉你他那里的天气。

某一周从周一到周五每天的活动分别是{读书，做清洁，散步，做清洁，散步}----这就是观察序列O，因为你可以观察得到。

从周一到周五的天气依次是{晴，兩，晴，晴，晴}----这就是状态序列S，状态序列是隐藏的，你不知道。

根据长期统计，某天晴的概率是0.6，兩的概率是0.4。则$\pi=[0.6,0.4]$。

从晴转晴的概率是0.7，从晴转兩的概率是0.3。从兩转晴的概率是0.4，从兩转兩的概率是0.6。则$A=\left[\begin{array}{cc}0.7&0.3\\0.4&0.6\end{array}\right]$。

天气晴时，散步的概率是0.4，看书的概率是0.3，做清洁的概率是0.3。天气兩时，散步的概率是0.1，看书的概率是0.4，做清洁的概率是0.5。则$B=\left[\begin{array}{cc}0.4&0.3&0.3\\0.4&0.1&0.5\end{array}\right]$。

该模型和实际情况有明显不符的地方：用一简单的状态转移矩阵A来表示状态的转移概率的前提是t时刻的状态只跟t-1时刻的状态有关，而实际上今天的天气跟过去几天的天气都有关系，而且跟过去几天的晴朗程度、兩量大小都有关系；混淆矩阵B认为今天的活动只跟今天的天气有关系，实际上今天的活动跟过去几天的活动也有关系，比如过去一周都没有打扫房间，那今天做清洁的概率就大大增加。

模型介绍完了。

评估问题

隐马尔可夫模型中包含一个评估问题：已知模型参数，计算某一特定输出序列的概率。通常使用forward算法解决。 

比如计算活动序列{读书，做清洁，散步，做清洁，散步}出现的概率，就属于评估问题。

如果穷举的话，观察序列会有2^5种，需要分别计算它们出现的概率，然后找出概率最大的。

穷举法中有很多重复计算，向前算法就是利用已有的结果，减少重复计算。

算法借助于一个矩阵Q[LEN][M]，其中M是所有状态的种数，Q[i][j]表示从第0天到第i天，满足观察序列，且第i天隐藏状态为Sj的所有可能的隐藏序列的概率之和。最终所求结果为Q[LEN-1][0]+...+Q[LEN-1][M-1]，即最后一天，所有隐藏状态下，分别满足观察序列的概率值之和。

比如Q[0][0]=p(第一天做卫生 且 第一天晴)=p(天晴)*p(做卫生|天晴)=P[0]*B[0][2]=0.6*0.3=0.18

Q[0][1]=p(第一天做卫生 且 第一天下雨)=p(下雨)*p(做卫生|下雨)=P[1]*B[1][2]=0.4*0.5=0.2

Q[1][0]=p(第一天做卫生 且 第二天晴 且 第二天做卫生)

=p(第一天做卫生 且 第二天晴)*p(天晴的情况下做卫生)

=p{ p(第一天做卫生 且 第一天晴)*p(从天晴转天晴)+p(第一天做卫生 且 第一天下雨)*p(从下雨转天晴) }*p(天晴的情况下做卫生)

={ Q[0][0]*A[0][0] + Q[0][1]*A[1][0] } * B[0][2]

Q[1][1]= ……

…… ……

可以看到计算Q矩阵的每i行时都用到了第i-1行的结果。

解码问题

解码问题是：已知模型参数，寻找最可能的能产生某一特定输出序列O(LEN)的隐含状态的序列。通常使用Viterbi算法解决。 

观察序列长度为LEN，则隐藏状态序列长度也是LEN，如果采用穷举法，就有M^LEN种可能的隐藏状态序列，我们要计算每一种隐藏状态到指定观察序列的概率，最终选择概率最大的。

穷举法中有很多重复计算，Viterbi算法就是利用已有的结果，减少重复计算。

跟评估问题非常相似，不同点在于评估算的是和，解码算的是最大值。

Viterbi算法主要就是在计算一个矩阵Q[LEN][M]，其中Q[i][j]表示从第0天到第i天，满足观察序列，且第i天隐藏状态为Sj的所有可能的隐藏序列的概率的最大值。另外还要建立一个矩阵Path[LEN][M]，用来记录状态序列中某一状态之前最可能的状态。

举个例子，假如指定观察序列是{读书，做卫生，散步，做卫生，散步}，求出现此观察序列最可能的状态序列是什么。

Q[0][0]=p(第一天读书 且 第一天晴)=p(天晴)*p(读书|天晴)

Path[0][0]=-1;

Q[0][1]=p(第一天读书 且 第一天下雨)=p(下雨)*p(读书|下雨)

Path[0][1]=-1;

关键是从第二天开始，Q[1][0]表示：满足“第一天读书 且 第二做卫生 且 第二天晴”的所有可能的隐藏序列的概率的最大值。那么满足“第一天读书 且 第二做清洁 且 第二天晴”的所有可能的隐藏序列有哪些呢？第二天是必须满足晴天的，第二天之前的状态可以任意变。则所有可能的隐藏序列就是“晴  晴”和“雨  晴”。实际上考虑第三天（及第三天以后）时，并不需要考虑“所有”可能的隐藏序列，而只需要考虑第二天的不同状态取值，这是因为马氏过程有无后效性--tm时刻所处状态的概率只和tm-1时刻的状态有关，而与tm-1时刻之前的状态无关。

Q[1][0]=

max{ p(第一天晴 且 第一天读书 且 第二天晴 且第二天做卫生) ，p(第一天下雨 且 第一天读书 且 第二天晴 且 第二天做卫生) }

=max{ p(第一天读书 且 第一天晴)*p(天晴转天晴)，p(第一天读书 且 第一天下雨)*p(下雨转天晴) } * p(做卫生|天晴)

=max{ Q[0][0]*A[0][0]，Q[0][1]*A[1][0] } * B[0][2]

假如Q[0][0]*A[0][0] < Q[0][1]*A[1][0]，则Path[1][0]=1；假如Q[0][0]*A[0][0] > Q[0][1]*A[1][0]，则Path[1][0]=0。

Q[1][1]= ……

…… ……

可以看到计算Q矩阵的每i行时都用到了第i-1行的结果。

下面给出两种算法的Java代码：

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

import xxzl.dm.utility.FileUtil;
import xxzl.dm.utility.Pair;
import xxzl.dm.utility.math.Smooth;

/**
 * 隐马尔可夫推断，包括评估问题和解码问题
 * 
 * @Author:zhangchaoyang
 * @Since:2015年3月29日
 * @Version:1.0
 */
public class HmmInference {

    /**
     * HMM的模型参数
     */
    private List<String> stateSet = new LinkedList<String>();// 状态值集合
    private List<String> observeSet = new LinkedList<String>();// 观察值集合
    private double[] stateProb;// 初始状态概率矩阵
    private double[][] stateTrans;// 状态转移矩阵
    private double[][] emission;// 发射矩阵
    private double[] minEmission;// 发射矩阵每一行的极小值（用于零值平滑）

    /**
     * 使用已标记好的训练样本，经过简单统计使用极大似然确定HMM的参数
     * 
     * @param tagFile
     *            文件格式：每行2列，第1列是观察值，第2列是状态值，不用序列之间用空行隔开。<br>
     *            /test/resources/corpus/wordcut.train是一个示例文件。
     */
    public void initParam(String tagFile) {
        Map<String, Integer> stateIndexMap = new HashMap<String, Integer>();// 状态值及其编号
        Map<String, Integer> observeIndexMap = new HashMap<String, Integer>();// 观察值及其编号
        int[] stateCount;// 状态值及其计数
        int[][] stateTransCount;// 状态转移计数矩阵
        int[][] confusionCount;// 混淆计数矩阵

        try {
            BufferedReader br = new BufferedReader(new FileReader(new File(
                    tagFile)));
            if (br.markSupported()) {
                br.mark(1024 * 1024 * 100);
            }
            String line = null;
            int stateTotal = 0;
            int observeTotal = 0;
            // 第一趟扫描文件，给stateIndexMap、observeIndexMap赋值
            while ((line = br.readLine()) != null) {
                String[] arr = line.split("\\s+");// 每行存储一个观察值、一个状态值，用空白符隔开。用空行隔开不同的序列。
                if (arr.length >= 2) {
                    String observe = arr[0];
                    String state = arr[1];
                    if (!observeIndexMap.containsKey(observe)) {
                        observeIndexMap.put(observe, observeTotal++);
                    }
                    if (!stateIndexMap.containsKey(state)) {
                        stateIndexMap.put(state, stateTotal++);
                    }
                }
            }
            if (br.markSupported()) {
                br.reset();
            } else {
                br.close();
                br = new BufferedReader(new FileReader(new File(tagFile)));
            }
            // System.out.println("state set:");
            // for (Entry<String, Integer> entry : stateIndexMap.entrySet()) {
            // System.out.println(entry.getValue() + ":" + entry.getKey());
            // }
            // 第二趟扫描文件，给stateTransCount、confusionCount、stateCount赋值
            for (int i = 0; i < stateTotal; i++) {
                stateSet.add("");
            }
            for (int i = 0; i < observeTotal; i++) {
                observeSet.add("");
            }
            stateTransCount = new int[stateTotal][];
            for (int i = 0; i < stateIndexMap.size(); i++) {
                stateTransCount[i] = new int[stateTotal];
            }
            confusionCount = new int[stateTotal][];
            for (int i = 0; i < stateTotal; i++) {
                confusionCount[i] = new int[observeTotal];
            }
            stateCount = new int[stateTotal];
            String preState = null;
            while ((line = br.readLine()) != null) {
                String[] arr = line.split("\\s+");
                if (arr.length >= 2) {
                    String observe = arr[0];
                    String state = arr[1];
                    int row = stateIndexMap.get(state);
                    int col = 0;
                    int oldCount = 0;
                    // if (observeIndexMap.containsKey(observe)) {
                    col = observeIndexMap.get(observe);
                    oldCount = confusionCount[row][col];
                    confusionCount[row][col] = oldCount + 1;
                    // }
                    stateCount[row] = stateCount[row] + 1;
                    if (preState == null) {
                        preState = state;
                    } else {
                        row = stateIndexMap.get(preState);
                        col = stateIndexMap.get(state);
                        oldCount = stateTransCount[row][col];
                        stateTransCount[row][col] = oldCount + 1;
                        preState = state;
                    }
                } else {
                    preState = null;
                }
            }
            br.close();
            // 给HMM基本参数赋值
            for (Entry<String, Integer> entry : stateIndexMap.entrySet()) {
                String state = entry.getKey();
                int index = entry.getValue();
                stateSet.set(index, state);
            }
            for (Entry<String, Integer> entry : observeIndexMap.entrySet()) {
                String observe = entry.getKey();
                int index = entry.getValue();
                observeSet.set(index, observe);
            }
            stateProb = calProbByCount(Smooth.GoodTuring(stateCount));
            // System.out.println("state initial prob:");
            // System.out.println(Arrays.toString(stateProb));
            stateTrans = new double[stateTransCount.length][];
            for (int i = 0; i < stateTransCount.length; i++) {
                stateTrans[i] = calProbByCount(stateTransCount[i]);// 计算状态转移概率时不作平滑，因为有些状态之间转移的概率就应该是0，如果平滑就会变成非0
            }
            // System.out.println("state transiation prob:");
            // for (int i = 0; i < stateTransProb.length; i++) {
            // System.out.println(Arrays.toString(stateTransProb[i]));
            // }
            emission = new double[confusionCount.length][];
            minEmission = new double[confusionCount.length];
            for (int i = 0; i < confusionCount.length; i++) {
                emission[i] = calProbByCount(Smooth
                        .GoodTuring(confusionCount[i]));
                double min = Double.MAX_VALUE;
                for (double ele : emission[i]) {
                    if (ele < min) {
                        min = ele;
                    }
                }
                minEmission[i] = min;
            }
        } catch (FileNotFoundException e) {
            e.printStackTrace();
        } catch (IOException e) {
            e.printStackTrace();
        }

    }

    /**
     * 采用向前算法解决评估问题：给定HMM的所有参数，评估一个观察序列出现的概率。
     * 
     * @param obs_seq
     * @return
     */
    public double estimate(List<String> obs_seq) {
        double rect = 0.0;
        int LEN = obs_seq.size();
        double[][] Q = new double[LEN][];
        // 状态的初始概率，乘上隐藏状态到观察状态的条件概率。
        Q[0] = new double[stateSet.size()];
        for (int j = 0; j < stateSet.size(); j++) {
            if (observeSet.contains(obs_seq.get(0))) {
                Q[0][j] = stateProb[j]
                        * emission[j][observeSet.indexOf(obs_seq.get(0))];
            } else {
                Q[0][j] = stateProb[j] * minEmission[j];
                System.err.println("观察值'" + obs_seq.get(0) + "'在已标记样本中未出现过");
            }
        }
        // 首先从前一时刻的每个状态，转移到当前状态的概率求和，然后乘上隐藏状态到观察状态的条件概率。
        for (int i = 1; i < LEN; i++) {
            Q[i] = new double[stateSet.size()];
            for (int j = 0; j < stateSet.size(); j++) {
                double sum = 0.0;
                for (int k = 0; k < stateSet.size(); k++) {
                    sum += Q[i - 1][k] * stateTrans[k][j];
                }
                if (observeSet.contains(obs_seq.get(i))) {
                    Q[i][j] = sum
                            * emission[j][observeSet.indexOf(obs_seq.get(i))];
                } else {
                    Q[i][j] = sum * minEmission[j];
                    System.err
                            .println("观察值'" + obs_seq.get(0) + "'在已标记样本中未出现过");
                }
            }
        }
        for (int i = 0; i < stateSet.size(); i++)
            rect += Q[LEN - 1][i];
        return rect;
    }

    /**
     * 采用viterbi进行解码：给定HMM的所有参数，给一个观察序列，评估最可能的状态序列是什么。
     * 
     * @param observe
     * @return
     */
    public Pair<Double, LinkedList<String>> viterbi(List<String> observe) {
        LinkedList<String> sta = new LinkedList<String>();
        int LEN = observe.size();
        int M = stateSet.size();
        double[][] Q = new double[LEN][];
        int[][] Path = new int[LEN][];
        Q[0] = new double[M];
        Path[0] = new int[M];
        for (int j = 0; j < M; j++) {
            if (observeSet.contains(observe.get(0))) {// 观察值在训练样本中未出现过，则概率设为0
                Q[0][j] = stateProb[j]
                        * emission[j][observeSet.indexOf(observe.get(0))];
            } else {
                Q[0][j] = stateProb[j] * minEmission[j] / 2;
                System.err.println("观察值'" + observe.get(0) + "'在已标记样本中未出现过");
            }
            Path[0][j] = -1;
        }
        for (int i = 1; i < LEN; i++) {
            Q[i] = new double[M];
            Path[i] = new int[M];
            for (int j = 0; j < M; j++) {
                double max = 0.0;
                int index = 0;
                for (int k = 0; k < M; k++) {
                    if (Q[i - 1][k] * stateTrans[k][j] > max) {
                        max = Q[i - 1][k] * stateTrans[k][j];
                        index = k;
                    }
                }
                if (observeSet.contains(observe.get(i))) {
                    Q[i][j] = max
                            * emission[j][observeSet.indexOf(observe.get(i))];
                } else {
                    Q[i][j] = max * minEmission[j] / 2;
                    System.err
                            .println("观察值'" + observe.get(0) + "'在已标记样本中未出现过");
                }
                Path[i][j] = index;
            }
        }
        // 找到最后一个时刻呈现哪种状态的概率最大
        double max = 0;
        int index = 0;
        for (int i = 0; i < M; i++) {
            if (Q[LEN - 1][i] > max) {
                max = Q[LEN - 1][i];
                index = i;
            }
        }
        sta.add(stateSet.get(index));
        // 动态规划，逆推回去各个时刻出现什么状态概率最大
        for (int i = LEN - 1; i > 0; i--) {
            index = Path[i][index];
            sta.add(stateSet.get(index));
        }
        // 把状态序列再顺过来
        Collections.reverse(sta);
        return Pair.of(max, sta);
    }

    public void baumWelch() {

    }

    public void initStateSet(String infile) {
        FileUtil.readLines(infile, stateSet);
    }

    public void initObserveSet(String infile) {
        FileUtil.readLines(infile, observeSet);
    }

    public void initStateProb(String infile) {
        assert stateSet != null;
        int stateCount = stateSet.size();
        assert stateCount > 0;
        stateProb = new double[stateCount];
        List<String> lines = new LinkedList<String>();
        FileUtil.readLines(infile, lines);
        assert lines.size() >= stateCount;
        for (int i = 0; i < stateCount; i++) {
            stateProb[i] = Double.parseDouble(lines.get(i));
        }
    }

    public void initStateTrans(String infile) {
        List<String> lines = new LinkedList<String>();
        FileUtil.readLines(infile, lines);
        stateTrans = new double[lines.size()][];
        for (int i = 0; i < lines.size(); i++) {
            String[] conts = lines.get(i).split("\\s+");
            stateTrans[i] = new double[conts.length];
            for (int j = 0; j < stateTrans[i].length; j++) {
                stateTrans[i][j] = Double.parseDouble(conts[j]);
            }
        }
    }

    public void initConfusion(String infile) {
        List<String> lines = new LinkedList<String>();
        FileUtil.readLines(infile, lines);
        emission = new double[lines.size()][];
        minEmission = new double[lines.size()];
        for (int i = 0; i < lines.size(); i++) {
            String[] conts = lines.get(i).split("\\s+");
            emission[i] = new double[conts.length];
            double min = Double.MAX_VALUE;
            for (int j = 0; j < emission[i].length; j++) {
                double ele = Double.parseDouble(conts[j]);
                emission[i][j] = ele;
                if (ele < min) {
                    min = ele;
                }
            }
            minEmission[i] = min;
        }
    }

    public List<String> getStateSet() {
        return stateSet;
    }

    public void setStateSet(List<String> stateSet) {
        this.stateSet = stateSet;
    }

    public List<String> getObserveSet() {
        return observeSet;
    }

    public void setObserveSet(List<String> observeSet) {
        this.observeSet = observeSet;
    }

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

    public void setStateProb(double[] stateProb) {
        this.stateProb = stateProb;
    }

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

    public void setStateTrans(double[][] stateTrans) {
        this.stateTrans = stateTrans;
    }

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

    public void setEmission(double[][] emission) {
        this.emission = emission;
        minEmission = new double[emission.length];
        for (int i = 0; i < emission.length; i++) {
            double min = Double.MAX_VALUE;
            for (double ele : emission[i]) {
                if (ele < min) {
                    min = ele;
                }
            }
            minEmission[i] = min;
        }
    }

    /**
     * 通过一组计数计算概率
     * 
     * @param countArr
     * @return
     */
    private double[] calProbByCount(double[] countArr) {
        double sum = 0.0;
        for (double count : countArr) {
            sum += count;
        }
        double[] prob = new double[countArr.length];
        for (int i = 0; i < countArr.length; i++) {
            prob[i] = countArr[i] / sum;
        }
        return prob;
    }

    /**
     * 通过一组计数计算概率
     * 
     * @param countArr
     * @return
     */
    private double[] calProbByCount(int[] countArr) {
        double sum = 0.0;
        for (double count : countArr) {
            sum += count;
        }
        double[] prob = new double[countArr.length];
        for (int i = 0; i < countArr.length; i++) {
            prob[i] = countArr[i] / sum;
        }
        return prob;
    }
}

 测试代码：

import java.util.ArrayList;
import java.util.LinkedList;
import java.util.List;

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

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

public class TestHmmInference {

    private static HmmInference hmm = new HmmInference();

    @BeforeClass
    public static void setup() {
        hmm.initParam("/Users/zhangchaoyang/OneDrive/msr_train.corp");
    }

    @Test
    public void testForward() {
        String sentence = "公布这些事实将对日记主人及其他被涉及者带来何种影响";
        List<String> obs_seq = str2List(sentence);
        double ratio = hmm.estimate(obs_seq);
        System.out.println("观察序列出现的概率：" + ratio);
    }

    @Test
    public void testViterbi() {
        String sentence = "公布这些事实将对日记主人及其他被涉及者带来何种影响";
        List<String> obs_seq = str2List(sentence);
        Pair<Double, LinkedList<String>> states = hmm.viterbi(obs_seq);
        System.out.println("最可能的状态序列：" + states.second + "，其概率为："
                + states.first);
        System.out.println("分词结果为：" + wordSeg(sentence, states.second));
    }

    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;
    }

    private String wordSeg(String sentence, List<String> tag) {
        StringBuilder sb = new StringBuilder();
        int len = sentence.length() <= tag.size() ? sentence.length() : tag
                .size();
        for (int i = 0; i < len; i++) {
            String word = sentence.substring(i, i + 1);
            sb.append(word);
            if (tag.get(i).equals("E") || tag.get(i).equals("S")) {
                sb.append("\t");
            }
        }
        return sb.toString();
    }
}