2. 程式人生 > >加權隨機取樣 (Weighted Random Sampling)

加權隨機取樣 (Weighted Random Sampling)

阿新 • • 發佈：2018-12-21

一個集合裡有 n 個元素，每個元素有不同的權重，現在要不放回地隨機抽取 m 個元素，每個元素被抽中的概率為元素的權重佔總權重的比例。要怎麼做呢？

簡單的解法

現在考慮只抽取一個元素，假設權重之和為 1。我們可以從 [0, 1] 中隨機得到一個權重，假設為 0.71，而後從第一個元素開始，不斷累加它們的權重，直到有一個元素的累加權重包含 0.71，則選取該元素。下面是個示意圖：

要選取 m 個元素，則可以按上面的方法先選取一個，將該元素從集合中去除，再反覆按上面的方法抽取剩餘的元素。這種方法的複雜度是 O(mn)，並且將元素從集合中刪除其實不太方便實現。

當然，最重要的是這個演算法需要多次遍歷資料，不適合用在流處理的場景中。

Algorithm A

  1. 對於集合 $V$ 中的元素 $v_i \in V$，選取均勻分佈的隨機數 $u_i = rand(0, 1)$ ，計算元素的特徵 $k_i = u_i^{(1/w_i)}$
  2. 將集合按 $k_i$ 排序，選取前 $m$ 大的元素。

演算法的正確性在作者 2006 年的論文 Weighted random sampling with a reservoir 裡給了詳細的證明。論文中給出了演算法的兩個變種 A-Res 與 A-ExpJ，它們都能在一次掃描中得到 m 個樣本。非常適合在流處理的場合中。

A-Res 演算法

A-Res(Algorithm A With a Reservoir) 是 Algorithm 的“蓄水池”版本，即維護含有 m

個元素的結果集，對每個新元素嘗試去替換結果集中權重最小的元素。步驟如下：

  1. 將集合 $V$ 的前 $m$ 個元素放入結果集合 $R$。
  2. 對於結果集裡的每個元素，計算特徵值 $k_i = u_i^{(1/w_i)}$，其中 $u_i = rand(0, 1)$
  3. 對 $i = m+1, m+2, \dots, n$ 重複步驟 4 ~ 6
    1. 將結果集中最小的特徵 $k$ 作為當前的閾值 $T$
    2. 對於元素 $v_i$，計算特徵 $k_i = u_i^{(1/w_i)}$，其中 $u_i = rand(0, 1)$
    3. 如果 $k_i > T$ 則將 $R$ 中擁有最小 $k$ 值的元素替換成 $v_i$。

論文證明了如果權重 $w_i$ 是一般連續分佈上的隨機變數，則上面的演算法中插入 $R$ 的次數為 $O(m \log(\frac{n}{m}))$。該演算法用 Python 實現如下：

import heapq
import random

def a_res(samples, m):
    """
    :samples: [(item, weight), ...]
    :k: number of selected items
    :returns: [(item, weight), ...]
    """

    heap = [] # [(new_weight, item), ...]
    for sample in samples:
        wi = sample[1]
        ui = random.uniform(0, 1)
        ki = ui ** (1/wi)

        if len(heap) < m:
            heapq.heappush(heap, (ki, sample))
        elif ki > heap[0][0]:
            heapq.heappush(heap, (ki, sample))

            if len(heap) > m:
                heapq.heappop(heap)

    return [item[1] for item in heap]

A-ExpJ 演算法

