3.初始状态概率 πi <script type="math/tex" id="MathJax-Element-20">\pi_i</script>的估计 π^i <script type="math/tex" id="MathJax-Element-21">\hat \pi_{i}</script>为S个样本中初始状态为 qi <script type="math/tex" id="MathJax-Element-22">q_i</script>的频率
由于监督学习需要使用训练数据，而人工标注训练数据往往代价很高，有时就会利用非监督学习的方法。

非监督学习问题

上述监督学习问题给定了大量一一对应的观察序列和隐藏序列，用最简单的概率统计方法就能求得转移矩阵，观测概率矩阵的频数，注意这里是频数而非概率。这部分的内容相对简单，但针对非监督学习问题时，由于多了隐藏变量，而系统的各种参数均未知，因此求解非监督学习问题时，就存在一定难度，本文用到的知识有极大似然估计，EM算法，基础概率论，如果对这些知识还不够熟悉的话，建议回到前言提到的链接，看完链接内容后，对理解Baum-Welch算法将大有帮助。

Baum-Welch算法
刚才提到了，非监督学习问题是为了计算模型参数 λ <script type="math/tex" id="MathJax-Element-273">\lambda</script>，使得在该参数下 P(O|λ) <script type="math/tex" id="MathJax-Element-274">P(O|\lambda)</script>的概率最大。这个问题便是我们的极大似然估计，但 P(O|λ) <script type="math/tex" id="MathJax-Element-275">P(O|\lambda)</script>并非孤立的存在，其背后与隐含状态相联系。这句话应该怎么理解呢，在海藻模型中，如我们观测到某一海藻序列 O={"dry","damp","soggy"} <script type="math/tex" id="MathJax-Element-276">O=\{"dry","damp","soggy"\}</script>，但是什么决定了海藻的湿度情况呢，很明显天气的因素占有很大的一部分，因此盲人在对海藻模型进行建模时，就把隐含的天气转移状态给考虑进去了。正如双硬币模型中，由于实习生b的失误，每组数据我们并不清楚是A掷的还是B掷的，遇到信息缺失的情况，就导致了用单纯的极大似然估计求导法是无法得到解析解的。

假设给定训练数据中包含S个长度为T的观测序列 {O1,O2,...,OS} <script type="math/tex" id="MathJax-Element-277">\{O_1,O_2,...,O_S\}</script>而没有对应的状态序列，目标是学习隐马尔可夫模型 λ=(A,B,π) <script type="math/tex" id="MathJax-Element-278">\lambda= (A,B,\pi)</script>的参数。我们将观测序列数据看作观测数据O，状态序列数据看作不可观测的隐数据 I <script type="math/tex" id="MathJax-Element-279">I</script>，那么隐马尔可夫模型事实上是一个含有隐变量的概率模型

1.确定完全数据的对数似然函数
所有观测数据写成 O=(o1,o2,...,oT) <script type="math/tex" id="MathJax-Element-281">O=(o_1,o_2,...,o_T)</script>，所有隐数据写成 I=(i1,i2,...,iT) <script type="math/tex" id="MathJax-Element-282">I=(i_1,i_2,...,i_T)</script>，完全数据是 (O,I)=(o1,o2,...,oT,i1,i2,...,iT) <script type="math/tex" id="MathJax-Element-283">(O,I)=(o_1,o_2,...,o_T,i_1,i_2,...,i_T)</script>。完全数据的对数似然函数是 logP(O,I|λ) <script type="math/tex" id="MathJax-Element-284">\log P(O,I|\lambda)</script>。

2.EM算法的E步：求Q函数 Q(λ,λ^) <script type="math/tex" id="MathJax-Element-285">Q(\lambda,\hat \lambda)</script>

3.EM算法的M步：极大化Q函数 Q(λ,λ^) <script type="math/tex" id="MathJax-Element-298">Q (\lambda,\hat \lambda)</script>求模型参数 A,B,π <script type="math/tex" id="MathJax-Element-299">A,B,\pi</script>
由于要极大化的参数在上式中单独地出现在3个项中，所以只需要对各项分别极大化。

