はじめに
『パターン認識と機械学習』の独学時のまとめです。一連の記事は「数式の行間埋め」または「R・Pythonでのスクラッチ実装」からアルゴリズムの理解を補助することを目的としています。本とあわせて読んでください。
この記事は、4.5.1項の内容です。ベイズロジスティック回帰のMAP推定と近似事後分布をR言語でします。
【数式読解編】
【前節の内容】
【他の節一覧】
【この節の内容】
・ベイズロジスティック回帰の実装:近似分布の導出とMAP推定
ベイズロジスティック回帰では、ロジスティック回帰に事前分布を導入します。しかし、ロジスティック回帰の事後分布は解析的に求められないため、ラプラス近似による近似事後分布を求めます。また、事後確率が最大となるMAP解を求めます。ニュートン-ラフソン法については4.3.3項、ラプラス近似については4.4.0項を参照してください。
利用するパッケージを読み込みます。
# 4.5.1項で利用するパッケージ library(tidyverse) library(mvnfast) library(gganimate)
推定過程をアニメーションで確認するのにgganimateパッケージを利用します。不要であれば省略してください。
・利用する関数の作成とデータの生成
基底関数の作成と、データの生成を行います。処理の流れは「【R】4.3.2-3：ロジスティック回帰の実装【PRMLのノート】 - からっぽのしょこ」と同じなので、詳しくはそちらを参照してください。
シグモイド基底関数の計画行列とクラスの境界線を引くための関数を作成します。
# ロジスティックシグモイド関数を作成 sigmoid <- function(x) { # ロジスティックシグモイド関数の計算 y <- 1 / (1 + exp(-x)) return(y) } # シグモイド基底関数の計画行列を作成 phi <- function(x_n) { # データを標準化 x_tilde_n <- (x_n - mean(x_n)) / sd(x_n) # マトリクスを初期化 y_n <- matrix(1, nrow = length(x_n), ncol = 2) # ロジスティックシグモイド関数による変換 y_n[, 2] <- sigmoid(x_tilde_n) return(y_n) } # 閾値の計算関数を作成:(シグモイド基底関数用) threshold <- function(w_m, x_vals) { # 回帰式の逆関数 x_tilde <- log(- w_m[1] / sum(w_m)) # 標準化の逆関数 x <- x_tilde * sd(x_vals) + mean(x_vals) return(x) }
シグモイド基底関数以外を使う場合は「【R】3.1.1：線形基底関数モデルの実装【PRMLのノート】 - からっぽのしょこ」も参照してください。
データ生成用のパラメータを指定して、人工的にデータを生成します。
# クラス割り当てのパラメータを指定:(クラス1となる確率) pi <- 0.4 # データ生成のパラメータを指定:(クラス0, クラス1) mu_k <- c(-10, 10) sigma_k <- c(5, 5) # データ数を指定 N <- 100 # 真のクラスを生成 t_n <- rbinom(n = N, size = 1, prob = pi) # データを生成 x_n <- rnorm(n = N, mean = mu_k[t_n+1], sd = sigma_k[t_n+1]) # 基底関数による変換 phi_x_nm <- phi(x_n) # 観測データをデータフレームに格納 data_df <- tidyr::tibble( x = x_n, t = t_n, class = as.factor(t_n) ) # 確認 head(data_df)
## # A tibble: 6 x 3 ## x t class ## <dbl> <int> <fct> ## 1 3.18 1 1 ## 2 7.01 1 1 ## 3 -12.2 0 0 ## 4 11.1 1 1 ## 5 -19.4 0 0 ## 6 5.33 1 1
生成した$N$個の観測データ(入力値$\mathbf{x} = \{x_1, \cdots, x_N\}$と目標値(各データのクラス)$\mathbf{t} = \{t_1, \cdots, t_N\}$)は、作図用にデータフレームに格納しておきます。
また、基底関数により非線形な変換を行った入力$\mathbf{x}$を$\boldsymbol{\Phi}$で表します。
・MAP推定と事後分布の近似
まずは、ニュートン-ラフソン法によりベイズロジスティック回帰のパラメータのMAP解を求めます。次に、MAP解を用いてラプラス近似による近似事後分布を求めます。これを指定した回数(収束するまで)繰り返します。
・処理の解説
まずは、繰り返し行う個々の処理を解説します。次に、それらをまとめて推定を行います。
モデルの重みパラメータと事前分布のパラメータの初期値を設定します。
# パラメータの次元数(基底関数の数)を設定:(固定) M <- 2 # 事前分布のパラメータの初期値を指定 m_m <- c(0, 0) # 平均 s_mm <- matrix(c(10, 0, 0, 10), nrow = M, ncol = M) # 分散共分散行列 s_inv_mm <- solve(s_mm) # 精度行列 # モデルのパラメータを初期化 w_logistic_m <- runif(n = M, min = -1, max = 1) # 確認 w_logistic_m
## [1] -0.4830551 0.9565946
この例では、ロジスティック回帰の重みパラメータ$\mathbf{w} = (w_0, w_1)^{\top}$の初期値を一様乱数で設定します。
$\mathbf{w}$の事前分布$p(\mathbf{w})$を、平均$\mathbf{m}_0$・分散共分散行列$\mathbf{S}_0$の2次元ガウス分布$\mathcal{N}(\mathbf{w} | \mathbf{m}_0, \mathbf{S}_0)$とします。超パラメータ$\mathbf{m}_0, \mathbf{S}_0$の下付き文字の0は、事前分布(初期値)を表す記号です。
変換した入力値の重み付き和を計算します。
# 重み付き和を計算 a_n <- (phi_x_nm %*% w_logistic_m) %>% as.vector() # 確認 w_logistic_m a_n[1:5]
## [1] -0.4830551 0.9565946 ## [1] 0.09365256 0.16535238 -0.20508344 0.23307428 -0.30894546
変換した入力値とロジスティック回帰のパラメータ$\mathbf{w}$の内積を計算します。
重み付き和をロジスティックシグモイド関数で変換します。
# ロジスティックシグモイド関数による変換 y_n <- sigmoid(a_n) # 確認 y_n[1:5]
## [1] 0.5233960 0.5412442 0.4489081 0.5580062 0.4233722
始めに作成したsigmoid()で計算します。
出力$\mathbf{y} = \{y_1, \cdots, y_N\}$の各項は、$0 < y_n < 1$の値をとり、対応する入力$x_n$のクラスが1である確率の推定値を表します。
重みパラメータを更新します。
# 中間変数を作成 r_nn <- diag(y_n) # 負の対数事後分布の勾配を計算 nabla_E_m <- (t(phi_x_nm) %*% (y_n - t_n) + s_inv_mm %*% (w_logistic_m - m_m)) %>% as.vector() # 負の対数事後分布のヘッセ行列を計算 h_mm <- s_inv_mm + t(phi_x_nm) %*% r_nn %*% phi_x_nm # 重みパラメータを更新 w_logistic_m <- (w_logistic_m - solve(h_mm) %*% nabla_E_m) %>% as.vector() # 確認 w_logistic_m
## [1] -2.448822 4.514220
次の式で$\mathbf{w}$を更新します。
収束した値(この例ではmax_iter_newtonに指定した回数更新した値)をMAP解$\mathbf{w}_{\mathrm{MAP}}$とします。
続いて、近似事後分布のパラメータを計算します。
# 事前分布のパラメータを更新 m_m <- w_logistic_m s_inv_mm <- h_mm s_mm <- solve(h_mm) # 確認 s_inv_mm s_mm
## [,1] [,2] ## [1,] 49.93012 25.96371 ## [2,] 25.96371 16.05245 ## [,1] [,2] ## [1,] 0.1260138 -0.2038186 ## [2,] -0.2038186 0.3919581
近似事後分布は、平均$\mathbf{m}_N$・分散共分散行列$\mathbf{S}_N$の2次元ガウス分布$\mathcal{N}(\mathbf{w} | \mathbf{m}_N, \mathbf{S}_N)$となります。超パラメータ$\mathbf{m}_N, \mathbf{S}_N$の下付き文字の$N$は、事後分布($N$個のデータから学習した値)を表す記号です。
パラメータは次の式で計算します。
最後に、負の対数事後確率密度を計算して、学習の推移を確認します。
# 負の対数事後確率を計算 y_n <- sigmoid(as.vector(phi_x_nm %*% w_logistic_m)) negative_log_likelihood <- -sum(dbinom(x = t_n, size = 1, prob = y_n, log = TRUE)) # 負の対数尤度関数 negative_log_prior <- -mvnfast::dmvn(X = w_logistic_m, mu = m_m, sigma = s_mm, log = TRUE) # 負の対数事前確率 # 確認 negative_log_likelihood negative_log_prior
## [1] 34.08891 ## [1] -0.5857354
次の式で計算します。ただし、$\mathbf{w}$に影響しない定数$\mathrm{const.}$は省略します。
以上が1回の試行で行う処理(計算)です。
次の試行では、推定した近似事後分布を事前分布として用いて、新たなMAP解と近似事後分布を求めます。この処理を収束するまで(この例ではmax_iter_laplaceに指定した回数)繰り返します。
・全体のコード
ニュートン-ラフソン法によりパラメータのMAP解を計算して、ラプラス近似による近似事後分布を計算します。
# パラメータの次元数(基底関数の数)を設定:(固定) M <- 2 # 事前分布のパラメータの初期値を指定 m_m <- c(0, 0) # 平均 s_mm <- matrix(c(10, 0, 0, 10), nrow = M, ncol = M) # 分散共分散行列 s_inv_mm <- solve(s_mm) # 精度行列 # モデルのパラメータを初期化 w_logistic_m <- runif(n = M, min = -1, max = 1) # 試行回数を指定 max_iter_laplace <- 500 max_iter_newton <- 100 # 推移の確認用の受け皿を作成 trace_w_mat <- matrix(NA, nrow = max_iter_laplace, ncol = M) # モデルの重みパラメータ trace_s_arr <- array(NA, dim = c(M, M, max_iter_laplace)) # 事前分布の分散共分散パラメータ trace_E_vec <- rep(NA, times = max_iter_laplace) # 負の対数事後確率 # ラプラス近似 for(i in 1:max_iter_laplace) { # ニュートン-ラフソン法による推定 for(j in 1:max_iter_newton) { # 重み付き和を計算 a_n <- (phi_x_nm %*% w_logistic_m) %>% as.vector() # ロジスティックシグモイド関数による変換 y_n <- sigmoid(a_n) # 中間変数を作成 r_nn <- diag(y_n) # 負の対数事後分布の勾配を計算 nabla_E_m <- (t(phi_x_nm) %*% (y_n - t_n) + s_inv_mm %*% (w_logistic_m - m_m)) %>% as.vector() # 負の対数事後分布のヘッセ行列を計算 h_mm <- s_inv_mm + t(phi_x_nm) %*% r_nn %*% phi_x_nm # 重みパラメータを更新 w_logistic_m <- (w_logistic_m - solve(h_mm) %*% nabla_E_m) %>% as.vector() } # 事前分布のパラメータを更新 m_m <- w_logistic_m s_inv_mm <- h_mm s_mm <- solve(h_mm) # 負の対数事後確率を計算 y_n <- sigmoid(as.vector(phi_x_nm %*% w_logistic_m)) negative_log_likelihood <- -sum(dbinom(x = t_n, size = 1, prob = y_n, log = TRUE)) # 負の対数尤度関数 negative_log_prior <- -mvnfast::dmvn(X = w_logistic_m, mu = m_m, sigma = s_mm, log = TRUE) # 負の対数事前確率 # 値を記録 trace_w_mat[i, ] <- w_logistic_m # モデルの重みパラメータ trace_s_arr[, , i] <- s_mm # 事前分布の分散共分散パラメータ trace_E_vec[i] <- negative_log_likelihood + negative_log_prior # 負の対数事後確率 # 途中経過を表示 #message("\r", "iter:", i, " (", round(i / max_iter_laplace * 100, 2), "%)", ", E=", trace_E_vec[i], appendLF = FALSE) }
試行回数max_iter_newton, max_iter_laplaceを指定して、パラメータw_logistic_m, m_m, s_mmの値を繰り返し更新します。
また、推移の確認用にパラメータをtrace_w_mat, trace_s_arrに、負の対数事後確率密度をtrace_E_vecに格納していきます。
・推定結果の可視化
推定結果と収束過程をグラフで確認します。
推定したパラメータを使って回帰曲線を計算します。
# 作図用のxの点を作成 x_vals <- seq(min(x_n) - 10, max(x_n) + 10, length.out = 1000) # MAP解によるモデルを計算 logistic_df <- tidyr::tibble( x = x_vals, # 入力値 a = phi(x_vals) %*% w_logistic_m %>% as.vector(), # 変換した入力の重み付き和 y = sigmoid(a) # 推定した確率値 ) # 確認 head(logistic_df)
## # A tibble: 6 x 3 ## x a y ## <dbl> <dbl> <dbl> ## 1 -36.6 -7.50 0.000554 ## 2 -36.5 -7.49 0.000559 ## 3 -36.5 -7.48 0.000563 ## 4 -36.4 -7.47 0.000568 ## 5 -36.3 -7.47 0.000572 ## 6 -36.2 -7.46 0.000577
未知の入力値(x軸の値)$x_{*}$に対する予測値(y軸の値)$y_{*}$を次の式で計算します。
x軸の値を作成してx_valsとします。x_valsの各要素に対して上の計算を行うことで、x軸の値ごとにクラス1となる確率(y軸の値)を求められます。
回帰曲線を作図します。
# MAP解によるモデルを作図 ggplot() + geom_point(data = data_df, aes(x = x, y = t, color = class)) + # 観測データ geom_line(data = logistic_df, aes(x = x, y = y)) + # 推定したモデル geom_vline(xintercept = threshold(w_logistic_m, x_vals), color = "#00A968", linetype = "dashed") + # クラスの閾値 labs(title = "Bayesian Logistic Regression", subtitle = paste0("iter:", max_iter_laplace, ", N=", N, ", threshold=", round(threshold(w_logistic_m, x_vals), 2), ", w=(", paste0(round(w_logistic_m, 2), collapse = ", "), ")"))
未知の入力$x_{*}$に対して、$y_{*} > 0.5$であればクラス1($t_{*} = 1$)、$y_{*} < 0.5$であればクラス0($t_{*} = 0$)に分類します。クラスの境界線($y = 0.5$となる$x$の値)は、始めに作成した関数threshold()で計算します。
事前分布を導入したことで、$\mathbf{w}$に対して制約を与えたことになります。そのため、4.3.2-3項のロジスティック回帰の結果よりもデータに対して緩やかにフィットしているのが分かります。
回帰曲線の推移をアニメーションで確認します。
・作図コード(クリックで展開)
各試行におけるモデルを計算します。
# フレームに利用する試行番号を指定 i_vec <- seq(1, max_iter_laplace, by = 1) #length(i_vec) # フレーム数 # 試行ごとのモデルを計算 trace_logistic_df <- tidyr::tibble() for(i in i_vec) { # i回目の値を取得 tmp_w_m <- trace_w_mat[i, ] # i回目のパラメータによるモデルを計算 tmp_logistic_df <- tidyr::tibble( x = x_vals, y = sigmoid(as.vector(phi(x_vals) %*% tmp_w_m)), threshold_x = threshold(tmp_w_m, x_vals), # y=0.5となるxの値 label = paste0( "iter:", i, ", E=", round(trace_E_vec[i], 2), ", threshold=", round(threshold(tmp_w_m, x_vals), 2), ", w=(", paste0(round(tmp_w_m, 2), collapse = ", "), ")" ) %>% as.factor() # フレーム切替用のラベル ) # 結果を結合 trace_logistic_df <- rbind(trace_logistic_df, tmp_logistic_df) # 途中経過を表示 #message("\r", "iter:", i, " (", round(i / max_iter_laplace * 100, 2), "%)", appendLF = "FALSE") } # 確認 head(trace_logistic_df)
## # A tibble: 6 x 4 ## x y threshold_x label ## <dbl> <dbl> <dbl> <fct> ## 1 -36.6 0.0104 0.919 iter:1, E=11.34, threshold=0.92, w=(-6.35, 11.93) ## 2 -36.5 0.0105 0.919 iter:1, E=11.34, threshold=0.92, w=(-6.35, 11.93) ## 3 -36.5 0.0106 0.919 iter:1, E=11.34, threshold=0.92, w=(-6.35, 11.93) ## 4 -36.4 0.0106 0.919 iter:1, E=11.34, threshold=0.92, w=(-6.35, 11.93) ## 5 -36.3 0.0107 0.919 iter:1, E=11.34, threshold=0.92, w=(-6.35, 11.93) ## 6 -36.2 0.0107 0.919 iter:1, E=11.34, threshold=0.92, w=(-6.35, 11.93)
各試行で求めたパラメータの値をtrace_w_matから取り出して回帰曲線を計算し、計算結果をデータフレームに追加していきます。
データへの当てはまり具合を確認するため、負の対数事後確率密度の値もtrace_E_vecから取り出します。
作図してgif画像として出力ます。
# MAP解によるモデルを作図 anime_model_graph <- ggplot() + geom_point(data = data_df, aes(x = x, y = t, color = class)) + # 観測データ geom_line(data = trace_logistic_df, aes(x = x, y = y)) + # 推定したモデル geom_vline(data = trace_logistic_df, aes(xintercept = threshold_x), color = "#00A968", linetype = "dashed") + # クラスの閾値 gganimate::transition_manual(label) + # フレーム xlim(c(min(x_vals), max(x_vals))) + # x軸の表示範囲 labs(title = "Bayesian Logistic Regression", subtitle = "{current_frame}") # gif画像に変換 gganimate::animate(anime_model_graph, nframe = length(i_vec), fps = 10)
徐々にデータに適合していくのを確認できます。
ちなみに、4.3.2-3項のロジスティック回帰では、max_iter_newton回の試行をそれぞれ1フレームとしました。ここでは、max_iter_newton解分の結果を1フレームとしています。
作図用の$\mathbf{w}$の点を作成して、近似事後分布を計算します。
# 作図用のwの点を作成 w_im <- expand.grid( seq(m_m[1] - 3 * sqrt(s_mm[1, 1]), m_m[1] + 3 * sqrt(s_mm[1, 1]), length.out = 500), # x軸の値 seq(m_m[2] - 3 * sqrt(s_mm[2, 2]), m_m[2] + 3 * sqrt(s_mm[2, 2]), length.out = 500) # y軸の値 ) %>% as.matrix() # 重みパラメータの近似事後分布を計算 posterior_df <- tidyr::tibble( w_0 = w_im[, 1], # x軸の値 w_1 = w_im[, 2], # y軸の値 density = mvnfast::dmvn(X = w_im, mu = m_m, sigma = s_mm) # 確率密度 ) # 確認 head(posterior_df)
## # A tibble: 6 x 3 ## w_0 w_1 density ## <dbl> <dbl> <dbl> ## 1 -10.3 17.6 0 ## 2 -10.3 17.6 0 ## 3 -10.3 17.6 0 ## 4 -10.3 17.6 0 ## 5 -10.3 17.6 0 ## 6 -10.3 17.6 0
$w_0, w_1$がとり得る値を作成して、全ての値の組み合わせとなるように$\mathbf{w}$の点を作成します。それぞれ平均$m_0, m_1$から標準偏差$\sqrt{s_{0,0}}, \sqrt{s_{1,1}}$の3倍を範囲とします。
多変量ガウス分布の確率密度は、mvnfastパッケージのdmvn()で計算できます。平均の引数muにm_mまたはw_logistic_mを、分散共分散行列の引数sigmaにs_mmを指定します。
近似事後分布のグラフを作成します。
# 重みパラメータのMAP解を格納 w_df <- tidyr::tibble( w_0 = w_logistic_m[1], # x軸の値 w_1 = w_logistic_m[2] # y軸の値 ) # 重みパラメータの近似事後分布を作図 ggplot() + geom_contour(data = posterior_df, aes(x = w_0, y = w_1, z = density, color = ..level..)) + # 近似事後分布 geom_point(data = w_df, aes(x = w_0, y = w_1), color = "red", shape = 4, size = 5) + # MAP解 labs(title = "Approximate Posterior Distribution", subtitle = paste0("iter:", max_iter_laplace, ", w=(", paste0(round(w_logistic_m, 2), collapse = ", "), ")", ", s=(", paste0(round(s_mm, 5), collapse = ", "), ")"), x = expression(w[0]), y = expression(w[1]), color = "density")
各軸の範囲を見ると分かるように、先の尖った分布になっています。
近似事後分布の推移もアニメーションで確認します。
・作図コード(クリックで展開)
# フレームに利用する試行回数を指定 i_vec <- seq(250, max_iter_laplace, by = 10) length(i_vec) # フレーム数 # 試行ごとに近似事後分布を計算 trace_w_df <- tidyr::tibble() trace_posterior_df <- tidyr::tibble() for(i in i_vec) { # i回目のパラメータを取得 tmp_w_m <- trace_w_mat[i, ] tmp_m_m <- trace_w_mat[i, ] tmp_s_mm <- trace_s_arr[, , i] # フレーム切替用のラベルを作成 label_text <- paste0( "iter:", i, ", E=", round(trace_E_vec[i], 2), ", w=(", paste0(round(tmp_w_m, 2), collapse = ", "), ")", ", s=(", paste0(round(tmp_s_mm, 5), collapse = ", "), ")" ) # i回目の重みパラメータを格納 tmp_w_df <- tidyr::tibble( w_0 = tmp_w_m[1], # x軸の値 w_1 = tmp_w_m[2], # y軸の値 label = as.factor(label_text) # フレーム切替用のラベル ) # 作図用のwの点を作成 w_im <- expand.grid( seq(tmp_m_m[1] - 3 * sqrt(tmp_s_mm[1, 1]), tmp_m_m[1] + 3 * sqrt(tmp_s_mm[1, 1]), length.out = 500), # x軸の値 seq(tmp_m_m[2] - 3 * sqrt(tmp_s_mm[2, 2]), tmp_m_m[2] + 3 * sqrt(tmp_s_mm[2, 2]), length.out = 500) # y軸の値 ) %>% as.matrix() # i回目の近似事後分布を計算 tmp_posterior_df <- tidyr::tibble( w_0 = w_im[, 1], # x軸の値 w_1 = w_im[, 2], # y軸の値 density = mvnfast::dmvn(X = w_im, mu = tmp_m_m, sigma = tmp_s_mm), # 確率密度 label = as.factor(label_text) # フレーム切替用のラベル ) # 結果を結合 trace_w_df <- rbind(trace_w_df, tmp_w_df) trace_posterior_df <- rbind(trace_posterior_df, tmp_posterior_df) # 途中経過を表示 message("\r", "iter:", i, " (", round(i / max_iter_laplace * 100, 2), "%)", appendLF = "FALSE") } # 重みパラメータの近似事後分布を作図 anime_posterior_graph <- ggplot() + geom_contour(data = trace_posterior_df, aes(x = w_0, y = w_1, z = density, color = ..level..)) + # 近似事後分布 geom_point(data = trace_w_df, aes(x = w_0, y = w_1), color = "red", shape = 4, size = 5) + # MAP解 gganimate::transition_manual(label) + # フレーム labs(title = "Approximate Posterior Distribution", subtitle = "{current_frame}", x = expression(w[0]), y = expression(w[1])) # gif画像に変換 gganimate::animate(anime_posterior_graph, nframe = length(i_vec), fps = 10)
左の図は、1から25回目の結果を1試行ずつ表示しています。右の図は、250から500(最後)回目の結果を10試行ずつ表示しています。
平均(MAP解)の周りの等高線の幅が狭くなっていく(分布が尖っていく)のが分かります。始めの頃の分布は、確率密度の値が小さいため描画されないことがあります。
負の対数事後確率密度の推移を作図します。
# 負の対数事後確率の推移を格納 E_df <- tidyr::tibble( iteration = 1:max_iter_laplace, E = trace_E_vec ) # 負の対数事後確率の推移を作図 ggplot(E_df, aes(x = iteration, y = E)) + geom_line() + labs(title = "Negative log Posterior", y = expression(E(w)))
試行を繰り返す度に値が小さくなっているのが分かります。
ロジスティック回帰の重みパラメータ(近似事後分布の平均パラメータ)と近似事後分布の精度パラメータの推移を作図します。
# モデルの重みパラメータの推移を格納 trace_w_df <- trace_w_mat %>% tidyr::as_tibble() %>% # データフレームに変換 dplyr::rename_with(.fn = ~0:(M-1), .cols = 1:M) %>% # 列名を付与 dplyr::mutate(iteration = 1:max_iter_laplace) %>% # 試行番号列を結合 tidyr::pivot_longer(cols = !iteration, names_to = "m", values_to = "value") # 縦型に変換 # モデルの重みパラメータの推移を作図 ggplot(trace_w_df, aes(x = iteration, y = value, color = m)) + geom_line() + labs(title = "w = m", y = expression(w[m])) # 事前分布の分散共分散パラメータの推移を格納 trace_s_df <- trace_s_arr %>% as.vector() %>% # ベクトルに変換 matrix(nrow = max_iter_laplace, ncol = M^2, byrow = TRUE) %>% # マトリクスに変換 dplyr::as_tibble() %>% # データフレームに変換 dplyr::rename_with(.fn = ~paste0("s_", rep(1:M, times = M), rep(1:M, each = M))) %>% # 列名を付与 dplyr::mutate(iteration = 1:max_iter_laplace) %>% # 試行番号列を結合 tidyr::pivot_longer(cols = !iteration, names_to = "idx", values_to = "value") # 縦型に変換 # 事前分布の分散共分散パラメータの推移を作図 ggplot(trace_s_df, aes(x = iteration, y = value, color = idx)) + geom_line() + #ylim(c(-0.5, 0.5)) + # y軸の表示範囲 labs(title = "S", y = expression(s[mm]))
各要素の値が収束しているのが分かります。
最後の近似事後分布の負の対数のグラフに、重みパラメータの推移を重ねて確認します。
# 作図用のwの点を作成 w_im <- expand.grid( seq(min(trace_w_mat[, 1]) - 5, max(trace_w_mat[, 1]) + 5, length.out = 500), seq(min(trace_w_mat[, 2]) - 5, max(trace_w_mat[, 2]) + 5, length.out = 500) ) %>% as.matrix() # 出力を計算 y_ni <- sigmoid(phi_x_nm %*% t(w_im)) # 負の対数事後確率を計算 negative_log_likelihood_i <- -colSums(log(y_ni^t_n) + log((1 - y_ni)^(1 - t_n))) # 負の対数尤度関数 negative_log_prior_i <- -mvnfast::dmvn(X = w_im, mu = m_m, sigma = s_mm, log = TRUE) # 負の対数事前確率 # 負の対数事後確率を格納 error_df <- tidyr::tibble( w_0 = w_im[, 1], # x軸の値 w_1 = w_im[, 2], # y軸の値 E = negative_log_likelihood_i + negative_log_prior_i ) # パラメータの推移をデータフレームに格納 trace_w_df <- trace_w_mat %>% dplyr::as_tibble() %>% # データフレームに変換 dplyr::rename_with(.fn = ~paste0("w_", 0:1), .cols = 1:M) %>% # 列名を付与 dplyr::mutate(iteration = 1:max_iter_laplace) # 試行回数を結合 # パラメータの推移を作図 ggplot() + geom_contour(data = error_df, aes(x = w_0, y = w_1, z = E, color = ..level..)) + # 誤差関数 #geom_contour(data = error_df, aes(x = w_0, y = w_1, z = E, color = ..level..), # breaks = seq(0, max(error_df[["E"]]) * 0.1, length.out = 10)) + # 誤差関数:(等高線を引く値を指定) geom_point(data = trace_w_df, aes(x = w_0, y = w_1), color = "#00A968") + # パラメータの推移 labs(title = "Negative Log Posterior", subtitle = paste0("iter:", max_iter_laplace), x = expression(w[0]), y = expression(w[1]), color = expression(E(w)))
(もう少しうまく可視化できると思ったのだが、一応)勾配に従ってパラメータが更新されているのが分かります。ただし、試行の度に誤差関数の勾配は変化します。
参考文献
- C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上下』,丸善出版,2012年.
おわりに
これはMAP推定ということでいいんですよね？本にそういう表現はないけど。あと、負の対数事後確率って呼んでいいんですか？誤差関数と呼んでもいいんですか？
4章もクライマックスですね、あと少し。
2021年12月27日は、Juice=Juiceの稲場愛香さんの24歳のお誕生日です！
サムネの真ん中の方です。来年もがんばりまなかん🌸゜｡(´ω ` )／🌙
【次節の内容】
年内にできるといいが・・・