A-Res 需要對每個元素產生一個隨機數，而生成高質量的隨機數有可能會有較大的效能開銷，，所以論文中給出了 A-ExpJ 演算法，能將隨機數的生成量從 $O(n)$ 減少到 $O(m\log(\frac{n}{m})))$。從步驟上看，很像我們最開始提出的簡單版本，設定一個閾值並跳過一些元素。具體步驟如下：

  1. 將集合 $V$ 的前 $m$ 個元素放入結果集合 $R$。
  2. 對於結果集裡的每個元素，計算特徵值 $k_i = u_i^{(1/w_i)}$，其中 $u_i = rand(0, 1)$
  3. 將 $R$ 中小最的特徵值記為閾值 $T_w$
  4. 對剩下的元素重複步驟 5 ~ 10
    1. 令 $r = rand(0, 1)$ 且 $X_w = \log( r )/\log(T_w)$
    2. 從當前元素 $v_c$ 開始跳過元素，直到遇到元素 $v_i$，滿足
    3. $w_c + w_{c+1} + \dots + w_{i-1} \lt X_w \le w_c + w_{c+1} + \dots + w_{i-1} + w_{i}$
    4. 使用 $v_i$ 替換 $R$ 中特徵值最小的元素。
    5. 令 $t_w = T_w^{w_i}$, $r2 = rand(t_w, 1)$, $v_i$ 的特徵 $k_i = r_2^{(1/w_i)}$
    6. 令新的閾值 $T_w$ 為此時 $R$ 中的最小特徵值。

Python 實現如下：

def a_expj(samples, m):
    """
    :samples: [(item, weight), ...]
    :k: number of selected items
    :returns: [(item, weight), ...]
    """

    heap = [] # [(new_weight, item), ...]

    Xw = None
    Tw = 0
    w_acc = 0
    for sample in samples:
        if len(heap) < m:
            wi = sample[1]
            ui = random.uniform(0, 1)
            ki = ui ** (1/wi)
            heapq.heappush(heap, (ki, sample))
            continue

        if w_acc == 0:
            Tw = heap[0][0]
            r = random.uniform(0, 1)
            Xw = math.log(r)/math.log(Tw)

        wi = sample[1]
        if w_acc + wi < Xw:
            w_acc += wi
            continue
        else:
            w_acc = 0

        tw = Tw ** wi
        r2 = random.uniform(tw, 1)
        ki = r2 ** (1/wi)
        heapq.heappop(heap)
        heapq.heappush(heap, (ki, sample))

    return [item[1] for item in heap]

驗證

我們用多次取樣的方式來嘗試驗證演算法的正確性。下面程式碼1中為 a、b、c 等元素賦予了不同的權重，取樣 10 萬次後計算被取樣的次數與元素 a 被取樣次數的比值。

overall = [('a', 10), ('b', 20), ('c', 50), ('d', 100), ('e', 200)]
def test_weighted_sampling(func, k):
    stat = {}
    for i in range(100000):
        sampled = func(overall, k)
        for item in sampled:
            if item[0] not in stat:
                stat[item[0]] = 0
            stat[item[0]] += 1
    total = stat['a']
    for a in stat:
        stat[a] = float(stat[a])/float(total)
    print(stat)

首先驗證 A-Res 演算法：

test_weighted_sampling(a_res, 1)
test_weighted_sampling(a_res, 2)
test_weighted_sampling(a_res, 3)
test_weighted_sampling(a_res, 4)
test_weighted_sampling(a_res, 5)

# output
{'e': 19.54951600893522, 'd': 9.864110201042442, 'c': 4.842889054355919, 'a': 1.0, 'b': 1.973566641846612}
{'b': 2.0223285486443383, 'e': 12.17949833260838, 'd': 8.95287806292591, 'c': 4.843410178338408, 'a': 1.0}
{'a': 1.0, 'e': 6.166443722530097, 'd': 5.597171794381808, 'b': 1.9579591056755208, 'c': 4.387922797630423}
{'b': 1.8358898492044953, 'e': 2.5878688779880092, 'c': 2.4081341327311896, 'd': 2.549897479820395, 'a': 1.0}
{'a': 1.0, 'd': 1.0, 'c': 1.0, 'b': 1.0, 'e': 1.0}

看到，在取樣一個元素時，b 被取樣到的次數約為 a 的 2 倍，而 e 則約為 20 倍，與overall 陣列中指定的權重一致。而取樣 5 個元素時，所有元素都會被選中。

同理驗證 A-ExpJ 演算法：

test_weighted_sampling(a_expj, 1)
test_weighted_sampling(a_expj, 2)
test_weighted_sampling(a_expj, 3)
test_weighted_sampling(a_expj, 4)
test_weighted_sampling(a_expj, 5)

