Metropolis-Hasting 算法 & 图上的Metropolis-Hasting Random Walk (MHRW)

背景介绍：

在概率与统计领域，往往需要根据已知的概率分布生成相应的采样集合，从而支持后续的统计分析与研究。有时候已知的概率分布比较简单，很容易生成一系列的采样集合，而有时候已获得的概率分布较为复杂，很难直接生成相应的采样集合。那么是否存在一种算法，不管给定的分布复杂与否，均可以统一的生成相应的采样集合呢？答案是肯定的。今天我们就将介绍一种普遍使用的算法：Metropolis-Hasting (MH) 算法。

MH算法名字中的Metropolis和Hasting分别代表两个人。其中，Nicholas Metropolis是1953年发表的论文Equation of State Calculations by Fast Computing Machines 中的第一作者。在这篇论文中为了讨论对称分布，首先提出了该算法。而在1970年，W. K. Hastings 将该算法扩展到一般情况，从而奠定了它的普适性。

Metropolis-Hasting 算法：

MH算法是一个基于Markov Chain Monte Carlo(MCMC) 的近似算法。它实现的原理是首先根据给定的任意一个概率分布，构造一个以该分布为静态分布的Markov Chain，然后执行该Markov Chain 到达收敛之后（每个点被访问的概率服从静态分布）开始采样，此时获得的采样集合近似认为是服从给定分布 π 的。

从上述的原理描述，我们不难看出，算法所面临的挑战在于如何通过给定的分布 π 来设计相应的Markov Chain？ 因为每一条Markov Chain都对应唯一的一个转移矩阵，所以该问题等价于如何根据 π 确定转移概率矩阵。

为了解决上述挑战，我们首先需要搞清楚什么样的转移概率矩阵对应的Markov Chain才会存在静态分布，以及是否该静态分布是唯一的？

1. 静态分布存在的一个充分不必要条件: 给定一条Markov Chain {Xt}t≥0，设其转移矩阵为P={p(i,j)}i,j∈S，其中S为状态空间，p(i,j)表示从转态i到状态j的一步转移概率。若存在一个分布π满足：p(i,j)π(i)=p(j,i)π(j)， 则π是该Markov Chain的一个静态分布。同时称该链是time-reversible的。

2. 静态分布唯一存在的充分必要条件：若一条Markov Chain 是非周期（aperiodic）和常回归的（positive recurrent）的，等价于该链是可遍历的（ergodicity），则可以证明该链存在唯一的静态分布，

上面的两个结论为我们如何根据分布π 设计转移矩阵指明了方向。换句话说，若我们能设计出满足上面两个条件的转移概率，其对应的Markov Chain一定存在唯一的静态分布，且该分布就是π。通过利用第一个结论，我们可以得到1式：

p(i,j)p(j,i)=π(j)π(i)

为了满足第二个结论，我们可以将每一步的转移概率分成两步：选举概率Q(i,j)和接受概率A(i,j)。即每一次转移时，首先从i的相邻状态中随机选择一个状态j 作为候选下一跳，然后以概率A(i,j)接受它作为真正的下一跳。因此，我们有2式：

p(i,j)=Q(i,j)A(i,j)

结合1式和2式，我们可以得到3式：

Q(i,j)A(i,j)Q(j,i)A(j,i)=π(j)π(i)

进而，我们可以得到4式：

A(i,j)=π(j)Q(j,i)π(i)Q(i,j)∗A(j,i)

记 T(j,i)=π(j)Q(j,i)π(i)Q(i,j)

从而我们有5式：A(i,j)=T(j,i)A(j,i)

结合T(j,i)和5式，不难得出，T(i,j)=1/T(j,i)和A(i,j)=F(T(i,j))，其中F满足F(x)=F(1/x)/x。由于F不唯一，我们可以选择其中一个解 F(x)=min{1,x}。最终，我们有：

A(i,j)=min{1,T(i,j)}=min{1,π(j)Q(j,i)π(i)Q(i,j)}

到这里，我们便根据给定的分布π，完成了对转移概率的设计，即对任意的 j≠i：

p(i,j)=Q(i,j)∗min{1,π(j)Q(j,i)π(i)Q(i,j)}

其中，p(i,i)=1−p(i,j)。

图上的Metropolis-Hasting Random Walk (MHRW)

作为一个具体的例子，接下来，我们将介绍在给定图上的基于MH算法的随机游走。

给定一个无向连通图G(V,E)，其中V是结点集合，E是边集合。对任意一个结点i，记N(i)为i的邻居集合，d(i)为邻居个数。MHRW算法假定选举概率 Q(i,j)=1/d(i)，即随机均匀的从i的邻居中选择一个。若假定给定的分布是均匀分布，即π(i)=1/|V|。则我们有对任意的 j≠i：

p(i,j)=1d(i)∗min{1,d(i)d(j)}

具体的算法伪代码如下：

v ← initial node
while stopping criterion not met do
Select node w uniformly at random from neighbors of v.
Generate uniformly at random a number 0≤p≤1.
if p < d(v)/d(w) then
v ← w.
else
Stay at v
end if
end while

参考文献：

[1]:Metropolis–Hastings algorithm -Wiki

[2]:Beyond Random Walk and Metropolis-Hastings Samplers: Why You Should Not Backtrack for Unbiased Graph Sampling

[3]: Walking in Facebook:A Case Study of Unbiased Sampling of OSNs

[4]:Equation of State Calculations by Fast Computing Machines