文章目錄

• 離散采樣Sampling
• • 問題定義
  • O(N)的方法
  • O(logN)的方法
  • O(1)的方法
  • Python代碼
  • 參考

離散采樣Sampling

問題定義

給定一個離散型随機變量的機率分布規律 P ( X = i ) = p i , i ∈ 1 , . . . N P(X=i)=p_i,i\in{1,...N} P(X=i)=pi​,i∈1,...N，希望設計一個方法能夠從該機率分布中進行采樣使得采樣結果盡可能服從機率分布P

O(N)的方法

想象随機事件依其機率的大小分布在一個長度為1的線段上。那麼我線上段中随機取一點，看看該點落在哪個事件對應的區間中，就取該區間對應的事件即可。

具體實作如下：

顯然線構造一個存儲累積事件機率的數組時間複雜度為O(N),每次進行線性查找的時間複雜度為O(N)

O(logN)的方法

上面在通過累加的方式構造出累積事件機率數組後，我們可以發現該數組滿足非遞減有序，對于有序數組的查找很容易想到使用二分查找進行優化，使用二分查找後的時間複雜度為 O ( l o g N ) O(logN) O(logN)

O(1)的方法

Alias方法整個機率分布轉化為一個1*N的矩形，對于每個事件i，轉換為對應矩形中的面積為 p i ∗ N ∑ k = 1 k = N p k \frac{p_i*N}{\sum_{k=1}^{k=N}p_k} ∑k=1k=N​pk​pi​∗N​

通過上述操作，一般會有某些位置面積大于1某些位置的面積小于1。我們通過将面積大于1的事件多出的面積補充到面積小于1對應的事件中，以確定每一個小方格的面積為1，同時，保證每一方格至多存儲兩個事件。

維護兩個數組accept和alias，accept數組中的accept[i]表示事件i占第i列矩形的面積的比例。

alias[i]表示第i列中不是事件i的另一個事件的編号

在進行采樣的時候，每次生成一個随機數 i ∈ [ 0 , N ) i∈[0,N) i∈[0,N),再生成一個随機數 r ∼ U n i f ( 0 , 1 ) r∼Unif(0,1) r∼Unif(0,1)，若 r < a c c e p t [ i ] r <accept[i] r<accept[i]，則表示接受事件i，否則，拒絕事件i傳回 a l i a s [ i ] alias[i] alias[i]

該算法對應的代碼如下，可以看到預處理alias table的時間複雜度仍為O(N)，而每次采樣産生時間的時間複雜度為O(1)。

Python代碼

import numpy as np
def gen_prob_dist(N):
    p = np.random.randint(0,100,N)
    return p/np.sum(p)



def create_alias_table(area_ratio):

    l = len(area_ratio)
    accept, alias = [0] * l, [0] * l
    small, large = [], []

    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):
    N = len(accept)
    i = int(np.random.random()*N)
    r = np.random.random()
    if r < accept[i]:
        return i
    else:
        return alias[i]

    

def simulate(N=100,k=10000,):
    

    truth = gen_prob_dist(N)
    
    area_ratio = truth*N
    accept,alias = create_alias_table(area_ratio)
    
    ans = np.zeros(N)
    for _ in range(k):
        i = alias_sample(accept,alias)
        ans[i] += 1
    return ans/np.sum(ans),truth

if __name__ == "__main__":
    alias_result,truth = simulate()

           

參考

1. Alias Method離散分布随機取樣
2. Darts, Dice, and Coins: Sampling from a Discrete Distribution
3. 【數學】時間複雜度O(1)的離散采樣算法—— Alias method/别名采樣方法
4. Alias:時間複雜度O(1)的離散采樣方法