はじめに

　『ベイズ推論による機械学習入門』(MLSシリーズ)の独学時のノートです。各種のモデルやアルゴリズムについて「数式・プログラム・図」を用いて解説します。
　本の補助として読んでください。

　この記事では、ベルヌーイ分布に対するベイズ推論をPythonでスクラッチ実装します。

【前節の内容】

【他の節の内容】

【この節の内容】

3.2.1 ベルヌーイ分布のベイズ推論の実装

　ベルヌーイモデル(Bernoulli model)に対するベイズ推論(Bayesian inference)を実装して、人工的に生成したデータを用いて、パラメータの学習と未知変数の予測を行う。ベルヌーイモデルでは、尤度関数をベルヌーイ分布(Bernoulli distribution)、事前分布をベータ分布(Beta distribution)とする。この記事では、Pythonを利用して実装する。
　ベルヌーイモデルについては「3.2.1：ベルヌーイモデルの生成モデルの導出【緑ベイズ入門のノート】 - からっぽのしょこ」、ベイズ推論については「3.2.1：ベルヌーイ分布のベイズ推論の導出【緑ベイズ入門のノート】 - からっぽのしょこ」、Rを利用する場合は「【R】3.2.1：ベルヌーイ分布のベイズ推論の実装【緑ベイズ入門のノート】 - からっぽのしょこ」を参照のこと。

　利用するライブラリを読み込む。

# ライブラリを読込
import numpy as np
from scipy.stats import beta
import matplotlib.pyplot as plt

ベイズ推論の実装

　まずは、ベルヌーイ分布に対するベイズ推論における一連の処理を確認する。
　生成モデル(ベルヌーイモデル)を設定して、モデルに従うデータ(トイデータ)を生成する。続いて、生成した観測データを用いて、事後分布の計算(パラメータの推定)を行う。さらに、事後分布のパラメータ(または観測データ)を用いて、予測分布の計算(未観測データの予測)を行う。

生成分布の設定

　データの生成分布(真の分布・ベルヌーイ分布) $p(x_n \mid \mu_{\mathrm{truth}}) = \mathrm{Bern}(x_n \mid \mu_{\mathrm{truth}})$ のパラメータ(真のパラメータ) $\mu_{\mathrm{truth}}$ を設定する。
　ベルヌーイ分布については「【Python】ベルヌーイ分布の作図 - からっぽのしょこ」を参照のこと。

　生成分布のパラメータ $\mu_{\mathrm{truth}}$ を設定する。

# 真のパラメータを指定
mu_truth = 0.25
print(mu_truth)

0.25

　ベルヌーイ分布の成功確率パラメータ(0から1の値) $0 \leq \mu_{\mathrm{truth}} \leq 1$ を指定する。
　 $\mu_{\mathrm{truth}}$ がベルヌーイモデルにおける真のパラメータであり、この値を求めるのがここでの目的(学習)である。

　生成分布の確率変数 $x$ の作図範囲を設定する。

# x軸の範囲を設定
x_min = 0 # (固定)
x_max = 1 # (固定)
print(x_min, x_max)

# x軸の値を作成
x_vec = np.arange(start=x_min, stop=x_max+1, step=1)
print(x_vec)

0 1
[0 1]

　0か1の2値変数 $x \in \{0, 1\}$ がとり得る範囲を設定する。

　生成分布の確率を計算する。

# 生成分布の確率を計算:式(2.16)
model_prob_vec = np.array([1.0-mu_truth, mu_truth])
print(model_prob_vec)

[0.75 0.25]

　 $x$ の値ごとに、ベルヌーイ分布に従う確率 $\mathrm{Bern}(x \mid \mu_{\mathrm{truth}})$ を計算する。
　この例では、失敗確率 $1-\mu$ と成功確率 $\mu$ を配列に格納している。

　生成分布のグラフを作成する。

# 生成分布のラベルを作成
model_param_lbl = f'$\\mu_{{truth}} = {mu_truth:.3g}$'

