Alias method 别名算法

以下内容来自转载:

因为需要用到Alias Sampling Method的方法,但是查了一下,发现没有找到靠谱的关于Alias Method的中文介绍,所以干脆自己写一个好了。 
关于Alias Method的介绍的比较好的是一个外国Blog:Darts, Dice, and Coins: Sampling from a Discrete Distribution,以下的介绍也主要参考这篇Blog里的算法。

问题:比如一个随机事件包含四种情况,每种情况发生的概率分别为: 12,13,112,112,问怎么用产生符合这个概率的采样方法。

最容易想到的方法

我之前有在【数学】均匀分布生成其他分布的方法中写过均匀分布生成其他分布的方法,这种方法就是产生0~1之间的一个随机数,然后看起对应到这个分布的CDF中的哪一个,就是产生的一个采样。比如落在0~ 12之间就是事件A,落在12~56之间就是事件B,落在56~1112之间就是事件C,落在1112~1之间就是事件D。 
但是这样的复杂度,如果用BST树来构造上面这个的话,时间复杂度为O(logN),有没有时间复杂度更低的方法。 

一个Naive的办法

一个Naive的想法如下:  
Alias method 别名算法 
1. 可以像上图这样采样,将四个事件排成4列:1~4,扔两次骰子,第一次扔1~4之间的整数,决定落在哪一列。 
Alias method 别名算法 
2. 如上如所示,将其按照最大的那个概率进行归一化。在1步中决定好哪一列了之后,扔第二次骰子,0~1之间的任意数,如果落在了第一列上,不论第二次扔几,都采样时间A,如果落在第二列上,第二次扔超过23则采样失败,重新采样,如果小于23则采样时间B成功,以此类推。 
3. 这样算法复杂度最好为O(1)最坏有可能无穷次,平均意义上需要采样O(N)

那怎么去改进呢?

Alias Method

这样Alias Method采样方法就横空出世了 
Alias method 别名算法 
还是如上面的那样的思路,但是如果我们不按照其中最大的值去归一化,而是按照其均值归一化。即按照1N(这里是$\frac{1}{4})归一化,即为所有概率乘以N,得到如下图: 
Alias method 别名算法 
其总面积为N,然后可以将其分成一个1*N的长方形,如下图:  
Alias method 别名算法 
将前两个多出来的部分补到后面两个缺失的部分中。 
先将1中的部分补充到4中: 
Alias method 别名算法 
这时如果,将1,2中多出来的部分,补充到3中,就麻烦了,因为又陷入到如果从中采样超过2个以上的事件这个问题中,所以Alias Method一定要保证:每列中最多只放两个事件 
所以此时需要讲1中的补满3中去: 
Alias method 别名算法 
再将2中的补到1中: 
Alias method 别名算法 
至此整个方法大功告成 
Alias Method具体算法如下: 
1. 按照上面说的方法,将整个概率分布拉平成为一个1*N的长方形即为Alias Table,构建上面那张图之后,储存两个数组,一个里面存着第i列对应的事件i矩形站的面积百分比【也即其概率】,上图的话数组就为Prab[2311313],另一个数组里面储存着第i列不是事件i的另外一个事件的标号,像上图就是Alias[2 NULL 1 1] 
2.产生两个随机数,第一个产生1~N 之间的整数i,决定落在哪一列。扔第二次骰子,0~1之间的任意数,判断其与Prab[i]大小,如果小于Prab[i],则采样i,如果大于Prab[i],则采样Alias[i]

这个算法是不是非常的精妙而且简洁,做到了O(1)事件复杂度的采样。但是还有一个问题是,如何去构建上面的第1步?即如何去拉平整个概率分布?这样预处理的时间复杂度又是多少呢?

Alias Method程序化构建

Naive方法

构建方法: 
1.找出其中面积小于等于1的列,如i列,这些列说明其一定要被别的事件矩形填上,所以在Prab[i]中填上其面积 
2.然后从面积大于1的列中,选出一个,比如j列,用它将第i列填满,然后Alias[i] = j,第j列面积减去填充用掉的面积。

以上两个步骤一直循环,直到所有列的面积都为1了为止。

存在性证明

那么Alias Table一定存在吗,如何去证明呢? 
要证明Alias Table一定存在,就说明上述的算法是能够一直运行下去,直到所有列的面积都为1了为止,而不是会中间卡住。 
一个直觉就是,这一定是可以一直运行下去的。上述方法每运行一轮,就会使得剩下的没有匹配的总面积减去1,在第n轮,剩下的面积为N-n,如果存在有小于1的面积,则一定存在大于1的面积,则一定可以用大于1的面积那部分把小于1部分给填充到1,这样就进入到了第n+1轮,最后一直到所有列面积都为1。 
更为严谨的证明见上面给出的那个Blog。

更快的构建方法

如果按照上面的方法去构建Alias Table,算法复杂度是O(n2)的,因为最多需要跑N轮,而每跑一轮的最大都需要遍历N次。一个更好的办法是用两个队列A,B去储存面积大于1的节点标号,和小于1的节点标号,每次从两个队列中各取一个,将大的补充到小的之中,小的出队列,再看大的减去补给之后,如果大于1,继续放入A中,如果等于1,则也出去,如果小于1则放入B中。 
这样算法复杂度为O(n)

至此Alias Method就讲完了,感觉还是一个非常精妙的方法,而且方法实现起来也非常的简单。值得学习。