はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」によって、「数式」と「プログラム」から理解するのが目標です。
省略してある内容等ありますので、本とあわせて読んでください。
この記事では、「尤度関数を精度が未知の多次元ガウス分布(多変量正規分布)」、「事前分布をウィシャート分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」の計算をR言語で実装します。
【数式読解編】
【前の節の内容】
【他の節の内容】
【この節の内容】
3.4.2 多次元ガウス分布の学習と予測：精度が未知の場合
尤度関数を多次元ガウス分布(Multivariate Gaussian Distribution)・多変量正規分布(Multivariate Normal Distribution)、事前分布をウィシャート分布(Wishart Distribution)とするモデルに対するベイズ推論を実装します。人工的に生成したデータを用いて、ガウス分布の精度パラメータを推定し、また未観測データに対する予測分布を求めます。
多次元ガウス分布については「多次元ガウス分布の定義式 - からっぽのしょこ」、ウィシャート分布については「ウィシャート分布の定義式 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用パッケージ library(tidyverse) library(mvnfast) library(gganimate)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。
magrittrパッケージのパイプ演算子%>%ではなく、ネイティブパイプ(ベースパイプ)演算子|>を使っています。%>%に置き換えても処理できます。
学習の推移をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。
ベイズ推論の実装
まずは、モデルを設定して、データを生成します。生成したデータを用いて、事後分布のパラメータを計算します。さらに、事後分布のパラメータを用いて、予測分布のパラメータを計算します。
生成分布の設定
データの生成分布(真の分布)として、多次元ガウス分布$\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})$を設定します。
ガウス分布のグラフ作成については「【R】2次元ガウス分布の作図 - からっぽのしょこ」、ウィシャート分布の可視化については「【R】ウィシャート分布の可視化 - からっぽのしょこ」を参照してください。
真の分布($D$次元ガウス分布)のパラメータ$\boldsymbol{\mu}, \boldsymbol{\Sigma} = \boldsymbol{\Lambda}^{-1}$を設定します。この例では2次元のグラフで表現するため、次元数を$D = 2$とします。パラメータの計算自体は次元数に関わらず行えます。
# 次元数を指定 D <- 2 # (既知の)平均ベクトルを指定 mu_d <- c(25, 50) # 真の分散共分散行列を指定 sigma_truth_dd <- matrix(c(900, -100, -100, 400), nrow = D, ncol = D) # 真の精度行列を計算 lambda_truth_dd <- solve(sigma_truth_dd) lambda_truth_dd
## [,1] [,2] ## [1,] 0.0011428571 0.0002857143 ## [2,] 0.0002857143 0.0025714286
平均ベクトル$\boldsymbol{\mu} = (\mu_1, \mu_2)^{\top}$をmu_dとして値を指定します。
分散共分散行列$\boldsymbol{\Sigma} = (\sigma_1^2, \sigma_{2,1}, \sigma_{1,2}, \sigma_2^2)$をsigma_truth_ddとして値を指定します。$\sigma_d$は$x_d$の標準偏差、$\sigma_d^2$は$x_d$の分散、$\sigma_{i,j} = \sigma_{j,i}$は$x_i, x_j$の共分散であり、$\sigma_d^2$は正の実数、$\sigma_{i,j}$は実数、また$\boldsymbol{\Sigma}$は正定値行列を満たす必要があります。
精度行列$\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{-1}$を計算してlambda_truth_ddとします。逆行列はsolve()で計算できます。lambda_truth_ddが真の値であり、この値を求めるのがここでの目的です。
真の分布の確率変数がとり得る値$\mathbf{x}$を作成します。
# グラフ用のxの値を作成 x_1_vec <- seq( mu_d[1] - sqrt(sigma_truth_dd[1, 1]) * 3, mu_d[1] + sqrt(sigma_truth_dd[1, 1]) * 3, length.out = 201 ) x_2_vec <- seq( mu_d[2] - sqrt(sigma_truth_dd[2, 2]) * 3, mu_d[2] + sqrt(sigma_truth_dd[2, 2]) * 3, length.out = 201 ) # グラフ用のxの点を作成 x_mat <- tidyr::expand_grid( x_1 = x_1_vec, x_2 = x_2_vec ) |> # 格子点を作成 as.matrix() # マトリクスに変換 head(x_mat)
## x_1 x_2 ## [1,] -65 -10.0 ## [2,] -65 -9.4 ## [3,] -65 -8.8 ## [4,] -65 -8.2 ## [5,] -65 -7.6 ## [6,] -65 -7.0
$x_1$(x軸)の値をx_1_vals、$x_2$(y軸)の値をx_2_valsとします。この例では、それぞれ平均を中心に標準偏差の3倍を範囲とします。
x_1_valsとx_2_valsの要素の全ての組み合わせ(格子状の点)をexpand_grid()で作成します。データフレームが出力されるので、as.matrix()でマトリクスに変換してx_matとします。x_matの各行が点$\mathbf{x} = (x_1, x_2)^{\top}$に対応します。
真の分布を計算して、作図用のデータフレームを作成します。
# 真の分布(多次元ガウス分布)を計算:式(2.72) model_dens_df <- tibble::tibble( x_1 = x_mat[, 1], # x軸の値 x_2 = x_mat[, 2], # y軸の値 dens = mvnfast::dmvn(X = x_mat, mu = mu_d, sigma = sigma_truth_dd) # 確率密度 ) model_dens_df
## # A tibble: 40,401 × 3 ## x_1 x_2 dens ## <dbl> <dbl> <dbl> ## 1 -65 -10 0.00000000549 ## 2 -65 -9.4 0.00000000611 ## 3 -65 -8.8 0.00000000680 ## 4 -65 -8.2 0.00000000756 ## 5 -65 -7.6 0.00000000839 ## 6 -65 -7 0.00000000931 ## 7 -65 -6.4 0.0000000103 ## 8 -65 -5.8 0.0000000114 ## 9 -65 -5.2 0.0000000126 ## 10 -65 -4.6 0.0000000140 ## # … with 40,391 more rows
x_matの行ごとに確率密度を計算します。多次元ガウス分布の確率密度は、mvnfastパッケージのdmvn()で計算できます。確率変数の引数Xにx_mat、平均ベクトルの引数muにmu_d、分散共分散行列の引数sigmaにsigma_truth_ddまたはlambda_truth_ddの逆行列を指定します。
真の分布のグラフを作成します。
# パラメータラベル用の文字列を作成 model_param_text <- paste0( "list(mu==(list(", paste(round(mu_d, 1), collapse = ", "), "))", ", Lambda==(list(", paste(round(lambda_truth_dd, 5), collapse = ", "), ")))" ) # 真の分布を作図 ggplot() + geom_contour(data = model_dens_df, mapping = aes(x = x_1, y = x_2, z = dens, color = ..level..)) + # 真の分布:(等高線図) #geom_contour_filled(data = model_dens_df, mapping = aes(x = x_1, y = x_2, z = dens, fill = ..level..), alpha = 0.8) + # 真の分布:(塗りつぶし等高線図) labs(title = "Multivariate Gaussian Distribution", subtitle = parse(text = model_param_text), color = "density", fill = "density", x = expression(x[1]), y = expression(x[2]))
geom_contour()またはgeom_contour_filled()で等高線グラフを描画します。
真のパラメータを求めることは、真の分布を求めることを意味します。
これまでのように、事前分布や事後分布の確率密度を計算してグラフで確認したいところですが、精度行列の分布(ウィシャート分布)をグラフ化するのは難しいです。そこで、精度行列の期待値(ウィシャート分布の期待値)の逆行列(分散共分散行列の期待値)によるマハラノビス距離の等高線グラフと、楕円形になるマハラノビス距離の等高線の軸によって可視化することにします。
比較のため、真の分散共分散行列を用いてマハラノビス距離とその軸のグラフを作成します。
真の分散共分散行列$\boldsymbol{\Sigma}$を用いて、マハラノビス距離を計算します。
# 真のΛによるマハラノビス距離を計算 model_delta_df <- tibble::tibble( x_1 = x_mat[, 1], # x軸の値 x_2 = x_mat[, 2], # y軸の値 delta = mahalanobis(x = x_mat, center = mu_d, cov = sigma_truth_dd) |> sqrt() # 距離 ) model_delta_df
## # A tibble: 40,401 × 3 ## x_1 x_2 delta ## <dbl> <dbl> <dbl> ## 1 -65 -10 4.65 ## 2 -65 -9.4 4.62 ## 3 -65 -8.8 4.60 ## 4 -65 -8.2 4.58 ## 5 -65 -7.6 4.56 ## 6 -65 -7 4.53 ## 7 -65 -6.4 4.51 ## 8 -65 -5.8 4.49 ## 9 -65 -5.2 4.46 ## 10 -65 -4.6 4.44 ## # … with 40,391 more rows
マハラノビス距離の2乗をmahalanobis()で計算します。変数の引数Xにx_mat、平均ベクトルの引数centerにmu_d、分散共分散行列の引数covにsigma_truth_ddまたはlambda_truth_ddを指定します。精度行列を使う場合はinverted = TRUEを指定します。
計算結果の平方根をとりマハラノビス距離を計算します。平方根はsqrt()で計算できます。
真の分散共分散行列$\boldsymbol{\Sigma}$の固有値と固有ベクトルを用いて、楕円(マハラノビス距離の等高線)の軸を計算します。
# 真のΛの固有値・固有ベクトルを計算 model_eigen <- eigen(sigma_truth_dd) model_lmd_d <- model_eigen[["values"]] model_u_dd <- model_eigen[["vectors"]] |> t() # 真のΛによる楕円の軸を計算 model_axis_df <- tibble::tibble( xend = mu_d[1] + model_u_dd[, 1] * sqrt(model_lmd_d), # x軸の値 yend = mu_d[2] + model_u_dd[, 2] * sqrt(model_lmd_d) # y軸の値 ) model_axis_df
## # A tibble: 2 × 2 ## xend yend ## <dbl> <dbl> ## 1 -4.77 55.7 ## 2 21.3 30.8
固有値$\lambda_d$と固有ベクトル$\mathbf{u}_d = (u_{d,1}, u_{d,2})^{\top}$をeigen()で計算します。リストが出力されるので、"values"で固有値(をまとめたベクトル)、"vectors"で固有ベクトル(をまとめたマトリクス)を取り出します。数式での成分と合わせるために転置しておきます。
軸番号(長軸か短軸か)を$i$、次元(x軸かy軸か)を$j$として、軸の先端の点を$\mu_i + u_{j,i} \sqrt{\lambda_j}$で計算します。
真の分散共分散行列によるマハラノビス距離の等高線とその軸のグラフを作成します。
# 等高線を引く値を指定 break_vals <- seq(0, ceiling(max(model_delta_df[["delta"]])), by = 0.5) # 真の分散共分散行列による距離と軸を作図 ggplot() + geom_contour(data = model_delta_df, aes(x = x_1, y = x_2, z = delta, color = ..level..), breaks = break_vals) + # 真のΛによる距離:(等高線図) #geom_contour_filled(data = model_delta_df, aes(x = x_1, y = x_2, z = delta, fill = ..level..), # breaks = break_vals, alpha = 0.8) + # 真のΛによる距離:(塗りつぶし等高線図) geom_segment(data = model_axis_df, mapping = aes(x = mu_d[1], y = mu_d[2], xend = xend, yend = yend, linetype = "model"), color = "red", size = 1, arrow = arrow(length = unit(10, "pt"))) + # 真のΛによる軸 scale_linetype_manual(breaks = "model", values = "solid", labels = "true model", name = "axis") + # (凡例表示用の黒魔術) coord_fixed(ratio = 1) + # アスペクト比 labs(title = "Mahalanobis Distance", subtitle = parse(text = model_param_text), color = "distance", fill = "distance", x = expression(x[1]), y = expression(x[2]))
geom_segment()で軸(線分)を描画します。線分の始点の引数x, yにmu_dの値、終点の引数xend, yendにxend, yend列を指定します。
データの生成
続いて、構築したモデルに従って観測データ$\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N\}$、$\mathbf{x}_n = (x_{n,1}, x_{n,2}, \cdots, x_{n,D})^{\top}$を生成します。
多次元ガウス分布の乱数生成については「【R】多次元ガウス分布の乱数生成 - からっぽのしょこ」を参照してください。
多次元ガウス分布に従う$N$個のデータをランダムに生成します。
# (観測)データ数を指定 N <- 100 # 多次元ガウス分布に従うデータを生成 x_nd <- mvnfast::rmvn(n = N, mu = mu_d, sigma = sigma_truth_dd) head(x_nd)
## [,1] [,2] ## [1,] 38.99642 55.58055 ## [2,] -20.83073 57.59026 ## [3,] 42.53096 41.83315 ## [4,] 20.72358 14.22715 ## [5,] 61.18807 24.50534 ## [6,] 21.36816 76.59592
多次元ガウス分布に従う乱数は、mvnfastパッケージのrmvn()で生成できます。データ数の引数nにN、平均ベクトルの引数muにmu_d、分散共分散行列の引数sigmaにsigma_truth_ddまたはlambda_truth_ddの逆行列を指定します。生成したN個のデータをx_ndとします。
観測データの散布図を真の分布と重ねて確認します。
# 観測データを格納 data_df <- tibble::tibble( x_1 = x_nd[, 1], # x軸の値 x_2 = x_nd[, 2] # y軸の値 ) # パラメータラベル用の文字列を作成 sample_param_text <- paste0( "list(mu==(list(", paste(round(mu_d, 1), collapse = ", "), "))", ", Lambda==(list(", paste(round(lambda_truth_dd, 5), collapse = ", "), "))", ", N==", N, ")" ) # 観測データの散布図を作成 ggplot() + geom_contour(data = model_dens_df, aes(x = x_1, y = x_2, z = dens, color = ..level..)) + # 真の分布:(等高線図) #geom_contour_filled(data = model_dens_df, aes(x = x_1, y = x_2, z = dens, fill = ..level..), alpha = 0.8) + # 真の分布:(塗りつぶし等高線図) geom_point(data = data_df, aes(x = x_1, y = x_2), color = "orange") + # 観測データ labs(title = "Multivariate Gaussian Distribution", subtitle = parse(text = sample_param_text), color = "density", fill = "density", x = expression(x[1]), y = expression(x[2]))
事前分布の設定
次は、尤度に対する共役事前分布を設定します。多次元ガウス分布の精度パラメータ$\boldsymbol{\Lambda}$の事前分布$p(\boldsymbol{\Lambda})$としウィシャート分布$\mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W})$を設定します。
$\boldsymbol{\Lambda}$の事前分布(ウィシャート分布)のパラメータ(超パラメータ)$\nu, \mathbf{W}$を設定します。
# Λの事前分布の自由度を指定 nu <- D # Λの事前分布の逆スケール行列を指定 w_dd <- diag(D) * 0.0001
ウィシャート分布の自由度$\nu$をnu、逆スケール行列$\mathbf{W} = (w_{1,1}, w_{2,1}, w_{1,2}, w_{2,2})$をw_ddとして値を指定します。$\nu$は$\nu > D - 1$、$\mathbf{W}$は正定値行列を満たす必要があります。
$\boldsymbol{\Lambda}$の事前分布の期待値(精度行列の期待値)を計算します。
# Λの事前分布の期待値を計算:式(2.89) E_lambda_dd <- nu * w_dd E_lambda_dd
## [,1] [,2] ## [1,] 2e-04 0e+00 ## [2,] 0e+00 2e-04
ウィシャート分布の期待値を次の式で計算します。
精度行列の期待値$\mathbb{E}[\boldsymbol{\Lambda}]$を用いて、マハラノビス距離を計算します。
# Λの期待値によるマハラノビス距離を計算 prior_delta_df <- tibble::tibble( x_1 = x_mat[, 1], # x軸の値 x_2 = x_mat[, 2], # y軸の値 delta = mahalanobis(x = x_mat, center = mu_d, cov = E_lambda_dd, inverted = TRUE) |> sqrt() # 距離 ) prior_delta_df
## # A tibble: 40,401 × 3 ## x_1 x_2 delta ## <dbl> <dbl> <dbl> ## 1 -65 -10 1.53 ## 2 -65 -9.4 1.53 ## 3 -65 -8.8 1.52 ## 4 -65 -8.2 1.52 ## 5 -65 -7.6 1.51 ## 6 -65 -7 1.51 ## 7 -65 -6.4 1.50 ## 8 -65 -5.8 1.50 ## 9 -65 -5.2 1.49 ## 10 -65 -4.6 1.49 ## # … with 40,391 more rows
真の分布のときと同様にして処理します。
精度行列の期待値$\mathbb{E}[\boldsymbol{\Lambda}]$の逆行列を用いて、楕円の軸を計算します。
# Λの期待値の固有値・固有ベクトルを計算 prior_eigen <- eigen(solve(E_lambda_dd)) prior_lmd_d <- prior_eigen[["values"]] prior_u_dd <- prior_eigen[["vectors"]] |> t() # Λの期待値による楕円の軸を計算 prior_axis_df <- tibble::tibble( xend = mu_d[1] + prior_u_dd[, 1] * sqrt(prior_lmd_d), # x軸の値 yend = mu_d[2] + prior_u_dd[, 2] * sqrt(prior_lmd_d) # y軸の値 ) prior_axis_df
## # A tibble: 2 × 2 ## xend yend ## <dbl> <dbl> ## 1 25 121. ## 2 -45.7 50
真の分布のときと同様にして処理します。
$\boldsymbol{\Lambda}$の事前分布の期待値(精度行列の期待値)によるマハラノビス距離の等高線とその軸のグラフを作成します。
# パラメータラベル用の文字列を作成 prior_param_text <- paste0( "list(nu==", nu, ", W==(list(", paste(w_dd, collapse = ", "), ")))" ) # Λの事前分布の期待値による距離と軸を作図 ggplot() + geom_contour(data = model_delta_df, mapping = aes(x = x_1, y = x_2, z = delta, color = ..level..), alpha = 0.5, linetype = "dashed") + # 真のΛによる距離 geom_contour(data = prior_delta_df, mapping = aes(x = x_1, y = x_2, z = delta, color = ..level..)) + # Λの期待値による距離 geom_segment(data = model_axis_df, mapping = aes(x = mu_d[1], y = mu_d[2], xend = xend, yend = yend, linetype = "model"), color = "red", size = 1, arrow = arrow(length = unit(10, "pt"))) + # 真のΛによる軸 geom_segment(data = prior_axis_df, mapping = aes(x = mu_d[1], y = mu_d[2], xend = xend, yend = yend, linetype = "prior"), color = "blue", size = 1, arrow = arrow(length = unit(10, "pt"))) + # Λの期待値による軸 scale_linetype_manual(breaks = c("prior", "model"), values = c("solid", "dashed"), labels = c("prior", "true model"), name = "axis") + # (凡例表示用の黒魔術) guides(linetype = guide_legend(override.aes = list(color = c("blue", "red"), size = c(0.5, 0.5), linetype = c("solid", "dashed")))) + # (凡例表示用の黒魔術) coord_fixed(ratio = 1) + # アスペクト比 labs(title = "Mahalanobis Distance", subtitle = parse(text = prior_param_text), color = "distance", x = expression(x[1]), y = expression(x[2]))
真の精度行列による等高線と軸を破線で重ねて描画します。
事後分布の計算
観測データ$\mathbf{X}$から精度パラメータ$\boldsymbol{\Lambda}$の事後分布$p(\boldsymbol{\Lambda} | \mathbf{X}) = \mathcal{W}(\boldsymbol{\Lambda} | \hat{\nu}, \hat{\mathbf{W}})$を求めます(精度パラメータ$\boldsymbol{\Lambda}$を分布推定します)。
観測データ$\mathbf{X}$を用いて、$\boldsymbol{\Lambda}$の事後分布(ウィシャート分布)のパラメータ$\hat{\nu}, \hat{\mathbf{W}}$を計算します。
# Λの事後分布の自由度を計算:式(3.116) nu_hat <- N + nu # Λの事後分布の逆スケール行列を計算:式(3.116) w_hat_dd <- solve((t(x_nd) - mu_d) %*% t(t(x_nd) - mu_d) + solve(w_dd)) nu_hat; w_hat_dd
## [1] 102 ## [,1] [,2] ## [1,] 1.026742e-05 2.188074e-06 ## [2,] 2.188074e-06 2.545089e-05
$\boldsymbol{\Lambda}$の事後分布のパラメータを次の式で計算します。
ただし、$\sum_{n=1}^N (\mathbf{x}_n - \mu) (\mathbf{x}_n - \mu)^{\top}$の計算について、$\tilde{x}_{n,d} = x_{n,d} - \mu_d$として$\tilde{\mathbf{X}}^{\top} \tilde{\mathbf{X}}$で計算しています。
$\boldsymbol{\Lambda}$の事後分布の期待値(精度行列の期待値)$\mathbb{E}[\boldsymbol{\Lambda}]$を計算します。
# Λの事後分布の期待値を計算:式(2.89) E_lambda_hat_dd <- nu_hat * w_hat_dd E_lambda_hat_dd
## [,1] [,2] ## [1,] 0.0010472763 0.0002231836 ## [2,] 0.0002231836 0.0025959910
事前分布のときと同様に計算してE_lambda_hat_ddとします。
精度行列の期待値$\mathbb{E}[\boldsymbol{\Lambda}]$を用いて、マハラノビス距離を計算します。
# Λの期待値によるマハラノビス距離を計算 posterior_delta_df <- tibble::tibble( x_1 = x_mat[, 1], # x軸の値 x_2 = x_mat[, 2], # y軸の値 delta = mahalanobis(x = x_mat, center = mu_d, cov = E_lambda_hat_dd, inverted = TRUE) |> sqrt() # 距離 ) posterior_delta_df
## # A tibble: 40,401 × 3 ## x_1 x_2 delta ## <dbl> <dbl> <dbl> ## 1 -65 -10 4.50 ## 2 -65 -9.4 4.48 ## 3 -65 -8.8 4.45 ## 4 -65 -8.2 4.43 ## 5 -65 -7.6 4.41 ## 6 -65 -7 4.38 ## 7 -65 -6.4 4.36 ## 8 -65 -5.8 4.34 ## 9 -65 -5.2 4.31 ## 10 -65 -4.6 4.29 ## # … with 40,391 more rows
事前分布のときと同様にして処理します。
精度行列の期待値$\mathbb{E}[\boldsymbol{\Lambda}]$の逆行列を用いて、楕円の軸を計算します。
# Λの期待値の固有値・固有ベクトルを計算 posterior_eigen <- eigen(solve(E_lambda_hat_dd)) posterior_lmd_d <- posterior_eigen[["values"]] posterior_u_dd <- posterior_eigen[["vectors"]] |> t() # Λの期待値による楕円の軸を計算 posterior_axis_df <- tibble::tibble( xend = mu_d[1] + posterior_u_dd[, 1] * sqrt(posterior_lmd_d), # x軸の値 yend = mu_d[2] + posterior_u_dd[, 2] * sqrt(posterior_lmd_d) # y軸の値 ) posterior_axis_df
## # A tibble: 2 × 2 ## xend yend ## <dbl> <dbl> ## 1 -6.07 54.4 ## 2 22.3 30.7
事前分布のときと同様にして処理します。
$\boldsymbol{\Lambda}$の事後分布の期待値(精度行列の期待値)によるマハラノビス距離の等高線とその軸のグラフを作成します。
# パラメータラベル用の文字列を作成 posterior_param_text <- paste0( "list(N==", N, ", hat(nu)==", nu_hat, ", hat(W)==(list(", paste(w_hat_dd, collapse = ", "), ")))" ) # Λの事後分布の期待値による距離と軸を作図 ggplot() + geom_contour(data = model_delta_df, aes(x = x_1, y = x_2, z = delta, color = ..level..), alpha = 0.5, linetype = "dashed") + # 真のΛによる距離 geom_contour(data = posterior_delta_df, aes(x = x_1, y = x_2, z = delta, color = ..level..)) + # Λの期待値による距離 geom_segment(data = model_axis_df, mapping = aes(x = mu_d[1], y = mu_d[2], xend = xend, yend = yend, linetype = "model"), color = "red", size = 1, arrow = arrow(length = unit(10, "pt"))) + # 真のΛによる軸 geom_segment(data = posterior_axis_df, mapping = aes(x = mu_d[1], y = mu_d[2], xend = xend, yend = yend, linetype = "posterior"), color = "blue", size = 1, arrow = arrow(length = unit(10, "pt"))) + # Λの期待値による軸 scale_linetype_manual(breaks = c("posterior", "model"), values = c("solid", "dashed"), labels = c("posterior", "true model"), name = "axis") + # (凡例表示用の黒魔術) guides(linetype = guide_legend(override.aes = list(color = c("blue", "red"), size = c(0.5, 0.5), linetype = c("solid", "dashed")))) + # (凡例表示用の黒魔術) coord_fixed(ratio = 1) + # アスペクト比 labs(title = "Mahalanobis Distance", subtitle = parse(text = posterior_param_text), color = "distance", x = expression(x[1]), y = expression(x[2]))
軸の方向が近付いていることから精度パラメータを推定できているのが分かります。
予測分布の計算
最後に、(観測データ$\mathbf{X}$から求めた)事後分布のパラメータ$\hat{\nu}, \hat{\mathbf{W}}$から未観測のデータ$\mathbf{x}_{*}$の予測分布$p(\mathbf{x}_{*} | \mathbf{X}) = \mathrm{St}(\mathbf{x}_{*} | \boldsymbol{\mu}_s, \hat{\boldsymbol{\Lambda}}_s, \hat{\nu}_s)$を求めます。予測分布は多次元スチューデントのt分布になります。
$\boldsymbol{\Lambda}$の事後分布のパラメータ$\hat{\nu}, \hat{\mathbf{W}}$を用いて、予測分布($D$次元スチューデントのt分布)のパラメータ$\boldsymbol{\mu}_s, \hat{\boldsymbol{\Lambda}}_s, \hat{\nu}_s$を計算します。
# 予測分布の位置ベクトルを計算:式(3.124') mu_s_d <- mu_d # 予測分布の逆スケール行列を計算:式(3.124') lambda_s_hat_dd <- (1 - D + nu_hat) * w_hat_dd # 予測分布の自由度を計算:式(3.124') nu_s_hat <- 1 - D + nu_hat mu_s_d; lambda_s_hat_dd; nu_s_hat
## [1] 25 50 ## [,1] [,2] ## [1,] 0.0010370089 0.0002209955 ## [2,] 0.0002209955 0.0025705402 ## [1] 101
予測分布のパラメータを次の式で計算します。
予測分布を計算します。
# 予測分布(多次元t分布)を計算:式(3.121) predict_dens_df <- tibble::tibble( x_1 = x_mat[, 1], # x軸の値 x_2 = x_mat[, 2], # y軸の値 dens = mvnfast::dmvt(X = x_mat, mu = mu_s_d, sigma = solve(lambda_s_hat_dd), df = nu_s_hat) # 確率密度 ) predict_dens_df
## # A tibble: 40,401 × 3 ## x_1 x_2 dens ## <dbl> <dbl> <dbl> ## 1 -65 -10 0.0000000230 ## 2 -65 -9.4 0.0000000252 ## 3 -65 -8.8 0.0000000275 ## 4 -65 -8.2 0.0000000300 ## 5 -65 -7.6 0.0000000327 ## 6 -65 -7 0.0000000356 ## 7 -65 -6.4 0.0000000388 ## 8 -65 -5.8 0.0000000422 ## 9 -65 -5.2 0.0000000459 ## 10 -65 -4.6 0.0000000499 ## # … with 40,391 more rows
x_matの値ごとに確率密度を計算します。多次元のスチューデントのt分布の確率密度は、mvnfastパッケージのdmvt()で計算できます。確率変数の引数Xにx_mat、位置ベクトルの引数muにmu_s_d、スケール行列の引数sigmaにlambda_s_hat_ddの逆行列、自由度の引数dfにnu_s_hatを指定します。lambda_s_hat_ddは逆スケール行列です。
予測分布のグラフを真の分布と重ねて作成します。
# パラメータラベル用の文字列を作成 predict_param_text <- paste0( "list(N==", N, ", mu[s]==(list(", paste0(mu_s_d, collapse = ", "), "))", ", hat(Lambda)[s]==(list(", paste(round(lambda_s_hat_dd, 5), collapse = ", "), "))", ", hat(nu)[s]==", nu_s_hat, ")" ) # 予測分布を作図 ggplot() + geom_contour(data = model_dens_df, aes(x = x_1, y = x_2, z = dens, color = ..level..), alpha = 0.5, linetype = "dashed") + # 真の分布 geom_contour(data = predict_dens_df, aes(x = x_1, y = x_2, z = dens, color = ..level..)) + # 予測分布:(等高線図) labs(title = "Multivariate Student's t Distribution", subtitle = parse(text = predict_param_text), color = "density", x = expression(x[1]), y = expression(x[2]))
真の分布の等高線を破線で示します。
観測データが増えると、予測分布が真の分布に近付きます。
以上で、ガウス分布のベイズ推論を実装できました。
学習推移の可視化
gganimateパッケージを利用して、事後分布と予測分布の推移をアニメーション(gif画像)で確認するためのコードです。
モデルの設定
真の分布のパラメータ$\boldsymbol{\mu}, \boldsymbol{\Lambda}$と、$\boldsymbol{\Lambda}$の事前分布のパラメータ(の初期値)$\nu, \mathbf{W}$、試行回数$N$を設定します。
# 次元数を指定 D <- 2 # (既知の)平均ベクトルを指定 mu_d <- c(25, 50) # 真の分散共分散行列を指定 sigma_truth_dd <- matrix(c(900, -100, -100, 400), nrow = D, ncol = D) # 真の精度行列を計算 lambda_truth_dd <- solve(sigma_truth_dd) # Λの事前分布の自由度を指定 nu <- D # Λの事前分布の逆スケール行列を指定 w_dd <- diag(D) * 0.00005 # データ数(試行回数)を指定 N <- 100
「ベイズ推論の実装」と同じく、パラメータを指定します。
・コード(クリックで展開)
真の分布の確率変数がとり得る値$\mathbf{x}$を作成します。
# グラフ用のxの値を作成 x_1_vec <- seq( mu_d[1] - sqrt(sigma_truth_dd[1, 1]) * 3, mu_d[1] + sqrt(sigma_truth_dd[1, 1]) * 3, length.out = 201 ) x_2_vec <- seq( mu_d[2] - sqrt(sigma_truth_dd[2, 2]) * 3, mu_d[2] + sqrt(sigma_truth_dd[2, 2]) * 3, length.out = 201 ) # グラフ用のxの点を作成 x_mat <- tidyr::expand_grid(x_1 = x_1_vec, x_2 = x_2_vec) |> # 格子点を作成 as.matrix() # マトリクスに変換
確率変数がとり得る値(x軸とy軸の値)を作成します。グラフが粗い場合や処理が重い場合は、mu_*_vec, x_*_vecの要素数(by引数の値)または値の間隔(length.out引数の値)を調整してください。
事後分布と予測分布のパラメータの計算(更新)処理について、2通りの方法を載せます。目的に応じて使い分けてください。
推論処理：for関数による処理
1つ目の処理方法では、for()を使って、1データずつ生成してパラメータの更新を繰り返し処理します。こちらの方が、前ステップで求めた事後分布(のパラメータ)を次ステップの事前分布(のパラメータ)として用いる逐次学習をイメージしやすいです。
・コード(クリックで展開)
超パラメータの初期値を使って、事前分布の期待値によるマハラノビス距離とその軸を計算します。
# Λの期待値を計算:式(2.89) E_lambda_dd <- nu * w_dd # Λの期待値によるマハラノビス距離を計算 anime_posterior_delta_df <- tibble::tibble( x_1 = x_mat[, 1], # x軸の値 x_2 = x_mat[, 2], # y軸の値 delta = mahalanobis(x = x_mat, center = mu_d, cov = E_lambda_dd, inverted = TRUE) |> sqrt(), # 距離 param = paste0( "N=", 0, ", nu=", nu, ", W=(", paste0(round(w_dd, 6), collapse = ", "), ")" ) |> as.factor() # フレーム切替用ラベル ) # Λの期待値の固有値・固有ベクトルを計算 prior_eigen <- eigen(solve(E_lambda_dd)) prior_lmd_d <- prior_eigen[["values"]] prior_u_dd <- prior_eigen[["vectors"]] |> t() # Λの期待値による楕円の軸を計算 anime_posterior_axis_df <- tibble::tibble( xend = mu_d[1] + prior_u_dd[, 1] * sqrt(prior_lmd_d), yend = mu_d[2] + prior_u_dd[, 2] * sqrt(prior_lmd_d), param = paste0( "N=", 0, ", nu=", nu, ", W=(", paste0(round(w_dd, 6), collapse = ", "), ")" ) |> as.factor() # フレーム切替用ラベル ) anime_posterior_delta_df; anime_posterior_axis_df
## # A tibble: 40,401 × 4 ## x_1 x_2 delta param ## <dbl> <dbl> <dbl> <fct> ## 1 -65 -10 1.08 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 2 -65 -9.4 1.08 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 3 -65 -8.8 1.08 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 4 -65 -8.2 1.07 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 5 -65 -7.6 1.07 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 6 -65 -7 1.07 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 7 -65 -6.4 1.06 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 8 -65 -5.8 1.06 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 9 -65 -5.2 1.06 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 10 -65 -4.6 1.05 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## # … with 40,391 more rows ## # A tibble: 2 × 3 ## xend yend param ## <dbl> <dbl> <fct> ## 1 25 150 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 2 -75 50 N=0, nu=2, W=(5e-05, 0, 0, 5e-05)
現在のパラメータを文字列結合し因子型に変換して、フレーム切替用のラベル列とします。
また、事前分布のパラメータを使って、予測分布のパラメータと予測分布を計算します。
# 初期値による予測分布のパラメータを計算:式(3.124) mu_s_d <- mu_d lambda_s_dd <- (nu - 1) * w_dd nu_s <- nu - 1 # 予測分布(多次元t分布)を計算:式(3.121) anime_predict_dens_df <- tibble::tibble( x_1 = x_mat[, 1], x_2 = x_mat[, 2], dens = mvnfast::dmvt(X = x_mat, mu = mu_s_d, sigma = solve(lambda_s_dd), df = nu_s), # 確率密度 param = paste0( "N=", 0, ", mu_s=(", paste0(mu_s_d, collapse = ", "), ")", ", lambda_s=(", paste(round(lambda_s_dd, 5), collapse = ", "), ")", ", nu_s=", nu_s ) |> as.factor() # フレーム切替用ラベル ) anime_predict_dens_df
## # A tibble: 40,401 × 4 ## x_1 x_2 dens param ## <dbl> <dbl> <dbl> <fct> ## 1 -65 -10 0.00000399 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 2 -65 -9.4 0.00000400 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 3 -65 -8.8 0.00000401 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 4 -65 -8.2 0.00000403 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 5 -65 -7.6 0.00000404 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 6 -65 -7 0.00000406 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 7 -65 -6.4 0.00000407 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 8 -65 -5.8 0.00000408 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 9 -65 -5.2 0.00000409 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 10 -65 -4.6 0.00000411 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## # … with 40,391 more rows
同様に、フレーム切替用のラベル列を作成します。
パラメータの更新処理をN回繰り返します。
# 観測データの受け皿を初期化 x_nd <- matrix(NA, nrow = N, ncol = D) # ベイズ推論 for(n in 1:N) { # 多次元ガウス分布に従うデータを生成 x_nd[n, ] <- mvnfast::rmvn(n = 1, mu = mu_d, sigma = sigma_truth_dd) |> as.vector() # Λの事後分布のパラメータを更新:式(3.116) nu <- 1 + nu w_dd <- solve((x_nd[n, ] - mu_d) %*% t(x_nd[n, ] - mu_d) + solve(w_dd)) # Λの期待値を計算:式(2.89) E_lambda_dd <- nu * w_dd # Λの期待値によるマハラノビス距離を計算 tmp_posterior_delta_df <- tibble::tibble( x_1 = x_mat[, 1], x_2 = x_mat[, 2], delta = mahalanobis(x = x_mat, center = mu_d, cov = E_lambda_dd, inverted = TRUE) |> sqrt(), # 距離 param = paste0( "N=", n, ", nu=", nu, ", W=(", paste0(round(w_dd, 6), collapse = ", "), ")" ) |> as.factor() # フレーム切替用ラベル ) # Λの期待値の固有値・固有ベクトルを計算 posterior_eigen <- eigen(solve(E_lambda_dd)) posterior_lmd_d <- posterior_eigen[["values"]] posterior_u_dd <- posterior_eigen[["vectors"]] |> t() # Λの期待値による楕円の軸を計算 tmp_posterior_axis_df <- tibble::tibble( xend = mu_d[1] + posterior_u_dd[, 1] * sqrt(posterior_lmd_d), yend = mu_d[2] + posterior_u_dd[, 2] * sqrt(posterior_lmd_d), param = paste0( "N=", n, ", nu=", nu, ", W=(", paste0(round(w_dd, 6), collapse = ", "), ")" ) |> as.factor() # フレーム切替用ラベル ) # 予測分布のパラメータを更新:式(3.124) #mu_s_d <- mu_d lambda_s_dd <- (nu - 1) * w_dd nu_s <- nu - 1 # 予測分布(多次元t分布)を計算:式(3.121) tmp_predict_dens_df <- tibble::tibble( x_1 = x_mat[, 1], x_2 = x_mat[, 2], dens = mvnfast::dmvt(X = x_mat, mu = mu_s_d, sigma = solve(lambda_s_dd), df = nu_s), param = paste0( "N=", n, ", mu_s=(", paste0(mu_s_d, collapse = ", "), ")", ", lambda_s=(", paste0(round(lambda_s_dd, 5), collapse = ", "), ")", ", nu_s=", nu_s ) |> as.factor() # フレーム切替用ラベル ) # 推論結果を結合 anime_posterior_delta_df <- rbind(anime_posterior_delta_df, tmp_posterior_delta_df) anime_posterior_axis_df <- rbind(anime_posterior_axis_df, tmp_posterior_axis_df) anime_predict_dens_df <- rbind(anime_predict_dens_df, tmp_predict_dens_df) # 動作確認 #print(paste0("n=", n, " (", round(n / N * 100, 1), "%)")) }
超パラメータに関して、$\hat{\nu}, \hat{\mathbf{W}}$に対応するnu_hat, w_hat_ddを新たに作るのではなく、nu, w_ddを繰り返し更新(上書き)していきます。これにより、事後分布のパラメータの計算式(3.116)の$N$と$\sum_{n=1}^N$の計算は、ループ処理によってN回繰り返し1とx_nd[n, ]を加えることで行います。n回目のループ処理のときには、n-1回分の1とx_nd[n, ]が既にnuとw_ddに加えられています。
更新したパラメータを使って事後分布と予測分布を計算して、それぞれ計算結果をanime_***_dfに結合していきます。
結果を確認します。
# 確認 head(x_nd); anime_posterior_delta_df; anime_posterior_axis_df; anime_predict_dens_df
## [,1] [,2] ## [1,] 30.384191 55.103412 ## [2,] 17.653001 5.724338 ## [3,] -6.928976 50.954658 ## [4,] 28.284898 42.151756 ## [5,] 48.905975 47.723875 ## [6,] 37.998124 33.089360 ## # A tibble: 4,080,501 × 4 ## x_1 x_2 delta param ## <dbl> <dbl> <dbl> <fct> ## 1 -65 -10 1.08 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 2 -65 -9.4 1.08 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 3 -65 -8.8 1.08 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 4 -65 -8.2 1.07 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 5 -65 -7.6 1.07 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 6 -65 -7 1.07 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 7 -65 -6.4 1.06 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 8 -65 -5.8 1.06 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 9 -65 -5.2 1.06 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 10 -65 -4.6 1.05 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## # … with 4,080,491 more rows ## # A tibble: 202 × 3 ## xend yend param ## <dbl> <dbl> <fct> ## 1 25 150 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 2 -75 50 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 3 -34.3 -6.25 N=1, nu=3, W=(5e-05, 0, 0, 5e-05) ## 4 81.2 -9.26 N=1, nu=3, W=(5e-05, 0, 0, 5e-05) ## 5 38.1 123. N=2, nu=4, W=(5e-05, -1e-06, -1e-06, 4.5e-05) ## 6 -44.6 62.5 N=2, nu=4, W=(5e-05, -1e-06, -1e-06, 4.5e-05) ## 7 45.6 113. N=3, nu=5, W=(4.7e-05, -1e-06, -1e-06, 4.5e-05) ## 8 -36.6 70.1 N=3, nu=5, W=(4.7e-05, -1e-06, -1e-06, 4.5e-05) ## 9 41.9 108. N=4, nu=6, W=(4.7e-05, -1e-06, -1e-06, 4.5e-05) ## 10 -31.9 66.5 N=4, nu=6, W=(4.7e-05, -1e-06, -1e-06, 4.5e-05) ## # … with 192 more rows ## # A tibble: 4,080,501 × 4 ## x_1 x_2 dens param ## <dbl> <dbl> <dbl> <fct> ## 1 -65 -10 0.00000399 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 2 -65 -9.4 0.00000400 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 3 -65 -8.8 0.00000401 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 4 -65 -8.2 0.00000403 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 5 -65 -7.6 0.00000404 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 6 -65 -7 0.00000406 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 7 -65 -6.4 0.00000407 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 8 -65 -5.8 0.00000408 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 9 -65 -5.2 0.00000409 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## 10 -65 -4.6 0.00000411 N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu… ## # … with 4,080,491 more rows
それぞれ「x_matの行数」掛ける「N+1」行のデータフレームになります。行数が増えすぎないように注意してください。
推論処理：tidyverseパッケージによる処理
2つ目の処理方法では、tidyverseパッケージの関数を使って、一度の処理でN+1回分の計算をします。こちらの方が、処理時間が短いです(が、こんな処理をしていいのかは知りません)。
・コード(クリックで展開)
ガウス分布に従う$N$個のデータを生成します。
# 多次元ガウス分布に従うデータを生成 x_nd <- mvnfast::rmvn(n = N, mu = mu_d, sigma = sigma_truth_dd) head(x_nd)
# [,1] [,2] ## [1,] 30.384191 55.103412 ## [2,] 17.653001 5.724338 ## [3,] -6.928976 50.954658 ## [4,] 28.284898 42.151756 ## [5,] 48.905975 47.723875 ## [6,] 37.998124 33.089360
事前分布を含めたN+1回分の事後分布のパラメータを計算します。
# 試行ごとにΛの事後分布のパラメータを計算 anime_posterior_param_df <- tibble::tibble( n = 0:N, # データ番号(試行回数) nu_hat = n + nu ) |> dplyr::group_by(n) |> # パラメータの計算用にグループ化 dplyr::mutate( tmp_x_lt = dplyr::if_else( n > 0, # 事前分布を除く true = t(t(x_nd) - mu_d)[0:n, , drop = FALSE] |> list(), # リストに格納 false = matrix(0, nrow = 1, ncol = D) |> list() # リストに格納 ), w_hat_lt = dplyr::if_else( n > 0, # 事前分布を除く true = solve(t(tmp_x_lt[[1]]) %*% tmp_x_lt[[1]] + solve(w_dd)) |> list(), # リストに格納 false = w_dd |> list() # リストに格納 ) ) |> dplyr::ungroup() # グループ化を解除 anime_posterior_param_df
## # A tibble: 101 × 4 ## n nu_hat tmp_x_lt w_hat_lt ## <int> <dbl> <list> <list> ## 1 0 2 <dbl [1 × 2]> <dbl [2 × 2]> ## 2 1 3 <dbl [1 × 2]> <dbl [2 × 2]> ## 3 2 4 <dbl [2 × 2]> <dbl [2 × 2]> ## 4 3 5 <dbl [3 × 2]> <dbl [2 × 2]> ## 5 4 6 <dbl [4 × 2]> <dbl [2 × 2]> ## 6 5 7 <dbl [5 × 2]> <dbl [2 × 2]> ## 7 6 8 <dbl [6 × 2]> <dbl [2 × 2]> ## 8 7 9 <dbl [7 × 2]> <dbl [2 × 2]> ## 9 8 10 <dbl [8 × 2]> <dbl [2 × 2]> ## 10 9 11 <dbl [9 × 2]> <dbl [2 × 2]> ## # … with 91 more rows
観測データ数(試行回数)を表す0からNまでのN+1個の整数を使って、事後分布のパラメータ$\hat{\nu}, \hat{\mathbf{W}}$それぞれ上手いこと計算してリストに格納します(もっといい処理方法があれば教えてください)。
N+1回分の精度行列の期待値を用いて、マハラノビス距離を計算します。
# 試行ごとにΛの期待値によるマハラノビス距離を計算 anime_posterior_delta_df <- tidyr::expand_grid( n = 0:N, # データ番号(試行回数) x_1 = x_1_vec, x_2 = x_2_vec, ) |> # 試行ごとに格子点を複製 dplyr::group_by(n) |> # 距離の計算用にグループ化 dplyr::mutate( delta = mahalanobis( x = x_mat, center = mu_d, cov = anime_posterior_param_df[["nu_hat"]][unique(n)+1] * anime_posterior_param_df[["w_hat_lt"]][[unique(n)+1]], # Λの期待値を計算:式(2.89) inverted = TRUE ) |> sqrt(), # 距離 param = paste0( "N=", unique(n), ", nu=", anime_posterior_param_df[["nu_hat"]][unique(n)+1], ", W=(", paste0(round(anime_posterior_param_df[["w_hat_lt"]][[unique(n)+1]], 6), collapse = ", "), ")" ) |> as.factor() # フレーム切替用ラベル ) |> dplyr::ungroup() anime_posterior_delta_df
## # A tibble: 4,080,501 × 5 ## n x_1 x_2 delta param ## <int> <dbl> <dbl> <dbl> <fct> ## 1 0 -65 -10 1.08 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 2 0 -65 -9.4 1.08 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 3 0 -65 -8.8 1.08 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 4 0 -65 -8.2 1.07 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 5 0 -65 -7.6 1.07 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 6 0 -65 -7 1.07 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 7 0 -65 -6.4 1.06 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 8 0 -65 -5.8 1.06 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 9 0 -65 -5.2 1.06 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 10 0 -65 -4.6 1.05 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## # … with 4,080,491 more rows
0からNまでの整数と、x軸とy軸の値x_1_vec, x_2_vecの全ての組み合わせをexpand_grid()で作成します。これにより、x_matに対応する行をN+1フレーム分に複製できます。
n列でグループ化することで、x_matごとに距離を計算できます。
確率変数(x_matの各行)と各試行(n列の値)の組み合わせごとに距離を計算します。
各試行における事後分布のパラメータをanime_posterior_param_dfから抽出して精度行列の期待値を計算して、マハラノビス距離を計算します。
また、パラメータごとにフレーム切替用のラベルを作成します。
因子型への変換処理では、無名関数function()の省略記法\()を使って、(\(引数){引数を使った具体的な処理})()としています。直前のパイプ演算子を%>%にすると、行全体(\(引数){処理})()を{}の中の処理(この例だとfactor(., levels = unique(.)))に置き換えられます(置き換えられるように引数名を.にしています)。
N+1回分の精度行列の期待値の逆行列を用いて、楕円の軸を計算します。
# 試行ごとにΛの期待値による楕円の軸を計算 anime_posterior_axis_df <- tidyr::expand_grid( n = 0:N, # データ番号(試行回数) name = c("y_1", "y_2") # 値の受け皿 ) |> # 試行ごとに受け皿を複製 dplyr::group_by(n) |> # 軸の計算用にグループ化 dplyr::mutate( # Λの期待値の固有値・固有ベクトルを計算 eigen_lt = (anime_posterior_param_df[["nu_hat"]][unique(n)+1] * anime_posterior_param_df[["w_hat_lt"]][[unique(n)+1]]) |> # Λの期待値を計算:式(2.89) solve() |> # 分散共分散行列に変換 eigen() |> list(), # リストに格納 # Λの期待値による楕円の軸を計算 xend = mu_d[1] + eigen_lt[[1]][["vectors"]][1, ] * sqrt(eigen_lt[[1]][["values"]]), # x軸の値 yend = mu_d[2] + eigen_lt[[1]][["vectors"]][2, ] * sqrt(eigen_lt[[1]][["values"]]), # y軸の値 param = paste0( "N=", unique(n), ", nu=", anime_posterior_param_df[["nu_hat"]][unique(n)+1], ", W=(", paste0(round(anime_posterior_param_df[["w_hat_lt"]][[unique(n)+1]], 6), collapse = ", "), ")" ) |> (\(.){factor(., levels = unique(.))})() # フレーム切替用ラベル ) |> dplyr::ungroup() # グループ化を解除 anime_posterior_axis_df
## # A tibble: 202 × 6 ## n name eigen_lt xend yend param ## <int> <chr> <list> <dbl> <dbl> <fct> ## 1 0 y_1 <eigen> 25 150 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 2 0 y_2 <eigen> -75 50 N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 3 1 y_1 <eigen> -34.3 -6.25 N=1, nu=3, W=(5e-05, 0, 0, 5e-05) ## 4 1 y_2 <eigen> 81.2 -9.26 N=1, nu=3, W=(5e-05, 0, 0, 5e-05) ## 5 2 y_1 <eigen> 38.1 123. N=2, nu=4, W=(5e-05, -1e-06, -1e-06, 4.5e-… ## 6 2 y_2 <eigen> -44.6 62.5 N=2, nu=4, W=(5e-05, -1e-06, -1e-06, 4.5e-… ## 7 3 y_1 <eigen> 45.6 113. N=3, nu=5, W=(4.7e-05, -1e-06, -1e-06, 4.5… ## 8 3 y_2 <eigen> -36.6 70.1 N=3, nu=5, W=(4.7e-05, -1e-06, -1e-06, 4.5… ## 9 4 y_1 <eigen> 41.9 108. N=4, nu=6, W=(4.7e-05, -1e-06, -1e-06, 4.5… ## 10 4 y_2 <eigen> -31.9 66.5 N=4, nu=6, W=(4.7e-05, -1e-06, -1e-06, 4.5… ## # … with 192 more rows
こちらも、anime_posterior_param_dfから事後分布のパラメータを抽出して精度行列の期待値の逆行列を計算して、楕円の軸を計算します。
また、anime_posterior_delta_dfのラベル列と対応するように、フレーム切替用のラベルを作成します。
同様に、初期値を含めたN+1回分の予測分布を計算します。
# 試行ごとに予測分布を計算 anime_predict_dens_df <- tidyr::expand_grid( n = 0:N, # データ番号(試行回数) x_1 = x_1_vec, x_2 = x_2_vec ) |> # 試行ごとに格子点を複製 dplyr::group_by(n) |> # 分布の計算用にグループ化 dplyr::mutate( # 予測分布のパラメータを計算:式(3.124') lambda_s_lt = ((1 - D + unique(n) + nu) * anime_posterior_param_df[["w_hat_lt"]][[unique(n)+1]]) |> list(), # リストに格納 nu_s = 1 - D + unique(n) + nu, # 予測分布(多次元t分布)を計算:式(2.72) dens = mvnfast::dmvt( X = x_mat, mu = mu_d, sigma = solve(lambda_s_lt[[1]]), df = unique(nu_s) ), # 確率密度 param = paste0( "N=", unique(n), ", mu_s=(", paste0(mu_d, collapse = ", "), ")", ", lambda_s=(", paste0(round(lambda_s_lt[[1]], 5), collapse = ", "), ")", ", nu_s=", nu_s ) |> (\(.){factor(., levels = unique(.))})() # フレーム切替用ラベル ) |> dplyr::ungroup() # グループ化を解除 anime_predict_dens_df
## # A tibble: 4,080,501 × 7 ## n x_1 x_2 lambda_s_lt nu_s dens param ## <int> <dbl> <dbl> <list> <dbl> <dbl> <fct> ## 1 0 -65 -10 <dbl [2 × 2]> 1 0.00000399 N=0, mu_s=(25, 50), lambda_… ## 2 0 -65 -9.4 <dbl [2 × 2]> 1 0.00000400 N=0, mu_s=(25, 50), lambda_… ## 3 0 -65 -8.8 <dbl [2 × 2]> 1 0.00000401 N=0, mu_s=(25, 50), lambda_… ## 4 0 -65 -8.2 <dbl [2 × 2]> 1 0.00000403 N=0, mu_s=(25, 50), lambda_… ## 5 0 -65 -7.6 <dbl [2 × 2]> 1 0.00000404 N=0, mu_s=(25, 50), lambda_… ## 6 0 -65 -7 <dbl [2 × 2]> 1 0.00000406 N=0, mu_s=(25, 50), lambda_… ## 7 0 -65 -6.4 <dbl [2 × 2]> 1 0.00000407 N=0, mu_s=(25, 50), lambda_… ## 8 0 -65 -5.8 <dbl [2 × 2]> 1 0.00000408 N=0, mu_s=(25, 50), lambda_… ## 9 0 -65 -5.2 <dbl [2 × 2]> 1 0.00000409 N=0, mu_s=(25, 50), lambda_… ## 10 0 -65 -4.6 <dbl [2 × 2]> 1 0.00000411 N=0, mu_s=(25, 50), lambda_… ## # … with 4,080,491 more rows
0からNまでの整数と、x軸とy軸の値x_1_vec, x_2_vecの全ての組み合わせを作成して、n列の値ごとに(x_matごとに)予測分布のパラメータを計算して、確率密度を計算します。
作図処理
事後分布と予測分布の推移をそれぞれアニメーションで可視化します。
・コード(クリックで展開)
観測データと対応するラベルをデータフレームに格納します。
# 観測データを格納 anime_data_df <- tibble::tibble( x_1 = c(NA, x_nd[, 1]), # x軸の値 x_2 = c(NA, x_nd[, 2]), # y軸の値 param = anime_posterior_delta_df[["param"]] |> unique() # フレーム切替用ラベル ) anime_data_df
## # A tibble: 101 × 3 ## x_1 x_2 param ## <dbl> <dbl> <fct> ## 1 NA NA N=0, nu=2, W=(5e-05, 0, 0, 5e-05) ## 2 30.4 55.1 N=1, nu=3, W=(5e-05, 0, 0, 5e-05) ## 3 17.7 5.72 N=2, nu=4, W=(5e-05, -1e-06, -1e-06, 4.5e-05) ## 4 -6.93 51.0 N=3, nu=5, W=(4.7e-05, -1e-06, -1e-06, 4.5e-05) ## 5 28.3 42.2 N=4, nu=6, W=(4.7e-05, -1e-06, -1e-06, 4.5e-05) ## 6 48.9 47.7 N=5, nu=7, W=(4.6e-05, -1e-06, -1e-06, 4.5e-05) ## 7 38.0 33.1 N=6, nu=8, W=(4.6e-05, 0, 0, 4.5e-05) ## 8 21.6 33.6 N=7, nu=9, W=(4.6e-05, 0, 0, 4.4e-05) ## 9 -0.117 38.7 N=8, nu=10, W=(4.4e-05, -1e-06, -1e-06, 4.4e-05) ## 10 -8.10 50.0 N=9, nu=11, W=(4.2e-05, -1e-06, -1e-06, 4.4e-05) ## # … with 91 more rows
事前分布には観測データが影響しないので欠損値NAとして、anime_posterior_delta_dfまたはanime_posterior_axis_dfのラベル列と合わせて格納します。
参考として、各試行における観測データを描画するのに使います。
事後分布の期待値によるマハラノビス距離とその軸の推移のアニメーションを作成します。
# 真のΛによるマハラノビス距離を計算 model_delta_df <- tibble::tibble( x_1 = x_mat[, 1], # x軸の値 x_2 = x_mat[, 2], # y軸の値 delta = mahalanobis(x = x_mat, center = mu_d, cov = lambda_truth_dd, inverted = TRUE) |> sqrt() # 距離 ) # 真のΛの固有値・固有ベクトルを計算 model_eigen <- eigen(sigma_truth_dd) model_lmd_d <- model_eigen[["values"]] model_u_dd <- model_eigen[["vectors"]] |> t() # 真のΛによる楕円の軸を計算 model_axis_df <- tibble::tibble( xend = mu_d[1] + model_u_dd[, 1] * sqrt(model_lmd_d), # x軸の値 yend = mu_d[2] + model_u_dd[, 2] * sqrt(model_lmd_d) # y軸の値 ) # 等高線を引く値を指定 break_vals <- seq(0, ceiling(max(model_delta_df[["delta"]])), by = 0.5) # Λの事後分布の期待値による距離と軸のアニメーションを作図 anime_posterior_graph <- ggplot() + geom_contour(data = model_delta_df, aes(x = x_1, y = x_2, z = delta, color = ..level..), breaks = break_vals, alpha = 0.5, linetype = "dashed") + # 真のΛによる距離 geom_contour(data = anime_posterior_delta_df, aes(x = x_1, y = x_2, z = delta, color = ..level..), breaks = break_vals) + # Λの期待値による距離 geom_segment(data = model_axis_df, mapping = aes(x = mu_d[1], y = mu_d[2], xend = xend, yend = yend, linetype = "model"), color = "red", size = 1, arrow = arrow(length = unit(10, "pt"))) + # 真のΛによる軸 geom_segment(data = anime_posterior_axis_df, mapping = aes(x = mu_d[1], y = mu_d[2], xend = xend, yend = yend, linetype = "posterior"), color = "blue", size = 1, arrow = arrow(length = unit(10, "pt"))) + # Λの期待値による軸 geom_point(data = anime_data_df, mapping = aes(x = x_1, y = x_2), color = "orange", size = 3) + # 観測データ gganimate::transition_manual(param) + # フレーム scale_linetype_manual(breaks = c("posterior", "model"), values = c("solid", "dashed"), labels = c("posterior", "true model"), name = "axis") + # (凡例表示用の黒魔術) guides(linetype = guide_legend(override.aes = list(color = c("blue", "red"), size = c(0.5, 0.5), linetype = c("solid", "dashed")))) + # (凡例表示用の黒魔術) coord_fixed(ratio = 1, xlim = c(min(x_1_vec), max(x_1_vec)), ylim = c(min(x_2_vec), max(x_2_vec))) + # アスペクト比 labs(title = "Mahalanobis Distance", subtitle = "{current_frame}", color = "distance", x = expression(x[1]), y = expression(x[2])) # gif画像を作成 gganimate::animate(anime_posterior_graph, nframes = N+1+10, end_pause = 10, fps = 10, width = 800, height = 600)
フレームの順番を示す列をtransition_manual()に指定して、animate()でgif画像を作成します。事前分布(初期値)を含むため、フレーム数はN+1です。
続いて、anime_predict_dens_dfのラベル列を使って観測データを格納します。
# 観測データを格納 anime_data_df <- tibble::tibble( x_1 = c(NA, x_nd[, 1]), # x軸の値 x_2 = c(NA, x_nd[, 2]), # y軸の値 param = anime_predict_dens_df[["param"]] |> unique() # フレーム切替用ラベル ) anime_data_df
## # A tibble: 101 × 3 ## x_1 x_2 param ## <dbl> <dbl> <fct> ## 1 NA NA N=0, mu_s=(25, 50), lambda_s=(5e-05, 0, 0, 5e-05), nu_s=1 ## 2 30.4 55.1 N=1, mu_s=(25, 50), lambda_s=(1e-04, 0, 0, 1e-04), nu_s=2 ## 3 17.7 5.72 N=2, mu_s=(25, 50), lambda_s=(0.00015, 0, 0, 0.00014), nu_s=3 ## 4 -6.93 51.0 N=3, mu_s=(25, 50), lambda_s=(0.00019, 0, 0, 0.00018), nu_s=4 ## 5 28.3 42.2 N=4, mu_s=(25, 50), lambda_s=(0.00024, 0, 0, 0.00023), nu_s=5 ## 6 48.9 47.7 N=5, mu_s=(25, 50), lambda_s=(0.00028, 0, 0, 0.00027), nu_s=6 ## 7 38.0 33.1 N=6, mu_s=(25, 50), lambda_s=(0.00032, 0, 0, 0.00031), nu_s=7 ## 8 21.6 33.6 N=7, mu_s=(25, 50), lambda_s=(0.00037, 0, 0, 0.00035), nu_s=8 ## 9 -0.117 38.7 N=8, mu_s=(25, 50), lambda_s=(4e-04, -1e-05, -1e-05, 4e-04), nu… ## 10 -8.10 50.0 N=9, mu_s=(25, 50), lambda_s=(0.00042, -1e-05, -1e-05, 0.00044)… ## # … with 91 more rows
各試行における観測データを描画するのに使います。
また、各試行までの観測データを格納します。
# 観測データを複製して格納 anime_alldata_df <- tidyr::expand_grid( frame = 1:N, # フレーム番号 n = 1:N # 試行回数 ) |> # 全ての組み合わせを作成 dplyr::filter(n < frame) |> # フレーム番号以前のデータを抽出 dplyr::mutate( x_1 = x_nd[n, 1], # x軸の値 x_2 = x_nd[n, 2], # y軸の値 param = unique(anime_predict_dens_df[["param"]])[frame+1] # フレーム切替用ラベル ) anime_alldata_df
## # A tibble: 4,950 × 5 ## frame n x_1 x_2 param ## <int> <int> <dbl> <dbl> <fct> ## 1 2 1 30.4 55.1 N=2, mu_s=(25, 50), lambda_s=(0.00015, 0, 0, 0.00014… ## 2 3 1 30.4 55.1 N=3, mu_s=(25, 50), lambda_s=(0.00019, 0, 0, 0.00018… ## 3 3 2 17.7 5.72 N=3, mu_s=(25, 50), lambda_s=(0.00019, 0, 0, 0.00018… ## 4 4 1 30.4 55.1 N=4, mu_s=(25, 50), lambda_s=(0.00024, 0, 0, 0.00023… ## 5 4 2 17.7 5.72 N=4, mu_s=(25, 50), lambda_s=(0.00024, 0, 0, 0.00023… ## 6 4 3 -6.93 51.0 N=4, mu_s=(25, 50), lambda_s=(0.00024, 0, 0, 0.00023… ## 7 5 1 30.4 55.1 N=5, mu_s=(25, 50), lambda_s=(0.00028, 0, 0, 0.00027… ## 8 5 2 17.7 5.72 N=5, mu_s=(25, 50), lambda_s=(0.00028, 0, 0, 0.00027… ## 9 5 3 -6.93 51.0 N=5, mu_s=(25, 50), lambda_s=(0.00028, 0, 0, 0.00027… ## 10 5 4 28.3 42.2 N=5, mu_s=(25, 50), lambda_s=(0.00028, 0, 0, 0.00027… ## # … with 4,940 more rows
フレーム番号と試行回数(データ番号)の全ての組み合わせを作成して、フレーム番号未満のデータを抽出します。フレーム切替用のラベルは、フレーム番号を使って抽出します。
各試行以前のデータを描画するのに使うので、フレーム番号が2からになります。
予測分布の推移のアニメーションを作成します。
# 真の分布(多次元ガウス分布)を計算:式(2.72) model_dens_df <- tibble::tibble( x_1 = x_mat[, 1], # x軸の値 x_2 = x_mat[, 2], # y軸の値 dens = mvnfast::dmvn(X = x_mat, mu = mu_d, sigma = sigma_truth_dd) # 確率密度 ) # 予測分布を作図 anime_predict_graph <- ggplot() + geom_contour(data = anime_predict_dens_df, aes(x = x_1, y = x_2, z = dens, color = ..level..)) + # 予測分布 geom_contour(data = model_dens_df, aes(x = x_1, y = x_2, z = dens, color = ..level..), alpha = 0.5, linetype = "dashed") + # 真の分布 geom_point(data = anime_alldata_df, aes(x = x_1, y = x_2), color ="orange", alpha = 0.5) + # 過去の観測データ geom_point(data = anime_data_df, aes(x = x_1, y = x_2), color ="orange", size = 3) + # 観測データ gganimate::transition_manual(param) + # フレーム labs(title = "Multivariate Student's t Distribution", subtitle = "{current_frame}", color = "density", x = expression(x[1]), y = expression(x[2])) # gif画像を作成 gganimate::animate(anime_predict_graph, nframes = N+1+10, end_pause = 10, fps = 10, width = 800, height = 600)
データが増えるにつれて、真の精度行列による軸に近付いていくを確認できます。
新たなデータによって精度(分布の広がり具合)が推移しているのを確認できます。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
mvnfastパッケージにも多次元ガウス分布関連の関数がありますが、書き換えるのが面倒だったので止めました。
- 2021/04/11：加筆修正しました。その際に数式読解編と記事を分割しました。
mvtnormパッケージからmvnfastパッケージに置き替えました。
- 2022/09/21：加筆修正しました。
寄り道して共分散行列のあれこれを勉強したおかげで満足のいく形で事後分布を可視化できて良かったです。期待値を使ってるので、まだ分布感を表現できてないけど。
【次の節の内容】
