马尔科夫蒙特卡洛算法(MCMC)

趁着周末，学习了此算法。一个重要的作用就是用来模拟目标分布的样本。下面看看具体情况。

---

1.名词解释
MCMC方法就是*构造合适的马尔科夫链进行抽样而使用蒙特卡洛方法进行积分计算,既然马尔科夫链可以收敛到平稳分布。我们可以建立一个以π为平稳分布的马尔科夫链，对这个链运行足够长时间之后，可以达到平稳状态。此时马尔科夫链的值就相当于在分布π(x)中抽取样本。利用马尔科夫链进行随机模拟的方法就是MCMC。

第一个MC: Monte Carlo(蒙特卡洛)。这个简单来说是让我们使用随机数（随机抽样）来解决计算问题。在MCMC中意味着：后验分布作为一个随机样本生成器，我们利用它来生成样本(simulation)，然后通过这些样本对一些感兴趣的计算问题（特征数，预测）进行估计。

第二个MC：Markov Chain(马尔科夫链)。第二个MC是这个方法的关键，因为我们在第一个MC中看到，我们需要利用后验分布生成随机样本，但后验分布太复杂，当这些样本独立时，利用大数定律样本均值会收敛到期望值。如果得到的样本是不独立的，那么就要借助于马尔科夫链进行抽样，利用Markov Chain的平稳分布这个概念实现对复杂后验分布的抽样。

2.马尔科夫链的定义如下：
设θ(t)是一个随机过程，如果它满足下面的性质：
p(θ(t+h)=xt+h|θ(s)=xs,s≤t)=p(θ(t+h)=xt+h|θ(t)=xt), 任意的h>0。

这个定义又称为马尔科夫性质，对一个马尔科夫链来说，未来状态只与当前t时刻有关，而与t时刻之前的历史状态无关（条件独立）。

马尔科夫链的一个很重要的性质是平稳分布。简单的说，主要统计性质不随时间而变的马尔科夫链就可以认为是平稳的。数学上有马氏链收敛定理，当步长n足够大时，一个非周期且任意状态联通的马氏链可以收敛到一个平稳分布π(x)。这个定理就是所有的MCMC方法的理论基础。

结论：一个Markov链可以由它的初始状态以及转移概率矩阵P完全确定。

3.什么是平稳分布？它和求极限概率分布有什么关系呢？

定义：Markov链有转移概率矩阵P，如果有一个概率分布{πi}满足,则称为这个Markov链的平稳分布。这个定义用矩阵形式写出来就是π*P=π.

 这个定义的含义：如果一个过程的初始状态X0有平稳分布π,我们可以知道对所有n，Xn有相同的分布π。再根据Markov性质可以得到，对任何k，有
 Xn,Xn+1,...,Xn+k的联合分布不依赖于n，显然这个过程是严格平稳的，平稳分布也由此得名！！

4.拒绝抽样

基本思想是，我们需要对一个分布f(x)进行采样，但是却很难直接进行采样，所以我们想通过另外一个容易采样的分布g(x)的样本，
用某种机制去除掉一些样本，从而使得剩下的样本就是来自与所求分布f(x)的样本。

具体的采样过程如下：

1. 对于g(x)进行采样得到一个样本xi, xi ~ g(x);
2. 对于均匀分布采样 ui ~ U(a,b);
3. 如果ui<= f(x)/[M*g(x)], 那么认为xi是有效的样本；否则舍弃该样本； （# 这个步骤充分体现了这个方法的名字：接受-拒绝）
4. 反复重复步骤1~3，直到所需样本达到要求为止。

示例：产生服从beta(2,7)的随机数。提议分布g取为均匀分布,常数M取为beta(2,7)的密度函数的最大值。

a <- 2
b <- 7
xmax <- (a-1)/(a+b-2)
dmax <- xmax^(a-1)*(1-xmax)^(b-1)*gamma(a+b)/(gamma(a)*gamma(b))
y <- runif(1000)
x <- na.omit(ifelse(runif(1000) <= dbeta(y,a,b)/dmax,y,NA))

查看ks检验:

z <- x[1:323]
ks.test(z,"pbeta",2,7)
#接受域图
plot(dbeta(x,2,7)~x,type = "