Metropolis-Hastings algorithm
1. 随机模拟的基本思想
假设我们有一个矩形区域(R)，面积为(S_0)。在此区域中，有一个不规则区域(M)，其面积(S)待求。
方法1：把不规则区域(M)划分为多个小的规则区域，由这些规则区域的面积总和(S')近似。
方法2：抓一把黄豆，均匀铺在(R)中，再统计落在不规则区域(M)中的黄豆比例。该比例可近似(frac{S}{S_0})。
方法2采用的是抽样（采样）解决问题的思想，即随机模拟。
因此，随机模拟的关键，在于如何将实际复杂问题，转化为抽样可以解决的问题。当然，这非常灵活，也考验我们的创造力。
我们接下来的核心问题，是利用计算机可直接实现的均匀抽样，借助随机模拟的方法，实现满足特定概率分布(p(x))的抽样。
2. 拒绝抽样
当我们的概率分布很简单时，如([0,1])上的均匀分布，抽样是很好实现的。现成的随机算法非常多。
但是，对于其他更复杂的分布，如果我们还想借助简单的均匀随机抽样，就需要换换思路了。
例如上图中的(p(z))（忽略波浪号），是较为复杂的概率密度函数。
我们先构造比较简单的、容易抽样的高斯分布(q(z))，乘一个常数(k)，使之满足：
拒绝抽样的方法是：
-
对高斯分布的样本进行采样，得到一个(z_0)，如图；
-
在([0,kq(z_0)])上均匀采样，得到(u_0)；
-
如果(u_0 > p(z))，就拒绝该采样结果；反之接受；
-
重复直到接受足够多的样本。
我们理解一下为什么拒绝抽样是可行的，这对我们之后运用随机模拟思想非常关键！
-
该简单分布与高斯分布的趋势大致相同，中间大两头小，因此取到中间的(z)的概率比较大；
-
通过(u_0)进一步筛选(z)，(p(z))较低时拒绝率较大，反之接受率较高，大致符合(p(z))分布要求。
一句话，通过简单抽样和设置接受率，“强迫”结果趋于复杂分布。
这样，我们就通过简单的高斯分布采样，以及计算机可直接实现的均匀采样，间接实现了复杂的(p(z))采样。
然而，在高维情况下，该方法的劣势很明显：
-
合适的简单分布(q(z))很难找到。
-
合适的(k)值很难找到。
这两个问题会导致拒绝率极高，计算效率很低。其根本缺陷在于：样本之间是独立无关的。
3. Metropolis-Hastings抽样
3.1 引入思想
结合马氏链知识（参考信息论相关书籍），我们知道：
对于遍历的马尔可夫链，当其转移次数足够多时，将会进入稳态分布(pi (x))，即各状态的出现概率趋于稳定。
进一步，假设变量(oldsymbol X)所有可能的取值，构成某遍历马氏链的状态空间。
则无论从什么状态出发，只要转移步数足够大，该马氏链将趋于稳定，即逐渐开始依次输出符合(p(x))分布的状态(X_i)。
这样，通过收集马氏链收敛后产生的(X_i)，我们就得到了符合分布(p(x))的样本。
此时稳态分布即(pi (x)=p(x))。
比如，我们希望抽样样本满足指数分布。此时，整数0到正无穷都是状态，但整数0的出现概率最大，随着数增大，出现概率越来越小。我们构造一个稳定输出服从指数分布状态的马氏链，则得到的样本服从指数分布。
这个绝妙的想法在1953年由Metropolis提出。为了研究粒子系统的平稳性质，Metropolis 考虑了物理学中常见的波尔兹曼分布的采样问题，首次提出了基于马氏链的蒙特卡罗方法（Metropolis算法），并在最早的计算机上编程实现。
Metropolis 算法是首个普适的采样方法，并启发了一系列 MCMC方法，所以人们把它视为随机模拟技术腾飞的起点。 Metropolis的这篇论文被收录在《统计学中的重大突破》中， Metropolis算法也被选为二十世纪的十个最重要的算法之一。
我们接下来介绍的MH算法是Metropolis 算法的一个常用的改进变种。
由于马氏链的收敛性质主要由转移矩阵(oldsymbol P)决定, 所以基于马氏链做采样的关键问题是：如何构造转移矩阵(oldsymbol P)，使平稳分布恰好是我们要的分布(p(x))。
3.2 理论基础：细致平稳条件
定理——细致平稳条件：
若非周期马氏链的转移矩阵(oldsymbol P)和分布(pi (x))满足：
则分布(pi (x))是马氏链的平稳分布。
证明：
这说明：
因此(pi (x))是平稳分布。
物理意义：
对于任意两个状态(i,j)，从状态(i)流失到状态(j)的概率质量，等于反向流动的概率质量，因此是状态概率是稳定的。
3.3 M-H算法实现
假设我们已经有了一个转移矩阵(oldsymbol Q)。一般情况下，该转移矩阵不满足细致平稳条件：
我们引入接受率(alpha(i,j))，满足：
其中：
-
(P(i))是状态(i)出现的概率，是假设马氏链满足复杂分布时输出状态的概率。
-
(oldsymbol Q)是初始转移概率矩阵，满足任意一个给定的简单分布，便于抽样即可。
显然，接受率最简单的构造方法是：
此时，细致平稳条件成立，该马氏链输出的状态满足稳态分布(p(x))！
综上，有几个要点必须要理解，有助于我们编程实现：
-
我们引入接受率，使得该马氏链以转移概率(oldsymbol Q_{ij}alpha(i,j))从状态(i)转移到状态(j)。
-
该转移概率的实现方式是：我们先按(oldsymbol Q_{ij})生成新状态，再按(alpha(i,j))的概率接受转移结果，乘积即为(oldsymbol Q_{ij}alpha(i,j))转移概率。
-
一定要理解这个接受率！在拒绝抽样中，如果(u_0 > p(z))，我们会放弃抽样结果，不纳入最终的统计。正是因为如此，抽样才会逼近复杂分布。而这里的接受是“接受转移”，如果(u>alpha)，意味着(t+1)时刻状态和(t)时刻相同，并且应该纳入统计而不是放弃结果！！
-
由于状态出现概率和转移概率满足细致平稳条件，因此状态出现概率即稳态分布概率。长期转移后，我们会得到想要的抽样分布效果。
图解转移过程：
要注意的是，当马氏链处在状态(i)时，我们并不知道如何选择下一个状态(j)。
我们在这里采用抽样的方法，并借助接受率，“强迫”转移过程逼近理想的遍历马氏链转移过程。
算法描述如下：
-
初始化马氏链状态(X_0=x_0)
-
从(t=0)开始，重复以下步骤：
-
假设该(t)时刻状态为(X_t=x_t)；
-
对(Q_{x_t,x})抽样，随机得到一个状态(x_{next})；
-
在([0,1])上均匀抽样，得到一个(u)；
-
若(u<alpha(x_t,x_{next})=P(x_{next})oldsymbol Q_{x_{next},x_t})，则接受转移(X_{t+1}=x_{next})；
-
否则不转移，即(X_{t+1}=x_t)。
-
3.4 算法升级
如果(alpha(x_t,x_{next}))太小，会导致转移很难发生，马氏链收敛过慢。
我们回到细致平稳条件：
我们希望在保持条件成立的前提下，使接受率尽量增大。
由于接受率的本质是概率，因此(alpha(i,j))和(alpha(j,i))中的较大者不应大于1。那么我们作下述改进即可：
解释：
-
如果(alpha(i,j))更大，那么应设为1。而第一项大于1，因此(alpha(i,j)=1)，结束。
-
如果(alpha(j,i))更大，那么为了等式成立，(alpha(i,j))必须等于buf。而buf(<1)，得逞～
综上，我们得到了升级版算法：
3.5 仿真实验
仿真目标：
用标准正态分布模拟指数分布（(lambda=0.1, mu=frac{1}{lambda}=10)）。
简化：
由于高斯分布是一个径向基函数（取值仅仅依赖于离原点距离的实值函数），因此该矩阵是转置对称的，(oldsymbol Q_{ij})=(oldsymbol Q_{ji})
代码实现：
function SamplefromExponential
clc;close all;clear;
% 随意设
data_size=100000; % 越大越准
start=1; % 非负即可
% 初始化
data=start;
index=1;
% 直接生成多个服从均值为10的指数分布的样本
data_dir=exprnd(10,[1,data_size]);
figure;
subplot(1,2,1)
histogram(data_dir,'BinWidth',1);
axis([0,100,0,data_size/10]);
% 随机模拟
while(index<data_size)
u=rand; % 均匀抽样
data_next=normrnd(data(index),1); % 简单分布为标准正态分布
buf=exppdf(data_next,10)/exppdf(data(index),10);
alpha=min(1,buf);
if u<alpha
data=[data data_next];
else
data=[data data(index)];
end
index = index+1;
end
subplot(1,2,2)
histogram(data,'BinWidth',1);
axis([0,100,0,data_size/10]);
结果：
