はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、3.2.3項の内容です。尤度関数をポアソン分布、事前分布をガンマ分布とした場合のパラメータの事後分布と未観測値の予測分布の計算をPythonで実装します。
省略してある内容等ありますので、本とあせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【他の節の内容】
【この節の内容】
・Pythonでやってみよう
人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。
利用するパッケージを読み込みます。
# 3.2.3項で利用するライブラリ import numpy as np import math # 対数ガンマ関数:lgamma() from scipy.stats import poisson, gamma, nbinom # ポアソン分布,ガンマ分布,負の二項分布 import matplotlib.pyplot as plt
この例では、ポアソン分布・負の二項分布の確率、ガンマ分布の確率密度を求めるのにガンマ関数$\Gamma(\cdot)$の計算を行います。ガンマ関数の計算には、mathライブラリの対数をとったガンマ関数lgamma()を使います。
または、SciPyライブラリのpoisson.pmf()・nbinom.pmf()でポアソン分布・負の二項分布の確率、gamma.pdf()でガンマ分布の確率密度を直接計算することもできます。
・モデルの構築
まずは、モデルの設定を行います。
尤度(ポアソン分布)$p(\mathbf{X} | \lambda)$のパラメータ$\lambda$を設定します。
# 真のパラメータを指定 lambda_truth = 4.0
$\lambda$をlambda_truthとして指定します。これが真のパラメータであり、この値を求めるのがここでの目的です。
尤度の確率を計算します。
# 作図用のxの値を設定 x_line = np.arange(4 * lambda_truth) # 尤度を計算:式(2.37) true_model = np.exp( x_line * np.log(lambda_truth) - [math.lgamma(x + 1) for x in x_line] - lambda_truth ) # 尤度を計算:SciPy ver true_model = poisson.pmf(k=x_line, mu=lambda_truth)
作図用に、ポアソン分布に従う変数$x_n$がとり得る非負の整数をnp.arange()で作成します。この例では、0からlambda_truthの4倍を範囲とします。
x_lineの各要素の値となる確率は、ポアソン分布の定義式
や、ポアソン分布の確率計算メソッドpoisson.pmf()で計算できます。
計算結果を確認しましょう。
# 確認 print(np.round(true_model, 3))
[0.018 0.073 0.147 0.195 0.195 0.156 0.104 0.06 0.03 0.013 0.005 0.002
0.001 0. 0. 0. ]
これを使って作図します。
尤度を作図します。
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 尤度を作図 plt.bar(x=x_line, height=true_model, color='purple') # 尤度 plt.xlabel('x') plt.ylabel('prob') plt.suptitle('Poisson Distribution', fontsize=20) plt.title('$\lambda=' + str(lambda_truth) + '$', loc='left') plt.show()
真のパラメータを求めることは、この真の分布を求めることを意味します。
・データの生成
続いて、構築した尤度に従って観測データ$\mathbf{X} = \{x_1, x_2, \cdots, x_N\}$を生成します。
ポアソン分布に従う$N$個のデータをランダムに生成します。
# (観測)データ数を指定 N = 50 # ポアソン分布に従うデータをランダムに生成 x_n = np.random.poisson(lam=lambda_truth, size=N)
生成するデータ数$N$をNとして、値を指定します。
ポアソン分布に従う乱数は、np.random.poisson()で生成できます。パラメータの引数lambdaにlambda_truth、試行回数の引数sizeにNを指定します。生成したN個のデータをx_nとします。
観測したデータ$\mathbf{X}$を確認しましょう。
# 観測データを確認 print(x_n[:10])
[3 5 3 3 0 4 1 4 6 1]
$x_n$は、非負の整数になっています。(Rのtable()のように簡単にカウントする関数ってないでしょうか？np.unique(x_n, return_counts=True)でカウントできると教えていただきました！)
$\mathbf{X}$をヒストグラムでも確認します。
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 観測データのヒストグラムを作図 plt.bar(x_line, [np.sum(x_n == x) for x in x_line]) # 観測データ plt.xlabel('x') plt.ylabel('count') plt.suptitle('Observation Data', fontsize=20) plt.title('$N=' + str(N) + ', \lambda=' + str(lambda_truth) + '$', loc='left') plt.show()
データ数が十分に大きいと、分布の形状が真の分布に近づきます。
・事前分布の設定
尤度に対する共役事前分布を設定します。
事前分布(ガンマ分布)のパラメータ(超パラメータ)を設定します。
# 事前分布のパラメータを指定 a = 1 b = 1
ガンマ分布のパラメータ$a,\ b$をそれぞれa, bとして、正の実数値を指定します。
事前分布の確率密度を計算します。
# 作図用のlambdaの値を設定 lambda_line = np.arange(0, 2 * lambda_truth, 0.001) # 事前分布を計算:式(2.56) ln_C_gam = a * np.log(b) - math.lgamma(a) # 正規化項(対数) prior = np.exp(ln_C_gam + (a - 1) * np.log(lambda_line) - b * lambda_line) # 確率密度 # 事前分布を計算:SciPy ver prior = gamma.pdf(x=lambda_line, a=a, scale=1 / b)
$\lambda$がとり得る0以上の値をnp.arange()で用意します。第3引数で間隔を指定できるので、グラフが粗かったり処理が重かったりする場合はこの値を調整してください。この例では、lambda_truthの2倍までとします。
lambda_lineの各要素の値に対して、確率密度を計算します。ガンマ分布の確率密度は、定義式
で計算します。ここで、$\Gamma(\cdot)$はガンマ関数であり、$\Gamma(x + 1) = x!$です。
ガンマ関数の計算はmath.gamma()で行えますが、値が大きくなると発散してしまします。そこで、対数をとったガンマ関数math.lgamma()で計算した後にnp.exp()で戻します。
ガンマ分布の確率密度は、gamma.pdf()でも計算できます。
計算結果は次のようになります。
# 確認 print(np.round(prior, 3))
[1. 0.999 0.998 ... 0. 0. 0. ]
np.log(0)がnanになります。
事前分布を作図します。
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 事前分布を作図 plt.plot(lambda_line, prior, color='purple') # 事前分布 plt.xlabel('$\lambda$') plt.ylabel('density') plt.suptitle('Gamma Distribution', fontsize=20) plt.title('$a=' + str(a) + ', b=' + str(b) + "$", loc='left') plt.show()
a, bの値を変更することで、ポアソン分布におけるパラメータと形状の関係を確認できます。
・事後分布の計算
観測データ$\mathbf{X}$からパラメータ$\lambda$の事後分布を求めます(パラメータ$\lambda$を分布推定します)。
観測データx_nを用いて事後分布のパラメータを計算します。
# 事後分布のパラメータを計算:式(3.38)
a_hat = np.sum(x_n) + a
b_hat = N + b
事後分布のパラメータは
で計算して、結果をa_hat, b_hatとします。
# 確認 print(a_hat) print(b_hat)
205
51
事前分布のパラメータ$a,\ b$に、それぞれ観測データの総和$\sum_{n=1}^N x_n$とデータ数$N$を加えています。
事後分布(ガンマ分布)の確率密度を計算します。
# 事後分布を計算 ln_C_gam = a_hat * np.log(b_hat) - math.lgamma(a_hat) # 正規化項(対数) posterior = np.exp( ln_C_gam + (a_hat - 1) * np.log(lambda_line) - b_hat * lambda_line ) # 確率密度 # 事前分布を計算:SciPy ver posterior = gamma.pdf(x=lambda_line, a=a_hat, scale=1 / b_hat)
更新した超パラメータa_hat, b_hatを用いて、事前分布のときと同様にして計算します。
計算結果は次のようになります。
# 確認 print(np.round(posterior, 5))
[0. 0. 0. ... 0. 0. 0.]
事後分布を作図します。
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 事後分布を作図 plt.plot(lambda_line, posterior, color='purple') # 事後分布 plt.vlines(x=lambda_truth, ymin=0, ymax=max(posterior), color='red', linestyle='--') # 真のパラメータ plt.xlabel('$\lambda$') plt.ylabel('density') plt.suptitle('Gamma Distribution', fontsize=20) plt.title('$N=' + str(N) + ', \hat{a}=' + str(a_hat) + ', \hat{b}=' + str(b_hat) + "$", loc='left') plt.show()
パラメータ$\lambda$の真の値付近をピークとする分布を推定できています。
・予測分布の計算
最後に、$\mathbf{X}$から未観測のデータ$x_{*}$の予測分布を求めます。
事後分布のパラメータa_hat, b_hat、または観測データx_nと事前分布のパラメータa, bを用いて予測分布(負の二項分布分布)のパラメータを計算します。
# 予測分布のパラメータを計算:式(3.44') r_hat = a_hat p_hat = 1 / (b_hat + 1)
予測分布のパラメータは
で計算して、結果をそれぞれr_hat, p_hatとします。
上の式だと、事後分布のパラメータa_hat, b_hatを使って計算できます。下の式だと、観測データx_nと事前分布のパラメータa, bを使って計算できます。
# 確認 print(r_hat) print(p_hat)
205
0.019230769230769232
$\hat{r},\ \hat{p}$は、$\mathbf{X}$から学習しているのが式からも分かります。
予測分布を計算します。
# 予測分布を計算:式(3.43) ln_C_NB = np.array([math.lgamma(x + r_hat) - math.lgamma(x + 1) for x in x_line]) - math.lgamma(r_hat) # 正規化項(対数) predict = np.exp(ln_C_NB + r_hat * np.log(1 - p_hat) + x_line * np.log(p_hat)) # 確率 # 予測分布を計算:SciPy ver predict = nbinom.pmf(k=x_line, n=r_hat, p=1 - p_hat) # 確率密度
x_lineの各様の値に対して、確率を計算します。負の二項分布の確率は、定義式
または負の二項分布の確率計算関数dnbinom()で計算します。
計算結果は次のようになります。
# 確認 print(np.round(predict, 3))
[0.019 0.074 0.146 0.193 0.193 0.156 0.105 0.061 0.031 0.014 0.006 0.002
0.001 0. 0. 0. ]
予測分布を尤度と重ねて作図します。
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 予測分布を作図 plt.bar(x=x_line, height=true_model, label='true', alpha=0.5, color='white', edgecolor='red', linestyle='--') # 真の分布 plt.bar(x=x_line, height=predict, label='predict', alpha=0.5, color='purple') # 予測分布 plt.xlabel('x') plt.ylabel('prod') plt.suptitle('Negative Binomial Distribution', fontsize=20) plt.title('$N=' + str(N) + ', \hat{r}=' + str(r_hat) + ', \hat{p}=' + str(np.round(p_hat, 3)) + '$', loc='left') plt.show()
観測データが増えると、予測分布が真の分布に近づきます。
・おまけ：推移の確認
animationモジュールを利用して、パラメータの推定値の推移のアニメーション(gif画像)を作成するためのコードです。
・コード(クリックで展開)
異なる点のみを簡単に解説します。
# 利用するライブラリ import numpy as np from scipy.stats import poisson, gamma, nbinom # ポアソン分布,ガンマ分布,負の二項分布 import matplotlib.pyplot as plt import matplotlib.animation as animation
・モデルの設定
# 真のパラメータを指定 lambda_truth = 4.0 # 事前分布のパラメータを指定 a = 1 b = 1 # 初期値による予測分布のパラメータを計算:式(3.44) r = a p = 1 / (b + 1)
・推論処理
試行ごとの結果をtrace_***に格納していきます。それぞれ初期値の結果を持つように作成しておきます。
# データ数(試行回数)を指定 N = 100 # 作図用の値を設定 lambda_line = np.arange(0, 2 * lambda_truth, 0.001) x_line = np.arange(4 * int(lambda_truth)) # 推移の記録用の受け皿を初期化 x_n = np.empty(N) trace_a = [a] trace_b = [b] trace_posterior = [gamma.pdf(x=lambda_line, a=a, scale=1 / b)] trace_r = [r] trace_p = [p] trace_predict = [nbinom.pmf(k=x_line, n=r, p=1 - p)] # ベイズ推論 for n in range(N): # ポアソン分布に従うデータを生成 x_n[n] = np.random.poisson(lam=lambda_truth, size=1) # 事後分布のパラメータを計算:式(3.38) a += x_n[n] b += 1 # 事後分布(ガンマ分布)を計算:式(2.56) trace_posterior.append(gamma.pdf(x=lambda_line, a=a, scale=1 / b)) # 確率密度 # 予測分布のパラメータを計算:式(3.44) r = a p = 1 / (b + 1) # 予測分布(負の二項分布)を計算:式(3.43) trace_predict.append(nbinom.pmf(k=x_line, n=r, p=1 - p)) # 確率 # 超パラメータを記録 trace_a.append(a) trace_b.append(b) trace_r.append(r) trace_p.append(p)
観測された各データによってどのように学習する(分布が変化する)のかを確認するため、for文で1データずつ処理します。よって、データ数Nがイタレーション数になります。
超パラメータに関して、$\hat{a},\ \hat{b}$に対応するa_hat, b_hatを新たに作るのではなく、a, bをイタレーションごとに更新(上書き)していきます。
それに伴い、事後分布のパラメータの計算式(3.38)の$\sum_{n=1}^N x_n$と$N$の計算は、ループ処理によってN回繰り返しx_n[n]と1を加えることで行います。n回目のループ処理のときには、n-1回分のx_n[n](bの場合は1)が既にaとbに加えられているわけです。
結果は次のようになります。
# 確認 print(trace_a[:10]) print(trace_b[:10]) print(trace_r[:10]) print(np.round(trace_p[:10], 3)) print(np.round(trace_posterior[:5], 2)) print(np.round(trace_predict[:5], 2))
[1, 5.0, 8.0, 9.0, 10.0, 16.0, 18.0, 19.0, 24.0, 27.0]
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
[1, 5.0, 8.0, 9.0, 10.0, 16.0, 18.0, 19.0, 24.0, 27.0]
[0.5 0.333 0.25 0.2 0.167 0.143 0.125 0.111 0.1 0.091]
[[1. 1. 1. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]]
[[0.5 0.25 0.13 0.06 0.03 0.02 0.01 0. 0. 0. 0. 0. 0. 0.
0. 0. ]
[0.13 0.22 0.22 0.17 0.11 0.07 0.04 0.02 0.01 0. 0. 0. 0. 0.
0. 0. ]
[0.1 0.2 0.23 0.19 0.13 0.08 0.04 0.02 0.01 0. 0. 0. 0. 0.
0. 0. ]
[0.13 0.24 0.24 0.18 0.11 0.06 0.03 0.01 0. 0. 0. 0. 0. 0.
0. 0. ]
[0.16 0.27 0.25 0.16 0.09 0.04 0.02 0.01 0. 0. 0. 0. 0. 0.
0. 0. ]]
(この値は下の画像の作成時とは別のときのものです。)
・事後分布の推移
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 作図処理を関数として定義 def update_posterior(n): # 前フレームのグラフを初期化 plt.cla() # nフレーム目の事後分布を作図 plt.plot(lambda_line, trace_posterior[n], color='purple') # 事後分布 plt.vlines(x=lambda_truth, ymin=0, ymax=np.nanmax(trace_posterior), color='red', linestyles='--') # 真のパラメータ plt.xlabel('$\lambda$') plt.ylabel('density') plt.suptitle('Gamma 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_3_Posterior.gif")
初期値(事前分布)を含むため、フレーム数はN + 1です。
・予測分布の推移
## 予測分布の推移をgif画像化 # 尤度を計算:式(2.37) true_model = poisson.pmf(k=x_line, mu=lambda_truth) # 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 作図処理を関数として定義 def update_predict(n): # 前フレームのグラフを初期化 plt.cla() # nフレーム目の予測分布を作図 plt.bar(x=x_line, height=true_model, label='true', alpha=0.5, color='white', edgecolor='red', linestyle='--') # 真の分布 plt.bar(x=x_line, height=trace_predict[n], label='predict', alpha=0.5, color='purple') # 予測分布 plt.xlabel('x') plt.ylabel('prob') plt.suptitle('Bernoulli Distribution', fontsize=20) plt.suptitle('Negative Binomial Distribution', fontsize=20) plt.title('$N=' + str(n) + ', \hat{r}=' + str(trace_r[n]) + ', \hat{p}=' + str(np.round(trace_p[n], 3)) + '$', loc='left') plt.ylim(0.0, np.nanmax(trace_predict)) # y軸の表示範囲 plt.legend() # gif画像を作成 predict_anime = animation.FuncAnimation(fig, update_predict, frames=N + 1, interval=100) predict_anime.save("ch3_2_3_Predict.gif")
(よく理解していないのでanimationモジュールの解説は省略...)
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
1年越しで気付けた理解しきれてなかったり勘違いしてたりなことが結構ある。
【次節の内容】
