からっぽのしょこ

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

【R】多次元ガウス分布の乱数生成

はじめに

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

 この記事では、R言語で多次元ガウス分布(多変量正規分布)の乱数を生成します。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

多次元ガウス分布の乱数生成

 多次元ガウス分布(Maltivariate Gaussian Distribution)・多変量正規分布(Maltivariate Normal Distribution)の乱数を生成します。多次元ガウス分布については「多次元ガウス分布の定義式 - からっぽのしょこ」を参照してください。

 利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)
library(mvnfast)
library(gganimate)

 この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。
 magrittrパッケージのパイプ演算子%>%ではなく、ベースパイプ(ネイティブパイプ)演算子|>を使っています。%>%に置き換えても処理できます。
 分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。

サンプリング

 まずは、ガウス分布の乱数を生成します。

 ガウス分布のパラメータ$\boldsymbol{\mu}, \boldsymbol{\Lambda}$とデータ数$N$を設定します。この例では、2次元のグラフで描画するため、次元数を$D = 2$とします。乱数の生成自体は次元数に関わらず行えます。

# 平均ベクトルを指定
mu_d <- c(6, 10)

# 分散共分散行列を指定
sigma_dd <- matrix(c(1, 0.6, 0.6, 1.5), nrow = 2, ncol = 2)

# データ数(サンプルサイズ)を指定
N <- 10000

 平均ベクトル$\boldsymbol{\mu} = (\mu_1, \mu_2)$、分散共分散行列$\boldsymbol{\Sigma} = (\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}$は正定値行列を満たす必要があります。

 ガウス分布に従う乱数を生成します。

# 多次元ガウス分布に従う乱数を生成
x_nd <- mvnfast::rmvn(n = N, mu = mu_d, sigma = sigma_dd)

# サンプルを格納
data_df <- tibble::tibble(
  x_1 = x_nd[, 1], # x軸の値
  x_2 = x_nd[, 2] # y軸の値
)
data_df
## # A tibble: 10,000 × 2
##      x_1   x_2
##    <dbl> <dbl>
##  1  4.83  9.25
##  2  5.64  8.85
##  3  4.82 10.0 
##  4  8.25 12.3 
##  5  6.25 11.9 
##  6  5.49 12.2 
##  7  5.95 11.6 
##  8  4.87  8.27
##  9  7.17 11.6 
## 10  7.70 11.0 
## # … with 9,990 more rows

 多次元ガウス分布の乱数は、mvnfastパッケージのrmvn()で生成できます。データ数の引数nN、平均の引数mumu_d、分散共分散行列の引数sigmasigma_ddを指定します。
 生成した値をサンプル$\mathbf{x}_n = (x_{n,1}, x_{n,2})$とします。

乱数の可視化

 続いて、生成した乱数のグラフを作成します。ガウス分布のグラフ作成については「【R】2次元ガウス分布の作図 - からっぽのしょこ」を参照してください。

 設定したパラメータに応じて、ガウス分布の確率変数の値$\mathbf{x}$を作成して、確率密度を計算します。

# xの値を作成
x_1_vals <- seq(
  from = mu_d[1] - sqrt(sigma_dd[1, 1]) * 3, 
  to = mu_d[1] + sqrt(sigma_dd[1, 1]) * 3, 
  length.out = 101
)
x_2_vals <- seq(
  from = mu_d[2] - sqrt(sigma_dd[2, 2]) * 3, 
  to = mu_d[2] + sqrt(sigma_dd[2, 2]) * 3, 
  length.out = 101
)

