首页 > 其他 > 详细

Alias sample(别名采样)

时间:2019-09-28 11:22:52      阅读:221      评论:0      收藏:0      [点我收藏+]

应用场景:加权采样,即按照随机事件出现的概率抽样

具体算法:

技术分享图片

举例如上,随机事件出现的概率依次是1/2,1/3,1/12,1/12;记随机事件的个数为N,则所有事件概率乘以N后概率为2,4/3,1/3,1/3;

记队列small,large分别存放小于1和大于1的事件下标(例子中small=[0,1],large=[2,3]);

记accept存放第i列对应的事件i矩形的面积百分比;alias存放第i列不是事件i的另外一个事件的标号;

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

上例中accept=[2/3,1,1/3,1/3],alias=[1,0,1,1],这里alias[1]的0是默认值,也可默认置为-1避免和事件0冲突;以上部分在源码create_alias_table函数中,时间复杂度是O(N);

随机采样1~N 之间的整数i,决定落在哪一列。随机采样0~1之间的一个概率值,如果小于accept[i],则采样i,如果大于accept[i],则采样alias[i];这一部分源码在alias_sample;

下面附上具体的源码:

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 = [], []
    area_ratio_ = np.array(area_ratio) * N
    for i, prob in enumerate(area_ratio_):
        if prob < 1.0:
            small.append(i)
        else:
            large.append(i)

    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 sample(别名采样)

原文:https://www.cnblogs.com/arachis/p/alias_sample.html

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!