MCMC 是马尔科夫链蒙特卡洛方法,英文是Markov Chain Monte Carlo。
首先说下什么是Monte Carlo方法。
MC方法的初衷是很多问题(比如求期望、积分)都不太好直接计算,所以MC通过随机采样的方式来近似求解这些问题。随机采样的方法很多,主要目的是生成一些服从目标分布(或者说近似分布)的样本,通过对这些采样的样本的分析来求解我们最终关心的问题。
而马尔科夫链这个方法又是另外一个方法。马尔科夫链是一个随机变量的序列,在这个序列上,当前时刻的随机变量如果只和上一个时刻的随机变量相关,不和历史上的序列相关,那么我们就说这个序列符合马尔科夫性。而具有这个性质的随机变量序列我们称之为马尔科夫链。马尔科夫链中比较重要的两个东西就是“状态转移概率分布”P和“状态概率分布”\(\pi(x)\)。
\[\pi_{n+1}(x)=\pi_n(x)P\]这两个分布决定了马尔科夫链中每个时刻的状态是啥,每个时刻怎么由当前·状态转移到下一个状态。如果存在某个状态分布,使得马尔科夫链怎么转移都还是在这个状态分布下,即: \(\pi(x)=\pi(x)P\)
那么这个状态分布我们称之为“平稳分布”。
至于马尔科夫链在什么条件下才会存在平稳分布后面再说。
假设马尔科夫链会收敛到平稳分布,那么我们可以构造这样一个马尔科夫链,使得其平稳分布就是我们想要求解的分布,然后再从这个马尔科夫链上采样,不就解决我们的问题了吗。这样,就诞生了MCMC方法。
那么怎样才能构造这样一个马尔科夫链,使得其平稳分布就是我们的目标分布呢?这就要介绍马尔科夫链的细致平稳方程(条件)。
\[\pi_ip_{ij}=\pi_jp_{ji} i,j=1,2,...\]根据细致平稳方程可以退出其符合平稳分布,所以满足细致平稳方程的分布就是马尔科夫链的平稳分布。
证明如下:
\[\sum_{i=1}^{\infty} \pi_ip_{ij} = \sum_{i=1}^{\infty} \pi_jp_{ji} = \pi_j \sum_{i=1}^{\infty} p_{ji}=\pi_j\]即 \(\pi P = \pi\)
那么我们面临的问题是如何根据平稳分布来构造概率转移矩阵,使得细致平稳方程成立。
在Metropolis-Hastings(M-H算法)中,状态转移矩阵是这么构建的。
我们已知一个比较容易进行采样的分布P,P也是另外一个马尔科夫链的状态转移矩阵。我们在P的基础上来构造状态转移矩阵。
因为有:
\[\pi(i)P(i,j)\pi(j)P(j,i)=\pi(j)P(j,i)\pi(i)P(i,j)\]所以 令
\[\alpha(i,j) = \pi (j)P(j,i)\] \[\alpha(j,i) = \pi (i)P(i,j)\]那么
\[\pi (i)P(i,j)\alpha(i,j)=\pi(j)P(j,i)\alpha(j,i)\]令
\[Q(i,j)=P(i,j)\alpha(i,j)\]那么
\[\pi(i)Q(i,j)=\pi(j)Q(j,i)\]显然,这个就是马尔科夫链的细致平稳方程。平稳分布是\(\pi(x)\),状态转移矩阵是Q,这也就是M-H算法中根据已知分布P构造出来的。
上面的P我们称之为建议分布,\(\alpha\)我们称之为接受率。实际执行时,我们从(0,1)的均匀分布中生成一个随机数,如果该随机数小于\(\alpha\),那么我们就接受由P采样的新状态,否则,就待在当前状态不动,即拒绝。
上面这个算法有一个缺点,就是\(\alpha\)可能很小,拒绝率高,导致采样效率比较低。M-H算法进行了改进。
\[\alpha(i,j) = min\{ \frac{\pi(j)P(j,i)}{\pi(i)P(i,j)},1\}\]这个改进相当于对原始的\(\alpha\)做了一个扩大,把其中一个扩大为1,提高了接受率。
但是,M-H采样还是有两个缺点,一是接受率在高维时计算量较大;二是有些数据的联合分布不好求,但是条件分布好求。这样,就有了最常用的Gibbs采样。
Gibbs采样每次只对随机变量中的一维进行采样,也就是对条件概率进行采样。可以证明,此时的接受率为1。
下面以一个二维随机变量来说明Gibbs采样。
假设二维概率分布p(x,y),考虑两点\(A(x_1,y_1),B(x_1,y_2)\)
有:
\[p(x_1, y_1) p(y_2 | x_1) = p(x_1) p(y_1 | x_1) p(y_2 | x_1) \\ p(x_1, y_2) p(y_1 | x_1) = p(x_1) p(y_2 | x_1) p(y_1 | x_1)\] \[p(x_1, y_1)p(y_2|x_1) = p(x_1, y_2)p(y_1|x_1)\]即 \(p(A)p(y_2|x_1) = p(B)p(y_1|x_1)\) .可以看到,在 \(x=x_1\) 时,如果以条件分布 \(p(y|x_1)\) 为转移概率,那么任意两点之间满足细致平稳条件。同样的,在 \(y=y_1\) 时细致平稳条件也成立
\[p(A)p(x_2|y_1) = p(C)p(x_1|y_1).\]所以我们构造如下的转移概率矩阵:
则有:
\[p(X)Q(X \rightarrow Y ) = p(Y )Q(Y \rightarrow X)\]上面是二维的情形。gibbs采样也可以推广到高维:
马尔科夫链方法就是定义一个马尔科夫链,使得最终的平稳分布是我们的目标分布。有了这个马尔科夫链之后,我们在这个马尔科夫链上进行随机游走,也就是采样。每个时刻都会产生一个样本。当步数足够多时,样本的分布就会趋近于平稳分布。在进入平稳分布之前的这个阶段称之为燃烧期(burn in)。那进入到平稳分布之后的样本就是我们需要的样本。我们可以基于这部分样本通过蒙特卡洛方法来求解我们关心的问题。
下面说下我在看相关章节的疑问:
1.拒绝采样中 目标采样的概率形式是否已知,C如何得到?
我理解应该是已知的。只有已知,才能确定是要拒绝还是接受。 但是C如何得到呢?
2.如何理解MCMC中的细致平稳方程
从公式上看,细致平稳方程是说从状态j到状态i的概率和从状态i到状态j的概率是相同的。所以,状态i上的概率是稳定的。
3.如何理解M-H算法中的接受率
参考:
1.https://www.cnblogs.com/pinard/p/6625739.html
2.http://blog.sciencenet.cn/blog-110554-1030316.html
3.《LDA数学八卦》