100囚人問題をJavaScriptでシミュレート


# 100囚人問題

[100 prisoners problem](https://en.wikipedia.org/wiki/100_prisoners_problem)(100囚人問題)というのは、

> 1から100まで連番の囚人番号が振られた100人の死刑囚(prisoner)に、最後のチャンスが与えられた。
> 
> - いま、100個の引き出し(drawer)のある棚には、1から100までの連番が書いてある紙が、1枚づつランダムに入っている。
> - 囚人は一人に付き50回まで、引き出しを開け中の番号を確認してから閉じることが許されている。
> - 囚人は一人づつ順にこの棚を確認し、100人全員が自分の番号が書かれた紙がある引き出しを当てることができたときのみ、生存することができる。
> - ただし、最初の一人が引き出しを開けた後は、囚人同士でコミュニケーションをとることは一切できない。
> 
> さて囚人たちにとって、最適な手段をとるにはどうすればよいか。

という問題です。

以降の説明は、上記リンク先のWikipedia記事にも書いてある内容への詳細な説明になります。
最後に、JavaScriptコードでのシミュレーションを行い、理想値に近い結果となることを確認するものとなります。


## 100囚人問題の選択戦略

一人が50個の引き出しをランダムに開けたら、100個の中から自分の番号が見つかる確率は1/2です。
このため、100人全員がランダムに開けたら、生存できる確率は1/2の100乗です。これは、30個以上0の続いたほぼゼロの確率です。
ここでは、このやり方を「ランダム選択」と呼ぶことにします。

ただし、以下の手順で引き出しを選んでいけば、生存できる確率がぐんと上がります。

0. 引き出しには連番を割り当てる(左上を1に、その右を2、というように、囚人皆で同じ番号を使う)
1. 各人は、自分の囚人番号と同じ番号の引き出しを最初に開ける
2. 引き出しの中の番号が、自分の囚人番号と同じなら、その引き出しの番号を覚えておく
3. 引き出しの中の番号が、自分の囚人番号ではないなら、その番号の引き出しを開け、2へ

このやり方を「巡回選択」と呼ぶことにします。そして、この巡回選択を取ることで、100囚人では31.18%で生存します。

## 4囚人問題の例

100囚人について考える前に、まず問題を小さくし、4囚人が2回引き出しを開けられる例で考えます。

4つの引き出し内の番号が取りうるパターンは4!=24個です。
引き出し番号順での中の番号がとりうる全パターンは、以下のものになります。

巡回選択で生存するパターン(10個):

- 1234

- 2134
- 3214
- 4231
- 1324
- 1432
- 1243

- 2143
- 3412
- 4321

巡回選択で生存できないパターン(14個):

- 1342
- 1423
- 3241
- 4213
- 2431
- 4132
- 2314
- 3124

- 2341
- 4123
- 2413
- 3142
- 3421
- 4312

生存パターンかどうかの違いについてはあとまわしにしますが、巡回選択での生存確率は、10/24 = 0.41666...の41.7%です。

これがランダム選択だと、1/2の4乗であり、1/16=0.0625の6.25%になります。
4囚人でも、巡回選択を選ぶことで、生存率は格段に上がることになります。

## 巡回選択での生存確率

この巡回選択をするというのは、自分の囚人番号から始め、引き出しの中の番号への**置換**を、50回まで適用することを行うことです。
この巡回選択では引き出し全体を、100要素の[置換群](https://ja.wikipedia.org/wiki/%E5%AF%BE%E7%A7%B0%E7%BE%A4)の一つになっている、とみなすものです。

この置換群は、いくつかの巡回置換の組として、分解することができます。
そして100囚人では、置換群の中の**最大の巡回置換の要素数**が50以下なら生存でき、51以上なら生存できない、ということになります。

4囚人の例では、取りうる引き出しのパターンは、

- 巡回1が4つの組(1): 1234
- 巡回2が1つ、巡回1が2つの組(4C2=6): 2134, 3214, 4231, 1324, 1432, 1243
- 巡回2が2つの組(3): 2143, 3412, 4321
- 巡回3が1つ、巡回1が1つの組(4x2=8): 1342, 1423, 3241, 4213, 2431, 4132, 2314, 3124
- 巡回4が1つの組(3!=6): 2341, 4123, 2413, 3142, 3421, 4312

になっています。このうち、巡回2までのパターンが生存できるパターンになっていたのでした。

---

以降は、囚人数を一般化したn囚人問題について扱います。

総要素数nによる置換群のうち、半分より多い要素数lの巡回がある組み合わせ数は、

- n個の中から巡回置換させるl個を選ぶ[組み合わせ数](https://ja.wikipedia.org/wiki/%E7%B5%84%E5%90%88%E3%81%9B_(%E6%95%B0%E5%AD%A6)): ${}_nC_l = \frac{n!}{l!(n-l)!}$
- 選んだl要素でつくる巡回のパターン数: $(l-1)!$
    - 巡回4の例を用いると、まず1を固定することで、巡回での置換を、1=>2=>3=>4=>1のようにならべると、間の234の[順列の個数](https://ja.wikipedia.org/wiki/%E9%A0%86%E5%88%97)分だけの巡回パターンがある
- 選んだl要素以外の残りの要素の間での置換がとりうる組み合わせ数: $(n-l)!$

をすべてかけた、$\frac{n!}{l}$ 個の組み合わせ数になります。

よって、n囚人問題での巡回選択で生存できる引き出しの組み合わせの個数は、引き出しが取りうる全パターン数$n!$から、半分より多い巡回のあるパターンをすべて差し引いた

```math
n! - \sum_{l=n/2+1}^{n} \frac{n!}{l}
```

であり、生存確率は、これを引き出しが取りうる全パターンの数$n!$で割った

```math
1 - \sum_{l=n/2+1}^{n} \frac{1}{l}
```

になります。

この式から、4囚人では1 - 1/3 - 1/4 = 5/12となるため、生存確率は0.4166...だったのでした。
また100囚人では、生存確率は1 - 1/51 - 1/52 - ... - 1/99 - 1/100であり、0.3118...となります。


## n囚人問題における巡回選択での生存確率の極限値

n囚人での巡回選択での生存確率は、

```math
1 - \sum_{l=\frac{n}{2}+1}^{n} \frac{1}{l}
```

です。

さらに、$\frac{1}{l} = \frac{1}{2l} + \frac{1}{2l} < \frac{1}{2l-1} + \frac{1}{2l}$であることから、たとえば1/3 < 1/5+1/6, 1/4 < 1/7+1/8のような関係になるため、4囚人での生存確率"1-1/3-1/4"より、8囚人での生存確率"1-1/5-1/6-1/7-1/8"のほうが小さくなります。

ただしこの生存確率は、nが増えるほど減っていく値ですが、nを増やしつづけても0まで減るわけではなく、ある[極限値に収束](https://ja.wikipedia.org/wiki/%E6%A5%B5%E9%99%90)する値になります。
以下の話は、この生存確率の極限値を求めるものです。

まず、この生存確率の式を変形し、

```math
\begin{align}
1 - \sum_{l=\frac{n}{2}+1}^{n} \frac{1}{l} &= 1 - \left(\sum_{l=1}^{n} \frac{1}{l} - \sum_{l=1}^{\frac{n}{2}} \frac{1}{l}\right) \\
  &= 1 - \sum_{l=1}^{n} \frac{1}{l} + \sum_{l=1}^{\frac{n}{2}} \frac{1}{l}
\end{align}
```
とします。この式でnの極限をとると、

```math
\lim_{n \to \infty}\left(1 - \sum_{l=1}^{n} \frac{1}{l} + \sum_{l=1}^{\frac{n}{2}} \frac{1}{l}\right)
```

となり、以降この極限の式について取り組みます。

この極限値の導出のために、[オイラー定数γ](https://ja.wikipedia.org/wiki/%E3%82%AA%E3%82%A4%E3%83%A9%E3%83%BC%E3%81%AE%E5%AE%9A%E6%95%B0)

```math
\gamma = \lim_{n \to \infty}\left(\sum_{l=1}^{n}\frac{1}{l} - \log{n}\right)
```

を使用します。(無限大に発散する[調和級数](https://ja.wikipedia.org/wiki/%E8%AA%BF%E5%92%8C%E7%B4%9A%E6%95%B0)と、同じく無限大に発散する[自然対数](https://ja.wikipedia.org/wiki/%E8%87%AA%E7%84%B6%E5%AF%BE%E6%95%B0)は、同じくらいゆっくり発散していて、それらの差がこのγ=0.5772...に近づく)

さらに$n=2m$とおき、このγを用いることで、

```math
\begin{align}
\lim_{n \to \infty}\left(\sum_{l=1}^{n/2}\frac{1}{l} - \log{n}\right) 
&= \lim_{m \to \infty}\left(\sum_{l=1}^{m}\frac{1}{l} - \log{2m}\right) \\
  &= \lim_{m \to \infty}\left(\sum_{l=1}^{m}\frac{1}{l} - (\log{2} + \log{m})\right) \\
  &= \lim_{m \to \infty}\left(\sum_{l=1}^{m}\frac{1}{l} - \log{m}\right) - \log{2} \\
  &= \gamma - \log{2} \\
\end{align}
```

という関係も導出できます。

これらを、先の極限の式へ適用すると、

```math
\begin{align}
\lim_{n \to \infty}\left(1 - \sum_{l=1}^{n} \frac{1}{l} + \sum_{l=1}^{\frac{n}{2}} \frac{1}{l}\right) 
&= \lim_{n \to \infty}\left(1 - \left(\sum_{l=1}^{n} \frac{1}{l} - \log{n}\right) + \left(\sum_{l=1}^{\frac{n}{2}} \frac{1}{l} - \log{n}\right)\right) \\
&= 1 - \gamma + (\gamma -\log{2}) \\
&= 1 - \log{2}
\end{align}
```

となって、n囚人問題の生存確率の**極限値は1-log2=0.3065...**となっていたのでした。

# JavaScriptでの100囚人問題のシミュレーション

以降は、100囚人問題の生存確率を、ランダム選択のときと巡回選択のときとを、シミュレーションによって求めてみます。
このシミュレーションでは、ランダムにシャッフルした引き出しを用意してランダム選択と巡回選択とでそれぞれ生存できるかどうかをチェックすることを、たとえば1万回繰り返して、そのうち生存できた総数を数えることで、それぞれの生存確率を計算します。

なお以下のコードでは、JavaScriptの配列インデックスに合わせて、囚人番号は0始まりにしています。

---

先に、JavaScriptの`Array`のためのユーティリティ関数を実装しておきます。

```js:[1]
const range = n => [...Array(n).keys()];
const shuffle = a => {
  for (let i = a.length - 1; i > 0; i--) { // Fisher-Yates
    const j = Math.random() * i >>> 0;
    [a[i], a[j]] = [a[j], a[i]];
  }
  return a;
};
const sample = (a, k) => {
  const n = a.length, idx = range(n), r = [];
  //console.assert(k < n);
  for (let i = 0; i < k; i++) {
    const j = Math.random() * (n - i) >>> 0;
    r.push(a[idx[j]]);
    idx[j] = idx[n - i - 1]; // overwrite chosen j with truncated last element
  }
  return r;
};
```

- `range`は、0からn-1まで順に入った配列を作ります
- `shuffle`は、[Fisher-Yates](https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A3%E3%83%83%E3%82%B7%E3%83%A3%E3%83%BC%E2%80%93%E3%82%A4%E3%82%A7%E3%83%BC%E3%83%84%E3%81%AE%E3%82%B7%E3%83%A3%E3%83%83%E3%83%95%E3%83%AB)でのインプレースなランダムシャッフルです
- `sample`は、aからランダムにk回ピックアップしていった結果を返す関数ですが、添字の配列でピックアップを行うことで、引数の配列は変化させないようにします

---

まず、ある引き出し(`drawers`)に対して、ランダム選択をとる関数`randomPick`の実装です。

```js:[2]
const randomPick = (prisoner, k, drawers) => {
  return sample(drawers, k).some(num => num === prisoner);
};
```

囚人番号`prisoner`と、引き出しを開けられる最大の回数`k`、番号入り引き出しの配列`drawers`を受け取り、自分の囚人番号が見つかったかどうかを返す関数です。

---

そして、`randomPick`と同じ引数リストをとる、巡回選択をとる関数`cyclicPick`の実装です。

```js:[3]
const cyclicPick = (prisoner, k, drawers) => {
  for (let num = prisoner, i = 0; i < k; i++) {
    if ((num = drawers[num]) === prisoner) return true;
  }
  return false;
};
```

`num`を`k`回置換させるまでに、囚人番号`prisoner`に到達できなければ、囚人番号を発見できなかったことになります。

---

以下は、囚人数`prisoners`に合わせたランダムな引き出しを用意し、`randomPick`と`cyclicPick`でそれぞれ`times`回シミュレーションして、成功数をカウントする関数`simulate`です。

```js:[4]
const challenge = (pick, drawers) => {
  const prisoners = drawers.length, k = Math.floor(prisoners / 2);
  return range(prisoners).every(prisoner => pick(prisoner, k, drawers));
};
const simulate = (prisoners, times) => {
  const drawers = range(prisoners);
  let successRandom = 0, successCyclic = 0;
  for (let i = 0; i < times; i++) {
    shuffle(drawers);
    successRandom += challenge(randomPick, drawers);
    successCyclic += challenge(cyclicPick, drawers);
  }
  return {successRandom, successCyclic};
};
```

最後は、シミュレーション設定(`prisoners = 100, times = 10000`)を与えてシミュレーションさせ、結果をアウトプットするメインコードです。

```js:[5]
{
  const prisoners = 100, times = 10000;
  const {successRandom, successCyclic} = simulate(prisoners, times);
  console.log(`${prisoners} Random Pick: ${successRandom}/${times}`); // => 0%
  console.log(`${prisoners} Cyclic Pick: ${successCyclic}/${times}`); // => 31%
}
```

このコードの実行結果は、

```shell-session
100 Random Pick: 0/10000
100 Cyclic Pick: 3138/10000
```

となりました。

巡回選択での生存確率は31.18%で、シミュレーションでも近い値が得られました。

## コード全体

<details><summary>https://gist.github.com/bellbind/5fff0abc164e5f81a56e92b2bfc52d60</summary><div>

```js:100prisoners.js
// Simulation of 100 prisoers problem strategies
// - https://en.wikipedia.org/wiki/100_prisoners_problem

// array utils
const range = n => [...Array(n).keys()];
const shuffle = a => {
  for (let i = a.length - 1; i > 0; i--) { // Fisher-Yates
    const j = Math.random() * i >>> 0;
    [a[i], a[j]] = [a[j], a[i]];
  }
  return a;
};
const sample = (a, k) => {
  const n = a.length, idx = range(n), r = [];
  //console.assert(k < n);
  for (let i = 0; i < k; i++) {
    const j = Math.random() * (n - i) >>> 0;
    r.push(a[idx[j]]);
    idx[j] = idx[n - i - 1]; // overwrite chosen j with truncated last element
  }
  return r;
};

// strategies
const randomPick = (prisoner, k, drawers) => {
  return sample(drawers, k).some(num => num === prisoner);
};
const cyclicPick = (prisoner, k, drawers) => {
  for (let num = prisoner, i = 0; i < k; i++) {
    if ((num = drawers[num]) === prisoner) return true;
  }
  return false;
};

// max size of cyclic permutation in drawers 
const cyclicMax = (drawers) => {
  const checked = new Set();
  let max = 0;
  for (let i = 0; i < drawers.length; i++) {
    if (checked.has(i)) continue;
    checked.add(i);
    let num = drawers[i], count = 1;
    while (!checked.has(num)) {
      checked.add(num);
      num = drawers[num], count += 1;
    }
    max = Math.max(max, count);
  }
  return max;
};

// simulation
const guard = (drawers, k) => new Proxy(drawers, {
  count: 0,
  get(target, key, recv) {
    const n = Number(key);
    if (Number.isInteger(n) && n >= 0) console.assert(this.count++ < k);
    return Reflect.get(target, key, recv);
  },
});
const challenge = (pick, drawers) => {
  const prisoners = drawers.length, k = Math.floor(prisoners / 2);
  return range(prisoners).every(prisoner => pick(prisoner, k, drawers));
  //return range(prisoners).every(prisoner => pick(prisoner, k, guard(drawers, k)));
};
const simulate = (prisoners, times) => {
  const drawers = range(prisoners);
  let successRandom = 0, successCyclic = 0;
  for (let i = 0; i < times; i++) {
    shuffle(drawers);
    successRandom += challenge(randomPick, drawers);
    successCyclic += challenge(cyclicPick, drawers);
    //console.assert(cyclicMax(drawers) <= Math.floor(prisoners / 2) === challenge(cyclicPick, drawers));
  }
  return {successRandom, successCyclic};
};

{
  const prisoners = 100, times = 100000;
  const {successRandom, successCyclic} = simulate(prisoners, times);
  console.log(`${prisoners} Random Pick: ${successRandom}/${times}`); // => 0%
  console.log(`${prisoners} Cyclic Pick: ${successCyclic}/${times}`); // => 31%
}

//[NOTE]
// cyclic on 4-prisoners: 41.67%, 8-prisoners: 36.55%, 65536-prisoners: 30.68%
// limit [n -> inf] prisoners: 1 - log(2) = 0.30685...
```

</div></details>

# リンク

- Rosetta Codeの100 prisoners: http://rosettacode.org/wiki/100_prisoners