Python Alias采样算法实现

Alias采样

离散分布：给你一个概率分布，是离散的，比如[1/2, 1/3, 1/12, 1/12]，代表某个变量属于事件A的概率为1/2， 属于事件B的概率为1/3，属于事件C的概率为1/12，属于事件D的概率为1/12。

离散分布的随机变量的取样问题：

一个随机事件包含四种情况，每种情况发生的概率分别为： 1/2, 1/3, 1/12, 1/12，问怎么用产生符合这个概率的采样方法。

时间复杂度为o(N)的采样算法

首先将其概率分布按其概率对应到线段上，像下图这样：

产生0~1之间的一个随机数，然后看起对应到线段的的哪一段，就产生一个采样事件。

比如落在0~ 1/2之间就是事件A，落在1/2~5/6之间就是事件B，落在5/6~11/12之间就是事件C，落在11/12~1之间就是事件D。

构建线段的时间复杂度为o(N)：

如果顺序查找线段的话，查找时间复杂度为O(N)，

如果二分查找的话，查找的时间复杂度为O(logN)。

时间复杂度O(1)的采样算法

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

2）构造Alias表

PS：Alias Method一定要保证：每列中最多只放两个事件；

下标:[0,1,2,3]

3）按照某种方法将整个概率分布拉平成为一个1*N的长方形即为Alias Table

储存两个数组:

1.存着第i列对应的事件i矩形所占的面积百分比 即概率，上图的话数组就为 Prob 。

2.存着第i列不是事件i的另外一个事件的标号，上图就是Alias [1, -1, 0, 0]。（-1即自身满足）

代码实现

# -*- coding: utf-8 -*-import numpy as np # 创建alias_tabledef create_alias_table(prob_val):    L = len(prob_val)    # 初始化两个数组    alias_prob = np.zeros(L)    # 储存概率    events_index = np.ones(L,dtype="int64") * -1 # 储存 下标/序号,-1表示自身足够      # 大的队列用于存储面积大于1的节点标号，小的队列存储面值小于1的节点标号    small_queue = []    large_queue = []     # 把 prob_val 均值归一化存储 并 把下标放到对应大/小队列    for index,prob in enumerate(prob_val):        alias_prob[index] = L * prob         if alias_prob[index] < 1.0:            small_queue.append(index)        else:            large_queue.append(index)     # 1. 每次从两个队列 各取一个，让 大的去补充小的，然后小的出small队列    # 2. 在看大的减去补给小的之后剩余的值    # 如果大于1，继续放到large队列；    # 如果恰好等于1，也出队列；    # 如果小于1加入small队列中；    while small_queue and large_queue:        small_index = small_queue.pop()        large_index = large_queue.pop()         # 因为 alias_index 中存的：另一个事件的标号，        # 那现在用大的概率补充小的概率，标号就要变成大的的事件的标号        events_index[small_index] = large_index        # 补充的原则是：大的概率要把小的概率 补满（补到概率为1），然后就是剩下的        alias_prob[large_index] = alias_prob[large_index] + alias_prob[small_index] - 1.0         # 判断补完后，剩余值得大小        if alias_prob[large_index] < 1.0:            small_queue.append(large_index)        elif alias_prob[large_index] > 1.0:            large_queue.append(large_index)     return alias_prob, events_index   if __name__ == '__main__':    alias_prob,events_index = create_alias_table([1/2, 1/3, 1/12, 1/12])    print(alias_prob,events_index)  # alias 采样def alias_smaple(alias_prob, events_index):    N = len(alias_prob)     # 第一个骰子，产生第一个1~N的随机数，决定落在哪一列    random_num1 = int(np.floor(np.random.rand()*N))    # 第二个骰子，产生-0~1之间的随机数，判断与accept_prob[random_num1]的大小    random_num2 = np.random.rand()     # 如果小于Prab[i]，则采样i，如果大于Prab[i]，则采样Alias[i]    print(random_num1)    print(random_num2,alias_prob[random_num1])    if random_num2 < alias_prob[random_num1]:        return random_num1    else:        return events_index[random_num1] if __name__ == '__main__':    alias_prob, events_index = create_alias_table([1/2, 1/3, 1/12, 1/12])    val = alias_smaple(alias_prob,events_index)    print(val)
免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删