はじめに

　機械学習で登場する確率分布について色々な角度から理解したいシリーズです。

　この記事は、「R言語 - Qiita Advent Calendar 2024 - Qiita」の6日目の記事です。
　この記事では、R言語を使って2次元ガウス分布のグラフを作成します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
2次元ガウス分布の作図：base plot版
- グラフの作成
参考文献
おわりに

2次元ガウス分布の作図：base plot版

　2次元ガウス分布(Two-dimensional Gaussian Distribution)・二変量正規分布(Bivariate Normal Distribution)のグラフを作成します。この記事では、base plot(組み込み関数)を用いて作図します。
　多次元ガウス分布(正規分布)については「多次元ガウス分布の定義式 - からっぽのしょこ」、ggplot2パッケージを用いる場合は「【R】2次元ガウス分布の作図 - からっぽのしょこ」を参照してください。

グラフの作成

　base plot を利用して、ガウス分布のグラフを作成します。この例では、2次元のグラフで描画するため、次元数が2の場合のみを扱います。
　確率密度の計算については「【R】多次元ガウス分布の計算 - からっぽのしょこ」を参照してください。

パラメータの設定

　パラメータを設定します。

# 次元数を設定:(固定)
D <- 2

# 平均ベクトルを指定
mu_d <- rep(0, times = D) # 0ベクトル
mu_d <- c(6, 4)

# 分散共分散行列を指定
sigma_dd <- diag(D) # 単位行列
sigma_dd <- c(9, 2.5, 2.5, 4) |> 
  matrix(nrow = D, ncol = D, byrow = TRUE)
mu_d; sigma_dd

[1] 6 4
     [,1] [,2]
[1,]  9.0  2.5
[2,]  2.5  4.0

　次元数 $D = 2$ として、平均ベクトル $\boldsymbol{\mu} = (\mu_1, \mu_2)$ 、分散共分散行列 $\boldsymbol{\sigma} = ( (\sigma_1^2, \sigma_{1,2}), (\sigma_{2,1}, \sigma_2^2) )$ を指定します。
　 $x_d$ の平均 $\mu_d$ は実数、 $x_d$ の分散 $\sigma_{d,d} = \sigma_d^2$ は正の実数、 $x_i, x_j\ (i \neq j)$ の共分散 $\sigma_{i,j} = \sigma_{j,i}$ は実数、また $\boldsymbol{\Sigma}$ は正定値行列を満たす必要があります。

　パラメータを数式で表示するための文字列を作成します。

# パラメータラベルを作成
param_label <- paste0(
  "list(", 
    "mu == {}(", paste0(mu_d, collapse = ", "), "), ", 
    "Sigma == {}(", 
      "{}(", paste0(sigma_dd[1, ], collapse = ", "), "), ", 
      "{}(", paste0(sigma_dd[2, ], collapse = ", "), ")", 
    ")", 
  ")"
)
param_label

[1] "list(mu == {}(6, 4), Sigma == {}({}(9, 2.5), {}(2.5, 4)))"

　expression() の記法を使ってギリシャ文字などの記号を使った数式を描画します。等号は ==、複数の(数式上の)変数を並べるのは {}(変数1, 変数2) と書きます。(プログラム上の)変数の値を使う場合は、文字列で指定し parse() のtext 引数に渡します。

変数の設定

　確率変数がとり得る値を作成します。

# 確率変数の範囲の調整値
sgm_num <- 3

# 確率変数を作成
x_1_vec <- seq(
  from = mu_d[1] - sgm_num * sqrt(sigma_dd[1, 1]), 
  to   = mu_d[1] + sgm_num * sqrt(sigma_dd[1, 1]), 
  length.out = 101
)
x_2_vec <- seq(
  from = mu_d[2] - sgm_num * sqrt(sigma_dd[2, 2]), 
  to   = mu_d[2] + sgm_num * sqrt(sigma_dd[2, 2]), 
  length.out = 101
)
x_1_vec[1:5]; x_2_vec[1:5]

[1] -3.00 -2.82 -2.64 -2.46 -2.28
[1] -2.00 -1.88 -1.76 -1.64 -1.52

　2次元空間における点 $\mathbf{x} = (x_1, x_2)$ として用いる2つの軸の値 $x_1, x_2$ を作成します。
　この例では、各軸( $d$ 番目の軸・ $x_d$ 軸)の範囲を平均 $\mu_d$ を中心に標準偏差 $\sigma_d$ の sgm_num 倍とします。
　各軸の値(の組み合わせによる格子点)の数が曲面の滑らかさ(ポリゴンの細かさ)に影響します。処理が重い場合は、点の数を調整してください。

確率密度の計算

