关于MCMC方法状态转移的问题

这是一个简短的文章，只讨论一个问题：在MCMC方法中，状态转移失败时，为什么接受上一个状态，而不是直接重新采样。
在探讨这个问题之前，你应该已经比较了解MCMC，如若不然，请移步刘建平大佬的博客看看MCMC的前世今生，地址如下：
MCMC详细讲解

刘大佬的文章是很详细的，唯有一点令人存疑。以M-H采样为例，它的步骤如下：

1）输入我们任意选定的马尔科夫链状态转移矩阵Q，平稳分布π(x)，设定状态转移次数阈值n1，需要的样本个数n2
2）从任意简单概率分布采样得到初始状态值x0
3）for t=0 to n1+n2−1:

a) 从条件概率分布Q(x|xt)中采样得到样本x∗
b) 从均匀分布采样u∼uniform[0,1]
c) 如果u<α(xt,x∗)=min{π(j)Q(j,i)/π(i)Q(i,j),1}, 则接受转移xt→x∗，即xt+1=x∗
d) 否则不接受转移，即xt+1=xt

在第（3）步的d中，如果不接受转移，则xt+1 = xt。
看到此处，我得小脑瓜嗡嗡的。为什么不是重新采样xt+1呢，而是让xt+1=xt？这合理吗？这是不是等于将xt重复了好多遍，这还是按照π的概率在采样吗？

关于这个问题，刘大佬也在回帖中进行了回复，乍一看，好像有那么点意思，但实际又并不是那么回事，直接贴出大佬的回复，这里不再评论。

下面给出我的分析。先说结论：采用xt+1 = xt可以，重新采样xt+1也可以，二者都行，但是前者的采样效率更高。
那么为什么呢？
先看公式，细致平稳的公式如下：

π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i)
α(i,j)=min{π(j)Q(j,i)/π(i)Q(i,j),1}

说人话就是， i 转移到 j 的概率等于 j 转移到 i 的概率。
当转移失败时，新的样本就等于原样本，那么这个概率是多少呢，更进一步， i 转移到 i 的概率是多少呢？是 Q(i,i)*a(i,i)再加上所有“转移失败的概率”之和。但是。但是，不管这个概率是多少， i 转移到 i 的概率永远等于 i 转移到 i 的概率。这句话没有错别字，就是 i 到 i 等于 i 到 i ，这不就是细致平稳的条件吗？所有它一定满足平稳的要求啊，会收敛到π分布。
那么现在破案了，MCMC的转移策略是没错的，但是明明一个样本老是double自己，它怎么还能是没错呢？我觉着可以这样理解：因为是所有样本都在double，并且对于在π中概率高（似然值高）的样本，别人转移到它的概率高，它转移出去的概率低（double几率高），在π中概率低的样本则相反。感性的认识就是，它们在按照π的分布double，概率高double多，概率低double少。妙蛙
那么不double行吗，也可以的，只不过会经常采样失败，样本数量的增加十分缓慢，半天才增加一点点（谭警官，只能拉一点点）。更细致一点（后面的说话有点随意，随缘理解吧，写公式太麻烦了），还可以看Q(i,j)*a(i,j)，Q是概率密度，a处处小于1，乘在一起积分一定小于1，所以大量miss是正常的，而采样失败则直接采纳上一个样本，恰恰弥补了积分不足1的那一部分。妙蛙
所以xt+1 = xt这一策略，显著提高了采样率。那么到底是不是这回事呢，实践是检验真理的唯一标准，直接跑代码。
先上图，后上代码。下图是用MCMC模拟采样，平稳分布π是正态分布，其中蓝色为标准的π分布，红色是抽到的样本的分布。在转移失败时，第一个图是采取xt+1 = xt的策略，第二个是直接对xt+1重新采样。可以看出，两个策略基本都很好地模仿了正态分布（当然看起来好像第一个更好，但是理论上应该是一样好才对）。

以下是代码，爱看不看，反正我不爱看。

import random
import math
from scipy.stats import norm
import matplotlib.pyplot as plt


def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=3, scale=2)
    return y

T = 50000
pi = [0 for i in range(T)]
sigma = 1
t = 0
while t < T-1:
    t = t + 1
    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None)
    alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))

    u = random.uniform(0, 1)
    if u < alpha:
        pi[t] = pi_star[0]
    else:
        # 两种策略的差别就在这两行
        
        # pi[t] = pi[t - 1]     # xt+1 = xt
        t -= 1                  # 重新采样xt+1


plt.scatter(pi, norm.pdf(pi, loc=3, scale=2))
num_bins = 50
plt.hist(pi, num_bins, normed=1, facecolor='red', alpha=0.7)
plt.show()

最后，其实还是有点心里不稳，如果有错，欢迎大佬指正。谢过。

（如需转载，请注明出处）