馬爾科夫鏈是描述狀態轉換的隨機過程，該過程具備“無記憶”的性質：即當前時刻$t$的狀態$s_t$的概率分佈只由前一時刻$t-1$的狀態$s_{t-1}$決定，與時間序列中$t-1$時刻之前的狀態無關。定義馬爾科夫鏈的轉移矩陣為$A$，有$$A_{ij}=p\left(s_{t}=j | s_{t-1}=i\right),\text{ }s_{t} | s_{t-1} \sim \operatorname{Discrete}\left(A_{s_{t-1}, :}\right)$$容易看出矩陣$A$的每行之和為1，給定一個起始狀態$s_1$(也可通過某個分佈產生)，可通過從上述分佈中抽樣生成序列$\left(s_{1}, s_2, \dots, s_{t}\right)$。

模型定義

隱馬爾科夫模型HMM假設觀測是從一個隱藏的馬爾科夫狀態序列生成的，如下圖所示：

• 不失一般性，假設共有$S$個離散狀態，即$s \in \{1,2,\dots,S\}$
• 不失一般性，假設共有$X$種觀測，即$x \in \{1,2,\dots,X\}$
• 定義初始的狀態分佈為$\vec{\pi}=(\pi_1,\pi_2,\dots,\pi_S)$，即$s_{1} \sim \operatorname{Discrete}(\vec{\pi})$
• 定義狀態$s$的轉移矩陣為$\mathcal{A}$，則$\mathcal{A}$為$S$行$S$列的矩陣，即$s_{i} |\left\{s_{i-1}=k^{\prime}\right\} \sim \operatorname{Discrete}\left(\mathcal{A}_{k^{\prime}, :}\right)$
• 定義觀測$x$的生成概率為$\mathcal{B}$，則$\mathcal{B}$為$S$行$X$列的矩陣，即$x_{i} |\left\{s_{i}=k^{\prime}\right\} \sim \operatorname{Discrete}\left(\mathcal{B}_{k^{\prime}, :}\right)$

則HMM的求解可歸納為以下兩個問題：

1. 訓練：已知觀測序列$\vec{o}=(x_1,x_2,\dots,x_T)$，使用最大似然估計計算引數$\vec{\pi},\mathcal{A},\mathcal{B}$，即$$\vec{\pi}_{\mathrm{ML}}, \mathcal{A}_{\mathrm{ML}}, \mathcal{B}_{\mathrm{ML}}=\arg \max _{\vec{\pi},\mathcal{A},\mathcal{B}} p\left(\vec{o} | \vec{\pi},\mathcal{A},\mathcal{B}\right)$$
2. 預測：已知觀測序列$\vec{o}$以及引數$\vec{\pi},\mathcal{A},\mathcal{B}$，估計生成觀測序列$\vec{o}$的最可能的隱藏狀態序列$\vec{s}=(s_{1}, \ldots, s_{T})$，即$$s_{1}, \ldots, s_{T}=\arg \max _{\vec{s}} p\left(\vec{s} | \vec{o}, \vec{\pi},\mathcal{A},\mathcal{B}\right)$$

引數估計

要解決問題1，首先看一下如何估計$p\left(\vec{o} | \vec{\pi},\mathcal{A},\mathcal{B}\right)$，根據概率公式，有$$\begin{aligned} p(\vec{o} | \vec{\pi}, \mathcal{A}, \mathcal{B}) &=\sum_{s_{1}=1}^{S} \cdots \sum_{s_{T}=1}^{S} p\left(\vec{o}, s_{1}, \ldots, s_{T} | \vec{\pi}, \mathcal{A},\mathcal{B}\right) \\ &=\sum_{s_{1}=1}^{S} \cdots \sum_{s_{T}=1}^{S} \pi_{s_1}\mathcal{B}_{s_1,x_1}\prod_{i=2}^{T}\mathcal{A}_{s_{i-1},s_i}\mathcal{B}_{s_i,x_i} \end{aligned}$$如果按上述公式直接進行計算，複雜度為$\mathcal{O}(TS^T)$，效率太低，考慮使用動態規劃的思想，將演算法複雜度降為$\mathcal{O}(TS^2)$，可以使用兩種方式：

1. Forward Algorithm

• 定義$\alpha_{t}(j)=p\left(x_{1}, x_{2} \ldots x_{t}, s_{t}=j | \vec{\pi}, \mathcal{A},\mathcal{B}\right),\text{ }t\in\{1,2,\cdots,T\},\text{ }j\in\{1,2,\cdots,S\}$，則演算法可以寫為：
  • Initialization: $\alpha_{1}(j)=\pi_j\mathcal{B}_{j,x_1}; \quad 1 \leq j \leq S$
  • Recursion: $\alpha_{t}(j)=\sum_{i=1}^{S} \alpha_{t-1}(i) \mathcal{A}_{i j} \mathcal{B}_{j,x_t} ; \quad 1 \leq j \leq S, 1<t \leq T$
  • Termination: $p\left(\vec{o} | \vec{\pi},\mathcal{A},\mathcal{B}\right)=\sum_{i=1}^{S} \alpha_{T}(i)$

