20201012 -
在学习《异常点检测》这本书的时候，在第十章的内容“离散数据的异常检测”记录中，涉及到隐马尔可夫模型（HMM）的学习，本篇文章具体记录HMM的学习过程。因为《异常点检测》书中关于这部分内容过于简短，本文主要学习文章[1]作为参考。
马尔可夫过程是一个随机过程，其未来状态和过去的状态有关，其中一阶的马尔可夫过程仅仅和上一时间的状态有关，如果状态空间是离散空间，该过程可以被称为马尔可夫链。
P
(
X
n
+
1
=
x
∣
X
1
=
x
1
,
X
2
,
=
x
2
,
.
.
.
,
X
n
=
x
n
)
=
P
(
X
n
+
1
=
x
∣
X
n
=
x
n
)
P(X_{n+1}=x|X_1=x_1,X_2,=x_2,...,X_n=x_n)=P(X_{n+1}=x|X_n=x_n)
P(Xn+1=x∣X1=x1,X2,=x2,...,Xn=xn)=P(Xn+1=x∣Xn=xn)
关于马尔可夫过程可以用矩阵来模拟其状态转移矩阵，该部分后续专门来研究，包括其一些平稳性的属性和其他的一些利用线性代数来解决其状态的内容。
对于一些实际的问题中，我们能够拥有一些观测值，但是却无法知晓内部的状态，也就是说存在一个潜在的变量或者说隐变量，举一个不是很恰当但是很像的例子，之前GMM模型中的内容，我们能够看到实际的观测值，但是我们不知道每个样本属于哪个高斯分布。这种包含隐藏状态的模型就是HMM。
那么对于HMM来说，在原始的马尔科夫过程的基础上，还包含了一个发射概率（emission probabilities），该概率表示基于某个内部状态，其出现某个观测值的可能性。那么总而言之，HMM将包含两个矩阵，一个是内部的状态转移矩阵，另一个是状态的观测值发射矩阵。
在HMM中，两个假设非常关键，下一个内部状态和当前的观测值仅仅依赖当前的内部状态，利用公式来描述，如下所示。
P
(
x
i
∣
x
1
,
x
2
,
.
.
.
,
x
i
−
1
)
=
P
(
x
i
∣
x
x
−
1
)
P
(
o
i
∣
x
i
,
x
2
,
.
.
.
,
x
i
−
1
)
=
P
(
o
i
∣
x
i
)
P(x_i|x_1,x_2,...,x_{i-1})=P(x_i|x_{x-1})\\ P(o_i|x_i,x_2,...,x_{i-1})=P(o_i|x_i)
P(xi∣x1,x2,...,xi−1)=P(xi∣xx−1)P(oi∣xi,x2,...,xi−1)=P(oi∣xi)
结合前面所说，正是这两个条件概率：内部状态的转移概率，外在观测值的概率；而完整的HMM就是这两个概率，再加上一个初始概率
π
\pi
π。通常来讲，HMM一般有三个应用方向：
下面的内容就专门来讲解这几个应用。
（最近一直在看该部分的理论推导，就是为什么要使用这种前向算法，一些文章中都没有具体说明，或者说非常简短，甚至于就是一句话，啊， 这个东西的计算度复杂。但这也是好事，就是我通过自己的学习把全概率公式，乘法公式，以及乘法公式都复习了一遍，虽然这东西都是基础，但是重新要自己的推导的时候，还是不是那么简单的。）
那么我按照我最近学习到的东西，仔细来捋一捋具体的内容。
首先再次明确HMM的两个性质，马尔可夫假设与观测值输出独立性。
I
n
t
e
r
n
a
l
S
t
a
t
e
:
P
(
q
i
∣
q
1
.
.
.
q
i
−
1
)
=
P
(
q
i
∣
q
i
−
1
)
O
u
t
p
u
t
I
n
d
e
p
e
n
d
e
n
c
e
:
P
(
o
i
∣
q
1
,
.
.
.
,
q
i
,
.
.
.
,
q
T
,
o
1
,
.
.
.
,
o
i
,
.
.
.
,
o
T
)
=
P
(
o
i
∣
q
i
)
Internal \quad State:P(q_i|q_1...q_{i-1})=P(q_i|q_{i-1})\\ Output \quad Independence: P(o_i|q_1,...,q_i,...,q_T,o_1,...,o_i,...,o_T)=P(o_i|q_i)
InternalState:P(qi∣q1...qi−1)=P(qi∣qi−1)OutputIndependence:P(oi∣q1,...,qi,...,qT,o1,...,oi,...,oT)=P(oi∣qi)
来看一下某个事件发生的概率，其中
Y
Y
Y代表某个序列值，
X
X
X代表内部状态。
P
(
Y
)
=
∑
x
P
(
Y
,
X
)
=
∑
x
P
(
Y
∣
X
)
P
(
X
)
P(Y)=\sum_x P(Y,X)=\sum_x P(Y|X)P(X)
P(Y)=x∑P(Y,X)=x∑P(Y∣X)P(X)
上述式子中，右边第一个是边缘分布律，通过遍历所有的内部状态空间
X
X
X，这个就是全概率公式，通过遍历整个状态空间来进行计算，注意，这里的全部状态必须能够填充事件
X
X
X的全部空间，而最后一个式子是乘法公式。那么如果要计算某个序列值的概率，可以通过遍历当前隐藏状态的所有状态来进行计算。但是这里还是存在这整个序列，依然是不好计算的。
P
(
y
1
,
y
2
,
.
.
.
,
y
t
∣
λ
)
=
∑
j
P
(
y
1
,
y
2
,
.
.
.
,
y
t
,
x
t
=
s
j
∣
λ
)
P(y_1,y_2,...,y_t|\lambda)=\sum_j P(y_1,y_2,...,y_t,x_t=s_j|\lambda)
P(y1,y2,...,yt∣λ)=j∑P(y1,y2,...,yt,xt=sj∣λ)
那么假设，我知道所有的观测序列
O
O
O和状态序列
Q
Q
Q，那么他们的状态值是多少呢？看下下面的公式。
P
(
O
,
Q
)
=
P
(
O
∣
Q
)
×
P
(
Q
)
=
∏
i
=
1
T
P
(
o
i
∣
q
i
)
×
∏
i
=
1
T
P
(
q
i
∣
q
i
−
1
)
P(O,Q)=P(O|Q)\times P(Q)=\prod_{i=1}^{T}P(o_i|q_i)\times\prod_{i=1}^{T}P(q_i|q_{i-1})
P(O,Q)=P(O∣Q)×P(Q)=i=1∏TP(oi∣qi)×i=1∏TP(qi∣qi−1)
上述公式就是简单利用了乘法公式和隐马尔可夫模型的特性来生成上述概率公式，因为这里他是一个连续的独立事件，可以通过乘法公式来进行展开，然后叠乘起来。但实际情况是，在估计观测值的时候，并知道具体的内部状态，仅仅有观测值。那么依然使用全概率公式。
P
(
O
)
=
∑
Q
P
(
O
,
Q
)
=
∑
Q
P
(
O
∣
Q
)
P
(
Q
)
P(O)=\sum_{Q} P(O,Q)=\sum_{Q}P(O|Q)P(Q)
P(O)=Q∑P(O,Q)=Q∑P(O∣Q)P(Q)
但是这种计算方式的话，需要计算全部的概率，也就是需要把所有的内部状态进行全排列，然后逐个利用前面的公式来进行计算，那么计算的隐藏序列个数将达到
N
T
N^T
NT个，这个量级是无法接受的。既然这样的话，那么可以利用前向算法或者后向算法的方式采用动态编程的思想来进行计算。
计算的方式就是，通过遍历前一个时间点的所有隐藏状态，然后利用叠加的方式计算计算出来现在每个隐藏状态的概率值，然后乘上此时的发射概率，来看一个比较形象的图片，该图片来自文章[2]
那么从公式的角度应该怎么写呢？
P
(
O
t
)
=
∑
i
P
(
O
t
,
Q
t
=
s
i
)
=
∑
i
P
(
O
t
∣
Q
t
=
s
i
)
P
(
Q
t
=
s
i
)
=
∑
i
P
(
O
t
∣
Q
t
=
s
i
)
∑
j
P
(
Q
t
=
s
i
,
Q
t
−
1
=
s
j
)
=
∑
i
P
(
O
t
∣
Q
t
=
s
i
)
∑
j
P
(
Q
t
=
s
i
∣
Q
t
−
1
=
s
j
)
P
(
Q
t
−
1
=
s
j
)
P(O_t)=\sum_iP(O_t, Q_t=s_i)=\sum_i P(O_t|Q_t=s_i)P(Q_t=s_i)\\ =\sum_iP(O_t|Q_t=s_i)\sum_j P(Q_t=s_i,Q_{t-1}=s_j)\\ =\sum_iP(O_t|Q_t=s_i)\sum_j P(Q_t=s_i|Q_{t-1}=s_j)P(Q_{t-1}=s_j)
P(Ot)=i∑P(Ot,Qt=si)=i∑P(Ot∣Qt=si)P(Qt=si)=i∑P(Ot∣Qt=si)j∑P(Qt=si,Qt−1=sj)=i∑P(Ot∣Qt=si)j∑P(Qt=si∣Qt−1=sj)P(Qt−1=sj)
上述公式是没错的，就是最后落实到了这个隐藏状态的概率上，但是怎么用递归的方式来表达出来呢？
（很多文章这里都是直接使用了递归公式，直接就是令这个东西得什么，虽然能知道这是正确的，但是没办法理解这个思考的过程）
那么既然这里我希望通过使用递归的方式，甚至于后面能够使用动态规划的思想来实现，那么必须能够写出来递归方程才行，要求解的内容是
P
(
O
1
:
t
)
=
∑
i
P
(
O
1
:
t
,
Q
t
=
s
i
)
P(O_{1:t})=\sum_i P(O_{1:t},Q_t=s_i)
P(O1:t)=∑iP(O1:t,Qt=si)
这里出现了这个公式，那么我能不能推导出来
P
(
O
1
:
t
−
1
)
P(O_{1:t-1})
P(O1:t−1)或者
P
(
O
1
:
t
−
1
,
Q
t
−
1
=
s
i
)
P(O_{1:t-1},Q_{t-1}=s_i)
P(O1:t−1,Qt−1=si)，也就是说他们是相同的方式呢？他们在整体的表现形式上是一样的。那么利用前向算法的思想。
P
(
O
1
:
t
,
Q
t
=
s
i
)
=
X
X
X
P
(
O
1
:
t
−
1
,
Q
t
−
1
=
s
j
)
P(O_{1:t},Q_t=s_i)=XXXP(O_{1:t-1},Q_{t-1}=s_j)
P(O1:t,Qt=si)=XXXP(O1:t−1,Qt−1=sj)
这里的
X
X
X
XXX
XXX是一个未知的内容，那么应该怎么写出来，下面来尝试推导一下。
P
(
O
t
,
Q
t
=
s
i
)
=
P
(
O
1
:
t
−
1
,
O
t
,
Q
t
=
s
i
)
=
P
(
O
t
∣
O
1
:
t
−
1
,
Q
t
=
s
i
)
P
(
O
1
:
t
−
1
,
Q
t
=
s
i
)
=
P
(
O
t
∣
Q
t
=
s
i
)
P
(
O
1
:
t
−
1
,
Q
t
=
s
i
)
=
P
(
O
t
∣
Q
t
=
s
i
)
∑
j
P
(
O
1
:
t
−
1
,
Q
t
=
s
i
,
Q
t
−
1
=
s
j
)
=
P
(
O
t
∣
Q
t
=
s
i
)
∑
j
P
(
O
1
:
t
−
1
,
Q
t
−
1
=
s
j
)
P
(
Q
t
=
s
i
∣
O
1
:
t
−
1
,
Q
t
−
1
=
s
j
)
=
P
(
O
t
∣
Q
t
=
s
i
)
∑
j
P
(
O
1
:
t
−
1
,
Q
t
−
1
=
s
j
)
P
(
Q
t
=
s
i
∣
Q
t
−
1
=
s
j
)
P(O_t,Q_t=s_i)=P(O_{1:t-1},O_t,Q_t=s_i)=\\ P(O_t|O_{1:t-1},Q_t=s_i)P(O_{1:t-1},Q_t=s_i)=\\ P(O_t|Q_t=s_i)P(O_{1:t-1},Q_t=s_i)=\\ P(O_t|Q_t=s_i)\sum_j P(O_{1:t-1},Q_t=s_i,Q_{t-1}=s_j)=\\ P(O_t|Q_t=s_i)\sum_j P(O_{1:t-1},Q_{t-1}=s_j)P(Q_t=s_i|O_{1:t-1},Q_{t-1}=s_j)=\\ P(O_t|Q_t=s_i)\sum_j P(O_{1:t-1},Q_{t-1}=s_j)P(Q_t=s_i|Q_{t-1}=s_j)
P(Ot,Qt=si)=P(O1:t−1,Ot,Qt=si)=P(Ot∣O1:t−1,Qt=si)P(O1:t−1,Qt=si)=P(Ot∣Qt=si)P(O1:t−1,Qt=si)=P(Ot∣Qt=si)j∑P(O1:t−1,Qt=si,Qt−1=sj)=P(Ot∣Qt=si)j∑P(O1:t−1,Qt−1=sj)P(Qt=si∣O1:t−1,Qt−1=sj)=P(Ot∣Qt=si)j∑P(O1:t−1,Qt−1=sj)P(Qt=si∣Qt−1=sj)
上述推导过程中，主要使用了全概率公式和乘法公式，同时结合HMM的性质得到。
令
a
i
j
a_{ij}
aij表示状态转移矩阵，而
b
j
b_j
bj表示观测值的发射矩阵，那么改写上述公式。
α
t
(
i
)
=
P
(
O
1
:
t
,
Q
t
=
s
i
)
α
t
(
i
)
=
b
i
(
Q
t
)
∑
j
a
j
i
α
t
−
1
(
j
)
\alpha_t(i)=P(O_{1:t},Q_t=s_i)\\ \alpha_t(i)=b_i(Q_t)\sum_j a_{ji}\alpha_{t-1}(j)
αt(i)=P(O1:t,Qt=si)αt(i)=bi(Qt)j∑ajiαt−1(j)
那么整体的计算过程如下：
I
n
i
t
i
a
l
i
z
e
:
α
1
(
j
)
=
π
j
b
j
(
Q
1
)
E
a
c
h
S
t
e
p
:
α
t
(
j
)
=
∑
i
=
1
N
α
t
−
1
(
i
)
a
i
j
b
j
(
y
t
)
F
i
n
a
l
R
e
s
u
l
t
:
P
(
Y
a
l
l
∣
λ
)
=
∑
i
=
1
N
α
T
(
i
)
Initialize: \qquad\alpha_1(j)=\pi_jb_j(Q_1)\\ Each Step:\qquad \alpha_t(j)=\sum_{i=1}^{N}\alpha_{t-1}(i)a_{ij}b_{j}(y_t)\\ Final Result:\qquad P(Y_{all}|\lambda)=\sum_{i=1}^{N}\alpha_T(i)
Initialize:α1(j)=πjbj(Q1)EachStep:αt(j)=i=1∑Nαt−1(i)aijbj(yt)FinalResult:P(Yall∣λ)=i=1∑NαT(i)
通过动态规划的方式，父问题利用子问题的解，然后得到最终的结果。
在很多文章中，都是直接利用了相应的结论，直接写出了递归的方程，这种方式并不是很好理解，一般通过自己推导出来的过程，记忆会更深刻。
从公式上来看，解码过程的公式跟求解似然概率的过程很像。仅仅是一个累加符号变成了max符号。但是我觉得，主要还是要思考这个过程中各个点是怎么样的关系来串联起来的。就是一个关联思考的过程，前面的地方是一种利用边缘概率或者说全概率公式将所有的情况都涵盖起来，然后计算出来，这里呢？目前来看我比较关注的一个点就是，这种仅仅max的方式，是不是有点像贪心算法？！每次仅仅选择最大的。
这部分首先参考了文章[3]，在一开始思考这部分问题的时候，就有点想不通为什么要用max这个操作呢？其实我觉得，他文章中的一段话阐述的非常好。我这里直接引用过来。
Like Baum-Welch algorithm, our training algorithm, the Viterbi algorithm is also a dynamic programming approach. As a reminder, the dynamic programming is an approach that you would reuse calculated result in the next calculation, and the act of reusing saves time! What’s more, it is indeed a natural choice for HMM, take the simple HMM as an example, in which any state depends on its first previous state. This statement, however, does not mean that there is no effect at all from the rest of the previous states to the current one, but their effects are all absorbed into the first previous state, and this first previous state becomes the only factor you need to consider when transiting to the next state.
来说一下其所传达的意思，虽然在公式中，某个内部状态仅仅和其前一内部状态有关，但是因为这里是要对整个序列进行预测，那么这就是一个联合概率，联合概率必然是要利用乘法的方式，你把整个所有的内部状态都进行一下内部状态的叠乘的过程。而max的过程，就是在前面的基础上，每次都选最大的，最后才会使整个内部状态序列最大的，其实，如果不加这个max的内容，假设你就是想得到某个内部状态序列的概率值，也是通过乘法公式的方式来实现。现在为了更好的理解。现在就需要把这部分的推导公式找出来，正如之前讲解的前向算法一样，他的完整概率推导公式是什么样子。
（我个人觉得，这部分内容，还是需要好好的思考，能够把这个公式计算的过程完全理解，而且我觉得，阅读这个东西的过程中，越来越觉得其实本质上需要思考这个动态规划的具体细节内容）
[1]Machine Learning — Hidden Markov Model (HMM)
[2]隐马尔可夫模型（三）——隐马尔可夫模型的评估问题(前向算法）
[3]Viterbi algorithm for prediction with HMM — Part 3 of the HMM series