　確率密度関数を設定します。

# 確率密度関数を指定
calc_dens <- function(x1, x2, params) {
  
  # パラメータを取得
  mu    <- params[["mu"]]
  sigma <- params[["sigma"]]
  
  # 確率密度を計算
  mvnfast::dmvn(X = cbind(x1, x2), mu = mu, sigma = sigma) # 多次元ガウス分布
}

　格子状の点に対して計算を行うために、確率密度の計算処理を自作関数として定義しておきます。
　多次元ガウス分布の確率密度は、mvnfast パッケージの dmvn() で計算できます。確率変数の引数 X に点 $\mathbf{x}$ を各行とするマトリクス、平均ベクトルの引数 mu に $\boldsymbol{\mu}$ 、分散共分散行列の引数 sigma に $\boldsymbol{\Sigma}$ を指定します。

　確率密度を計算します。

# 確率密度を計算
param_lt <- list(mu = mu_d, sigma = sigma_dd) # パラメータを格納
dens_mat <- outer(X = x_1_vec, Y = x_2_vec, FUN = calc_dens, params = param_lt)
dens_mat[1:5, 1:5]

             [,1]         [,2]         [,3]         [,4]         [,5]
[1,] 5.081938e-05 5.757903e-05 6.495423e-05 7.295559e-05 8.158640e-05
[2,] 5.757903e-05 6.535632e-05 7.386165e-05 8.311099e-05 9.311207e-05
[3,] 6.495423e-05 7.386165e-05 8.362549e-05 9.426846e-05 1.058040e-04
[4,] 7.295559e-05 8.311099e-05 9.426846e-05 1.064590e-04 1.197034e-04
[5,] 8.158640e-05 9.311207e-05 1.058040e-04 1.197034e-04 1.348401e-04

　outer() を使って x_1_vec, x_2_vec それぞれの要素の組み合わせ(格子点)ごとに確率密度を計算します。
　dens_mat の行数は x_1_vec の要素数(x軸の点の数)、列数は x_2_vec の要素数(y軸の点の数)に対応します。

等高線図の作成

　等高線図を作成します。

# 2次元ガウス分布の等高線を作図
contour(
  x = x_1_vec, y = x_2_vec, z = dens_mat, 
  #nlevels = 20, 
  xlab = expression(x[1]), ylab = expression(x[2]), 
  main = "2D Gaussian Distribution", 
  sub = parse(text = param_label), 
  asp = 1
)

　等高線図は contour() で描画できます。
　等高線の値を指定する場合は levels 引数、等高線の数を指定する場合は nlevels 引数を使います。
　asp 引数にx軸・y軸のアスペクト比を指定できます。

　塗りつぶし等高線図を作成します。

# 2次元ガウス分布の塗りつぶしを作図
filled.contour(
  x = x_1_vec, y = x_2_vec, z = dens_mat, 
  #nlevels = 20, 
  xlab = expression(x[1]), ylab = expression(x[2]), 
  main = "2D Gaussian Distribution", 
  sub = parse(text = param_label), 
  key.title = title(main = "density"), 
  asp = 1
)

　色付けした等高線図は filled.contour() で描画できます。

　グラデーションの配色を設定します。

# 等高線の数を指定
lvl_num <- 13

# カラーパレットを指定
col.pal <- colorRampPalette(colors = c("blue", "green", "yellow", "red")) # 配色を指定

# カラーマップを作成
#col_map_vec <- cm.colors(n = lvl_num-1)
#col_map_vec <- terrain.colors(n = lvl_num-1)
col_map_vec <- col.pal(n = lvl_num-1)
col_map_vec[1:5]

[1] "#0000FF" "#0045B9" "#008B73" "#00D02E" "#17FF00"

　グラデーションに用いる色をcolorRampPalette() に指定します。配色に応じた色(カラーコード)を出力する関数が出力されます。
　グラデーションの分割数(階級数)を指定して、色ベクトルを作成します。　cm.colors() などの組み込みのカラーパレットを使うこともできます。

　グラデーションの階級を設定します。

# グラデーションの範囲を設定
u <- 0.005
dens_min <- 0
dens_max <- ceiling(max(dens_mat) /u)*u # u単位で切り上げ

# グラデーションの階級を設定
dens_level_vec <- seq(from = dens_min, to = dens_max, length.out = lvl_num)
dens_max; dens_level_vec[1:5]

[1] 0.03
[1] 0.0000 0.0025 0.0050 0.0075 0.0100

　グラデーション(z軸の値)の下限・上限を指定して、グラデーションの分割数で等間隔に区切ります。

　任意の配色の塗りつぶし等高線図を作成します。

