Chainerチュートリアル問6.4(ベイズの定理)の解法3選と解説

問6.4(ベイズの定理)

一様な確率で1から6までの数が出るサイコロがあります。 サイコロを何回か連続で振る試行を考えます。サイコロを振った後、直前に出た数より大きい数が出たときに「ラッキーである」と呼ぶことにします。 例えばサイコロを5回振って 2,5,3,4,4 と順に出た場合、2回目と4回目はラッキーですが、それ以外はラッキーではありません。
いま、サイコロを 10 回振って一度もラッキーではなかったとします。このとき、次にサイコロを振ってラッキーになる確率はいくらになるでしょうか。プログラムを使って分析してみて下さい。

Chainerチュートリアル https://tutorials.chainer.org/ja/Exercise_Step_01.html より引用

私は、プログラミング初心者であり、この問題に対するアプローチが全くわからなかった。
そこで情報を調べてみると、解答が書かれたサイトを発見。
コードの読解にかなり時間と体力を要したので、備忘録として記すことにする。

解法1:modとn進法の考え方を使った方法

import time
StartTime = time.time()

n = 11 # (n-1)回目まで全てアンラッキーで、n回目でラッキーになる
n1_unlucky = 0
n_lucky = 0
deme =[0] * n # 出目のリストを用意

for i in range(6 ** n):
    temp = i # temp = 一時的な値
    for j in range(n):
        deme[j] = temp % 6 + 1 # *1
        if 0 < j < n -1 and deme[j] > deme[j-1]:#ラッキーが途中で出たらストップ
            break
        temp = temp // 6 # *2

    if j == n - 1: # (n-1)回ラッキーが出なかった場合
        n1_unlucky += 1
        if deme[j] > deme[j-1]: # n回目がラッキー
            n_lucky += 1

print("Time={}".format(time.time() - StartTime))
print(n_lucky / n1_unlucky)
print(n_lucky,'/',n1_unlucky)

from fractions import Fraction
print(Fraction(n_lucky,n1_unlucky))

terarail質問「chainer チュートリアル　ベイズの定理」 https://teratail.com/questions/267662 より引用(一部書き換えあり)

解説

まず初期値である
・n1_unlucky
・n_lucky
・deme
を設定(demeはリスト形式で)。

次に、for文で0～6のn乗-1、つまり全ての場合の数で回している。
tempという一時的な値を起き、サイコロを投げる回数nでまたfor文を回す。

*1と*2の意味を理解するのに1時間ほどかかったが、これは余りによる分類(mod)と桁数を１つ下げる(n進法)の考え方を使い、MECE(もれなくダブりなく)を実現している。
どういうことか、わかりやすいように10進法で説明する。

*1と*2のわかりやすい説明

10進法が最もわかりやすいため、6面サイコロを10面サイコロに変え、計算式中の6もすべて10に変えて説明する。
ここで、たとえば temp = 2345という数に対する、jに対するfor文を考える。

まずj=0、つまりdeme[0]は、temp % 10 + 1 となり、10で割ったときの余り+1、つまり「一の位の数+1」を表す。
そしてtempに temp // 10、つまり10で割ったときの商を再代入するのだが、これは一の位を除外することを意味する(2345÷10=234余り5であるため)。deme[0]で使った一の位は、もう用無しだということだ。というよりここで用無しにしておかないと、あとでダブりが出てくる。

さて、ここまでくれば理解できるだろう。あとはこれらの操作を繰り返していくと、2345という数に対するdemeのリストは[6, 5, 4, 3, 1, 1, 1, 1, 1, 1, 1]ということになる(n=11とした)。
つまり、数とリストが1対1に全て対応するのだ。これでMECEが作り出せていることになる。

今10進法で考えたが、サイコロは6面なので、6進法バージョンにしてあげたのが、コードの意味になる。

その後は、分母となるn1_unluckyに1を加え、n回目がラッキーだった場合のみ分子にも1を加えれば、確率が求められる。

解法1の弱点

n=8ぐらいまでは処理が少なくてちゃんと機能するのだが、n=9ぐらいからかかる時間が増えていき、n=11にもなると処理が終わるまで約2時間ほどかかってしまう。
気になる人は、nに色んな値を代入し、ぜひ試してみてほしい。

そこで、計算量の問題を解決したのが、次の解法2である。

解法2:計算量を抑え処理を迅速にした方法

import time
StartTime = time.time()

from collections import defaultdict
# ラッキー判定
def lucky(p,c):
    return p < c

# 探索関数。引数はそれまでのサイコロの出目のリスト
def search(l):
    global cnt
    p = l[-1]  # リストの一番最後の出目
    # 10投して全部ノーラッキーだったら総数にカウントする
    # リストの要素数が10だったら条件成立
    if len(l) == N:
        #
        cnt += 1
        cnt_d[p] += 1
        return 
    # それ以外の場合は次のサイコロ試行。１～６
    for c in range(1,7):
        # ラッキーだったらそこで探索終了
        if lucky(p,c):
            pass
        else:
            # ノーラッキーだったら引数のリストに今の出目を加えて次の試行へ
            search(l + [c])

cnt = 0  # ノーラッキー総数
cnt_d = defaultdict(int)  # 最後の出目別ノーラッキー数
N = 10   # サイコロ投擲回数
total = 6 ** N  # 10回のサイコロ投擲総組み合わせ数

# 第1投目を初期値として与えて探索開始
for i in range(1,7):
    search([i])

#10回振って一度もラッキーしない確率
print(cnt,'/',total)
# うち最後の目別のカウント
print(cnt_d)

