Alias 采样

参考：【数学】时间复杂度O(1)的离散采样算法—— Alias method/别名采样方法
原文：关于数学理论

对于问题：
比如一个随机事件包含四种情况，每种情况发生的概率分别为： $\frac{1}{2},\frac{1}{3},\frac{1}{12},\frac{1}{12}$，问怎么用产生符合这个概率的采样方法。

两次掷骰子：第一次决定所属列，第二次决定事件概率。

算法流程

1. 事件N中情况概率按照均值归一化，即每种情况概率乘N。
   

2. 构造Alias表。
   Alias Method一定要保证：每列中最多只放两个事件；
   

3. 按照某种方法将整个概率分布拉平成为一个1*N的长方形即为Alias Table，然后储存两个数组，一个里面存着第i列对应的事件i矩形所占的面积百分比，即概率，上图的话数组就为 $Prob[\frac{2}{3}, 1, \frac{1}{3}, \frac{1}{3}]$。另一个数组里面储存着第i列不是事件i的另外一个事件的标号，上图就是Alias[2，NULL，1，1]。

Alias Table构造方法

用两个队列A,B去储存面积大于1的节点标号，和小于1的节点标号，每次从两个队列中各取一个，将大的补充到小的之中，小的出队列，再看大的减去补给之后，如果大于1，继续放入A中，如果等于1，则也出去，如果小于1则放入B中。算法复杂度为O(n)。

python代码

class AliasSampling:
    """Helper function for performing sampling."""
    def __init__(self, prob):
        self.n = len(prob)
        self.U = np.array(prob) * self.n
        self.K = [i for i in range(len(prob))]
        overfull, underfull = [], []
        for i, U_i in enumerate(self.U):
            if U_i > 1:
                overfull.append(i)
            elif U_i < 1:
                underfull.append(i)
        while len(overfull) and len(underfull):
            i, j = overfull.pop(), underfull.pop()
            self.K[j] = i
            self.U[i] = self.U[i] - (1 - self.U[j])
            if self.U[i] > 1:
                overfull.append(i)
            elif self.U[i] < 1:
                underfull.append(i)

    def sampling(self, n=1):
        x = np.random.rand(n)
        i = np.floor(self.n * x)
        y = self.n * x - i
        i = i.astype(np.int32)
        res = [i[k] if y[k] < self.U[i[k]] else self.K[i[k]] for k in range(n)]
        if n == 1:
            return res[0]
        else:
            return res