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)的采样算法---Alias

(1)alias分为两步

  • 做表
  • 根据表采样

(2)做表

将概率分布的每个概率乘上N,画出柱状图。N为事件数量。

 

其总面积为N,可以看出某些位置面积大于1某些位置的面积小于1,

将面积大于1的事件多出的面积补充到面积小于1对应的事件中,以确保每一个小方格的面积为1,同时,保证每一方格至多存储两个事件,这样我们就能看到一个1*N的矩形啦。

表里面有两个数组,一个数组存储的是事件i占第i列矩形的面积的比例,即

另一个数组存储第i列中不是事件i的另一个事件的编号,即

做表的时间复杂度是o(N)。

(3)采样

产生两个随机数,

第一个产生一个1~N 之间的整数 i,决定落在哪一列。

第二个产生一个0~1之间的随机数,判断其与Prab[i]大小,

如果小于Prab[i],则采样 i,(表示若其小于事件i占第i列矩形的面积的比例,则表示接受事件i

如果大于Prab[i],则采样Alias[i](否则,接收第i列中不是事件i的另一个事件。)

这种方式采样的时间复杂度为 o(1) ,且对应事件i的概率,完全对应原来的概率分布。

计算:选择第一列的概率为1/4,则选择第一列的蓝色部分的概率为1/4*2/3 = 1/6,蓝色部分还有第三、第四列,则蓝色(事件A)的总概率为1/4 * 2/3 * 3 = 1/2。【第一次产生随机数的概率为1/4,第二次产生随机数选择蓝色部分的概率为1/4 *(2/3*3) = 1/2】。--------对应原来的概率分布。

四、代码

1、制作表:(create_alias_table函数,O(N))

(1)概率分布 area_ratio 的每个概率乘上N

(2)small,large分别存放小于1和大于1的事件下标。

(3)每次从small,large中各取一个,将大的补充到小的之中,小的出队列,再看大的减去补给之后,如果大于1,继续放入large中,如果等于1,则也出去,如果小于1则放入small中。

 

(4)获取accept、alias

accept存放第i列对应的事件i矩形的面积百分比;

alias存放第i列不是事件i的另外一个事件的标号;

上例中accept=[2/3,1,1/3,1/3],alias=[2,0,1,1],这里alias[1]的0是默认值,也可默认置为-1避免和事件0冲突;

2、采样:(alias_sample函数,O(1))

随机采样1~N 之间的整数i,决定落在哪一列。

随机采样0~1之间的一个概率值,

如果小于accept[i],则采样i,

如果大于accept[i],则采样alias[i];

import numpy as np
def create_alias_table(area_ratio):
    """
    area_ratio[i]代表事件i出现的概率
    :param area_ratio: sum(area_ratio)=1
    :return: accept,alias
    """
    N = len(area_ratio)
    accept, alias = [0] * N, [0] * N
    small, large = [], []
  ###------(1)概率 * N ----- area_ratio_
= np.array(area_ratio) * N

###------(2)获取small 、large -----
for i, prob in enumerate(area_ratio_): if prob < 1.0: small.append(i) else: large.append(i)
  ###------(3)修改柱状图 ----- (4)获取accept和alias -----
while small and large: small_idx, large_idx = small.pop(), large.pop() accept[small_idx] = area_ratio_[small_idx] alias[small_idx] = large_idx area_ratio_[large_idx] = area_ratio_[large_idx] - \ (1 - area_ratio_[small_idx]) if area_ratio_[large_idx] < 1.0: small.append(large_idx) else: large.append(large_idx)
while large: large_idx = large.pop() accept[large_idx] = 1 while small: small_idx = small.pop() accept[small_idx] = 1 return accept, alias def alias_sample(accept, alias): """ :param accept: :param alias: :return: sample index """ N = len(accept) i = int(np.random.random()*N) r = np.random.random() if r < accept[i]: return i else: return alias[i]

 

 

参考文献:

Alias method:时间复杂度O(1)的离散采样方法

时间复杂度O(1)的离散采样算法—— Alias method/别名采样方法

Alias sample(别名采样)

posted on 2020-04-22 00:01  吱吱了了  阅读(4484)  评论(0编辑  收藏  举报

导航