（1）上式中的第一项可以写成：

（2）上式中的第二项可以写成

（3）上式中的第三项可以写成

算法实际物理含义

EM算法中M步的各公式的难点在于如何求得这些概率，如 aij <script type="math/tex" id="MathJax-Element-319">a_{ij}</script>该公式分子分母上的联合概率如何计算。其实在我看来，对隐马尔可夫模型中的各种概率计算最后均是映射到节点上去做计算。当然，我们先来观察由EM算法推导出的参数计算公式。

观察式子 aij <script type="math/tex" id="MathJax-Element-320">a_{ij}</script>和 bj(k) <script type="math/tex" id="MathJax-Element-321">b_j(k)</script>，你会发现不管是分子，还是分母，它们都是概率计算，只不过对应的一些状态不一样。具体以 aij <script type="math/tex" id="MathJax-Element-322">a_{ij}</script>举例，如在 aij <script type="math/tex" id="MathJax-Element-323">a_{ij}</script>的分母上计算式子 P(O,it=i,it+1=j|λ^) <script type="math/tex" id="MathJax-Element-324">P(O,i_t=i,i_{t+1}=j |\hat\lambda)</script>，仔细想想，我们在计算什么的时候，有遇到过类似的式子？其实在阐述隐马尔可夫模型的第一个概率计算问题时，我们就做过类似的求解。概率计算是为了计算 P(O|λ) <script type="math/tex" id="MathJax-Element-325">P(O|\lambda)</script>的概率，但我们是把式子扩展为 ∑IP(O,I|λ) <script type="math/tex" id="MathJax-Element-326">\sum_I P(O,I|\lambda)</script>进行计算的。即我们需要在任何隐藏状态序列下求出 P(O|λ) <script type="math/tex" id="MathJax-Element-327">P(O|\lambda)</script>的概率。由此我们用前向算法和后向算法来求解该问题，很好的把概率计算问题，映射到了物理节点上去做计算，并且借助物理节点存储中间变量的特性大大的简化了算法的复杂度。

同样地， P(O,it=i,it+1=j|λ^) <script type="math/tex" id="MathJax-Element-328">P(O,i_t=i,i_{t+1}=j |\hat\lambda)</script>不就可以看成是对

于是我们就有了前向后向算法中 ξt(i,j) <script type="math/tex" id="MathJax-Element-332">\xi_t(i,j)</script>的定义了，它的定义式为：

还记得前向算法和后向算法是如何定义中间节点的嘛，为了计算 P(O|λ) <script type="math/tex" id="MathJax-Element-343">P(O|\lambda)</script>，我们给每一个t时刻的隐含状态节点定义了实际的物理含义，即把它们命名为 αt(i) <script type="math/tex" id="MathJax-Element-344">\alpha_t(i)</script>和 βt(i) <script type="math/tex" id="MathJax-Element-345">\beta_t(i)</script>，两个中间变量分别从两边进行有向边加权和有向边汇聚，形成一种递归结构，并且由此不断传播至两端，对任意t=1时刻，和t=T时刻，分别进行累加就能求得 P(O|λ) <script type="math/tex" id="MathJax-Element-346">P(O|\lambda)</script>，我们还举出了一个小例子，来论证前向算法和后向算法只要满足有向边加权和有向边汇聚就能得到算法的一致性，今天我们根据前向后向算法做进一步的例子推广，从而真正理解 ξt(i,j) <script type="math/tex" id="MathJax-Element-347">\xi_t(i,j)</script>的物理含义。

书中利用前向概率和后向概率的定义可以将观测序列概率 P(O|λ) <script type="math/tex" id="MathJax-Element-348">P(O|\lambda)</script>统一写成

假设我们从t =T-1时刻开始递推，那么上述式子，把t=T-1代入得

我们把问题用节点图表示

从A=1,B=1，向上传播得C=840元，从C=1，向下传播得A+B=840元。而 P(O|λ) <script type="math/tex" id="MathJax-Element-361">P(O|\lambda)</script>这个式子告诉我们，从中间某一层节点，从C出发到达该层和从A,B出发到达该层，等到的两个中间变量进行相乘，且对该层所有节点进行累加完毕后，获得的资源总数不变，即840元。不信，我们分别对第二层和第三层的节点计算一遍，如下图

