隐马尔可夫链简介
定义 1 [隐马尔可夫链 (HMM)]:
一个双离散随机过程 (double discrete stochastic process) { C t } , { X t } 被称作隐马尔可夫链,若它满足以下两个条件:
⎩ ⎨ ⎧ ( 1 ) P ( C t = c t ∣ C 1 t − 1 = c 1 t − 1 ) = P ( C t = c t ∣ C t − 1 = c t − 1 ) ( 2 ) P ( X t = x t ∣ C 1 t = c 1 t , X 1 t − 1 = x 1 t − 1 ) = P ( X t = x t ∣ C t = c t ) (HMM 性质 1) (HMM 性质 2)
使用简便记号,隐马尔可夫链的定义也可简写为
{ p ( c t ∣ c 1 t − 1 ) = p ( c t ∣ c t − 1 ) p ( x t ∣ c 1 t , x 1 t − 1 ) = p ( x t ∣ c t )
在隐马尔可夫链的定义中,{ C t } 被称作状态过程 (state process),它是不可观察的(隐藏的),HMM 性质 1 说明其是一个马尔可夫链;{ X t } 被称为观测值过程 (observation process),HMM 性质 2 说明它仅与当前状态有关,而与之前的状态和观测值均无关。
在此基础上,我们还可以定义加强版(改进版)的隐马尔可夫链,如下所示。
定义 2 [加强版的隐马尔可夫链 (Revised HMM)]:
( 1 ) p ( c t ∣ x 1 t − 1 , c 1 t − 1 ) = p ( c t ∣ c t − 1 ) ( 加强版的 HMM 性质 1)
( 2 ) p ( x t ∣ x 1 , … , x t − 1 , x t + 1 , … , x T , c 1 , … , c t ) = p ( x t ∣ c t ) ( 加强版的 HMM 性质 2)
加强版的 HMM 性质 1 说明,给定第 t − 1 个状态变量,第 t 个状态变量与其他所有变量均不相关;加强版的 HMM 性质 2 说明,给定第 t 个状态变量,第 t 个观测值与其他所有变量均不相关。
在后续部分中,我们将不再详细区分 HMM 性质与加强版的 HMM 性质。总的来说,在隐马尔可夫链框架下,一共两个随机过程,{ C t } 与 { X t } ;{ C t } 是马尔可夫链,给定 C t − 1 ,C t 与其他所有变量(过去的观测值和过去的状态)均无关,{ X t } 是观测值过程,给定 C t ,X t 与其他所有变量(之前和之后的观测值、过去的状态)均无关。
三个基本问题
应用隐马尔可夫链面临三个基本问题 :
当参数给定时,求出观测值序列发生的概率,即
P ( X 1 T = x 1 T ∣ θ )
通常对于这个问题我们使用向前/向后算法 (Forward/Backward Approach)
找出一个隐藏的状态序列 c 1 , c 2 , … , c T ,其能最好的解释观测值序列 x 1 , x 2 , … , x T
这个问题也称为解码 (decoding),通常使用维特比算法 (Viterbi Algorithm)
找出能使观测值序列发生概率最大的参数,即
θ ^ = arg θ max P ( X 1 T = x 1 T )
此问题即参数估计 ,也称为学习 (study) 或者训练 (train),是三个基本问题中最复杂的,通常使用期望最大化算法 (Expectation Maximization Method or called Baum–Welch Method)。
隐马尔可夫链的若干基本性质
我们首先定义一个符号
p i ( x ) = P ( X t = x ∣ C t = i )
可见 p i ( x ) 是给定 C t 在 t 时刻位于状态 i 时,X t 的条件概率。它也被称为状态相关分布 (state-dependent distribution),或者发射分布 (emission distribution)。形象地理解,可以把观测值 X t 当作状态 C t “发射”的“信号”(signal)。
在此基础上,定义
P ( x ) = p 1 ( x ) p 2 ( x ) ⋱ p m ( x ) m × m
P ( x ) 被称为对角发射矩阵 (diagonal emission matrix),它的第 i 个对角线元素即是 p i ( x ) 。
定理 1 [单变量分布 (univariate distribution)]:
P ( X t = x ) = i = 1 ∑ m u i ( t ) ⋅ p i ( x )
写成矩阵形式为
P ( X t = x ) = u ( t ) ⋅ P ( x ) ⋅ 1 T
证明:
P ( X t = x ) = i = 1 ∑ m P ( X t = x ∣ C t = i ) ⋅ P ( C t = i ) (基于定理 1.1 ,全概率公式)
= i = 1 ∑ m u i ( t ) ⋅ p i ( x )
上式也可以用矩阵运算写成
P ( X t = x ) = [ u 1 ( t ) , u 2 ( t ) , … , u m ( t ) ] 1 × m ⋅ p 1 ( x ) p 2 ( x ) ⋱ p m ( x ) m × m ⋅ 1 1 ⋮ 1 m × 1
= u ( t ) ⋅ P ( x ) ⋅ 1 T
此外,我们知道,u ( t ) = u ( 1 ) ⋅ A t − 1 ,因此,上述定理的矩阵形式也可写为
P ( X t = x ) = u ( 1 ) ⋅ A t − 1 ⋅ P ( x ) ⋅ 1 T = π ⋅ A t − 1 ⋅ P ( x ) ⋅ 1 T
若马尔可夫链是稳态的,u ( t ) ≡ δ ,则上式又可简化为
P ( X t = x ) = δ ⋅ P ( x ) ⋅ 1 T
定理 2:
p ( x t , x t + k , c t , c t + k ) = p ( c t ) ⋅ p c t ( x t ) ⋅ a c t , c t + k ( k ) ⋅ p c t + k ( x t + k )
证明:
p ( x t , x t + k , c t , c t + k ) = p ( x t , x t + k ∣ c t , c t + k ) ⋅ p ( c t , c t + k ) = p ( x t ∣ c t , c t + k ) ⋅ p ( x t + k ∣ x t , c t , c t + k ) ⋅ p ( c t + k ∣ c t ) ⋅ p ( c t ) = p ( x t ∣ c t ) ⋅ p ( x t + k ∣ c t + k ) ⋅ p ( c t ) ⋅ a c t , c t + k ( k ) = p ( c t ) ⋅ p c t ( x t ) ⋅ a c t , c t + k ( k ) ⋅ p c t + k ( x t + k ) (取条件于 c t , c t + k ) (基于引理 1 及取条件于 c t ) (基于 HMM 性质 2 )
证毕
定理 2 可以直观地理解如下:隐马尔可夫链在时刻 t 和 t + k ,状态变量取值为 c t 和 c t + k ,观测变量取值为 x t 和 x t + k 的概率,等于在 t 时刻,状态变量取值于 c t 的概率,乘以“发射信号”x t 的概率,乘以状态变量 k 步转移到 c t + k 的概率(此时时刻为 t + k ),再乘以“发射信号”x t + k 的概率。
定理 3 [双变量分布 (Bivariate Distribution)]:
P ( X t = v , X t + k = w ) = i = 1 ∑ m j = 1 ∑ m u i ( t ) ⋅ p i ( v ) ⋅ a ij ( k ) ⋅ p j ( w )
矩阵形式可写为
P ( X t = v , X t + k = w ) = u ( t ) ⋅ P ( v ) ⋅ A k ⋅ P ( w ) ⋅ 1 T
证明:
P ( X t = v , X t + k = w ) = P ( X t = v , X t + k = w , i ⋃ ( C t = i ) , j ⋃ ( C t + k = j ) ) = i = 1 ∑ m j = 1 ∑ m P ( X t = v , X t + k = w , C t = i , C t + k = j ) = i = 1 ∑ m j = 1 ∑ m P ( C t = i ) ⋅ p i ( v ) ⋅ a ij ( k ) ⋅ p j ( w ) = i = 1 ∑ m j = 1 ∑ m u i ( t ) ⋅ p i ( v ) ⋅ a ij ( k ) ⋅ p j ( w ) (基于定理 2 )
上式也可用矩阵运算写为
P ( X t = v , X t + k = w ) = [ u 1 ( t ) , … , u m ( t ) ] 1 × m ⋅ p 1 ( v ) p 2 ( v ) ⋱ p m ( v ) m × m ⋅ a 11 ( k ) ⋮ a m 1 ( k ) ⋯ ⋱ ⋯ a 1 m ( k ) ⋮ a mm ( k ) m × m ⋅ p 1 ( w ) p 2 ( w ) ⋱ p m ( w ) m × m ⋅ 1 1 ⋮ 1 m × 1
= u ( t ) ⋅ P ( v ) ⋅ A ( k ) ⋅ P ( w ) ⋅ 1 T
= u ( t ) ⋅ P ( v ) ⋅ A k ⋅ P ( w ) ⋅ 1 T (基于定理 2 A ( k ) = A k )
证毕
若马尔可夫链是稳态的,则 u ( t ) ≡ δ ,上式又可简化为
P ( X t = v , X t + k = w ) = δ ⋅ P ( v ) ⋅ A k ⋅ P ( w ) ⋅ 1 T
定理 4:
P ( X 1 T = x 1 T , C 1 T = c 1 T ) = π c 1 ⋅ t = 1 ∏ T − 1 a c t , c t + 1 ⋅ t = 1 ∏ T p c t ( x t ) = π c 1 ⋅ t = 2 ∏ T a c t − 1 , c t ⋅ t = 1 ∏ T p c t ( x t )
证明:
P ( X 1 T = x 1 T , C 1 T = c 1 T ) = p ( x 1 T , c 1 T ) = p ( x 1 T ∣ c 1 T ) ⋅ p ( c 1 T ) (取条件于 c 1 T )
p ( c 1 T ) = p ( c 2 T ∣ c 1 ) ⋅ p ( c 1 ) = π c 1 ⋅ p ( c 2 ∣ c 1 ) ⋅ p ( c 3 T ∣ c 1 , c 2 ) = π c 1 ⋅ a c 1 , c 2 ⋅ p ( c 3 T ∣ c 2 ) = π c 1 ⋅ a c 1 , c 2 ⋅ p ( c 3 ∣ c 2 ) ⋅ p ( c 4 T ∣ c 2 , c 3 ) ⋮ = π c 1 ⋅ a c 1 , c 2 ⋅ a c 2 , c 3 ⋯ a c T − 1 , c T = π c 1 ⋅ t = 1 ∏ T − 1 a c t , c t + 1 = π c 1 ⋅ t = 2 ∏ T a c t − 1 , c t (取条件于 c 1 ) (基于引理 1 ) (基于 HMM 性质 1 ) (基于引理 1 ) (重复上述过程)
p ( x 1 T ∣ c 1 T ) = p ( x 1 ∣ c 1 T ) ⋅ p ( x 2 T ∣ c 1 T , x 1 ) = p ( x 1 ∣ c 1 ) ⋅ p ( x 2 T ∣ c 2 T ) = p c 1 ( x 1 ) ⋅ p ( x 2 ∣ c 2 T ) ⋅ p ( x 3 T ∣ c 2 T , x 2 ) = p c 1 ( x 1 ) ⋅ p ( x 2 ∣ c 2 ) ⋅ p ( x 3 T ∣ c 3 T ) = p c 1 ( x 1 ) ⋅ p c 2 ( x 2 ) ⋅ p ( x 3 ∣ c 3 T ) ⋅ p ( x 4 T ∣ c 3 T , x 3 ) ⋮ = p c 1 ( x 1 ) ⋅ p c 2 ( x 2 ) ⋯ p c T ( x T ) = t = 1 ∏ T p c t ( x t ) (基于引理 1 ) (基于 HMM 性质 2 ) (基于引理 1 ) (基于 HMM 性质 2 ) (基于引理 1 ) (重复上述过程)
结合上述两步的结果,可得
P ( X 1 T = x 1 T , C 1 T = c 1 T ) = π c 1 ⋅ t = 1 ∏ T − 1 a c t , c t + 1 ⋅ t = 1 ∏ T p c t ( x t ) = π c 1 ⋅ t = 2 ∏ T a c t − 1 , c t ⋅ t = 1 ∏ T p c t ( x t )
证毕
定理 4 可以直观地理解如下:出现状态序列 c 1 , c 2 , … , c T 和观测值序列 x 1 , x 2 , … , x T 的概率,等于初始时刻状态位于 c 1 的概率(右边第 1 项),乘以状态从 c 1 转移到 c 2 的概率,再乘以状态从 c 2 转移到 c 3 的概率,以此类推,直到乘以状态从 c T − 1 转移到 c T 的概率(右边第 2 项),最后还要乘以在每个状态 c t 下“发射信号”x t 的概率(右边第 3 项)。
从定理 4 出发,我们也可以导出观测值序列 x 1 , x 2 , … , x T 出现的概率。
定理 5:
P ( X 1 T = x 1 T ) = c 1 , c 2 , … , c T = 1 ∑ m P ( X 1 T = x 1 T , C 1 T = c 1 T )
证明:
P ( X 1 T = x 1 T ) = P ( X 1 T = x 1 T , c 1 = 1 ⋃ m ( C 1 = c 1 ) ) = c 1 = 1 ∑ m P ( X 1 T = x 1 T , C 1 = c 1 ) = c 1 = 1 ∑ m P ( X 1 T = x 1 T , C 1 = c 1 , c 2 ⋃ ( C 2 = c 2 ) ) = c 1 = 1 ∑ m c 2 = 1 ∑ m P ( X 1 T = x 1 T , C 1 = c 1 , C 2 = c 2 ) ⋮ = c 1 = 1 ∑ m c 2 = 1 ∑ m ⋯ c T = 1 ∑ m P ( X 1 T = x 1 T , C 1 = c 1 , C 2 = c 2 , … , C T = c T ) = c 1 , c 2 , … , c T = 1 ∑ m P ( X 1 T = x 1 T , C 1 T = c 1 T ) (重复上述步骤)
证毕
隐马尔可夫链的似然函数
隐马尔可夫链的似然函数是个很重要的概念,在参数估计、条件分布、预测等问题中有着非常重要的应用。此外,后面介绍的前向(后向)概率等许多参量也可以用似然函数来表示。幸运的是,不但正常情况下隐马尔可夫链的似然函数有明晰的解析表达式,在个别数据缺失的情况下,依然可以得到似然函数的解析表达式。
正常情况下隐马尔可夫链的似然函数
定理 6 [无缺失数据时 HMM 的似然函数]:
若隐马尔可夫链共有 T 个观测值,x 1 , x 2 , … , x T ,则其似然函数可表示为
L = π P ( x 1 ) A P ( x 2 ) ⋯ A P ( x T ) ⋅ 1 T = π P ( x 1 ) ⋅ s = 2 ∏ T A P ( x s ) ⋅ 1 T
证明:
L = P ( X 1 T = x 1 T ) = c 1 , … , c T = 1 ∑ m P ( X 1 T = x 1 T , C 1 T = c 1 T ) = c 1 , … , c T = 1 ∑ m π c 1 ⋅ t = 2 ∏ T a c t − 1 , c t ⋅ t = 1 ∏ T p c t ( x t ) = c 1 , … , c T = 1 ∑ m π c 1 ⋅ a c 1 , c 2 ⋅ a c 2 , c 3 ⋯ a c T − 1 , c T ⋅ p c 1 ( x 1 ) ⋅ p c 2 ( x 2 ) ⋅ p c 3 ( x 3 ) ⋯ p c T ( x T ) = π P ( x 1 ) A P ( x 2 ) A P ( x 3 ) ⋯ A P ( x T ) ⋅ 1 T = π P ( x 1 ) ⋅ s = 2 ∏ T A P ( x s ) ⋅ 1 T (基于定理 5 ) (基于定理 4 )
从矩阵运算的观点来看,LHS(似然函数)是一个实数,RHS 是 ( 1 × m ) ⋅ ( m × m ) ⋯ ( m × m ) ⋅ ( m × 1 ) 的矩阵乘法,其结果 ( 1 × 1 ) 也是一个实数。 证毕
特别地,如果马尔可夫链(隐马尔可夫链的状态过程)是稳态的,则有 π = δ = δ ⋅ A ,因此上式可简化为
L = δ ⋅ A P ( x 1 ) ⋅ s = 2 ∏ T A P ( x s ) ⋅ 1 T = δ ⋅ s = 1 ∏ T A P ( x s ) ⋅ 1 T
数据缺失情况下隐马尔可夫链的似然函数
在个别数据缺失的情况下,隐马尔可夫链的似然函数也有类似的表达式。
定理 7 [有缺失数据时 HMM 的似然函数]:
若观察值 x 3 , x 5 , x 6 缺失,则似然函数可表示为
L − ( 3 , 5 , 6 ) = π P ( x 1 ) ⋅ A P ( x 2 ) ⋅ A 2 P ( x 4 ) ⋅ A 3 P ( x 7 ) ⋯ A P ( x T ) ⋅ 1 T
其中 L − ( 3 , 5 , 6 ) 代表 x 3 , x 5 , x 6 缺失时的似然函数。
证明:
L − ( 3 , 5 , 6 ) = P ( X 1 = x 1 , X 2 = x 2 , X 4 = x 4 , X 7 = x 7 , X 8 T = x 8 T ) = 除 c 3 , c 5 , c 6 之外的所有 c t ∑ π c 1 ⋅ a c 1 , c 2 ⋅ a c 2 , c 4 ( 2 ) ⋅ a c 4 , c 7 ( 3 ) ⋅ a c 7 , c 8 ⋯ a c T − 1 , c T ⋅ p c 1 ( x 1 ) ⋅ p c 2 ( x 2 ) ⋅ p c 4 ( x 4 ) ⋅ p c 7 ( x 7 ) ⋅ p c 8 ( x 8 ) ⋯ p c T ( x T ) = 除 c 3 , c 5 , c 6 之外的所有 c t ∑ π c 1 ⋅ p c 1 ( x 1 ) ⋅ a c 1 , c 2 ⋅ p c 2 ( x 2 ) ⋅ a c 2 , c 4 ( 2 ) ⋅ p c 4 ( x 4 ) ⋅ a c 4 , c 7 ( 3 ) ⋅ p c 7 ( x 7 ) ⋅ a c 7 , c 8 ⋅ p c 8 ( x 8 ) ⋯ a c T − 1 , c T ⋅ p c T ( x T ) = π P ( x 1 ) ⋅ A P ( x 2 ) ⋅ A ( 2 ) P ( x 4 ) ⋅ A ( 3 ) P ( x 7 ) ⋅ A P ( x 8 ) ⋯ A P ( x T ) ⋅ 1 T = π P ( x 1 ) ⋅ A P ( x 2 ) ⋅ A 2 P ( x 4 ) ⋅ A 3 P ( x 7 ) ⋅ A P ( x 8 ) ⋯ A P ( x T ) ⋅ 1 T (基于定理 2 A ( k ) = A k )
证毕
评注: 对于缺失数据的似然函数,有一个简便的记忆方法。若 x t 缺失,则对应的对角发射矩阵(diagonal emission matrix),P ( x t ) ,就被替换为单位阵 I ;等价的,对于所有的 i = 1 , 2 , … , m ,p i ( x t ) ≡ 1 。
例: x 3 缺失
则 P ( x 3 ) → I ,似然函数的变化为
⋯ P ( x 2 ) A P ( x 3 ) A P ( x 4 ) ⋯ → ⋯ P ( x 2 ) ⋅ A ⋅ I ⋅ A P ( x 4 ) ⋯ = ⋯ P ( x 2 ) A 2 P ( x 4 ) ⋯
例: x 5 , x 6 缺失
则 P ( x 5 ) → I , P ( x 6 ) → I ,似然函数变化为
⋯ P ( x 4 ) A P ( x 5 ) A P ( x 6 ) A P ( x 7 ) ⋯ → ⋯ P ( x 4 ) ⋅ A ⋅ I ⋅ A ⋅ I ⋅ A ⋅ P ( x 7 ) = ⋯ P ( x 4 ) A 3 P ( x 7 )
数据缺失情况下的隐马尔可夫链似然函数在后文提到的分布预测等问题中有很多应用。
两类隐马尔可夫链
我们将介绍两类常见的隐马尔可夫链,泊松-隐马尔可夫链 (Poisson-HMM) 和正态-隐马尔可夫链 (Normal-HMM)。它们都是从“发射信号”这个角度来划分的。若发射概率是泊松分布的,则称这种 HMM 为泊松-隐马尔可夫链,正态-隐马尔可夫链与此同理。类似地,也可定义二项-隐马尔可夫链 (Binomial-HMM) 或者伽马-隐马尔可夫链 (Gamma-HMM) 等。
定义 3 [泊松-隐马尔可夫链 (Poisson-HMM)]:
若隐马尔可夫链共有 m 个状态,T 个观测值,即 i = 1 , 2 , … , m , t = 1 , 2 , … , T ,在每个状态 i 下,发射概率(emission probability)是泊松分布的并且强度为 λ i ,即
p i ( x t ) = P ( X t = x t ∣ C t = i ) = x t ! e − λ i ⋅ λ i x t
则称此 HMM 为泊松-隐马尔可夫链。
在泊松-隐马尔可夫链框架下,参数为 θ = [ π , A , λ ] ,其中
⎩ ⎨ ⎧ π = [ π 1 , π 2 , … , π m ] 1 × m A = a 11 ⋮ a m 1 ⋯ ⋱ ⋯ a 1 m ⋮ a mm m × m λ = [ λ 1 , λ 2 , … , λ m ] 1 × m (初始分布向量) (转移概率矩阵) (发射参数向量)
与泊松-隐马尔可夫链类似,我们也可以定义正态-隐马尔可夫链,区别仅在于在每个状态 i 下,发射概率为正态分布而非泊松分布。
定义 4 [正态-隐马尔可夫链 (Normal-HMM)]:
若隐马尔可夫链共有 m 个状态,T 个观测值,即 i = 1 , 2 , … , m , t = 1 , 2 , … , T ,在每个状态 i 下,发射概率是正态分布的并且参数为 μ i 和 σ i ,即
p i ( x t ) = P ( X t = x t ∣ C t = i ) = 2 π σ i 1 e − ( x t − μ i ) 2 / ( 2 σ i 2 )
则称此 HMM 为正态-隐马尔可夫链。
在正态-隐马尔可夫链框架下,参数为 θ = [ π , A , μ , σ ] ,其中
⎩ ⎨ ⎧ π = [ π 1 , π 2 , … , π m ] 1 × m A = a 11 ⋮ a m 1 ⋯ ⋱ ⋯ a 1 m ⋮ a mm m × m μ = [ μ 1 , μ 2 , … , μ m ] 1 × m σ = [ σ 1 , σ 2 , … , σ m ] 1 × m (初始分布向量) (转移概率矩阵) (发射参数向量) (发射参数向量)
向前/向后算法
我们曾提到过,隐马尔可夫链的核心是三个基本问题,其中第一个问题就是计算观测值出现的概率,P ( X 1 T = x 1 T ∣ θ ) 。回答这个问题需使用向前/向后算法,而这个算法的核心是前向/后向概率。此外,前向/后向概率也是其他重要算法如后面提到的期望最大化算法的基础。
前向概率与向前算法
定义 5 [前向概率 (Forward Probability)]:
α i ( t ) = P ( X 1 t = x 1 t , C t = i ) ( t = 1 , 2 , … , T ; i = 1 , 2 , … , m )
可见前向概率是一个联合概率,它是出现观测值 x 1 , x 2 , … , x t 且状态在 t 时刻位于 i 的概率。
仿照前面行向量 u ( t ) , π 的定义,我们也可以定义前向概率的向量表示 (vector representation),
α ( t ) = [ α 1 ( t ) , α 2 ( t ) , … , α m ( t ) ] 1 × m
对于前向概率的向量表示,我们有如下的定理
定理 8:
α ( t ) = π P ( x 1 ) A P ( x 2 ) ⋯ A P ( x t ) = π P ( x 1 ) ⋅ s = 2 ∏ t A P ( x s )
直观地看,上式的左边是一个 1 × m 的行向量,右边是 ( 1 × m ) ⋅ ( m × m ) ⋯ ( m × m ) 的矩阵乘法,结果也是一个 1 × m 的行向量。
定理 9 [前向概率的初值条件 (Initial Condition of Forward Probability)]:
α i ( 1 ) = π i ⋅ p i ( x 1 ) ( i = 1 , 2 , … , m )
矩阵(向量)形式为
α ( 1 ) = π P ( x 1 )
证明:
α i ( 1 ) = P ( X 1 = x 1 , C 1 = i ) = P ( X 1 = x 1 ∣ C 1 = i ) ⋅ P ( C 1 = i ) = π i ⋅ p i ( x 1 ) (取条件于 C 1 )
矩阵形式下,可写为
[ α 1 ( 1 ) , α 2 ( 1 ) , … , α m ( 1 ) ] 1 × m = [ π 1 , π 2 , … , π m ] 1 × m ⋅ p 1 ( x 1 ) p 2 ( x 1 ) ⋱ p m ( x 1 ) m × m
也即
α ( 1 ) = π P ( x 1 )
证毕
定理 10 [前向概率的递归公式 (Recursion of Forward Probability)]:
α j ( t + 1 ) = ( i = 1 ∑ m α i ( t ) ⋅ a ij ) ⋅ p j ( x t + 1 ) ( t = 1 , 2 , … , T − 1 ; j = 1 , 2 , … , m )
矩阵形式为
α ( t + 1 ) = α ( t ) ⋅ A ⋅ P ( x t + 1 ) ( t = 1 , 2 , … , T − 1 )
证明:
(1) 实数形式 (scalar form)
α j ( t + 1 ) = P ( X 1 t + 1 = x 1 t + 1 , C t + 1 = j ) = P ( X 1 t + 1 = x 1 t + 1 , C t + 1 = j , i ⋃ ( C t = i ) ) = i ∑ P ( X 1 t + 1 = x 1 t + 1 , C t = i , C t + 1 = j ) = i ∑ P ( X t + 1 = x t + 1 , C t + 1 = j ∣ X 1 t = x 1 t , C t = i ) ⋅ P ( X 1 t = x 1 t , C t = i ) α i ( t ) (取条件于 X 1 t , C t ) = i ∑ P ( X t + 1 = x t + 1 , C t + 1 = j ∣ C t = i ) ⋅ α i ( t ) = i ∑ α i ( t ) ⋅ P ( C t + 1 = j ∣ C t = i ) ⋅ P ( X t + 1 = x t + 1 ∣ C t = i , C t + 1 = j ) (基于引理 1 ) = i ∑ α i ( t ) ⋅ a ij ⋅ P ( X t + 1 = x t + 1 ∣ C t + 1 = j ) = ( i ∑ α i ( t ) ⋅ a ij ) ⋅ p j ( x t + 1 ) (基于 HMM 性质 1 、 2 ) (基于 HMM 性质 2 )
(2) 矩阵形式 (matrix form)
α ( t + 1 ) = π P ( x 1 ) ⋅ s = 2 ∏ t + 1 A P ( x s ) (基于定理 8 )
= π P ( x 1 ) ⋅ s = 2 ∏ t A P ( x s ) ⋅ A P ( x t + 1 )
= π P ( x 1 ) ⋅ s = 2 ∏ t A P ( x s ) α ( t ) ⋅ A P ( x t + 1 )
(3) 实际上,我们也可以从矩阵形式得到实数形式
矩阵形式 α ( t + 1 ) = α ( t ) ⋅ A ⋅ P ( x t + 1 ) 说明
[ α 1 ( t + 1 ) , … , α j ( t + 1 ) , … , α m ( t + 1 ) ] 1 × m = [ α 1 ( t ) , … , α j ( t ) , … , α m ( t ) ] 1 × m ⋅ a 11 a 21 ⋮ a m 1 ⋯ ⋯ ⋱ ⋯ a 1 j a 2 j ⋮ a mj ⋯ ⋯ ⋱ ⋯ a 1 m a 2 m ⋮ a mm m × m ⋅ p 1 ( x t + 1 ) ⋱ p j ( x t + 1 ) ⋱ p m ( x t + 1 ) m × m = [ i ∑ α i ( t ) ⋅ a i 1 , … , i ∑ α i ( t ) ⋅ a ij , … , i ∑ α i ( t ) ⋅ a im ] 1 × m ⋅ p 1 ( x t + 1 ) ⋱ p j ( x t + 1 ) ⋱ p m ( x t + 1 ) m × m
等式左边的第 j 个元素为 α j ( t + 1 )
等式右边的第 j 个元素为 ( i ∑ α i ( t ) ⋅ a ij ) ⋅ p j ( x t + 1 )
矩阵相等意味着矩阵每个对应元素均相等,因此,矩阵形式可以得出
α j ( t + 1 ) = ( i ∑ α i ( t ) ⋅ a ij ) ⋅ p j ( x t + 1 )
(也即是实数形式) 证毕
前向概率的一个重要应用是隐马尔可夫链的似然函数可以用它来表示,如下述定理所示。
定理 11 [似然函数的前向概率表达]:
L = P ( X 1 T = x 1 T ) = i = 1 ∑ m α i ( T ) = α ( T ) ⋅ 1 T
证明:
L = P ( X 1 T = x 1 T ) = P ( X 1 T = x 1 T , i ⋃ ( C T = i ) ) = i = 1 ∑ m P ( X 1 T = x 1 T , C T = i ) α i ( T ) = i = 1 ∑ m α i ( T ) = [ α 1 ( T ) , α 2 ( T ) , … , α m ( T ) ] 1 × m ⋅ 1 1 ⋮ 1 m × 1 = α ( T ) ⋅ 1 T
实际上,由定理 6,我们有
L = π P ( x 1 ) A P ( x 2 ) ⋯ A P ( x T ) ⋅ 1 T
此外,由定理 8,我们又有
α ( T ) = π P ( x 1 ) ⋅ s = 2 ∏ T A P ( x s ) = π P ( x 1 ) A P ( x 2 ) ⋯ A P ( x T )
因此,很明显的
L = α ( T ) ⋅ 1 T
对于前向概率的论述到此可以暂时告一段落,我们将给出利用前向概率求全部观测值概率(似然函数)的向前算法。
向前算法(forward approach)
向前算法用来求概率 P ( X 1 T = x 1 T ∣ θ ) ,也就是似然函数的值。它是一个从初值开始,经过迭代最终到达终值的过程,这就是所谓的“向前”。
α i ( 1 ) = π i ⋅ p i ( x 1 ) ( i = 1 , 2 , ⋯ , m )
或者等价的
α ( 1 ) = π P ( x 1 )
α j ( t + 1 ) = ( i = 1 ∑ m α i ( t ) ⋅ a ij ) ⋅ p j ( x t + 1 ) ( t = 1 , 2 , ⋯ , T − 1 , j = 1 , 2 , ⋯ , m )
或者等价的
α ( t + 1 ) = α ( t ) ⋅ A ⋅ P ( x t + 1 ) ( t = 1 , 2 , ⋯ , T − 1 )
L = P ( X 1 T = x 1 T ) = i = 1 ∑ m α i ( T ) = α ( T ) ⋅ 1 T
后向概率与向后算法
定义 6 [后向概率(Backward Probability)]:
β i ( t ) = P ( X t + 1 T = x t + 1 T ∣ C t = i ) ( t = 1 , 2 , ⋯ , T − 1 , i = 1 , 2 , ⋯ , m )
与前向概率是一个联合概率不同,后向概率是一个条件概率,它是给定 t 时刻状态位于 i 的条件下,出现 x t + 1 , x t + 2 , ⋯ , x T 的概率。
与前向概率类似,我们也可以定义后向概率的向量表示,即
β ( t ) = [ β 1 ( t ) , β 2 ( t ) , ⋯ , β m ( t ) ] 1 × m
定理 12:
β ( t ) T = A P ( x t + 1 ) A P ( x t + 2 ) ⋯ A P ( x T ) ⋅ 1 T = [ s = t + 1 ∏ T A P ( x s ) ] ⋅ 1 T ( t = 1 , 2 , ⋯ , T − 1 )
直观地来看,上式的左边为 β ( t ) T ,是一个 m × 1 的列向量,上式的右边是 ( m × m ) ⋅ ( m × m ) ⋯ ( m × 1 ) 的矩阵运算,结果也是一个 m × 1 的列向量。
定理 13 [后向概率的终值条件(Terminal Condition of Backward Probability)]:
β i ( T ) = 1 ( i = 1 , 2 , ⋯ , m )
矩阵形式为
β ( T ) = 1 1 × m
定理 14:
β i ( t ) = j ∑ P ( X t + 1 T = x t + 1 T ∣ C t + 1 = j ) ⋅ a ij
证明:
RHS = j ∑ P ( X t + 1 T = x t + 1 T ∣ C t = i , C t + 1 = j ) ⋅ P ( C t + 1 = j ∣ C t = i ) = j ∑ P ( C t = i , C t + 1 = j ) P ( X t + 1 T = x t + 1 T , C t = i , C t + 1 = j ) ⋅ P ( C t = i ) P ( C t = i , C t + 1 = j ) (基于 HMM 性质 2 ) = j ∑ P ( C t = i ) P ( X t + 1 T = x t + 1 T , C t = i , C t + 1 = j ) = j ∑ P ( X t + 1 T = x t + 1 T , C t + 1 = j ∣ C t = i ) = P ( X t + 1 T = x t + 1 T , j ⋃ ( C t + 1 = j ) C t = i ) = P ( X t + 1 T = x t + 1 T ∣ C t = i ) = β i ( t ) = LHS
从定理 14 出发,我们可以得到后向概率的递归公式。
定理 15 [后向概率的递归公式(Recursion of Backward Probability)]:
β i ( t ) = j = 1 ∑ m p j ( x t + 1 ) ⋅ β j ( t + 1 ) ⋅ a ij ( t = T − 1 , ⋯ , 2 , 1 , i = 1 , 2 , ⋯ , m )
矩阵形式可写为
β ( t ) T = A ⋅ P ( x t + 1 ) ⋅ β ( t + 1 ) T ( t = T − 1 , ⋯ , 2 , 1 )
证明:
(1) 实数形式(scalar form)
β i ( t ) = j ∑ P ( X t + 1 T = x t + 1 T ∣ C t + 1 = j ) ⋅ a ij (基于定理 14 ) = j ∑ a ij ⋅ P ( X t + 1 = x t + 1 ∣ C t + 1 = j ) ⋅ P ( X t + 2 T = x t + 2 T ∣ X t + 1 = x t + 1 , C t + 1 = j ) (基于引理 1 ) = j ∑ a ij ⋅ p j ( x t + 1 ) ⋅ P ( X t + 2 T = x t + 2 T ∣ C t + 1 = j ) (基于 HMM 性质 2 ) = j ∑ a ij ⋅ p j ( x t + 1 ) ⋅ β j ( t + 1 )
(2) 矩阵形式(matrix form)
β ( t ) T = [ s = t + 1 ∏ T A P ( x s ) ] ⋅ 1 T (基于定理 12 )
= A P ( x t + 1 ) ⋅ [ s = t + 2 ∏ T A P ( x s ) ] ⋅ 1 T = A ⋅ P ( x t + 1 ) ⋅ β ( t + 1 ) T
(3) 与前向概率的递归公式类似,这里我们也可以从矩阵形式得到实数形式。
矩阵形式 β ( t ) T = A ⋅ P ( x t + 1 ) ⋅ β ( t + 1 ) T 意味着
β 1 ( t ) ⋮ β i ( t ) ⋮ β m ( t ) m × 1 = a 11 ⋮ a i 1 ⋮ a m 1 ⋯ ⋱ ⋯ ⋱ ⋯ a 1 j ⋮ a ij ⋮ a mj ⋯ ⋱ ⋯ ⋱ ⋯ a 1 m ⋮ a im ⋮ a mm m × m ⋅ p 1 ( x t + 1 ) ⋱ p j ( x t + 1 ) ⋱ p m ( x t + 1 ) m × m ⋅ β 1 ( t + 1 ) ⋮ β j ( t + 1 ) ⋮ β m ( t + 1 ) m × 1
= a 11 ⋮ a i 1 ⋮ a m 1 ⋯ ⋱ ⋯ ⋱ ⋯ a 1 j ⋮ a ij ⋮ a mj ⋯ ⋱ ⋯ ⋱ ⋯ a 1 m ⋮ a im ⋮ a mm m × m ⋅ p 1 ( x t + 1 ) β 1 ( t + 1 ) ⋮ p j ( x t + 1 ) β j ( t + 1 ) ⋮ p m ( x t + 1 ) β m ( t + 1 ) m × 1
LHS 的第 i 个元素为 β i ( t )
RHS 的第 i 个元素为
a i 1 ⋅ p 1 ( x t + 1 ) β 1 ( t + 1 ) + ⋯ + a ij ⋅ p j ( x t + 1 ) β j ( t + 1 ) + ⋯ + a im ⋅ p m ( x t + 1 ) β m ( t + 1 ) = j = 1 ∑ m a ij ⋅ p j ( x t + 1 ) ⋅ β j ( t + 1 )
矩阵相等意味着 LHS 与 RHS 的每一个元素均相等,即
β i ( t ) = j = 1 ∑ m p j ( x t + 1 ) ⋅ β j ( t + 1 ) ⋅ a ij
也就是实数形式。
与前向概率类似,隐马尔可夫链的似然函数也可以用后向概率来表示。
定理 16 [似然函数的后向概率表达]:
L = P ( X 1 T = x 1 T ) = i = 1 ∑ m π i ⋅ p i ( x 1 ) ⋅ β i ( 1 )
矩阵形式可写为
L = π ⋅ P ( x 1 ) ⋅ β ( 1 ) T
证明:
1. 实数形式
L = P ( X 1 T = x 1 T ) = P ( X 1 T = x 1 T , i ⋃ ( C 1 = i ) ) = i ∑ P ( X 1 T = x 1 T , C 1 = i ) = i ∑ P ( X 1 T = x 1 T ∣ C 1 = i ) ⋅ P ( C 1 = i ) = i ∑ π i ⋅ P ( X 1 = x 1 ∣ C 1 = i ) ⋅ P ( X 2 T = x 2 T ∣ C 1 = i , X 1 = x 1 ) (取条件于 C 1 ) = i ∑ π i ⋅ p i ( x 1 ) ⋅ P ( X 2 T = x 2 T ∣ C 1 = i ) (基于引理 1 ) = i = 1 ∑ m π i ⋅ p i ( x 1 ) ⋅ β i ( 1 ) (基于 HMM 性质 2 )
2. 矩阵形式
以上的实数形式可用矩阵运算写为
L = [ π 1 , π 2 , ⋯ , π m ] 1 × m ⋅ p 1 ( x 1 ) p 2 ( x 1 ) ⋱ p m ( x 1 ) m × m ⋅ β 1 ( 1 ) β 2 ( 1 ) ⋮ β m ( 1 ) m × 1 = π ⋅ P ( x 1 ) ⋅ β ( 1 ) T
实际上,由定理 6,我们有
L = π P ( x 1 ) A P ( x 2 ) ⋯ A P ( x T ) ⋅ 1 T
由定理 12,我们又有
β ( 1 ) T = [ s = 2 ∏ T A P ( x s ) ] ⋅ 1 T = A P ( x 2 ) ⋯ A P ( x T ) ⋅ 1 T
因此,很明显的有
L = π ⋅ P ( x 1 ) ⋅ β ( 1 ) T
似然函数不但可以用单独的前向概率或者后向概率来表示,也可以把前向和后向概率结合在一起来表达似然函数,基于下面的两个定理。
定理 17:
P ( X 1 T = x 1 T , C t = i ) = α i ( t ) ⋅ β i ( t ) ( t = 1 , 2 , ⋯ , T , i = 1 , 2 , ⋯ , m )
证明:
P ( X 1 T = x 1 T , C t = i ) = P ( X t + 1 T = x t + 1 T ∣ X 1 t = x 1 t , C t = i ) ⋅ P ( X 1 t = x 1 t , C t = i ) (取条件于 C t , X 1 t ) = α i ( t ) ⋅ P ( X t + 1 T = x t + 1 T ∣ C t = i ) = α i ( t ) ⋅ β i ( t ) (基于 HMM 性质 2 )
定理 18 [似然函数的前向/后向概率表达]:
L = i = 1 ∑ m α i ( t ) ⋅ β i ( t )
矩阵形式为
L = α ( t ) ⋅ β ( t ) T
证明:
L = P ( X 1 T = x 1 T ) = P ( X 1 T = x 1 T , i ⋃ ( C t = i ) ) = i ∑ P ( X 1 T = x 1 T , C t = i ) = i ∑ α i ( t ) ⋅ β i ( t ) (基于定理 17 )
利用矩阵运算,上式可写为
L = [ α 1 ( t ) , α 2 ( t ) , ⋯ , α m ( t ) ] 1 × m ⋅ β 1 ( t ) β 2 ( t ) ⋮ β m ( t ) m × 1 = α ( t ) ⋅ β ( t ) T
一个更直观地证明上述矩阵形式的方法是,由定理 6 可得
L = α ( t ) π P ( x 1 ) A P ( x 2 ) ⋯ A P ( x t ) ⋅ β ( t ) T A P ( x t + 1 ) ⋯ A P ( x T ) ⋅ 1 T = α ( t ) ⋅ β ( t ) T (基于定理 8 和定理 12 )
仿照前向概率的部分,这里我们先给出利用后向概率求全部观测值概率(似然函数)的向后算法的框架,然后给出向后算法的细节及伪代码实现。附录中的 Matlab 程序把向前和向后算法统一起来,在一个函数中加以实现。
向后算法(backward approach)
向后算法也是用来求概率 P ( X 1 T = x 1 T ∣ θ ) ,也就是似然函数的值。它是一个从终值开始,经过迭代最终到达初值的过程,这就是所谓的“向后”。
β i ( T ) = 1 ( i = 1 , 2 , ⋯ , m )
或者等价的
β ( T ) = 1
β i ( t ) = j = 1 ∑ m a ij ⋅ p j ( x t + 1 ) ⋅ β j ( t + 1 ) ( t = T − 1 , ⋯ , 2 , 1 , i = 1 , 2 , ⋯ , m )
或者等价的
β ( t ) T = A ⋅ P ( x t + 1 ) ⋅ β ( t + 1 ) T ( t = T − 1 , ⋯ , 2 , 1 )
L = P ( X 1 T = x 1 T ) = i = 1 ∑ m π i ⋅ p i ( x 1 ) ⋅ β i ( 1 )
或者等价的
L = π ⋅ P ( x 1 ) ⋅ β ( 1 ) T
其他数值参量
到现在为止,我们已经介绍了前向/后向概率和向前/向后算法,也就是说,隐马尔可夫链的第一个基本问题已经回答完毕了。不过,在回答第二个、第三个基本问题之前,我们还是希望再论述一下其他相关的数值参量(other Desired Quantities)。一方面,这些数值参量与向前/向后概率有很密切的关系;另一方面,它们也是后面介绍的 Baum-Welch 算法的基础。
γ 与 g
定义 7:
γ i ( t ) = P ( C t = i ∣ X 1 T = x 1 T ) ( t = 1 , 2 , ⋯ , T , i = 1 , 2 , ⋯ , m )
从定义 7 可以看出,γ i ( t ) 是给定所有观测值 x 1 , x 2 , ⋯ , x T 的条件下,状态变量在 t 时刻位于 i 的条件概率。之前我们曾经定义过 u i ( t ) = P ( C t = i ) 来表示状态在 t 时刻位于 i 的无条件概率,因此,γ i ( t ) 可以被看作 u i ( t ) 的条件估计(conditional guess/estimation),给定全部观测值这个条件。
与之前的做法类似,这里我们也用 γ ( t ) 来作为它的向量表示
γ ( t ) = [ γ 1 ( t ) , γ 2 ( t ) , ⋯ , γ m ( t ) ] 1 × m
定理 19:
γ i ( t ) = L α i ( t ) ⋅ β i ( t )
证明:
γ i ( t ) = P ( C t = i ∣ X 1 T = x 1 T ) = P ( X 1 T = x 1 T ) P ( X 1 T = x 1 T , C t = i ) = L α i ( t ) ⋅ β i ( t ) (分子基于定理 17 )
接下来我们定义 g
定义 8
g i = t = 1 ∑ T − 1 γ i ( t ) ( i = 1 , 2 , ⋯ , m )
很自然的,它的向量表示为
g = [ g 1 , g 2 , ⋯ , g m ] = [ t = 1 ∑ T − 1 γ 1 ( t ) , t = 1 ∑ T − 1 γ 2 ( t ) , ⋯ , t = 1 ∑ T − 1 γ m ( t ) ] 1 × m
g i 可不仅仅是把 T − 1 个 γ i ( t ) 相加这么简单,它有很重要的性质
定理 20:
g i 是
(1) 给定所有观测值的条件下,状态变量位于 i 的期望次数;
(2) 给定所有观测值的条件下,状态变量从 i 转移出去的期望次数(在最后时刻 T 无转移)
证明:
构造示性函数 1 i ( t ) 如下
1 i ( t ) = { 1 0 (给定所有观测值,若状态变量在 t 时刻位于 i ) (其他情况)
则有
g i = t = 1 ∑ T − 1 γ i ( t ) = t = 1 ∑ T − 1 E [ 1 i ( t )] = E [ t = 1 ∑ T − 1 1 i ( t ) ]
因此,g i 是位于 i 的期望次数,也是从 i 转移出去的期望次数,在给定所有观测值的条件下。
ε 与 h
定义 9:
ε i , j ( t ) = P ( C t = i , C t + 1 = j ∣ X 1 T = x 1 T ) ( t = 1 , 2 , ⋯ , T − 1 , i , j = 1 , 2 , ⋯ , m )
从上述定义可知,ε i , j ( t ) 是给定全部观测值 x 1 , x 2 , ⋯ , x T 的条件下,状态变量在 t 时刻位于 i ,在 t + 1 时刻转移到 j 的条件概率。与 γ i ( t ) 的情况类似,它也可被视为相应的无条件概率的条件估计。
ε i , j ( t ) 也可以用 α i ( t ) , β i ( t ) 等参量来表示,不过我们先要给出下面这个定理。
定理 21:
P ( X 1 T = x 1 T , C t = i , C t + 1 = j ) = α i ( t ) ⋅ a ij ⋅ p j ( x t + 1 ) ⋅ β j ( t + 1 )
证明:
P ( X 1 T = x 1 T , C t = i , C t + 1 = j ) = P ( X t + 1 T = x t + 1 T , C t + 1 = j ∣ X 1 t = x 1 t , C t = i ) ⋅ P ( X 1 t = x 1 t , C t = i ) (取条件于 X 1 t , C t ) = α i ( t ) ⋅ P ( X t + 1 T = x t + 1 T , C t + 1 = j ∣ C t = i ) = α i ( t ) ⋅ P ( C t + 1 = j ∣ C t = i ) ⋅ P ( X t + 1 T = x t + 1 T ∣ C t = i , C t + 1 = j ) (基于 HMM 性质 1 、 2 ) = α i ( t ) ⋅ a ij ⋅ P ( X t + 1 T = x t + 1 T ∣ C t + 1 = j ) (基于引理 1 ) = α i ( t ) ⋅ a ij ⋅ P ( X t + 1 = x t + 1 ∣ C t + 1 = j ) ⋅ P ( X t + 2 T = x t + 2 T ∣ C t + 1 = j , X t + 1 = x t + 1 ) (基于 HMM 性质 2 ) = α i ( t ) ⋅ a ij ⋅ p j ( x t + 1 ) ⋅ P ( X t + 2 T = x t + 2 T ∣ C t + 1 = j ) (基于引理 1 ) = α i ( t ) ⋅ a ij ⋅ p j ( x t + 1 ) ⋅ β j ( t + 1 ) (基于 HMM 性质 2 )
定理 22:
ε i , j ( t ) = L α i ( t ) ⋅ a ij ⋅ p j ( x t + 1 ) ⋅ β j ( t + 1 ) ( t = 1 , 2 , ⋯ , T − 1 , i , j = 1 , 2 , ⋯ , m )
证明:
ε i , j ( t ) = P ( C t = i , C t + 1 = j ∣ X 1 T = x 1 T ) = P ( X 1 T = x 1 T ) P ( X 1 T = x 1 T , C t = i , C t + 1 = j ) = L α i ( t ) ⋅ a ij ⋅ p j ( x t + 1 ) ⋅ β j ( t + 1 ) (基于定理 21 )
最后一个定义的数值参量是 h ij ,即
定义 10:
h ij = t = 1 ∑ T − 1 ε i , j ( t ) ( i , j = 1 , 2 , ⋯ , m )
与前面的向量表示类似,这里我们也可以写出 h ij 的矩阵表示(Matrix Representation)
H = h 11 ⋮ h m 1 ⋯ ⋱ ⋯ h 1 m ⋮ h mm m × m
与 g i 类似,h ij 也有一个很重要的性质。
定理 23:
h ij 是在给定所有观测值的条件下,状态变量从 i 转移到 j 的期望次数
证明:
仍引入示性函数
1 i , j ( t ) = { 1 0 (状态变量 t 时刻位于 i , t + 1 时刻转移到 j ,给定所有观测值) (其他情况)
因此
h ij = t = 1 ∑ T − 1 ε i , j ( t ) = t = 1 ∑ T − 1 E [ 1 i , j ( t )] = E [ t = 1 ∑ T − 1 1 i , j ( t ) ]
故 h ij 是在给定所有观测值的条件下,状态变量从 i 转移到 j 的期望次数。
期望最大化算法
我们将先跳过第二个问题,直接回答第三个问题,即参数估计。参数估计是隐马尔可夫链三个问题中最难的也是最核心的,解决这个难题一般用期望最大化算法(Expectation Maximization, EM),在隐马尔可夫链框架下也称为 Baum-Welch 算法。
期望最大化算法的基本思想
期望最大化算法可视为极大似然估计方法的拓展,主要用来解决含有不可观测变量的参数估计问题。假设在我们的统计模型中,有两组变量,x 1 , x 2 , ⋯ , x T 是可观测变量(observed),y 1 , y 2 , ⋯ , y T 是不可观测变量(unobserved or hidden),则 { x 1 , x 2 , ⋯ , x T ; y 1 , y 2 , ⋯ , y T } 被称为完整数据(complete data),而完整数据的似然函数(complete data likelihood function)为
L ( θ ∣ x 1 , x 2 , ⋯ , x T ; y 1 , y 2 , ⋯ , y T ) = p ( x 1 , x 2 , ⋯ , x T ; y 1 , y 2 , ⋯ , y T ∣ θ ) = p ( x 1 T ∣ θ ) ⋅ p ( y 1 T ∣ x 1 T , θ )
评注: 当含有不可观测变量时,似然函数是随机的,因为 y 1 , y 2 , ⋯ , y T 是随机数据(不知道确切数值),故我们可以把似然函数看成当 x 1 , x 2 , ⋯ , x T 给定时,关于 y 1 , y 2 , ⋯ , y T 的函数。
相应的,
L ( θ ∣ x 1 , x 2 , ⋯ , x T ) = p ( x 1 T ∣ θ )
被称为不完整数据的似然函数(incomplete data likelihood function)。
在完整数据的似然函数基础上,我们定义完整数据的对数似然函数(complete data log-likelihood function, CDLL),这个 CDLL 是我们今后要处理的主要对象
log L ( θ ∣ x 1 T , y 1 T ) = log p ( x 1 T , y 1 T ∣ θ ) = log [ p ( x 1 T ∣ θ ) ⋅ p ( y 1 T ∣ x 1 T , θ ) ]
当然,这个 CDLL 也是一个随机变量,因为 y 1 , y 2 , ⋯ , y T 是随机的。它也仍可被视为关于 y 1 , y 2 , ⋯ , y T 的函数,当 x 1 , x 2 , ⋯ , x T 给定时。
期望最大化算法(EM)如其名所示,包含两个步骤:求期望(E-Step)与最大化(M-Step)。下面我们分别加以论述。
E-Step
E-Step 的核心思想是求完整数据对数似然函数(CDLL)的条件期望。这个条件期望是对于随机数据 y 1 , y 2 , ⋯ , y T 而言的,给定观测值 x 1 , x 2 , ⋯ , x T 和目前的参数估计值 θ ~ 。
定义 11 [CDLL 的条件期望,Q 函数]:
Q ( θ , θ ) = E [ log p ( x 1 T , y 1 T ∣ θ ) x 1 T , θ ]
注: 在 Q 函数中
x 1 T 是常数,在计算条件期望时使用
θ 是一个正常变量,有条件分布 p ( y 1 T ∣ x 1 T , θ ) 或 f ( y 1 T ∣ x 1 T , θ ~ )
y 1 T 是随机数据,有条件分布 p ( y 1 T ∣ x 1 T , θ ~ )
因此,我们可以计算 Q 函数
Q ( θ , θ ) = E [ log p ( x 1 T , y 1 T ∣ θ ) x 1 T , θ ]
M-Step
M-Step 主要是求使得 Q 函数(条件期望)最大化的那个 θ ^ ,用数学语言表示就是
θ ^ = arg θ max Q ( θ , θ ~ )
E-Step 与 M-Step 交替进行,直到满足一定条件时为止,故 EM 算法本质上是一个数值迭代的算法。
Baum-Welch 算法
Baum-Welch 算法是在隐马尔可夫链环境下期望最大化算法的特殊形式。在泊松—隐马尔可夫链(Poisson-HMM)环境下,θ = [ π , A , λ ] ;在正态—隐马尔可夫链(Normal-HMM)环境下,θ = [ π , A , μ , σ ] 。隐马尔可夫链的完整数据对数似然函数(CDLL)为
log L ( θ ∣ x 1 T , c 1 T ) = log p ( x 1 T , c 1 T ∣ θ )
这里 c 1 , ⋯ , c T 是随机数据,相当于上文中提到的不可观测变量 y 1 , ⋯ , y T 。
基于定理 4,有
p ( x 1 T , c 1 T ) = π c 1 ⋅ t = 1 ∏ T − 1 a c t , c t + 1 ⋅ t = 1 ∏ T p c t ( x t )
故 HMM 的 CDLL 可写为
log π c 1 + t = 1 ∑ T − 1 log a c t , c t + 1 + t = 1 ∑ T log p c t ( x t )
因此在隐马尔可夫链环境下,Q 函数可表示为
Q = E [ log p ( x 1 T , c 1 T ∣ θ ) ⋅ p ( c 1 T ∣ x 1 T , θ ) ] = c 1 , ⋯ , c T ∑ log p ( x 1 T , c 1 T ∣ θ ) ⋅ p ( c 1 T ∣ x 1 T , θ ) = CDLL c 1 , ⋯ , c T ∑ log π c 1 ⋅ p ( c 1 T ∣ x 1 T , θ ) + CDLL c 1 , ⋯ , c T ∑ ( t = 1 ∑ T − 1 log a c t , c t + 1 ) ⋅ p ( c 1 T ∣ x 1 T , θ ) + CDLL c 1 , ⋯ , c T ∑ ( t = 1 ∑ T log p c t ( x t ) ) ⋅ p ( c 1 T ∣ x 1 T , θ ~ ) = ( a ) + ( b ) + ( c ) [ 基于式 (4.3)]
可见,Q 函数可以被分解为三项之和,而这三项我们又可以分别处理。
(1) 对于第一项,我们可仅关注于第 1 个时间层,即仅考虑 t = 1 时的边缘分布,因此
( a ) = i = 1 ∑ m log π i ⋅ P ( C 1 = i ∣ X 1 T = x 1 T ) = i = 1 ∑ m log π i ⋅ γ i ( 1 )
(2) 对于第二项,我们仅关注在每个时间层,从 i 到 j 的转移,
( b ) = i = 1 ∑ m j = 1 ∑ m t = 1 ∑ T − 1 log a ij ⋅ P ( C t = i , C t + 1 = j ∣ X 1 T = x 1 T ) = i = 1 ∑ m j = 1 ∑ m t = 1 ∑ T − 1 log a ij ⋅ ε i , j ( t )
(3) 对于第三项,我们仅关注在每个时间层,在所有状态下的发射,
= i = 1 ∑ m t = 1 ∑ T log p i ( x t ) ⋅ P ( C t = i ∣ X 1 T = x 1 T ) = i = 1 ∑ m t = 1 ∑ T log p i ( x t ) ⋅ γ i ( t )
故在隐马尔可夫链环境下,我们可以把 Q 函数简化为
Q ( θ , θ ~ ) = i = 1 ∑ m log π i ⋅ γ i ( 1 ) + i = 1 ∑ m j = 1 ∑ m t = 1 ∑ T − 1 log a ij ⋅ ε i , j ( t ) + i = 1 ∑ m t = 1 ∑ T log p i ( x t ) ⋅ γ i ( t ) = ( A ) + ( B ) + ( C )
Baum-Welch 算法仍是分为 E-Step 和 M-Step 两大步骤(其本质上仍是 EM 算法)。
E-Step
给定观测值序列 x 1 , x 2 , … , x T ,给定当前的参数估计值 ( θ = [ π , A , λ ] or θ = [ π , A , μ , σ ~ ]) 的情况下,计算 α i ( t ) , β i ( t ) , p i ( x t ) , γ i ( t ) , ε i , j ( t ) , g i , h ij , L 等参量值,则 ( A ) , ( B ) , ( C ) 就可以算出,继而 Q 函数的值(条件期望)也可以算出,这也就是 E-Step 的意义。
M-Step
求使 Q 函数(条件期望)最大化的参数值 ( θ ^ = [ π ^ , A ^ , λ ^ ] or θ ^ = [ π ^ , A ^ , μ ^ , σ ^ ]) 。因为我们已经把 Q 函数分解为三个部分,且可以观察到,( A ) 仅与 π (初始分布)相关,( B ) 仅与 A (转移矩阵)相关,( C ) 仅与 p i ( x t ) (发射参数)相关,则我们可以分别最优化 ( A ) , ( B ) , ( C ) ,继而得到 Baum-Welch 公式。
首先我们对 ( A ) 进行最优化,得到 π ^ i 的最优估计量
定理 24 [初始分布的最优估计]:
π ^ i = γ i ( 1 )
证明:
最优化模型为
⎩ ⎨ ⎧ max i = 1 ∑ m log π i ⋅ γ i ( 1 ) s.t. i = 1 ∑ m π i = 1
其拉格朗日函数 (Lagrange Function) 为
L = i = 1 ∑ m log π i ⋅ γ i ( 1 ) + λ ( 1 − π 1 − π 2 − ⋯ − π m )
∂ π i ∂ L = 0
⇒ π i γ i ( 1 ) − λ = 0
⇒ π ^ i = λ γ i ( 1 )
此外
i ∑ π i = 1
⇒ i ∑ π ^ i = 1
⇒ i ∑ λ γ i ( 1 ) = 1
⇒ λ = i ∑ γ i ( 1 )
同时我们注意到 γ i ( 1 ) 是一个(条件)概率,有
i ∑ γ i ( 1 ) = i ∑ P ( C 1 = i ∣ X 1 T = x 1 T ) = P ( i ⋃ ( C 1 = i ) ∣ X 1 T = x 1 T ) = 1
因此 λ = ∑ i γ i ( 1 ) = 1
综合以上信息,我们可以得出
π ^ i = γ i ( 1 )
接下来我们对 ( B ) 进行最优化,得到 a ij 的最优估计量。
定理 25 [转移概率的最优估计]:
a ^ ij = t = 1 ∑ T − 1 γ i ( t ) t = 1 ∑ T − 1 ε i , j ( t ) = g i h ij
证明:
可以知道,马尔可夫链转移概率的最优估计是
a ^ ij = k = 1 ∑ m f ik f ij
其中 f ij 是从状态 i 到状态 j 的转移次数,k = 1 ∑ m f ik 是从状态 i 转移出去的次数之和。
由定理 23 可知,h ij [ 即 t = 1 ∑ T − 1 ε i , j ( t ) ] 是给定所有观测值的条件下,从 i 转移到 j 的期望次数;由定理 20 又可以知道,g i [ 即 t = 1 ∑ T − 1 γ i ( t ) ] 是给定所有观测值的条件下,从 i 转移出去的期望次数。
因此,转移概率的最优估计是
a ^ ij = t = 1 ∑ T − 1 γ i ( t ) t = 1 ∑ T − 1 ε i , j ( t ) = g i h ij
证毕
最后我们对 ( C ) 进行最优化,得到发射参数的最优估计量。对于不同的隐马尔可夫链而言,前两个最优估计量 ( π ^ i , a ^ ij ) 均是相同的,但不同的 HMM 其发射分布不同,因此其发射参数的最优估计量也并不相同。对于泊松-隐马尔可夫链和正态-隐马尔可夫链来说,其发射参数的最优估计量有明晰的解析表达式。
定理 26 [泊松-隐马尔可夫链的发射参数最优估计]:
λ ^ i = t = 1 ∑ T γ i ( t ) t = 1 ∑ T γ i ( t ) ⋅ x t
证明:
在 Poisson-HMM 环境下,其发射分布是
p i ( x t ) = x t ! e − λ i ⋅ λ i x t
我们需要最大化
L = i = 1 ∑ m t = 1 ∑ T log p i ( x t ) ⋅ γ i ( t ) = i ∑ t ∑ ( − λ i + x t ⋅ log λ i − log x t ! ) ⋅ γ i ( t )
因此
∂ λ i ∂ L = 0
⇒ t ∑ ( − 1 + λ i x t ) ⋅ γ i ( t ) = 0
⇒ λ i t ∑ x t ⋅ γ i ( t ) = t ∑ γ i ( t )
⇒ λ ^ i = t = 1 ∑ T γ i ( t ) t = 1 ∑ T γ i ( t ) ⋅ x t
注: 这里的分母是 t = 1 ∑ T γ i ( t ) ,不是 g i ,因为 g i = t = 1 ∑ T − 1 γ i ( t ) 。
定理 27 [正态-隐马尔可夫链的发射参数最优估计]:
μ ^ i σ ^ i 2 = t = 1 ∑ T γ i ( t ) t = 1 ∑ T γ i ( t ) ⋅ x t = t = 1 ∑ T γ i ( t ) t = 1 ∑ T γ i ( t ) ⋅ ( x t − μ ^ i ) 2
证明:
在 Normal-HMM 环境下,其发射分布是
p i ( x t ) = 2 π σ i 1 e − 2 σ i 2 ( x t − μ i ) 2
我们需要最大化
L = i = 1 ∑ m t = 1 ∑ T log p i ( x t ) ⋅ γ i ( t ) = i ∑ t ∑ log [ ( 2 π ) − 2 1 ⋅ ( σ i 2 ) − 2 1 ⋅ e − 2 σ i 2 ( x t − μ i ) 2 ] ⋅ γ i ( t ) = i ∑ t ∑ ( − 2 1 log 2 π − 2 1 log σ i 2 − 2 σ i 2 ( x t − μ i ) 2 ) ⋅ γ i ( t )
分别对 μ i , σ i 求偏导并令偏导数为 0 可得
∂ μ i ∂ L ⇒ t ∑ γ i ( t ) ⋅ ( − 2 σ i 2 1 ⋅ 2 ( x t − μ i ) ⋅ ( − 1 ) ) ⇒ t ∑ γ i ( t ) ⋅ ( x t − μ i ) ⇒ t ∑ γ i ( t ) ⋅ x t ⇒ μ ^ i = 0 = 0 = 0 = ( t ∑ γ i ( t ) ) ⋅ μ i = t = 1 ∑ T γ i ( t ) t = 1 ∑ T γ i ( t ) ⋅ x t
(这里 μ ^ i 的表达式与 Poisson-HMM 中 λ ^ i 相同)
∂ σ ^ i 2 ∂ L ⇒ i ∑ t ∑ γ i ( t ) ⋅ ( − 2 1 ⋅ σ ^ i 2 1 − 2 1 ( x t − μ i ) 2 ⋅ ( − 1 ) ⋅ ( σ ^ i 2 ) 2 1 ) ⇒ i ∑ t ∑ γ i ( t ) ⋅ ( σ ^ i 2 ( x t − μ i ) 2 − 1 ) ⇒ i ∑ t ∑ γ i ( t ) ⋅ [ ( x t − μ i ) 2 − σ ^ i 2 ] ⇒ σ ^ i 2 = 0 = 0 = 0 = 0 = t = 1 ∑ T γ i ( t ) t = 1 ∑ T γ i ( t ) ⋅ ( x t − μ ^ i ) 2
证毕
在最后,我们给出使用期望最大化算法估计泊松-隐马尔可夫链和正态-隐马尔可夫链参数的 Baum-Welch 公式。
泊松-隐马尔可夫链的 Baum-Welch 公式
π ^ i a ^ ij λ ^ i = γ i ( 1 ) = t = 1 ∑ T − 1 γ i ( t ) t = 1 ∑ T − 1 ε i , j ( t ) = g i h ij = t = 1 ∑ T γ i ( t ) t = 1 ∑ T γ i ( t ) ⋅ x t
正态-隐马尔可夫链的 Baum-Welch 公式
π ^ i a ^ ij μ ^ i σ ^ i 2 = γ i ( 1 ) = t = 1 ∑ T − 1 γ i ( t ) t = 1 ∑ T − 1 ε i , j ( t ) = g i h ij = t = 1 ∑ T γ i ( t ) t = 1 ∑ T γ i ( t ) ⋅ x t = t = 1 ∑ T γ i ( t ) t = 1 ∑ T γ i ( t ) ⋅ ( x t − μ ^ i ) 2
若估计的是正态分布的标准差,则调整为
σ ^ i = t = 1 ∑ T γ i ( t ) t = 1 ∑ T γ i ( t ) ⋅ ( x t − μ ^ i ) 2
维特比算法
现在我们来回答隐马尔可夫链的第二个基本问题,解码。所谓全局解码 (Global Decoding),就是找到最优的状态序列 s 1 , s 2 , … , s T ,其能最好的解释观测到的序列 x 1 , x 2 , … , x T 。从数学上讲,我们希望最大化
P ( C 1 T = c 1 T ∣ X 1 T = x 1 T )
等价的,也就是最大化
P ( C 1 T = c 1 T , X 1 T = x 1 T )
定义 12 [最优状态序列 (Optimal State Sequence)]:
若 s 1 , s 2 , … , s T 是最优状态序列,则
s 1 , s 2 , … , s T = arg c 1 , … , c T max P ( X 1 T = x 1 T , C 1 T = c 1 T )
寻找最优序列一般使用维特比算法 (Viterbi Algorithm),我们先给出 V j ( t ) 的定义及其递归公式
定义 13:
V i ( t ) V i ( 1 ) = c 1 , … , c t − 1 max P ( C 1 t − 1 = c 1 t − 1 , C t = i , X 1 t = x 1 t ) ( t = 2 , 3 , … , T ) = P ( C 1 = i , X 1 = x 1 ) = π i ⋅ p i ( x 1 )
定理 28** [V j ( t ) 的递归公式 (recursion)]:
V j ( t ) = ( i max V i ( t − 1 ) ⋅ a ij ) ⋅ p j ( x t ) ( t = 2 , 3 , … , T , i = 1 , 2 , … , m )
证明:
V j ( t ) = c 1 , … , c t − 1 max P ( C 1 t − 1 = c 1 t − 1 , C t = j , X 1 t = x 1 t ) = i max c 1 , … , c t − 2 max P ( C 1 t − 2 = c 1 t − 2 , C t − 1 = i , C t = j , X 1 t = x 1 t ) = i max c 1 , … , c t − 2 max P ( X t = x t , C t = j ∣ X 1 t − 1 = x 1 t − 1 , C 1 t − 2 = c 1 t − 2 , C t − 1 = i ) ⋅ P ( C 1 t − 2 = c 1 t − 2 , C t − 1 = i , X 1 t − 1 = x 1 t − 1 ) = i max c 1 , … , c t − 2 max P ( X t = x t , C t = j ∣ C t − 1 = i ) ⋅ P ( C 1 t − 2 = c 1 t − 2 , C t − 1 = i , X 1 t − 1 = x 1 t − 1 ) (取条件于 X 1 t − 1 = x 1 t − 1 , C 1 t − 2 = c 1 t − 2 , C t − 1 = i ) = i max c 1 , … , c t − 2 max P ( X t = x t ∣ C t = j ) ⋅ P ( C t = j ∣ C t − 1 = i ) ⋅ P ( C 1 t − 2 = c 1 t − 2 , C t − 1 = i , X 1 t − 1 = x 1 t − 1 ) (基于 H MM 性质 1 、 2 ) = i max [ p j ( x t ) ⋅ a ij ] ⋅ c 1 , … , c t − 2 max P ( C 1 t − 2 = c 1 t − 2 , C t − 1 = i , X 1 t − 1 = x 1 t − 1 ) (基于引理 1 ) = i max [ p j ( x t ) ⋅ a ij ] ⋅ V i ( t − 1 ) = [ i max V i ( t − 1 ) ⋅ a ij ] ⋅ p j ( x t ) = [ i max V i ( t − 1 ) ⋅ a ij ] ⋅ P ( X t = x t ∣ C t = j ) (基于 H MM 性质 2 ) = [ i max V i ( t − 1 ) ⋅ a ij ] ⋅ p j ( x t )
我们引入 V j ( t ) 的目的是给出如下计算最优序列的公式,即维特比公式。
定理 29 [维特比公式 (Viterbi Formula)]:
⎩ ⎨ ⎧ s T = arg max i V i ( T ) s t = arg max i ( V i ( t ) ⋅ a i , s t + 1 ) ( t = T − 1 , … , 2 , 1 )
证明:
(1) 找出 s T
c 1 , … , c T max P ( X 1 T = x 1 T , C 1 T = c 1 T ) = i max c 1 , … , c T − 1 max P ( X 1 T = x 1 T , C 1 T − 1 = c 1 T − 1 , C T = i ) = i max V i ( T )
因此,找出 s T 等同于找出 i 来最大化 V i ( T ) ,即
s T = arg i max V i ( T ) , 即 i max V i ( T ) = V s T ( T )
(2) 找出 s T − 1
V s T ( T ) = [ i max V i ( T − 1 ) ⋅ a i , s T ] ⋅ p s T ( x T ) (基于定理 28 )
注意到现在 s T 已经被找出,因此 p s T ( x T ) 是固定的。故找出 s T − 1 等同于找到 i 来最大化 V i ( T − 1 ) ⋅ a i , s T ,即
s T − 1 = arg i max [ V i ( T − 1 ) ⋅ a i , s T ]
找出 s T − 1 后
V s T ( T ) = [ V s T − 1 ( T − 1 ) ⋅ a s T − 1 , s T ] ⋅ p s T ( x T )
(3) 找出 s T − 2
V s T − 1 ( T − 1 ) = [ i max V i ( T − 2 ) ⋅ a i , s T − 1 ] ⋅ p s T − 1 ( x T − 1 ) (基于定理 28 )
现在 s T − 1 已经被找出,因此 p s T − 1 ( x T − 1 ) 是已知的。故找出 s T − 2 等同于找到 i 来最大化 V i ( T − 2 ) ⋅ a i , s T − 1 ,即
s T − 2 = arg i max [ V i ( T − 2 ) ⋅ a i , s T − 1 ]
找出 s T − 2 后
V s T − 1 ( T − 1 ) = [ V s T − 2 ( T − 2 ) ⋅ a s T − 2 , s T − 1 ] ⋅ p s T − 1 ( x T − 1 )
(4) 基于同样的逻辑,可以得出
s t = arg i max ( V i ( t ) ⋅ a i , s t + 1 ) ( t = T − 1 , … , 2 , 1 )
证毕