# 2次元ガウス分布の等高線を作図
filled.contour(
  x = x_1_vec, y = x_2_vec, z = dens_mat, 
  zlim = c(dens_min, dens_max), 
  levels = dens_level_vec, col = col_map_vec, 
  plot.axes = {
    contour(x = x_1_vec, y = x_2_vec, z = dens_mat, levels = dens_level_vec, add = TRUE)
  }, 
  xlab = expression(x[1]), ylab = expression(x[2]), 
  main = "2D Gaussian Distribution", 
  sub = parse(text = param_label), 
  key.title = title(main = "density"), 
  asp = 1
)

　filled.contour() の col 引数に色ベクトルを指定します。
　plot.axes 引数に contour() を指定して等高線を描画できます。

曲面図の作成

　曲面図を作成します。

# 2次元ガウス分布の曲面を作図
persp(
  x = x_1_vec, y = x_2_vec, z = dens_mat, 
  border = "black", lwd = 1, 
  ticktype = "detailed", nticks = 10, 
  xlab = "x_1", ylab = "x_2", zlab = "density", 
  main = "2D Gaussian Distribution", 
  sub = parse(text = param_label), 
  scale = FALSE, expand = 200, 
  theta = 30, phi = 30
)

　曲面図は persp() で描画できます。x軸・y軸線に垂直な線の位置はそれぞれ x_1_vec, x_2_vec の値に対応します。
　デフォルト( scale = TRUE, expand = 1 )では立方体に収まるように各軸の縮尺が変更されます。scale 引数を TRUE にすると各軸の範囲に応じて縮尺が決まります。expand 引数はz軸のサイズ(倍率)を指定できます。
　theta 引数に水平方向の表示角度、phi 引数に垂直方向の表示角度を指定できます。
　軸ラベルに expression() は使えません。

　グラデーションの配色を設定します。

# グラデーションの分割数を指定
col_num <- 100

# カラーパレットを指定
col.pal <- colorRampPalette(colors = c("blue", "green", "yellow", "red")) # 配色を指定

# カラーマップを作成
#col_map_vec <- cm.colors(n = col_num)
#col_map_vec <- terrain.colors(n = col_num)
col_map_vec <- col.pal(n = col_num)
col_map_vec[1:5]

[1] "#0000FF" "#0007F7" "#000FEF" "#0017E7" "#001EE0"

　「等高線図の作成」と同様に処理します。

　色の割り当て用のz軸の値を計算します。

# 点の数を取得
x_1_num <- nrow(dens_mat)
x_2_num <- ncol(dens_mat)

# ポリゴンごとの平均確率密度を計算
dens_bar_mat <- 0.25 * (
  dens_mat[-x_1_num, -x_2_num] + dens_mat[-1, -x_2_num] + dens_mat[-x_1_num, -1] + dens_mat[-1, -1]
) # (左上 + 左下 + 右上 + 右下) / 4 
dens_bar_mat[1:5, 1:5]

             [,1]         [,2]         [,3]         [,4]         [,5]
[1,] 5.783344e-05 6.543781e-05 7.372061e-05 8.269126e-05 9.235083e-05
[2,] 6.543781e-05 7.417628e-05 8.371665e-05 9.407389e-05 1.052536e-04
[3,] 7.372061e-05 8.371665e-05 9.465536e-05 1.065587e-04 1.194382e-04
[4,] 8.269126e-05 9.407389e-05 1.065587e-04 1.201765e-04 1.349461e-04
[5,] 9.235083e-05 1.052536e-04 1.194382e-04 1.349461e-04 1.518057e-04

　各ポリゴンを形成する4点の平均確率密度を計算します。

　グラデーションの階級を設定します。

# グラデーションの範囲を設定
u <- 0.005
dens_min <- 0
dens_max <- ceiling(max(dens_mat) /u)*u # u単位で切り上げ

# グラデーションの階級を設定
dens_level_vec <- seq(from = dens_min, to = dens_max, length.out = col_num+1)
dens_max; dens_level_vec[1:5]

[1] 0.03
[1] 0.0000 0.0003 0.0006 0.0009 0.0012

　「等高線図の作成」と同様に処理します。

　z軸の値に応じて色を割り当てます。

# ポリゴンごとの色を割り当て
dens_range_vec <- cut(x = dens_bar_mat, breaks = dens_level_vec)
dens_range_vec[1:5]

(0,0.0003] (0,0.0003] (0,0.0003] (0,0.0003] (0,0.0003]
100 Levels: (0,0.0003] (0.0003,0.0006] (0.0006,0.0009] ... (0.0297,0.03]

　cut() を使って各ポリゴンの平均確率密度に対応する階級ベクトルを作成します。階級ごとに因子レベルが設定されます。

