给定一个的概率分布 p(x) , 我们希望产生服从该分布的样本。前面介绍过一些随机采样算法（如拒绝采样、重要性采样）可以产生服从特定分布的样本，但是这些采样算法存在一些缺陷（如难以选取合适的建议分布，只适合一元随机变量等）。下面将介绍一种更有效的随机变量采样方法：MCMC 和 Gibbs采样，这两种采样方法不仅效率更高，而且适用于多元随机变量的采样。
如果一定条件下，马尔可夫链可以收敛到平稳分布。一个绝妙的想法是：如果能构造一个转移矩阵为 P 的马尔可夫链，是其平稳分布刚好是 p(x) 。那我们就可以从任何一个初始状态 x0 出发，沿着马尔可夫链转移，得到一个转移序列 x0 , x1 , x2 ,…, xn , xn+1 , xT , 如果马尔可夫链在第 n 布已经收敛，我们就得到了 p(x) 的样本 xn,xn+1…xT , 过程如下所示：
设置 t=0 生成一个随机状态 u ， 设定初始状态 x0=u重复t=t+1
根据转移概率 p(xt|xt−1) , 采样得到 u
设置 xt=u直到 t=T
这种方法在1953年被 Metropolis 想到了，为了研究粒子系统的平稳性质，Metropolis 考虑了常见的冷静的背包分布的采样问题，首次提出了基于马尔可夫链蒙特卡罗的采样方法，即 Metropolis 算法，并在计算机上进行了编程实现。Metropolis 算法是首个普适的采样方法，并启发了一系列MCMC方法，所以人们把它作为随机模拟技术腾飞的起点。Metropolis 的这篇论文被收录在《统计学中的重大突破》，Metropolis 算法当选为二十世纪十个最重要的算法之一。
Metropolis 算法假定目标概率分布为 p(x) , 我们已经有一个转移矩阵为 Q （q(i,j) 表示从状态 i 转移到状态 j 的概率）的马尔可夫链。通常情况下不满足马尔可夫链的细致平稳条件：
所以 p(x) 通常不是这个马尔可夫链的平稳分布。我们可否对该外向的向日葵链做一个改造，使得细致平稳条件成立呢？譬如，我们引入一个接受率 α(i,j) , 我们希望：
p(i)q(i,j)α(i,j)=p(j)q(j,i)a(j,i)
取什么样的 α(i,j) 才能让上式成立呢？最简单的，根据对称性，我们可以取
α(i,j)=p(j)q(j,i) α(j,i)=p(i)q(i,j)
此时细致平稳条件成立：
p(i)q(i,j)α(i,j)q′(i,j)=p(j)q(j,i)α(j,i)q′(j,i)
于是我们将原来具有转移矩阵 Q 的一个普通的外向的向日葵链，改造成了具有转移矩阵 Q′ （ q′(i,j) 表示从状态 i 转移到状态 j 的概率）的外向的向日葵链。而此外向的向日葵链的平稳分布刚好就是目标概率分布 p(x) ！
接受率 α(i,j) 的理解
由于 α(i,j)=p(j)q(j,i)≤1 ，因此
细致平稳条件的直观含义如下图所示：
一般介绍MCMC的资料中没有说明状态自我转移永远是满足细致平稳条件的，自己转给自己，当然转入和转出相等。 因此可以增加自我转移的概率 q′(i,i) ，使 ∑j∈Sq′(i,j)=1 成立 。Metropolis算法细致平稳条件如下图所示：
从以上分析可以发现， α(i,j) 表示在原来转移矩阵为 Q 的外向的向日葵链中，从状态 i 以概率 q(i,j) 转移到状态 j 的时候，我们以 α(i,j) 接受这个转移。而以概率 q(i,j)(1−α(i,j)) 拒绝这个转移（进行自我转移）。因此称 α 为接受率。
Metropolis算法流程假定目标分布为 p(x) ，我们已经有了一个转移矩阵为 Q （对应元素为 q(i,j) )，经典MCMC算法流程如下所示
以上分析中 p(x) , q(x|y) 说的都是离散的情形，事实上即便这两个分布式连续的，以上算法仍然有效（此时 p(x) 表示概率密度， q(x|y) 表示条件概率密度），因此可以生成服从连续的概率分布 p(x) 的样本。
注：MCMC算法是以于马尔可夫链为基础的。马尔可夫链通常研究的是时间和状态都是离散的随机过程。而连续分布实际上对应的状态是连续的。那么MCMC算法可以对连续分布进行采样的原因又是什么？？
Metropolis-Hastings 算法以上的Metropolis采样算法已经能漂亮地工作了，不过它有个小问题：外向的向日葵链在转移的过程中的接受率 α(i,j) 可能偏小，采样过程中容易原地踏步（自我转移），拒绝大量的跳转。这会导致外向的向日葵链收敛到稳态分布需要很长的时间。有没有办法提升接受率呢？
假设在Metropolis算法中 α(i,j)=0.1 ， α(j,i)=0.2 ，满足细致平稳条件，下式成立：
将上式两边同时扩大5倍
p(i)q(i,j)×0.5=p(j)q(j,i)×1
看！我们将接受率 α(i,j) ， α(j,i) 放大了5倍，而细致平稳条件仍然成立！这启发我们可以把接受率 α(i,j) ， α(j,i) 同比例放大，最优的情况是使得两个数中较大的一个放大到1（因为接受率要求小于等于1）。所以我们可以取
α(i,j)=min{p(j)q(j,i)p(i)q(i,j),1}
于是，通过对Metropolis算法中的接受率的微小改造，我们就得到了如下教科书中最常见的Metropolis-Hastings 算法。
注：对未归一化的概率分布 p(x)=p~(x)Xp ( Xp 未知)。MH算法仍然适用。只需将 Algorithm 6 中的接受率可变换为：
例子
采用 MH算法生成服从Cauchy分布的样本。MH算法采样的过程如下所示：
当前的状态为 xt
利用建议分布（这里使用正态分布）进行采样得到 y
若接受转移，则 xt+1=y
matlab代码如下：
%%M-H算法T = 500; %the maximum number of iterationssigma = 1;%the deviation of normal proposal densityx_min=-10; x_max=10 %define a range for starting valuesx = zeros(1,T);%init storage space for our samplesseed=1; rand('state',seed); randn('state',seed); %random seedx(1) = unifrnd(x_min,x_max);%%Start samplingt=1;while t<T %Iterate until we have T samples t=t+1; %propose a new value for theta using a normal proposal density y = normrnd(x(t-1),sigma); %Calculate the acceptance ratio alpha = min([1,(cauchy(y)*normpdf(x(t-1),y,sigma))/(cauchy(x(t-1))*normpdf(y,x(t-1),sigma))]); %draw a uniform deviate from [0,1] u=rand; %accept this proposal? if u<alpha x(t)=y; else x(t)=x(t-1); endend%% Display histogram of our samplesfigure(1); clf;subplot(3,1,1);nbins=200;thetabins = linspace(x_min,x_max,nbins);counts = hist(x,thetabins);bar(thetabins,counts/sum(counts),'k');xlim([x_min x_max]);xlabel('x'); ylabel('p(x)');%% Overlay the theoretical densityy=cauchy(thetabins);hold on;plot(thetabins,y/sum(y),'r--','LineWidth',3);set(gca,'YTick',[]);%% Display history of our samplessubplot(3,1,2:3);stairs(x,1:T,'k-');ylabel('t');xlabel('x');set(gca,'YDir','reverse');xlim([x_min x_max]);实验结果：
随机采样方法整理与讲解（MCMC、Gibbs Sampling等
MCMC方法的理解
Computational Statistics with Matlab