はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
この記事では、R言語で多次元ガウス分布(多変量正規分布)から多次元ガウス分布を生成します。
【前の内容】
【他の記事一覧】
【この記事の内容】
多次元ガウス分布から確率分布の生成
多次元ガウス分布(Multivariate Gaussian Distribution)・多変量正規分布(Multivariate Normal Distribution)から多次元ガウス分布を生成します。多次元ガウス分布については「多次元ガウス分布の定義式 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用パッケージ library(tidyverse) library(mvnfast) library(patchwork)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。
また、patchworkパッケージの演算子を使うため、patchworkも読み込む必要があります。
magrittrパッケージのパイプ演算子%>%ではなく、ベースパイプ(ネイティブパイプ)演算子|>を使っています。%>%に置き換えても処理できます。
生成分布の設定
まずは、パラメータの生成分布としてガウス分布を設定して、ガウス分布の平均ベクトルを生成(サンプリング)します。生成分布を$\mathcal{N}(\boldsymbol{\mu} | \boldsymbol{\mu}_{\mathrm{gen}}, \boldsymbol{\Sigma}_{\mathrm{gen}})$、生成された分布を$\mathcal{N}(\mathbf{x} | \boldsymbol{\mu}, \boldsymbol{\Sigma})$で表すことにします。2次元ガウス分布のグラフ作成については「【R】2次元ガウス分布の作図 - からっぽのしょこ」を参照してください。
ガウス分布(生成分布)のパラメータ$\boldsymbol{\mu}_{\mathrm{gen}}, \boldsymbol{\Sigma}_{\mathrm{gen}}$とサンプルサイズ$N$を設定します。この例では、2次元のグラフで可視化するため、次元数を$D = 2$とします。分布の生成自体は次元数に関わらず行えます。
# (生成分布の)平均ベクトルを指定 mu_gen_d <- c(4, -2) # (生成分布の)分散共分散行列を指定 sigma_gen_dd <- matrix(c(1, 0.4, 0.4, 1.5), nrow = 2, ncol = 2) # 分布の数(サンプルサイズ)を指定 N <- 9
平均ベクトル$\boldsymbol{\mu}_{\mathrm{gen}} = (\mu_1, \mu_2)$、分散共分散行列$\boldsymbol{\Sigma}_{\mathrm{gen}} = (\sigma_1^2, \sigma_{2,1}, \sigma_{1,2}, \sigma_2^2)$とデータ数(パラメータ数)$N$を指定します。
平均$\mu_d$は実数、分散$\sigma_d^2$は正の実数、共分散$\sigma_{i,j} = \sigma_{j,i}$は実数、また$\boldsymbol{\Sigma}_{\mathrm{gen}}$は正定値行列を満たす必要があります。
ガウス分布に従う乱数を生成します。
# 多次元ガウス分布の平均ベクトルμを生成 mu_nd <- mvnfast::rmvn(n = N, mu = mu_gen_d, sigma = sigma_gen_dd) head(mu_nd)
## [,1] [,2] ## [1,] 3.900304 0.212285 ## [2,] 4.765049 -3.439237 ## [3,] 3.733442 -2.433367 ## [4,] 4.397805 -2.560559 ## [5,] 4.035782 -1.944171 ## [6,] 4.599251 -1.476575
多次元ガウス分布の乱数は、mvnfastパッケージのrmvn()で生成できます。データ数の引数nにN、平均の引数muにmu_gen_d、分散共分散行列の引数sigmaにsigma_gen_ddを指定します。
生成した値を、ガウス分布(生成された分布)の平均ベクトル$\boldsymbol{\mu}_n = (\mu_{n,1}, \mu_{n,2})$のN個のサンプルmu_ndとします。
生成したパラメータ(乱数)をデータフレームに格納します。
# パラメータを格納 param_df <- tibble::tibble( mu_1 = mu_nd[, 1], # x軸の値 mu_2 = mu_nd[, 2] # y軸の値 ) |> # 値を格納 dplyr::arrange(mu_1, mu_2) |> # 作図時の配置調整用に並べ替え dplyr::mutate( parameter = paste0("(", round(mu_1, 2), ", ", round(mu_2, 2), ")") |> factor() # 色分け用ラベル ) param_df
## # A tibble: 9 × 3 ## mu_1 mu_2 parameter ## <dbl> <dbl> <fct> ## 1 2.75 -1.04 (2.75, -1.04) ## 2 3.73 -2.43 (3.73, -2.43) ## 3 3.90 0.212 (3.9, 0.21) ## 4 4.04 -1.94 (4.04, -1.94) ## 5 4.40 -2.56 (4.4, -2.56) ## 6 4.60 -1.48 (4.6, -1.48) ## 7 4.77 -3.44 (4.77, -3.44) ## 8 4.92 -3.39 (4.92, -3.39) ## 9 5.26 -1.27 (5.26, -1.27)
mu_ndの各列(各次元)を列とするデータフレームを作成します。
グラフを分割して描画する際にx軸の値が小さい順に並ぶように、行を並べ替えてラベル列を作成します。
ガウス分布の期待値$\mathbb{E}[\boldsymbol{\mu}]$を計算します。
# パラメータの期待値を計算 E_mu_d <- mu_gen_d
$\boldsymbol{\mu}$の期待値$\mathbb{E}[\boldsymbol{\mu}] = \boldsymbol{\mu}_{\mathrm{gen}}$をE_mu_dとします。
これは、サンプリングされたパラメータとその分布の基準を示すのに使います。
設定したパラメータに応じて、ガウス分布の確率変数の値$\boldsymbol{\mu}$を作成して、確率密度を計算します。
# 確率変数μの値を作成 mu_1_vals <- seq( from = mu_gen_d[1] - sqrt(sigma_gen_dd[1, 1]) * 3, to = mu_gen_d[1] + sqrt(sigma_gen_dd[1, 1]) * 3, length.out = 100 ) mu_2_vals <- seq( from = mu_gen_d[2] - sqrt(sigma_gen_dd[2, 2]) * 3, to = mu_gen_d[2] + sqrt(sigma_gen_dd[2, 2]) * 3, length.out = 100 ) # 確率変数μの点を作成 mu_mat <- tidyr::expand_grid( mu_1 = mu_1_vals, mu_2 = mu_2_vals ) |> # 格子点を作成 as.matrix() # マトリクスに変換 # 多次元ガウス分布を計算 gaussian_generator_df <- tibble::tibble( mu_1 = mu_mat[, 1], # x軸の値 mu_2 = mu_mat[, 2], # y軸の値 density = mvnfast::dmvn(X = mu_mat, mu = mu_gen_d, sigma = sigma_gen_dd) # 確率密度 ) gaussian_generator_df
## # A tibble: 10,000 × 3 ## mu_1 mu_2 density ## <dbl> <dbl> <dbl> ## 1 1 -5.67 0.000156 ## 2 1 -5.60 0.000178 ## 3 1 -5.53 0.000203 ## 4 1 -5.45 0.000230 ## 5 1 -5.38 0.000260 ## 6 1 -5.30 0.000293 ## 7 1 -5.23 0.000329 ## 8 1 -5.15 0.000367 ## 9 1 -5.08 0.000408 ## 10 1 -5.01 0.000452 ## # … with 9,990 more rows
$\mu_1$(x軸)の値をmu_1_vals、$\mu_2$(y軸)の値をmu_2_valsとします。この例では、それぞれ平均を中心に標準偏差の3倍を範囲とします。
mu_1_valsとmu_2_valsの要素の全ての組み合わせ(格子状の点)をexpand_grid()で作成します。データフレームが出力されるので、as.matrix()でマトリクスに変換してmu_matとします。mu_matの各行が点$\boldsymbol{\mu} = (\mu_1, \mu_2)$に対応します。
パラメータの値を数式で表示するための文字列を作成します。
# パラメータラベルを作成:(数式表示用) generator_param_text <- paste0( "list(", "mu[paste(g,e,n)]==group('(', list(", paste0(mu_gen_d, collapse = ", "), "), ')')", ", Sigma[paste(g,e,n)]==group('(', list(", paste0(sigma_gen_dd, collapse = ", "), "), ')')", ")" ) generator_param_text
## [1] "list(mu[paste(g,e,n)]==group('(', list(4, -2), ')'), Sigma[paste(g,e,n)]==group('(', list(1, 0.4, 0.4, 1.5), ')'))"
ギリシャ文字などの記号を使った数式を表示する場合は、expression()の記法を使います。等号は"=="、複数の(数式上の)変数を並べるには"list(変数1, 変数2)"とします。(プログラム上の)変数の値を使う場合は、parse()のtext引数に指定します。
生成分布(ガウス分布)とパラメータのサンプルを描画します。
# 2次元ガウス分布を作図 gaussian_generator_graph <- ggplot() + geom_contour_filled(data = gaussian_generator_df, mapping = aes(x = mu_1, y = mu_2, z = density, fill = ..level..), alpha = 0.8) + # パラメータの生成分布 geom_point(mapping = aes(x = E_mu_d[1], y = E_mu_d[2]), color = "red", size = 6, shape = 4) + # パラメータの期待値 geom_point(data = param_df, mapping = aes(x = mu_1, y = mu_2, color = parameter), alpha = 0.8, size = 6) + # パラメータのサンプル labs(title ="Multivariate Gaussian Distribution", subtitle = parse(text = generator_param_text), color = expression(mu), fill = "density", x = expression(mu[1]), y = expression(mu[2])) gaussian_generator_graph
期待値をバツ印で示します。
以上で、生成分布を設定して、パラメータを生成しました。次は、パラメータのサンプルを用いて、ガウス分布を作図します。
分布の作図：多次元ガウス分布
次に、生成した値をガウス分布の平均ベクトルとして利用します。
ガウス分布(生成された分布)の分散共分散行列$\boldsymbol{\Sigma}$を設定します。
# (生成された分布の)分散共分散行列を指定 sigma_dd <- matrix(c(1, 0, 0, 1), nrow = 2, ncol = 2)
生成分布のときと同様に、値を指定します。
生成された分布の確率変数の値$\mathbf{x}$を作成します。
# 確率変数xの値を作成 x_1_vals <- mu_1_vals x_2_vals <- mu_2_vals # 確率変数xの点を作成 x_mat <- mu_mat
この例では、作図時に対応関係が分かりやすいように、$\mathbf{x} = (x_1, x_2)$の値としてmu_matの値を使います。
生成分布の期待値(生成された分布の平均パラメータの期待値)を用いて、ガウス分布を計算します。
# パラメータの期待値による多次元ガウス分布を計算 E_gaussian_df <- tibble::tibble( x_1 = x_mat[, 1], # x軸の値 x_2 = x_mat[, 2], # y軸の値 density = mvnfast::dmvn(X = x_mat, mu = E_mu_d, sigma = sigma_dd) # 確率密度 ) E_gaussian_df
## # A tibble: 10,000 × 3 ## x_1 x_2 density ## <dbl> <dbl> <dbl> ## 1 1 -5.67 0.00000207 ## 2 1 -5.60 0.00000271 ## 3 1 -5.53 0.00000353 ## 4 1 -5.45 0.00000458 ## 5 1 -5.38 0.00000590 ## 6 1 -5.30 0.00000756 ## 7 1 -5.23 0.00000963 ## 8 1 -5.15 0.0000122 ## 9 1 -5.08 0.0000154 ## 10 1 -5.01 0.0000193 ## # … with 9,990 more rows
mean引数にE_mu_dを指定します。
パラメータのサンプルごとにガウス分布を計算します。
# パラメータのサンプルごとに多次元ガウス分布を計算 gaussian_sample_df <- tidyr::expand_grid( n = 1:N, # データ番号 x_1 = x_1_vals, x_2 = x_2_vals ) |> # パラメータごとに確率変数xを複製 dplyr::mutate( mu_1 = mu_nd[n, 1], mu_2 = mu_nd[n, 2], ) |> # パラメータのサンプルを抽出 dplyr::arrange(mu_1, mu_2, x_1, x_2) |> # 作図時の配置調整用に並べ替え dplyr::group_by(n) |> # 分布の計算用にグループ化 dplyr::mutate( density = mvnfast::dmvn( X = x_mat, mu = unique(cbind(mu_1, mu_2)), sigma = sigma_dd ) # 確率密度 ) |> dplyr::ungroup() |> # グループ化を解除 dplyr::mutate( parameter = paste0("(", round(mu_1, 2), ", ", round(mu_2, 2), ")") |> (\(.){factor(., levels = unique(.))})() # 色分け用ラベル ) gaussian_sample_df
## # A tibble: 90,000 × 7 ## n x_1 x_2 mu_1 mu_2 density parameter ## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> ## 1 8 1 -5.67 2.75 -1.04 0.000000751 (2.75, -1.04) ## 2 8 1 -5.60 2.75 -1.04 0.00000106 (2.75, -1.04) ## 3 8 1 -5.53 2.75 -1.04 0.00000148 (2.75, -1.04) ## 4 8 1 -5.45 2.75 -1.04 0.00000206 (2.75, -1.04) ## 5 8 1 -5.38 2.75 -1.04 0.00000284 (2.75, -1.04) ## 6 8 1 -5.30 2.75 -1.04 0.00000391 (2.75, -1.04) ## 7 8 1 -5.23 2.75 -1.04 0.00000535 (2.75, -1.04) ## 8 8 1 -5.15 2.75 -1.04 0.00000729 (2.75, -1.04) ## 9 8 1 -5.08 2.75 -1.04 0.00000986 (2.75, -1.04) ## 10 8 1 -5.01 2.75 -1.04 0.0000133 (2.75, -1.04) ## # … with 89,990 more rows
パラメータ番号とする1からNの整数と確率変数x_1_vals, x_2_valsの要素の全ての組み合わせをexpand_grid()を作成します。これにより、サンプルごとにx_matを複製できます。
パラメータ列nでグループ化することで、x_matごとに確率密度を計算できます。
dmvn()のX引数にx_mat、mu引数にmu_ndの値、sigma引数にsigma_ddを指定します。
param_dfと対応するようにラベル列を作成します。
ここまでで、作図に用いるデータを作成できました。ここからは作図に関する処理を行います。
パラメータラベル用の文字列を作成します。
# パラメータラベルを作成:(数式表示用) sample_param_text <- paste0( "list(", "Sigma==group('(', list(", paste0(sigma_dd, collapse = ", "), "), ')')", ")" ) sample_param_text
## [1] "list(Sigma==group('(', list(1, 0, 0, 1), ')'))"
生成分布のときと同様に、文字列を指定します。
N個のガウス分布のグラフを分割して描画します。
# サンプルによる2次元ガウス分布を作図:分割して表示 gaussian_sample_graph <- ggplot() + geom_contour_filled(data = gaussian_sample_df, mapping = aes(x = x_1, y = x_2, z = density, fill = ..level..), alpha = 0.8) + # サンプルによる分布 geom_point(data = param_df, mapping = aes(x = mu_1, y = mu_2, color = parameter), size = 6, shape = 4, show.legend = FALSE) + # サンプルによる分布の期待値 facet_wrap(. ~ parameter, dir = "v", labeller = label_bquote(mu==.((as.character(parameter))))) + # グラフの分割 coord_cartesian(xlim = c(min(mu_1_vals), max(mu_1_vals)), ylim = c(min(mu_2_vals), max(mu_2_vals))) + # 軸の表示範囲 labs(title ="Multivariate Gaussian Distribution", subtitle = parse(text = sample_param_text), color = expression(mu), fill = "density", x = expression(x[1]), y = expression(x[2])) gaussian_sample_graph
facet_wrap()に列を指定すると、その列の値ごとにグラフを分割して描画できます。
パラメータの生成分布(ガウス分布)と生成された分布(ガウス分布)を並べて描画します。
# グラフを並べて描画 gaussian_generator_graph + gaussian_sample_graph
patchworkパッケージの+演算子を使うと左右に並べて描画できます。
続いて、グラフを重ねて表示することで、生成された分布の位置関係を比較しやすくします。
作図用に確率密度(z軸の値)の最大値を計算します。
# 確率密度の最大値を計算 max_dens <- mvnfast::dmvn(X = E_mu_d, mu = E_mu_d, sigma = sigma_dd) max_dens
## [1] 0.1591549
$\mathbf{x} = \boldsymbol{\mu}$の時の確率密度が最大になります。
パラメータラベル用の文字列を作成します。
# パラメータラベルを作成:(数式表示用) sample_param_text <- paste0( "list(", "E(mu)==group('(', list(", paste0(E_mu_d, collapse = ", "), "), ')')", ", Sigma==group('(', list(", paste0(sigma_dd, collapse = ", "), "), ')')", ")" ) sample_param_text
## [1] "list(E(mu)==group('(', list(4, -2), ')'), Sigma==group('(', list(1, 0, 0, 1), ')'))"
N+1個のガウス分布のグラフを重ねて描画します。
# サンプルによる2次元ガウス分布を作図:重ねて表示 gaussian_sample_graph <- ggplot() + geom_contour(data = E_gaussian_df, mapping = aes(x = x_1, y = x_2, z = density), breaks = max_dens*exp(-0.5), color = "red", linetype ="dashed") + # 期待値による分布 geom_point(mapping = aes(x = E_mu_d[1], y = E_mu_d[2]), color = "red", size = 6, shape = 4) + # 期待値による分布の期待値 geom_contour(data = gaussian_sample_df, mapping = aes(x = x_1, y = x_2, z = density, color = parameter), breaks = max_dens*exp(-0.5), show.legend = FALSE) + # サンプルによる分布 geom_point(data = param_df, mapping = aes(x = mu_1, y = mu_2, color = parameter), size = 6, shape = 4, show.legend = FALSE) + # サンプルによる分布の期待値 coord_cartesian(xlim = c(min(mu_1_vals), max(mu_1_vals)), ylim = c(min(mu_2_vals), max(mu_2_vals))) + # 軸の表示範囲 labs(title ="Multivariate Gaussian Distribution", subtitle = parse(text = sample_param_text), color = expression(mu), x = expression(x[1]), y = expression(x[2])) gaussian_sample_graph
geom_contour()のbreaks引数に値を指定すると、その値で等高線を引きます。この例では、確率密度が最大値の$\exp(-\frac{1}{2})$倍を指定します。
期待値による分布を破線で示します。
生成分布と生成された分布を並べて描画します。
# グラフを並べて描画 gaussian_generator_graph + gaussian_sample_graph + patchwork::plot_layout(guides = "collect") # 凡例の位置をまとめる
(定義通りですが)左図におけるパラメータのサンプルの点と、右図における分布の確率密度が最大となる点(期待値)が対応しているのが分かります。
この記事では、多次元ガウス分布からの分布生成を確認しました。
参考文献
- C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上』,丸善出版,2012年.
おわりに
スピンオフの分散共分散行列シリーズまで誕生した多次元ガウス分布シリーズもこれで一旦終わりです。いつかまた統計量の導出編などでお会いしましょう。
【次の内容】
