设观察序列为 $X$，隐藏序列为 $Z$，且隐藏序列有 $K$ 种可选状态。设隐藏状态的初始概率分布为 $p(z_1 \mid \pi)$，其中 $\pi_k$ 是 $z_1$ 取第 $k$ 个状态（$S_k$）的概率，即 $p(z_1 = S_k) = \pi_k$。
设隐藏状态的转换概率矩阵为 $A$，即 $p(z_{i+1} = S_j \mid z_i = S_k) = A_{kj}$。
设隐藏状态到观察状态的概率分布为 $p(x_i \mid z_i, \theta)$。如果观察序列有 $D$ 种状态，用矩阵表示的话，$\theta$ 是个 $K \times D$ 的矩阵（注意，这里假设 $p(x_i \mid z_i, \theta)$ 是离散的多项式分布，后面的推导也基于这个假设）。
下面看如何用 EM 算法求得模型参数 $\pi, A, \theta$。
先看隐藏状态和观察状态的联合概率分布：
\begin{align} p(X,Z \mid \pi, A, \theta) & = p(x_1, x_2, \dots, x_N, z_1, z_2, \dots, z_N \mid \pi, A, \theta)\nonumber\\ & = p(z_1 \mid \pi)\prod_{n=2}^N{p(z_n \mid z_{n-1}, A)}\prod_{m=1}^N{p(x_m \mid z_m, \theta)}. \end{align}
上面是说，在知道模型参数的情况下，可以通过这个公式得出一组数据的概率。这里的“一组数据”是指一对儿 $(X, Z)$。
反过来，如果手里有成对儿的 $(X, Z)$，也就是完整的数据，模型参数是可以通过最大似然（ML）（或者最大后验概率）方法求得的：
\begin{align} argmax_{\pi, A, \theta}\log{{p(X, Z \mid \pi, A, \theta)}}. \end{align}
但是，在只知道部分数据（$X$），而另一部分数据（$Z$）未知的情况下如何求模型参数呢？这里的难题是，不知道 $Z$ 就无法根据 $p(X, Z)$ 得到概率，也就无法利用优化算法求解。
（以下为书写方便，用 $\Phi$ 代替 $\pi, A, \theta$）
一个思路是，既然 $Z$ 未知，那就把 $Z$ 所有可能的情况都考虑进来。但是 $Z$ 都有哪些情况呢？因为不知道模型参数，也就无法知道 $Z$ 的分布。
这就是两难的局面：
此时，EM to the rescure：既然不能同时得到 $Z$ 和 模型参数，那就迭代交替计算它们。而要开始这个迭代就需要一方先有值，这里一般先随机初始化模型参数。然后迭代的步骤是：
上面第二步中具体求解的公式是： \begin{align} argmax_{\Phi}\sum_{Z}{p(Z \mid X, \Phi^{old})\log{{p(X, Z \mid \Phi)}}}. \end{align} 注意，在当前步骤里 $\Phi^{old}$ 是不变的，跟 $\Phi$ 不一样。$\Phi$ 是待求解的变量。
这个迭代过程在模型参数不再变化（或者 $p(X \mid \Phi) = \sum_Z{p(X, Z \mid \Phi)}$ 不再变化）的时候停止。
下面看如何把这个方法应用在 HMM 上。
上面式子中的 $\sum_{Z}{p(Z \mid X, \Phi^{old})\log{{p(X, Z \mid \Phi)}}}$ 是关于 $\log{p(X, Z \mid \Phi)}$ 的期望，就是 EM 中 E。这个 E 如何表现为关于 $\pi, A, \theta$ 的函数？
把上面的等式（1）代入进来可以看到： \begin{align} \sum_{Z}{p(Z \mid X, \Phi^{old})\log{{p(X, Z \mid \Phi)}}} & = \sum_{Z}{p(Z \mid X, \Phi^{old})\left(\log{p(z_1 \mid \pi)} + \sum_{n=2}^N{\log{p(z_n \mid z_{n-1}, A)}} + \sum_{m=1}^N{\log{p(x_m \mid z_m, \theta)}}\right)}\nonumber\\ & = \sum_{Z}{p(Z^{old})\log{p(z_1 \mid \pi)}} + \sum_{Z}{p(Z^{old})\sum_{n=2}^N{\log{p(z_n \mid z_{n-1}, A)}}} + \sum_{Z}{p(Z^{old})\sum_{m=1}^N{\log{p(x_m \mid z_m, \theta)}}} \end{align}
（上面用 $p(Z^{old})$ 代替了 $p(Z \mid X, \Phi^{old})$）
先看右手边求和的第一项 $\sum_{Z}{p(Z^{old})\log{p(z_1 \mid \pi)}}$。这是把所有可能的 $Z$ 都列出来，然后加和：
\begin{align} \sum_{Z}{p(Z^{old})\log{p(z_1 \mid \pi)}} & = p(z_1 = S_1, z_2 = S_1, \dots, z_N = S_1)\log{p(z_1 = S_1 \mid \pi)} \\ & + p(z_1 = S_1, z_2 = S_1, \dots, z_N = S_2)\log{p(z_1 = S_1 \mid \pi)} \nonumber \\ & + \cdots \nonumber \\ & + p(z_1 = S_1, z_2 = S_1, \dots, z_N = S_K)\log{p(z_1 = S_1 \mid \pi)} \nonumber \\ & + p(z_1 = S_1, z_2 = S_2, \dots, z_N = S_1)\log{p(z_1 = S_1 \mid \pi)} \nonumber \\ & + \cdots \nonumber \\ & + p(z_1 = S_1, z_2 = S_K, \dots, z_N = S_K)\log{p(z_1 = S_1 \mid \pi)} \nonumber \\ & + p(z_1 = S_2, z_2 = S_1, \dots, z_N = S_1)\log{p(z_1 = S_2 \mid \pi)} \nonumber \\ & + \cdots \nonumber \\ & + p(z_1 = S_2, z_2 = S_K, \dots, z_N = S_K)\log{p(z_1 = S_2 \mid \pi)} \nonumber \\ & + \cdots \nonumber \\ & + p(z_1 = S_K, z_2 = S_1, \dots, z_N = S_1)\log{p(z_1 = S_K \mid \pi)} \nonumber \\ & + \cdots \nonumber \\ & + p(z_1 = S_K, z_2 = S_K, \dots, z_N = S_K)\log{p(z_1 = S_K \mid \pi)}. \end{align}
可以看到 $log$ 中的 $z_1$ 只有 $K$ 种变化，而前面的 $p(Z^{old})$ 有 $K^N$ 种变化。可以把 $K$ 种 $log$ 提取公因子：
\begin{align} \sum_{Z}{p(Z^{old})\log{p(z_1 \mid \pi)}} & = \sum_{z_2,\dots,z_N}{p(z_1 = S_1, z_2, \dots, z_N)\log{p(z_1 = S_1 \mid \pi)}}\nonumber\\ & + \sum_{z_2,\dots,z_N}{p(z_1 = S_2, z_2, \dots, z_N)\log{p(z_1 = S_2 \mid \pi)}}\nonumber\\ & + \cdots\nonumber\\ & + \sum_{z_2,\dots,z_N}{p(z_1 = S_K, z_2, \dots, z_N)\log{p(z_1 = S_K \mid \pi)}}. \end{align}
而上式中的 $\sum_{z_2,\dots,z_N}{p(z_1, z_2, \dots, z_N)}$ 就是 $z_1$ 的边缘分布（Marginal Distribution）$p(z_1 \mid X, \Phi^{old})$。设 $\gamma(z_1) = \sum_{z_2,\dots,z_N}{p(z_1, z_2, \dots, z_N)}$，则上式等于：
\begin{align} \sum_{Z}{p(Z^{old})\log{p(z_1 \mid \pi)}} & = \gamma(z_1 = S_1)\log{p(z_1 = S_1 \mid \pi)}\nonumber\\ & + \gamma(z_1 = S_2)\log{p(z_1 = S_2 \mid \pi)}\nonumber\\ & + \cdots\nonumber\\ & + \gamma(z_1 = S_K)\log{p(z_1 = S_K \mid \pi)}. \end{align}
再看 $log$ 项，其中 $p(z_1 = S_1 \mid \pi)$ 就是 $\pi_1$，$p(z_1 = S_K \mid \pi)$ 就是 $\pi_K$。所以上式进一步等于：
\begin{align} \sum_{Z}{p(Z^{old})\log{p(z_1 \mid \pi)}} & = \gamma(z_1 = S_1)\log{\pi_1}\nonumber\\ & + \gamma(z_1 = S_2)\log{\pi_2}\nonumber\\ & + \cdots\nonumber\\ & + \gamma(z_1 = S_K)\log{\pi_K}\nonumber\\ & = \sum_{k=1}^K{\gamma(z_{1k})\log{\pi_k}}. \end{align}
上式中用 $\gamma(z_{1k})$ 表示 $\gamma(z_{1} = S_k)$。
这样等式（4）右边第一项表示成了关于 $\pi$ 的函数。下面看右边第二项。
第二项同样可以使用展开然后提取公因子的方式转化为关于 $A$ 的函数。只是需要注意，$p(z_n \mid z_{n-1})$ 涉及两个 $z$，所以提取公因子后得到的是相邻两个 $z$ 的边缘分布。
定义 $\xi(z_{n-1},z_n) = p(z_{n-1}, z_n \mid X, \Phi^{old})$，并且注意到 $\log{p(z_n = S_k \mid z_{n-1} = S_j)} = log{A_{jk}}$，等式（4）右边第二项等于：
$$ \sum_{Z}{p(Z^{old})\sum_{n=2}^N{\log{p(z_n \mid z_{n-1}, A)}}} = \sum_{n=2}^N{\sum_{j=1}^K{\sum_{k=1}^K{\xi(z_{(n-1)j},z_{nk})\log{A_{jk}}}}}. $$
同样的方法可以得出等式（4）右边第三项等于：
$$ \sum_{Z}{p(Z^{old})\sum_{m=1}^N{\log{p(x_m \mid z_m, \theta)}}} = \sum_{n=1}^N{\sum_{k=1}^K{\sum_{d=1}^D{\gamma(z_{nk})\log{\theta_{kd}}}}}. $$
现在等式（4）可以重新表示为：
\begin{equation} \sum_{Z}{p(Z \mid X, \Phi^{old})\log{{p(X, Z \mid \Phi)}}} = \sum_{k=1}^K{\gamma(z_{1k})\log{\pi_k}} + \sum_{n=2}^N{\sum_{j=1}^K{\sum_{k=1}^K{\xi(z_{(n-1)j},z_{nk})\log{A_{jk}}}}} + \sum_{n=1}^N{\sum_{k=1}^K{\sum_{d=1}^D{\gamma(z_{nk})\log{\theta_{kd}}}}}. \end{equation}
不过要计算 $argmax_{\pi, A, \theta}$ 依然缺少关于 $Z$ 的两类边缘分布 $\gamma(z_{nk})$ 和 $\xi(z_{(n-1)j}, z_{nk})$，下面看如何计算。
sum-product 是树形的概率图中计算变量边缘分布的通用算法，自然也可以应用在 HMM 中。sum-product 利用变量之间的条件无关性（conditional-independent）来避免重复计算，把看起来指数时间的复杂度恢复成多项式时间。
下面用个简单的例子说明一下。
假设要求解的边缘分布是 $p(z_3) = \sum_{z_1,z_2,z_4,z_5}{P(z_1, z_2, z_3, z_4, z_5)}$，并且 $P$ 可以分解成 $P(z_1, z_2, z_3, z_4, z_5) = p_1(z_1, z_2)p_2(z_2, z_3)p_3(z_3, z_4)p_4(z_4, z_5)$。
如果简单的尝试每组可能的 $Z$ 的话，加和的量是 $K^5$（假设每个 $z_i$ 有 $K$ 种状态）。但是利用 $P$ 可分解的事实可以把上式重新组织成
$$p(z_3) = \left(\left(\sum_{z_1}{p_1(z_1,z_2)}\right)\sum_{z_2}{p_2(z_2,z_3)}\right)\left(\sum_{z_4}{p_3(z_3,z_4)}\left(\sum_{z_5}{p_4(z_4,z_5)}\right)\right)$$.
这是一个因式分解的过程。但是其中 $\sum_{z_1}{p_1(z_1,z_2)}$ 并不是一个固定的值，而是关于 $z_2$ 的函数。设
$$f_1(z_2) = \sum_{z_1}{p_1(z_1,z_2)}$$
，上式简化为 $$ p(z_3) = \left(\sum_{z_2}{p_2(z_2,z_3)f_1(z_2)}\right)\left(\sum_{z_4}{p_3(z_3,z_4)}\left(\sum_{z_5}{p_4(z_4,z_5)}\right)\right). $$
上面 $\sum_{z_2}{p_2(z_2,z_3)f_1(z_2)}$ 是关于 $z_3$ 的函数。类似的设
$$f_2(z_3) = \sum_{z_2}{p_2(z_2,z_3)f_1(z_2)}$$
，$p(z_3)$ 可以进一步简化为
$$ p(z_3) = f_2(z_3)\left(\sum_{z_4}{p_3(z_3,z_4)}\left(\sum_{z_5}{p_4(z_4,z_5)}\right)\right). $$
$z_4, z_5$ 也做相应的操作，$p(z_3)$ 最终简化为
$$ p(z_3) = f_2(z_3)f_4(z_3) $$
，其中 $f_4(z_3) = \sum_{z_4}{p_3(z_3,z_4)f_5(z_4)}$ 以及 $f_5(z_4) = \sum_{z_5}{p_4(z_4,z_5)}$。
那如此分解之后的加和项有多少呢？
可以看到计算 $f_1(z_2)$ 时对每个 $z_2$ 有 $K$ 次加和，所以 $f_1$ 会有 $K^2$ 次加和。相应的 $f_2(z_3)$ 对每个 $z_3$ 也是 $K$ 次加和（假设需要知道每个 $z_3$ 的边缘分布），$f_2$ 同样是 $K^2$ 次加和。$f_4, f_5$ 一样，所以一共是 $4K^2$ 次加和（实际是 $(N-1)K^2$ 次）。
所以这个方法把随变量个数成指数级增长的时间复杂度降为了与变量个数成线性相关。
上面的例子是树形概率图中最简单的形式——链式（后面将看到 HMM 正是链式结构）。计算 $p(z_3)$ 的步骤可以是：
前两步是正向的 $f_1 \rightarrow f_2$ ，下面两步是反向的 $f_4 \leftarrow f_5$，所以针对 HMM 的 sum-product 算法又叫 forward-backward。
上面不是 sum-product 的完整介绍，因为并没有涉及真正的树形的情况，比如一个变量出现在三个或以上因子里或者一个因子含有三个或以上变量的情况。但是总体思路都是提取公因子，减少计算量。详细的可以参考 PRML 8.4 小节。
回过头来看，待计算的边缘分布是 $\gamma(z_{nk}) = p(z_n = S_k \mid X, \Phi^{old})$ （后面将省略 $\Phi^{old}$）。但是这个分布不能直接利用 forward-backward 算法，解决方法是先计算 $p(z_n = S_k, X)$，然后再计算 $\gamma{z_{nk}}$。因为
\begin{align} \gamma(z_{nk}) & = \frac{p(z_n = S_k, X)}{\sum_{j=1}^K{p(z_n = S_j, X)}}\nonumber\\ & = \frac{p(z_n = S_k \mid X)p(X)}{\sum_{j=1}^K{p(z_n = S_j \mid X)p(X)}}\nonumber\\ & = \frac{p(z_n = S_k \mid X)}{\sum_{j=1}^K{p(z_n = S_j \mid X)}}\nonumber\\ & = p(z_n = S_k \mid X). \end{align}
然后看如何计算 $p(z_n, X)$：
\begin{align} p(z_n, X) &= \sum_{z_1, z_2, \dots, z_{n-1}, z_{n+1}, \dots, z_N}p(z_1,z_2,\dots,z_N, X)\nonumber\\ &=\sum_{z_1, z_2, \dots, z_{n-1}, z_{n+1}, \dots, z_N}\left(p(z_1)p(x_1 \mid z_1)\right)\left(p(z_2 \mid z_1)p(x_2 \mid z_2)\right)\cdots\left(p(z_N \mid z_{N-1})p(x_N \mid z_N)\right)\nonumber\\ &=\sum_{z_1}p(z_1)p(x_1 \mid z_1)p(z_2 \mid z_1)p(x_2 \mid z_2)\sum_{z_2}p(z_3\mid z_2)p(x_3 \mid z_2)\cdots\sum_{z_{N-1}}p(z_N\mid z_{N-1})p(x_N \mid z_{N-1}). \end{align}
设：
\begin{align} h(z_1) &= p(z_1)p(x_1 \mid z_1)\nonumber\\ g_n(z_{n}, z_{n+1}) &= p(z_{n+1} \mid z_{n})p(x_{n+1} \mid z_{n+1}) \end{align}
$p(z_n, X)$ 表示为：
\begin{align} p(z_n, X) = \left(\sum_{z_1}p(z_1)h_1(z_1)g_1(z_1,z_2)\sum_{z_2}g_2(z_2, z_3)\cdots\sum_{z_{n-1}}g_{n-1}(z_{n-1},z_n)\right)\left(\sum_{z_n}g_n(z_n,z_{n+1})\cdots\sum_{z_{N-1}}g_{N-1}(z_{N-1}, Z_N)\right). \end{align}
这个结构与上一节中的例子相似，可以使用 forward-backward 算法计算。之后再用公式（11）计算 $\gamma(r_{nk})$。
第二类边缘分布是 $\xi(z_{n-1}, z_n) = p(z_{n-1}, z_n \mid X, \Phi^{old})$。这个也是先计算 $p(z_{n-1}, z_n, X)$ （省略 $\Phi^{old}$）：
\begin{align} p(z_{n-1}, z_n, X) = \left(\sum_{z_1}p(z_1)h_1(z_1)g_1(z_1,z_2)\sum_{z_2}g_2(z_2, z_3)\cdots\sum_{z_{n-2}}g_{n-2}(z_{n-2},z_{n-1})\right)\sum_{z_{n-1}}g_{n-1}(z_{n-1},z_n)\left(\sum_{z_n}g_n(z_n,z_{n+1})\cdots\sum_{z_{N-1}}g_{N-1}(z_{N-1}, Z_N)\right)\nonumber. \end{align}
上式中左边括号内用 forward 算法计算，右边括号内容用 backward 算法计算，然后再与 $\sum_{z_{n-1}}g_{n-1}(z_{n-1},z_n)$ 相乘。有了 $p(z_{n-1}, z_n, X)$ 之后再计算每一对儿 $p(z_{n-1} = S_i, z_{n} = S_j, X)$ 在 $\sum_{t=1}^K{\sum_{u=1}^K{p(z_{n-1} = S_t, z_{n} = S_u, X)}}$ 中的比例即可，$p(X)$ 会被消掉。
前面给出了 $\gamma(z_{nk})$ 和 $\xi(z_{n-1},z_n)$ 两个边缘分布，现在式子（10）就完整了。剩下的就是寻找 $argmax_{\pi,A,\theta}$。
这个优化问题可以简单的利用拉格朗日乘子法求解。下面给出对 $\pi_i$ 的推导过程。
式子（10）中只有第一项跟 $\pi$ 有关，并且约束条件是 $\sum_i^K{\pi_i} = 1$。引入拉格朗日乘子 $\lambda$ 得：
\begin{align} \sum_{k=1}^K{\gamma(z_{1k})\log{\pi_k}} - \lambda(\sum_k^K{\pi_k} - 1). \end{align}
对每个 $\pi_k$ 求导并等于 $0$，得：
\begin{align} \gamma(r_{11})\frac{1}{\pi_1} - \lambda &= 0\\ \gamma(r_{12})\frac{1}{\pi_2} - \lambda &= 0\\ \cdots\nonumber\\ \gamma(r_{1K})\frac{1}{\pi_K} - \lambda &= 0\\ \end{align}
各自乘以 $\pi_k$ 得 $\gamma(r_{1k}) - \lambda\pi_k = 0$，然后 $1\dots K$ 的式子相加，得：
$$ \sum_{k=1}^K{r_{1k}} = \lambda(\sum_{k=1}^K{\pi_k}). $$
因为约束条件 $\sum_i^K{\pi_i} = 1$ 所以，上式是 $\lambda = \sum_{k=1}^K{r_{1k}}$ 。带入（16）（17）（18）等等，得：
$$ \pi_k = \frac{\gamma(r_{1k})}{\sum_{j=1}^K{\gamma(r_{1j})}}. $$
用类似的方法可以求得 $A$ 和 $\theta$ 的最优解：
\begin{align} A_{jk} &= \frac{\sum_{n=2}^N{\xi(z_{(n-1)j,z_{nk}})}}{\sum_{l=1}^K{\sum_{n=2}^N{\xi(z_{(n-1)j}, z_{nl})}}}\\ \theta_{ki} &= \frac{\sum_{n=1}^N{\gamma(z_(nk))x_{ni}}}{\sum_{n=1}^N{\gamma(z_(nk))}} \end{align}
（其中，$x_{ni}$ 表示第 $n$ 个 $x$ 是否是第 $i$ 个状态，等于就是 1，否则是 0，一共 D 个状态）。
注意，$\theta$ 的最优解表达式是基于 $p(x_i\mid z_i)$ 服从离散的多项式分布这个假设导出的。其他的分布会有不同的最优解表达式。
上面给出了针对 HMM 的 EM 算法推导过程。具体实现可以参考 HMM-import.ipynb 中的 HMMem 函数（其中支持对多组观察序列的处理）。
需要注意的是，EM 算法不保证得到全局最优解，需要多次尝试不同的参数初始值。