· 開発

長さのわからないストリームから k 個を公平に選ぶ: リザーバーサンプリング(Algorithm R)

全体の件数がわからないデータを一度だけ読みながら k 個を無作為に選ぶリザーバーサンプリング。置き換え確率が k/i になる理由、どの要素も k/n で残る証明、k 個の組み合わせの均等性、よくある範囲の誤りをまとめました。

開発アルゴリズム

数GBあるログから無作為に100行だけ取り出したいとします。行数はわかりません。先に最後まで数えてからもう一度読んで選ぶとファイルを二度読むことになり、標準入力やネットワークのストリームのように読み直せない入力ではそれもできません。すべてをメモリに載せてから混ぜる方法は、入力が大きくなるほど負担が増えます。

リザーバーサンプリング(reservoir sampling)はこういう場面のための方法です。入力は一度だけ読み、メモリは k 個分の配列ひとつしか使いません。入力がどこで終わっても、その時点で配列に入っている k 個が無作為な標本になります。

Algorithm R

NIST の Dictionary of Algorithms and Data Structures(DADS)は、このアルゴリズムを次のように定義しています。

  1. 最初の k 個はそのまま大きさ k の配列(リザーバー)に入れる。
  2. j 番目の要素(j > k)が来たら、1 以上 j 以下の整数 M を無作為に選ぶ。
  3. M ≤ k なら配列の M 番目を j 番目の要素で置き換える。そうでなければ捨てる。

j 番目の要素がリザーバーに入る確率は k/j で、入る場合にどの枠を追い出すかは k 個の枠から均等に決まります。ひとつの M で両方を同時に決めているわけです。DADS は、この方法が Knuth の Algorithm 3.4.2R をメモリ上で処理するようにした変形であり、Knuth はこれを Alan G. Waterman によるものとしている、と書いています。この記事でも略して Algorithm R と呼びます。

JavaScript で書くと次のようになります。添字は0から始まるので、M の代わりに0以上 seen 未満の m を引き、m < k なら置き換えます。

// 0 以上 n 未満の整数を均等に引く (1 ≤ n ≤ 2^32)
function randomBelow(n) {
  const limit = Math.floor(2 ** 32 / n) * n;
  const buf = new Uint32Array(1);
  do {
    crypto.getRandomValues(buf);
  } while (buf[0] >= limit);
  return buf[0] % n;
}

function reservoirSample(iterable, k) {
  if (!Number.isInteger(k) || k < 0) {
    throw new RangeError("k は 0 以上の整数である必要があります");
  }
  const reservoir = [];
  if (k === 0) return reservoir;

  let seen = 0;
  for (const item of iterable) {
    seen++;
    if (seen <= k) {
      reservoir.push(item);
      continue;
    }
    if (seen > 2 ** 32) throw new RangeError("randomBelow の範囲を超えました");
    const m = randomBelow(seen); // 0 ≤ m < seen
    if (m < k) reservoir[m] = item;
  }
  return reservoir;
}

randomBelow は32ビットの乱数のうち floor(2^32 / n) × n 以上の値を捨てて引き直します。残る値の個数が n の倍数になるので、% n の結果は0から n−1 までそれぞれ同じ個数ずつ対応します。この関数は n が 2^32 を超えると使えないため、要素がそれより多くなりうる入力では、もっと広い範囲の整数乱数が必要です。

要素ひとつあたりの処理は定数なので時間は O(n)、追加のメモリは O(k) です。乱数は k 番目より後の要素ごとにひとつずつ引くので、k ≥ 1 なら合計 max(n − k, 0) 回です。入力が k 個以下なら一度も引かず、k = 0 なら入力を読む前に返るので、やはり0回です。

どの要素も k/n の確率で残る理由

入力が n 個で終わり、1 ≤ k ≤ n だとします。以下はこの記事による計算です。(k = 0 なら関数は空の配列を返すので、n > 0 のときどの要素も残る確率は0です。n = 0 なら確率を考える要素がありません。)

まず、j 番目の要素(j > k)を処理する段階で、リザーバーの特定の枠が追い出される確率を考えます。M がちょうどその枠の番号でなければならないので 1/j です。その枠にいた要素がこの段階を生き残る確率は 1 − 1/j = (j − 1)/j です。

最初の k 個のうちのひとつは確率1で入り、そのあと k+1 番目から n 番目までのすべての段階を生き残る必要があります。

(k/(k+1)) × ((k+1)/(k+2)) × … × ((n−1)/n) = k/n

分子と分母が順に約分され、k/n だけが残ります。

**i 番目の要素(i > k)**は、まず確率 k/i で入り、i+1 番目から n 番目までの段階を生き残る必要があります。生き残る確率は同じ約分で i/n になります。

(k/i) × (i/n) = k/n

どの要素も、最終的なリザーバーに残る確率は k/n で同じです。入力がどこで終わってもこれが成り立つので、n を前もって知る必要はありません。

要素ごとの確率が同じでも、それだけでは足りない

各要素が k/n で選ばれるというだけでは、「公平に選んだ」と言うには足りません。たとえば4個の要素を {1, 2} と {3, 4} の二組に分け、コインを投げてどちらかの組を選ぶと、各要素が選ばれる確率はどれも 1/2 です。それでも {1, 3} のような組み合わせは決して出てきません。

