ランダムサンプリング

Takami Torao #RandomSampling #WRS #ReservoirSampling
  • このエントリーをはてなブックマークに追加

概要

ランダムサンプリング (RS; random sampling) または無作為抽出、乱択とは、ある集合から無作為に要素を選択すること。社会科学分野での統計のために利用されることも多いが、このページではソフトウェアアーキテクチャに焦点を当てている。

  1. 概要
  2. 特徴
  3. 単一のサンプリングアルゴリズム
    1. 累積和法
    2. Alias 法
      1. ロビンフッド法
      2. 正方ヒストグラムの生成
      3. サンプリングの実行
      4. Alias 法サンプリングの実装
  4. 複数のサンプリングアルゴリズム
    1. Naive サンプリング
    2. Reservoir サンプリング法
      1. Algorithm A
      2. Algorithm A-Res
      3. Algorithm A-ExpJ
    3. アルゴリズム比較
  5. 参考文献

特徴

集合の各要素が「選ばれやすさ」に関連するパラメータを持っている場合は重み付きランダムサンプリング (WRS; weighted random sampling) あるいは非等価ランダムサンプリング (unequal random sampling) と呼ばれる。WRS に対して重みを持たないランダムサンプリングを一様ランダムサンプリング (uniformed random sampling) と呼ぶ。分散システムにおいてマシンの相対的な性能でタスク分配の優先度を重み付けしたり、Web 広告での関連度や出稿料金でランダムに広告を選択する場合などでは WRS が利用される。

母集団から続けてサンプリングを行うとき、次の要素を選択する前に抽出した要素を元の母集団に戻す方法を復元サンプリング (sampling with replacement) と呼ぶ。復元サンプリングではそれぞれの抽選は事象として独立しており、\(n\) 個の母集団に対して \(m\) のサンプリングを試行したときに各要素が何回選ばれているかを確率変数とする離散確率分布を多項分布で表すことができる。同一の要素が複数回選択される可能性がある (したがって \(m\) 回の試行で選択される要素はかならず \(m\) 個以下となる)。一般的に復元サンプリングは母集団が非常に大きく重複選択の起きる可能性が無視できる場合や、後述するように選択機会の公平性が優先される場合に適用される。

次の要素を選択する前に抽出された要素を母集団には戻さずまだ選択されていない要素のみで抽選を続ける方法を非復元サンプリング (sampling without replacement) と呼ぶ (これらは壺問題における復元と同じ意味)。非復元サンプリングではサンプリングの対象となる母集団は \(m\) 回のサンプリングまでに毎回変化することから、それぞれのサンプリングは事象として独立しておらず、サンプリングごとに選択確率が変化する。詳細は非復元サンプリングの公平性で説明している。

\(n\) 個の母集団に \(k\) 個の「当たり」が含まれているとき、\(m\) 回のサンプリングで引き当てる当たりの個数を確率変数とする離散確率分布を、復元サンプリングでは二項分布、非復元サンプリングでは超幾何分布で表すことができる。これは一定数の故障ノードが含まれる集合から WRS で合意のための委員会クラスタを選出するときに \(2f+1\) や \(3f+1\) の前提が破られる確率分布などに使用することができる。

母集団が巨大で全ての要素を検討するには非効率であったり、母集団の全ての要素を検討する必要がない場合は、サンプリングの前に部分母集団を作成する多段サンプリング (multi-stage sampling)層化サンプリング (stratified sampling) といった体系的サンプリング (systematic sampling) の手法が使われる。それらに対し、母集団に対して単純にサンプリングを行う手法を単純サンプリング (simple sampling) と呼ぶ。

ランダムサンプリングは乱数を使用した非決定論的な標本抽出手法である。もし結果に非決定性が不要であれば (重み付き) ラウンドロビンも検討することができる。

単一のサンプリングアルゴリズム

母集団から無作為に 1 つの要素を選択するアルゴリズムについて。

累積和法

累積和法 (cumulative sum) は Algorithm D [1] とも呼ばれるシンプルな WRS アルゴリズムである。母集団に含まれている要素を直列に並べ、一様乱数が重み累積のどの要素の区間に含まれるかで判定する (これは累積離散確率分布から逆変換サンプリングで乱数を生成する方法と同じである)。個人的に重みの円グラフを高速回転させ矢を放つ Fig 1 で示すような回転ダーツを用いて説明することが多い。

Fig 1. Weighted Random Sampling by cumulative sum

以下は累積和法による重み付き非復元サンプリングで選択された要素のインデックスを参照するサンプルコードである。

use rand;

use crate::wrs::Weighted;

pub fn sampling<T: Weighted>(population: &[&T]) -> usize {
  let total_weight = population.iter().map(|e| e.weight()).sum::<f64>();
  let threshold = total_weight * rand::random::<f64>();
  let mut cumsum = 0.0;
  for (i, e) in population.iter().enumerate() {
    if cumsum <= threshold && threshold < cumsum + e.weight() {
      return i;
    }
    cumsum += e.weight();
  }
  panic!()
}

