R语言与Markov Chain Monte Carlo(MCMC)方法学习笔记(1)

       蒙特卡洛方法被誉为20世纪最伟大的十大算法之一。它由美国拉斯阿莫斯国家实验室的三位科学家John von Neumann, Stan Ulam 和 Nick Metropolis于1946年提出。

       蒙特卡洛算法之所以那么有名，我的理解就是它利用随机模拟给出了一个十分普遍的求解许多问题近似解的办法。一个十分形象的例子是：在广场上画一个边长一米的正方形，在正方形内部随意用粉笔画一个不规则的形状，现在要计算这个不规则图形的面积，怎么计算?蒙特卡洛(Monte Carlo)方法告诉我们，均匀的向该正方形内撒N（N 是一个很大的自然数）个黄豆，随后数数有多少个黄豆在这个不规则几何形状内部，比如说有M个，那么，这个奇怪形状的面积便近似于M/N，N越大，算出来的值便越精确。在这里我们要假定豆子都在一个平面上，相互之间没有重叠。(撒黄豆只是一个比喻。)我们在《R语言与分类算法的绩效评估》一文计算AUC时用的就是这个方法。

        要说道Monte Carlo模拟，那么我们就得从随机数的产生开始说。

一、常见的抽样方法

        常见的抽样方法有许多，如直接抽样法（逆变换法）、拒绝抽样法、重要抽样法等。由于逆变换法过于简单，我们在这里就不再讨论了，我们先来看看拒绝抽样。

接受-拒绝抽样（Acceptance-Rejectionsampling)

       拒绝抽样的算法也十分简单，我们之所以会花一定的篇幅介绍它，是因为它是MCMC方法的一个基础。拒绝抽样的基本思想是，我们需要对一个分布f(x)进行采样，但是却很难直接进行采样，所以我们想通过另外一个容易采样的分布g(x)的样本，用某种机制去除掉一些样本，从而使得剩下的样本就是来自与所求分布f(x)的样本。

      它有几个条件：1）对于任何一个x，有f(x)<=M*g(x);也就是说我们的初步采样是必须包括进一步取样的全体的 2) g(x)容易采样；3) 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）的密度函数的最大值。显然这是满足拒绝抽样的几个条件。

对应的R代码给出如下：

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)

        One-sampleKolmogorov-Smirnov test

data:  z

D = 0.0317, p-value = 0.9026

alternative hypothesis: two-sided

        可见的确是生成了分布为beta（2,7）的随机数。我们可以用图来说明这个办法：

       可见g(x)如果与f(x)的形状相差较大时，效率是比较低的，本次运行的