由此，我们明确了 ξt(i,j) <script type="math/tex" id="MathJax-Element-362">\xi_t(i,j)</script>的物理含义，它只是从t时刻出发，由前向算法导出的中间节点 Si <script type="math/tex" id="MathJax-Element-363">S_i</script>和从t+1时刻出发，由后向导出的中间节点 Sj <script type="math/tex" id="MathJax-Element-364">S_j</script>，且节点 Si和Sj <script type="math/tex" id="MathJax-Element-365">S_i和S_j</script>中间还有一条加权有向边的关系 aijbj(ot+1) <script type="math/tex" id="MathJax-Element-366">a_{ij}b_j(o_{t+1})</script>，所以我们得

算法（Baum-Welch算法）

输入：观测数据 O=(o1,o2,...,oT) <script type="math/tex" id="MathJax-Element-375">O=(o_1,o_2,...,o_T)</script>
输出：隐马尔可夫模型参数
（1）初始化
对 n=0 <script type="math/tex" id="MathJax-Element-376">n=0</script>，选取 a(0)ij,bj(k)(0),π(0)i <script type="math/tex" id="MathJax-Element-377">a_{ij}^{(0)},b_j(k)^{(0)},\pi_i^{(0)}</script>，得到模型 λ(0)=(A(0),B(0),π(0)) <script type="math/tex" id="MathJax-Element-378">\lambda^{(0)} = (A^{(0)},B^{(0)},\pi^{(0)})</script>
（2）递推，对 n=1,2,..., <script type="math/tex" id="MathJax-Element-379">n=1,2,...,</script>

a(n+1)ij=∑T−1t=1ξt(i,j)∑T−1t=1γt(i)

<script type="math/tex; mode=display" id="MathJax-Element-380">a_{ij}^{(n+1)} = \frac{\sum_{t=1}^{T-1}\xi_t(i,j)}{\sum_{t=1}^{T-1}\gamma_t(i)}</script>

b(n+1)j(k)=∑T−1t=1,ot=vkγt(j)∑Tt=1γt(j)

<script type="math/tex; mode=display" id="MathJax-Element-381">b_{j}^{(n+1)}(k) = \frac{\sum_{t=1,o_t=v_k}^{T-1}\gamma_t(j)}{\sum_{t=1}^{T}\gamma_t(j)}</script>

π(n+1)i=γ1(i)

<script type="math/tex; mode=display" id="MathJax-Element-382">\pi_i^{(n+1)}=\gamma_1(i)</script>
（3）终止，得到模型参数 λ(n+1)=(A(n+1),B(n+1),π(n+1)) <script type="math/tex" id="MathJax-Element-383">\lambda^{(n+1) } =(A^{(n+1)},B^{(n+1)},\pi^{(n+1)})</script>

理论分析总算完毕了，简单总结一下前向后向算法，首先隐马尔可夫模型参数的估计问题是一个隐藏变量的极大似然估计，因此我们用到了EM算法来解决上述参数估计问题，从EM算法中，求得Q函数，从而能够对Q函数进行求偏导，得到极大似然函数的极值，求偏导算出了参数估计的公式，与先前 λ^ <script type="math/tex" id="MathJax-Element-384">\hat\lambda</script>参数产生了关系，并进一步需要计算大量联合概率，而联合概率的计算巧妙的使用了节点图的各种性质，用中间变量降低了节点计算的复杂度，导出了对计算有帮助的定义，方便参数求解。

Code Time

Baum-Welch算法的Python实现

接着前文hmm.py和test.py文件继续添加。

1.在hmm.py中添加baum_welch_train算法

