再帰を用いた無作為抽出標本の生成のコード

Question

MITのアルゴリズムイントロダクション第３版の5章の練習問題5.3-7の質問です。
この問題は無作為抽出標本を生成することを証明する問題です。
下の擬似コードが{1,2,3,.....,n}から0＜＝m＜＝nであるm部分集合Sを無作為抽出することの証明です。

具体的な数字を当てはめて動きを追うことはできるのですが、感覚的に無作為抽出になることも理解できず、もちろんループ不変式を用いた帰納的な証明もできません。
再帰のサイクルの途中段階では最終的なループより小さい集合から抽出しているし、特に6行目は理解できません。
後もしこのテキストの解答が閲覧できるサイトがあれば教えていただければ幸いです。

Comments

[swordone]

RANDOM(1, n)というのは、1以上n以下の整数をランダムに得るということで合っていますか？

[macuserdesu]

合ってます！

[episteme]

↑それとも1以上n"未満"?

[macuserdesu]

1以上n以下です。nも含みます！

Answer 1

無作為抽出した場合、このm部分集合に特定の要素が属する確率を数学的に求めるとm/n(※)なので、すべての要素に対して、部分集合に属する確率がm/nであれば証明ができることになります。

(※)全パターンはnCm、特定の要素(xとする)が属するパターンは、xを除いた(n-1)個から、xとともに属する(m-1)個を選ぶ組み合わせに等しいので、(n-1)C(m-1)=(m/n)×nCm。よって確率は{(m/n)×nCm}/nCm = m/n

再帰の繰り返しにより、(0, n-m)に到達し、これは空集合を返します。そして(1, n-m+1)から要素の選定が始まります。これから、n-m+1以下の要素aについては、RANDOMでaが1回でも選ばれれば部分集合に入ります。この確率を直に求めるのは難しいので、「aが1回も選ばれない」確率を求めます。
(1, n-m+1)の時、aが選ばれない確率は(n-m)/(n-m+1)
(2, n-m+2)の時、aが選ばれない確率は(n-m+1)/(n-m+2)
(3, n-m+3)の時、aが選ばれない確率は(n-m+2)/(n-m+3)
…
(m-1, n-1)の時、aが選ばれない確率は(n-2)/(n-1)
(m, n)の時、aが選ばれない確率は(n-1)/n
ここで再帰が終了します。aが1回も選ばれない確率は、これらの確率をすべて掛けて(n-m)/nとなります。
したがって、aが1回以上選ばれる確率は、1-(n-m)/n = m/nとなります。

n-m+2以上の要素bについては、b=n-m+1+iとしたとき、(i, b-1)までに選ばれることはありませんが、(i+1, b)のときだけは、RANDOMでb自身が選ばれるときだけではなく、それまでに選んだi個の要素のいずれかが選ばれたときに、部分集合に入ります。したがって、先ほどと同様にbが1回も選ばれない確率を考えると、
(i+1, b)の時、bが選ばれない確率は(b-i-1)/b
(i+2, b+1)の時、bが選ばれない確率はb/(b+1)
(i+3, b+2)の時、bが選ばれない確率は(b+1)/(b+2)
…
(m-1, n-1)の時、bが選ばれない確率は(n-2)/(n-1)
(m, n)の時、bが選ばれない確率は(n-1)/n
ここで再帰が終了します。bが1回も選ばれない確率は、これらの確率をすべて掛けて(b-i-1)/nとなります。b=n-m+1+iなので、この確率は(n-m)/nです。
したがって、bが部分集合に入る確率は、1-(n-m)/n = m/nとなります。

以上より、どの要素も部分集合に入る確率がm/nであることが示せました。

直感的には、番号が後の方は最初に選ばれるチャンスがない代わりに、あとで大きなチャンスがあるために釣り合っている、といった感じでしょうか。

Comments

[episteme]

そか、無作為抽出になってることを証明するんだから、
どれも同確率で現れることを証明せにゃならんのか。

[swordone]

まあ、明らかにフィッシャー・イェーツのアルゴリズムを応用すれば簡単だと思うんだけど…

[episteme]

n個の玉の入った袋からm個取り出すだけですからねー

Answer 2

1. S = RANDOM_SAMPLE(m-1, n-1) は n を含まない。
2. i = RANDOM(1, n) が既に S に含まれていたら n を追加する。(1.より重複しない)
3. i が S に含まれていないなら i を追加する。(当然重複しない)

...だから S には m個の 1～n が重複なく現れる。

[追記] とんでもなく勘違いしてました。そっと閉じてください orz

で、書いてみた.

#include <iostream>
#include <set>
#include <random>

std::set<int> random_sample(int m, int n) {
  if ( m == 0 ) { 
    return std::set<int>();
  } else {
    std::set<int> S = random_sample(m-1, n-1);
    std::random_device rnd;
    std::uniform_int_distribution<> dist(1,n);
    int i = dist(rnd);
    if ( S.find(i) != S.end() ) {
      S.insert(n);
    } else {
      S.insert(i);
    }
    return S;
  }
}

int main() {
  int hist[10];
  std::fill_n(hist, 10, 0);
  for ( int i = 0; i < 10000; ++i ) {
    for ( int item : random_sample(3, 10)) { ++hist[item-1]; }
  }
  for ( int item : hist ) {
   std::cout << item << ' ';
  }
  std::cout << std::endl;
}

Comments

[swordone]

「無作為抽出」である、つまり各要素が確率的な偏りなく抽出される、という点が疑問なのだと思います。

[episteme]

面目ねぇ... お詫びの実装を追記しました。