Markov Chain Monte Carlo(MCMC) 方法

Monte Carlo 方法

假设我们要求一个原函数并不明确的函数\(f(x)\)的在某个区间\([a,b]\)上的积分

\(\theta = \int_{a}^b f(x)dx\)

因为\(f(x)\)的原函数不知道，所以无法用牛顿-莱布尼茨公式计算。这里采用一种称为monte carlo的方法来模拟近似求解，它的思想如下，首先将待求的式子化为

\(\theta = \int_a^b f(x)dx = \int_a^b \frac{f(x)}{p(x)}p(x)dx=E_{x\sim p(x)}[\frac{f(x)}{p(x)}]\)

这个式子将积分看做是随机变量\(\frac{f(x)}{p(x)}\)的期望，再由强大数定律有

\(\theta = E_{x\sim p(x)}[\frac{f(x)}{g(x)}] = \lim_{n \rightarrow \infty}\frac1n\sum_{i=0}^n \frac{f(x_i)}{p(x_i)}\quad (1)\)

上述式子成立要求\(\frac{f(x_i)}{p(x_i)}\)相互独立同分布(期望存在)。那么根据这个式子我们就可以独立地产生一组服从\(p(x)\)这个密度函数的随机数，然后按照上述等式最右边的式子来计算得到\(\theta\)的估计值. 这种思想就是Monte Carlo方法的核心思想。

现在一个问题就是如何产生一组服从特定密度函数\(p(x)\)的相互独立的随机数，而MCMC方法可以解决这个问题
PS: 实际上这里的独立性要求其实不需要，因为可以证明以\(P(x)\)为平稳概率的马氏链满足上面的(1)式(也就是说尽管不一定有\(\frac{f(x_i)}{p(x_i)}\),可以证明(1)也同样成立)，所以上面的问题可 进一步改为

现在的问题就是如何产生一组服从特定密度函数\(p(x)\)的随机数，MCMC方法可以解决这个问题

Markov Chain Mote Carlo(MCMC)

要用MCMC方法，必须要找到一个平稳分布是\(\pi(i)\)的马氏链，更为具体一点就是 通过已知的平稳分布\(\pi(i)\)来确定一个马氏链转移概率\(p(i,j)\)(马氏链除了定义状态以外就是定义转移概率了)，使得该马氏链在这个转移概率下经过长时间转移后有平稳分布\(\pi(i)\)。目前我们可以知道的平稳分布\(\pi(i)\)与转移概率\(p(i,j)\)的关系是下式

\(\pi(i)p(i,j) = \pi(j)p(i,j)\)

这是时逆马氏链的条件，鉴于我们的后面的讨论都基于这个式子就可以知道MCMC方法构造的马氏链一定是时逆马氏链(时逆马氏链一定有稳态分布)。我们从这个式子出发，对任意一个以\(Q(i,j)\)为转移概率的马氏链来说上式并不一定成立，但是我们让他们强行成立的话有

\(\pi(i)Q(i,j)\alpha(i,j) = \pi(j)Q(j,i)\alpha(j,i)\)

其中\(\alpha(i,j) = \pi(j)Q(j,i)\ ;\alpha(j,i) = \pi(i)Q(i,j)\), 所以我们可以看到实际上构造的马氏链是以\(Q(i,j)\alpha(i,j)\)为转移概率的马氏链，这里的\(\alpha(i,j)\)又称为接受率。所以下面我们在进行抽样的时候状态之间的转移应该包括与\(Q(i,j)\)对应的过程和与\(\alpha(i,j)\)对应的过程，具体的MCMC算法的过程如下

1. 输入一个马氏链的转移概率\(Q(i,j)\)以及期望的平稳分布\(\pi(i)\), 当然还有状态转移所需要的次数\(n_T\)以及要采集随机数的个数\(n_r\)

2. 按照一定的概率分布随机产生一个状态\(x_0\)

3. 进入循环，循环次数为(\(n_T+n_r\)), 每次循环做以下的一些事情

   a. 从条件概率分布\(Q(i,j)\)中采样得到\(x^*\)（这里采样可以使用接受-拒绝采样来得到，可见

   (https://www.cnblogs.com/pinard/p/6625739.html)

   b. 以\(\alpha(i,j)\)的概率来接受\(x^*\)(生成一个[0,1]之间的均匀分布的随机数\(u\), 如果\(u<\alpha(i,j)\)就接受\(x^*\),即\(x_{n+1} = x^*\);否则就拒绝，\(x_{n+1} = x_n\))

Metropolis-Hasting 算法

M-H算法是针对MCMC算法里面\(\alpha(i,j)\)太小而导致马氏链需要很长的一段时间才收敛的问题而提出的。在MCMC算法中\(\alpha(i,j)\)是介于[0,1]之间的，而在M-H算法中\(\alpha(i,j)\)或者\(\alpha(j,i)\)中有一个为1。因为

\(\pi(i)Q(i,j)\alpha(i,j) = \pi(j)Q(j,i)\alpha(j,i)\)

等式两边同时扩大\(\frac{1}{\alpha(i,j)}\)倍或者同时扩大\(\frac{1}{\alpha(j,i)}\)倍等式不变，for example,两边同时扩大\(\frac1{\alpha(j,i)}\)倍, 就有

\(\pi(i)Q(i,j)\frac{\alpha(i,j)}{\alpha(j,i)} = \pi(j)Q(j,i) \rightarrow\pi(i)Q(i,j)\alpha_{MH}(i,j) = \pi(j)Q(j,i)\alpha_{MH}(j,i)\)

而新的\(\alpha_{MH}(i,j) = \frac{\alpha(i,j)}{\alpha(j,i)} = \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)}\)，而新的\(\alpha_{MH}(j,i) = 1\). 但是应该记住的是不管是不是新的，接受率都不能大于1，所以归纳一下可以得到\(\alpha_{MH}(i,j) = \min\{\frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\}\), 这样M-H算法构造马氏链的转移概率就是\(P(i,j) = Q(i,j)\alpha_{MH}(i,j)\)。具体的M-H算法采取一组服从某个概率分布\(\pi(i)\)的随机数如下所示

1. 输入一个马氏链的转移概率\(Q(i,j)\)以及期望的平稳分布\(\pi(i)\), 当然还有状态转移所需要的次数\(n_T\)以及要采集随机数的个数\(n_r\)

2. 按照一定的概率分布随机产生一个状态\(x_0\)

3. 进入循环，循环次数为(\(n_T+n_r\)), 每次循环做以下的一些事情

   a. 从条件概率分布\(Q(i,j)\)中采样得到\(x^*\)

   b. 以\(\alpha_{MH}(i,j)\)的概率来接受\(x^*\)

Gibbs 算法

Gibbs算法主要应用于多维随机向量的采样。具体可见 https://www.cnblogs.com/pinard/p/6645766.html

PS: 其中要求联合稳态分布的条件分布，可以用联合分布除以对应的边缘分布，而边缘分布可以通过对所有其他变量积分(连续)或者对所有其他变量的可能情况求和得到(离散)。