n個あるうちからk個取り出したい。R だと sample 関数を使うだけなのだが、その仕組みが気になる。
重複を許す場合(with-replacement)は簡単で、1からnまでの一様乱数をk個作るだけ。計算量はO(k)。
重複を許さない場合はどうする? ソースコードを見る前に考えてみる。
まず、１つ目の要素を k/n の確率で選択する。これを選択した場合は、２つ目の要素を (k - 1)/(n - 1) の確率で選択する。１つ目の要素が棄却された場合は、2つ目の要素を k/(n - 1) の確率で選択する。このとき、2つ目の要素も k/n の確率で選択されることは明らか。3つ目以降の要素についても同様に進めればよい。すると、すべての要素について判定しないといけないので、O(n)の時間がかかる。これは嬉しくない。
乱数を１つ作るたびに、要素が１つ抽出されるようにしたい。問題は、重複抽出を認めないので、抽出が進むにつれてどんどん要素が脱落していって、乱数 i を作っても、「アタマから数えて i 個目の残っている要素」を容易には探し出せないこと。アタマから見ていくのは論外。Fenwick木とかを使えばlogオーダになるのは知っているが、そんなにややこしい問題なのか、これ?
というあたりで、答えを見てみる。R のソースコードで sample を grep すると、http://svn.r-project.org/R/trunk/src/main/random.c が見つかる。少し追いかけて、SampleNoReplace がそれだと分かる。
なんと! 要素を抽出するたびに、その脱落した場所に最後の要素を持ってきているのだ。そうすると、前半は常に「残っている要素」で充填されていることになるので、一様乱数でサンプリングできる。
言われてみれば当たり前だが、すこし感激したのでメモ。
注意: R の sample の実装だと、最初に添字を別の配列にコピーしているので、結局 O(n) かかっている。乱数発生は k 回で済むんだけど……
