ココイチのスプーンの当選確率をベイズ推定する
はじめに
最近ココイチでスプーンがあたるキャンペーンをやっていますね。 ベイズをやっていると、こういうときに手持ちの情報から当選確率をどこまで推定できるか気になりますね。 ということで、5回ココイチに行って2つのスプーンをゲットしてきました。 今回はこの情報を用いて当選確率のベイズ推定を行っていきたいと思います。
当選したスプーン。我が家には炊飯器がないため使用機会はない。
当選確率のベイズ推定
目的変数は当たる、当たらないの二値判定なので、よくあるコイン投げの例題と同じように計算ができます。
下記のサイトを参考にしてベイズ推定の計算を実施し、事後確率の算出を行いました。
http://arduinopid.web.fc2.com/P19.html
数式の説明は上記の記事や、探せばいくらでも見つかるのでここでは割愛します。
３回目くらいの更新まで手計算でやりましたが、 途中でめんどくさくなったので後はpythonに任せます。
#スプーンが当たるか否か data=[0,0,0,1,1] #事前分布(一様) pd=1 #積分するシンボル theta = symbols('θ') X=np.linspace(0,1,100) #事後分布の逐次計算 for num,hit in enumerate(data): if hit==1: f=integrate(pd*theta,(theta,0,1)) k=1/f pd=k*pd*theta else: f=integrate(pd*(1-theta),(theta,0,1)) k=1/f pd=k*pd*(1-theta) plt.plot(X,[pd.subs(theta, i) for i in X]) plt.show()
１回ごとの当選確率の事後分布の変化を見ていきます。
１回目 はずれ
２回目 はずれ
0の確率密度が伸びました。
３回目 はずれ
0の確率密度がさらに伸びました。
４回目 当たり
ようやく当たりが出たため、0の確率密度が0になり、 それらしい分布になってきました。
５回目 当たり
5回中２回中当たりが出ているので、確率密度が最大になるのは0.4なのですが、 だいぶ分布の裾野が広く、私のもつ情報だけだと本当の確率が0.2とか0.6とかでも全然不思議ではないですね。
データを増やす
私の持っているデータだけだと流石に少ないので、もう少し情報を増やします。
検索したら135回挑戦し、23回当選、112回外れという情報の記事を見つけたので、 その数字を使ってもう少し推定を進めたいと思います。
上記の条件でベイズ推定を行った時の、10回ごとの事後分布の変化を以下に示します。
回数を重ねるごとに確率分布の幅が狭まっていく様子がわかります。 135サンプルでベイズ推定を行うことで、本当の確率はかなりの割合で10~25%の間に収まることがわかります。
ちなみに135サンプル実行後の事後分布の式は以下のような結構つらい感じになります。
MCMC
よく使うMCMCでも同じような値になるかを確認してみます。
#MCMC import pystan from pystan import StanModel # Stanコード model = """ data{ int<lower=0> N; int<lower=0, upper=1> y[N]; } parameters{ real<lower=0, upper=1> theta; } model{ for(n in 1:N){ y[n] ~ bernoulli(theta); } } """ stan_data = {'N': len(data),'y': data} sm = pystan.StanModel(model_code=model) %time fit = sm.sampling(data=stan_data, iter=1000, chains=3)
MCMC実行後の事後分布を確認します。
# サンプル列を抽出 la = fit.extract(permuted=True) pd.Series(la["theta"]).plot(kind="kde", secondary_y=True) #カーネル密度推定
逐次計算をしたときとほぼ同じ分布を得ることができました。
まとめ
ベイズ推定を用いてココイチのスプーンの当選確率の推定を行いました。 135サンプルの情報を用いてベイズ推定を行うことで、本来の確率が10~25%にいる可能性が高いことが推定できました。 ベイズはこのように情報が少ない場合でも推定を行うことができ、 値の確からしさを知ることができるのが良いですね。
