珠玉のプログラミングのサンプリング問題のでかいケース

珠玉のプログラミング 本質を見抜いたアルゴリズムとデータ構造

あれ？復刊されてたのか！私が持ってるのはピアソンのやつ。

ちょっと前のCodeIQの問題を（題意に反してちゃんと）と解こうとしていて、この本のサンプリング問題を思い出しました。
この章は母集団からn個の値をソート済みで重複なく取り出す安い方法の話。

N = 10
M = 3

select = M
remaining = N

while select > 0
  if rand < select.fdiv(remaining)
    puts N - remaining
    select -= 1
  end
  remaining -= 1
end

ループのたびに確率を変えていくのが面白いところ！

件の問題は母集団がものすごく大きいケースなんだけど、本の通りやるとものすごく時間がかかるので、もっと安い方法が必要です。
たぶん、つぎの出現までの間隔の方を乱数で求められれば大丈夫なはず。

最近、統計の本の勉強会してたのでなんか上手くできそうな気がしてやってみた。

Head First Statistics ―頭とからだで覚える統計の基本

母集団がN個、アタリがM個、x番目に最初のアタリが出る確率を元のアルゴリズムを使って計算するとx-1回ハズレてx回目にアタリなので、こんなの。

各分数の分子はM個先の分母と同じで、約分できそう。

Mが1のとき

2のとき

3のとき

これが確率密度関数なのかな。積分して逆関数を作れば、N個中M個がアタリのときの最初のアタリの分布の乱数が作れます。たぶん、次の関数みたいになる。

N = 10000
M = 5
BIN = 10

def generate_0(select, remaining)
  ary = []
  (remaining - 1).downto(0) do |x|
    next if rand > select.fdiv(x)
    select -= 1
    ary.unshift(x)
    break if select <= 0
  end
  ary
end

def generate_1(select, size)
  ary = []
  tail = 0
  while select > 0
    tail = (1 - (1 - rand) ** (1.quo(select))) * (size - tail - select) + tail
    ary << tail
    select -= 1
  end
  ary
end

histogram = M.times.collect {[0] * BIN}

100000.times do
  generate_1(M, N).zip(histogram) do |v, h|
    h[v / (N / BIN)] += 1
  end
end

histogram.each do |line|
  puts line.join("\t")
end

generate_0が元のアルゴリアズム、generate_1が今回作ったアルゴリズム。前者がO(N)、後者がO(M)です。

# 元のアルゴリズム
real	1m44.655s
user	1m44.066s
sys	0m0.088s

# 今回のアルゴリズム
real	0m0.844s
user	0m0.825s
sys	0m0.008s

ちなみに分布のヒストグラムはこんな感じです。

Nの大きさによらず、重複なくアタリを選ぶことができました。もともと解きたかった問題はNが5垓くらいでしたから、元のアルゴリズムでは無理でした。でもなー、FloatにBigNumを乗じても5垓種類の値は出ない感じするなー。でも満足したー。