単純な累積和法は線形探索となることから \(O(n)\) とあまり効率的ではないが、母集団に変化がないのであれば平衡二分探索木 (AVL 木) や B-Tree のような \(O(\log n)\) の構造を持つ探索木を実装してインデックスとして使用することで、繰り返しのサンプリングを高速化することができる。

Alias 法

Alias 法 (Walker's alias method) は離散確率分布からランダムサンプリングを行うための効率的なアルゴリズム。正方ヒストグラム (square histogram, エイリアス) [2] を作成し、横方向と縦方向のたかだか 2 回の一様ランダムサンプリングで WRS を行うことができる。正方ヒストグラムの作成に \(O(n)\) を要するが、サンプリングは \(O(1)\) と非常に効率的であるため、要素やそれらの重み付けに変化がない母集団に対してサンプリングを繰り返して行うケースに向いている。

ロビンフッド法

正方ヒストグラムはロビンフッド法 (Robin Hood method) で作成することができる。これは平均より高いものから低いものへ、低いほうが平均と一致する量の富を移動する操作を、貧富の差がなくなるまで (全てが平均と等しい富を持つまで) 繰り返す方法である (Fig 2 参照)。

各要素 \(i \in \{1,2,\ldots,n\}\) に対応する重みを \(w_i\) としたとき、最終的に全ての \(w_i\) はそれらの平均 (\(\ref{sqhg_average}\)) と等しくなる。\[ \begin{equation} a = \frac{\sum_i w_i}{n} \label{sqhg_average} \end{equation} \] 富の移動 1 回の操作で必ず富の低いほうが平均と等しくなり、また最後の 1 回は必ず両者が平均と等しくなる。したがってロビンフッド法は最大でも \(n-1\) の操作で完了する。富の移動によって富者が貧者となる可能性がある点に注意 (例えば Fig 2 の 2)。

Fig 2. Robin Hood 法による正方ヒストグラムの作成手順

正方ヒストグラムの生成

Alias 法ではロビンフッド法に基づいて以下の手順で \(K\) と \(U\) を生成する。

  1. \(K_i=i\) で初期化する。
  2. \(w_j \lt a\) となるような \(j\) が存在するかぎり以下の 3 ステップを繰り返す。
    1. \(w_i \gt a\) となるような \(i\) を適当に選択する (\(j\) が存在するなら \(i\) も必ず存在する)。以下の 2 ステップは \(i\) から \(j\) へ \(a - w_j\) の量を移動する意味の操作である。
    2. \(K_j=i\) を設定する。
    3. \(w_i \to w_i - (a - w_j) = w_i + w_j - a\) へ再設定する。
  3. すべての区分 \(d_i\) に対してしきい値 \(U_i = \frac{w_i}{a}\) を計算する。

この操作によって 1 つまたは 2 つの領域を持つ等しい高さ \(a\) の区分 \(d_i\) で構成されるヒストグラムができあがる (Fig 3)。区分 \(d_i\) の下の領域は要素 \(i\) の確率域、上の領域は要素 \(K_i\) の確率域を表している。そして \(U_i\) は区分 \(d_i\) における \(i\) と \(K_i\) の確率域の境界を表している。

Fig 3. 正方ヒストグラム

\(a=1\) となるように (つまり \(\sum_i w_i = n\) となるように) 各重み \(w_i\) を設計すると \(w_i=U_i\) となるため操作 3. は不要になる。また誤差がなく \(U_i = 1\) と正確に一致させることができる場合はその区分の \(K_i\) が参照されることもないため操作 1. も省略可能である。

操作 2.1. の判断があるため正方ヒストグラムの作成は浮動小数点の誤差に非常に敏感である。\(a\) の算出に除算が用いられることから、多くのケースでは浮動小数点の加算と比較が多用されることに注意しなければならない。特に、極端に大きな重みを持つ要素が存在するケースでは操作 2.3. の繰り返しで桁落ち、情報落ちの誤差の累積が大きくなる。これは、\(w_j \lt a\) となるような \(j\) が存在しなくなったにも関わらず、\(w_i \gt a\) となるような \(i\) が存在する (ように見える) といった現象で顕在化する。

ラフな対処としては、\(w_j \lt a\) となるような \(j\) が存在しなくなった時点で \(w_i \gt a\) に見えている全ての \(i\) を \(w_i = a\) とする方法が考えられる。あるいはより正確に、計算機イプシロン \(\epsilon \simeq 10^{-15}\) (倍精度浮動小数点の場合)、ループカウンタ \(m\in\{1,\ldots\}\) として、操作 2.1. の判断を \(w_i/a-1 \gt m\epsilon\)、\(a/w_j-1 \gt m\epsilon\) のように累積分の機械イプシロンの幅を想定して許容範囲を設けることで対処することもできる。

サンプリングの実行

サンプリング操作は、ランダムに選んだ区分 \(d_i\) に対して一様乱数が上部と下部のどちらの確率域を指しているかでサンプルを決定する。以下の手順では 1 回の一様乱数生成でその 2 つを行っている。

  1. \(0 \le x \lt 1\) となるような一様乱数 \(x\) を生成する。
  2. \(i = \lfloor n x \rfloor + 1\), \(y = n x + 1 - i\) とする (つまり \(i \in \{1,2,\ldots,n\}、0 \le y \lt 1\))。
  3. 区分 \(d_i\) において \(y \lt U_i\) であれば下の確率域の \(i\)、そうでなければ上の確率域の \(K_i\) をサンプルとする。

Alias 法サンプリングの実装

以下は Alias 法に基づく WRS の実装例である。ここでは \(w_i \gt a\) となるものを overfull, \(w_i \lt a\) となるものを underfull となるように分類している。

複数のサンプリングアルゴリズム

trials = 0 ; rand() × 0 , exp() × 0 , log() × 0

Naive サンプリング

ある母集団から複数の標本を選択する最も単純な方法は、1 つの標本を選択するランダムサンプリング (抽選) を既定の個数が選択されるまで繰り返すことである。

この方法での復元サンプリングは各要素ごとに何回選択されたかの回数が生成される。これは選択回数を確率変数とした離散確率分布 (多項分布) と等価である。

一方で、少数の母集団に対する復元サンプリングは重複選択の発生する頻度が高く、それによって規定数からの差異が無視できないなくなることがある。シミュレーションなどで WRS を扱う場合は重複選択なく正確に \(m\) 個が選択されることが重要であるケースが多いことから非復元サンプリングが利用される。非復元サンプリングは \(m\) 回の抽選に対して必ず \(m\) 個の標本を生成するが、それぞれの抽選は独立しておらず、要素数の減少により抽選が進むにつれて当選の期待値が高くなる。

Reservoir サンプリング法

Reservoir サンプリングは任意のデータストリームから非復元サンプリングで特定個数のサンプルを抽出するときに利用できるアルゴリズムである。

Algorithm A

Algorithm A と呼ばれている非復元 WRS アルゴリズム [1] は母集団から式 (\(\ref{reservoir_key}\)) で示すキーが最も大きい \(m\) 個の要素を選択する。\[ \begin{equation} k_i = u_i^{\frac{1}{w_i}} \label{reservoir_key} \end{equation} \] ここで \(u_i\) は要素ごとに割り振られた 0 以上 1 未満の一様乱数を表している。具体的な実装は以下のようになる。

val samples = population
  .sortBy { e => - math.pow(math.random, 1 / e.weight()) }
  .take(m)

A アルゴリズムは非常に単純で見通しが良く簡易実装に向いているが、母集団のサイズと同じ回数 \(O(n)\) の乱数生成とべき乗演算を行わなければならないため、結果的に累積和法より遅くなることも多い。

しかし \(k_i\) の計算は要素ごとに独立しており、重みの総和も必要ないことから、ランダムサンプリングの分散実行や、母集団全体をメモリ上に置くことが現実的ではない巨大なデータ集合に対して 1 パスで WRS を行うことができるという利点を持つ (このときに使用する一時保存領域を Reservoir (貯水池の意味) と呼ぶ)。近年では Apache Spark の RDD のような抽象化ストリームに基づいたデータ処理プラットフォームと親和性の高い WRS アルゴリズムとして注目を集めている。

Reservoir サンプリングでは「Reservoir 内の最小の \(k_{\rm min}\) より大きな \(k\) を持つ要素を見つけたら、\(k_{\rm min}\) を持つ要素と置き換える」という操作が頻繁に行われるため、しばしば最小の要素が常に先頭に配置されるヒープ、優先キュー、ソートキューといったデータ構造を使用する。

Algorithm A-Res

Algorithm A-Res (A-Reservoir) はサイズ \(m\) の Reservoir 領域を使用し、母集団の中で \(k_i\) の大きい \(m\) 個のみを保持する Algorithm A の改良バージョンである。以下のアルゴリズムで処理を行う。

  1. 先頭の \(m\) 個の要素に対して \(k\) を計算し Reservoir に挿入する。
  2. \(m+1\) 個目以降の要素に対して \(k\) を計算し、Reservoir に含まれる要素で最も小さい \(k\) より大きければ、その最も小さい \(k\) を持つ要素と置き換える。
  3. 最終的に Reservoir に残った要素が非復元 WRS によって得られた標本である。
use rand;
use std::collections::BinaryHeap;

use crate::wrs::a::Sample;
use crate::wrs::Weighted;

pub fn sampling<'s, T: Weighted>(population: &[&'s T], m: usize) -> Vec<&'s T> {
  assert!(m <= population.len());
  let mut reservoir = BinaryHeap::new();  // in general, priority queue
  for i in 0..m {
    let k = rand::random::<f64>().powf(1.0 / population[i].weight());
    reservoir.push(Sample { value: population[i], key: k });
  }
  for i in m..population.len() {
    let k = rand::random::<f64>().powf(1.0 / population[i].weight());
    if k > reservoir.peek().unwrap().key {
      reservoir.pop();  // remove a sample with the smallest key
      reservoir.push(Sample { value: population[i], key: k });
    }
  }
  assert!(reservoir.len() == m);
  reservoir.iter().map(|s| s.value).collect()
}