def baum_welch_train(self, observations, criterion=0.05):
    n_states = self.A.shape[0]
    # 观察序列的长度T
    n_samples = len(observations)

    done = False
    while not done:
        # alpha_t(i) = P(o_1,o_2,...,o_t,q_t = s_i | hmm)
        # Initialize alpha
        # 获得所有前向传播节点值 alpha_t(i)
        alpha = self._forward(observations)

        # beta_t(i) = P(o_t+1,o_t+2,...,o_T | q_t = s_i , hmm)
        # Initialize beta
        # 获得所有后向传播节点值 beta_t(i)
        beta = self._backward(observations)

        # 计算 xi_t(i,j) -> xi(i,j,t)
        xi = np.zeros((n_states, n_states, n_samples - 1))
        # 在每个时刻
        for t in range(n_samples - 1):
            # 计算P(O | hmm)
            denom = sum(alpha[:, -1])
            for i in range(n_states):
                # numer[1,:] = 行向量，alpha[i,t]=实数，slef.A[i,:] = 行向量
                # self.B[:,observations[t+1]].T = 行向量,beta[:,t+1].T = 行向量
                numer = alpha[i, t] * self.A[i, :] * self.B[:, observations[t + 1]].T * beta[:, t + 1].T
                xi[i, :, t] = numer / denom

            # 计算gamma_t(i) 就是对j进行求和
            gamma = np.sum(xi, axis=1)
            # need final gamma elements for new B
            prod = (alpha[:, n_samples - 1] * beta[:, n_samples - 1]).reshape((-1, 1))
            # 合并T时刻的节点
            gamma = np.hstack((gamma, prod / np.sum(prod)))
            # 列向量
            newpi = gamma[:, 0]
            newA = np.sum(xi, 2) / np.sum(gamma[:, :-1], axis=1).reshape((-1, 1))
            newB = np.copy(self.B)

            # 观测状态数
            num_levels = self.B.shape[1]
            sumgamma = np.sum(gamma, axis=1)
            for lev in range(num_levels):
                mask = observations == lev
                newB[:, lev] = np.sum(gamma[:, mask], axis=1) / sumgamma

            if np.max(abs(self.pi - newpi)) < criterion and \
                            np.max(abs(self.A - newA)) < criterion and \
                            np.max(abs(self.B - newB)) < criterion:
                done = 1
            self.A[:], self.B[:], self.pi[:] = newA, newB, newpi

2.在hmm.py中添加模拟序列生成函数

def simulate(self, T):
    def draw_from(probs):
        # np.random.multinomial 为多项式分布，1为实验次数，类似于投掷一枚骰子，丢出去是几，probs每个点数的概率，均为1/6
        # 给定行向量的概率，投掷次数为1次，寻找投掷的点数
        return np.where(np.random.multinomial(1, probs) == 1)[0][0]

    observations = np.zeros(T, dtype=int)
    states = np.zeros(T, dtype=int)
    states[0] = draw_from(self.pi)
    observations[0] = draw_from(self.B[states[0], :])
    for t in range(1, T):
        states[t] = draw_from(self.A[states[t - 1], :])
        observations[t] = draw_from(self.B[states[t], :])
    return observations, states

回到海藻模型，我们可以用这样一串代码完成Baum-Welch算法的训练，并且评估其准确率。

3.在test.py中添加测试代码

# 参数估计
observations_data, states_data = h.simulate(100)
guess = hmm.HMM(array([[0.33, 0.33, 0.34],
                       [0.33, 0.33, 0.34],
                       [0.33, 0.33, 0.34]]),
                array([[0.25, 0.25, 0.25, 0.25],
                       [0.25, 0.25, 0.25, 0.25],
                       [0.25, 0.25, 0.25, 0.25]]),
                array([0.7, 0.15, 0.15])
                )
guess.baum_welch_train(observations_data)
# 预测问题
states_out = guess.state_path(observations_data)[1]
p = 0.0
for s in states_data:
    if next(states_out) == s: p += 1

print(p / len(states_data))

经过多次测试，本算法的预测准确率在0.3~0.5。可见隐马尔可夫模型的参数估计的准确率还没有到令人满意的程度。

参考文献

1. EM算法及其推广学习笔记
2. 隐马尔可夫学习笔记（一）
3. 李航. 统计学习方法[M]. 北京：清华大学出版社，2012