# output
{'e': 19.78311444652908, 'c': 4.915572232645403, 'd': 9.840900562851782, 'a': 1.0, 'b': 1.9838649155722325}
{'e': 11.831543244771057, 'c': 4.709157716223856, 'b': 1.9720180893159978, 'd': 8.75183719615602, 'a': 1.0}
{'d': 5.496249062265567, 'c': 4.280007501875469, 'e': 6.046324081020255, 'b': 1.9321080270067517, 'a': 1.0}
{'a': 1.0, 'd': 2.5883654175335105, 'c': 2.440760540383957, 'e': 2.62591841571643, 'b': 1.8787559581808126}
{'a': 1.0, 'd': 1.0, 'c': 1.0, 'b': 1.0, 'e': 1.0}

與 A-Res 的結果類似。

小結

文章中介紹了 A-Res 與 A-ExpJ 兩種演算法，按照步驟用 Python 實現了一個簡單的版本，最後用取樣的方式驗證了演算法的正確性。

加權隨機取樣本身不難，但如果需要在一次掃描中完成就不容易了。難以想像上面的演算法直到 2006 年才提出。演算法本身如此之簡單，也讓不感嘆數學與概率的精妙。

參考

  1. 1.修改自 https://blog.xingwudao.me/2017/09/26/sampling/ ↩

相關推薦

加權隨機取樣 (Weighted Random Sampling)

一個集合裡有 n 個元素，每個元素有不同的權重，現在要不放回地隨機抽取 m 個元素，每個元素被抽中的概率為元素的權重佔總權重的比例。要怎麼做呢？ 簡單的解法 現在考慮只抽取一個元素，假設權重之和為 1。我們可以從 [0, 1] 中隨機得到一個權重，假設為 0.71，而後從第一個元素開始，不斷累加它們的權重，直

隨機取樣和隨機模擬：吉布斯取樣Gibbs Sampling

吉布斯取樣演算法詳解 為什麼要用吉布斯取樣 什麼是sampling? sampling就是以一定的概率分佈，看發生什麼事件。舉一個例子。甲只能E：吃飯、學習、打球，時間T：上午、下午、晚上，天氣W：晴朗、颳風、下雨。現在要一個sample，這個sample可以是：打

隨機模擬的基本思想和常用取樣方法（sampling）

通常，我們會遇到很多問題無法用分析的方法來求得精確解，例如由於式子特別，真的解不出來； 一般遇到這種情況，人們經常會採用一些方法去得到近似解（越逼近精確解越好，當然如果一個近似演算法與精確解的接近程度能夠通過一個式子來衡量或者有上下界，那麼這種近似演算法比較好，因為人們可以

Python-- 有權重的隨機選擇， Weighted Random Choice