DADS は、Algorithm R ではありうる k 個の組み合わせがどれも同じ確率で出ると書いています。これも帰納法で確かめられます。ここもこの記事による計算で、以下には k で割る段階があるため 1 ≤ k ≤ n を仮定します。k = 0 ならありうる組み合わせは空集合ひとつだけで、上のコードも常に空の配列を返すので、別に証明することはありません。t 個(t ≥ k)を処理した後、最初の t 個から作れる k 個の組み合わせ C(t, k) 通りがすべて確率 1/C(t, k) だと仮定し、t+1 番目の要素を処理します。以下では C(t+1, k) = C(t, k) × (t+1)/(t+1−k) を使います。

  • 新しい要素を含まない組み合わせ S は、直前のリザーバーが S で、新しい要素が捨てられた場合にだけ出ます。確率は (1/C(t, k)) × (t+1−k)/(t+1) = 1/C(t+1, k) です。
  • 新しい要素を含む組み合わせ S について、新しい要素を除いた k−1 個を T とします。直前のリザーバーは T にもうひとつの要素 x を加えたものでなければならず、x の候補は t − (k−1) 個あります。新しい要素が入り(k/(t+1))、x の枠が選ばれる(1/k)必要があるので、全体の確率は (t−k+1) × (1/C(t, k)) × (1/(t+1)) = 1/C(t+1, k) です。

最初の k 個を入れた直後はありうる組み合わせがひとつしかないので仮定は成り立ち、一段階進めても保たれます。したがって最後まで読んだとき、C(n, k) 通りの組み合わせはすべて同じ確率です。

よくある誤りと境界条件

範囲をひとつ間違える。 0から数えるループでは、i は (i+1) 番目の要素です。ここで randomBelow(i + 1) ではなく randomBelow(i) を使うと、入る確率が k/(i+1) ではなく k/i になります。1 ≤ k < n とすると、i = k、つまり k+1 番目の要素では randomBelow(k) が必ず k 未満を返すので、この要素は必ず入ります。同じ前提のもと、前の節と同じ方法で最後まで計算すると(この記事による計算)、最初の k 個が残る確率は (k−1)/(n−1)、それ以外の要素は k/(n−1) になります。最初の k 個だけが不利になる、偏った標本です。k = 0 の場合、この誤ったコードは最初の要素で randomBelow(0) を呼びます。コメントに書いた範囲(1 ≤ n)の外なので、確率を論じる以前に不正な呼び出しです。n ≤ k なら誤った呼び出しまで到達せず、上の確率の式も当てはまりません。上のコードでは seen を先に増やしてから randomBelow(seen) を呼び、1から数えた個数をそのまま使うことでこの混乱を避けています。

入力が k 個より少ない場合。 上の関数は受け取った要素をすべて入れ、k 個より短い配列を返します。長さのわからないストリームではこの場合を前もって防げないので、エラーとするか短い結果を受け入れるかは呼び出し側で決める必要があります。ちなみに Python の random.sample(population, k) は、k が母集団より大きいと ValueError を送出します。

k が0以下の場合。 上のコードは k = 0 なら入力を読まずに空の配列を返し、負の数や整数でない値は RangeError で拒否します。別の方針でもかまいませんが、どちらにするかは明示しておくほうが安全です。負の k をそのままにすると seen <= k が一度も真にならず、黙って空の結果が返りがちです。

結果の順序。 リザーバーに入る要素の集合は均等ですが、配列の中の並び順まで無作為になるわけではありません。入力が k 個以下なら入力順のままで、最初の k 個のうち最後まで残った要素は元の位置にとどまります。順序に意味がある用途なら、結果を Fisher–Yates でもう一度シャッフルしてください。

乱数の品質。 上の証明は、M が1から j まで正確に均等に出ることを前提にしています。浮動小数点の乱数に j を掛けて切り捨てる方法は、j が大きくなるにつれて精度の限界を確かめる必要があります。Python のドキュメントは、標準の生成器であるメルセンヌ・ツイスタは決定的なので暗号用途にはまったく適さず、セキュリティ用途には secrets モジュールを使うよう案内しています。標本を予測されると困る場面なら、アルゴリズムより先に乱数源を見直すべきです。

長さがわかるとき、もっと速くしたいとき

すべての要素がすでにメモリにあり件数もわかっているなら、リザーバーサンプリングは要りません。Fisher–Yates シャッフルを k 回だけ回すか、Python なら random.sample を使えば済みます。ただし random.sample は母集団としてシーケンスを受け取り、Python 3.11 からは集合を自動でリストに変換しなくなりました。ストリームに使うには結局すべてをリストにする必要があるので、一度読んで終わりにしたい場面には向きません。

逆に n が k よりずっと大きいと、ほとんどの要素は乱数をひとつ引いたうえで捨てられます。DADS は Jeffrey Scott Vitter の1985年の論文 “Random Sampling with a Reservoir”(ACM TOMS 11(1))を紹介し、その Algorithm X、Y、Z は要素ごとに読み飛ばすかどうかを決める代わりに、いくつ読み飛ばすかをまとめて決めることで、使う乱数をずっと少なくしていると説明しています。乱数の生成コストが問題になるなら、調べてみる価値のある方向です。

参考資料

いずれも2026-10-13に確認しました。

ブログへ