はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「Rで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は3.3.1項の内容です。尤度関数・事前分布を1次元ガウス分布とした場合の平均パラメータの事後分布を導出し、学習した事後分布を用いた予測分布を導出します。またその推論過程をR言語で実装します。
省略してる内容等ありますので、本と併せて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てば幸いです。
【前節の内容】
【他の節一覧】
【この節の内容】
3.3.1 平均が未知の場合
・事後分布の計算
ガウス分布に従うと仮定するN個の観測データ$\boldsymbol{X} = \{x_1, x_2, \cdots, x_N\}$を用いて、平均パラメータ$\mu$の事後分布を求めていく。
観測モデル$p(\boldsymbol{X} | \mu)$を用いて、観測データ$\boldsymbol{X}$によって学習した$\mu$の事後分布$p(\mu | \boldsymbol{X})$はベイズの定理より
となる。分母の$p(\boldsymbol{X})$は$\mu$に影響しないため省略して、比例関係にのみ注目する。省略した部分については、最後に正規化することで対応できる。
次にこの分布の具体的な形状を明らかにしていく。対数をとって指数部分の計算を分かりやすくして進めると
【途中式の途中式】
- 式(3.49)の下から2行目の式を用いている。
- 対数をとると積が和に変わる。また$\mu$と関係のない項($- \ln p(\boldsymbol{X})$)を$\mathrm{const.}$とおく。
- $\mathcal{N}(x_n | \mu, \lambda^{-1}),\ \mathcal{N}(\mu | m, \lambda_{\mu}^{-1})$に、1次元ガウス分布の定義式(2.64)を用いて具体的な確率分布を代入する。
- それぞれ対数をとる。$\ln \sqrt{\frac{1}{x}} = \ln x^{-\frac{1}{2}} = -\frac{1}{2} \ln x$、また$\ln \exp(x) = x$である。
- $-\frac{1}{2}$を括り出し、$\mu$と関係のない項を$\mathrm{const.}$に含める。
- 2乗部分を展開する。
- $\sum_{n=1}^N$に関する括弧を展開し、$\mu$と関係のない項を$\mathrm{const.}$に含める。
- $\mu^2,\ \mu$の項をそれぞれまとめて式を整理する。
となる。
これまでと同様に、この式が(対数をとった)1次元ガウス分布の形(2.65)をしていることを確認したい。そこで
とおき、平方完成を行う。
後の項は$\mu$と無関係であるため$\mathrm{const.}$に含めて
とおくと、事後分布$p(\mu | \boldsymbol{X})$は
平均$\hat{\mu}$、精度$\hat{\lambda}_{\mu}^{-1}$の1次元ガウス分布であることが確認できる。
また事後分布のパラメータの計算式(3.53)、(3.54)によって、ハイパーパラメータ$m,\ \lambda_{\mu}$が更新される。
・予測分布の計算
続いて未観測のデータ$x_{*}$に対する予測分布を求めていく。
(観測データによる学習を行っていない)事前分布$p(\mu)$を用いて、パラメータ$\mu$を周辺化することで予測分布$p(x_{*})$となる。
この計算を直接するのではなく、より簡単に計算したい。そこで予測分布$p(x_{*})$と事前分布$p(\mu)$の関係を考えてみる。ベイズの定理より
という関係が成り立つことが分かる。両辺で対数をとり、予測分布に関して整理すると
となる。$x_{*}$と関係のない$\ln p(\mu)$を$\mathrm{const.}$とおいている。
つまり式(3.57)の計算によっても予測分布を求められることが分かる。
$p(\mu | x_{*})$は、1つのデータ$x_{*}$が与えられた下での$\mu$の条件付き分布である。つまり$p(\mu | x_{*})$は、N個の観測データが与えられた下での事後分布$p(\mu | \boldsymbol{X})$の導出と同様の手順で求められる。
従って事後分布のパラメータ(3.53)、(3.54)を用いると、$p(\mu | x_{*})$は$N = 1$より
平均$m(x_{*})$、精度$\lambda + \lambda_{\mu}$を持つ1次元ガウス分布となる。ただし
とおく。
$p(\mu | x_{*})$が求まったので、予測分布$p(x_{*})$の計算に戻る。
式(3.58)を式(3.57)に代入して、$x_{*}^2,\ x_{*}$について整理すると
となる。適宜$x_{*}$と関係のない項を$\mathrm{const.}$に含めている。
この式について
とおき、平方完成を行うと
となる。後の項は$x_{*}$と無関係であるため$\mathrm{const.}$に含めて
とおくと、予測分布$p(x_{*})$は
平均$\mu_{*}$、精度$\lambda_{*}$の1次元ガウス分布となることが分かる。
ちなみに精度パラメータ$\lambda_{*}$の逆数である分散$\sigma^2 = \lambda_{*}^{-1}$に注目すると
である。
精度パラメータ$\lambda_{*}$を構成するパラメータ$\lambda_{\mu},\ m$について、事後分布パラメータ(3.53)、(3.54)を用いると
また
となり、(観測データ$\boldsymbol{X}$によって学習した)事後分布$p(\mu | \boldsymbol{X})$を用いた予測分布$p(x_{*} | \boldsymbol{X})$
が得られる。
・Rでやってみよう
各確率分布に従いランダムに生成したデータを用いて、パラメータを推定してみましょう。
# 利用パッケージ library(tidyverse)
利用するパッケージを読み込みます。
・事後分布
## パラメータの初期値を指定 # 観測モデルのパラメータ mu_truth <- 25 lambda <- 0.01 # 事前分布のパラメータ m <- 20 lambda_mu <- 0.001 # 試行回数 N <- 50
観測モデルの平均パラメータ$\mu$をmu_truth
とします。この値を推定するのが目的です。
また観測モデルの精度パラメータ$\lambda$をlambda
とします。この値は推定の中で更新されません。
続いて事前分布のパラメータ(ハイパーパラメータの初期値)を指定します。平均$m$をm
、精度$\lambda_{\mu}$をlambda_mu
とします。
生成するデータ数$N$をN
とします。
$\lambda,\ \lambda_{\mu}$はどちらも分散の逆数であるため、値が大きくなるほど散らばり具合が小さくなります。
# ガウス分布に従うデータを生成 x_n <- rnorm(n = N, mean = mu_truth, sd = sqrt(lambda^(-1)))
ガウス分布に従う乱数を発生させる関数rnorm()
を使って、ランダムにデータを生成します。
試行回数の引数n
にはN
、平均の引数mean
にはmu_truth
を指定します。
しかし標準偏差の引数sd
にはlambda
のままでは指定できません。そこで標準偏差$\sigma = \sqrt{\lambda_{\mu}^{-1}}$に置き換えるため、sqrt(lambda^(-1))
とします。
サンプルを確認してみましょう。
# 観測データを確認 summary(x_n) ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 3.958 19.825 25.161 25.193 30.541 44.753
このデータを用いて事後分布のパラメータを計算します。
# 事後分布のパラメータを計算 lambda_mu_hat <- N * lambda + lambda_mu m_hat <- (lambda * sum(x_n) + lambda_mu * m) / lambda_mu_hat
事後分布のパラメータの計算式(3.53)、(3.54)の計算を行い、$\hat{\lambda}_{\mu},\ \hat{m}$をそれぞれlambda_mu_hat, m_hat
とします。
m_hat
の計算には、事前分布のパラメータlambda_mu
と事後分布のパラメータlambda_mu_hat
の両方を使うので注意しましょう。
# 事後分布を計算 posterior_df <- tibble( mu = seq( m_hat - 3 * sqrt(lambda_mu_hat^(-1)), m_hat + 3 * sqrt(lambda_mu_hat^(-1)), by = 0.01 ), # 作図用の値 density = dnorm(mu, mean = m_hat, sd = sqrt(lambda_mu_hat^(-1))) # 確率密度 )
$\mu$が取り得る値をseq()
を使って値を用意します。
ただしmu_truth
の値によって作図に必要な範囲が変わります。ここでは$\mu$を中心に標準偏差$\sigma$の3倍をプロットすることにします。そこでseq()
の第1引数(最小値)、第2引数(最大値)に、mu_truth
から3 * sqrt(lambda^(-1))
を引いた値と足した値をそれぞれ指定します。
mu
の各値の確率密度は、ガウス分布の定義式(2.64)で計算できます。しかし計算式が複雑なので、ここでは手を抜いてみましょう。
ガウス分布に従う確率密度を返す関数dnorm()
を使います。引数についてはrnorm()
と同様です。ただしここでは事前分布のパラメータm
とlambda_mu
を用います。
前節で扱った確率分布についてもこのような関数が用意されています。
計算結果は作図用にデータフレームにまとめておきます。推定結果を確認してみましょう。
head(posterior_df) ## # A tibble: 6 x 2 ## mu density ## <dbl> <dbl> ## 1 20.9 0.00314 ## 2 21.0 0.00320 ## 3 21.0 0.00327 ## 4 21.0 0.00334 ## 5 21.0 0.00341 ## 6 21.0 0.00349
このデータフレームを用いてグラフを描きます。
# 作図 ggplot(posterior_df, aes(mu, density)) + geom_line(color = "#56256E") + # 折れ線グラフ geom_vline(aes(xintercept = mu_truth), color = "red", linetype = "dashed") + # 垂直線 labs(title = "Gaussian Distribution", subtitle = paste0("N=", N, ", m_hat=", round(m_hat, 2), ", lambda_mu_hat=", lambda_mu_hat), x = expression(mu)) # ラベル
ggplot2
パッケージを利用してグラフを作成します。
推定値は折れ線グラフgeom_line()
を用いて描きます。
パラメータ$\mu$の真の値を垂直線geom_vline()
で示します。
m_hat
はr m_hat
のように細かい値となるので、round()
で丸めておきます。
・予測分布
# 予測分布のパラメータを計算 lambda_hat <- lambda * lambda_mu_hat / (lambda + lambda_mu_hat) mu_hat <- m_hat
予測分布のパラメータの計算式の計算を行い、$\hat{\lambda}_{*},\ \hat{\mu}_{*}$をそれぞれlambda_hat, mu_hat
とします。
lambda_hat <- (N * lambda + lambda_mu) * lambda / ((N + 1) * lambda + lambda_mu) mu_hat <- (lambda * sum(x_n) + lambda_mu * m) / (N * lambda + lambda_mu)
このように、観測データx_n
と観測モデルと事前分布のパラメータlambda, lambda_mu
を使って計算することもできます。
# 予測分布の計算 predict_df <- tibble( x = seq( mu_hat - 3 * sqrt(lambda_hat^(-1)), mu_hat + 3 * sqrt(lambda_hat^(-1)), by = 0.01 ), # 作図用の値 density = dnorm(x, mean = mu_hat, sd = sqrt(lambda_hat^(-1))) # 確率密度 )
$x_{*}$の分布も事後分布のときと同様にして計算します。
推定結果を確認してみましょう。
head(predict_df) ## # A tibble: 6 x 2 ## x density ## <dbl> <dbl> ## 1 -5.12 0.000439 ## 2 -5.11 0.000440 ## 3 -5.10 0.000441 ## 4 -5.09 0.000443 ## 5 -5.08 0.000444 ## 6 -5.07 0.000445
こちらもggplot2
パッケージを利用して作図します。
# 作図 ggplot(predict_df, aes(x, density)) + geom_line(color = "#56256E") + # 折れ線グラフ geom_vline(aes(xintercept = mu_truth), color = "red", linetype = "dashed") + # 垂直線 labs(title = "Gaussian Distribution", subtitle = paste0("N=", N, ", m_hat=", round(mu_hat, 2), ", lambda_mu_hat=", round(lambda_hat, 2)), x = expression(mu)) # ラベル
・おまけ
gganimate
パッケージを利用して、パラメータの推定値の推移のgif画像を作成するためのコードです。
・コード(クリックで展開)
# 利用パッケージ library(tidyverse) library(gganimate) ## パラメータの初期値を指定 # 観測モデルのパラメータ mu_truth <- 25 lambda <- 0.01 # 事前分布のパラメータ m <- 20 lambda_mu <- 0.001 # 試行回数 N <- 50 # 事前分布を計算 posterior_df <- tibble( mu = seq( mu_truth - 2 * sqrt(lambda^(-1)), mu_truth + 2 * sqrt(lambda^(-1)), by = 0.01 ), # 作図用の値 density = dnorm(mu, mean = m, sd = sqrt(lambda_mu^(-1))), # 確率密度 N = 0 # 試行回数 ) # 初期値による予測分布のパラメータを計算 lambda_star <- lambda * lambda_mu / (lambda + lambda_mu) mu_star <- m # 初期値による予測分布を計算 predict_df <- tibble( x = seq( mu_truth - 3 * sqrt(lambda^(-1)), mu_truth + 3 * sqrt(lambda^(-1)), by = 0.01 ), # 作図用の値 density = dnorm(x, mean = mu_star, sd = sqrt(lambda_star^(-1))), # 確率密度 N = 0 # 試行回数 ) # パラメータを推定 x_n <- rep(0, N) # 受け皿 for(n in 1:N){ # ガンマ分布に従うデータを生成 x_n[n] <- rnorm(n = 1, mean = mu_truth, sd = sqrt(lambda^(-1))) # パラメータを更新 lambda_mu_old <- lambda_mu lambda_mu <- 1 * lambda + lambda_mu m <- (lambda * x_n[n] + lambda_mu_old * m) / lambda_mu # 事後分布を計算 tmp_posterior_df <- tibble( mu = seq( mu_truth - 2 * sqrt(lambda^(-1)), mu_truth + 2 * sqrt(lambda^(-1)), by = 0.01 ), # 作図用の値 density = dnorm(mu, mean = m, sd = sqrt(lambda_mu^(-1))), # 確率密度 N = n # 試行回数 ) # 予測分布のパラメータを計算 lambda_star <- lambda * lambda_mu / (lambda + lambda_mu) mu_star <- m # 予測分布を計算 tmp_predict_df <- tibble( x = seq( mu_truth - 3 * sqrt(lambda^(-1)), mu_truth + 3 * sqrt(lambda^(-1)), by = 0.01 ), # 作図用の値 density = dnorm(x, mean = mu_star, sd = sqrt(lambda_star^(-1))), # 確率密度 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(mu, density)) + geom_line(color = "#56256E") + # 折れ線グラフ geom_vline(aes(xintercept = mu_truth), color = "red", linetype = "dashed") + # 垂直線 transition_manual(N) + # フレーム labs(title = "Gaussian Distribution", subtitle = "N= {current_frame}", x = expression(mu)) # ラベル # 描画 animate(posterior_graph) ## 予測分布 predict_graph <- ggplot(predict_df, aes(x, density)) + geom_line(color = "#56256E") + # 折れ線グラフ geom_vline(aes(xintercept = mu_truth), color = "red", linetype = "dashed") + # 垂直線 transition_manual(N) + # フレーム labs(title = "Gaussian Distribution", subtitle = "N= {current_frame}", x = expression(x)) # ラベル # 描画 animate(predict_graph)
異なる点のみを簡単に解説します。
各データによってどのように学習する(推定値が変化する)のかを確認するため、こちらはforループで1データずつ処理します。
よって生成するデータ数として設定したN
がイタレーション数になります。
パラメータの推定値について、$\hat{m},\ \hat{\lambda}_{\mu}$に対応するm_hat, lambda_mu_hat
を新たに作るのではなく、m, lambda_mu
をイタレーションごとに更新(上書き)していきます。
それに伴い、事後分布のパラメータの計算式(3.54)、(3.53)の$\sum_{n=1}^N$の計算は、forループによってN回繰り返しx_n[n]
を加えることで行います。n回目のループ処理のときには、n-1回分のx_n[n]
(lambda_mu
の場合はlambda
)が既にm
とlambda_mu
に加えられているわけです。
$\hat{m}$の計算には、更新前の$\lambda_{\mu}$と更新後の$\hat{\lambda}_{\mu}$の両方が必要です。そこで更新前のlambda_mu
の値をlambda_mu_old
として残しておきます。そしてlambda_mu
を更新してから、lambda_mu_old
とlambda_mu
を使ってm
を更新します。
作図範囲を固定するため(更新されない)観測モデルのパラメータを使います。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
ここまでは1日1項ペースで進んでいましたが、次とその次でかなり詰まっております…
2020/03/05:加筆修正しました。
【次節の内容】