モンモール問題

$n$名でプレゼント交換を一回行ったとき、誰一人も自分のプレゼントに当たらない確率$Pn$は？

P_{n}=1-\frac{1}{1 !}+\frac{1}{2 !}-\frac{1}{3 !}+\cdots+(-1)^{n} \frac{1}{n !}=\sum_{k=0}^{n} \frac{(-1)^{k}}{k !}

いきなりですが上式が答えになります。何かに似ていますね。
そうですテイラー近似式です。

e^{x} \fallingdotseq 1+\frac{1}{1 !} x+\frac{1}{2 !} x^{2}+\frac{1}{3 !} x^{3}+\cdots+\frac{1}{n !} x^{n}

$x=-1$ を代入すると近似できます。

P_{n} \fallingdotseq e^{-1}=0.36787944 \cdots \fallingdotseq 0.3679

実際に$n$に値を入れてみると次のようになります。

$n$	1	2	3	4	5	6	7	$\cdots$	10	50	100
$P_{n}$	0	0.5	0.3333	0.375	0.3667	0.3681	0.3679	$\cdots$	0.3679	0.3679	0.3679

$n \geqq 5$ ですでに四捨五入すれば$37\%$です。なかなか速く収束するようです。
とここまで偉そうに書いてきましたが$P_{n}$の導出がないのを不思議に思った方がおられると思います。
鋭いです。私は$P_{n}$がなぜ一番上の式で求められるかわからないからです😭~~（プリント丸写しだし）~~。
また$P_{n}$の式の証明が東工大の入試で出たこともあるようです~~（言い訳）~~。

ここからが本題です。
理論がわからなくても模擬実験で確かめることができます。
シミュレーションってやつですね。
というこでモンテカルロ法を使って本当に$e^{-1}$で近似できるのかを確かめます。

モンテカルロ法

私も全然詳しくないのですが軽く説明させてもらうと
計算機で疑似乱数を用いて何回もシミュレートすれば大数の法則で実際の値に近いデータを得られるよってことです。たぶん
ざっくりですみません。
円周率を求めるのが有名ですね。

では実際のコードを見てみましょう。言語はCです。

montmort_montecarlo.c

/* モンモール問題
n名でプレゼント交換を一回行ったとき,誰一人も自分のプレゼントに当たらない確率Pnは？
この問題はnが十分大きいとき1/e(約0.3679)に収束する。
それをモンテカルロ法で確かめる。
*/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define NUM_TRI 100000000  // 試行回数

void shuffle(int ary[], int size)  // 配列をごちゃまぜにする
{
    for (int i=0; i<size; i++) {
        int j = rand()%size;
        int tmp = ary[i];
        ary[i] = ary[j];
        ary[j] = tmp;
    }
}

double monte_carlo(int n)  // nは参加人数
{
    int failure = 0;
    int present[n];
    for (int i=0; i<n; i++) {
        present[i] = i;
    }

    for (int i=0; i<NUM_TRI; i++) {
        shuffle(present, n);
        for (int j=0; j<n; j++) {
            if (present[j] == j) {   // 自分のプレゼントが
                failure++;           // 自分のとこに来たので
                break;               // 失敗
            }
        }
    }
    return 1 - (double)failure / NUM_TRI;  // 余事象で求める
}

int main(void)
{
    srand((unsigned int)time(NULL));

    printf("n=1の時   :Pn=%f\n", monte_carlo(1));
    printf("n=2の時   :Pn=%f\n", monte_carlo(2));
    printf("n=3の時   :Pn=%f\n", monte_carlo(3));
    printf("n=4の時   :Pn=%f\n", monte_carlo(4));
    printf("n=5の時   :Pn=%f\n", monte_carlo(5));
    printf("n=6の時   :Pn=%f\n", monte_carlo(6));
    printf("n=7の時   :Pn=%f\n", monte_carlo(7));
    printf("n=8の時   :Pn=%f\n", monte_carlo(8));
    printf("n=9の時   :Pn=%f\n", monte_carlo(9));
    printf("n=10の時  :Pn=%f\n", monte_carlo(10));
    printf("n=50の時  :Pn=%f\n", monte_carlo(50));
    printf("n=100の時 :Pn=%f\n", monte_carlo(100));

    return 0;
}

実行結果

$ ./monmort_montecarlo
n=1の時   :Pn=0.000000
n=2の時   :Pn=0.499987
n=3の時   :Pn=0.333385
n=4の時   :Pn=0.375038
n=5の時   :Pn=0.366645
n=6の時   :Pn=0.368129
n=7の時   :Pn=0.367885
n=8の時   :Pn=0.367863
n=9の時   :Pn=0.367930
n=10の時  :Pn=0.367870
n=50の時  :Pn=0.367992
n=100の時 :Pn=0.367885

どうやら$e^{-1}$に近づくのは本当らしいですね。

モンモール問題をモンテカルロ法で確かめる（C言語）

モンモール問題

モンテカルロ法