这部分要讲解马尔科夫链排在前三的应用。
Google \text{Google} Google
PageRank是谷歌的核心技术。用于实现搜索引擎得到的网页重要性排序。
Web Pages → Page Ranking → Searching → Searching Engines \text{Web Pages} \rightarrow \text{Page Ranking} \rightarrow \text{Searching} \rightarrow \text{Searching Engines} Web Pages→Page Ranking→Searching→Searching Engines
PageRank是1998-2000年由Brein和Page做出的工作
这可以看做是一个有向图求解极限转移概率的问题。之前讨论的是无向图的极限概率转移。
对于无向图，只有度的概念，而有向图分成了入度和出度。
之所以用有向图模型是因为，我们可以把网页看做是图的结点，把节点之间的连接，看做是网页之间的连接
Web Pages → Nodes of Graph Link → Edges of Graph \text{Web Pages} \rightarrow \text{Nodes of Graph} \\ \text{Link} \rightarrow \text{Edges of Graph} Web Pages→Nodes of GraphLink→Edges of Graph
一开始，人们考虑用入度来表示结点的重要性。因为如果一个网页很重要，必然会有很多的网站引用他。
Importance ⇔ in degree \text{Importance} \Leftrightarrow \text{in degree} Importance⇔in degree
但是只用入度来描述是不够的，因为这个有向图模型中C和D的入度都是2，但是可以看出来，D更加的重要。
Brin和Page提出了自己的观点，认为网页节点的重要性应该等同于极限分布。当经过n步转移之后，各个节点被访问的情况就趋于平稳了，就能够判断哪个节点更加的重要了，就是访问节点次数除以总的访问次数。
Importance ⇔ Limit Distribution \text{Importance} \Leftrightarrow \text{Limit Distribution} Importance⇔Limit Distribution
但是，求这个模型的极限分布是有问题的，因为这个概率转移矩阵很稀疏，只有个别点是有值的，大部分点是没有值的。因此会存在分布极限存在性问题，只有不可约+非周期才能够保证分布极限的存在。这个稀疏的矩阵这两点都难以保证。
因此，Brin和Pgae的策略是，在原来的一步转移概率矩阵上做了调整，增加上一个全一矩阵，得到了一个新的一步转移概率矩阵。
P ~ = ( 1 − α ) P + α n ( 1 . . . 1 . . . . . . . . . 1 . . . 1 ) Google Matrix \widetilde{P} = (1- \alpha)P + \frac{\alpha}{n} \begin{pmatrix} 1 & ...&1 \\ ... & ...&... \\ 1&...& 1 \end{pmatrix} \\ \text{Google Matrix} P =(1−α)P+nα⎝⎛1...1.........1...1⎠⎞Google Matrix
这个模型可以从这些角度进行分析
π = π P ~ \pi = \pi \widetilde{P} π=πP
即求矩阵P的左特征值。
Principal Left Eigenvector of P ~ ⇒ Q R \text{Principal Left Eigenvector of } \widetilde{P} \Rightarrow QR Principal Left Eigenvector of P ⇒QR
求特征值可以使用QR分解的方法
MCMC是离散马尔科夫蒙特卡洛的缩写
MCMC ⇒ Markov Chain Monte Carlo \text{MCMC} \Rightarrow \text{Markov Chain Monte Carlo} MCMC⇒Markov Chain Monte Carlo
蒙特卡洛技术是指，我们想要做某种仿真，仿真中有很多的随机效应，这样的随机效应一定是由某种随机变量描述的。
这个随机变量的分布可以很简单，比如用高斯分布、均匀分布等模型来描述，这个分布也可能很复杂，可能是不连续的、不光滑的，即使知道了分布，也很那求均值。
因为，我们很难通过分布来描述一个随机变量，最终就只能通过仿真而不是解析的方式来描述这些样本
Simulation → Random Effect → Random Variables → Distribution → Simulation \text{Simulation} \rightarrow \text{Random Effect} \rightarrow \text{Random Variables} \rightarrow \text{Distribution} \rightarrow \text{Simulation} Simulation→Random Effect→Random Variables→Distribution→Simulation
所谓仿真，就是产生一组伪随机样本,这些伪随机样本可以产生某种分布
G e n e r a t e Z 1 , . . . , Z n ∼ f ( x ) Pseudo-Random Sample Generate \quad Z_{1},...,Z_{n} \sim f(x) \\ \text{Pseudo-Random Sample} GenerateZ1,...,Zn∼f(x)Pseudo-Random Sample
我们对这些分布并不了解，但是我们可以转化为对样本的研究。
比如我们需要计算期望，但是积分非常难算，我们就可以用样本均值替代。
E ( Z ) = ∫ − ∞ + ∞ x f ( x ) d x → 1 n ∑ k = 1 n Z k E(Z) = \int_{-\infty}^{+\infty} xf(x)dx \\ \rightarrow \frac{1}{n} \sum_{k=1}^n Z_k E(Z)=∫−∞+∞xf(x)dx→n1k=1∑nZk
甚至我们可以通过直方图的方法，估计分布情况。
如果我们有足够多的样本，甚至就不需要分布了，因为我们对分布没有办法进行解析计算。
因此，蒙特卡洛模拟在随机变量仿真上，具有非常重要的地位。蒙特卡洛模拟由下面的人发明
蒙特卡洛的核心思想，就是把对分布的解析研究，转换为对样本的研究。
具体目的是，如何做到，随便给一个分布，就能够得到服从这个分布的样本。
蒙特卡洛的困难在于，有了样本是可以得到对分布的认知的，但是现在要做的是逆向操作，有了分布如何去获取样本
Z 1 , . . . , Z n → f ( x ) f ( x ) ? → Z 1 , . . , Z n Z_1,...,Z_n \rightarrow f(x) \\ f(x) ?\rightarrow Z_1,..,Z_n Z1,...,Zn→f(x)f(x)?→Z1,..,Zn
同时，我们希望这种方法是通用的，一般来说，对某种具体的分布的采样策略是已知的，但是如果没有通用的策略，如果出现了一种新的分布，就没有办法进行仿真。
Universal Pseudo-Random Sampler \text{Universal Pseudo-Random Sampler} Universal Pseudo-Random Sampler
我们希望有一个通用的伪随机采样方法，能够得到任意分布的样本。
马氏链是对这个问题最好的解决方案。P是一步转移概率矩阵。bi是初始分布。我们假设是不可约的，非周期的。
P = ( P i j ) i j b i = P ( Z 0 = i ) Irreducible + Non-Perdic ⇒ π = π P P = (P_{ij})_{ij} \\ b_i = P(Z_0 =i) \\ \text{Irreducible + Non-Perdic} \Rightarrow \pi = \pi P P=(Pij)ijbi=P(Z0=i)Irreducible + Non-Perdic⇒π=πP
现在我们假设极限分布是已知的，就是我们的目标分布。虽然目标分布可能是连续的，不过不本质，连续空间马氏链也有求解极限分布的方法。
我们假设极限分布就是我们已知的需要采样的分布。然后我们反过来求一步转移概率矩阵。
这是个非适定问题。因为P有n2个元素，π有n个元素，也就是用n个方程求解n2个未知数
假如是能够求解到一步转移概率的话，就能产生随机样本了。因为马氏链的转移规律已经知道了，我们只需要让马氏链运行起来即可。然后抛弃前面的样本，因为前面的马氏链还没有收敛。然后每次run一步都能够得到一个样本轨道上的点，每一步都是服从样本分布的。
这是迄今为止，最好的，通用的产生随机样本的方法。
需要分三步来求解。
Detailed Balance \text{Detailed Balance} Detailed Balance
我们假设有极限分布π，还有一步转移概率P，我们要求解的方程如下
π = ( π 0 , π 1 , . . . , π n ) P = ( P i j ) i j π = π P \pi = (\pi_0,\pi_1,...,\pi_n) \quad P = (P_{ij})_{ij} \\ \pi = \pi P π=(π0,π1,...,πn)P=(Pij)ijπ=πP
但是这个方程不好验证，我们需要找一个充分条件，当满足充分条件的时候，这个方程成立。这个条件就叫做细致平衡
π i P i j = π j P j i ⇒ π = π P \pi_i P_{ij} = \pi_j P_{ji} \Rightarrow \pi = \pi P πiPij=πjPji⇒π=πP
细致指的是，具体到两个状态之间的转移关系。
平衡指的是，我们假设π是不同结点上的某种物质的量，最终平衡就反应了不同结点上有不同物质的量，那这个方程直观解释就是，从i转移到j的量和从j转移到i的量是一样的，转移量是一样的，就是平衡了。
我们来验证这个事情，我们对第j个方程用细致平衡，来看等式两边是否成立,πj就是π与P第j列的乘法
π j = ( π P ) j = ∑ i π i P i j = ∑ i π j P j i = π j ∑ i P j i = π j ∗ 1 = π j \pi_j=(\pi P)_j = \sum_i \pi_i P_{ij} = \sum_i \pi_j P_{ji} \\ =\pi_j \sum_i P_{ji} = \pi_j *1 = \pi_j πj=(πP)j=i∑πiPij=i∑πjPji=πji∑Pji=πj∗1=πj
因此细致平衡成立，并且是方程的充分条件，细致平衡成立，方程就成立。
其实MCMC的一步转移矩阵P开始是任取的,但是一般不满足细致平衡条件
∀ P ( P r o p o s a l ) π i P i j = π j P j i \forall P(Proposal) \\ \pi_i P_{ij} \cancel= \pi_j P_{ji} ∀P(Proposal)πiPij= πjPji
不过对P稍加调整，就能够满足细致平衡条件
P ~ i j = P i j m i n ( 1 , π j P j i π i P i j ) π i P i j ~ = π j P j i ~ \widetilde{P}_{ij} = P_{ij} min(1,\frac{\pi_j P_{ji}}{\pi_i P_{ij}}) \\ \pi_i \widetilde{P_{ij}} = \pi_j \widetilde{P_{ji}} P ij=Pijmin(1,πiPijπjPji)πiPij =πjPji
我们进行证明
π i P ~ i j = { π j P j i π i P i j ≥ π j P j i π i P i j π i P i j < π j P j i \pi_i \widetilde{P}_{ij} = \begin{cases} \pi_j P_{ji} & \pi_i P_{ij} \geq \pi_j P_{ji}\\ \pi_iP_{ij} & \pi_i P_{ij} < \pi_j P_{ji} \end{cases} πiP ij={πjPjiπiPijπiPij≥πjPjiπiPij<πjPji
π j P ~ j i = { π j P j i π i P i j ≥ π j P j i π i P i j π i P i j < π j P j i \pi_j \widetilde{P}_{ji} = \begin{cases} \pi_j P_{ji} & \pi_i P_{ij} \geq \pi_j P_{ji}\\ \pi_iP_{ij} & \pi_i P_{ij} < \pi_j P_{ji} \end{cases} πjP ji={πjPjiπiPijπiPij≥πjPjiπiPij<πjPji
则
π i P i j ~ = π j P j i ~ \pi_i \widetilde{P_{ij}} = \pi_j \widetilde{P_{ji}} πiPij =πjPji
这是metropolis和hastings的工作
Metropolis ( 1953 ) Hastings ( 1970 ) \text{Metropolis}(1953) \\ \text{Hastings}(1970) Metropolis(1953)Hastings(1970)
一步转移矩阵经过无穷多次迭代收敛就是我们需要的目标矩阵，收敛以后得到的样本轨道，就是我们需要的样本。但是收敛之前的样本要丢弃。
后来人们发现P矩阵有很多，需要好好挑选，因为不同矩阵的收敛次数不一样，可能有些200次就达到了收敛，有些要20000次才能达到收敛。
开始MCMC没有得到广泛应用，因为计算量比较大，人们觉得只是小样本量的模拟，但是后来随机计算机性能的提高，MCMC成为了通用性的模拟随机分布的重要工具。
先解释一下什么是隐马尔科夫
Hidden Markov Model \text{Hidden Markov Model} Hidden Markov Model
隐马尔科夫模型是对马尔科夫模型的一种扩展。隐马尔科夫包括了状态转移模型和状态观测模型两部分，状态转移模型用来描述内部规律的变化，但是很多时候运动规律是不可获取的，就会需要有一个观测模型，对状态进行观测。
Sequential ⇒ Markov \text{Sequential} \Rightarrow \text{Markov} Sequential⇒Markov
隐马尔科夫有两个转移方程，一个是状态间的转移，一个是状态到观测的转移
States P i j = P ( S n = j ∣ S n − 1 = i ) Observations B i j ( n ) = P ( O n = j ∣ S n = i ) \text{States} \\ P_{ij} = P(S_n = j|S_{n-1} = i) \\ \text{Observations} B_{ij}(n) = P(O_n=j|S_n = i) StatesPij=P(Sn=j∣Sn−1=i)ObservationsBij(n)=P(On=j∣Sn=i)
从状态转移关系可以看出，我们已经做了平稳性假定，因为这个概率与时刻n没有关系
隐马尔科夫典型应用就是语音识别
Speech Recognition \text{Speech Recognition} Speech Recognition
具体原理是，每个字的发音是一种系统响应，然后我们会对声音进行观测。隐马尔科夫可以把两个效应都描述出来。一个是发音1和2的区别，一个是每个人发音1的区别。
隐马尔科夫模型能够同时描述内在变化规律和外在表现规律。
我们现在假设，我们已经通过各种训练，得到语音识别需要的每个不同发音的状态转移规律和状态观测规律，下面我们要做两种计算。
我们假设现在有一些观测到的数据,我们可以通过分析这组观测数据在哪一个发音的马氏链中出现的概率大，来判断这是哪个音。
O 1 , O 2 , . . . , O n O_1,O_2,...,O_n O1,O2,...,On
P ( O 1 , . . . , O n ) P(O_1,...,O_n) P(O1,...,On)
用全概率公式展开
P ( O 1 , . . . , O n ) = ∑ S 1 , . . , S n P ( O 1 , . . . , O n ∣ S 1 , . . . , S n ) P ( S 1 , . . . , S n ) = ∑ S 1 , . . , S n ( ∏ k = 1 n P ( O k ∣ S k ) ) ∏ k = 1 n P ( S k ∣ S k − 1 ) P(O_1,...,O_n) = \sum_{S_1,..,S_n} P(O_1,...,O_n|S_1,...,S_n)P(S_1,...,S_n) \\ = \sum_{S_1,..,S_n} (\prod_{k=1}^n P(O_k|S_k)) \prod_{k=1}^n P(S_k|S_{k-1}) P(O1,...,On)=S1,..,Sn∑P(O1,...,On∣S1,...,Sn)P(S1,...,Sn)=S1,..,Sn∑(k=1∏nP(Ok∣Sk))k=1∏nP(Sk∣Sk−1)
前一个是观测概率B，后一个是转移概率P。我们依次算一下，就能够得到最好的结果。这里面的B和P都是事先训练好的。
如果直接这么求，计算量太大，重复计算量太高。可以利用马尔科夫性简化计算
我们设一个中间变量
λ i ( k ) = P ( O 1 , O 2 , . . . , O k , S k = i ) λ j ( k + 1 ) = P ( O 1 , O 2 , . . . , O k , O k + 1 , S k + 1 = j ) \lambda_i(k) = P(O_1,O_2,...,O_k,S_k = i)\\ \lambda_j(k+1) = P(O_1,O_2,...,O_k,O_{k+1},S_{k+1} = j) λi(k)=P(O1,O2,...,Ok,Sk=i)λj(k+1)=P(O1,O2,...,Ok,Ok+1,Sk+1=j)
找一下递推关系
λ j ( k + 1 ) = P ( O 1 , O 2 , . . . , O k , O k + 1 , S k + 1 = j ) = ∑ i P ( O 1 , O 2 , . . . , O k , O k + 1 , S k + 1 = j , S k = i ) = ∑ i P ( O k + 1 , S k + 1 = j ∣ O 1 , . . . O k , S k = i ) P ( O 1 , . . . , O k , S k = i ) = ∑ i P ( O k + 1 , S k + 1 = j ∣ S k = i ) λ i ( k ) \lambda_j(k+1) = P(O_1,O_2,...,O_k,O_{k+1},S_{k+1} = j) \\ = \sum_iP(O_1,O_2,...,O_k,O_{k+1},S_{k+1} = j,S_k=i) \\ = \sum_i P(O_{k+1},S_{k+1}=j|O_1,...O_k,S_k = i)P(O_1,...,O_k,S_k=i) \\ = \sum_i P(O_{k+1},S_{k+1}=j|S_k=i) \lambda_i(k) λj(k+1)=P(O1,O2,...,Ok,Ok+1,Sk+1=j)=i∑P(O1,O2,...,Ok,Ok+1,Sk+1=j,Sk=i)=i∑P(Ok+1,Sk+1=j∣O1,...Ok,Sk=i)P(O1,...,Ok,Sk=i)=i∑P(Ok+1,Sk+1=j∣Sk=i)λi(k)
可以根据马尔科夫性忘记前面的状态。再用全概率公式变形
λ j ( k + 1 ) = ∑ i P ( O k + 1 , S k + 1 = j ∣ S k = i ) λ i ( k ) = ∑ i P ( O k + 1 ∣ S k + 1 = j , S k = i ) P ( S k + 1 = j ∣ S k = i ) λ i ( k ) = P ( O k + 1 ∣ S k + 1 = j ) ∑ i λ i ( k ) P ( S k + 1 = j ∣ S k = i ) = B j O k + 1 ( k + 1 ) ( ∑ i λ i ( k ) P i j ) \lambda_j(k+1) =\sum_i P(O_{k+1},S_{k+1}=j|S_k=i) \lambda_i(k) \\ = \sum_i P(O_{k+1}|S_{k+1}=j,S_k=i)P(S_{k+1}=j|S_k=i) \lambda_i(k) \\ = P(O_{k+1}|S_{k+1}=j)\sum_i \lambda_i(k) P(S_{k+1}=j|S_k=i) \\ = B_{jO_{k+1}}(k+1)(\sum_i \lambda_i(k)P_{ij}) λj(k+1)=i∑P(Ok+1,Sk+1=j∣Sk=i)λi(k)=i∑P(Ok+1∣Sk+1=j,Sk=i)P(Sk+1=j∣Sk=i)λi(k)=P(Ok+1∣Sk+1=j)i∑λi(k)P(Sk+1=j∣Sk=i)=BjOk+1(k+1)(i∑λi(k)Pij)
我们要求的概率可以表示为
P ( O 1 , . . . , O n ) = ∑ i P ( O 1 , . . . , O n , S n = i ) = ∑ i λ i ( n ) P(O_1,...,O_n) = \sum_i P(O_1,...,O_n,S_n=i) = \sum_{i}\lambda_i(n) P(O1,...,On)=i∑P(O1,...,On,Sn=i)=i∑λi(n)
这样做运算量可以节约很多的数量级，得到的这个结果叫做前向递推
Forward Recursion \text{Forward Recursion} Forward Recursion
前向递推的方法使得小计算量的语音识别在嵌入式系统上也能得到简单实现。
第二个计算，我们要计算状态的条件概率的最大值
m a x S 1 , . . . , S n P ( S 1 , . . . , S n ∣ O 1 , . . . , O n ) max_{S_1,...,S_n} P(S_1,...,S_n|O_1,...,O_n) maxS1,...,SnP(S1,...,Sn∣O1,...,On)
就是在已知某个观测的前提下，计算这个状态后验概率的最大值。
m
a
x
S
1
,
.
.
.
,
S
n
P
(
S
1
,
.
.
.
,
S
n
∣
O
1
,
.
.
.
,
O
n
)
⇔
m
a
x
S
1
,
.
.
.
,
S
n
P
(
S
1
,
.
.
.
,
S
n
,
O
1
,
.
.
.
,
O
n
)
max_{S_1,...,S_n} P(S_1,...,S_n|O_1,...,O_n) \Leftrightarrow max_{S_1,...,S_n} P(S_1,...,S_n,O_1,...,O_n)
maxS1,...,SnP(S1,...,Sn∣O1,...,On)⇔maxS1,...,SnP(S1,...,Sn,O1,...,On)
如果用贝叶斯公式的话，P(O1,…On)会出现在分母上，这是个常数，不需要关心，因此可以转化为求右边式子的最大值。
有一个标准解法，就是维特比迭代
Viterbi \text{Viterbi} Viterbi
μ i ( k ) = m a x S 1 , . . . , S k − 1 P ( O 1 , . . . , O k , S k = i , . . . , S 1 ) μ j ( k + 1 ) = m a x S 1 , . . . , S k P ( O 1 , . . . , O k + 1 , S k + 1 = j , . . . , S 1 ) = m a x S 1 , . . . , S k P ( O k + 1 , S k + 1 ∣ O k , . . . , O 1 , S k , . . . , S 1 ) P ( O k , . . . , O 1 , S k , . . . S 1 ) = m a x S 1 , . . . , S k P ( O k + 1 ∣ S k + 1 ) P ( S k + 1 = j ∣ S k = i ) P ( O k , . . . , O 1 , S k = i , . . . S 1 ) = m a x i P ( O k + 1 ∣ S k + 1 = j ) m a x S 1 , . . . , S k − 1 P ( S k + 1 = j ∣ S k = i ) P ( O k , . . . , O 1 , S k = i , . . . S 1 ) = m a x i ( B j , O k + 1 ( k + 1 ) P i j μ i ( k ) ) \mu_i(k) = max_{S_1,...,S_{k-1}}P(O_1,...,O_k,S_k=i,...,S_1) \\ \mu_j(k+1) = max_{S_1,...,S_{k}} P(O_1,...,O_{k+1},S_{k+1}=j,...,S_1) \\ =max_{S_1,...,S_{k}} P(O_{k+1},S_{k+1}|O_k,...,O_1,S_k,...,S_1) P(O_k,...,O_1,S_k,...S_1) \\ = max_{S_1,...,S_{k}} P(O_{k+1}|S_{k+1})P(S_{k+1}=j|S_k=i) P(O_k,...,O_1,S_k=i,...S_1) \\ = max_{i} P(O_{k+1}|S_{k+1}=j)max_{S_1,...,S_{k-1}}P(S_{k+1}=j|S_k=i) P(O_k,...,O_1,S_k=i,...S_1) \\ = max_i( B_{j,O_{k+1}(k+1)}P_{ij} \mu_i(k)) μi(k)=maxS1,...,Sk−1P(O1,...,Ok,Sk=i,...,S1)μj(k+1)=maxS1,...,SkP(O1,...,Ok+1,Sk+1=j,...,S1)=maxS1,...,SkP(Ok+1,Sk+1∣Ok,...,O1,Sk,...,S1)P(Ok,...,O1,Sk,...S1)=maxS1,...,SkP(Ok+1∣Sk+1)P(Sk+1=j∣Sk=i)P(Ok,...,O1,Sk=i,...S1)=maxiP(Ok+1∣Sk+1=j)maxS1,...,Sk−1P(Sk+1=j∣Sk=i)P(Ok,...,O1,Sk=i,...S1)=maxi(Bj,Ok+1(k+1)Pijμi(k))
维特比的递推和前向递推具有相似的结构
μ j ( k + 1 ) = m a x i B j , O k + 1 ( k + 1 ) ( P i j μ i ( k ) ) \mu_j(k+1) = max_i B_{j,O_{k+1}(k+1)}(P_{ij} \mu_i(k)) \\ μj(k+1)=maxiBj,Ok+1(k+1)(Pijμi(k))
λ j ( k + 1 ) = ∑ i B j O k + 1 ( k + 1 ) ( λ i ( k ) P i j ) \lambda_j(k+1)=\sum_i B_{jO_{k+1}}(k+1)( \lambda_i(k)P_{ij}) \\ λj(k+1)=i∑BjOk+1(k+1)(λi(k)Pij)
其中
μ i ( k ) = m a x S 1 , . . . , S k − 1 P ( O 1 , . . . , O k , S k = i , . . . , S 1 ) λ i ( k ) = P ( O 1 , O 2 , . . . , O k , S k = i ) \mu_i(k) = max_{S_1,...,S_{k-1}}P(O_1,...,O_k,S_k=i,...,S_1) \\ \lambda_i(k) = P(O_1,O_2,...,O_k,S_k = i) μi(k)=maxS1,...,Sk−1P(O1,...,Ok,Sk=i,...,S1)λi(k)=P(O1,O2,...,Ok,Sk=i)
