確率とモンテカルロシミュレーション

問題
Aさんには三人の子供がいる
「日曜日生まれの男の子はいるか？」と訊くと「いる」と答えた
Aさんの子供が三人とも男の子である確率は？

最近２ちゃんまとめで上の問題やってみて、
これをモンテカルロでやってみたら結果は理論解に収束するのかな？
って思ったのでやってみた

　モンテカルロシミュレーションとは

シミュレーションや数値計算を乱数を用いて行う手法の総称。元々は、中性子が物質中を動き回る様子を探るために考案された。
強みは基礎的な理論さえ把握していれば、一定のアプローチに従って幅広い問題を解決できるという点

by Wiki

まぁ簡単にいえば、乱数つかって事象を表現して、それをループ文で何回も試し大数の法則でまとめ上げる手法

例：コインを投げた時の裏表の確率
表=1,裏=0として、0か1の乱数を生成。これを10万回(試行回数)繰り返し、1(表)が出たときの数をカウントしておく
コインが表になる確率は、「カウント数 / 10万(試行回数)」で求まる。

　理論解（ベイズの定理）

初見で1/4って思って恥じ書いたｗｗ三人のうち、一人が男確定だから、残り二人が男の確立Pは下の通り \[ P=\frac{1}{2}\cdot\frac{1}{2} =\frac{1}{4} \] (￣ー￣) ﾄﾞﾔｯ!
はい。日曜日を軽視してましたーｗｗ

これはどう見ても事後確率だから、ベイズの定理を使って解いてくのが定石

事象$A$：3人とも男が生まれる
事象$B$：少なくとも1人は日曜日に生まれの男がいる
　　とおくと、求める確率$P(A|B)$はベイズの定理より \[ P(A|B)=\frac{P(A \land B)}{P(B)} \]
$P(B)$は「3人とも日曜生まれの男でない」場合の余事象なので \[ P(B)=1-\frac{ 13^3}{14^3}=\frac{547}{2744} \]
$P(A \land B)$は3人とも男で、うち少なくとも一人は日曜日に生まれてる \[ P(A \land B)= _3C_0\left(\frac{1}{14}\right)^3\left(\frac{6}{14}\right)^0+ _3C_1\left(\frac{1}{14}\right)^2\left(\frac{6}{14}\right)^1+ _3C_2\left(\frac{1}{14}\right)^1\left(\frac{6}{14}\right)^2\\ =\frac{127}{2744} \] よって \[ P(A|B)=\frac{P(A \land B)}{P(B)}=\frac{127}{547} \]

　数値解（モンテカルロ法）

今回はCで書いてみた
誕生曜日、性別を持つ構造体をつくり、事象Bを満たす三兄弟を乱数使って生成
三兄弟が全員男である場合をカウントしていく。
最後に、試行回数(ループ数)とカウンターの比率がP(A|B)の確率となる

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

#define SEX_BOY 0
#define SEX_GIRL 1

#define WEEK_SUNDAY    0
#define WEEK_MONDAY    1
#define WEEK_TUESDAY   2
#define WEEK_WEDNESDAY 3
#define WEEK_THURSDAY  4
#define WEEK_FRIDAY    5
#define WEEK_SATURDAY  6

#define LOOP_MAX 10000000

typedef struct {
    int sex;
    int week;
} human_t;

main() {
    srand((unsigned)time(NULL));
    long unsigned loop_count, check_count = 0;
    const human_t target_child = {SEX_BOY, WEEK_SUNDAY};
    human_t children[3];

    for (loop_count = 0; loop_count < LOOP_MAX; loop_count++) {

        // 事象B：少なくとも一人日曜生まれの男がいる条件を成立させる
        int check_flg = -1;
        while(check_flg != 0) {
            int i;
            for (i = 0; i < 3; i++) {
                human_t child = {rand() % 2, rand() % 7};
                children[i] = child;
                if (memcmp(&target_child, &child, sizeof(human_t)) == 0) {
                    check_flg = 0;
                }
            }
        }

        // 事象A：全員男の場合を調べる
        if( (children[0].sex == SEX_BOY)
              &&(children[1].sex == SEX_BOY)
              &&(children[2].sex == SEX_BOY)) {
            check_count++;
        }
    }
    double reslt = check_count/(double)LOOP_MAX;
    printf("reslt = %f",reslt);
}

結果は

reslt = 0.232030

理論解を小数点表記にすると、　127/547 = 0.232175503
よって、誤差は0.0006
いい感じ（＾ｑ＾）

★ちなみに★

事象Bのロジックについて。。。。
ソースの31～42行目で生成してる三兄弟は
三兄弟を乱数で生成 → その中に日曜生まれ男がいるかチェック → いれば三兄弟を指標として採用。いなければ作り直し

// 事象B：少なくとも一人日曜生まれの男がいる条件を成立させる
int check_flg = -1;
while(check_flg != 0) {
    int i;
    for (i = 0; i < 3; i++) {
        human_t child = {rand() % 2, rand() % 7};
        children[i] = child;
        if (memcmp(&target_child, &child, sizeof(human_t)) == 0) {
            check_flg = 0;
        }
    }
}

となってる。
これが面倒で
三兄弟を乱数で生成、その後、兄弟の誰かを日曜生まれ男に差し替える

// 事象B：少なくとも一人日曜生まれの男がいる条件を成立させる
human_t children[3] = {
    {rand()%2, rand()%7},
    {rand()%2, rand()%7},
    {rand()%2, rand()%7}
};
// 三人のうち誰かを日曜生まれ男に差し替え
children[rand()%3] = target_child;

とか
長男を日曜生まれ男で固定し、残り二人を乱数で生成

// 事象B：少なくとも一人日曜生まれの男がいる条件を成立させる
human_t children[3] = {
    {SEX_BOY, WEEK_SUNDAY}, // 日曜男
    {rand()%2, rand()%7},
    {rand()%2, rand()%7}
};

などをすると、数値解析結果は

• 三兄弟を乱数で生成、その後、兄弟の誰かを日曜生まれ男に差し替える
reslt = 0.249997
• 長男を日曜生まれ男で固定し、残り二人を乱数で生成
reslt = 0.249940

僕が理論解でミスった値、1/4に収束していることがわかる
これは事後確率であることがプログラム上で示せていないため
つまり、ちょっと難しいが、モンテカルロシミュレーションでも事後確率と普通の確率は区別して求めることができる！

　まとめ

モンテカルロシミュレーションで事後確率を示すには、
前提条件の定義さえしっかりプログラムできていれば求めることができる！