关于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()
最后,其实还是有点心里不稳,如果有错,欢迎大佬指正。谢过。
(如需转载,请注明出处)