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)}]=\underset{n\to \mathrm{\infty }}{lim}\frac{1}{n}\sum _{i=0}^n\frac{f(x_i)}{p(x_i)}(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^{\ast }$（这里采样可以使用接受-拒绝采样来得到，可见

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

   b. 以$\alpha (i,j)$的概率来接受$x^{\ast }$(生成一个[0,1]之间的均匀分布的随机数$u$, 如果$u<\alpha (i,j)$就接受$x^{\ast }$,即$x_{n+1}=x^{\ast }$;否则就拒绝，$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,两边同时扩大$\frac{1}{\alpha (j,i)}$倍, 就有

$\pi (i)Q(i,j)\frac{\alpha (i,j)}{\alpha (j,i)}=\pi (j)Q(j,i)\to \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^{\ast }$

   b. 以${\alpha }_{MH}(i,j)$的概率来接受$x^{\ast }$

Gibbs 算法

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

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