はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、3.2.1項の内容です。尤度関数をベルヌーイ分布、事前分布をベータ分布とした場合のパラメータの事後分布と未観測値の予測分布の計算をPythonで実装します。
省略してある内容等ありますので、本とあせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【他の節の内容】
【この節の内容】
3.2.1 ベルヌーイ分布の学習と予測
人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。
利用するパッケージを読み込みます。
# 3.2.1項で利用するライブラリ import numpy as np import math # 対数ガンマ関数:lgamma() #from scipy.stats import beta # ベータ分布 import matplotlib.pyplot as plt
この例では、ベータ分布の確率密度を求めるのにガンマ関数$\Gamma(\cdot)$の計算を行います。mathライブラリの対数をとったガンマ関数lgamma()を使います。または、SciPyライブラリのbetaを使って直接ベータ分布の確率密度を計算することもできます。
・モデルの構築
まずは、モデルの設定を行います。
尤度(ベルヌーイ分布)$p(x_n | \mu)$のパラメータ$\mu$を設定します。
# 真のパラメータを指定 mu_truth = 0.25
$x_n = 1$となる確率$\mu$をmu_truthとして、$0 \leq \mu \leq 1$の値を指定します。これが真のパラメータであり、この値を求めるのがここでの目的です。
尤度の確率を計算します。
# x軸の値を設定 x_point = np.array([0, 1]) # 尤度(ベルヌーイ分布)を計算 true_model = np.array([1 - mu_truth, mu_truth])
$x_n$がとり得る値0と1をNumpy配列x_pointとして作成しておきます。
それに対応した確率も作成します。ここでは簡単に、$x_n = 0$となる確率$1 - \mu$と合わせて、NumPy配列に格納します。
ベルヌーイ分布の定義式
でも計算できます。ただし$a^0 = 1$なので、$x = 0$のとき$\mu^0 (1 - \mu)^1 = 1 - \mu$、$x = 1$のとき$\mu^1 (1 - \mu)^0 = \mu$になります。よって簡単に$1 - \mu$だけを計算して、true_modelとしています。
結果を確認しましょう。
# 確認 print(true_model)
[0.75 0.25]
これを使って作図します。
尤度を描画します。
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 尤度を作図 plt.bar(x=x_point, height=true_model, color='pruple') # 尤度 plt.xlabel('x') plt.ylabel('prob') plt.xticks(ticks=x_point, labels=x_point) # x軸目盛 plt.suptitle('Bernoulli Distribution', fontsize=20) plt.title('$\mu=' + str(mu_truth) + '$', loc='left') plt.ylim(0.0, 1.0) plt.show()
真のパラメータを求めることは、この真の分布を求めることを意味します。
・データの生成
続いて、構築した尤度に従って観測データ$\mathbf{X} = \{x_1, x_2, \cdots, x_N\}$を生成します。
ベルヌーイ分布に従う$N$個のデータをランダムに生成します。
# データ数を指定 N = 50 # (観測)データを生成 x_n = np.random.binomial(n=1, p=mu_truth, size=N)
生成するデータ数$N$をNとして、値を指定します。
ベルヌーイ分布に従う乱数は、二項分布に従う乱数生成関数np.random.binomial()のn引数を1にすることで生成できます。また、確率の引数pにmu_truth、試行回数の引数sizeにNを指定します。生成したN個のデータをx_nとします。
観測したデータ$\mathbf{X}$を確認しましょう。
# 確認 print(N - np.sum(x_n)) print(np.sum(x_n))
40
10
$x_n = 0$の数は$N - \sum_{n=1}^N x_n$、$x_n = 1$の数は$\sum_{n=1}^N x_n$です。
$\mathbf{X}$をヒストグラムでも確認します。
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 観測データのヒストグラムを作図 plt.bar(x=x_point, height=[N - np.sum(x_n), np.sum(x_n)]) # 観測データ plt.xlabel('x') plt.ylabel('count') plt.xticks(ticks=x_point, labels=x_point) # x軸目盛 plt.suptitle('Observation Data', fontsize=20) plt.title('$N=' + str(N) + ', \mu=' + str(mu_truth) + '$', loc='left') plt.show()
データ数が十分に大きいと、分布の形状が真の分布に近づきます。
・事前分布の設定
尤度に対する共役事前分布を設定します。
事前分布(ベータ分布)のパラメータ(超パラメータ)を設定します。
# 事前分布のパラメータを指定 a = 1.0 b = 1.0
ベータ分布のパラメータ$a,\ b$をそれぞれa, bとして、$a > 0,\ b > 0$の値を指定します。
事前分布の確率密度を計算します。
# x軸の値を設定 mu_line = np.arange(0.0, 1.001, 0.001) # 事前分布(ベータ分布)を計算 ln_C_beta = math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b) # 正規化項 prior = np.exp(ln_C_beta) * mu_line**(a - 1) * (1 - mu_line)**(b - 1) #prior = beta.pdf(x=mu_line, a=a, b=b)
np.arange(0, 1)で、$\mu$がとり得る0から1までの値を用意します。第3引数で間隔を指定できるので、グラフが粗かったり処理が重かったりする場合はこの値を調整してください。
mu_lineの各要素に対して、確率密度を計算します。ベータ分布の確率密度は、定義式
で計算します。ここで、$\Gamma(\cdot)$はガンマ関数です。
ガンマ関数の計算はmath.gamma()で行えますが、値が大きくなると発散してしまします。そこで、対数をとったガンマ関数math.lgamma()で計算した後にnp.exp()を使います。
ベータ分布の確率密度は、beta.pdf()でも計算できます。
計算結果は次のようになります。
# 確認 print(prior)
[1. 1. 1. ... 1. 1. 1.]
事前分布を描画します。
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 事前分布を作図 plt.plot(mu_line, prior, color='purple') plt.xlabel('$\mu$') plt.ylabel('density') plt.suptitle('Beta Distribution', fontsize=20) plt.title('a=' + str(a) + ', b=' + str(b), loc='left') plt.show()
a, bの値を変更することで、ベータ分布におけるパラメータと形状の関係を確認できます。
・事後分布の計算
観測データ$\mathbf{X}$からパラメータ$\mu$の事後分布を求めます(パラメータ$\mu$を分布推定します)。
観測データx_nを用いて事後分布のパラメータを計算します。
# 事後分布のパラメータを計算
a_hat = np.sum(x_n) + a
b_hat = N - np.sum(x_n) + b
事後分布のパラメータは
で計算して、結果をa_hat, b_hatとします。
# 確認 print(a_hat) print(b_hat)
11.0
41.0
事前分布のパラメータ$a,\ b$に、それぞれ$x_n = 1,\ x_n = 0$の数を加えています。
事後分布(ベータ分布)の確率密度を計算します。
# 事後分布(ベータ分布)の確率密度を計算 ln_C_beta = math.lgamma(a_hat + b_hat) - math.lgamma(a_hat) - math.lgamma(b_hat) posterior = np.exp(ln_C_beta) * mu_line**(a_hat - 1) * (1 - mu_line)**(b_hat - 1) #posterior = beta.pdf(x=mu_line, a=a_hat, b=b_hat)
更新した超パラメータa_hat, b_hatを用いて、事前分布のときと同様にして計算します。
計算結果は次のようになります。
# 確認 print(posterior)
[0.00000000e+000 5.03334242e-019 4.95174824e-016 ... 5.64601707e-097
5.18670837e-109 0.00000000e+000]
事後分布を描画します。
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 事後分布を作図 plt.plot(mu_line, posterior, color='purple') # 事後分布 plt.vlines(x=mu_truth, ymin=0, ymax=max(posterior), linestyles='--', color='red') # 真のパラメータ plt.xlabel('$\mu$') plt.ylabel('density') plt.suptitle('Beta Distribution', fontsize=20) plt.title('$N=' + str(N) + ', \hat{a}=' + str(a_hat) + ', \hat{b}=' + str(b_hat) + '$', loc='left') plt.show()
パラメータ$\mu$の真の値付近をピークとする分布を推定できています。
・予測分布の計算
最後に、$\mathbf{X}$から未観測のデータ$x_{*}$の予測分布を求めます。
事後分布のパラメータa_hat, b_hat、または観測データx_nと事前分布のパラメータa, bを用いて予測分布(ベルヌーイ分布)のパラメータを計算します。
# 予測分布のパラメータを計算 mu_hat_star = a_hat / (a_hat + b_hat) #mu_hat_star = (np.sum(x_n) + a) / (N + a + b)
予測分布のパラメータの計算式
の結果をmu_hat_starとします。
上の式だと、事後分布のパラメータa_hat, b_hatを使って計算できます。下の式だと、観測データx_nと事前分布のパラメータa, bを使って計算できます。
# 確認 print(mu_hat_star)
0.21153846153846154
$\hat{\mu}_{*}$は、$x_{*} = 1$となる確率を表し、$\mathbf{X}$から学習しているのが式からも分かります。
予測分布を計算します。
# 予測分布(ベルヌーイ分布)を計算 predict = np.array([1 - mu_hat_star, mu_hat_star])
尤度のときと同様に処理できます。
計算結果は次のようになります。
# 確認 print(predict)
[0.78846154 0.21153846]
予測分布を尤度と重ねて描画します。
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 予測分布を作図 plt.bar(x=x_point, height=true_model, label='true', alpha=0.5, color='white', edgecolor='red', linestyle='dashed') # 真のモデル plt.bar(x=x_point, height=predict, label='predict', alpha=0.5, color='purple') # 予測分布 plt.xlabel('x') plt.ylabel('prob') plt.suptitle('Bernoulli Distribution', fontsize=20) plt.title('$N=' + str(N) + ', \hat{\mu}_{*}=' + str(np.round(mu_hat_star, 2)) + '$', loc='left') plt.ylim(0.0, 1.0) plt.legend() plt.show()
観測データが増えると、予測分布が真の分布に近づきます。(本当は塗りつぶしの色を白にするのではなく、枠をそのままで中だけ透過したい。)
・おまけ：推移の確認
animationモジュールを利用して、パラメータの推定値の推移のアニメーション(gif画像)を作成するためのコードです。
・コード(クリックで展開)
# 利用するライブラリ import numpy as np from scipy.stats import beta # ベータ分布 import matplotlib.pyplot as plt import matplotlib.animation as animation
異なる点のみを簡単に解説します。
# 真のパラメータを指定 mu_truth = 0.4 # 事前分布のパラメータを指定 a = 1.0 b = 1.0 # 初期値による予測分布のパラメーターを計算 mu_star = a / (a + b) # 事後分布のx軸の値を作成 mu_line = np.arange(0.0, 1.001, 0.001)
試行ごとの結果をtrace_***に格納していきます。それぞれ初期値の結果を持つように作成しておきます。
# データ数(試行回数)を指定 N = 100 # 推移の記録用の受け皿を初期化 x_n = np.empty(N) trace_a = [a] trace_b = [b] trace_mu = [mu_star] trace_posterior = [beta.pdf(x=mu_line, a=a, b=b)] trace_predict = [[1 - mu_star, mu_star]] # 推論処理 for n in range(N): # (観測)データを生成 x_n[n] = np.random.binomial(n=1, p=mu_truth, size=1) # 事後分布のパラメータを計算 a += x_n[n] b += 1 - x_n[n] # 事後分布を計算 trace_posterior.append(beta.pdf(x=mu_line, a=a, b=b)) # 予測分布のパラメータを計算 mu_star = a / (a + b) # 予測分布を計算 trace_predict.append([1 - mu_star, mu_star]) # 値を記録 trace_a.append(a) trace_b.append(b) trace_mu.append(mu_star)
観測された各データによってどのように学習する(分布が変化する)のかを確認するため、for文で1データずつ処理します。よって、データ数Nがイタレーション数になります。
パラメータの推定値に関して、$\hat{a},\ \hat{b}$に対応するa_hat, b_hatを新たに作るのではなく、a, bをイタレーションごとに更新していきます。
それに伴い、事後分布のパラメータの計算式(3.15)の$\sum_{n=1}^N$の計算は、forループによってN回繰り返しx_n[n]を加えることで行います。n回目のループ処理のときには、n-1回分のx_n[n](bの場合は1 - x_n[n])が既にaとbに加えられているわけです。
結果は次のようになります。
# 確認 print(trace_a[:10]) print(trace_b[:10]) print(np.round(trace_mu[:10], 2)) print(np.round(trace_posterior[:10], 2)) print(np.round(trace_predict[:10], 2))
[1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0]
[1.0, 2.0, 2.0, 3.0, 4.0, 5.0, 5.0, 6.0, 7.0, 8.0]
[0.5 0.33 0.5 0.4 0.33 0.29 0.38 0.33 0.3 0.27]
[[1. 1. 1. ... 1. 1. 1. ]
[2. 2. 2. ... 0. 0. 0. ]
[0. 0.01 0.01 ... 0.01 0.01 0. ]
...
[0. 0. 0. ... 0. 0. 0. ]
[0. 0. 0. ... 0. 0. 0. ]
[0. 0. 0. ... 0. 0. 0. ]]
[[0.5 0.5 ]
[0.67 0.33]
[0.5 0.5 ]
[0.6 0.4 ]
[0.67 0.33]
[0.71 0.29]
[0.62 0.38]
[0.67 0.33]
[0.7 0.3 ]
[0.73 0.27]]
(こんな感じになります。次の画像の作成時とは別の値です。)
・事後分布の推移
## 事後分布の推移をgif画像化 # 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 作図処理を関数として定義 def update_posterior(n): # 前フレームのグラフを初期化 plt.cla() # nフレーム目の事後分布を作図 plt.plot(mu_line, trace_posterior[n], color='purple') # 事後分布 plt.vlines(x=mu_truth, ymin=0, ymax=np.nanmax(trace_posterior), linestyles='--', color='red') # 真のパラメータ plt.xlabel('$\mu$') plt.ylabel('density') plt.suptitle('Beta Distribution', fontsize=20) plt.title('$N=' + str(n) + ', \hat{a}=' + str(trace_a[n]) + ', \hat{b}=' + str(trace_b[n]) + '$', loc='left') # gif画像を作成 posterior_anime = animation.FuncAnimation(fig, update_posterior, frames=N + 1, interval=100) posterior_anime.save("ch3_2_1_Posterior.gif")
・予測分布の推移
## 予測分布の推移をgif画像化 # x軸の点を作成 x_point = np.array([0, 1]) # 画像サイズを指定 fig = plt.figure(figsize=(12, 8)) # 作図処理を関数として定義 def update_predict(n): # 前フレームのグラフを初期化 plt.cla() # nフレーム目の予測分布を作図 plt.bar(x=x_point, height=[1 - mu_truth, mu_truth], alpha=0.5, color='white', edgecolor='red', linestyle='dashed', label='true') # 真のモデル plt.bar(x=x_point, height=trace_predict[n], alpha=0.5, color='purple', label='predict') # 予測分布 plt.xlabel('x') plt.ylabel('prob') plt.xticks(ticks=x_point, labels=x_point) # x軸目盛 plt.suptitle('Bernoulli Distribution', fontsize=20) plt.title('$N=' + str(n) + ', \hat{\mu}_{*}=' + str(np.round(trace_mu[n], 2)) + '$', loc='left') plt.ylim(0.0, 1.0) plt.legend() # gif画像を作成 predict_anime = animation.FuncAnimation(fig, update_predict, frames=N + 1, interval=100) predict_anime.save("ch3_2_1_Predict.gif")
(よく理解していないので、animationの解説は省略...)
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
Rで実装編の記事を書いてから丁度1年が経ちました。あれからPythonも使えるようになりました！折角なので他のも実装していきます。
2021年2月19日は、モーニング娘。'21の森戸知沙希さんの21歳のお誕生日です！おめでとうございます！
動画はカントリー・ガールズ時代のものです。ちなみにサムネの方ではなくバスケットボール持ってる方です！
【次節の内容】
