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(別名采樣)