Alias method/别名采样方法证明和实现

最新推荐文章于 2022-05-13 21:29:49 发布

原创

最新推荐文章于 2022-05-13 21:29:49 发布 · 1.5k 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#采样

本文介绍了Alias method，一种时间复杂度为O(1)的离散采样算法。通过数学证明，展示了该算法如何将不同概率分布的样本分配到体积为1的矩形容器中，每个容器最多使用两个概率项。采用归纳法，当样本数量较少时算法的正确性显而易见，随着样本数量增加，通过调整概率分布确保算法的正确执行。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

算法介绍：【数学】时间复杂度O(1)的离散采样算法—— Alias method/别名采样方法

先上代码


class alias_sampling(object):
    def __init__(self,prob):
        self.n = len(prob)
        self.area_ratio = prob * self.n
        self.__create_alias_table()

    def __create_alias_table(self):
        low_indices,hight_indices = [],[]
        self.accept,self.alias = [0] * self.n,[0] * self.n
        for i,ratio in enumerate(self.area_ratio):
            if ratio < 1.0:
                low_indices.append(i)
            else:
                hight_indices.append(i)

        while low_indices and hight_indices:
            low_id,hight_id = low_indices.pop(),hight_indices[-1]
            self.accept[low_id] = self.area_ratio[low_id]
            self.alias[low_id] = hight_id
            self.area_ratio[hight_id] -= (1 - self.area_ratio[low_id])

            if self.area_ratio[hight_id] < 1.0:
                hight_indices.pop()
                low_indices.append(hight_id)

        for hight_id in hight_indices:
            self.accept[hight_id] = 1.
        for low_id in low_indices:
            self.accept[low_id] = 1.
        del self.area_ratio

    def get_sample(self):
        id = np.random.randint(0,self.n)
        if self.accept[id] > np.random.uniform(0,1):
            return id
        return self.alias[id]

    def get_samples(self,n):
        for i in range(n):
            yield self.get_sample()

def test(n=100,k=1000000):

    def gen_prop_dist(n):
        p = np.random.uniform(0,10,n)
        return p / np.sum(p)

    p = gen_prop_dist(n)
    sampling = alias_sampling(p)
    p_sample = np.zeros(n)

    for sample in sampling.get_samples(k):
        p_sample[sample] += 1

    p_sample /= k

    print(np.linalg.norm(p-p_sample,2))


if __name__ == '__main__':
    for i in range(10):
        test()

假设概率分布为 ${p_i\}_{i=0}^{n-1}$