1.解决的问题
计算机怎么在任意给定的概率分布P上采样?首先可以想到把它拆成两步:
(1)首先等概率的从采样区间里取一个待定样本x,并得到它的概率为p(x)
(2)然后在均匀分布U[0,1]上取一个值,如果这个值高于了样本出现的概率p(x),那么就舍弃掉这个待定样本,如果低于或等于这个值就保留。
不断重复这个步骤,我们就能得到满足概率分布P的样本X。
这个方法看似有效,但存在一个致命的问题,那就是慢,因为概率p的均值,随着采样的区域增大而减小【如果P是多维概率分布,那么随着维度不断增大,那么这个值将趋近于0,也就是维度灾难】,那么就意味着这个样本被保留下来的概率非常小。
那么怎么改进?
我们把上述的第(2)步改一下,不再用均匀U[0,1]进行采样,而是用U[0,max(p)],而选取的规则变为,如果从均匀分布U[0,max(p)]里取得的值u高于p(x),那么就保留该样本。
这样能够保证,至少有一个可能出现的样本在概率区间接受率为1。
这种改进比第一种的命中率要高上不少,如果p是均匀分布,那么采样命中率将是100%,当然如果采样是均匀分布,那么我直接取就好了,没必要多这么一步。
那么在此之上还有什么改进呢?
你自然能想到,我们不再用均匀分布U来计算样本是否可以保留,我们找一个好实现的概率分布Q,这个Q要和分布P相似,我们放缩该概率密度函数,使得min(q)>max(p)。
这样一来,命中率就又提高了,而且可以预想到,如果概率分布q=p那么采样命中率就会是100%,但这自然不可能。
那么有没有一种方法可以让采样的命中率是100%呢?
有,就是MCMC,在MCMC中,所有样本都是有效采样。而原理则是用频率取代概率。
2.MCMC
假设我们手里有个小球,每隔一秒就会变一种颜色(包括当前颜色),一段时间后,我们统计每种颜色的时间,就可以得到,各颜色出现的概率。
而这就是MCMC的原理,我们要做的是借助马尔可夫链,实现这样一个小球。
我们给小球输入一个指令,给定它每种颜色出现的概率,也就是P。
这时候小球内部有个状态机会来实现这个概率,首先会按照颜色的数量建立状态,同时每种状态之间会有线路连接,且每种状态之间都是相互可达的,即连通。
这时候有钢珠在每种状态里游走,钢珠到了哪种状态,哪种状态代表的灯就会亮起。
好了,如何让小球的灯能够以某种稳定的频率让灯亮起。
研究发现,只要让两两状态之间转换的概率与对方状态的概率乘积相等。
即P(A)Q(A->B)=P(B)Q(B->A),这就是马尔科夫链细致平衡条件,A和B代表任意两个状态。Q(A-B)是状态A到状态B的转移概率。
这时候,我们只要跟随钢珠在各状态之间的游走,记录下状态值即样本值,就能完成采样。
具体步骤如下,
(1)记录钢珠初始值状态X1。
(2)随机选择一个目标状态。
(3)查看抵达目标状态的转移概率。
(4)从U[0,1]中采一个值,如果值大于转移概率,那么状态不变,如果小于或等于转移概率,那么状态转化为目标状态。
重复如上步骤,就能得到一个服从P分布的样本了。
但这个算法虽然没有无效采样,但状态A转移到状态的转移概率过低,那么就会长时间卡在这一种状态。那么有没有方法改进呢。
当然有,那就是等比例缩放平衡方程两边的转移概率,P(A)Q(A->B)=P(B)Q(B->A)。
这里有个限制,那就是不能破坏转移概率小于或等于1概率特性,同时也不能破坏平衡方程本身。
那么我们能找到的一个缩放系数就是MAX[P(B->A),P(A->B)],两边同时除以该系数。
则有P(A)MIN[1,Q(A->B)/Q(B->A)] = P(B)MIN[1,Q(A->B)/Q(B->A)]
这样就能保证至少有一边可以必定转化为另一边。
到这里MCMC就算比较完整了。
问题就是Q怎么选。
这里有两种方案,一个就是Q任取,只要满足转移矩阵的要求,然后补一个系数配平方程
即: P(A)Q(A->B)*[P(B)Q(B->A)]=P(B)Q(B->A)*[P(A)Q(A->B)],
缩放后有:
P(A)Q(B->A)*MIN[1,P(B)Q(B->A)/P(B)Q(B->A)] = P(B)Q(B->A)*MIN[1,P(A)Q(A->B)/P(B)Q(B->A)]
另一种那就是选Q(A->B) = P(B),Q(B->A)=P(A)。
前者是HM的思想,配平的系数的使用方法类比前一节的转移接受率。
后者是Gibbs的思想。
具体实现有很多细节,可以看下面推荐的教学视频,这里推荐先看P6
机器学习-白板推导系列-蒙特卡洛方法5(Monte Carlo Method)- 再回首_哔哩哔哩_bilibili