コイン投げを例にベイズ推定 試行毎に事後分布をみる

はじめに

最近ベイズ統計が熱いと聞いて、勉強せねばと思い「完全独習 ベイズ統計学入門」を買いました。
途中まではわかりやすいのですが、数式を極力使わないで図をもって説明するスタイルなので、連続分布がでてきた途端にわかりにくく感じました。
そこで勉強のため、簡単な例としてコイン投げを題材に、1回投げるたびに事後分布を見てみようと思います。
MCMC法などはまだよく理解していないので、確率密度関数を細切れにして計算します。

参考

小島寛之「完全独習 ベイズ統計学入門」ダイヤモンド社, 2015年11月

コイン投げを例題とした、上記の書籍と同様の図を用いた説明です。
北数教 第108回数学教育実践研究会 平成31年1月26日(土)「ベイズ推定の思想」

こちらは同じ題材ですが、pyMCを使ってMCMC法で計算した発展的な内容です。
私もpyMC使えるようになりたいです。
コインの裏表の事後分布を試行毎に見てみる(PyMC)

「完全独習 ベイズ統計学入門」の初めの章にあたります。図が多用されていてわかりやすいです。
5分でスッキリ理解するベイズ推定

ベイズ推定

ベイズ統計の強みは、
1. データが少なくても推測でき、データが多くなるほど正確になるという性質
2. 入ってくる情報に瞬時に反応して、自動的に推測をアップデートする学習機能
にあるそうです¹。

ベイズ推定は、ベイズの定理が基礎になっています。
$$p(\theta|x)=p(\theta) \frac{p(x|\theta)}{\int p(x|\theta)p(\theta)d\theta}$$
ここで$\theta$はパラメータ、$p(\theta)$は事前分布、$p(x|\theta)$は$\theta$であるときの$x$の条件付き確率(尤度)、$p(\theta|x)$は事後分布です。
コイン投げの例で言えば、$\theta$は表か裏、$p(\theta)$は$\theta$である確率分布、$x$はコインが表か裏かのデータで、分母は$x$が得られる確率を表します。
元々の$\theta$が起きる確率は事前確率$p(\theta)$としていましたが、$x$が起きるという情報を得た元では事後確率$p(\theta|x)$になります。これをベイズ更新と言います。

情報を得て更新された事後分布$p(\theta|x)$を再び事前分布$p(\theta)$として、次々に更新していくことができます。掛け算の順番を入れ替えても答えは変わらないので、コイン投げの例でいうと、$x=${表, 表, 裏, 裏}でも$x=${表, 裏, 表, 裏}でも最終的な事後分布は同じになります。これをベイズ推定の逐次合理性と言います。今までのデータの結果が事前分布に反映されているので、ある種の学習機能と言えます²。

プログラム

0. モジュール

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Times New Roman"      #全体のフォントを設定
plt.rcParams["xtick.direction"] = "in"               #x軸の目盛線を内向きへ
plt.rcParams["ytick.direction"] = "in"               #y軸の目盛線を内向きへ
plt.rcParams["xtick.minor.visible"] = True           #x軸補助目盛りの追加
plt.rcParams["ytick.minor.visible"] = True           #y軸補助目盛りの追加
plt.rcParams["xtick.major.width"] = 1.5              #x軸主目盛り線の線幅
plt.rcParams["ytick.major.width"] = 1.5              #y軸主目盛り線の線幅
plt.rcParams["xtick.minor.width"] = 1.0              #x軸補助目盛り線の線幅
plt.rcParams["ytick.minor.width"] = 1.0              #y軸補助目盛り線の線幅
plt.rcParams["xtick.major.size"] = 10                #x軸主目盛り線の長さ
plt.rcParams["ytick.major.size"] = 10                #y軸主目盛り線の長さ
plt.rcParams["xtick.minor.size"] = 5                 #x軸補助目盛り線の長さ
plt.rcParams["ytick.minor.size"] = 5                 #y軸補助目盛り線の長さ
plt.rcParams["font.size"] = 14                       #フォントの大きさ
plt.rcParams["axes.linewidth"] = 1.5                 #囲みの太さ

1. 一様な事前分布

前もってなんの知識もないときは、とりあえず一様な事前分布を仮定します。これを理由不十分の原則と言うそうです。
まずはこの原則に従ってやってみます。

N = 1000 #計算分割数

p_x_front = np.linspace(0.0,1.0,N) #p(x|表)
p_x_back = 1.0 - p_x_front #p(x|裏)

prior_dist = np.ones(N) #p(θ)事前分布

fig,axes = plt.subplots(figsize=(8,6))
axes.plot(p_x_front,prior_dist,label="Prior distribution")
axes.scatter(np.sum(prior_dist*p_x_front/N),0.5,s=100,label="expected value")
axes.legend(loc="best")
axes.set_xlabel(r"$p(\theta)$")
axes.set_ylabel("Probability density")
axes.set_ylim(-0.1,2.5)
axes.set_xticks(np.linspace(0.0,1.0,11))

縦軸は確率密度です。積分すると1になります。丸プロットは期待値です。期待値の縦軸の値は見やすいところにくるようにしただけです。

def update(p_x_theta,p_theta):
    return (p_x_theta * p_theta) / np.sum(p_x_theta * p_theta / N)

fig,axes = plt.subplots(figsize=(8,6))
dist_front1_back0 = update(p_x_front,prior_dist)
axes.plot(p_x_front,dist_front1_back0,label="front=1, back=0")
axes.scatter(np.sum(dist_front1_back0*p_x_front/N),0.5,s=100)