# xの点を作成
x_mat <- tidyr::expand_grid(
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # 格子点を作成
  as.matrix() # マトリクスに変換

# 多次元ガウス分布を計算
dens_df <- tibble::tibble(
  x_1 = x_mat[, 1], # x軸の値
  x_2 = x_mat[, 2], # y軸の値
  density = mvnfast::dmvn(X = x_mat, mu = mu_d, sigma = sigma_dd) # 確率密度
)
dens_df
## # A tibble: 10,201 × 3
##      x_1   x_2  density
##    <dbl> <dbl>    <dbl>
##  1     3  6.33 0.000355
##  2     3  6.40 0.000399
##  3     3  6.47 0.000447
##  4     3  6.55 0.000499
##  5     3  6.62 0.000554
##  6     3  6.69 0.000612
##  7     3  6.77 0.000673
##  8     3  6.84 0.000736
##  9     3  6.91 0.000801
## 10     3  6.99 0.000869
## # … with 10,191 more rows

 $x_1$(x軸)の値をx_1_vals、$x_2$(y軸)の値をx_2_valsとします。この例では、それぞれ平均を中心に標準偏差の3倍を範囲とします。
 x_1_valsx_2_valsの要素の全ての組み合わせ(格子状の点)をexpand_grid()で作成します。データフレームが出力されるので、as.matrix()でマトリクスに変換してx_matとします。x_matの各行が点$\mathbf{x} = (x_1, x_2)$に対応します。
 多次元ガウス分布の確率密度は、dmvn()で計算できます。確率変数の引数Xx_mat、平均の引数mumu_d、共分散の引数sigmasigma_ddを指定します。

 サンプルの散布図を確率分布と重ねて作成します。

# パラメータラベルを作成:(数式表示用)
math_text <- paste0(
  "list(", 
  "mu==group('(', list(", paste0(mu_d, collapse = ", "), "), ')')", 
  ", Sigma==group('(', list(", paste0(sigma_dd, collapse = ", "), "), ')')", 
  ")"
)

# サンプルの散布図を作成
ggplot() + 
  geom_point(data = data_df, mapping = aes(x = x_1, y = x_2), 
             color = "orange", alpha = 0.3) + # サンプル
  geom_contour(data = dens_df, mapping = aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 元の分布
  labs(title ="Maltivariate Gaussian Distribution", 
       subtitle = parse(text = math_text), 
       color = "density", 
       x = expression(x[1]), y = expression(x[2]))

2次元ガウス分布の乱数の散布図


 サンプルの度数のヒートマップを作成します。

# サンプルの度数のヒートマップを作成
ggplot() + 
  geom_bin_2d(data = data_df, mapping = aes(x = x_1, y = x_2, fill = ..count..), 
              alpha = 0.8) + # サンプル
  geom_contour(data = dens_df, mapping = aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 元の分布
  scale_fill_distiller(palette = "Spectral") + # 塗りつぶしの色
  scale_color_distiller(palette = "Spectral") + # 等高線の色
  labs(title = "Maltivariate Gaussian Distribution", 
       subtitle = parse(text = math_text), 
       fill = "frequency", color = "density", 
       x = expression(x[1]), y = expression(x[2]))

2次元ガウス分布の乱数の度数のヒートマップ

 geom_bin_2d()で2次元の変数に対するヒートマップを描画できます。タイルの色の引数fill..count..を指定すると度数、..density..を指定すると密度に応じて色付けされます。

 密度の等高線図を作成します。

# サンプルの密度の等高線図を作成
ggplot() + 
  geom_density_2d_filled(data = data_df, mapping = aes(x = x_1, y = x_2, fill = ..level..), 
                         alpha = 0.8) + # サンプル
  geom_contour(data = dens_df, mapping = aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 元の分布
  scale_color_viridis_c(option = "D") + # 等高線の色
  labs(title = "Maltivariate Gaussian Distribution", 
       subtitle = parse(text = math_text), 
       fill = "density", color = "density", 
       x = expression(x[1]), y = expression(x[2]))

2次元ガウス分布の乱数の密度の等高線図

 geom_density_2d_filled()で密度の等高線図を描画できます。

 データ数が十分に増えると、ヒストグラムの形が分布の形に近付きます。

乱数と分布の関係をアニメーションで可視化

 次は、サンプルサイズとヒストグラムの形状の関係をアニメーションで確認します。

1データずつ可視化

 1データずつ乱数を生成するアニメーションを作成します。

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

 データ数を指定して、サンプルを生成します。

# データ数(フレーム数)を指定
N <- 100

# 多次元ガウス分布に従う乱数を生成
x_nd <- mvnfast::rmvn(n = N, mu = mu_d, sigma = sigma_dd)
head(x_nd)
##          [,1]      [,2]
## [1,] 6.709243  7.818683
## [2,] 4.155555  8.257296
## [3,] 5.502608  9.615666
## [4,] 7.670199 11.452761
## [5,] 5.087102  9.606752
## [6,] 4.831511  8.720070

 x_nn番目の要素を、n番目のデータ(n回目にサンプリングされた値)とみなします。アニメーションのn番目のフレームでは、n個のサンプルx_nd[1:n, ]のグラフを描画します。

 サンプルをデータフレームに格納します。

# サンプルを格納
anime_data_df <- tibble::tibble(
  n = 1:N, # データ番号
  x_1 = x_nd[, 1], # x軸の値
  x_2 = x_nd[, 2], # y軸の値
  parameter = paste0("mu=(", paste0(mu_d, collapse = ", "), "), Sigma=(", paste0(sigma_dd, collapse = ", "), "), N=", n) |> 
    factor(levels = paste0("mu=(", paste0(mu_d, collapse = ", "), "), Sigma=(", paste0(sigma_dd, collapse = ", "), "), N=", 1:N)) # フレーム切替用ラベル
)
anime_data_df
## # A tibble: 100 × 4
##        n   x_1   x_2 parameter                                 
##    <int> <dbl> <dbl> <fct>                                     
##  1     1  6.71  7.82 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=1 
##  2     2  4.16  8.26 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=2 
##  3     3  5.50  9.62 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=3 
##  4     4  7.67 11.5  mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=4 
##  5     5  5.09  9.61 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=5 
##  6     6  4.83  8.72 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=6 
##  7     7  4.99 11.5  mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=7 
##  8     8  6.44  8.16 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=8 
##  9     9  6.95 10.5  mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=9 
## 10    10  3.37  7.69 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=10
## # … with 90 more rows

 $\mathbf{x}$のサンプルx_ndをデータフレームに格納して、フレーム切替用のラベルを作成します。ラベルが文字列型だと文字列の基準で順序が決まるので、因子型にしてサンプリング回数に応じたレベル(順序)を設定します。
 このデータフレームは、各試行におけるサンプルを描画するのに使います。

 サンプリング回数ごとに、それまでのサンプルを持つデータフレームを作成します。

# サンプルを複製して格納
anime_freq_df <- tidyr::expand_grid(
  frame = 1:N, # フレーム番号
  n = 1:N # データ番号
) |> # 全ての組み合わせを作成
  dplyr::filter(n <= frame) |> # 各試行までのサンプルを抽出
  dplyr::mutate(
    x_1 = x_nd[n, 1], 
    x_2 = x_nd[n, 2], 
    parameter = paste0("mu=(", paste0(mu_d, collapse = ", "), "), Sigma=(", paste0(sigma_dd, collapse = ", "), "), N=", frame) |> 
      factor(levels = paste0("mu=(", paste0(mu_d, collapse = ", "), "), Sigma=(", paste0(sigma_dd, collapse = ", "), "), N=", 1:N)) # フレーム切替用ラベル
  )
anime_freq_df
## # A tibble: 5,050 × 5
##    frame     n   x_1   x_2 parameter                                
##    <int> <int> <dbl> <dbl> <fct>                                    
##  1     1     1  6.71  7.82 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=1
##  2     2     1  6.71  7.82 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=2
##  3     2     2  4.16  8.26 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=2
##  4     3     1  6.71  7.82 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=3
##  5     3     2  4.16  8.26 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=3
##  6     3     3  5.50  9.62 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=3
##  7     4     1  6.71  7.82 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=4
##  8     4     2  4.16  8.26 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=4
##  9     4     3  5.50  9.62 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=4
## 10     4     4  7.67 11.5  mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=4
## # … with 5,040 more rows

 フレーム番号とデータ番号(それぞれ1からNの整数)の全ての組み合わせをexpand_grid()で作成して、フレーム番号以下のデータ番号を抽出します。
 データ番号をインデックスとして使って、x_ndから対応するサンプルを抽出します。
 フレーム番号をデータ番号として、フレーム切替用のラベルを作成します。
 このデータフレームは、ヒートマップなどのグラフを描画するのに使います。

 サンプルの散布図のアニメーション(gif画像)を作成します。

# 散布図のアニメーションを作図
anime_freq_graph <- ggplot() + 
  geom_point(data = anime_freq_df, mapping = aes(x = x_1, y = x_2), 
             color = "orange", alpha = 0.5, size = 3) + # n個のサンプル
  geom_contour(data = dens_df, mapping = aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 元の分布
  geom_point(data = anime_data_df, mapping = aes(x = x_1, y = x_2), 
             color = "orange", size = 6) + # n番目のサンプル
  gganimate::transition_manual(parameter) + # フレーム
  labs(title ="Maltivariate Gaussian Distribution", 
       subtitle = "{current_frame}", 
       color = "density", 
       x = expression(x[1]), y = expression(x[2]))

# gif画像を作成
gganimate::animate(anime_freq_graph, nframes = N+10, end_pause = 10, fps = 10, width = 800, height = 600)

 transition_manual()にフレームの順序を表す列を指定します。この例では、因子型のラベルのレベルの順に描画されます。
 animate()のフレーム数の引数nframesにデータ数(サンプルサイズ)、フレームレートの引数fpsに1秒当たりのフレーム数を指定します。fps引数の値が大きいほどフレームが早く切り替わります。ただし、値が大きいと指定した通りに動作しません。

 サンプルの度数のヒートマップのアニメーションを作成します。

# メッシュ用の値を設定
x_1_breaks <- seq(from = min(x_1_vals), to = max(x_1_vals), length.out = 30)
x_2_breaks <- seq(from = min(x_2_vals), to = max(x_2_vals), length.out = 30)

# 度数のヒートマップのアニメーションを作図
anime_freq_graph <- ggplot() + 
  geom_bin_2d(data = anime_freq_df, mapping = aes(x = x_1, y = x_2, fill = ..count..), 
              breaks = list(x = x_1_breaks, y = x_2_breaks), 
              alpha = 0.8) + # n個のサンプル
  geom_contour(data = dens_df, mapping = aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 元の分布
  geom_point(data = anime_data_df, mapping = aes(x = x_1, y = x_2), 
             color = "orange", size = 6) + # n番目のサンプル
  gganimate::transition_manual(parameter) + # フレーム
  scale_fill_distiller(palette = "Spectral") + # 塗りつぶしの色
  scale_color_distiller(palette = "Spectral") + # 等高線の色
  labs(title = "Maltivariate Gaussian Distribution", 
       subtitle = "{current_frame}", 
       fill = "frequency", color = "density", 
       x = expression(x[1]), y = expression(x[2]))

# gif画像を作成
gganimate::animate(anime_freq_graph, nframes = N+10, end_pause = 10, fps = 10, width = 800, height = 600)

 集計範囲(メッシュ)を固定するために、geom_bin_2d()breaks引数に区切り位置を指定します。x軸とy軸を分けて指定する場合は、リストにして渡します。この例では、mu_vals, lambda_valsそれぞれの最小値から最大値の範囲を30分割します。(指定しなくてもよしなにしてくれるっぽい?)

 サンプルの密度の等高線のアニメーションを作成します。

# 除去するフレーム数を指定
n_min <- 10

# (データが少ないと密度を計算できないため)最初のフレームを除去
tmp_data_df <- anime_data_df|> 
  dplyr::filter(n > n_min) |> # 始めのデータを削除
  dplyr::mutate(
    parameter = parameter |> 
      as.character() |> 
      (\(.){factor(., levels = unique(.))})() # レベルを再設定
  )
tmp_freq_df <- anime_freq_df |> 
  dplyr::filter(frame > n_min) |> # 始めのデータを削除
  dplyr::mutate(
    parameter = parameter |> 
      as.character() |> 
      (\(.){factor(., levels = unique(.))})() # レベルを再設定
  )

# 密度の等高線図のアニメーションを作図
anime_freq_graph <- ggplot() + 
  geom_density_2d_filled(data = tmp_freq_df, mapping = aes(x = x_1, y = x_2, fill = ..level..), 
                         alpha = 0.8) + # n個のサンプル
  geom_contour(data = dens_df, mapping = aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 元の分布
  geom_point(data = tmp_data_df, mapping = aes(x = x_1, y = x_2), 
             color = "orange", size = 6) + # n番目のサンプル
  gganimate::transition_manual(parameter) + # フレーム
  scale_color_viridis_c(option = "D") + # 等高線の色
  labs(title = "Maltivariate Gaussian Distribution", 
       subtitle = "{current_frame}", 
       fill = "density", color = "density", 
       x = expression(x[1]), y = expression(x[2]))

# gif画像を作成
gganimate::animate(anime_freq_graph, nframes = N-n_min+10, end_pause = 10, fps = 10, width = 800, height = 600)

 データが1つだと密度を計算できず、またデータが少ないと密度が大きくなり(z軸の最大値が大きくなり)アニメーション全体でのグラデーション(等高線)が分かりにくくなるので、始めのフレームに対応する行を取り除いておきます。
 データフレームからデータを除いても、因子レベルの情報がアニメーションのフレームに影響してしまうので、因子の情報を再設定する必要があります。


2次元ガウス分布の乱数の散布図

2次元ガウス分布の乱数の度数のヒートマップ

2次元ガウス分布の乱数の密度の等高線図


複数データずつ可視化

 続いて、フレームごとに複数データを生成します。

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

 データ数とフレーム数を指定します。

# データ数を指定
N <- 10000

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

# 1フレーム当たりのデータ数を計算
n_per_frame <- N %/% frame_num
n_per_frame
## [1] 100

 データ数Nとフレーム数frame_numを指定して、1フレーム当たりのデータ数n_per_frameを計算します。

 サンプルを生成して、フレームごとにn_per_frame個ずつサンプルが増えるデータフレームを作成します。

# 多次元ガウス分布に従う乱数を生成
x_nd <- mvnfast::rmvn(n = N, mu = mu_d, sigma = sigma_dd)

# サンプルを複製して格納
anime_freq_df <- tidyr::expand_grid(
  frame = 1:frame_num, # フレーム番号
  n = 1:N # データ番号
) |> # 全ての組み合わせを作成
  dplyr::filter(n <= frame*n_per_frame) |> # 各フレームで利用するサンプルを抽出
  dplyr::mutate(
    x_1 = x_nd[n, 1], 
    x_2 = x_nd[n, 2], 
    parameter = paste0("mu=(", paste0(mu_d, collapse = ", "), "), Sigma=(", paste0(sigma_dd, collapse = ", "), "), N=", frame*n_per_frame) |> 
      factor(levels = paste0("mu=(", paste0(mu_d, collapse = ", "), "), Sigma=(", paste0(sigma_dd, collapse = ", "), "), N=", 1:frame_num*n_per_frame)) # フレーム切替用ラベル
  )
anime_freq_df
## # A tibble: 505,000 × 5
##    frame     n   x_1   x_2 parameter                                  
##    <int> <int> <dbl> <dbl> <fct>                                      
##  1     1     1  5.19  8.36 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=100
##  2     1     2  6.16 12.3  mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=100
##  3     1     3  6.19 10.2  mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=100
##  4     1     4  7.74 11.3  mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=100
##  5     1     5  4.30  8.45 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=100
##  6     1     6  6.15  7.93 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=100
##  7     1     7  4.60  9.75 mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=100
##  8     1     8  6.61 10.6  mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=100
##  9     1     9  7.21 11.5  mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=100
## 10     1    10  6.28 10.7  mu=(6, 10), Sigma=(1, 0.6, 0.6, 1.5), N=100
## # … with 504,990 more rows

 1データずつのときと同様に、作図用のデータを作成します。こちらは、フレームごとに「フレーム番号」掛ける「1フレーム当たりのデータ数」番目までのサンプルを抽出します。

 散布図のアニメーションを作成します。

# 散布図のアニメーションを作図
anime_freq_graph <- ggplot() + 
  geom_point(data = anime_freq_df, mapping = aes(x = x_1, y = x_2), 
             color = "orange", alpha = 0.3, size = 3) + # サンプル
  geom_contour(data = dens_df, mapping = aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 元の分布
  gganimate::transition_manual(parameter) + # フレーム
  labs(title ="Maltivariate Gaussian Distribution", 
       subtitle = "{current_frame}", 
       color = "density", 
       x = expression(x[1]), y = expression(x[2]))

# gif画像を作成
gganimate::animate(anime_freq_graph, nframes = frame_num+10, end_pause = 10, fps = 10, width = 800, height = 600)

 1データずつのときと同様に処理します。

 サンプルの度数のヒートマップのアニメーションを作成します。

# メッシュ用の値を設定
x_1_breaks <- seq(from = min(x_1_vals), to = max(x_1_vals), length.out = 30)
x_2_breaks <- seq(from = min(x_2_vals), to = max(x_2_vals), length.out = 30)

# 度数のヒートマップのアニメーションを作図
anime_freq_graph <- ggplot() + 
  geom_bin_2d(data = anime_freq_df, mapping = aes(x = x_1, y = x_2, fill = ..count..), 
              breaks = list(x = x_1_breaks, y = x_2_breaks), 
              alpha = 0.8) + # サンプル
  geom_contour(data = dens_df, mapping = aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 元の分布
  gganimate::transition_manual(parameter) + # フレーム
  scale_fill_distiller(palette = "Spectral") + # 塗りつぶしの色
  scale_color_distiller(palette = "Spectral") + # 等高線の色
  labs(title = "Maltivariate Gaussian Distribution", 
       subtitle = "{current_frame}", 
       fill = "frequency", color = "density", 
       x = expression(x[1]), y = expression(x[2]))

# gif画像を作成
gganimate::animate(anime_freq_graph, nframes = frame_num+10, end_pause = 10, fps = 10, width = 800, height = 600)


 サンプルの密度の等高線のアニメーションを作成します。

# 密度の等高線図のアニメーションを作図
anime_freq_graph <- ggplot() + 
  geom_density_2d_filled(data = anime_freq_df, mapping = aes(x = x_1, y = x_2, fill = ..level..), 
                         alpha = 0.8) + # サンプル
  geom_contour(data = dens_df, mapping = aes(x = x_1, y = x_2, z = density, color = ..level..)) + # 元の分布
  gganimate::transition_manual(parameter) + # フレーム
  scale_color_viridis_c(option = "D") + # 等高線の色
  labs(title = "Maltivariate Gaussian Distribution", 
       subtitle = "{current_frame}", 
       fill = "density", color = "density", 
       x = expression(x[1]), y = expression(x[2]))

# gif画像を作成
gganimate::animate(anime_freq_graph, nframes = frame_num+10, end_pause = 10, fps = 10, width = 800, height = 600)


2次元ガウス分布の乱数の散布図

2次元ガウス分布の乱数の度数のヒートマップ

2次元ガウス分布の乱数の密度の等高線図

 サンプルが増えるに従って、元の分布に近付くのを確認できます。

 この記事では、ガウス分布の乱数生成を確認しました。次は、ガウス分布から確率分布を生成します。

参考文献

  • C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上』,丸善出版,2012年.

おわりに

 色々寄り道したけど本筋に戻ってこれました。

【次の内容】

www.anarchive-beta.com