　塗りつぶし曲面図を作成します。

# 2次元ガウス分布の曲面を作図
persp(
  x = x_1_vec, y = x_2_vec, z = dens_mat, 
  col = col_map_vec[dens_range_vec], #shade = 0.05, 
  border = "black", lwd = 1, 
  ticktype = "detailed", nticks = 10, 
  xlab = "x_1", ylab = "x_2", zlab = "density", 
  main = "2D Gaussian Distribution", 
  sub = parse(text = param_label), 
  scale = FALSE, expand = 200, 
  theta = 30, phi = 30
)

　persp() の col 引数に各ポリゴンの色を指定できます。ポリゴンごとに、階級の因子レベルを添字として用いて色ベクトルからカラーコードを取り出して、色を割り当てます。

　回転するアニメーションでグラフを確認します。

作図コード(クリックで展開)

# 追加パッケージ
library(magick)

# 一時書き出しフォルダを指定
dir_path <- "tmp_folder"

# フレーム数を指定
frame_num <- 360

# 表示角度を作成
h_vals <- seq(from = 0, to = 360, length.out = frame_num+1)[1:frame_num] # 水平方向

# アングルごとに作図
for(frame_i in 1:frame_num) {
  
  # 表示角度を取得
  h <- h_vals[frame_i]
  
  # 書き出し開始
  file_path <- paste0(dir_path, "/", stringr::str_pad(frame_i, width = nchar(frame_num), pad = "0"), ".png")
  png(filename = file_path, width = 800, height = 800, units = "px")
  
  # 2次元ガウス分布を作図
  persp(
    x = x_1_vec, y = x_2_vec, z = dens_mat, 
    col = col_map_vec[dens_range_vec], #shade = 0.05, 
    border = "black", lwd = 1, 
    ticktype = "detailed", nticks = 10, 
    xlab = "x_1", ylab = "x_2", zlab = "density", 
    main = "2D Gaussian Distribution", 
    sub = parse(text = param_label), 
    scale = FALSE, expand = 200, 
    theta = h, phi = 30
  )
  
  # 書き出し終了
  dev.off() 
  
  # 途中経過を表示
  message("\r", frame_i, " / ", frame_num, appendLF = FALSE)
}

# 動画を作成
paste0(dir_path, "/", stringr::str_pad(1:frame_num, width = nchar(frame_num), pad = "0"), ".png") |> # ファイルパスを作成
  magick::image_read() |> # 画像ファイルを読込
  magick::image_animate(fps = 1, dispose = "previous") |> # gif画像を作成
  magick::image_write_video(path = "2d_gaussian_surf.mp4", framerate = 30) -> tmp_path # mp4ファイルを書出

　(base plot系の関数で直接アニメーションにすることもできるっぽいですが、1枚(フレーム)ずつ図を書き出して magick パッケージでやっつけました。)

　この記事では、base plotを用いてグラフを作成しました。次の記事では、ggplot2パッケージを用いて作成します。

Enjoy!

参考文献

R／RStudioでやさしく学ぶプログラミングとデータ分析【ダウンロード特典付き！】

作者:掌田津耶乃
マイナビ出版

Amazon

　3Dプロットについては扱っていませんが、base plotの色々なグラフについて解説されています。

おわりに

　もうアドカレの季節ですか……私がこの時期に扱うネタのレベルが年々下がっているのが気になってますが、今年は今のところ2つ投稿する予定です。
　そういえば今年こそは1人アドカレをしようと秋ごろにネタの仕込みをしていたのですが、別の準備に手を取られて忘れていました。

　そしてブログを開設して7年目の1つ目の記事でした。
　私はR歴7年半くらいのはずですが、tidyverseネイティブで生きてきたのでbase Rなあれこれはほとんど触ってきませんでした。近頃はmagrittrパイプよりbaseパイプを使うようにしていますが。
　特にきっかけがあったわけではないですが、(ネタもないことだし)試しにbase plotを使ってみました。拡張(?)派生(?)パッケージを使うと装飾などが充実するようですが、(キリがないことですし)今回は組み込みパッケージ縛りでやっています。使ってみて作図処理が早いという印象でした。
　でもやっぱりggplot2で3Dプロットを描き(書き)たーい。

　さて2024年12月6日は、モーニング娘。の10期メンバーの石田亜佑美さんの卒業の日です。

　投稿時点ではまだ卒コンは始まっていませんが、おめでとうございます！！！！
　あぁついにこの日が、10期メンバーのいないモーニング娘。になるんですね。あゆみんのリーダー姿を見たかったなんて言わないよ絶対…。

【次の内容】

からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

【R】2次元ガウス分布の作図【base plot】