A-Res は A と同様に \(k\) を計算するために比較的コストの高い乱数生成と指数演算が母集団のサイズ \(n\) と同じ回数必要である。

Algorithm A-ExpJ

Algorithm A-ExpJ (A-Exponential Jump) は A-Res の計算量を大幅に削減した改良バージョンである。A-Res において新しい要素が Reservoir に入るまでスキップされる要素の重みの合計を \(w_s\) とすると、\(w_s\) は指数分布に従う乱数とみなすことができる。したがって、要素ごとに乱数や \(k\) を計算する代わりに、指数分布に従う乱数 \(w_s\) を生成し、該当する要素までジャンプすることで代用可能という理論に基づいている。

use rand;
use std::collections::BinaryHeap;

use crate::wrs::a::Sample;
use crate::wrs::Weighted;

pub fn sampling<'s, T: Weighted>(population: &[&'s T], m: usize) -> Vec<&'s T> {
  assert!(m <= population.len());
  let mut reservoir = BinaryHeap::new();  // in general, priority queue
  for i in 0..m {
    let k = rand::random::<f64>().powf(1.0 / population[i].weight());
    reservoir.push(Sample { value: population[i], key: k });
  }
  let mut x = rand::random::<f64>().ln() / reservoir.peek().unwrap().key.ln();
  for i in 0..m {
    x -= population[i].weight();
    if x <= 0.0 {
      let t = reservoir.peek().unwrap().key.powf(population[i].weight());
      let r = (t + (1.0 - t) * rand::random::<f64>()).powf(1.0 / population[i].weight());
      reservoir.pop();
      reservoir.push(Sample { value: population[i], key: r });
      x += rand::random::<f64>().ln() / reservoir.peek().unwrap().key.ln();
    }
  }
  assert!(reservoir.len() == m);
  reservoir.iter().map(|s| s.value).collect()
}

