書記が数学やるだけ#587 ベータ二項モデルのベイズ統計

鈴華書記（Writer Rinka）

2022年11月21日 06:18

実際に ベータ二項モデルで実践してみる。

問題

1回も当たりが出なかったからといって，当たりくじはないと断言できるのか？

説明

解答

標本データしか考慮しない最尤推定量では当然0になるが，事前分布を考慮することで〜5%くらいであると予測できる。

実装例を以下に示す。今回はPyMC3をメインに用いる（注意：可視化にはArviZを用いることがある）。

import pymc3 as pm
import arviz as az
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
plt.style.use('seaborn-darkgrid')

ベータ分布のグラフを示す。事前分布に比べ，事後分布では下方修正されていることもわかる。

#ベータ分布（青：事前分布，赤：事後分布）
x = np.linspace(0, 1, 1000)
f, ax = plt.subplots(sharex=True, sharey=True)
a = 2
b1 = 20
b2 = 40
y1 = stats.beta(a, b1).pdf(x)
y2 = stats.beta(a, b2).pdf(x)

ax.plot(x, y1)
ax.set_xlim([0,0.25])
ax.legend(fontsize=12)

ax.plot(x, y2)
ax.set_xlim([0,0.25])
ax.legend(fontsize=12)

plt.figure()

実際にモデルを構築していく。

#くじ引き，20回引いて当たりは0回
data = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]

#ベータ二項モデル
with pm.Model() as model:
   #一様事前分布
   theta = pm.Beta('theta', alpha=2, beta=20)
   #標本
   y = pm.Bernoulli('y', p=theta, observed=data)
   #サンプリング
   trace = pm.sample(1000)

結果の要約を示す。