はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、3.3.2項の内容です。「尤度関数を精度が未知の1次元ガウス分布(正規分布)」、「事前分布をガンマ分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」の計算をPythonで実装します。
省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【他の節の内容】
【この節の内容】
・Pythonでやってみよう
人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。
利用するライブラリを読み込みます。
# 3.3.2項で利用するライブラリ import numpy as np import math # 対数ガンマ関数:lgamma() from scipy.stats import norm, gamma, t # 1次元ガウス分布, ガンマ分布, 1次元スチューデントのt分布 import matplotlib.pyplot as plt
この例では、ガンマ分布・1次元スチューデントのt分布の確率密度を求めるのにガンマ関数$\Gamma(\cdot)$の計算を行います。ガンマ関数の計算には、mathライブラリの対数をとったガンマ関数lgamma()を使います。
または、SciPyライブラリのnorm.pdf()・gamma.pdf()・t.pdf()で1次元ガウス分布・ガンマ分布・1次元スチューデントのt分布の確率密度を直接計算することもできます。
・モデルの構築
まずは、モデルを設定します。
尤度(ガウス分布)のパラメータを設定します。
# 真のパラメータを指定 mu = 25 lambda_truth = 0.01
平均パラメータ$\mu$をmuとして値を指定します。これは与えられている値として使います。
精度パラメータ$\lambda$をlambda_truthとして値を指定します。精度は分散の逆数なので、値が大きいほど散らばり具合が小さくなり、逆数の平方根が標準偏差$\sigma = \sqrt{\lambda^{-1}}$になります。これが真のパラメータであり、この値を求めるのがここでの目的です。
尤度の確率密度を計算します。
# 作図用のxの値を設定 x_line = np.linspace( mu - 4 * np.sqrt(1 / lambda_truth), mu + 4 * np.sqrt(1 / lambda_truth), num=1000 ) # 尤度を計算:式(2.64) ln_C_N = - 0.5 * (np.log(2 * np.pi) - np.log(lambda_truth)) # 正規化項(対数) true_model = np.exp(ln_C_N - 0.5 * lambda_truth * (x_line - mu)**2) # 確率密度 # 尤度を計算:SciPy ver #true_model = norm.pdf(x=x_line, loc=mu, scale=np.sqrt(1 / lambda_truth)) # 確率密度
作図用に、ガウス分布に従う変数$x_n$がとり得る値をnp.linspace()で作成してx_lineとします。この例では、平均値を中心に標準偏差の4倍を範囲とします。np.linspace()を使うと指定した要素数で等間隔に切り分けます。np.arange()を使うと切り分ける間隔を指定できます。処理が重い場合は、この値を調整してください。
x_lineの値ごとに確率密度を計算します。1次元ガウス分布の確率密度は、対数をとった定義式
で計算して、最後にnp.exp()をします。scipy.statsのガウス分布の確率密度関数norm.pdf()でも計算できます。
計算結果を確認しましょう。
# 確認 print(x_line[:10]) print(true_model[:10])
[-15. -14.91991992 -14.83983984 -14.75975976 -14.67967968
-14.5995996 -14.51951952 -14.43943944 -14.35935936 -14.27927928]
[1.33830226e-05 1.38182046e-05 1.42666228e-05 1.47286482e-05
1.52046611e-05 1.56950518e-05 1.62002199e-05 1.67205753e-05
1.72565380e-05 1.78085384e-05]
これを使って作図します。
尤度を作図します。
# 尤度を作図 plt.figure(figsize=(12, 9)) plt.plot(x_line, true_model, color='purple') # 尤度 plt.xlabel('x') plt.ylabel('density') plt.suptitle('Gaussian Distribution', fontsize=20) plt.title('$\mu=' + str(mu) + ', \lambda=' + str(lambda_truth) + '$', loc='left') plt.show()
真のパラメータを求めることは、この真の分布を求めることを意味します。
・データの生成
続いて、構築したモデルに従って観測データ$\mathbf{X} = \{x_1, x_2, \cdots, x_N\}$を生成します。
ガウス分布に従う$N$個のデータをランダムに生成します。
# (観測)データ数を指定 N = 50 # ガウス分布に従うデータを生成 x_n = np.random.normal(loc=mu, scale=np.sqrt(1 / lambda_truth), size=N)
生成するデータ数$N$をNとして値を指定します。
ガウス分布に従う乱数は、np.random.normal()で生成できます。平均の引数locにmu、標準偏差の引数scaleにnp.sqrt(1 / lambda_truth)、試行回数の引数sizeにNを指定します。生成したN個のデータをx_nとします。
観測したデータを確認しましょう。
# 観測データを確認 print(x_n[:5])
[15.1945144 26.17510772 33.79452528 26.72806556 32.47875526]
観測データをヒストグラムでも確認します。
# 観測データのヒストグラムを作図 plt.figure(figsize=(12, 9)) plt.hist(x=x_n, bins=50) # 観測データ plt.xlabel('x') plt.ylabel('count') plt.suptitle('Gaussian Distribution', fontsize=20) plt.title('$N=' + str(N) + ', \mu=' + str(mu) + ', \sigma=' + str(np.sqrt(1 / lambda_truth)) + '$', loc='left') plt.show()
データ数が十分に大きいと、分布の形状が真の分布に近づきます。
・事前分布の設定
次に、尤度に対する共役事前分布を設定します。
$\lambda$の事前分布(ガンマ分布)のパラメータ(超パラメータ)を設定します。
# lambdaの事前分布のパラメータを指定 a = 1 b = 1
ガンマ分布のパラメータ$a,\ b$をそれぞれa, bとして正の実数値を指定します。
$\lambda$の事前分布の確率密度を計算します。
# 作図用のlambdaの値を設定 lambda_line = np.linspace(0, 4 * lambda_truth, num=1000) # lambdaの事前分布を計算:式(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) # 確率密度 # lambdaの事前分布を計算:SciPy ver #prior = gamma.pdf(x=lambda_line, a=a, scale=1 / b) # 確率密度
作図用に、ガンマ分布に従う変数$\lambda$がとり得る値をnp.linspace()で作成してlambda_lineとします。この例では、0から真の値の4倍までを範囲とします。
lambda_lineの値ごとに確率密度を計算します。ガンマ分布の確率密度は、対数をとった定義式
で計算して、最後にnp.exp()をします。ここで、$\Gamma(\cdot)$はガンマ関数です。
ガンマ関数の計算はmath.gamma()で行えますが、値が大きくなると発散してしまします。そこで、対数をとったガンマ関数math.lgamma()で計算した後にnp.exp()で戻します。
ガンマ分布の確率密度は、scipy.statsのgamma.pdf()でも計算できます。
計算結果は次のようになります。
# 確認 print(lambda_line[:10]) print(prior[:10])
[0.00000000e+00 4.00400400e-05 8.00800801e-05 1.20120120e-04
1.60160160e-04 2.00200200e-04 2.40240240e-04 2.80280280e-04
3.20320320e-04 3.60360360e-04]
[ nan 0.99995996 0.99991992 0.99987989 0.99983985 0.99979982
0.99975979 0.99971976 0.99967973 0.9996397 ]
np.log(0)が-Infになります。
$\lambda$の事前分布を作図します。
# lambdaの事前分布を作図 plt.figure(figsize=(12, 9)) plt.plot(lambda_line, prior, label='prior', color='purple') # lambdaの事前分布 plt.xlabel('$\lambda$') plt.ylabel('density') plt.suptitle('Gamma Distribution', fontsize=20) plt.title('$a=' + str(a) + ', b=' + str(b) + '$', loc='left') plt.show()
(lambda_lineの範囲をもっと広くして)a, bの値を変更することで、ガウス分布におけるパラメータと形状の関係を確認できます。
・事後分布の計算
観測データ$\mathbf{X}$からパラメータ$\lambda$の事後分布を求めます(パラメータ$\lambda$を分布推定します)。
観測データx_nを用いて、$\lambda$の事後分布(ガンマ分布)のパラメータを計算します。
# lambdaの事後分布のパラメータを計算:式(3.69) a_hat = 0.5 * N + a b_hat = 0.5 * np.sum((x_n - mu)**2) + b
$\lambda$の事後分布のパラメータは
で計算して、結果をa_hat, b_hatとします。
# 確認 print(a_hat) print(b_hat)
26.0
2366.5072354060935
求めたパラメータを使って、$\lambda$の事後分布の確率密度を計算します。
# lambdaの事後分布の計算:式(2.56) 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) # 確率密度 # lambdaの事前分布を計算:SciPy ver #posterior = gamma.pdf(x=lambda_line, a=a_hat, scale=1 / b_hat) # 確率密度
更新した超パラメータa_hat, b_hatを用いて、事前分布のときと同様にして計算します。
計算結果は次のようになります。
# 確認 print(posterior[:10])
[0.00000000e+00 3.60885189e-48 1.10145652e-40 2.52986438e-36
3.05783511e-33 7.36228787e-31 6.38840371e-29 2.74100925e-27
7.02356462e-26 1.21400102e-24]
$\lambda$の事後分布を作図します。
# lambdaの事後分布を作図 plt.figure(figsize=(12, 9)) plt.plot(lambda_line, posterior, label='posterior', color='purple') # lambdaの事後分布 plt.vlines(x=lambda_truth, ymin=0, ymax=max(posterior), label='$\lambda_{truth}$', color='red', linestyle='--') # 真のlambda plt.xlabel('$\lambda$') plt.ylabel('density') plt.suptitle('Gamma Distribution', fontsize=20) plt.title('$N=' + str(N) + ', a=' + str(a_hat) + ', b=' + str(np.round(b_hat, 1)) + '$', loc='left') plt.legend() plt.show()
パラメータ$\lambda$の真の値付近をピークとする分布を推定できています。
・予測分布の計算
最後に、観測データ$\mathbf{X}$から未観測のデータ$x_{*}$の予測分布を求めます。
$\lambda$の事後分布のパラメータa_hat, b_hat、または観測データx_nと$\lambda$の事前分布のパラメータa, bを用いて、予測分布(スチューデントのt分布)のパラメータを計算します。
# 予測分布のパラメータを計算:式(3.79') mu_s = mu lambda_s_hat = a_hat / b_hat nu_s_hat = 2 * a_hat # 予測分布のパラメータを計算:式(3.79') #lambda_s_hat = (N + 2 * a) / (np.sum((x_n - mu)**2) + 2 * b) #nu_s_hat = N + 2 * a
予測分布のパラメータは
で計算して、結果をmu_s, lambda_s_hat, nu_s_hatとします。$\mu_s$は平均、$\lambda_s$は何？、$\nu_s$は自由度です。
それぞれ上の式だと、事後分布のパラメータa_hat, b_hatを使って計算できます。下の式だと、観測データx_nと事前分布のパラメータa, bを使って計算できます。
# 確認 print(mu_s) print(lambda_s_hat) print(nu_s_hat)
25
0.010986655612543856
52.0
$\mathbf{X}$から$\hat{\lambda}_s,\ \hat{\nu}_s$を学習しているのが式からも分かります。
求めたパラメータを使って、予測分布を計算します。
# 予測分布を計算:式(3.76) ln_C_St = math.lgamma(0.5 * (nu_s_hat + 1)) - math.lgamma(0.5 * nu_s_hat) # 正規化項(対数) ln_term1 = 0.5 * np.log(lambda_s_hat / np.pi / nu_s_hat) ln_term2 = - 0.5 * (nu_s_hat + 1) * np.log(1 + lambda_s_hat / nu_s_hat * (x_line - mu_s)**2) predict = np.exp(ln_C_St + ln_term1 + ln_term2) # 確率密度 # 予測分布を計算:SciPy ver #predict = t.pdf(x=x_line, df=nu_s_hat, loc=mu_s, scale=np.sqrt(1 / lambda_s_hat))
x_lineの値ごとに確率密度を計算します。1次元スチューデントのt分布の確率密度は、対数をとった定義式
で計算して、最後にnp.exp()をします。この$\pi$は円周率です。
ガンマ分布の確率密度は、scipy.statsのt.pdf()でも計算できます。
計算結果は次のようになります。
# 確認 print(predict[:10])
[1.85237888e-05 1.90268207e-05 1.95429919e-05 2.00726288e-05
2.06160654e-05 2.11736434e-05 2.17457124e-05 2.23326300e-05
2.29347623e-05 2.35524834e-05]
予測分布を尤度と重ねて作図します。
# 予測分布を作図 plt.figure(figsize=(12, 9)) plt.plot(x_line, predict, label='predict', color='purple') # 予測分布 plt.plot(x_line, true_model, label='true', color='red', linestyle='--') # 真の分布 plt.xlabel('x') plt.ylabel('density') plt.suptitle("Student's t Distribution", fontsize=20) plt.title('$N=' + str(N) + ', \mu_s=' + str(mu_s) + ', \hat{\lambda}_s=' + str(np.round(lambda_s_hat, 3)) + ', \hat{\\nu}_s=' + str(nu_s_hat) + '$', loc='left') plt.legend() plt.show()
観測データが増えると、予測分布が真の分布に近づきます。
・おまけ：アニメーションで推移の確認
animationモジュールを利用して、事後分布と予測分布の推移をアニメーション(gif画像)で確認するためのコードです。
・コード(クリックで展開)
異なる点のみを簡単に解説します。
# 利用するライブラリ import numpy as np from scipy.stats import norm, gamma, t # 1次元ガウス分布, ガンマ分布, 1次元スチューデントのt分布 import matplotlib.pyplot as plt import matplotlib.animation as animation
・モデルの設定
# 真のパラメータを指定 mu = 25 lambda_truth = 0.01 # lambdaの事前分布のパラメータを指定 a = 1 b = 1 # 初期値による予測分布のパラメータを計算:式(3.79) mu_s = mu lambda_s = a / b nu_s = 2 * a # 作図用のxの値を設定 x_line = np.linspace( mu - 4 * np.sqrt(1 / lambda_truth), mu + 4 * np.sqrt(1 / lambda_truth), num=1000 ) # 作図用のlambdaの値を設定 lambda_line = np.linspace(0, 4 * lambda_truth, num=1000)
・推論処理
各試行の結果をリストに格納していく必要があります。事後分布をtrace_posterior、予測分布をtrace_predict、パラメータについてもそれぞれtrace_***として、初期値の結果を持つように作成しておきます。
# データ数(試行回数)を指定 N = 100 # 推移の記録用の受け皿を初期化 x_n = np.empty(N) trace_a = [a] trace_b = [b] trace_posterior = [gamma.pdf(x=lambda_line, a=a, scale=1 / b)] trace_mu_s = [mu_s] trace_lambda_s = [lambda_s] trace_nu_s = [nu_s] trace_predict = [t.pdf(x=x_line, df=nu_s, loc=mu_s, scale=np.sqrt(1 / lambda_s))] # ベイズ推論 for n in range(N): # ガウス分布に従うデータを生成 x_n[n] = np.random.normal(loc=mu, scale=np.sqrt(1 / lambda_truth), size=1) # lambdaの事前分布のパラメータを更新:式(3.69) a += 0.5 b += 0.5 * (x_n[n] - mu)**2 # lambdaの事前分布(ガンマ分布)を計算:式(2.56) trace_posterior.append(gamma.pdf(x=lambda_line, a=a, scale=1 / b)) # 予測分布のパラメータを更新:式(3.79) mu_s = mu lambda_s = a / b nu_s = 2 * a # 予測分布(スチューデントのt分布)を計算:式(3.76) trace_predict.append( t.pdf(x=x_line, df=nu_s, loc=mu_s, scale=np.sqrt(1 / lambda_s)) ) # 超パラメータを記録 trace_a.append(a) trace_b.append(b) trace_mu_s.append(mu_s) trace_lambda_s.append(lambda_s) trace_nu_s.append(nu_s)
観測された各データによってどのように学習する(分布が変化する)のかを確認するため、for文で1データずつ処理します。よって、データ数Nがイタレーション数になります。
超パラメータに関して、$\hat{a},\ \hat{b}$に対応するa_hat, b_hatを新たに作るのではなく、a, bをイタレーションごとに更新(上書き)していきます。
それに伴い、事後分布のパラメータの計算式(3.69)の$\sum_{n=1}^N (x_n - \mu)^2$と$\frac{N}{2}$の計算は、ループ処理によってN回繰り返し(x_n[n] - mu)**2と0.5を加えることで行います。n回目のループ処理のときには、n-1回分の(x_n[n] - mu)**2と0.5が既にaとbに加えられているわけです。
・事後分布の推移
# 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 作図処理を関数として定義 def update_posterior(n): # 前フレームのグラフを初期化 plt.cla() # nフレーム目のlambdaの事後分布を作図 plt.plot(lambda_line, trace_posterior[n], label='posterior', color='purple') # lambdaの事後分布 plt.vlines(x=lambda_truth, ymin=0, ymax=np.nanmax(trace_posterior), label='$\lambda_{trurh}$', color='red', linestyle='--') # 真のlambda 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(np.round(trace_b[n], 1)) + '$', loc='left') plt.legend() # gif画像を作成 posterior_anime = animation.FuncAnimation(fig, update_posterior, frames=N + 1, interval=100) posterior_anime.save("ch3_3_2_Posterior.gif")
初期値(事前分布)を含むため、フレーム数の引数framesはN + 1です。
・予測分布の推移
# 尤度を計算:式(2.64) true_model = norm.pdf(x=x_line, loc=mu, scale=np.sqrt(1 / lambda_truth)) # 確率密度 # 画像サイズを指定 fig = plt.figure(figsize=(12, 9)) # 作図処理を関数として定義 def update_predict(n): # 前フレームのグラフを初期化 plt.cla() # nフレーム目の予測分布を作図 plt.plot(x_line, trace_predict[n], label='predict', color='purple') # 予測分布 plt.plot(x_line, true_model, label='true', color='red', linestyle='--') # 真の分布 plt.xlabel('x') plt.ylabel('density') plt.suptitle("Student's t Distribution", fontsize=20) plt.title('$N=' + str(n) + ', \hat{\mu}_s=' + str(trace_mu_s[n]) + ', \hat{\lambda}_s=' + str(np.round(trace_lambda_s[n], 3)) + ', \hat{\\nu}_s=' + str(trace_nu_s[n]) + '$', loc='left') plt.ylim(0, np.nanmax(trace_predict)) plt.legend() # gif画像を作成 predict_anime = animation.FuncAnimation(fig, update_predict, frames=N + 1, interval=100) predict_anime.save("ch3_3_2_Predict.gif")
(よく理解していないので、animationの解説は省略…とりあえずこれで動きます……)
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
Python歴もそろそろ1年になります。こうやってスクラッチで組むことしかしてないので、未だにNumPyとPyplotしかほとんど触ったことがない。
【次節の内容】