A-ExpJ も A-Res と同様に母集団全体をメモリ上に保持する必要がなくデータストリームとして扱うことができる。

アルゴリズム比較

以下はそれぞれの非復元 WRS アルゴリズムのコスト理論値と、上記のサンプル Rust コードを実際に Core i7 2.7GHz macOS 10.14 で実行した結果である。

累積和法 A-Res A-ExpJ
最悪計算量 \(O(m n)\) \(O(n)\) \(O(n)\)
乱数生成 \(O(m)\) \(O(n)\) \(O(m \log \frac{n}{m})\)
べき乗演算 \(0\) \(O(n)\) \(O(m \log \frac{n}{m})\)
\(m=5,n=100\) 524 [nsec] 4,277 [nsec] 723 [nsec]
\(m=25,n=100\) 1,769 [nsec] 4,991 [nsec] 2,331 [nsec]
\(m=25,n=1000\) 9,090 [nsec] 40,357 [nsec] 2,138 [nsec]
\(m=25,n=10000\) 96,336 [nsec] 374,878 [nsec] 2,146 [nsec]
\(m=250,n=1000\) 76,089 [nsec] 43,913 [nsec] 16,592 [nsec]

累積和法は標本数 \(m\) や母集団 \(n\) が小さいケースでは速いがそれらの増加に伴って実行時間も大きく増加する傾向にある。A-ExpJ は \(m\) や \(n\) の増加に対する影響が小さく、結果からはスケーラビリティが高い傾向が見られる。全体的な傾向として、A-Res は \(n\) の影響が大きく、A-ExpJ は \(m\) の影響が大きく、累積和法はその両方の影響が大きく出ている。

参考文献

  1. Efraimidis P, Spirakis P. Weighted Random Sampling. In: Kao MY. (eds) Encyclopedia of Algorithms. Springer, New York, NY (日本語訳)
  2. G. Marsaglia, W. W. Tsang, and J. Wang. Fast Generation of Discrete Random Variables, Journal of Statistical Software, 11(3):1-11, 2004. (日本語訳)
  3. 汪 金芳, 標本調査法入門, 平成 17 年5 月 16 日