2. Backward Algorithm

• 定義$\beta_{t}(i)=p\left(x_{t+1}, x_{t+2} \ldots x_{T} | s_{t}=i, \vec{\pi}, \mathcal{A},\mathcal{B}\right),\text{ }t\in\{1,2,\cdots,T\},\text{ }i\in\{1,2,\cdots,S\}$，則演算法可以寫為：
  • Initialization: $\beta_{T}(i)=1; \quad 1 \leq i \leq S$
  • Recursion: $\beta_{t}(i)=\sum_{j=1}^{S} \mathcal{A}_{i j} \mathcal{B}_{j,x_{t+1}} \beta_{t+1}(j) ; \quad 1 \leq i \leq S, 1 \leq t < T$
  • Termination: $p\left(\vec{o} | \vec{\pi},\mathcal{A},\mathcal{B}\right)=\sum_{j=1}^{S} \pi_{j} \mathcal{B}_{j,x_1} \beta_{1}(j)$

使用EM演算法求解$\vec{\pi}_{\mathrm{ML}}, \mathcal{A}_{\mathrm{ML}}, \mathcal{B}_{\mathrm{ML}}$，關於EM演算法的具體介紹可參考文章EM演算法和高斯混合模型GMM介紹。具體到該問題，又被稱為Forward-Backward Algorithm(i.e., Baum-Welch Algorithm)：

• 假設引數的當前估計值$\mathcal{\Lambda}^*=[\vec{\pi}^*,\mathcal{A}^*,\mathcal{B}^*]$，則有
  • E-step: 定義$q(\vec{s})=p(\vec{s} | \vec{o}, \mathcal{\Lambda}^*)$，則$\mathcal{L}(\mathcal{\Lambda})=\mathbb{E}_{q}[\ln p(\vec{o}, \vec{s} |  \mathcal{\Lambda})]$。容易看出$$\ln p(\vec{o}, \vec{s} | \mathcal{\Lambda})= \underbrace{\ln \pi_{s_1}}_{\text{ initial state }} +\sum_{t=2}^{T} \underbrace{\ln \mathcal{A}_{s_{t-1}, s_t}}_{\text { Markov chain }} + \sum_{t=1}^{T} \underbrace{\ln \mathcal{B}_{s_t, x_t}}_{\text { observations }}$$因此$\mathcal{L}$的形式可寫為$$\mathcal{L}(\mathcal{\Lambda})=\sum_{k=1}^{S} \gamma_{1}(k) \ln \pi_{k}+\sum_{t=2}^{T} \sum_{j=1}^{S} \sum_{k=1}^{S} \xi_{t}(j, k) \ln \mathcal{A}_{j, k}+\sum_{t=1}^{T} \sum_{k=1}^{S} \gamma_{t}(k) \ln \mathcal{B}_{k, x_{t}}$$其中$\gamma_t(k)=p(s_t=k | \vec{o}, \mathcal{\Lambda}^*),\text{ }\xi_t(j,k)=p(s_{t-1}=j,s_t=k | \vec{o}, \mathcal{\Lambda}^*); \quad 1 \leq t \leq T$，由貝葉斯公式可知$$\begin{cases} \gamma_t(k)=\frac{p(\vec{o}, s_t=k | \mathcal{\Lambda}^*)}{p(\vec{o} | \mathcal{\Lambda}^*)}=\frac{\alpha_t(k)\beta_t(k)}{\sum_{m=1}^S\alpha_t(m)\beta_t(m)} \\ \xi_t(j,k)=\frac{p(\vec{o}, s_{t-1}=j, s_t=k | \mathcal{\Lambda}^*)}{p(\vec{o} | \mathcal{\Lambda}^*)}=\frac{\alpha_{t-1}(j)\mathcal{A}_{jk}^*\mathcal{B}_{k,x_t}^*\beta_t(k)}{\sum_{m=1}^S\alpha_t(m)\beta_t(m)} \end{cases}$$注意此時的$\alpha_t(.)$以及$\beta_t(.)$是從引數的當前估計值$\mathcal{\Lambda}^*$計算出來的
  • M-step: 更新引數的估計值$\mathcal{\Lambda}^*=\arg\max_{\mathcal{\Lambda}}\mathcal{L}(\mathcal{\Lambda})$，有$$\pi_{k}^*=\frac{\gamma_{1}(k)}{\sum_{j=1}^S \gamma_{1}(j)}, \quad \mathcal{A}_{jk}^*=\frac{\sum_{t=2}^{T} \xi_{t}(j, k)}{\sum_{t=2}^{T} \sum_{l=1}^{S} \xi_{t}(j, l)}, \quad \mathcal{B}_{kv}^*=\frac{\sum_{t=1}^{T} \gamma_{t}(k) I\left(x_{t}=v\right)}{\sum_{t=1}^{T} \gamma_{t}(k)}$$在實際應用中使用的不僅僅是1個序列，假設有$N$個序列，每個序列的長度為$T_n,\text{ }n\in\{1,2,\cdots,N\}$，則每個序列可以計算出自己的$\gamma_t,\xi_t$，記為$\gamma_t^{(n)},\xi_t^{(n)}$，更新公式變為$$\pi_{k}^*=\frac{\sum_{n=1}^{N} \gamma_{1}^{(n)}(k)}{\sum_{n=1}^{N} \sum_{j=1}^S \gamma_{1}^{(n)}(j)}, \quad \mathcal{A}_{jk}^*=\frac{\sum_{n=1}^{N} \sum_{t=2}^{T_{n}} \xi_{t}^{(n)}(j, k)}{\sum_{n=1}^{N} \sum_{t=2}^{T_{n}} \sum_{l=1}^{S} \xi_{t}^{(n)}(j, l)}, \quad \mathcal{B}_{kv}^*=\frac{\sum_{n=1}^{N} \sum_{t=1}^{T_{n}} \gamma_{t}^{(n)}(k) I\left(x_{t}^{(n)}=v\right)}{\sum_{n=1}^{N} \sum_{t=1}^{T_{n}} \gamma_{t}^{(n)}(k)}$$

