はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「Rで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は3.3.2項の内容です。尤度関数・事前分布を1次元ガウス分布とした場合の精度パラメータの事後分布を導出し、学習した事後分布を用いた予測分布を導出します。またその推論過程をR言語で実装します。
省略してる内容等ありますので、本と併せて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てば幸いです。
【前節の内容】
【他の節一覧】
【この節の内容】
3.3.2 精度が未知の場合
・事後分布の計算
ガウス分布に従うと仮定するN個の観測データ$\boldsymbol{X} = \{x_1, x_2, \cdots, x_N\}$を用いて、精度パラメータ$\lambda$の事後分布を求めていく。
観測モデル$p(\boldsymbol{X} | \lambda)$を用いて、観測データ$\boldsymbol{X}$によって学習した$\lambda$の事後分布$p(\lambda | \boldsymbol{X})$はベイズの定理より
となる。分母の$p(\boldsymbol{X})$は$\lambda$に影響しないため省略して、比例関係にのみ注目する。省略した部分については、最後に正規化することで対応できる。
次にこの分布の具体的な形状を明らかにしていく。対数をとって指数部分の計算を分かりやすくして進めると
【途中式の途中式】
- 式(3.66)の下から2行目の式を用いている。
- 対数をとると積が和に変わる。また$\mu$に影響しない項($- \ln p(\boldsymbol{X})$)を$\mathrm{const.}$とおく。
- $\mathcal{N}(x_n | \mu, \lambda^{-1}),\ \mathrm{Gam}(\lambda | a, b)$に、それぞれ対数をとった1次元ガウス分布の定義式(2.65)、対数をとったガンマ分布の定義式(2.58)を用いて具体的な確率分布を代入する。
- $\sum_{n=1}^N$に関する括弧を展開し、$\lambda$と関係のない項を$\mathrm{const.}$に含める。$- \ln x^{-1} = -(-1) \ln x = \ln x$である。
- $\ln \lambda,\ \lambda$の項をそれぞれまとめて式を整理する。
となる。
式(3.67)について
とおくと、事後分布$p(\lambda | \boldsymbol{X})$は
パラメータ$\hat{a},\ \hat{b}$を持つガンマ分布であることが分かる。
また事後分布のパラメータの計算式(3.69)によって、ハイパーパラメータ$a,\ b$が更新される。
・予測分布の計算
続いて未観測のデータ$x_{*}$に対する予測分布を求めていく。
(観測データによる学習を行っていない)事前分布$p(\lambda)$を用いて、パラメータ$\lambda$を周辺化することで予測分布$p(x_{*})$となる。
3.3.1項のときと同様に、この計算を直接するのではなく予測分布$p(x_{*})$と事前分布$p(\lambda)$の関係を考えてみる。ベイズの定理より
という関係が成り立つことが分かる。両辺で対数をとり、予測分布に関して整理すると
となる。$x_{*}$と関係のない$\ln p(\lambda)$を$\mathrm{const.}$とおいている。
つまり式(3.72)の計算によっても予測分布を求められることが分かる。
$p(\lambda | x_{*})$は、1つのデータ$x_{*}$が与えられた下での$\lambda$の条件付き分布である。つまり$p(\lambda | x_{*})$は、N個の観測データが与えられた下での事後分布$p(\lambda | \boldsymbol{X})$の導出と同様の手順で求められる。
従って事後分布のパラメータ(3.69)を用いると、$p(\lambda | x_{*})$は$N = 1$より
パラメータ$\frac{1}{2} + a,\ b(x_{*})$を持つガンマ分布となる。ただし
とおく。
$p(\lambda | x_{*})$が求まったので、予測分布$p(x_{*})$の計算に戻る。
式(3.73)を式(3.72)に代入すると
となる。5行目から6行目では、$a \ln x y = a (\ln x + \ln y) = a \ln x + a \ln y$の式変形を行っている。また適宜$x_{*}$と関係のない項を$\mathrm{const.}$にまとめている。
この式の形状の確率分布はスチューデントのt分布と呼ばれる。
この式について、$\frac{1}{2 b} = \frac{a}{2 a b}$として
とおくと、予測分布$p(x_{*})$は
パラメータ$\mu_s,\ \lambda_s,\ \nu_s$を持つスチューデントのt分布となることが分かる。
予測分布のパラメータ$\lambda_s,\ \nu_s$を構成する事前分布のパラメータ$a,\ b$について、事後分布のパラメータ(3.69)を用いると
また
となり、(観測データ$\boldsymbol{X}$によって学習した)事後分布$p(\lambda | \boldsymbol{X})$を用いた予測分布$p(x_{*} | \boldsymbol{X})$
が得られる。
・Rでやってみよう
各確率分布に従いランダムに生成したデータを用いて、パラメータを推定してみましょう。
# 利用パッケージ library(tidyverse)
利用するパッケージを読み込みます。
・事後分布
## パラメータの初期値を指定 # 観測モデルのパラメータ lambda_truth <- 0.01 mu <- 25 # 事前分布のパラメータ a <- 1 b <- 1 # 試行回数 N <- 50
観測モデルの精度パラメータ$\lambda$をlambda_truthとします。この値を推定するのが目的です。
観測モデルの平均パラメータ$\mu$をmuとします。この値は推定の中で更新されません。
続いて事前分布のパラメータ(ハイパーパラメータの初期値)を指定します。$a$をa、$b$をbとします。
生成するデータ数$N$をNとします。
$\lambda$は分散の逆数であるため、値が大きくなるほど散らばり具合が小さくなります。
# ガウス分布に従うデータを生成 x_n <- rnorm(n = N, mean = mu, sd = sqrt(lambda_truth^(-1)))
ガウス分布に従う乱数を発生させる関数rnorm()を使って、ランダムにデータを生成します。
試行回数の引数nにはN、平均の引数meanにはmuを指定します。
しかし標準偏差の引数sdにはlambda_truthのままでは指定できません。そこで標準偏差$\sigma = \sqrt{\lambda^{-1}}$に置き換えるため、sqrt(lambda_truth^(-1))とします。
サンプルを確認してみましょう。
# 観測データを確認 summary(x_n) ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## -1.801 15.410 25.974 24.359 32.162 43.662
このデータを用いて事後分布のパラメータを計算します。
# 事後分布のパラメータを計算 a_hat <- N / 2 + a b_hat <- sum((x_n - mu)^2) / 2 + b
事後分布のパラメータの計算式(3.69)の計算を行い、$\hat{a},\ \hat{b}$をそれぞれa_hat, b_hatとします。
# 事後分布を計算 posterior_df <- tibble( lambda = seq(0, 3 * lambda_truth, by = 0.00001), # 作図用の値 C_G = a_hat * log(b_hat) - lgamma(a_hat), # 正規化項(対数) density = exp(C_G + (a_hat - 1) * log(lambda) - b_hat * lambda) # 確率密度 )
$\lambda$が取り得る値をseq()を使って用意します。
ただしlambda_truthの値によって必要な範囲が変わるので、ここでは0からlambda_truthの3倍までとします。そこでseq()の第1引数(最小値)に0、第2引数(最大値)に3 * mu_truthを指定します。
lambdaの各値に対してガンマ分布の定義式(2.57)の計算を行い、確率密度を求めます。ただし値が大きくなるとgamma()で計算できなくなるため、対数をとって計算することにします。なので最後にexp()で値を戻します。
計算結果は作図用にデータフレームにまとめておきます。推定結果を確認してみましょう。
head(posterior_df) ## # A tibble: 6 x 3 ## lambda C_G density ## <dbl> <dbl> <dbl> ## 1 0 149. 0. ## 2 0.00001 149. 6.32e-61 ## 3 0.00002 149. 2.06e-53 ## 4 0.00003 149. 5.06e-49 ## 5 0.00004 149. 6.53e-46 ## 6 0.00005 149. 1.68e-43
このデータフレームを用いてグラフを描きます。
# 作図 ggplot(posterior_df, aes(lambda, density)) + geom_line(color = "#56256E") + # 折れ線グラフ geom_vline(aes(xintercept = lambda_truth), color = "red", linetype = "dashed") + # 垂直線 labs(title = "Gamma Distribution", subtitle = paste0("N=", N, ", a_hat=", a_hat, ", b_hat=", round(b_hat, 2)), x = expression(lambda)) # ラベル
ggplot2パッケージを利用してグラフを作成します。
推定値は折れ線グラフgeom_line()を用いて描きます。
パラメータ$\lambda$の真の値を垂直線geom_vline()で示します。
事前分布のパラメータ$\hat{b}$はr b_hatのように細かい値となるので、round()で丸めておきます。
・予測分布
# 予測分布のパラメータを計算 mu_s <- mu lambda_s_hat <- a_hat / b_hat nu_s_hat <- 2 * a_hat
予測分布のパラメータの計算式の計算を行い、$\mu_s,\ \hat{\lambda}s,\ \hat{\mu}s$をそれぞれmu_s, lambda_s_hat, nu_s_hatとします。
lambda_s_hat <- (N + 2 * a) / (sum((x_n - mu)^2) + 2 * b) nu_s_hat <- N + 2 * a
このように、観測データx_nと事前分布のパラメータa, bを使って計算することもできます。
# 予測分布を計算 predict_df <- tibble( x = seq(-mu_s, 3 * mu_s, by = 0.01), # 作図用の値 C_St = lgamma((nu_s_hat + 1) / 2) - lgamma(nu_s_hat / 2), # 正規化項(対数) term1 = log(lambda_s_hat / pi / nu_s_hat) / 2, term2 = - (nu_s_hat + 1) / 2 * log(1 + lambda_s_hat / nu_s_hat * (x - mu_s)^2), density = exp(C_St + term1 + term2) # 確率密度 )
$x_{*}$が取り得る値をseq()を使って用意します。
ただしmu_sの値によって作図に必要な範囲が変わります。ここでは$-\mu_s$から$3 \mu_s$までとします。必要に応じて調整してください。特にこのままですと、mu_sが負の値をとるとエラーとなります。
xの各値に対してスチューデントのt分布の定義式(3.76)の計算を行い、確率密度を求めます。ちなみにこのpiは、前節で出てきたパラメータではなく円周率のことです。ご注意ください。
推定結果を確認してみましょう。
head(predict_df) ## # A tibble: 6 x 5 ## x C_St term1 term2 density ## <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 -25 1.62 -4.90 -9.51 0.00000279 ## 2 -25.0 1.62 -4.90 -9.51 0.00000280 ## 3 -25.0 1.62 -4.90 -9.50 0.00000280 ## 4 -25.0 1.62 -4.90 -9.50 0.00000281 ## 5 -25.0 1.62 -4.90 -9.50 0.00000282 ## 6 -25.0 1.62 -4.90 -9.49 0.00000283
こちらもggplot2パッケージを利用して作図します。
# 作図 ggplot(predict_df, aes(x, density)) + geom_line(color = "#56256E") + # 折れ線グラフ labs(title = "Student's t Distribution", subtitle = paste0("N=", N, ", mu_s=", mu_s, ", lambda_s_hat=", round(lambda_s_hat, 2), ", nu_s_hat=", nu_s_hat)) # ラベル
・おまけ
gganimateパッケージを利用して、パラメータの推定値の推移のgif画像を作成するためのコードです。
・コード(クリックで展開)
# 利用パッケージ library(tidyverse) library(gganimate) ## パラメータの初期値を指定 # 観測モデルのパラメータ lambda_truth <- 0.01 mu <- 25 # 事前分布のパラメータ a <- 1 b <- 1 # 試行回数 N <- 50 # 事前分布を計算 posterior_df <- tibble( lambda = seq(0, 3 * lambda_truth, by = 0.00001), # 作図用の値 C_G = a * log(b) - lgamma(a), # 正規化項(対数) density = exp(C_G + (a - 1) * log(lambda) - b * lambda), # 確率密度 N = 0 # 試行回数 ) # 初期値による予測分布のパラメータを計算 mu_s <- mu lambda_s <- a / b nu_s <- 2 * a # 初期値による予測分布を計算 predict_df <- tibble( x = seq(-mu_s, 3 * mu_s, by = 0.01), # 作図用の値 C_St = lgamma((nu_s + 1) / 2) - lgamma(nu_s / 2), # 正規化項(対数) term1 = log(lambda_s / pi / nu_s) / 2, term2 = - (nu_s + 1) / 2 * log(1 + lambda_s / nu_s * (x - mu_s)^2), density = exp(C_St + term1 + term2), # 確率密度 N = 0 # 試行回数 ) # パラメータを推定 x_n <- rep(0, N) # 受け皿 for(n in 1:N){ # ガウス分布に従うデータを生成 x_n[n] <- rnorm(n = 1, mean = mu, sd = sqrt(lambda_truth^(-1))) # ハイパーパラメータを更新 a <- 1 / 2 + a b <- (x_n[n] - mu)^2 / 2 + b # 事後分布を計算 tmp_posterior_df <- tibble( lambda = seq(0, 3 * lambda_truth, by = 0.00001), # 作図用のlambda C_G = a * log(b) - lgamma(a), # 正規化項(対数) density = exp(C_G + (a - 1) * log(lambda) - b * lambda), # 確率密度 N = n # 試行回数 ) # 予測分布のパラメータを更新 mu_s <- mu lambda_s <- a / b nu_s <- 2 * a # 予測分布を計算 tmp_predict_df <- tibble( x = seq(-mu_s, 3 * mu_s, by = 0.01), # 作図用の値 C_St = lgamma((nu_s + 1) / 2) - lgamma(nu_s / 2), # 正規化項(対数) term1 = log(lambda_s / pi / nu_s) / 2, term2 = - (nu_s + 1) / 2 * log(1 + lambda_s / nu_s * (x - mu_s)^2), density = exp(C_St + term1 + term2), # 確率密度 N = n # 試行回数 ) # 結果を結合 posterior_df <- rbind(posterior_df, tmp_posterior_df) predict_df <- rbind(predict_df, tmp_predict_df) } # 観測データを確認 summary(x_n) ## 事後分布 # 作図 posterior_graph <- ggplot(posterior_df, aes(lambda, density)) + geom_line(color = "#56256E") + # 折れ線グラフ geom_vline(aes(xintercept = lambda_truth), color = "red", linetype = "dashed") + # 垂直線 transition_manual(N) + # フレーム labs(title = "Gamma Distribution", subtitle = "N= {current_frame}", x = expression(lambda)) # ラベル # 描画 animate(posterior_graph) ## 予測分布 # 作図 predict_graph <- ggplot(predict_df, aes(x, density)) + geom_line(color = "#56256E") + # 折れ線グラフ transition_manual(N) + # フレーム labs(title = "Student's t Distribution", subtitle = "N= {current_frame}") # ラベル # 描画 animate(predict_graph)
異なる点のみを簡単に解説します。
各データによってどのように学習する(推定値が変化する)のかを確認するため、こちらはforループで1データずつ処理します。
よって生成するデータ数として設定したNがイタレーション数になります。
パラメータの推定値について、$\hat{a},\ \hat{b}$に対応するa_hat, b_hat新たに作るのではなく、a, bをイタレーションごとに更新(上書き)していきます。
それに伴い、事後分布のパラメータの計算式(3.69)の$\sum_{n=1}^N$の計算は、forループによってN回繰り返し(x_n[n] - mu)^2を加えることで行います。n回目のループ処理のときには、n-1回分の(x_n[n] - mu)^2(aの場合は1 / 2)が既にaとbに加えられているわけです。
(微妙なできだけど初期値が高すぎる所為だと思われる…。)
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
式(3.75)を導出できない。難しい訳ではないのですができませんでした…。どこがおかしいのか分かる方はぜひ教えてください。(アドバイスいただいて解けました！説明文を勘違いしてました。ありがとうございます！)
ちなみに、次でも同様に式を整理できない部分があり現在止まっております。(こっちは丁寧に書き直したらできました。)
それとは別に予測分布の方もRでやってみたいと思うので、次の項を読み(書き)終えたら一旦全編修正作業に移る予定です。(追加しました！)
2020/03/05：加筆修正しました。
【次節の内容】