Reservoir sampling - 巨人の肩に乗ってニヤニヤするブログ

$N$ 件のデータから $k$ 件をランダムサンプリングする場合のアルゴリズムで、
$N$ が未知数であっても、 $N$ を求める必要がなく、
1passでランダムサンプリングができる。
今まで $N$ を求めるためにかかっていた時間を家事や育児にあてることができる。

概要

$N$ 件のデータを $S[0],S[1], \cdots ,S[N-1]$ とする。
このアルゴリズムでは、Reservoir（貯水池）を $R[0],R[1], \cdots ,R[k-1]$ として、
そこに $k$ 件のデータを随時貯めておく。
Reservoirは処理の過程でどんどん更新されていき、
最終的にはこれがサンプリング結果となる。

アルゴリズム

アルゴリズムは基本的に２つのステップから構成される。

ステップ①
最初の $k$ 件をとりあえずReservoirにぶっこんでおく。
$\displaystyle R[i ] \left( 0\leq i < k \right)$ にデータ $S[i]$ を代入。

ステップ②
データ $S[i] \left( k \leq i < N \right)$ については $\frac{k}{i+1}$ の確率でReservoirの保持内容を更新する。
$[0,i]$ の整数値をとる一様乱数 $j$ を生成し、 $j < k$ の場合、 $R[j]$ に $S[i]$ を代入。

解釈

本当にランダムサンプリングされているのか気になるところ。
２つの場合に分けてデータ $S [ i ]$ が選ばれる確率を考えてみる。

データ $S [ i ] \left( i < k \right)$ がReservoirに残る確率
$\hspace{3em} = R [ i ]$ が更新されない確率
$\hspace{3em} = \frac{k}{k+1} \frac{k+1}{k+2} \cdots \frac{N-1}{N}$
$\hspace{3em} = \frac{k}{N}$

データ $S [ i ] \left( k \leq i \right)$ がReservoirに残る確率
$\hspace{3em} = R [ j ]$ に $S [ i ]$ が代入される確率 $\times R [ j ]$ が更新されない確率
$\hspace{3em} = \frac{k}{i} \frac{i}{i+1} \frac{i+1}{i+2} \cdots \frac{N-1}{N}$
$\hspace{3em} = \frac{k}{N}$

実装例

行数が未知のテキストファイルから $k$ 行をReservoir samplingにて取り出す例。
サンプリング結果はテキストファイルとして出力。

import random

def random_sampling(inputFile, k=10000, header=True):
    R = []
    columns = ''
    with open(inputFile) as fin:
        if header: columns = fin.next()
        for i in range(k):
            R.append(fin.next())
        i = k
        for line in fin:
            j = random.randint(0,i)
            if j<k: R[j] = line
            i += 1
    outputFile = inputFile.split('.')[0] + '_random' + str(k) + '.' \
                    + inputFile.split('.')[1]
    with open(outputFile,"w") as fout:
        if header: fout.write(columns)
        for line in R:
            fout.write(line)

参考

Reservoir sampling - Wikipedia, the free encyclopedia