はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、3.5節の内容です。線形回帰モデルの「観測モデルを1次元ガウス分布(1次元正規分布)」、「事前分布を多次元ガウス分布(多変量正規分布)」とした場合の「パラメータの事後分布」と「未観測値の予測分布」の計算をRで実装します。
省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【数式読解編】
【他の節一覧】
【この節の内容】
・Rでやってみよう
人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。
利用するパッケージを読み込みます。
# 利用パッケージ library(tidyverse) library(mvtnorm)
mvtnormパッケージは、多次元ガウス分布に関するパッケージです。多次元ガウス分布に従う乱数生成関数rmvnorm()を使います。
・観測モデルの設定
この例では、ある入力$x_n$を$\mathbf{x}_n = (1, x_n, x_n^2, \cdots, x_n^{M-1})$とした$M$次元ベクトルを入力値として用いて、$y_n$を出力します。
そこでまずは、$x_n$から$\mathbf{x}_n$を作成する関数を定義しておきます。
# {x^(m-1)}ベクトル作成関数を定義 x_vector <- function(x_n, M) { # 受け皿を作成 x_nm <- matrix(NA, nrow = length(x_n), ncol = M) # m乗を計算 for(m in 1:M) { x_nm[, m] <- x_n^(m-1) } return(x_nm) }
$x^0 = 1$であることを利用して、$m$番目(列目)の項を$x_n^{m-1}$で計算します。
$N$個のデータ$\{x_1, x_2, \cdots, x_N\}$の場合は、それぞれ$M$次元にした$N \times M$の行列
として扱うことにします($M \times N$のマトリクスとした方が実装(内積の処理)が楽ですが、これまでの章での扱いと合わせます。多くの計算でt()で転置することになります)。
試しに、1から5の値をx_vector()で4次元ベクトルにしてみます。
# 確認 x_vector(1:5, 4)
## [,1] [,2] [,3] [,4] ## [1,] 1 1 1 1 ## [2,] 1 2 4 8 ## [3,] 1 3 9 27 ## [4,] 1 4 16 64 ## [5,] 1 5 25 125
1列目は全て1になり、2列目が元の値です。
続いて、真のモデルを確認していきます。
入力$\mathbf{x} = (1, x, x^2, \cdots, x^{M-1})$とパラメータ$\mathbf{w} = (w_1, w_2, \cdots, w_M)$の内積(重み付き和)
を真のモデルとします。
実際に計算して、グラフで確認しましょう。
観測モデルのパラメータ$\mathbf{w}$を設定します。
# 真の次元数を指定 M_truth <- 4 # 真のパラメータを生成 w_truth_m <- sample(x = seq(-1, 1, by = 0.1), size = M_truth) # 確認 w_truth_m
## [1] 0.8 0.6 0.4 -0.1
この例では、真のパラメータの次元数をM_truthとして値を指定して、$\mathbf{w}$の値を-1から1の範囲でランダムに生成することにします。作成した真のパラメータをw_truth_mとします。
作図用に、$x_n$がとり得る値を作成します。
# 作図用のx軸の値を作成 x_line <- seq(-3, 3, by = 0.01) # 作図用のxをM次元に拡張 x_truth_mat <- x_vector(x_line, M_truth)
$x_n$がとり得る値(x軸の値)をseq()で作成してx_lineとします。この例では、-3から3の範囲をグラフにすることにします(範囲を広くすると$x_n^{M-1}$の影響が強くなり、パラメータに関わらず似た形のグラフになってつまんないです)。
また、x_lineの各要素をx_vector()でM_truth次元ベクトルに変換してx_truth_matとします。
# 確認 x_truth_mat[1:5, ]
## [,1] [,2] [,3] [,4] ## [1,] 1 -3.00 9.0000 -27.00000 ## [2,] 1 -2.99 8.9401 -26.73090 ## [3,] 1 -2.98 8.8804 -26.46359 ## [4,] 1 -2.97 8.8209 -26.19807 ## [5,] 1 -2.96 8.7616 -25.93434
2列目がx_lineの要素ですね。
真のモデルの出力を計算します。
y_line <- t(w_truth_m) %*% t(x_truth_mat) %>% as.vector()
作図用に、真のモデルの入出力をデータフレームに格納します。
# 真のモデルをデータフレームに格納 model_df <- tibble( x = x_line, y = y_line ) # 確認 head(model_df)
## # A tibble: 6 x 2 ## x y ## <dbl> <dbl> ## 1 -3 5.3 ## 2 -2.99 5.26 ## 3 -2.98 5.21 ## 4 -2.97 5.17 ## 5 -2.96 5.12 ## 6 -2.95 5.08
ggplot2で作図するにはデータフレームを渡す必要があります。
真のモデルを作図します。
# 真のモデルを作図 ggplot(model_df, aes(x = x, y = y)) + geom_line() + labs(title = "Observation Model", subtitle = paste0("w=(", paste0(w_truth_m, collapse = ', '), ")"))
これがノイズを含めない観測モデルです。
・データの生成
実際に観測される出力値$y_n$にはノイズ成分$\epsilon_n$が含まれているとします。
ノイズ成分は、平均$0$、分散$\sigma^2 = \lambda^{-1}$の1次元のガウス分布に従うと仮定します。分散の逆数$\lambda = \frac{1}{\sigma^2}$を精度と呼びます。
ノイズ成分のパラメータを設定します。
# ノイズ成分の標準偏差を指定 sigma <- 1.5 # ノイズ成分の精度を計算 lambda <- 1 / sigma^2 # 確認 lambda
## [1] 0.4444444
1次元ガウス分布の標準偏差$\sigma$をsigmaとして値を指定し、精度を計算してlambdaとします。これは作図時に標準偏差を利用するためです。
試しに、ノイズ成分を1つ生成してみます。
# ノイズを1つ生成 rnorm(n = 1, mean = 0, sd = sqrt(1 / lambda))
## [1] 0.7615452
1次元ガウス分布に従う乱数は、rnorm()で生成できます。データ数の引数nにはデータ数、平均の引数meanに0、標準偏差の引数sdにsqrt(1 / lambda))を指定します。精度$\lambda$を標準偏差$\sigma = \sqrt{\frac{1}{\lambda}}$に戻しています(勿論sigmaを指定すればいいのですが、式(3.142)に合わせておきます)。
以上で観測モデルのパラメータを設定できたので、入力値$\mathbf{X}$を作成して、出力値$\mathbf{y} = \{y_1, y_2, \cdots, y_N\}$を求めます。
$N$個の入力値を生成します。
# (観測)データ数を指定 N <- 50 # 入力値を生成 x_n <- sample(seq(min(x_line), max(x_line), by = 0.01), size = N, replace = TRUE) # 入力値をM次元に拡張 x_truth_nm <- x_vector(x_n, M_truth)
生成するデータ数$N$をNとして値を指定します。
データの生成に関して分布の仮定をおいていないため、sample()で一様な確率に従いランダムに値を生成します。ただし$x_n$がとり得る範囲は、グラフの描画範囲x_lineにしておきます。生成したN個のデータをx_nとします。
また、x_nの各要素をx_vector()でM_truth次元ベクトルに変換してx_truth_nmとします。
1次元ガウス分布に従う$N$個のノイズ成分を生成して、出力値$\mathbf{y}$を計算します。
# ノイズ成分を生成 epsilon_n <- rnorm(n = N, mean = 0, sd = sqrt(1 / lambda)) # 出力値を計算:式(3.141) y_n <- (t(w_truth_m) %*% t(x_truth_nm) + epsilon_n) %>% as.vector()
観測モデル(3.141)に従い、出力値$\mathbf{y}$を計算してy_nとします。
作図用に、観測データ$\mathbf{y},\ \mathbf{X}$をデータフレームに格納します。
# 観測データをデータフレームに格納 sample_df <- tibble( x = x_n, y = y_n ) # 確認 head(sample_df)
## # A tibble: 6 x 2 ## x y ## <dbl> <dbl> ## 1 1.27 2.76 ## 2 -2.02 2.51 ## 3 -1.2 1.80 ## 4 0.45 0.533 ## 5 2.19 3.90 ## 6 2.95 4.66
観測データの散布図を真のモデルのグラフと重ねて確認します。
# 観測データの散布図を作成 ggplot() + geom_point(data = sample_df, aes(x, y)) + # 観測データ geom_line(data = model_df, aes(x, y), color = "blue") + # 真のモデル labs(title = "Observation Data", subtitle = paste0("N=", N, ", w=(", paste0(round(w_truth_m, 1), collapse = ', '), ")", ", sigma=", round(sqrt(1 / lambda), 1)))
以上が観測モデルに関する処理です。
・事前分布の設定
次に、観測モデルのパラメータ$\mathbf{w}$の事前分布$p(\mathbf{w})$として、多次元ガウス分布$\mathcal{N}(\mathbf{w} | \mathbf{m}, \boldsymbol{\Lambda}^{-1})$を設定します。
$\mathbf{w}$の事前分布のパラメータ(超パラメータ)を設定します。
# 事前分布の次元数を指定 M <- 5 # 事前分布の平均を指定 m_m <- rep(0, M) # 事前分布の精度行列を指定 sigma_mm <- diag(M) * 10 lambda_mm <- solve(sigma_mm^2) # 確認 lambda_mm
## [,1] [,2] [,3] [,4] [,5] ## [1,] 0.01 0.00 0.00 0.00 0.00 ## [2,] 0.00 0.01 0.00 0.00 0.00 ## [3,] 0.00 0.00 0.01 0.00 0.00 ## [4,] 0.00 0.00 0.00 0.01 0.00 ## [5,] 0.00 0.00 0.00 0.00 0.01
事前分布の次元数$M$をMとして値を指定します。(ところで入力値は$x_n$と$\mathbf{x}_n$のどちらの状態で得られるものなのでしょうかね。$M$次元ベクトル$\mathbf{x}_n$が手元にある場合は、真の$M$を事前分布の次元数とします。この例では、真の次元数が分からない$x_n$しか手元にないものとして、任意の値を指定することにします。)
多次元ガウス分布の平均パラメータ$\mathbf{m} = (m_1, \cdots, m_M)$をm_m、標準偏差$\sigma_{i,i}$と相関係数$\sigma_{i,j}$の行列$(\sigma_{1,1}, \cdots, \sigma_{M,M})$をsigma_mmとして値を指定します。sigma_mmから分散共分散行列$\boldsymbol{\Sigma} = (\sigma_{1,1}^2, \cdots, \sigma_{M,M}^2)$の逆行列である精度行列(精度パラメータ)$\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{-1}$を計算してlambda_mmとします。逆行列はsolve()で計算します。lambda_mmを直接指定してもかまいません。diag()は単位行列を作成する関数です。
また、入力値に関するデータx_n, x_lineをM次元ベクトルに変換します。
# 入力値をM次元に拡張 x_nm <- x_vector(x_n, M) # 作図用のxをM次元に拡張 x_mat <- x_vector(x_line, M)
それぞれx_nm, x_lineとします。
事前分布から$\mathbf{w}$をサンプリングしたモデルを確認しましょう。
# 事前分布からサンプリングしたwを用いたモデルを比較 prior_df <- tibble() for(i in 1:5) { # サンプルサイズを指定 # パラメータを生成 prior_w_m <- mvtnorm::rmvnorm(n = 1, mean = m_m, sigma = solve(lambda_mm)) %>% as.vector() # 入出力値をデータフレームに格納 tmp_df <- tibble( x = x_line, y = as.vector(t(prior_w_m) %*% t(x_mat)), smp_num = as.factor(i) # サンプル番号 ) # 結果を結合 prior_df <- rbind(prior_df, tmp_df) }
真のパラメータはsample()で生成しましたが、事前分布からのサンプリングでは仮定した$M$次元ガウス分布から生成します。多次元ガウス分布に従う乱数は、mvtnorm::rmvnorm()で生成できます。データ数の引数nに1、平均の引数meanにm_m、分散共分散行列の引数sigmaにsolve(lambda_mm)を指定します。
生成したパラメータprior_w_mを用いて出力を計算して、データフレームに格納します。
サンプリングした$\mathbf{w}$を用いたモデルを作図します。
# 事前分布からサンプリングしたパラメータによるモデルを作図 ggplot() + geom_line(data = prior_df, aes(x = x, y = y, color = smp_num)) + # サンプリングしたwを用いたモデル geom_line(data = model_df, aes(x =x, y = y), color = "blue", linetype = "dashed") + # 真のモデル ylim(min(model_df[["y"]]) - 3 * sigma, max(model_df[["y"]]) + 3 * sigma) + # y軸の表示範囲 labs(title = "Sampling from Prior Distribution", subtitle = paste0("m=(", paste0(m_m, collapse = ", "), ")"))
この時点では当然観測モデルとは程遠いグラフになります。
・事後分布の計算
観測データ$\mathbf{y},\ \mathbf{X}$からパラメータ$\mathbf{w}$の事後分布$p(\mathbf{w} | \mathbf{y}, \mathbf{X})$を求めます(パラメータ$\mathbf{w}$を分布推定します)。事後分布は多次元ガウス分布$\mathcal{N}(\mathbf{w} | \hat{\mathbf{m}}, \hat{\boldsymbol{\Lambda}}^{-1})$になります。
観測データy_n、x_nmを用いて、$\mathbf{w}$の事後分布のパラメータを計算します。
# 事後分布のパラメータを計算:式(3.148) lambda_hat_mm <- lambda * t(x_nm) %*% x_nm + lambda_mm m_hat_m <- solve(lambda_hat_mm) %*% (lambda * t(t(y_n) %*% x_nm) + lambda_mm %*% matrix(m_m)) %>% as.vector()
事後分布パラメータを
を計算して、結果をlambda_hat_mm、m_hat_mとします。ただし効率よく処理するために、$\sum_{n=1}^N \mathbf{x}_n \mathbf{x}_n^{\top}$の計算を$\mathbf{X}^{\top} \mathbf{X}$で、$\sum_{n=1}^N y_n \mathbf{x}_n$の計算を$\mathbf{y} \mathbf{X}$で計算しています。
# 確認 lambda_hat_mm; m_hat_m
## [,1] [,2] [,3] [,4] [,5] ## [1,] 22.23222 9.72000 53.27987 60.28596 288.3090 ## [2,] 9.72000 53.28987 60.28596 288.30899 399.3105 ## [3,] 53.27987 60.28596 288.31899 399.31048 1915.8348 ## [4,] 60.28596 288.30899 399.31048 1915.84481 2774.4289 ## [5,] 288.30899 399.31048 1915.83481 2774.42892 14009.2688 ## [1] 0.95775178 0.75846328 0.50341487 -0.08903103 -0.02894226
事後分布からも$\mathbf{w}$をサンプリングしてモデルを確認しましょう。
# 事後分布からサンプリングしたwを用いたモデルを比較 posterior_df <- tibble() for(i in 1:5) { # サンプルサイズを指定 # パラメータを生成 posterior_w_m <- mvtnorm::rmvnorm(n = 1, mean = m_hat_m, sigma = solve(lambda_hat_mm)) %>% as.vector() # 入出力値をデータフレームに格納 tmp_df <- tibble( x = x_line, y = as.vector(t(posterior_w_m) %*% t(x_mat)), smp_num = as.factor(i) # サンプル番号 ) # 結果を結合 posterior_df <- rbind(posterior_df, tmp_df) }
更新した超パラメータm_m, lambda_hat_mmを用いて、事前分布のときと同様に処理します。
サンプリングした$\mathbf{w}$を用いたモデルを作図します。
# 事後分布からサンプリングしたパラメータによるモデルを作図 ggplot() + geom_line(data = posterior_df, aes(x = x, y = y, color = smp_num)) + # サンプリングしたwを用いたモデル geom_line(data = model_df, aes(x = x, y = y), color = "blue", linetype = "dashed") + # 真のモデル geom_point(data = sample_df, aes(x, y)) + # 観測データ ylim(min(model_df[["y"]]) - 3 * sigma, max(model_df[["y"]]) + 3 * sigma) + # y軸の表示範囲 labs(title = "Sampling from Posterior Distribution", subtitle = paste0("N=", N, ", m_hat=(", paste0(round(m_hat_m, 2), collapse = ", "), ")"))
真のモデルに近い形になっているのを確認できます。
パラメータ$\mathbf{w}$の期待値($\mathbf{w}$の事後分布の期待値)$\mathbb{E}[\mathbf{w}]$を確認しましょう。
# 真のパラメータを確認 w_truth_m # パラメータの期待値を確認 round(m_hat_m, 2)
## [1] 0.8 0.6 0.4 -0.1 ## [1] 0.96 0.76 0.50 -0.09 -0.03
パラメータの期待値は多次元ガウス分布の期待値(2.76)より、$\mathbb{E}[\mathbf{w}] = \hat{\mathbf{m}}$です。
真のパラメータw_truth_mとパラメータの期待値m_hat_mを比べると、真の値に近い値を推定できているのが分かります。真のパラメータにはない5次元目についても0に近い値になっています。
・予測分布の計算
最後に、観測データ$\mathbf{y},\ \mathbf{X}$と新規の入力値$\mathbf{x}_{*}$から未観測値$y_{*}$の予測分布$p(y_{*} | \mathbf{x}_{*}, \mathbf{y}, \mathbf{X})$を求めます。予測分布は1次元ガウス分布$\mathcal{N}(y_{*} | \hat{\mu}_{*}, \hat{\lambda}_{*})$になります。
$\mathbf{w}$の事後分布パラメータを用いて、予測分布のパラメータを計算します。
# 予測分布のパラメータを計算:式(3.155') mu_star_hat_line <- t(m_hat_m) %*% t(x_mat) %>% as.vector() sigma2_star_hat_line <- 1 / lambda + x_mat %*% solve(lambda_hat_mm) %*% t(x_mat) %>% diag()
x_lineの要素ごとに、予測分布パラメータ
を計算します。
入力(x軸)の各値に対する出力(y軸)の平均$\hat{\mu}_{*}$をmu_star_hat_lineとします。$\hat{\mu}_{*}$の計算式を見ると、真のパラメータ$\mathbf{w}$の代わりにパラメータの期待値$\mathbb{E}[\mathbf{w}] = \hat{\mathbf{m}}$を用いて、ノイズのない出力の計算$\mathbb{E}[\mathbf{w}]^{\top} \mathbf{x}_{*}$をしているのが分かります。また、1次元ガウス分布の期待値(2.66)より、出力の期待値$\mathbb{E}[y_{*}] = \hat{\mu}_{*}$なので、$\mathbf{w}$の期待値を用いて$y_{*}$の期待値を求めているとも言えます。
また、精度の逆数は分散$\hat{\sigma}_{*}^2 = \hat{\lambda}_{*}^{-1}$なのでsigma2_star_hat_lineとします。
# 確認 mu_star_hat_line[1:5]; sigma2_star_hat_line[1:5]
## [1] 3.272610 3.257184 3.241707 3.226182 3.210610 ## [1] 4.375009 4.310422 4.247926 4.187465 4.128984
計算結果をデータフレームに格納します。
# 予測分布をデータフレームに格納 predict_df <- tibble( x = x_line, E_y = mu_star_hat_line, # 予測分布の期待値 minus_sigma_y = E_y - sqrt(sigma2_star_hat_line), # 予測分布の期待値-sigma plus_sigma_y = E_y + sqrt(sigma2_star_hat_line) # 予測分布の期待値+sigma )
作図用に、y軸の平均mu_star_hat_lineをE_y、平均から標準偏差(分散の平方根)sqrt(sigma2_star_hat_line)を足し引きしたものをそれぞれplus_sigma_y、minus_sigma_yとして、データフレームを作成します。
# 確認 head(predict_df)
## # A tibble: 6 x 4 ## x E_y minus_sigma_y plus_sigma_y ## <dbl> <dbl> <dbl> <dbl> ## 1 -3 3.27 1.18 5.36 ## 2 -2.99 3.26 1.18 5.33 ## 3 -2.98 3.24 1.18 5.30 ## 4 -2.97 3.23 1.18 5.27 ## 5 -2.96 3.21 1.18 5.24 ## 6 -2.95 3.19 1.18 5.21
予測分布を真のモデルと重ねて作図します。
# 予測分布を作図 ggplot() + geom_line(data = predict_df, aes(x = x, y = E_y), color = "orange") + # 予測分布の期待値 geom_ribbon(data = predict_df, aes(x = x, ymin = minus_sigma_y, ymax = plus_sigma_y), fill = "#00A968", alpha = 0.1, color = "#00A968", linetype = "dotted") + # 予測分布の標準偏差 geom_line(data = model_df, aes(x = x, y = y), color = "blue", linetype = "dashed") + # 真のモデル geom_point(data = sample_df, aes(x = x, y = y)) + # 観測データ ylim(min(model_df[["y"]]) - 3 * sigma, max(model_df[["y"]]) + 3 * sigma) + # y軸の表示範囲 labs(title = "Predictive Distribution", subtitle = paste0("N=", N, ", m_hat=(", paste0(round(m_hat_m, 2), collapse = ", "), ")"), y = "y")
標準偏差の範囲をgeom_ribbon()で表示します。
観測データが増えると予測分布の平均(オレンジの実線)が真のモデル(青の破線)に近付き、また精度の範囲(緑の塗りつぶし)が狭まります。
x軸の各点$x_{*}$において、y軸方向に平均$\hat{\mu}_{*}$、精度$\hat{\lambda}_{*}$の1次元ガウス分布になっています。
・過学習の確認
いくつか値を変えて試してみましょう。
まずは、データ数$N$を減らして次元数$M$を増やした場合を見てみましょう。真の次元数4に対して、$N = 10$、$M = 10$とします。
観測モデル(青の破線)への当てはまりが悪くなり、また精度も悪くなっています(標準偏差の範囲が広くなっています)。パラメータ数(次元数)が増えたことでモデル(オレンジの実線)の表現力が上がり、データへの当てはまりが過剰に良くなっています。
次は、x軸の描画範囲を広くして、観測データのない範囲に対する予測分布を確認してみましょう。$N = 50$、$M = 10$として、またx_nの作成時の範囲をx_lineの範囲よりも狭くします。
$M$の値が真の値よりも大きいと、観測データの範囲を越えたとたんに予測できていません。
真の値と同じ$M = 4$とします。
$M$が真の値と一致していると、予測ができていますね。m_hat_mの値がw_truth_mに近付いています。
・おまけ：アニメーションで推移の確認
gganimateパッケージを利用して、gif画像を作成するコードです。
・コード(クリックで展開)
異なる点のみを簡単に解説します。
# 利用するパッケージ library(tidyverse) library(mvtnorm) library(gganimate)
# {x^(m-1)}ベクトル作成関数を定義 x_vector <- function(x_n, M) { # 受け皿を作成 x_nm <- matrix(NA, nrow = length(x_n), ncol = M) # m乗を計算 for(m in 1:M) { x_nm[, m] <- x_n^(m-1) } return(x_nm) }
・観測データ数が予測精度へ与える影響
まずは、観測データの増加が予測精度へ与える影響を確認します。
・コード(クリックで展開)
・観測モデルの設定
# 真の次元数を指定 M_truth <- 4 # ノイズ成分の標準偏差を指定 sigma <- 1.5 # ノイズ成分の精度を計算 lambda <- 1 / sigma^2 # パラメータを生成 w_truth_m <- sample(x = seq(-1, 1, by = 0.1), size = M_truth) # 作図用のx軸の値を作成 x_line <- seq(-3, 3, by = 0.01) # 観測モデルをデータフレームに格納 model_df <- tibble( x = x_line, y = t(w_truth_m) %*% t(x_vector(x_line, M_truth)) %>% as.vector() )
・事前分布の設定
# 事前分布の次元数を指定 M <- 5 # 事前分布のパラメータを指定 m_m <- rep(0, M) lambda_mm <- diag(M) * 0.01 # 作図用のxをM次元に拡張 x_mat <- x_vector(x_line, M) # 初期値による予測分布のパラメータを計算:式(3.155) mu_star_line <- t(m_m) %*% t(x_mat) %>% as.vector() sigma2_star_line <- 1 / lambda + x_mat %*% solve(lambda_mm) %*% t(x_mat) %>% diag() # 初期値による予測分布をデータフレームに格納 predict_df <- tibble( x = x_line, E_y = mu_star_line, minus_sigma_y = E_y - sqrt(sigma2_star_line), plus_sigma_y = E_y + sqrt(sigma2_star_line), label = as.factor( paste0( "N=", 0, ", m=(", paste(m_m, collapse = ", "), ")" ) ) )
各試行の結果を同じデータフレームに格納していく必要があります。予測分布をpredict_dfとして、初期値の結果を持つように作成しておきます。
・推論処理
# (観測)データ数(試行回数)を指定 N <- 100 # 受け皿を初期化 x_n <- rep(NA, N) y_n <- rep(NA, N) x_df <- tibble( x = NA, y = NA, label = as.factor( paste0("N=", 0, ", m=(", paste(m_m, collapse = ", "), ")") ) ) # ベイズ推論 for(n in 1:N) { # 入力値を生成 x_n[n] <- sample(seq(min(x_line), max(x_line), by = 0.01), size = 1, replace = TRUE) # 出力値を計算:式(3.141) y_n[n] <- t(w_truth_m) %*% t(x_vector(x_n[n], M_truth)) + rnorm(n = 1, mean = 0, sd = sqrt(1 / lambda)) # 入力値をM次元に拡張 x_1m <- x_vector(x_n[n], M) # 事後分布のパラメータを更新:式(3.148) old_lambda_mm <- lambda_mm lambda_mm <- lambda * t(x_1m) %*% x_1m + lambda_mm m_m <- solve(lambda_mm) %*% (lambda * y_n[n] * t(x_1m) + old_lambda_mm %*% matrix(m_m)) %>% as.vector() # 予測分布のパラメータを更新:式(3.155') mu_star_line <- t(m_m) %*% t(x_mat) %>% as.vector() sigma2_star_line <- 1 / lambda + x_mat %*% solve(lambda_mm) %*% t(x_mat) %>% diag() # n回目までの観測データをデータフレームに格納 tmp_x_df <- tibble( x = x_n[1:n], y = as.vector(y_n[1:n]), label = as.factor( paste0( "N=", n, ", m_hat=(", paste(round(m_m, 2), collapse = ", "), ")" ) ) ) # 結果を結合 x_df <- rbind(x_df, tmp_x_df) # 予測分布を計算 tmp_predict_df <- tibble( x = x_line, E_y = mu_star_line, minus_sigma_y = E_y - sqrt(sigma2_star_line), plus_sigma_y = E_y + sqrt(sigma2_star_line), label = as.factor( paste0( "N=", n, ", m_hat=(", paste(round(m_m, 2), collapse = ", "), ")" ) ) ) # 結果を結合 predict_df <- rbind(predict_df, tmp_predict_df) }
観測された各データによってどのように学習する(分布が変化する)のかを確認するため、for()で1データずつ処理します。よって、データ数Nがイタレーション数になります。
一度の処理で事後分布のパラメータを計算するのではなく、事前分布(1ステップ前の事後分布)に対して繰り返し観測データの情報を与えることでパラメータを更新(上書き)していきます。
それに伴い、事後分布のパラメータの計算式(3.148)の$\sum_{n=1}^N \mathbf{x}_n$の計算は、ループ処理によって$N$回繰り返しx_nd[n, ]を加えることで行います。$n$回目のループ処理のときには、$n-1$回分のx_nd[n, ]が既にlambda_mmとm_mに加えられているわけです。
ただし、事後分布のパラメータの計算において、更新前(1ステップ前)のパラメータ$\boldsymbol{\Lambda}$を使うため、old_lambda_mmとして値を一時的に保存しておきます。
・作図処理
# 予測分布を作図 predict_graph <- ggplot() + geom_line(data = predict_df, aes(x = x, y = E_y), color = "orange") + # 予測分布の期待値 geom_ribbon(data = predict_df, aes(x = x, ymin = minus_sigma_y, ymax = plus_sigma_y), fill = "#00A968", alpha = 0.1, color = "#00A968", linetype = "dotted") + # 予測分布の標準偏差 geom_line(data = model_df, aes(x = x, y = y), color = "blue", linetype = "dashed") + # 真のモデル geom_point(data = x_df, aes(x = x, y = y)) + # 観測データ gganimate::transition_manual(label) + # フレーム ylim(min(model_df[["y"]]) - 3 * sigma, max(model_df[["y"]]) + 3 * sigma) + # y軸の表示範囲 labs(title = "Predictive Distribution", subtitle = "{current_frame}", y = "y") # gif画像を作成 gganimate::animate(predict_graph, nframes = N, fps = 10)
各フレームの順番を示す列をgganimate::transition_manual()に指定します。初期値(事前分布)を含むため、フレーム数はN + 1です。
(わざわざ学習してるっぽい振る舞いをさせていますが、これが良くないことは先ほど確認した通りです。)
・事前分布の次元数が予測精度へ与える影響
次に、同じ観測データを用いたとき、パラメータの数(事前分布の次元数)の違いが予測精度へ与える影響を確認します。
・コード(クリックで展開)
・観測モデルの設定
# 真の次元数を指定 M_truth <- 4 # ノイズ成分の標準偏差を指定 sigma <- 1.5 # ノイズ成分の精度を計算 lambda <- 1 / sigma^2 # パラメータを生成 w_truth_m <- sample(x = seq(-1, 1, by = 0.1), size = M_truth) # 作図用のx軸の値を作成 x_line <- seq(-3, 3, by = 0.01) # 観測モデルをデータフレームに格納 model_df <- tibble( x = x_line, y = t(w_truth_m) %*% t(x_vector(x_line, M_truth)) %>% as.vector() )
・観測データの生成
# (観測)データ数(試行回数)を指定 N <- 10 # 入力値を生成 x_n <- sample(seq(min(x_line), max(x_line), by = 0.01), size = N, replace = TRUE) # 出力値を計算:式(3.141) y_n <- (t(w_truth_m) %*% t(x_vector(x_n, M_truth)) + rnorm(n = N, mean = 0, sd = sqrt(1 / lambda))) %>% as.vector() y_n # 観測データをデータフレームに格納 x_df <- tibble( x = x_n, y = y_n )
・推論処理
# 事前分布の次元数の最大値(試行回数)を指定 M_max <- 15 # 受け皿を作成 predict_df <- tibble() # ベイズ推論 for(m in 1:M_max) { # 事前分布のパラメータをm次元に初期化 m_m <- rep(0, m) lambda_mm <- diag(m) * 0.01 # 入力値をm次元に拡張 x_nm <- x_vector(x_n, m) # 作図用のxをm次元に拡張 x_mat <- x_vector(x_line, m) # 事後分布のパラメータを計算:式(3.148) lambda_hat_mm <- lambda * t(x_nm) %*% x_nm + lambda_mm m_hat_m <- solve(lambda_hat_mm) %*% (lambda * t(t(y_n) %*% x_nm) + lambda_mm %*% matrix(m_m)) %>% as.vector() # 予測分布のパラメータを計算:式(3.155) mu_star_hat_line <- t(m_hat_m) %*% t(x_mat) %>% as.vector() sigma2_star_hat_line <- 1 / lambda + x_mat %*% solve(lambda_hat_mm) %*% t(x_mat) %>% diag() # 予測分布を計算 tmp_predict_df <- tibble( x = x_line, E_y = mu_star_hat_line, minus_sigma_y = E_y - sqrt(sigma2_star_hat_line), plus_sigma_y = E_y + sqrt(sigma2_star_hat_line), label = as.factor( paste0( "N=", N, ", M=", m, ", m_hat=(", paste(round(m_hat_m, 2), collapse = ", "), ")" ) ) ) # 結果を結合 predict_df <- rbind(predict_df, tmp_predict_df) }
こちらはM_maxをイタレーション数として、推論処理をM_max回繰り返します。
イタレーションごとに事前分布の次元数(パラメータ数)が変わるので、事前分布のパラメータを初期化します。また、x_nとx_lineの次元数も更新します。
・作図処理
# 予測分布を作図 predict_graph <- ggplot() + geom_line(data = predict_df, aes(x = x, y = E_y), color = "orange") + # 予測分布の期待値 geom_ribbon(data = predict_df, aes(x = x, ymin = minus_sigma_y, ymax = plus_sigma_y), fill = "#00A968", alpha = 0.1, color = "#00A968", linetype = "dotted") + # 予測分布の標準偏差 geom_line(data = model_df, aes(x = x, y = y), color = "blue", linetype = "dashed") + # 真のモデル geom_point(data = x_df, aes(x = x, y = y)) + # 観測データ gganimate::transition_manual(label) + # フレーム ylim(min(model_df[["y"]]) - 3 * sigma, max(model_df[["y"]]) + 3 * sigma) + # y軸の表示範囲 labs(title = "Predictive Distribution", subtitle = "{current_frame}", y = "y") # gif画像を作成 gganimate::animate(predict_graph, nframes = M_max, fps = 10)
$M$が大きくなるにしたがって、観測データに過剰適合していくのが確認できます。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
実装すると数式の間違いや理解が足りていなかったところに気付けてとても良いです。まぁ理解を深める工程の分先へは進めませんがね。
そういえばこれまでのモデルでは、予測分布のパラメータの計算に入力$\mathbf{x}_{*}$自体を使わなかったので気にもならなかったのですが、1つのデータに対して予測するので$\sum_{n=1}^N$の計算が不要みたいな流れでしたが、じゃあ複数のデータ(x_mline)の場合はどうなるんだ？結果の分布(グラフ)がちゃんとそれっぽくなってるんでこのままでいいんだろうなとして考えると、事後分布では観測した全てのデータの影響を考慮するため$N$個のデータとして扱う必要があるが、予測分布においては未知のデータ間の影響を考慮しようがないため独立に計算すればいいのだろうと理解しました。いかがでしょうか？
勉強をするほど疑問が出てくるのは、(うまく回っているうちは)良いことなんだと思ってます。
そんなただでさえ思ってるより進捗が出せない状況なのですが、この内容をPythonでも実装してみることにしました。
最後に昨日公開されたモーニング娘。'20さんの新曲です！ぜひBGMに♪
エビエビエビデンス🦐
【次節の内容】
- 2021/04/19：加筆修正しました。
今回の例だと、-3から3までを0.01刻みで新規の入力$x_{*}$として、そのそれぞれの出力$y_{*}$を分布で求めている。つまり、601個のy軸方向に広がる1次元ガウス分布を並べたのが記事中の予測分布の図。(自己解決)