dist_front1_back1 = update(p_x_back,dist_front1_back0)
axes.plot(p_x_front,dist_front1_back1,label="front=1, back=1")
axes.scatter(np.sum(dist_front1_back1*p_x_front/N),0.5,s=100)

dist_front2_back1 = update(p_x_front,dist_front1_back1)
axes.plot(p_x_front,dist_front2_back1,label="front=2, back=1")
axes.scatter(np.sum(dist_front2_back1*p_x_front/N),0.5,s=100)

dist_front2_back2 = update(p_x_back,dist_front2_back1)
axes.plot(p_x_front,dist_front2_back2,label="front=2, back=2")
axes.scatter(np.sum(dist_front2_back2*p_x_front/N),0.2,s=100)

axes.legend(loc="best")
axes.set_xlabel(r"$p(\theta)$")
axes.set_ylabel("Probability density")
axes.set_ylim(-0.1,2.5)
axes.set_xticks(np.linspace(0.0,1.0,11))

update関数はベイズの定理そのものです。コインが表か裏かの情報を得るたびに事後分布が変化します。また、データが多くなるほど正確になるという性質がオレンジプロットと赤プロットを比べるとわかります。

表100回、裏100回でたときの事後分布を見てみます。

def updates(front,back):
    p_theta = prior_dist[:]

    for i in range(front):
        p_theta = update(p_x_front,p_theta)

    for i in range(back):
        p_theta = update(p_x_back,p_theta)

    return p_theta

p_theta = updates(100,100)
fig,axes = plt.subplots(figsize=(8,6))
axes.plot(p_x_front,p_theta)
axes.scatter(np.sum(p_theta*p_x_front/N),0.5,s=100)
axes.set_xlabel(r"$p(\theta)$")
axes.set_ylabel("Probability density")
axes.set_ylim(-0.1,12.0)
axes.set_xticks(np.linspace(0.0,1.0,11))

だいぶ分散が小さくなりました。期待値は0.5です。

2. 偏った事前分布

ベイズ推定には主観を反映させることができます。今度は表が出やすいだろうと思っているとして、事前分布に$\alpha=3,\beta=2$のベータ分布を用います。

a = 3.0
b = 2.0
x = np.linspace(0.0,1.0,N)
prior_dist = x**(a-1.0) * (1.0-x)**(b-1.0)
prior_dist /= np.sum(prior_dist) / N

fig,axes = plt.subplots(figsize=(8,6))
axes.plot(p_x_front,prior_dist,label="Prior distribution")
axes.scatter(np.sum(prior_dist*p_x_front/N),0.5,s=100,label="expected value")
axes.legend(loc="best")
axes.set_xlabel(r"$p(\theta)$")
axes.set_ylabel("Probability density")
axes.set_ylim(-0.1,2.5)
axes.set_xticks(np.linspace(0.0,1.0,11))

期待値はおよそ0.6です。

fig,axes = plt.subplots(figsize=(8,6))
dist_front1_back0 = update(p_x_front,prior_dist)
axes.plot(p_x_front,dist_front1_back0,label="front=1, back=0")
axes.scatter(np.sum(dist_front1_back0*p_x_front/N),0.5,s=100)

dist_front1_back1 = update(p_x_back,dist_front1_back0)
axes.plot(p_x_front,dist_front1_back1,label="front=1, back=1")
axes.scatter(np.sum(dist_front1_back1*p_x_front/N),0.5,s=100)

dist_front2_back1 = update(p_x_front,dist_front1_back1)
axes.plot(p_x_front,dist_front2_back1,label="front=2, back=1")
axes.scatter(np.sum(dist_front2_back1*p_x_front/N),0.5,s=100)

dist_front2_back2 = update(p_x_back,dist_front2_back1)
axes.plot(p_x_front,dist_front2_back2,label="front=2, back=2")
axes.scatter(np.sum(dist_front2_back2*p_x_front/N),0.2,s=100)

axes.legend(loc="best")
axes.set_xlabel(r"$p(\theta)$")
axes.set_ylabel("Probability density")
axes.set_ylim(-0.1,2.5)
axes.set_xticks(np.linspace(0.0,1.0,11))

一様な事前分布を用いたときと比べると、表が出やすいだろうという主観に引っ張られて期待値が高くなっています。

再び、表100回、裏100回でたときの事後分布を見てみます。

p_theta = updates(100,100)
fig,axes = plt.subplots(figsize=(8,6))
axes.plot(p_x_front,p_theta)
axes.scatter(np.sum(p_theta*p_x_front/N),0.5,s=100)
axes.set_xlabel(r"$p(\theta)$")
axes.set_ylabel("Probability density")
axes.set_ylim(-0.1,12.0)
axes.set_xticks(np.linspace(0.0,1.0,11))

期待値が0.502となり、僅かに0.5より大きくなりましたが、データの数が増えてくると段々と初めの主観の影響が薄れ、正確な分布を得られるようになりました。

まとめ

ベイズ推定の大まかな流れが掴めました。私はMCMC法を理解していないので、次はRとStanではじめる ベイズ統計モデリングによるデータ分析入門を読んで、今回より発展的な内容について書けたらいいなと思います。

---

1. 小島寛之「完全独習 ベイズ統計学入門」ダイヤモンド社, 2015年11月, p.6 ↩

2. 小島寛之「完全独習 ベイズ統計学入門」ダイヤモンド社, 2015年11月, p.145 ↩