# 分母は10投ノーラッキー数ごとに11投目の出目数の乗算
fb = cnt * 6
# 分子はラッキー総数
fi = 0
# ラッキー数計算
for i in range(1,7):
    fi += (cnt_d[i]) * (6-i)

# 10投目までノーラッキーのとき、11投目にラッキーが出る確率
print("Time={}".format(time.time() - StartTime))
print(fi / fb)
print(fi,'/',fb)
from fractions import Fraction
print(Fraction(fi, fb))

terarail質問「chainer チュートリアル　ベイズの定理」 https://teratail.com/questions/267662 より引用(ほぼそのまま引用)

defaultdictを使い、10回振って一度もラッキーが起こらなかったときの最後の出目をカウントできる。これがこの解法の肝。

解説

たとえば最後の出目が1であれば、次にラッキーになるのは2,3,4,5,6の5通り。よって、最後の出目が1の場合の数に5を書ければ、場合の数が導き出せる。
最後の出目が2であれば、次にラッキーになるのは3,4,5,6の4通りなので、4をかける。
それが、fi += (cnt_d[i]) * (6-i) の意味である。

解法1では、for文をrange(6 ** n)の範囲で回していたので、計算量が莫大になってしまったが、解法2ではrange(1,7)で収まるように工夫している。
また、luckyだった場合はそこで探索を終了しているため、それでも計算量をかなり抑えられている。
N=10の場合でも、分子の場合の数は3003通りしかなく、分子の場合の数しか計算しないアルゴリズムになっているため、非常に高速に答えを求めることができるのだ。

実験してみたところ、
・N=30…Time=1.6741032600402832
・N=40…Time=9.471542835235596（分子の数1221759）
・N=50…Time=30.463029861450195（分子の数3478761）
と、N=50ぐらいまではスムーズに動くことが確認できた。

解法3:数学的に考えてみたら、すごい簡単な式になった

ここからは完全にオリジナルの内容である。
解法1,2をいきなり見てもプログラミング初心者である私には理解ができなかったので、まずは数学的にアプローチを試みた。
以下の画像が、私の試行錯誤した結果だ。

結論、n回目まで全てアンラッキーだったとき、n+1回目にラッキーが起こる確率P(B|A)は、なんと
$\displaystyle \frac{5n}{6(n+1)}$
で表せるのだ。
なので、コードは以下のようになる。

# n回目までアンラッキーだったとき、(n+1)回目でラッキーが出る確率
n = 10

fi = 5 * n
fb = 6 * (n + 1)

print(fi / fb)
print(fi,'/',fb)
from fractions import Fraction
print(Fraction(fi, fb))

解説

画像に記載の通り、まず私はn=2,n=3,n=4の場合をそれぞれ考えた。
まず条件の方の確率である「n回目まですべてアンラッキーである確率」について解説すると、たとえばn=3のとき、最初に2が出たら、
・(2,2,2)
・(2,2,1)
・(2,1,1)
の3通りとなる。

ここで気づいてほしいのが、最初の2を除いた3組(2,2),(2,1),(1,1)が、ちょうどn=2のときに最初に1が出た場合と2が出た場合の和になっているということだ。
n=3のときに最初に3が出た場合も、n=2のときに最初に1が出た場合+2が出た場合+3が出た場合と等しくなる。
よって、和の記号であるΣを重ねて、分子を表現することができるのだ。

またここで、三重和の公式を使う。
$\displaystyle \sum_{j=1}^{m} \sum_{i=1}^{j} \sum_{k=1}^{i} k = \frac{1}{4!}m(m+1)(m+2)(m+3)$
この公式の証明は、以下の参考文献を見てほしい。
３重 Σ を組合せ論的に求めるhttp://blog.livedoor.jp/ddrerizayoi/archives/56202181.html

すると、今回はサイコロの問題なので、m=6であり、n回アンラッキーの確率P(A)は、
$\displaystyle P(A)=\frac{6\times \frac{1}{n!}・6・(6+1)・…(6+n-1)}{6^{n+1}}$
と求めることができた。

次にn+1回目がラッキーになる同時確率P(A∧B)だが、これにも法則性があり、分子がちょうど
5×(n回目まですべてアンラッキーの分子の最後の数)+4×(n回目まですべてアンラッキーの分子の最後から2番目の数)+…
というふうになっている。これをシグマを使って表し、シグマの公式を当てはめて求めると、
$\displaystyle P(A \land B)=\frac{\sum_{k=1}^{5} \frac{k}{(n-1)!}・(k+1)・(k+2)・…(k+n-1)}{6^{n+1}} =\frac{ \frac{5}{(n-1)!・(n+1)}・(5+1)・(5+2)・…(5+n)}{6^{n+1}}$
と求めることができる。

さて、ここまで計算すれば、あとは約分をしていけばOK。
実際、$(5+1)=6$ , $(5+n)=(6+n-1)$なので、大半がスッキリと約分ができ、結論で示したように
$P(B|A)=\displaystyle \frac{5n}{6(n+1)}$
と、ものすごくきれいな数式に帰着するのだ。

総括

解法1は非常にシンプルな考え方で計算を試みているが、その分計算量が大きくなりすぎてn=10だとかなり処理に負荷がかかる。
一方で解法3(私の解法)は、数学的アプローチが完了したため計算は全くいらないが、実装するまでかなりの時間を要した。
最終的にバランスが取れているのは、解法2の方法であると思う。
なので解法2のようなコードが書けるように、日々努力を積み重ねていきたいと思う。

にしても、数学的にこんなにもきれいな数式($P(B|A)=\frac{5n}{6(n+1)}$)で表せるなんて、取り掛かる前は全く思わなかった。非常に良問である。
もしかしたら確率漸化式等を使えば、もっと簡単に解けるのかな？

以上、何かの参考になれば幸いです。