蒙特卡罗方法

MCMC是（Markov Chain Monte Carlo）缩写，中文马尔科夫链蒙特卡罗。

蒙特卡罗方法

Monte Carlo思想

最早的蒙特卡罗方法是为了求各或积分问题，比如

[ heta=int_a^bf(x)dx ]

使用蒙特卡罗方法求得积分的近似值。在([a,b])区间上随机采样一个近似点(x_0)，然后用(f(x_0))代表区间上所有(f(x))的值，这样近似求解的结果为

[(b-a)f(x_0) ]

进一步，采用([a,b])区间上的n个点近似代替，则近似结果为

[frac{b-a}{n}sum_{i=0}^{n-1}f(x_i) ]

此时假设区间是均匀分布的，如不是均匀分布，则结果差很多。

我们假设(x)在区间([a,b])上的概率分布函数为(p(x))，那么我们定积分可以表示为

[ heta=int_a^bf(x)dx=int_a^bfrac{f(x)}{p(x)}p(x)dxapproxfrac{1}{n}sum_{i=0}^{n-1}frac{f(x_i)}{p(x_i)} ]

上式右边即代入了蒙特卡罗方法。可以看出，假设([a,b])间是均匀分布的，那么(p(x_i)=frac{1}{b-a})，代入后可得

[frac{1}{n}sum_{i=0}^{n-1}frac{f(x_i)}{{1}/{(b-a)}}=frac{b-a}{n}sum_{i=0}^{n-1}f(x_i) ]

所以，使用蒙特卡罗方法的关键是得到x的概率分布p(x)。

概率分布采样

如果得到了x的概率分布，我们可以基于此去采样。对于均匀分布，常用的采样方法即uniform(0,1)，通过线性同余发生器，可以很容易得到伪随机样本。地于其他常见的分布，也可以通过函数转换而得。

接受－拒绝采样

对于不常见的分布，一个可行的办法是接受－拒绝采样方法。即然p(x)太复杂没有办法直接采样，那么我们选择一个可采样的分布比如高斯分布，然后按照一定的方法拒绝某些样本，以达到近似p(x)的目的。

采样过程中，设定一个常用的概率分布函数q(x)，以及一个常数k，使得p(x)总在kq(x)的下方。采用均匀分布随机采样的方法，如果这此采样落在了p(x)的上方，则拒绝此次采样，否则接受，直到采够n个样本，最后结果为

[frac{1}{n}sum_{i=1}^{n-1}frac{f(z_i)}{p(z_i)} ]

整个过程就是用q(x)来模拟p(x)的概率分布的过程。

蒙特卡罗总结

蒙特卡罗主要解决一些不常见分布的问题，但对于一些多维分布或是条件分布的问题无法使用。

马尔科夫链

Markov Chain主要思想

马尔科夫链，它假设某一时刻的状态转移概率只信赖于它的前一个状态。在一个系统中，只要求出任意两个状态之间的转换概率，这个马尔科夫模型就定了。

马尔科夫链状态转移矩阵的性质

马尔科夫模型的状态转移矩阵收敛到稳定的概率分布与我们的初始概率分布无关。

基于马尔科夫链采样

我们达到平稳分布后，就很容易进行采样。

马尔科夫采样过程：

1. 输入马尔科夫状态转移矩阵，设定状态转移次数(n_1)和需要采样的次数(n_2)

2. 从任意简单概率分布采样得到初始状态值(x_0)

3. for t=0 to (n_1+n_2-1)：从条件概率分布(p(x|x_t))中采样得到样本(x_{t+1})

   样本((x_{n_1},x_{n_1+1},...,x_{n_1+n_2+1}))即为我们需要的平稳分布的样本集。

马尔科夫链采样总结

我们得到马尔科夫状态转移矩阵，就可以采样得到样本集，然后进行蒙特卡罗模拟。

关键问题是又转换到如何得到马尔科夫状态转移矩阵，下面介绍MCMC方法解决该问题。

MCMC采样

马尔科夫链细致平衡条件

如果非周期马尔科夫链的状态转移矩阵P和概率分布(pi(x))对所有的i,j满足

[pi(i)P(i,j)=pi(j)P(j,i) ]

则称概率分布(pi(x))是状态转移矩阵P的平稳分布。由平稳分布条件有：

[sum_{i=1}^{infty}pi(i)p(i,j)=sum_{i=1}^{infty}pi(j)P(j,i)=pi(j)sum_{i=1}^{infty}P(j,x)=pi(j) ]

将上式用矩阵表示为

[pi P=pi ]

也就是说，只要我们找到了可以使概率分布(pi(x))满足细致平稳条件的矩阵P即可。如何找P呢

随机找一个状态转移矩阵Q，很难满足细致平稳条件。

[pi(i)Q(i,j) epi(j)Q(j,i) ]

MCMC采样

对细致平稳条件进行改造

[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) ]

代入上式两边即可约为衡等式，即可使原不等式变为等式

[pi(i)Q(i,j)＝pi(j)Q(j,i) ]

即有

[P(i,j)=Q(i,j)alpha(i,j) ]

也就是说，我们的目标矩阵P，可以通过任意一个矩阵Q，乘以(alpha(i,j))得到，我们称(alpha(i,j))为接受率。

MCMC过程如下：

1. 输入矩阵Q，平稳分布(pi)，状态转移次数阈值(n_1)，需要样本数(n_2)
2. 从任意简单概率分布采样得到样本(x_0)
3. for t=0 to (n_1+n_2+1):
   1. 从条件概率分布(Q(x|x_t))中采样得到样本x
   2. 从均匀分布采样u~uniform[0,1]
   3. 如果(u<alpha(x_t,x_*)=pi(x_*)Q(x_*,x_t))则接受转移(x_t o x_*)即(x_{t+1}=x_*)
   4. 否则不接受转移，即(x_{t+1}=x_t)

MH采样

略