import random def windex(lst):     '''an attempt to make a random.choose() function that makes weighted choices     accepts a list

隨機取樣池化--S3Pool: Pooling with Stochastic Spatial Sampling

本文將常規池化看作兩個步驟： 1）以步長為1在特徵圖上滑動池化視窗，尺寸大小基本保持不變， leaves the spatial resolution intact 2）以一種 uniform 和 deterministic 的方式進行降取樣 我們認為

機器學習方法(八)：隨機取樣方法整理（MCMC、Gibbs Sampling等）

轉載請註明出處：Bin的專欄，http://blog.csdn.net/xbinworld 本文是對參考資料中多篇關於sampling的內容進行總結+搬運，方便以後自己翻閱。其實參考資料中的資料寫的比我好，大家可以看一下！好東西多分享！PRML的第11章也是sampling，

【編程珠璣】【第一章】生成隨機數、隨機取樣的問題

當前 rand 可用 生成 奇數 sel 浪費 print 運行時 一、利用隨機數函數生成隨機數 問題1（《編程珠璣》習題12.1後半段）： 給定一個rand()，可以產生從0到RAND_MAX的隨機數，其中RAND_MAX很大（常見值:16位int能表示的最大整數3276

隨機森林（Random Forest）--- 轉載

市場營銷 ssi -o afr actual 所有 很好 struct 驗證 1 什麽是隨機森林？ 　　作為新興起的、高度靈活的一種機器學習算法，隨機森林（Random Forest，簡稱RF）擁有廣泛的應用前景，從市場營銷到醫療保健保險，既可以用來做市場營銷模擬的建模，統

Python-隨機模塊-random

star abcdefg 模塊 abcd IV AC operator sample pre random 　　生成隨機變量 生成 [0, 1) 隨機數　　　　　　　　　　　　.random import random random_number = random.

python生成隨機密碼：random庫

要求： （1）使用random庫，採用0x1010作為隨機種子。 （2）密碼由26個字母的大小寫，10個數字字元和！@#￥%……&*等8個字元組成。 （3）每個密碼長度固定在10個字元。 （4）程式執行每次長生20個密碼，每個密碼一行。 （5）每次長生的20個密碼儲存

關於亂序（shuffle）與隨機取樣（sample）的一點探究

shuffle shuffle的意思就是讓序列亂序，本質上就是讓序列裡面的每一個元素等概率分佈在序列的任何位置。、 使用random.shuffle(iterable), 這樣會直接改變iterable的順序，需要注意的是random.shuffle()函式沒有返回值。 random模組

JavaScript隨機數函式,隨機整數函式,random()

<script> document.write(parseInt(10*Math.random())); 　//輸出0～10之間的隨機整數 document.write(Math.floor(Math.random()*10+1)); 　//輸出1～10之間的隨機整數 fun

3. 集成學習（Ensemble Learning）隨機森林（Random Forest）

總結 子節點 clas 支持向量機 2個 最終 分類算法 容易 oot 1. 前言 相信看了之前關於集成學習的介紹，大家對集成學習有了一定的了解。本文在給大家介紹下遠近聞名的隨機森林（RF）算法。 隨機森林是集成學習中可以和梯度提升樹GBDT分庭抗禮的算法，尤其是它可以很方

3. 整合學習（Ensemble Learning）隨機森林（Random Forest）

1. 前言 相信看了之前關於整合學習的介紹，大家對整合學習有了一定的瞭解。本文在給大家介紹下遠近聞名的隨機森林（RF）演算法。 隨機森林是整合學習中可以和梯度提升樹GBDT分庭抗禮的演算法，尤其是它可以很方便的並行訓練，在如今大資料大樣本的的時代很有誘惑力。 2. 隨機森林原理 隨機森林是Baggin

ArcMap中解決時態資料隨機取樣的動畫生成問題

在做ArcMap操作練習的時候遇到這樣的問題，根據已有的颱風資料（時態資料）來生成其動畫。 資料：2016/09/10  --  2016/09/16 的時態颱風資料 要素已根據時間唯一值進行渲染。 生成時間動畫： ①啟用時間，並設定時間步長間隔為1天。

[轉]隨機取樣一致性（RANSAC）

作者：桂。 時間：2017-04-25  21:05:07 連結：http://www.cnblogs.com/xingshansi/p/6763668.html  前言 仍然是昨天的問題，別人問到最小二乘、霍夫變換、RANSAC在直線擬合上的區別。昨天梳理

隨機函式模組random

random.random(),隨機指定0-1之間的浮點數,含左不含右 import random v = random.random() print(v) 輸出: 0.3250622831488351 random.uniform(1,5)隨機指定1-5之間的浮點數,左右都

python--隨機函式（random uniform randint randrange shuffle sample）

random() random()方法：返回隨機生成的一個實數，它在[0,1)範圍內 運用random()方法的語法： import random #random()方法不能直接訪問，需要匯入random模組，然後通過random靜態物件呼叫該方法 random.random rando

Random Walk 隨機漫步理論 Random Walk Theory 隨即漫步應用

目錄 馬爾可夫鏈 更高的緯度 異常擴散 應用 變種 在圖表上 也可以看看 參考 參考書目 外部連結 Random walk 文章來源：https://en.wikipedia.org/wiki/Random_

隨機搜尋（Random Searching）演算法概述

個人理解，隨機搜尋是一種在巨大資料規模下執行一個耗時上無法接受的程式的優化方法。 它可以用以對一個搜尋演算法施展優化的前提是： 1、資料規模大，精確的結果難以在一定時間計算出。 2、結果的些許的不精確能夠被接受。 3、求取的結果是最優化（optimization）問題，有