genkuroki/100 prisoners problem.ipynb

## 100 prisoners problem.ipynb
{
  "cells": [
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "# 百囚人問題\n\n黒木玄\n\n2017-09-03, 2019-02-24\n\n私は次の問題を『数学セミナー』誌2012.06号に掲載された岩沢宏和氏の「続確率パズルの迷宮」の連載第5回目で知った。ピーター・ウィンクラー著『とっておきの数学パズル』(日本評論社)にこの問題が載っているらしい。\n\n>問題：100人の囚人の名札が100個の別々の木箱のに入れられていて、それらの箱は、ある部屋のテーブルの上に1列に並べられている。囚人たちは、1人ひとりその部屋に連れていかれる。その際、最も多くて50個の箱まで中を覗いてよいのだが、その間、部屋のいっさいのものを原状のままにしておくとともに、部屋から出た後には他の囚人たちといっさい連絡をとってはならないことになっている。\n\n>囚人たちには、事前に作戦を立てる機会が与えられているが、実際、それは絶対に必要である。なぜなら、100人が全員、自分の名札を見つけてこない限り、後で全員が処刑されてしまうからである。\n\n>囚人たちの助かる確率が30%を超えるような作戦を考え出してほしい。\n\nウィンクラー氏によるコメント\n\n>囚人たちがそれぞれ無作為に50個の箱を覗いていったら、彼らが生き残れる確率は、嫌になるほど低い\n\n>$$(1/2)^{100} \\risingdotseq 0.0000000000000000000000000000008$$\n\n>という値である。下手をすればもっとひどいことになる――たとえば、全員が同じ50個の箱を見ることにしておいたら、生き残れる確率は0である。だから、30%を超えようなどというのはとてもばかげた相談に見えるかもしれない――そう、その通りだが、そう見えるのはむしろ問題を正しく理解している証拠だ。その無理難題に挑戦してほしい。\n\n以上は『数学セミナー』誌2012.06号に掲載された岩沢宏和氏の記事からの引用である。岩沢氏の記事では確率の上限がどうなるかに関する証明は書いてあるが、この問題の答えは載っていない。\n\n上のパズルは百囚人問題(100 prisoners problem)となどと呼ばれ、wikipedia に項目が設けられるほど有名になっている。\n\nhttps://en.wikipedia.org/wiki/100_prisoners_problem\n\nこのパズルはとてもよい数学パズルだと思ったので、私は mathtodon で紹介することにした。\n\nhttps://mathtod.online/@genkuroki/563079\n\nこのノートブックの目的は、Julia言語を使って、百囚人問題関係の計算をしてみせることである。\n\n**2019年2月24日補足:** 以上で紹介した問題は連載をまとめた次の本のp.156以降でも紹介されている.\n\n* 岩沢宏和, <a href=\"https://www.amazon.co.jp/dp/4535787476\">確率パズルの迷宮</a>, 日本評論社, 2014, vi+288 pages."
    },
    {
      "metadata": {
        "toc": true
      },
      "cell_type": "markdown",
      "source": "<h1>目次<span class=\"tocSkip\"></span></h1>\n<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#全員が生き残れる確率\" data-toc-modified-id=\"全員が生き残れる確率-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>全員が生き残れる確率</a></span></li><li><span><a href=\"#置換を扱うパッケージを使ってみる\" data-toc-modified-id=\"置換を扱うパッケージを使ってみる-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>置換を扱うパッケージを使ってみる</a></span></li><li><span><a href=\"#生き残りのモンテカルロシミュレーション\" data-toc-modified-id=\"生き残りのモンテカルロシミュレーション-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>生き残りのモンテカルロシミュレーション</a></span></li><li><span><a href=\"#確率の計算\" data-toc-modified-id=\"確率の計算-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>確率の計算</a></span><ul class=\"toc-item\"><li><span><a href=\"#補題\" data-toc-modified-id=\"補題-4.1\"><span class=\"toc-item-num\">4.1&nbsp;&nbsp;</span>補題</a></span></li><li><span><a href=\"#生き残れる確率\" data-toc-modified-id=\"生き残れる確率-4.2\"><span class=\"toc-item-num\">4.2&nbsp;&nbsp;</span>生き残れる確率</a></span></li></ul></li></ul></div>"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "maxnworkers = 8\nusing Distributed\nif nworkers() < maxnworkers\n    @show addprocs(maxnworkers - nworkers())\nend\n@show workers()\n\n@everywhere using Permutations\n@everywhere using Permutations: randperm",
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "text": "addprocs(maxnworkers - nworkers()) = [2, 3, 4, 5, 6, 7, 8]\nworkers() = [2, 3, 4, 5, 6, 7, 8]\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## 全員が生き残れる確率\n\n実は生き残れる確率は\n\n$$\\displaystyle\n1 - \\sum_{k=51}^{100} \\frac{1}{k} = 0.3118278206898\\cdots\n$$\n\nになる。この確率で生き残れるようにする方法は後述。"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "n = 50\n@time 1 - sum(k->1/k,(n+1):2n)",
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": "  0.145020 seconds (284.53 k allocations: 14.394 MiB, 5.77% gc time)\n",
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "execution_count": 2,
          "data": {
            "text/plain": "0.31182782068980486"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "全体の人数が $2n=2\\times 10^9$ の同様の問題でも全員が生き残れる確率はそう変わらない。"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "n = 10^9\n@time 1 - sum(k->1/k,(n+1):2n)",
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": "  2.509592 seconds (95.34 k allocations: 4.844 MiB)\n",
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "execution_count": 3,
          "data": {
            "text/plain": "0.30685281969005473"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "$n\\to\\infty$ での極限で全員が生き残れる確率は $1-\\log 2$ になる。"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "1 - log(2)",
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 4,
          "data": {
            "text/plain": "0.3068528194400547"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## 置換を扱うパッケージを使ってみる\n\n$N$ 次の置換とは数 $1,2,\\ldots,N$ の各々をそれらのどれかに重複無しに対応させる函数(写像)のことである。例えば $10$ 次の置換 $\\sigma$ を\n\n$$\\begin{aligned}\n&\n\\sigma(1)=9, \\quad\n\\sigma(2)=4, \\quad\n\\sigma(3)=3, \\quad\n\\sigma(4)=6, \\quad\n\\sigma(5)=1, \\quad\n\\\\ &\n\\sigma(6)=5, \\quad\n\\sigma(7)=8, \\quad\n\\sigma(8)=10, \\quad\n\\sigma(9)=2, \\quad\n\\sigma(10)=7\n\\end{aligned}$$\n\nと定めることができる。数学の教科書ではこれをよく\n\n$$\n\\sigma = \\begin{pmatrix}\n1 & 2 & 3 & 4 & 5 & 6 & 7 & 8  & 9 & 10 \\\\\n9 & 4 & 3 & 6 & 1 & 5 & 8 & 10 & 2 & 7 \\\\\n\\end{pmatrix}\n$$\n\nのように書いたりする。これは上段の数を下段の数に対応させる函数という意味である。\n\nJulia言語では配列で表わされた置換をランダムに生成する函数として、 `randperm(N)` が用意されている。上の $\\sigma$ は配列によって\n\n$$\n[9, 4, 3, 6, 1, 5, 8, 10, 2, 7]\n$$\n\nと表わされる。\n\n任意の置換は以下のような方法によって巡回置換(その意味は以下の説明を読めばわかるだろう)の積で表わすことができる。例として上の置換 $\\sigma$ を用いよう。\n\n$\\sigma$ は1を9に移し、さらに9を2に移し、2を4に移す。これを次々に繰り返すと1から10までの数は10個しかないので、いつかは1に戻る。$\\sigma$ を1から10までの数に次々に作用させた結果を追跡すると、$\\sigma$ は\n\n    1→9→2→4→6→5→1\n    3→3\n    7→8→10→7\n\nと数を移すことがわかる。1→9→2→4→6→5→1のように数を巡回的に移し、他の数を動かさない置換を (1,9,2,4,6,5) と書き、長さ6の巡回置換と呼ぶ。ループの先頭はどれに選んでも巡回置換(写像)としては同じものになる。たとえば\n\n$$(1,9,2,4,6,5) = (2,4,6,5,1,9).$$\n\n他のループについても同様で (3) は長さ1の巡回置換を表わし、(7,8,10)は長さ3の巡回置換を表わす。これらの写像としての合成(この場合の合成の順番は可換)によって $\\sigma$ は\n\n$$\n\\sigma = (1,9,2,4,6,5)(3)(7,8,10)\n$$\n\nと表わされる。このとき、置換 $\\sigma$ は長さ1,3,6の巡回置換を含むと言うことにしよう。\n\nPermutations パッケージを使うと、`rp = randperm(N)` で生成した配列としての置換を `p = Permutation(rp)` によって巡回置換の積に分解された置換 `p` に変換できる。そして、置換 `p` に含まれる巡回置換たちを `c = cycles(p)` によって配列の配列として取り出せる。\n\n次の問題について考えよう。\n\n>$N$ 次の置換 $\\sigma$ をランダムに等確率で生成するとき、$\\sigma$ に含まれる巡回置換の長さの最大値が $N/2$ 以下になる確率を求めよ。\n\n以下ではJulia言語で書かれた `maxlensim(niters; N=100)` 函数によってモンテカルロシミュレーションを実行している。理論的に $N=100$ の場合にその確率は\n\n$$\\displaystyle\n1 - \\sum_{k=51}^{100} \\frac{1}{k} = 0.3118278206898\\cdots\n$$\n\nになり、モンテカルロシミュレーションの結果も概ねこの値になっていることを確認できるはずである。"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "@show rp = randperm(10);",
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "text": "rp = randperm(10) = [8, 10, 6, 1, 9, 2, 3, 4, 5, 7]\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "p = Permutation(rp)",
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 6,
          "data": {
            "text/plain": "(1,8,4)(2,10,7,3,6)(5,9)"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "c = cycles(p)",
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 7,
          "data": {
            "text/plain": "3-element Array{Array{Int64,1},1}:\n [1, 8, 4]       \n [2, 10, 7, 3, 6]\n [5, 9]          "
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "maximum(length.(c))",
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 8,
          "data": {
            "text/plain": "5"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "function maxlensim(niters; N=100)\n    n = div(N,2)\n    s = @distributed (+) for i in 1:niters\n        ifelse(maximum(length.(cycles(Permutation(randperm(N))))) <= n, 1, 0)\n    end\n    return s\nend",
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 9,
          "data": {
            "text/plain": "maxlensim (generic function with 1 method)"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "niters = 10^6\n@time maxlensim(niters)/niters",
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "text": "  5.570878 seconds (1.94 M allocations: 98.949 MiB, 0.42% gc time)\n",
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "execution_count": 10,
          "data": {
            "text/plain": "0.312074"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## 生き残りのモンテカルロシミュレーション\n\nランダムに生成された $N$ 次の置換が $N/2$ より長い巡回置換を含まない確率は30%より大きい。\n\nこのことを使えば、囚人たちが助かる確率を30%より大きくできる。\n\n囚人たちは全員に1から100までの番号を重複無しにランダムに割り振り、囚人たちは一直線に並んでいる木箱たちに1から100までの番号を重複無しにランダムに割り振る。これによって、囚人たちは木箱の中に名札が入っている状況を100次の置換とみなすことができる。そして、その置換に含まれる巡回置換の長さが50以下ならば、以下のようにして全員が生き残ることができる。\n\n$k$ 番の囚人は最初に $k$ 番の箱を最初に開け、出て来た名札が $l$ 番の囚人のものならばその次に $l$ 番の木箱を開け、同様に繰り返すことにする。$k$ 番の囚人が50箱以内に自分の名札を発見できればその囚人は成功したことになる。\n\n以下では、本当に30%より大きな確率で囚人たちが生き残ることができるかどうかをモンテカルロシミュレーションで確認している。"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "@everywhere function liveordie(card::Array{Int64,1})\n    N = length(card)\n    n = div(N,2)\n    c = 0\n    for i in 1:N\n        c = card[i]\n        k = 1\n        while c != i\n            c = card[c]\n            k += 1\n            if k > n\n                return false # DIE\n            end\n        end\n    end\n    return true # LIVE\nend",
      "execution_count": 11,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "function liveordiesim(niters; N=100)\n    s = @distributed (+) for i in 1:niters\n        ifelse(liveordie(randperm(N)), 1, 0)\n    end\n    return s\nend",
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 12,
          "data": {
            "text/plain": "liveordiesim (generic function with 1 method)"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "niters = 10^6\n@time liveordiesim(niters)/niters",
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "stream",
          "text": "  0.829628 seconds (119.57 k allocations: 5.886 MiB, 0.53% gc time)\n",
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "execution_count": 13,
          "data": {
            "text/plain": "0.310969"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "collapsed": true
      },
      "cell_type": "markdown",
      "source": "## 確率の計算"
    },
    {
      "metadata": {
        "collapsed": true
      },
      "cell_type": "markdown",
      "source": "### 補題\n\n$2k > N$ のとき $N$ 次の置換 $\\sigma$ が長さ $k$ の巡回置換を含む確率は $1/k$ である。\n\n**証明**　$2k > N$ であると仮定する。 $N$ 次の置換で長さ $k$ の巡回置換を含むもの全体の個数を $a_k$ と書くことにしよう。 $2k > N$ という仮定より、そのような置換が含む長さ $k$ の巡回置換は1つだけなので、そのような置換は長さ $k$ の巡回置換($N(N-1)\\cdots(N-k+1)/k$ 個存在する)とそこに含まれない数達の置換($(N-k)!$ 個存在する)の積で一意的に表わされる。$a_k$ はそれらの個数の積になるので、\n\n$$\na_k = N!/k.\n$$\n\nゆえに、$N$ 次の置換が長さ $k$ の巡回置換を含む確率は\n\n$$\na_k/N! = 1/k\n$$\n\nになる。 q.e.d."
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "### 生き残れる確率\n\n$2n$ 次の置換に含まれる巡回置換の長さの最大値が $n$ 以下になる確率は\n\n$$\\displaystyle\n1 - \\sum_{k=n+1}^{2n}\\frac{1}{k}\n$$\n\nである。\n\n**証明**　上の補題より、ただちに得られる。 q.e.d."
    }
  ],
  "metadata": {
    "_draft": {
      "nbviewer_url": "https://gist.github.com/fa348364bec065eccab2d7aa55dd929e"
    },
    "gist": {
      "id": "fa348364bec065eccab2d7aa55dd929e",
      "data": {
        "description": "Julia/Statistics/百囚人問題/100 prisoners problem.ipynb",
        "public": true
      }
    },
    "kernelspec": {
      "name": "julia-1.1",
      "display_name": "Julia 1.1.0",
      "language": "julia"
    },
    "language_info": {
      "file_extension": ".jl",
      "name": "julia",
      "mimetype": "application/julia",
      "version": "1.1.0"
    },
    "toc": {
      "nav_menu": {
        "height": "31px",
        "width": "252px"
      },
      "number_sections": true,
      "sideBar": true,
      "skip_h1_title": true,
      "title_cell": "目次",
      "title_sidebar": "目次",
      "toc_cell": true,
      "toc_position": {},
      "toc_section_display": "block",
      "toc_window_display": false
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}
	{
	"cells": [
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "# 百囚人問題\n\n黒木玄\n\n2017-09-03, 2019-02-24\n\n私は次の問題を『数学セミナー』誌2012.06号に掲載された岩沢宏和氏の「続確率パズルの迷宮」の連載第5回目で知った。ピーター・ウィンクラー著『とっておきの数学パズル』(日本評論社)にこの問題が載っているらしい。\n\n>問題：100人の囚人の名札が100個の別々の木箱のに入れられていて、それらの箱は、ある部屋のテーブルの上に1列に並べられている。囚人たちは、1人ひとりその部屋に連れていかれる。その際、最も多くて50個の箱まで中を覗いてよいのだが、その間、部屋のいっさいのものを原状のままにしておくとともに、部屋から出た後には他の囚人たちといっさい連絡をとってはならないことになっている。\n\n>囚人たちには、事前に作戦を立てる機会が与えられているが、実際、それは絶対に必要である。なぜなら、100人が全員、自分の名札を見つけてこない限り、後で全員が処刑されてしまうからである。\n\n>囚人たちの助かる確率が30%を超えるような作戦を考え出してほしい。\n\nウィンクラー氏によるコメント\n\n>囚人たちがそれぞれ無作為に50個の箱を覗いていったら、彼らが生き残れる確率は、嫌になるほど低い\n\n>$$(1/2)^{100} \\risingdotseq 0.0000000000000000000000000000008$$\n\n>という値である。下手をすればもっとひどいことになる――たとえば、全員が同じ50個の箱を見ることにしておいたら、生き残れる確率は0である。だから、30%を超えようなどというのはとてもばかげた相談に見えるかもしれない――そう、その通りだが、そう見えるのはむしろ問題を正しく理解している証拠だ。その無理難題に挑戦してほしい。\n\n以上は『数学セミナー』誌2012.06号に掲載された岩沢宏和氏の記事からの引用である。岩沢氏の記事では確率の上限がどうなるかに関する証明は書いてあるが、この問題の答えは載っていない。\n\n上のパズルは百囚人問題(100 prisoners problem)となどと呼ばれ、wikipedia に項目が設けられるほど有名になっている。\n\nhttps://en.wikipedia.org/wiki/100_prisoners_problem\n\nこのパズルはとてもよい数学パズルだと思ったので、私は mathtodon で紹介することにした。\n\nhttps://mathtod.online/@genkuroki/563079\n\nこのノートブックの目的は、Julia言語を使って、百囚人問題関係の計算をしてみせることである。\n\n2019年2月24日補足: 以上で紹介した問題は連載をまとめた次の本のp.156以降でも紹介されている.\n\n* 岩沢宏和, <a href=\"https://www.amazon.co.jp/dp/4535787476\">確率パズルの迷宮</a>, 日本評論社, 2014, vi+288 pages."
	},
	{
	"metadata": {
	"toc": true
	},
	"cell_type": "markdown",
	"source": "<h1>目次<span class=\"tocSkip\"></span></h1>\n<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#全員が生き残れる確率\" data-toc-modified-id=\"全員が生き残れる確率-1\"><span class=\"toc-item-num\">1  </span>全員が生き残れる確率</a></span></li><li><span><a href=\"#置換を扱うパッケージを使ってみる\" data-toc-modified-id=\"置換を扱うパッケージを使ってみる-2\"><span class=\"toc-item-num\">2  </span>置換を扱うパッケージを使ってみる</a></span></li><li><span><a href=\"#生き残りのモンテカルロシミュレーション\" data-toc-modified-id=\"生き残りのモンテカルロシミュレーション-3\"><span class=\"toc-item-num\">3  </span>生き残りのモンテカルロシミュレーション</a></span></li><li><span><a href=\"#確率の計算\" data-toc-modified-id=\"確率の計算-4\"><span class=\"toc-item-num\">4  </span>確率の計算</a></span><ul class=\"toc-item\"><li><span><a href=\"#補題\" data-toc-modified-id=\"補題-4.1\"><span class=\"toc-item-num\">4.1  </span>補題</a></span></li><li><span><a href=\"#生き残れる確率\" data-toc-modified-id=\"生き残れる確率-4.2\"><span class=\"toc-item-num\">4.2  </span>生き残れる確率</a></span></li></ul></li></ul></div>"
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "maxnworkers = 8\nusing Distributed\nif nworkers() < maxnworkers\n @show addprocs(maxnworkers - nworkers())\nend\n@show workers()\n\n@everywhere using Permutations\n@everywhere using Permutations: randperm",
	"execution_count": 1,
	"outputs": [
	{
	"output_type": "stream",
	"text": "addprocs(maxnworkers - nworkers()) = [2, 3, 4, 5, 6, 7, 8]\nworkers() = [2, 3, 4, 5, 6, 7, 8]\n",
	"name": "stdout"
	}
	]
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "## 全員が生き残れる確率\n\n実は生き残れる確率は\n\n$$\\displaystyle\n1 - \\sum_{k=51}^{100} \\frac{1}{k} = 0.3118278206898\\cdots\n$$\n\nになる。この確率で生き残れるようにする方法は後述。"
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "n = 50\n@time 1 - sum(k->1/k,(n+1):2n)",
	"execution_count": 2,
	"outputs": [
	{
	"output_type": "stream",
	"text": " 0.145020 seconds (284.53 k allocations: 14.394 MiB, 5.77% gc time)\n",
	"name": "stdout"
	},
	{
	"output_type": "execute_result",
	"execution_count": 2,
	"data": {
	"text/plain": "0.31182782068980486"
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "全体の人数が $2n=2\\times 10^9$ の同様の問題でも全員が生き残れる確率はそう変わらない。"
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "n = 10^9\n@time 1 - sum(k->1/k,(n+1):2n)",
	"execution_count": 3,
	"outputs": [
	{
	"output_type": "stream",
	"text": " 2.509592 seconds (95.34 k allocations: 4.844 MiB)\n",
	"name": "stdout"
	},
	{
	"output_type": "execute_result",
	"execution_count": 3,
	"data": {
	"text/plain": "0.30685281969005473"
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "$n\\to\\infty$ での極限で全員が生き残れる確率は $1-\\log 2$ になる。"
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "1 - log(2)",
	"execution_count": 4,
	"outputs": [
	{
	"output_type": "execute_result",
	"execution_count": 4,
	"data": {
	"text/plain": "0.3068528194400547"
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "## 置換を扱うパッケージを使ってみる\n\n$N$ 次の置換とは数 $1,2,\\ldots,N$ の各々をそれらのどれかに重複無しに対応させる函数(写像)のことである。例えば $10$ 次の置換 $\\sigma$ を\n\n$$\\begin{aligned}\n&\n\\sigma(1)=9, \\quad\n\\sigma(2)=4, \\quad\n\\sigma(3)=3, \\quad\n\\sigma(4)=6, \\quad\n\\sigma(5)=1, \\quad\n\\\\ &\n\\sigma(6)=5, \\quad\n\\sigma(7)=8, \\quad\n\\sigma(8)=10, \\quad\n\\sigma(9)=2, \\quad\n\\sigma(10)=7\n\\end{aligned}$$\n\nと定めることができる。数学の教科書ではこれをよく\n\n$$\n\\sigma = \\begin{pmatrix}\n1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\\\\n9 & 4 & 3 & 6 & 1 & 5 & 8 & 10 & 2 & 7 \\\\\n\\end{pmatrix}\n$$\n\nのように書いたりする。これは上段の数を下段の数に対応させる函数という意味である。\n\nJulia言語では配列で表わされた置換をランダムに生成する函数として、 `randperm(N)` が用意されている。上の $\\sigma$ は配列によって\n\n$$\n[9, 4, 3, 6, 1, 5, 8, 10, 2, 7]\n$$\n\nと表わされる。\n\n任意の置換は以下のような方法によって巡回置換(その意味は以下の説明を読めばわかるだろう)の積で表わすことができる。例として上の置換 $\\sigma$ を用いよう。\n\n$\\sigma$ は1を9に移し、さらに9を2に移し、2を4に移す。これを次々に繰り返すと1から10までの数は10個しかないので、いつかは1に戻る。$\\sigma$ を1から10までの数に次々に作用させた結果を追跡すると、$\\sigma$ は\n\n 1→9→2→4→6→5→1\n 3→3\n 7→8→10→7\n\nと数を移すことがわかる。1→9→2→4→6→5→1のように数を巡回的に移し、他の数を動かさない置換を (1,9,2,4,6,5) と書き、長さ6の巡回置換と呼ぶ。ループの先頭はどれに選んでも巡回置換(写像)としては同じものになる。たとえば\n\n$$(1,9,2,4,6,5) = (2,4,6,5,1,9).$$\n\n他のループについても同様で (3) は長さ1の巡回置換を表わし、(7,8,10)は長さ3の巡回置換を表わす。これらの写像としての合成(この場合の合成の順番は可換)によって $\\sigma$ は\n\n$$\n\\sigma = (1,9,2,4,6,5)(3)(7,8,10)\n$$\n\nと表わされる。このとき、置換 $\\sigma$ は長さ1,3,6の巡回置換を含むと言うことにしよう。\n\nPermutations パッケージを使うと、`rp = randperm(N)` で生成した配列としての置換を `p = Permutation(rp)` によって巡回置換の積に分解された置換 `p` に変換できる。そして、置換 `p` に含まれる巡回置換たちを `c = cycles(p)` によって配列の配列として取り出せる。\n\n次の問題について考えよう。\n\n>$N$ 次の置換 $\\sigma$ をランダムに等確率で生成するとき、$\\sigma$ に含まれる巡回置換の長さの最大値が $N/2$ 以下になる確率を求めよ。\n\n以下ではJulia言語で書かれた `maxlensim(niters; N=100)` 函数によってモンテカルロシミュレーションを実行している。理論的に $N=100$ の場合にその確率は\n\n$$\\displaystyle\n1 - \\sum_{k=51}^{100} \\frac{1}{k} = 0.3118278206898\\cdots\n$$\n\nになり、モンテカルロシミュレーションの結果も概ねこの値になっていることを確認できるはずである。"
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "@show rp = randperm(10);",
	"execution_count": 5,
	"outputs": [
	{
	"output_type": "stream",
	"text": "rp = randperm(10) = [8, 10, 6, 1, 9, 2, 3, 4, 5, 7]\n",
	"name": "stdout"
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "p = Permutation(rp)",
	"execution_count": 6,
	"outputs": [
	{
	"output_type": "execute_result",
	"execution_count": 6,
	"data": {
	"text/plain": "(1,8,4)(2,10,7,3,6)(5,9)"
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "c = cycles(p)",
	"execution_count": 7,
	"outputs": [
	{
	"output_type": "execute_result",
	"execution_count": 7,
	"data": {
	"text/plain": "3-element Array{Array{Int64,1},1}:\n [1, 8, 4] \n [2, 10, 7, 3, 6]\n [5, 9] "
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "maximum(length.(c))",
	"execution_count": 8,
	"outputs": [
	{
	"output_type": "execute_result",
	"execution_count": 8,
	"data": {
	"text/plain": "5"
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "function maxlensim(niters; N=100)\n n = div(N,2)\n s = @distributed (+) for i in 1:niters\n ifelse(maximum(length.(cycles(Permutation(randperm(N))))) <= n, 1, 0)\n end\n return s\nend",
	"execution_count": 9,
	"outputs": [
	{
	"output_type": "execute_result",
	"execution_count": 9,
	"data": {
	"text/plain": "maxlensim (generic function with 1 method)"
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "niters = 10^6\n@time maxlensim(niters)/niters",
	"execution_count": 10,
	"outputs": [
	{
	"output_type": "stream",
	"text": " 5.570878 seconds (1.94 M allocations: 98.949 MiB, 0.42% gc time)\n",
	"name": "stdout"
	},
	{
	"output_type": "execute_result",
	"execution_count": 10,
	"data": {
	"text/plain": "0.312074"
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "## 生き残りのモンテカルロシミュレーション\n\nランダムに生成された $N$ 次の置換が $N/2$ より長い巡回置換を含まない確率は30%より大きい。\n\nこのことを使えば、囚人たちが助かる確率を30%より大きくできる。\n\n囚人たちは全員に1から100までの番号を重複無しにランダムに割り振り、囚人たちは一直線に並んでいる木箱たちに1から100までの番号を重複無しにランダムに割り振る。これによって、囚人たちは木箱の中に名札が入っている状況を100次の置換とみなすことができる。そして、その置換に含まれる巡回置換の長さが50以下ならば、以下のようにして全員が生き残ることができる。\n\n$k$ 番の囚人は最初に $k$ 番の箱を最初に開け、出て来た名札が $l$ 番の囚人のものならばその次に $l$ 番の木箱を開け、同様に繰り返すことにする。$k$ 番の囚人が50箱以内に自分の名札を発見できればその囚人は成功したことになる。\n\n以下では、本当に30%より大きな確率で囚人たちが生き残ることができるかどうかをモンテカルロシミュレーションで確認している。"
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "@everywhere function liveordie(card::Array{Int64,1})\n N = length(card)\n n = div(N,2)\n c = 0\n for i in 1:N\n c = card[i]\n k = 1\n while c != i\n c = card[c]\n k += 1\n if k > n\n return false # DIE\n end\n end\n end\n return true # LIVE\nend",
	"execution_count": 11,
	"outputs": []
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "function liveordiesim(niters; N=100)\n s = @distributed (+) for i in 1:niters\n ifelse(liveordie(randperm(N)), 1, 0)\n end\n return s\nend",
	"execution_count": 12,
	"outputs": [
	{
	"output_type": "execute_result",
	"execution_count": 12,
	"data": {
	"text/plain": "liveordiesim (generic function with 1 method)"
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "niters = 10^6\n@time liveordiesim(niters)/niters",
	"execution_count": 13,
	"outputs": [
	{
	"output_type": "stream",
	"text": " 0.829628 seconds (119.57 k allocations: 5.886 MiB, 0.53% gc time)\n",
	"name": "stdout"
	},
	{
	"output_type": "execute_result",
	"execution_count": 13,
	"data": {
	"text/plain": "0.310969"
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"collapsed": true
	},
	"cell_type": "markdown",
	"source": "## 確率の計算"
	},
	{
	"metadata": {
	"collapsed": true
	},
	"cell_type": "markdown",
	"source": "### 補題\n\n$2k > N$ のとき $N$ 次の置換 $\\sigma$ が長さ $k$ の巡回置換を含む確率は $1/k$ である。\n\n証明　$2k > N$ であると仮定する。 $N$ 次の置換で長さ $k$ の巡回置換を含むもの全体の個数を $a_k$ と書くことにしよう。 $2k > N$ という仮定より、そのような置換が含む長さ $k$ の巡回置換は1つだけなので、そのような置換は長さ $k$ の巡回置換($N(N-1)\\cdots(N-k+1)/k$ 個存在する)とそこに含まれない数達の置換($(N-k)!$ 個存在する)の積で一意的に表わされる。$a_k$ はそれらの個数の積になるので、\n\n$$\na_k = N!/k.\n$$\n\nゆえに、$N$ 次の置換が長さ $k$ の巡回置換を含む確率は\n\n$$\na_k/N! = 1/k\n$$\n\nになる。 q.e.d."
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "### 生き残れる確率\n\n$2n$ 次の置換に含まれる巡回置換の長さの最大値が $n$ 以下になる確率は\n\n$$\\displaystyle\n1 - \\sum_{k=n+1}^{2n}\\frac{1}{k}\n$$\n\nである。\n\n証明　上の補題より、ただちに得られる。 q.e.d."
	}
	],
	"metadata": {
	"_draft": {
	"nbviewer_url": "https://gist.github.com/fa348364bec065eccab2d7aa55dd929e"
	},
	"gist": {
	"id": "fa348364bec065eccab2d7aa55dd929e",
	"data": {
	"description": "Julia/Statistics/百囚人問題/100 prisoners problem.ipynb",
	"public": true
	}
	},
	"kernelspec": {
	"name": "julia-1.1",
	"display_name": "Julia 1.1.0",
	"language": "julia"
	},
	"language_info": {
	"file_extension": ".jl",
	"name": "julia",
	"mimetype": "application/julia",
	"version": "1.1.0"
	},
	"toc": {
	"nav_menu": {
	"height": "31px",
	"width": "252px"
	},
	"number_sections": true,
	"sideBar": true,
	"skip_h1_title": true,
	"title_cell": "目次",
	"title_sidebar": "目次",
	"toc_cell": true,
	"toc_position": {},
	"toc_section_display": "block",
	"toc_window_display": false
	}
	},
	"nbformat": 4,
	"nbformat_minor": 2
	}