随机模型
详细内容可以参见 Github 博客
模型分为确定性模型和随机模型。而随机模型主要分为以下三类
- 概率模型
- 马氏链模型
- 回归模型(单独总结)
概率模型
主要用到概率的运算,概率分布期望方差等基本知识。以轧刀中的浪费为例。
轧钢中的浪费
问题为:已知成品材的规定长度为 l l l 以及粗轧后的钢材均方差 σ \sigma σ ,确定粗轧后的钢材长度均值为 m m m ,使得当轧机以 m m m 粗轧,再进行精轧得到成品时浪费最少。
最简单的思路便是把两部分的损失加起来:
W = ∫ l ∞ ( x − l ) p ( x ) d x + ∫ − ∞ l x p ( x ) d x W=\int_l^{\infty}(x-l)p(x)dx+\int_{-\infty}^lxp(x)dx W=∫l∞(x−l)p(x)dx+∫−∞lxp(x)dx
利用概率密度函数的性质(积分为1),期望定义,上式可改为:
(*) W = m − l ∫ l ∞ p ( x ) d x = m − l P W=m-l\int_l^{\infty}p(x)dx =m-lP \tag{*} W=m−l∫l∞p(x)dx=m−lP(*)
但事实这并不合理,应该考虑没得到一根钢材平均浪费了多少而不是每轧一个浪费了多少。
粗轧 N N N 根浪费了 m N − l P N mN-lPN mN−lPN 长度,但是只有 P N PN PN 成了成品,所以
J = m N − l P N P N = m P − l J=\frac{mN-lPN}{PN}=\frac{m}{P}-l J=PNmN−lPN=Pm−l
最小化 J J J 即可,其中
P ( m ) = ∫ l ∞ p ( x ) d x p ( x ) = 1 2 π σ e − ( x − m ) 2 2 σ 2 P(m)=\int_l^{\infty}p(x)dx \\ p(x)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-m)^2}{2\sigma^2}} P(m)=∫l∞p(x)dxp(x)=2πσ1e−2σ2(x−m)2
马氏链模型
在一个离散时间集合 T = { 0 , 1 , 2 , ⋯   } T=\{0,1,2,\cdots\} T={0,1,2,⋯} 和一个有限或可列无穷的状态空间 S = { 1 , 2 , ⋯   } S=\{1,2,\cdots\} S={1,2,⋯}上,一个随机过程在任一时刻从一个状态以一定的概率向其他状态转移(或保持原状态不变)。记 X n X_n Xn 为时刻 n n n 时时刻过程所处的状态,假定:
- 在时刻0, 过程所处的状态 X 0 X_0 X0 是 S S S 上的一个随机变量;
- 在任一时刻 n n n,给定 X 0 , X 1 , ⋯ X n X_0,X_1,\cdots X_n X0,X1,⋯Xn时, X n + 1 X_{n+1} Xn+1 的条件分布只与 X n X_n Xn 有关,而与 X 0 , X 1 , ⋯ X n − 1 X_0,X_1,\cdots X_{n-1} X0,X1,⋯Xn−1无关。
满足上述条件的随机过程为马尔可夫链,简称马氏链。通常有以下两种类型:
正则链
从任意状态出发都可以经有限次转移到达任意状态
马氏链是正则链的充要条件是:存在正整数 N N N ,使得 P N > 0 \boldsymbol{P}^N>0 PN>0 。即每个元素都大于零
正则链存在唯一的极限状态概率 w = ( w 1 , w 2 , ⋯   , w n ) \boldsymbol{w}=(w_1,w_2,\cdots,w_n) w=(w1,w2,⋯,wn) ,使得 a n → w a_n \rightarrow \boldsymbol{w} an→w 。稳态状态与初始状态无关,满足 w P = w , ∑ w i = 1 \boldsymbol{w}\boldsymbol{P}=\boldsymbol{w},\sum w_i=1 wP=w,∑wi=1
从状态
i
i
i 经
n
n
n 次转移第一次到
j
j
j 的概率称为首达概率,记作
f
i
j
(
n
)
f_{ij}(n)
fij(n) ,则转移的平均次数:
μ
i
j
=
∑
i
=
1
∞
n
f
i
j
(
n
)
\mu_{ij}=\sum_{i=1}^{\infty}nf_{ij}(n)
μij=i=1∑∞nfij(n)
特别的,针对正则链, μ i i = 1 / w i \mu_{ii}=1/w_i μii=1/wi
吸收链
转移概率 p i i = 1 p_{ii}=1 pii=1 的状态 i i i 称为吸收状态。如果某一个马氏链包含只招一个吸收状态并且从非吸收状态可经有限次转移达到吸收状态,则称之为吸收链。
吸收链的转移矩阵可以化为简单的标准形式( r r r 个吸收态全部提前)
P = [ I r × r O R Q ] P = \left[ {\begin{array}{c} {{I_{r \times r}}}&O\\ R&Q \end{array}} \right] P=[Ir×rROQ]
上述中 I − Q I-Q I−Q 可逆。且逆阵 M = ( I − Q ) − 1 = ∑ s = 0 ∞ Q s M=(I-Q)^{-1}=\sum \limits_{s=0}^{\infty}Q^s M=(I−Q)−1=s=0∑∞Qs 。记元素全为 1 的向量 e = ( 1 , 1 , ⋯   , 1 ) e=(1,1,\cdots,1) e=(1,1,⋯,1) 。则 y = M e y=Me y=Me 的第 i i i 个分量是第 i i i 个非吸收状态出发,被某个吸收状态吸收的平均次数。首达概率矩阵 F = M R F=MR F=MR
隐马尔科夫模型 HMM
网络结构
HMM (Hidden Markov Model)是关于时序的概率模型,描述由一个隐藏的马尔科夫链生成不可观测的状态随机序列,再由各个状态生成观测随机序列的过程.
就是马尔科夫链本身状态不可观测,但是每个不可观测状态 z z z 又会生成一个可观测状态 x x x 。
定义
参数:
- Q = { q 1 , q 2 , ⋯   , q N } Q=\{q_1,q_2,\cdots,q_N\} Q={q1,q2,⋯,qN} 为可能的状态集合 , N N N 是可能的状态数
- V = { v 1 , v 2 , ⋯   , v M } V=\{v_1,v_2,\cdots,v_M\} V={v1,v2,⋯,vM} 为可能的观测集合, M M M 是观测集合的大小
- I = { i 1 , i 2 , ⋯   , i T } I=\{i_1,i_2,\cdots,i_T\} I={i1,i2,⋯,iT} 为长度为 T T T 的状态序列(不可测), i i ∈ Q i_i \in Q ii∈Q
- O = { o 1 , o 2 , ⋯   , o T } O=\{o_1,o_2,\cdots,o_T\} O={o1,o2,⋯,oT} 为对应的观测序列。 o i ∈ V o_i \in V oi∈V
HMM的三要素
HMM由初始概率分布 π \pi π (向量)、状态转移概率分布 A A A (矩阵) 以及观测概率分布 B B B (矩阵) 确定. π \pi π 和 A A A 决定状态序列, B B B 决定观测序列。因此, HMM可以用三元符号表示, 称为HMM的三要素:
λ = ( A , B , π ) \lambda=(A,B,\pi) λ=(A,B,π)
A A A 是转移概率矩阵, A = [ a i j ] N × N A=[a_{ij}]_{N\times N} A=[aij]N×N ,其中 a i j a_{ij} aij 是时刻 t t t 处于状态 q i q_i qi 的条件下时刻 t + 1 t+1 t+1 转移至 q j q_j qj 的概率
a i j = P ( i t + 1 = q j ∣ i t = q i ) , i = 1 , 2 , ⋯   , N ; j = 1 , 2 , ⋯   , N a_{ij}=P(i_{t+1}=q_j\mid i_t=q_i), \quad i=1,2,\cdots,N;j=1,2,\cdots,N aij=P(it+1=qj∣it=qi),i=1,2,⋯,N;j=1,2,⋯,N
B B B 是观测转移概率矩阵: 其中, b j ( k ) b_j(k) bj(k) 是在时刻 t t t 处于状态 q j q_j qj 的条件下生成观测 v k v_k vk的概率:
B = [ b j ( k ) ] N × M b j ( k ) = P ( o t = v k ∣ i t = q j ) , k = 1 , 2 , ⋯   , M ; j = 1 , 2 , ⋯   , N B=[b_{j}(k)]_{N\times M} \\ b_j(k)=P(o_t=v_k\mid i_t=q_j),\quad k=1,2,\cdots,M;j=1,2,\cdots,N B=[bj(k)]N×Mbj(k)=P(ot=vk∣it=qj),k=1,2,⋯,M;j=1,2,⋯,N
π \pi π 是初始状态概率向量: π = ( π ) N × 1 \pi=(\pi)_{N \times1} π=(π)N×1 。
π i = P ( i 1 = q i ) , i = 1 , 2 , ⋯   , N \pi_i=P(i_1=q_i),\quad i=1,2,\cdots,N πi=P(i1=qi),i=1,2,⋯,N
一个HMM模型仅有以上三项确定,而与观测序列无关,因为本身就是随机的
HMM的两个假设
齐次马尔可夫性假设 :任意时刻 t t t 的状态, 只依赖于其前一刻的状态, 与其他时刻的状态及观测无关, 也与时刻 t t t 无关.
P ( i t ∣ i t − 1 , o t − 1 , ⋯   , i 1 , o 1 ) = P ( i t ∣ i t − 1 ) , t = 1 , 2 , ⋯   , T P(i_t\mid i_{t-1},o_{t-1},\cdots,i_1,o_1)=P(i_t\mid i_{t-1}), \quad t=1,2,\cdots,T P(it∣it−1,ot−1,⋯,i1,o1)=P(it∣it−1),t=1,2,⋯,T
观测独立性假设 : 任何时刻的观测只依赖于该时刻的马尔科夫链状态. 与其他观测及状态无关.
P ( o t ∣ i T , o T , ⋯   , i t , o t , ⋯   , i 1 , o 1 ) = P ( o t ∣ i t ) , t = 1 , 2 , ⋯   , T P(o_t \mid i_{T},o_{T},\cdots,i_{t},o_{t},\cdots,i_1,o_1)=P(o_t\mid i_{t}), \quad t=1,2,\cdots,T P(ot∣iT,oT,⋯,it,ot,⋯,i1,o1)=P(ot∣it),t=1,2,⋯,T
HMM的三个基本问题
实现可参考:隐马尔科夫模型(HMM)及其Python实现
数学原理参考书籍:《统计学习方法》
概率计算问题——前向后向算法(动态规划)
给定模型 λ = ( A , B , π ) \lambda=(A,B,\pi) λ=(A,B,π) 和观测序列 O = ( o 1 , o 2 , ⋯   , o T ) O=(o_1,o_2,\cdots,o_T) O=(o1,o2,⋯,oT) ,计算模型 λ \lambda λ 下观测序列 O O O 出现的概率 P ( O ∣ λ ) P(O \mid\lambda) P(O∣λ)
显然给定模型可以直接计算,即对所有的可能状态序列下的该 观测序列求和
P ( O ∣ λ ) = ∑ I P ( O ∣ I , λ ) P ( I ∣ λ ) = ∑ i 1 , i 2 , ⋯   , i T ⎵ N × N × ⋯ × N π i 1 a i 1 i 2 b i 1 ( o 1 ) ⋯ a i T − 1 i T b i T ( o T ) P(O\mid\lambda)=\sum_IP(O \mid I,\lambda)P(I \mid \lambda) \\ =\sum\limits_{\underbrace {i_1,i_2, \cdots ,i_T}_{N \times N \times \cdots \times N}} {\pi_{i_1}a_{i_1i_2}b_{i_1}(o_1)\cdots a_{i_{T-1}i_T}b_{i_T}(o_T)} P(O∣λ)=I∑P(O∣I,λ)P(I∣λ)=N×N×⋯×N i1,i2,⋯,iT∑πi1ai1i2bi1(o1)⋯aiT−1iTbiT(oT)
但是上述算法的复杂度达到了 O ( T N T ) O(TN^T) O(TNT) ,不可取。
前向概率:给定隐马尔科夫模型 λ \lambda λ ,定义为到时刻 t t t 部分观测序列为 o 1 , o 2 , ⋯   , o t o_1,o_2,\cdots,o_t o1,o2,⋯,ot 且状态为 q i q_i qi 的概率
α t ( i ) = P ( o 1 , o 2 , ⋯   , o t ∣ i t = q i , λ ) \alpha_t(i)=P(o_1,o_2,\cdots,o_t \mid i_t=q_i,\lambda) αt(i)=P(o1,o2,⋯,ot∣it=qi,λ)后向概率:给定隐马尔科夫模型 λ \lambda λ ,定义为到时刻 t t t 状态为 q i q_i qi 的条件下,从 t + 1 t+1 t+1 到 T T T 的部分观测序列为 o t + 1 , o t + 2 , ⋯   , o T o_{t+1},o_{t+2},\cdots,o_T ot+1,ot+2,⋯,oT 的概率。
β t ( i ) = P ( o t + 1 , o t + 2 , ⋯   , o T ∣ i t = q i , λ ) \beta_t(i)=P(o_{t+1},o_{t+2},\cdots,o_T \mid i_t=q_i,\lambda) βt(i)=P(ot+1,ot+2,⋯,oT∣it=qi,λ)
本处给出前向算法的步骤:(递归实现)
-
给定初值: α 1 ( i ) = π i b i ( o 1 ) , i = 1 , 2 , ⋯   , N \alpha_1(i)=\pi_ib_i(o_1),i=1,2,\cdots,N α1(i)=πibi(o1),i=1,2,⋯,N
-
递推:对于 t = 1 , 2 , ⋯   , T − 1 t=1,2,\cdots,T-1 t=1,2,⋯,T−1
α t + 1 ( i ) = [ ∑ j = 1 N α t ( j ) a j i ] b i ( o t + 1 ) \alpha_{t+1}(i)=[\sum_{j=1}^N\alpha_t(j)a_{ji}]b_i(o_{t+1}) αt+1(i)=[j=1∑Nαt(j)aji]bi(ot+1)
- 最终: P ( O ∣ λ ) = ∑ i = 1 N α T ( i ) P(O \mid \lambda)=\sum \limits_{i=1}^N\alpha_T(i) P(O∣λ)=i=1∑NαT(i)
显然通过前向概率的引入作为中间变量从而实现了递归,大大降低了算法复杂度
学习问题——Baum-Welch算法
已知观测序列 O = o 1 , o 2 , ⋯   , o T O=o_1,o_2,\cdots,o_T O=o1,o2,⋯,oT 。估计模型 λ = ( A , B , π ) \lambda=(A,B,\pi) λ=(A,B,π) ,使得 P ( O ∣ λ ) P(O \mid \lambda) P(O∣λ) 最大,即用极大似然法估计参数
有监督的学习算法:极大似然法
假设有监督意味着还知道每一个状态序列对应的状态序列 I = i 1 , i 2 , ⋯   , i T I=i_1,i_2,\cdots,i_T I=i1,i2,⋯,iT 。
相当于给定了完全数据 { ( o 1 , i 1 ) , ( o 2 , i 2 ) , ⋯   , ( o T , i T ) } \{(o_1,i_1),(o_2,i_2),\cdots,(o_T,i_T)\} {(o1,i1),(o2,i2),⋯,(oT,iT)} ,具体估计如下:
转移概率的估计:
设样本时刻 t t t 处于状态 i i i 的,时刻 t + 1 t+1 t+1 处于状态 j j j 的频数为 A i j A_{ij} Aij ,那么
α ^ i j = A i j ∑ j = 1 N A i j , i = 1 , 2 , ⋯   , N ; j = 1 , 2 , ⋯   , N \widehat{\alpha}_{ij}=\frac{A_{ij}}{\sum \limits_{j=1}^NA_{ij}} ,\quad i=1,2,\cdots,N;j=1,2,\cdots,N α ij=j=1∑NAijAij,i=1,2,⋯,N;j=1,2,⋯,N
似乎只对到达的末端进行求和算频率
观测概率的估计:
设样本中状态为 j j j 并观测为 k k k 的频数为 $B_{jk} $ ,那么
b ^ j ( k ) = B j k ∑ k = 1 M B j k , j = 1 , 2 , ⋯   , N ; k = 1 , 2 , ⋯   , M \widehat{b}_{j}(k)=\frac{B_{jk}}{\sum \limits_{k=1}^MB_{jk}},\quad j=1,2,\cdots,N;k=1,2,\cdots,M b j(k)=k=1∑MBjkBjk,j=1,2,⋯,N;k=1,2,⋯,M
初始状态的估计: π i \pi_i πi 即样本中初始状态为 q i q_i qi 的频率
从中便可以看到极大似然估计的美,基本与样本的均值形式等同,如样本服从指数分布,则参数 μ \mu μ 便是样本均值。
然而监督学习代价往往比较高,有时便会利用以下无监督学习方法。
无监督的学习算法:Baum-Weich算法
此时我们仅仅有一堆观测数据,这种情况下我们可以使用 EM 算法,将状态变量视为隐变量。使用 EM 算法学习HMM的参数的算法称为 BW 算法
步骤如下:
确定完全数据的对数似然函数 log P ( O , I ∣ λ ) \log P(O,I\mid \lambda) logP(O,I∣λ)
EM 算法的 E 步:计算 Q Q Q 函数
Q Q Q 函数是完全数据的对数似然函数在给定的观测 O O O 与当前参数 λ ^ \widehat{\lambda} λ 下对未观测数据 I I I 的条件概率分布的期望,即
Q ( λ , λ ^ ) = E I [ log P ( O , I ∣ λ ) ∣ O , λ ^ ] = ∑ I log P ( O , I ∣ λ ) P ( I ∣ O , λ ^ ) P ( O , I ∣ λ ) = π i 1 a i 1 i 2 b i 1 ( o 1 ) ⋯ a i T − 1 i T b i T ( o T ) Q(\lambda,\widehat{\lambda})=E_I[\log P(O,I\mid\lambda)\mid O,\widehat{\lambda}]=\sum_{I}\log P(O,I\mid \lambda)P(I\mid O,\widehat{\lambda}) \\ P(O,I\mid\lambda)=\pi_{i_1}a_{i_1i_2}b_{i_1}(o_1)\cdots a_{i_{T-1}i_T}b_{i_T}(o_T) Q(λ,λ )=EI[logP(O,I∣λ)∣O,λ ]=I∑logP(O,I∣λ)P(I∣O,λ )P(O,I∣λ)=πi1ai1i2bi1(o1)⋯aiT−1iTbiT(oT)
可以理解为计算给定初始参数,针对每一种可能的状态序列 I I I 下对数似然函数的期望,最大化后的 λ \lambda λ 优于 λ ^ \widehat{\lambda} λ
EM 算法的 M 步:最大化上述 Q Q Q 函数估计 λ \lambda λ 作为下一次迭代初始参数,具体推倒见《统计学习方法》,直接给出参数估计式
(*) α i j = ∑ t = 1 T − 1 ξ t ( i , j ) ∑ t = 1 T − 1 γ t ( i ) b j ( k ) = ∑ t = 1 , o t = v k T γ t ( j ) ∑ t = 1 T γ t ( j ) π i = γ 1 ( i ) \alpha_{ij}=\frac{\sum \limits_{t=1}^{T-1}\xi_t(i,j)}{\sum \limits_{t=1}^{T-1}\gamma_t(i)} \\ b_j(k)=\frac{\sum \limits_{t=1,o_t=v_k}^{T}\gamma_t(j)}{\sum \limits_{t=1}^{T}\gamma_t(j)} \tag{*} \\ \pi_i=\gamma_1(i) αij=t=1∑T−1γt(i)t=1∑T−1ξt(i,j)bj(k)=t=1∑Tγt(j)t=1,ot=vk∑Tγt(j)πi=γ1(i)(*)
其中
γ t ( i ) = P ( i t = q i ∣ O , λ ) = P ( i t = q i , O ∣ λ ) P ( O ∣ λ ) = α t ( i ) β t ( i ) ∑ k = 1 N α t ( j ) β t ( j ) \gamma_t(i)=P(i_t=q_i\mid O,\lambda) \\ =\frac{P(i_t=q_i,O\mid \lambda)}{P(O \mid \lambda)} \\ =\frac{\alpha_t(i)\beta_t(i)}{\sum \limits_{k=1}^N\alpha_t(j)\beta_t(j)} γt(i)=P(it=qi∣O,λ)=P(O∣λ)P(it=qi,O∣λ)=k=1∑Nαt(j)βt(j)αt(i)βt(i)
表示给定模型参数和所有观测,时刻 t t t 处于状态 q i q_i qi 的概率。 α , β \alpha,\beta α,β 分别为前向和后向概率。
ξ t ( i , j ) = α t ( i ) a i j b j ( o t + 1 ) β t + 1 ( j ) ∑ i = 1 N ∑ j = 1 N α t ( i ) a i j b j ( o t + 1 ) β t + 1 ( j ) \xi_t(i,j)=\frac{\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}{\sum \limits_{i=1}^N \sum\limits_{j=1}^N\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)} ξt(i,j)=i=1∑Nj=1∑Nαt(i)aijbj(ot+1)βt+1(j)αt(i)aijbj(ot+1)βt+1(j)
表示给定模型和所有观测,时刻 t t t 处于状态 q i q_i qi ,而时刻 t + 1 t+1 t+1 处于状态 q j q_j qj 的概率。
下面给出算法流程图
输入:观测数据 O = ( o 1 , o 2 , ⋯   , o T ) O=(o_1,o_2,\cdots,o_T) O=(o1,o2,⋯,oT);
输出:隐马尔科夫模型参数
- 参数初始化:对 n = 0 n=0 n=0 ,选取 α i j ( 0 ) , b j ( k ) ( 0 ) , π i ( 0 ) \alpha_{ij}^{(0)},b_j(k)^{(0)},\pi_i^{(0)} αij(0),bj(k)(0),πi(0) 得到模型 λ ( 0 ) = ( A ( 0 ) , B ( 0 ) , π ( 0 ) ) \lambda^{(0)}=(A^{(0)},B^{(0)},\pi^{(0)}) λ(0)=(A(0),B(0),π(0))
- 递推,对 $n=1,2,\cdots $ ,按(*)式求得下一次迭代模型参数,即 ( * )式右边方程式由 λ ( i ) \lambda^{(i)} λ(i) 计算,得到 λ ( i + 1 ) \lambda^{(i+1)} λ(i+1) 的模型参数。
- 终止
∣ ∣ λ ( i + 1 ) − λ ( i ) ∣ ∣ < ε ||\lambda^{(i+1)}-\lambda^{(i)}||<\varepsilon ∣∣λ(i+1)−λ(i)∣∣<ε
得到模型参数 λ = λ ( n + 1 ) \lambda=\lambda^{(n+1)} λ=λ(n+1)
论文中绘制算法流程图时上述公式都要写出来
预测算法——维特比(Viterbi)算法
也称为解码问题。已知观测序列 O = o 1 , o 2 , ⋯   , o T O=o_1,o_2,\cdots,o_T O=o1,o2,⋯,oT 和模型 λ = ( A , B , π ) \lambda=(A,B,\pi) λ=(A,B,π) ,求给定观测序列使得条件概率 P ( I ∣ O ) P(I \mid O) P(I∣O) 最大的状态序列 I I I 。相当于预测出一个观测序列对应的状态向量
维特比算法实际是用动态规划解隐马尔科夫模型预测问题,即用动态规划求概率最大路径(最优路径),这是一个状态路径对应一个状态序列。