# 生成分布を作図
fig, ax = plt.subplots(
    figsize=(8, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Bernoulli distribution', fontsize=20)
ax.bar(
    x=x_vec, height=model_prob_vec, 
    color='red', 
    label='true model', zorder=10
) # 真の分布
ax.set_xticks(ticks=x_vec) # x軸目盛
ax.set_xlabel('$x$')
ax.set_ylabel('probability')
ax.set_title(model_param_lbl, loc='left') # パラメータラベル
ax.legend()
ax.grid(zorder=0)
plt.show()

　真の分布(ベルヌーイ分布)を赤色の棒グラフで示す。

　真のパラメータ $\mu_{\mathrm{truth}}$ を求めることは、真の分布 $\mathrm{Bern}(x_n \mid \mu_{\mathrm{truth}})$ を求めることを意味する。

データの生成

　設定した生成分布(ベルヌーイ分布) $p(x_n \mid \mu_{\mathrm{truth}}) = \mathrm{Bern}(x_n \mid \mu_{\mathrm{truth}})$ に従うデータ(観測データ) $\mathbf{X}$ を作成する。
　ベルヌーイモデルのデータ生成については「準備中」を参照のこと。

　生成分布からデータ $\mathbf{X}$ を生成する。

# データ数を指定
N = 100

# 観測データを生成
x_n = np.random.binomial(n=1, p=mu_truth, size=N)
print(x_n[:10])

[0 0 0 1 0 0 0 0 1 0]

　データ数 $N$ を指定して、ベルヌーイ分布に従う乱数 $\mathbf{X} = \{x_1, \cdots, x_N\}$ を生成する。
　二項分布の乱数生成関数 np.random.binomial() の試行回数の引数 n に 1、成功確率の引数 p に $\mu_{\mathrm{truth}}$ 、サンプルサイズの引数 size に $N$ を指定する。

　観測データ $\mathbf{X}$ を集計する。

# 度数を集計
obs_freq_vec = np.array([N-np.sum(x_n), np.sum(x_n)])
print(obs_freq_vec)

# 相対度数を計算
obs_relfreq_vec = obs_freq_vec / N
print(obs_relfreq_vec)

[75 25]
[0.75 0.25]

　 $x$ の値ごとに、観測データ x_n に含まれる要素数をカウントして、度数 $N_x$ と相対度数 $\frac{N_x}{N}$ を求める。
　この例では、失敗回数 $N_0 = N - N_1$ と成功回数 $N_1 = \sum_{n=1}^N x_n$ を配列に格納している。

　観測データ $\mathbf{X}$ のグラフを作成する。

# 観測データのラベルを作成
obs_param_lbl  = f'$N = {N}, '
obs_param_lbl += f'\\mu_{{truth}} = {mu_truth:.3g}$'

# 相対度数軸の範囲を設定
u = 0.1
relfreq_max = np.max(obs_relfreq_vec)
relfreq_max = max(relfreq_max, model_prob_vec.max()) # 生成確率と比較
relfreq_max = np.ceil(relfreq_max /u)*u # u単位で切り上げ

# 観測データを作図
fig, ax = plt.subplots(
    figsize=(8, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Bernoulli distribution', fontsize=20)
ax2 = ax.twinx() # 第2軸の設定用

ax.bar(
    x=x_vec, height=model_prob_vec, 
    facecolor='none', edgecolor='red', linewidth=1.0, linestyle='--', 
    label='true model\n(generation distribution)', zorder=10
) # 生成分布
ax.bar(
    x=x_vec, height=obs_relfreq_vec, 
    color='hotpink', alpha=0.5, 
    label='observation data', zorder=11
) # 観測データ

ax.set_xticks(ticks=x_vec) # x軸目盛
relfreq_vals = ax.get_yticks()  # 相対度数
freq_vals    = relfreq_vals * N # 度数
ax2.set_yticks(ticks=freq_vals, labels=[f'{y:.1f}' for y in freq_vals]) # 度数軸目盛
ax.set_xlabel('$x$')
ax.set_ylabel('probability, relative frequency')
ax2.set_ylabel('frequency')
ax.set_title(obs_param_lbl, loc='left') # パラメータラベル
ax.legend()
ax.grid(zorder=0)
ax.set_ylim(ymin=0.0, ymax=relfreq_max) # (目盛の共通化用)
ax2.set_ylim(ymin=0.0, ymax=relfreq_max*N) # (目盛の共通化用)

plt.show()

　 $x$ の値ごとの生成確率(生成分布)を赤色(白抜き)の棒グラフ(破線)、観測データ $\mathbf{X}$ の度数 $N_x$ (相対度数 $\frac{N_x}{N}$ )を桃色(塗りつぶし)の棒グラフ(実線)で示す。

　データ数 $N$ が十分に大きいと、観測データ $\mathbf{X}$ のヒストグラムの形状が生成分布 $\mathrm{Bern}(x_n \mid \mu_{\mathrm{truth}})$ に近付く。

事前分布の設定

　パラメータ $\mu$ の事前分布(ベータ分布) $p(\mu \mid a, b) = \mathrm{Beta}(\mu \mid a, b)$ のパラメータ(超パラメータ) $a, b$ を設定する。
　ベータ分布については「準備中」を参照のこと。

　事前分布のパラメータ $a, b$ を設定する。

# 事前分布のパラメータを指定
a = 1.0
b = 1.0
print(a)
print(b)

1.0
1.0

　ベータ分布の形状パラメータ(正の値) $a \gt 0, b \gt 0$ を指定する。

　事前分布の確率変数 $\mu$ の作図範囲を設定する。

# μ軸の範囲を設定
mu_min = 0.0
mu_max = 1.0
print(mu_min, mu_max)

# μ軸の値を作成
mu_vec = np.linspace(start=mu_min, stop=mu_max, num=1001)
print(mu_vec[:5])

0.0 1.0
[0.    0.001 0.002 0.003 0.004]

　0から1の値の変数 $\mu \in (0, 1)$ がとり得る範囲を設定する。

　事前分布の確率密度を計算する。

# 事前分布の確率密度を計算:式(2.41)
prior_dens_vec = beta.pdf(x=mu_vec, a=a, b=b)
print(prior_dens_vec[:5])

[1. 1. 1. 1. 1.]

　 $\mu$ の値ごとに、ベータ分布に従う確率密度 $\mathrm{Bern}(\mu \mid a, b)$ を計算する。
　ベータ分布の確率密度関数 sp.stats.beta.pdf() の確率変数の引数 x に $\mu$ 、形状パラメータの引数 a, b に $a, b$ を指定する。

　事前分布のグラフを作成する。

# 事前分布のラベルを作成
prior_param_lbl  = f'$\\mu_{{truth}} = {mu_truth:.3g}, '
prior_param_lbl += f'a = {a:.3g}, b = {b:.3g}$'

# 事前分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Beta distribution', fontsize=20)
ax2 = ax.twiny() # 第2軸の設定用

ax.axvline(
    x=mu_truth, 
    color='red', linewidth=1.0, linestyle='--', 
    label='true parameter', zorder=10
) # 真のパラメータ
ax.plot(
    mu_vec, prior_dens_vec, 
    color='purple', linewidth=1.0, 
    label='posterior distribution', zorder=11
) # 事前分布

ax2.set_xticks(ticks=[mu_truth], labels=['$\\mu_{truth}$']) # パラメータラベル
ax.set_xlabel('$\\mu$')
ax.set_ylabel('density')
ax.set_title(prior_param_lbl, loc='left') # パラメータラベル
ax.legend()
ax.grid(zorder=0)
ax.set_xlim(xmin=mu_min, xmax=mu_max) # (目盛の共通化用)
ax2.set_xlim(xmin=mu_min, xmax=mu_max) # (目盛の共通化用)

plt.show()

　真のパラメータを赤色の垂直線(破線)、事前分布(ベータ分布)を紫色の曲線(実線)で示す。

　真のパラメータ(真の分布のパラメータ) $\mu_{\mathrm{truth}}$ と、パラメータ $\mu$ の事前分布 $\mathrm{Beta}(\mu \mid a, b)$ の位置関係を図で確認する。

　この例では、無情報事前分布として、全ての確率密度が等しい分布を設定している。

事後分布の計算

　観測データ $\mathbf{X}$ から、パラメータ $\mu$ の事後分布(ベータ分布) $p(\mu \mid \mathbf{X}, a, b) = \mathrm{Beta}(\mu \mid \hat{a}, \hat{b})$ のパラメータ(超パラメータ) $\hat{a}, \hat{b}$ を求める(真のパラメータ $\mu_{\mathrm{truth}}$ を分布推定する)。

　事後分布のパラメータ $\hat{a}, \hat{b}$ を計算する。

# 事後分布のパラメータを計算:式(3.15)
a_hat = np.sum(x_n) + a
b_hat = N - np.sum(x_n) + b
print(a_hat)
print(b_hat)

26.0
76.0

　事後分布のパラメータは、次の式で計算できる。

$\displaystyle \begin{aligned} \hat{a} &= \sum_{n=1}^N x_n + a \\ \hat{b} &= N - \sum_{n=1}^N x_n + b \end{aligned} \tag{3.15}$

　事後分布の確率密度を計算する。

# 事前分布の確率密度を計算:式(2.41)
posterior_dens_vec = beta.pdf(x=mu_vec, a=a_hat, b=b_hat)
print(posterior_dens_vec[:5])

[0.00000000e+00 2.27237102e-50 7.07307257e-43 1.65666946e-38
 2.04182158e-35]

　事前分布のときと同様にして、ベータ分布に従う確率密度 $\mathrm{Beta}(\mu \mid \hat{a}, \hat{b})$ を計算する。

　事後分布のグラフを作成する。

# 事後分布のラベルを作成
posterior_param_lbl  = f'$N = {N}, '
posterior_param_lbl += f'\\mu_{{truth}} = {mu_truth:.3g}, '
posterior_param_lbl += f'\\hat{{a}} = {a_hat:.3g}, \\hat{{b}} = {b_hat:.3g}$'

# 事後分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Beta distribution', fontsize=20)
ax2 = ax.twiny() # 第2軸の設定用

ax.axvline(
    x=mu_truth, 
    color='red', linewidth=1.0, linestyle='--', 
    label='true parameter', zorder=10
) # 真のパラメータ
ax.plot(
    mu_vec, posterior_dens_vec, 
    color='purple', linewidth=1.0, 
    label='posterior distribution', zorder=11
) # 事後分布

ax2.set_xticks(ticks=[mu_truth], labels=['$\\mu_{truth}$']) # パラメータラベル
ax.set_xlabel('$\\mu$')
ax.set_ylabel('density')
ax.set_title(posterior_param_lbl, loc='left') # パラメータラベル
ax.legend()
ax.grid(zorder=0)
ax.set_xlim(xmin=mu_min, xmax=mu_max) # (目盛の共通化用)
ax2.set_xlim(xmin=mu_min, xmax=mu_max) # (目盛の共通化用)

plt.show()

　真のパラメータを赤色の垂直線(破線)、事後分布(ベータ分布)を紫色の曲線(実線)で示す。

　データ数 $N$ が十分に大きいと、パラメータ $\mu$ の事後分布 $\mathrm{Beta}(\mu \mid \hat{a}, \hat{b})$ のピークが真のパラメータ $\mu_{\mathrm{truth}}$ に近付く。

予測分布の計算

　観測データ $\mathbf{X}$ から、未観測データ $x_{*}$ の予測分布(ベルヌーイ分布) $p(x_{*} \mid \mathbf{X}, a, b) = \mathrm{Bern}(x_{*} \mid \hat{\mu}_{*})$ のパラメータ $\hat{\mu}_{*}$ を求める。

　予測分布のパラメータ $\hat{\mu}_{*}$ を計算する。

# 事後分布のパラメータにより予測分布のパラメータを計算:式(3.19')
mu_star_hat = a_hat / (a_hat + b_hat)
print(mu_star_hat)

0.2549019607843137

# 事後分布のパラメータにより予測分布のパラメータを計算:式(3.19')
mu_star_hat = (np.sum(x_n) + a) / (N + a + b)
print(mu_star_hat)

0.2549019607843137

　予測分布のパラメータは、事後分布のパラメータ $\hat{a}, \hat{b}$ または観測データ $\mathbf{X}$ を用いて、次の式で計算できる。

$\displaystyle \begin{aligned} \hat{\mu}_{*} &= \frac{\hat{a}}{\hat{a} + \hat{b}} \\ &= \frac{\sum_{n=1}^N x_n + a}{N + a + b} \end{aligned} \tag{3.19'}$

　予測分布の確率を計算する。

# 予測分布の確率を計算:式(2.16)
predict_prob_vec = np.array([1.0-mu_star_hat, mu_star_hat])
print(predict_prob_vec)

[0.74509804 0.25490196]

　生成分布のときと同様にして、ベルヌーイ分布に従う確率 $\mathrm{Bern}(x \mid \hat{\mu}_{*})$ を計算する。

　予測分布のグラフを作成する。

# 予測分布のラベルを作成
predict_param_lbl  = f'$N = {N}, '
predict_param_lbl += f'\\mu_{{truth}} = {mu_truth:.3g}, '
predict_param_lbl += f'\\hat{{\\mu}}_{{*}} = {mu_star_hat:.3g}$'

# 予測分布を作図
fig, ax = plt.subplots(
    figsize=(8, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Bernoulli distribution', fontsize=20)
ax.bar(
    x=x_vec, height=model_prob_vec, 
    facecolor='none', edgecolor='red', linewidth=1.0, linestyle='--', 
    label='true model', zorder=10
) # 真の分布
ax.bar(
    x=x_vec, height=predict_prob_vec, 
    color='purple', alpha=0.5, 
    label='predict distribution', zorder=11
) # 予測分布
ax.set_xticks(ticks=x_vec) # x軸目盛
ax.set_xlabel('$x$')
ax.set_ylabel('probability')
ax.set_title(predict_param_lbl, loc='left') # パラメータラベル
ax.legend()
ax.grid(zorder=0)
plt.show()

　真の分布(ベルヌーイ分布)を赤色(白抜き)の棒グラフ(破線)、予測分布(ベルヌーイ分布)を紫色(塗りつぶし)の棒グラフ(実線)で示す。

　データ数 $N$ が十分に大きいと、未観測データ $x_{*}$ の予測分布 $\mathrm{Bern}(x_{*} \mid \hat{\mu}_{*})$ の形状が真の分布 $\mathrm{Bern}(x_{n} \mid \mu_{\mathrm{truth}})$ に近付く。

　以上で、ポアソンモデルのベイズ推論を実装した。

スポンサードリンク

学習の推移

　次は、ベルヌーイ分布に対するベイズ推論を図で確認する。
　データ数を増やして分布の変化をアニメーションで確認する。
　作図コードについては「Suyama-Bayes/code/bernoulli_model/bayesian_inference/plot_parameter_updates.py at master · anemptyarchive/Suyama-Bayes · GitHub」を参照のこと。

データ数と分布の関係

　データ数 $N$ を変化させたときの事後分布 $p(\mu \mid \mathbf{X}, a, b)$ の変化をアニメーションにする。

　 $n$ 個のデータから求めた( $n$ 回更新した)事後分布(ベータ分布) $\mathrm{Beta}(\mu \mid a^{(n)}, b^{(n)})$ を紫色の曲線(実線)、 $n$ 番目の観測データ $x_n$ に対応する位置 $\mu = x_n$ を桃色の点で示す。

　「ベイズ推論の実装」では、 $N$ 個(複数)のデータ $\mathbf{X}$ を用いて、事後分布のパラメータ $\hat{a}, \hat{b}$ を一括更新した。
　 $n$ 番目(1つ)のデータ $x_n$ を用いて逐次更新する場合、 $n$ 回目の事後分布のパラメータ $a^{(n)}, b^{(n)}$ は、次の式で計算できる。

$\displaystyle \begin{aligned} a^{(n)} &\leftarrow x_n + a^{(n-1)} \\ b^{(n)} &\leftarrow 1 - x_n + b^{(n-1)} \end{aligned}$

　超パラメータの初期値(1回目の更新における事前分布のパラメータ)を $a^{(0)}, b^{(0)}$ として、 $a^{(n-1)}, b^{(n-1)}$ は、 $n-1$ 回目の更新における事後分布のパラメータであり、また $n$ 回目の更新における事前分布のパラメータを表す。

　データ数 $N$ が大きくなるのに応じて、パラメータ $\mu$ の事後分布 $\mathrm{Beta}(\mu \mid \hat{a}, \hat{b})$ のピークが真のパラメータ $\mu_{\mathrm{truth}}$ に近付いていくのを確認できる。

　データ数 $N$ を変化させたときの予測分布 $p(x_{*} \mid \mathbf{X}, a, b)$ の変化をアニメーションにする。

　 $n$ 個のデータから求めた( $n$ 回更新した)予測分布(ベルヌーイ分布) $\mathrm{Bern}(x_{*} \mid \mu_{*}^{(n)})$ を紫色(塗りつぶし)の棒グラフ(実線)、 $n$ 番目の観測データ $x_n$ を桃色の点で示す。

　「ベイズ推論の実装」では、 $N$ 個(複数)のデータ $\mathbf{X}$ を用いて、予測分布のパラメータ $\hat{\mu}_{*}$ を一括更新した。
　 $n$ 番目(1つ)のデータ $x_n$ を用いて逐次更新する場合、 $n$ 回目の予測分布のパラメータ $\mu_{*}^{(n)}$ は、次の式で計算できる。

$\displaystyle \mu_{*}^{(n)} = \frac{a^{(n)}}{a^{(n)} + b^{(n)}}$

　超パラメータの初期値 $a^{(0)}, b^{(0)}$ を用いて求めたパラメータを $\mu_{*}^{(0)}$ で表す。

　データ数 $N$ が大きくなるのに応じて、未観測データ $x_{*}$ の予測分布 $\mathrm{Bern}(x_{*} \mid \hat{\mu}_{*})$ の形状が真の分布 $\mathrm{Bern}(x_{n} \mid \mu_{\mathrm{truth}})$ に近付いていくのを確認できる。

観測データと分布の関係

　データ数 $N$ の変化による観測データと事後分布、予測分布の関係をアニメーションにする。

　真の分布の期待値 $\mathbb{E}[x_n]$ と真のパラメータ $\mu_{\mathrm{truth}}$ の各図(分布)に対応する位置を赤色の垂直線(破線)、観測データの標本平均 $\bar{x}$ を桃色の垂直線(破線)、事後分布の期待値 $\mathbb{E}[\mu]$ と予測分布の期待値 $\mathbb{E}[x_{*}]$ を紫色の垂直線(破線)で示す。

　生成分布(ベルヌーイ分布)、事後分布(ベータ分布)、予測分布(ベルヌーイ分布)の期待値は、それぞれ次の式で計算できる。

$\displaystyle \begin{align} \mathbb{E}_{\mathrm{Bern}(x_n \mid \mu_{\mathrm{truth}})}[x] &= \mu_{\mathrm{truth}} \\ \mathbb{E}_{\mathrm{Beta}(\mu \mid \hat{a}, \hat{b})}[\lambda] &= \frac{\hat{a}}{\hat{a} + \hat{b}} \\ &= \hat{\mu}_{*} \\ \mathbb{E}_{\mathrm{Bern}(x_{*} \mid \hat{\mu}_{*})}[x] &= \hat{\mu}_{*} \end{align}$

　真の分布 $\mathrm{Bern}(x_n \mid \mu_{\mathrm{truth}})$ や真のパラメータ $\mu_{\mathrm{truth}}$ と、観測データ $\mathbf{X}$ 、事後分布 $\mathrm{Beta}(\mu \mid \hat{a}, \hat{b})$ 、予測分布 $\mathrm{Bern}(x_{*} \mid \hat{\mu}_{*})$ 、またそれぞれの統計量の対応関係を確認できる。

　以上で、ベルヌーイモデルのベイズ推論における学習推移を確認した。

　この記事では、ベルヌーイ分布に対するベイズ推論を実装した。次の記事では、カテゴリ分布を扱う。