イカサマコインである確率をベイズの定理で計算する【python】

#問題

コインが二枚ある
一枚は1/2の確率で表が出るコイン
もう一枚はイカサマで70%が表になるコイン
イカサマのコインかどうかを当てられたら勝ち

今追加で50回試行したとき35回表だった

さて、この時イカサマコインといえるか？

#古典統計の場合

イカサマコインではない場合、50回中25回表になるはず。
比率の差の検定を行うと、
カイ二乗値は4.16
p値は0.04
危険率5%で検定すると、「イカサマでないコインを投げたとして35回表が出る」ことは起きても4%くらいだからイカサマでないコインと言い切ることはできない
という解釈になる

「言い切れない」とか曖昧な結果でなく、
イカサマコインであるかどうかというような1,0の確率が欲しい

#ベイズ統計では

P(X|Y) = \frac{P(Y|X)P(X)}{P(Y)}

この場合
二項分布

\begin{align}
 {}_n C_{x} p^x (1-p)^{n-x}
\end{align}

を尤度P(Y|X)
として扱う

事前確率P(X)はイカサマコインである確率、もしくは公正なコインである確率を置く

今回はイカサマコインも公正なコインもどちらも等しく選ばれるものとして
P(X)=0.5
として考えておく。

起こりうる事象は、
イカサマコインが選ばれて35回表がでる
公正コインが選ばれて35回表がでる

これを計算式で表すと

公正

\begin{align}
 {}_{50} C_{35} (0.5)^{35} (1-0.5)^{50-35}=0.0019
\end{align}

イカサマ

\begin{align}
 {}_{50} C_{35} (0.7)^{35} (1-0.7)^{50-35}=0.1223
\end{align}

となる。
それぞれの値に事前確率としてイカサマコインが選ばれる確率をかける
これだけでは確率として考えられないので、両者の値を足して周辺分布とする

公正

\frac{0.1223×0.5}{0.1223×0.5+0.0019×0.5}=0.016

イカサマ

\frac{0.0019×0.5}{0.1223×0.5+0.0019×0.5}=0.984

import math
import numpy as np
import random

def combinations_count(n, r):
    return math.factorial(n) // (math.factorial(n - r) * math.factorial(r))

loaded_prior=0.5
loaded_heads=0.7
loaded_tails=0.3

fair_prior=1-loaded_prior
fair_heads=0.5
fair_tails=0.5

n=50
x=35
n_x=n-x

def coin(n,x,n_x,lp,lh,lt,fp,fh,ft):
    loaded_likelihood=combinations_count(n, x)*(loaded_heads**x)*(loaded_tails**n_x)
    fair_likelihood=combinations_count(n, x)*(fair_heads**x)*(fair_tails**n_x)
    marginal = loaded_likelihood*loaded_prior + fair_likelihood*fair_prior
    load_prob=(loaded_likelihood*loaded_prior)/marginal
    fair_prob=(fair_likelihood*fair_prior)/marginal
    return load_prob, fair_prob

load_prob,fair_prob= coin(n=n,x=x,n_x=n_x,lp=loaded_prior,lh=loaded_heads,lt=loaded_tails,fp=fair_prior,fh=fair_heads,ft=fair_tails)
print("loaded coin "+str(round(load_prob,3))+"%, fair coin "+str(round(fair_prob,3))+"%")

loaded coin 0.984%, fair coin 0.016%

こうして計算すると、イカサマコインである確率は事前には0.5だったが0.98まで上昇した。
ちなみに33回表が出たら90%を超える

try_time = 50

for o in range(try_time):
    n_x=try_time-o
    lo_p,fa_p=coin(n=try_time,x=o,n_x=n_x,
                          lp=loaded_prior,lh=loaded_heads,lt=loaded_tails,
                          fp=fair_prior,fh=fair_heads,ft=fair_tails)
    print("head "+str(o)+"  loaded coin "+str(round(lo_p,3))+"%, fair coin "+str(round(fa_p,3))+"%")

head 0  loaded coin 0.0%, fair coin 1.0%
head 1  loaded coin 0.0%, fair coin 1.0%
head 2  loaded coin 0.0%, fair coin 1.0%
head 3  loaded coin 0.0%, fair coin 1.0%
head 4  loaded coin 0.0%, fair coin 1.0%
head 5  loaded coin 0.0%, fair coin 1.0%
head 6  loaded coin 0.0%, fair coin 1.0%
head 7  loaded coin 0.0%, fair coin 1.0%
head 8  loaded coin 0.0%, fair coin 1.0%
head 9  loaded coin 0.0%, fair coin 1.0%
head 10  loaded coin 0.0%, fair coin 1.0%
head 11  loaded coin 0.0%, fair coin 1.0%
head 12  loaded coin 0.0%, fair coin 1.0%
head 13  loaded coin 0.0%, fair coin 1.0%
head 14  loaded coin 0.0%, fair coin 1.0%
head 15  loaded coin 0.0%, fair coin 1.0%
head 16  loaded coin 0.0%, fair coin 1.0%
head 17  loaded coin 0.0%, fair coin 1.0%
head 18  loaded coin 0.0%, fair coin 1.0%
head 19  loaded coin 0.0%, fair coin 1.0%
head 20  loaded coin 0.0%, fair coin 1.0%
head 21  loaded coin 0.0%, fair coin 1.0%
head 22  loaded coin 0.001%, fair coin 0.999%
head 23  loaded coin 0.002%, fair coin 0.998%
head 24  loaded coin 0.005%, fair coin 0.995%
head 25  loaded coin 0.013%, fair coin 0.987%
head 26  loaded coin 0.029%, fair coin 0.971%
head 27  loaded coin 0.065%, fair coin 0.935%
head 28  loaded coin 0.14%, fair coin 0.86%
head 29  loaded coin 0.275%, fair coin 0.725%
head 30  loaded coin 0.469%, fair coin 0.531%
head 31  loaded coin 0.674%, fair coin 0.326%
head 32  loaded coin 0.828%, fair coin 0.172%
head 33  loaded coin 0.918%, fair coin 0.082%
head 34  loaded coin 0.963%, fair coin 0.037%
head 35  loaded coin 0.984%, fair coin 0.016%
head 36  loaded coin 0.993%, fair coin 0.007%
head 37  loaded coin 0.997%, fair coin 0.003%
head 38  loaded coin 0.999%, fair coin 0.001%
head 39  loaded coin 0.999%, fair coin 0.001%
head 40  loaded coin 1.0%, fair coin 0.0%
head 41  loaded coin 1.0%, fair coin 0.0%
head 42  loaded coin 1.0%, fair coin 0.0%
head 43  loaded coin 1.0%, fair coin 0.0%
head 44  loaded coin 1.0%, fair coin 0.0%
head 45  loaded coin 1.0%, fair coin 0.0%
head 46  loaded coin 1.0%, fair coin 0.0%
head 47  loaded coin 1.0%, fair coin 0.0%
head 48  loaded coin 1.0%, fair coin 0.0%
head 49  loaded coin 1.0%, fair coin 0.0%

#以上