• 迭代上述步驟直到收斂為止

隱藏狀態序列估計

為了解決問題2，仍使用動態規劃的思想(Viterbi Algorithm)，定義$v_{t}(j)=\max _{s_{1}, \ldots, s_{t-1}} p\left(s_{1}, \ldots s_{t-1}, x_{1}, x_{2}, \ldots x_{t}, s_{t}=j | \mathcal{\Lambda}\right)$，此外還要定義$b_t(j)$用來儲存$v_{t}(j)$對應的$s_{t-1}$，則演算法可以寫為

• Initialization: $$v_1(j)=\pi_j\mathcal{B}_{j,x_1},\text{ }b_1(j)=0; \quad 1 \leq j \leq S$$
• Recursion: $$\begin{aligned} v_{t}(j) &=\max _{i\in\{1,2,\cdots,S\}} v_{t-1}(i) \mathcal{A}_{i j} \mathcal{B}_{j,x_t} ; \quad 1 \leq j \leq S, 1<t \leq T \\ b_{t}(j) &= \underset{i\in\{1,2,\cdots,S\}}{\arg \max } v_{t-1}(i)  \mathcal{A}_{i j} \mathcal{B}_{j,x_t} ; \quad 1 \leq j \leq S, 1<t \leq T \end{aligned}$$
• Termination: $$\begin{aligned} & \text { The best score: }  P^*=\max _{\vec{s}} p\left(\vec{s},\vec{o} | \mathcal{\Lambda} \right)=\max _{i\in\{1,2,\cdots,S\}} v_{T}(i) \\ & \text { The start of backtrace: } s_{T}^*=\underset{i\in\{1,2,\cdots,S\}}{\operatorname{argmax}} v_{T}(i) \end{aligned}$$ 
• Backtrace: $$s_{t}^*=b_{t+1}(s_{t+1}^*); \quad 1 \leq t<T$$則$s_1^*,s_2^*,\ldots,s_T^*$即為$\arg\max _{\vec{s}} p\left(\vec{s},\vec{o} | \mathcal{\Lambda} \right)$，容易看出$\arg\max _{\vec{s}} p\left(\vec{s},\vec{o} | \mathcal{\Lambda} \right)=\arg\max _{\vec{s}} p\left(\vec{s} | \vec{o}, \mathcal{\Lambda} \right)$

下面介紹一個使用HMM的簡單例子：假設有兩個骰子，一個未經處理（fair, 記為0），一個經過處理（loaded, 記為1）。每次投擲時使用前一次投擲的骰子，或者使用另一個骰子，觀測序列是多次投擲所得到的骰子點數序列。假設真實的引數如下圖所示，目標是通過觀測序列估計真實的引數，進而估計出每次投擲使用的是哪個骰子。

最終的結果如下圖所示，右圖表示使用Viterbi演算法得到的最可能的隱藏狀態序列（即每次投擲使用了哪種骰子），注意左圖進行簡單地四捨五入後的結果不等於右圖，因為左圖並沒有考慮狀態之間的關聯。

&n