問題
コインが二枚ある
一枚は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%