背景
贝叶斯统计推断是一个非常艺术的东西,他的先验结果(prior knowledge)是一种专家的经验,而贝叶斯公式给出的后验分布是贝叶斯推断的基本工具。每个大学生都应该学过数理统计,其中提及的贝叶斯统计推断都是一些基于简单先验和简单后验的结果。
但当我们一旦碰到复杂的先验(比如高维参数、复杂先验分布),我们用贝叶斯公式得出的后验分布将变得非常复杂,在计算上会非常困难。为此,先人提出了MCMC算法方便我们可以对任何后验分布进行计算或推断。其思想就是其名字:两个MC。
第一个MC: Monte Carlo(蒙特卡洛)。这个简单来说是让我们使用随机数(随机抽样)来解决计算问题。在MCMC中意味着:后验分布作为一个随机样本生成器,我们利用它来生成样本(simulation),然后通过这些样本对一些感兴趣的计算问题(特征数,预测)进行估计。
第二个MC:Markov Chain(马尔科夫链)。第二个MC是这个方法的关键,因为我们在第一个MC中看到,我们需要利用后验分布生成随机样本,但后验分布太复杂,一些package中根本没有相应的随机数生成函数(如rnorm(),rbinom()
)怎么办?答案是我们可以利用Markov Chain的平稳分布这个概念实现对复杂后验分布的抽样。
Markov Chain
为了能顺利阐述MCMC算法,这里就简单地讲一下所涉及到的马尔科夫链概念。因为都是我个人的理解,这里所讲不敢说都是正确的,希望各位童鞋能够独立思考。
1.什么是Markov Chain
随机过程${X_n},X_n$的状态空间为有限集或可列集,比如当$X_n=i$,就称过程在n处于状态i。
定义:如果对任何一列状态$i_0,i_1,...,i_{n-1},i,j$,及对任何n≥0,随机过程${X_n}$满足Markov性质:
$$P{X_{n+1}=j|X_0=i_0,...X_{n-1}=i_{n-1},X_n=i}\=P{X_{n+1}=j|X_n=i}$$
转移概率
$P{x_{n+1}=j|X_n=i}$成为Markov链的一步转移概率$P_{ij}^{n,n+1}$,当这个概率与n无关,则记之为$P_{ij}$
转移概率矩阵P
这个矩阵的元素就是一步转移概率$P_{ij}$
结论
一个Markov链可以由它的初始状态以及转移概率矩阵P完全确定。(证明略,自行百度或翻书)
2.n步转移概率
所谓n步转移概率就是从状态i走n步正好到状态j的概率,我们记为$P_{ij}^{(n)}$。
利用概率分割的思想,由基础概率论中全概率公式可以得到$$P_{ij}{(n)}=\Sigma_{k=0}{\infty}P_{ik}P_{kj}^{(n-1)}$$ 写成矩阵形式就是$$P^{(n)}=P\times P^{(n-1)}$$
进一步推广,我们就推出了著名的Chapman-Kolmogorov方程:$$P_{ij}{(n+m)}=\Sigma_{k=0}{\infty}P_{ik}{(n)}P_{kj}{(m)}$$写成矩阵形式就是$$P{(m+n)}=P{(m)}\times P^{(n)}$$
3.Markov链的极限定理和平稳分布
现在我们已经有了n步转移概率的概念,一个很简单的想法就是如果n趋向无穷会怎么样?这个问题也是后面极限分布以及平稳分布的来源,更是MCMC算法中第二个MC的关键。
要回答这个问题首先要掌握几个关键概念,我先列出来,如果不熟悉的可以自行百度或翻书:
- 互达性($i \leftrightarrow j$)
- 周期性(d(i))
- 常返与瞬过($f_{ii}=1$)
- 常返时($\T_i$)
3.1互达性、周期性
几个重要结论(证明略,自行百度,或者call me)
- $i \leftrightarrow j,\Rightarrow d(i)=d(j)$
- $若存在d(i)<\infty,则存在N,对所有n>N,有P_{ii}^{(nd(i))}>0$
- $P_{ji}^{(m)}>0,\Rightarrow 存在N,对所有n>N,有P_{ii}^{(m+nd(i))}>0$
- 对于非周期不可约Markov链的转移概率矩阵P,存在N,当n≥N时,$P^{(n)}>0$
3.2常返、瞬过、常返时
引入重要概率$f_{ij}{(n)}$表示从i出发在n步**首次**到达j的概率,我们约定$f_{ii}{(0)}=0$
$定义f_{ij}=\Sigma_{n=1}{\infty}f_{ij}{(n)}$表示从i出发最终到j的概率。
$定义:如果f_{ii}=1,则状态i是常返的,否则是瞬过的$
$我们记T_i为首次返回i的步长,且\mu_i=ET_i$
$若\mu_i=\infty,则称状态i是零常返,否则正常返$
3.3 极限定理
是时候回答上面那个问题了!(摘自网络共享PPT)
另外对于一个正常返、非周期的状态(也称为遍历ergodic),我们有结论:$$对所有i\leftrightarrow j,有limP_{ji}{(n)}=limP_{ii}{(n)}=\frac{1}{\mu_i}$$
到这里,我们可以想象了,n步转移概率矩阵最终的极限形式应该是由相同的行向量组成的!!!我们把那个行向量称为极限概率分布。 但是,问题还没完,要计算极限概率,就要根据极限定理,求出很多东西(如$f_{ii}^{(n)}$),实际中并不方便。所以,我们要引入平稳分布这个概念。最关键的部分来了!!!!
3.3 平稳分布
什么是平稳分布?它和求极限概率分布有什么关系呢?
定义:Markov链有转移概率矩阵P,如果有一个概率分布${\pi_i }满足\pi_j =\Sigma_{i=0}^{\infty}\pi_i P_{ij}$,则称为这个Markov链的平稳分布。这个定义用矩阵形式写出来就是π*P=π.
这个定义内容很丰富:如果一个过程的初始状态$X_0$有平稳分布π,我们可以知道对所有n,$X_n$有相同的分布π。再根据Markov性质可以得到,对任何k,有$X_n,X_{n+1},...,X_{n+k}$的联合分布不依赖于n,显然这个过程是严格平稳的,平稳分布也由此得名!!
重要定理
若一个markov链中所有的状态都是遍历的,则对所有i,j有$limP_{ij}{(n)}=\pi_j存在,且\pi={\pi_j,j≥0}就是平稳分布$。反之,拖一个不可约Markov链只存在一个平稳分布,且这个Markov链的所有状态都是遍历的,则这个平稳分布就是该Markov链的极限概率分布。$$limP_{ij}{(n)}=\pi_j$$
上面这条定理就给出了通过求平稳分布来求Markov链极限概率分布的简单方法。到这里我对第二个MC的解读也就结束了。下一期我们就要开始具体讲讲怎么应用MCMC算法对后验分布进行抽样和推断。