からっぽのしょこ

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

【R】ウィシャート分布の可視化

はじめに

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

 この記事では、R言語でウィシャート分布の乱数を用いたマハラノビス距離の等高線とその軸により可視化します。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

ウィシャート分布の可視化

 ウィシャート分布(Wishart Distribution)の確率密度をグラフで表現するのは難しいので、期待値と乱数を用いたマハラノビス距離(Mahalanobis Distance)のグラフにより可視化します。ウィシャート分布については「ウィシャート分布の定義式 - からっぽのしょこ」、マハラノビス距離については「分散共分散行列とユークリッド距離・マハラノビス距離の関係の導出 - からっぽのしょこ」、固有ベクトルについては「2.3.0:分散共分散行列と固有値・固有ベクトルの関係の導出【PRMLのノート】 - からっぽのしょこ」、楕円の軸については「【R】2次元ガウス分布の作図 - からっぽのしょこ」を参照してください。

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

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

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

期待値による距離と軸の計算

 まずは、ウィシャート分布のパラメータを設定して、分散共分散行列(確率変数の逆行列)の期待値によるマハラノビス距離の楕円と軸を計算します。

 ウィシャート分布のパラメータ$\nu, \mathbf{W}$とサンプルサイズ$N$を設定します。この例では、2次元のグラフで描画するため、次元数を$D = 2$とします。

# 次元数を指定
D <- 2

# 自由度を指定
nu <- D + 8

# 逆スケール行列を指定
w_dd <- matrix(c(5, 0.6, 0.6, 12), nrow = D, ncol = D)

# パラメータの数(サンプルサイズ)を指定
N <- 100

 自由度$\nu$と逆スケール行列$\mathbf{W}$、サンプルサイズ$N$を指定します。
 自由度は$\nu > D - 1$、逆スケール行列は$D \times D$の正定値行列を満たす必要があります。

 分散共分散行列の期待値$\mathbb{E}[\boldsymbol{\Sigma}]$を計算します。

# 分散共分散行列の期待値を計算
E_sigma_dd <- solve(nu * w_dd)
E_sigma_dd
##              [,1]         [,2]
## [1,]  0.020120724 -0.001006036
## [2,] -0.001006036  0.008383635

 ウィシャート分布の期待値(精度行列$\boldsymbol{\Lambda}$の期待値)$\mathbb{E}[\boldsymbol{\Lambda}] = \nu \mathbf{W}$の逆行列をとり、分散共分散行列の期待値$\mathbb{E}[\boldsymbol{\Sigma}] = (\nu \mathbf{W})^{-1}$を求めます。

 分散共分散行列に対応する変数の値$\mathbf{x}$を作成します。

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

# xの点を作成
x_mat <- tidyr::expand_grid(
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # 格子点を作成
  as.matrix() # マトリクスに変換
head(x_mat)
##             x_1        x_2
## [1,] -0.4255426 -0.2746866
## [2,] -0.4255426 -0.2691374
## [3,] -0.4255426 -0.2635881
## [4,] -0.4255426 -0.2580389
## [5,] -0.4255426 -0.2524897
## [6,] -0.4255426 -0.2469405

 $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)$に対応します。

 分散共分散行列の期待値$\mathbb{E}[\boldsymbol{\Sigma}]$によるマハラノビス距離を計算します。

# 期待値によるマハラノビス距離(楕円)を計算
E_delta_df <- tibble::tibble(
  x_1 = x_mat[, 1], # x軸の値
  x_2 = x_mat[, 2], # y軸の値
  delta = mahalanobis(x = x_mat, center = 0, cov = E_sigma_dd) |> 
    sqrt() # マハラノビス距離
)
E_delta_df
## # A tibble: 10,000 × 3
##       x_1    x_2 delta
##     <dbl>  <dbl> <dbl>
##  1 -0.426 -0.275  4.42
##  2 -0.426 -0.269  4.37
##  3 -0.426 -0.264  4.33
##  4 -0.426 -0.258  4.29
##  5 -0.426 -0.252  4.24
##  6 -0.426 -0.247  4.20
##  7 -0.426 -0.241  4.16
##  8 -0.426 -0.236  4.12
##  9 -0.426 -0.230  4.07
## 10 -0.426 -0.225  4.03
## # … with 9,990 more rows

 マハラノビス距離は、mahalanobis()で計算できます。変数の引数Xx_mat、平均ベクトルの引数center0(またはc(0, 0))、分散共分散行列の引数covE_sigma_ddを指定します。

 分散共分散行列の期待値$\mathbb{E}[\boldsymbol{\Sigma}]$の固有値と固有ベクトルを計算します。

# 期待値の固有値・固有ベクトルを計算
res_eigen  <- eigen(E_sigma_dd)
E_lambda_d <- res_eigen[["values"]]
E_u_dd     <- res_eigen[["vectors"]] |> 
  t()
E_lambda_d; E_u_dd
## [1] 0.020206332 0.008298028
##             [,1]        [,2]
## [1,] -0.99639908  0.08478722
## [2,] -0.08478722 -0.99639908

 固有値$\lambda_1, \lambda_2$と固有ベクトル$\mathbf{u}_1 = (u_{1,1}, u_{1,2}), \mathbf{u}_2 = (u_{2,1}, u_{2,2})$をeigen()で計算します。リストが出力されるので、"values"で固有値(をまとめたベクトル)、"vectors"で固有ベクトル(をまとめたマトリクス)を取り出します。数式での成分と合わせるために転置しておきます。

 固有値と固有ベクトルを用いて、楕円の軸を計算します。

# 期待値による楕円の長軸を計算
E_axis_df <- tibble::tibble(
  xend = E_u_dd[1, 1] * sqrt(E_lambda_d[1]), 
  yend = E_u_dd[1, 2] * sqrt(E_lambda_d[1])
)
E_axis_df
## # A tibble: 1 × 2
##     xend   yend
##    <dbl>  <dbl>
## 1 -0.142 0.0121

 軸番号を$i$、次元を$j$として、軸の先端の点を$u_{j,i} \sqrt{\lambda_j}$で計算します。ここでは、長い方の軸($i = 1$)のデータフレームに格納します。

 ここまでは、ウィシャート分布の期待値に用いて計算しました。続いて、ウィシャート分布のサンプルに用いて計算します。

マハラノビス距離(楕円)の作図

 次は、ウィシャート分布から精度行列をサンプリングして、サンプルごとにマハラノビス距離の楕円と軸を計算します。乱数生成については「【R】ウィシャート分布から確率分布の生成 - からっぽのしょこ」を参照してください。

 ウィシャート分布に従う乱数を生成します。

# 分散共分散行列を生成
lambda_ddn <- rWishart(n = N, df = nu, Sigma = w_dd)
lambda_ddn[, , 1:3]
## , , 1
## 
##           [,1]      [,2]
## [1,]  35.01167 -23.27527
## [2,] -23.27527  64.35598
## 
## , , 2
## 
##           [,1]       [,2]
## [1,] 66.733981   7.994431
## [2,]  7.994431 115.121494
## 
## , , 3
## 
##           [,1]       [,2]
## [1,] 81.253923   6.525797
## [2,]  6.525797 232.926304

 ウィシャート分布の乱数は、rWishart()で生成できます。データ数の引数nN、自由度の引数dfnu、逆スケール行列の引数Sigmaw_ddを指定します。
 生成した値をガウス分布の精度行列$\boldsymbol{\Lambda}_n = (\lambda_{1,1,n}, \cdots, \lambda_{D,D, n})$のN個のサンプルlambda_ddnとします。(記号が被っていますが、$\lambda_{i,j,n}$はウィシャート分布の乱数(精度行列)の各成分です。)

 精度行列のサンプル$\boldsymbol{\Lambda}_n$ごとにマハラノビス距離を計算します。

# サンプルごとにマハラノビス距離(楕円)を計算
res_delta_df <- tidyr::expand_grid(
  n = 1:N, # データ番号
  x_1 = x_1_vals, 
  x_2 = x_2_vals 
) |> # パラメータごとに格子点を複製
  dplyr::group_by(n) |> # 距離の計算用にグループ化
  dplyr::mutate(
    delta = mahalanobis(x = x_mat, center = 0, cov = lambda_ddn[, , unique(n)], inverted = TRUE) |> 
      sqrt() # マハラノビス距離
  ) |> 
  dplyr::ungroup() # グループ化を解除
res_delta_df
## # A tibble: 1,000,000 × 4
##        n    x_1    x_2 delta
##    <int>  <dbl>  <dbl> <dbl>
##  1     1 -0.426 -0.275  2.40
##  2     1 -0.426 -0.269  2.38
##  3     1 -0.426 -0.264  2.36
##  4     1 -0.426 -0.258  2.35
##  5     1 -0.426 -0.252  2.33
##  6     1 -0.426 -0.247  2.32
##  7     1 -0.426 -0.241  2.30
##  8     1 -0.426 -0.236  2.29
##  9     1 -0.426 -0.230  2.28
## 10     1 -0.426 -0.225  2.27
## # … with 999,990 more rows

 パラメータ番号(1からNの整数)と変数x_1_vals, x_2_valsの要素の全ての組み合わせをexpand_grid()を作成します。これにより、サンプルごとにx_matを複製できます。
 パラメータ列nでグループ化することで、x_matごとにマハラノビス距離を計算できます。
 mahalanobis()x引数にx_matcenter引数に0cov引数にlambda_ddnn番目の値(マトリクス)を指定します。cov引数に精度行列を指定する場合は、inverted引数にTRUEを指定します。

 精度行列のサンプル$\boldsymbol{\Lambda}_n$ごとのマハラノビス距離の等高線図を作成します。

# 等高線を引く距離を指定
delta <- 1

# サンプルごとにマハラノビス距離(楕円)を作図
ggplot() + 
  geom_contour(data = E_delta_df, mapping = aes(x = x_1, y = x_2, z = delta), 
               breaks = delta, color ="red", size = 1, linetype = "dashed") + # 期待値による分布
  geom_contour(data = res_delta_df, mapping = aes(x = x_1, y = x_2, z = delta, color = factor(n)), 
               breaks = delta, alpha = 0.5, show.legend = FALSE) + # サンプルによる分布
  coord_fixed(ratio = 1, xlim = c(min(x_1_vals), max(x_1_vals)), ylim = c(min(x_2_vals), max(x_2_vals))) + # アスペクト比
  labs(title = "Mahalabobis Distance", 
       subtitle = parse(text = paste0("list(nu==", nu, ", W==(list(", paste0(round(w_dd, 2), collapse = ", "), ")), N==", N, ")")), 
       x = expression(x[1]), y = expression(x[2]))

ウィシャート分布の乱数を用いたマハラノビス距離の等高線図

 geom_contour()breaks引数に値を指定すると、その値で等高線を引きます。
 期待値による楕円を赤色の破線で示します。

 それぞれ楕円形になります。

 精度行列のサンプル$\boldsymbol{\Lambda}_n$ごとに固有値と固有ベクトルを計算して、楕円の軸を計算します。

# サンプルごとに楕円の長軸を計算
res_axis_df <- tibble::tibble(
  n = 1:N # データ番号
) |> 
  dplyr::group_by(n) |> # 軸の計算用にグループ化
  dplyr::mutate(
    lambda_1 = lambda_ddn[, , n] |> 
      solve() |> 
      eigen() |> 
      (\(.){.[["values"]][1]})(), # 固有値
    u_11 = lambda_ddn[, , n] |> 
      solve() |> 
      eigen() |> 
      (\(.){t(.[["vectors"]])[1, 1]})(), # 固有ベクトルの成分
    u_12 = lambda_ddn[, , n] |> 
      solve() |> 
      eigen() |> 
      (\(.){t(.[["vectors"]])[1, 2]})(), # 固有ベクトルの成分
    xend = u_11 * sqrt(lambda_1), # 終点のx軸の値
    yend = u_12 * sqrt(lambda_1)  # 終点のy軸の値
  ) |> 
  dplyr::ungroup() # グループ化を解除
res_axis_df
## # A tibble: 100 × 6
##        n lambda_1   u_11     u_12   xend      yend
##    <int>    <dbl>  <dbl>    <dbl>  <dbl>     <dbl>
##  1     1   0.0451 -0.876 -0.483   -0.186 -0.103   
##  2     2   0.0153 -0.987  0.159   -0.122  0.0196  
##  3     3   0.0123 -0.999  0.0429  -0.111  0.00477 
##  4     4   0.0292  0.621  0.784    0.106  0.134   
##  5     5   0.0211 -0.938  0.345   -0.136  0.0501  
##  6     6   0.0161 -0.989 -0.146   -0.125 -0.0185  
##  7     7   0.0868 -0.999 -0.0374  -0.294 -0.0110  
##  8     8   0.0255 -0.793  0.609   -0.127  0.0973  
##  9     9   0.0318 -0.975  0.223   -0.174  0.0397  
## 10    10   0.0190 -1.00  -0.00472 -0.138 -0.000650
## # … with 90 more rows

 上手いことして(もっといいやり方があれば教えてください)、固有値と固有ベクトルの値の列を作成して、軸の先端の点を$u_{j,i}^{(n)} \sqrt{\lambda_j^{(n)}}$で計算します。($\lambda_j^{(n)}$はウィシャート分布の$n$番目の乱数$\boldsymbol{\Lambda}_n$の逆行列$\boldsymbol{\Sigma}_n$の固有値です。)

 精度行列のサンプル$\boldsymbol{\Lambda}_n$ごとの楕円の軸を描画します。

# サンプルごとに楕円の長軸を作図
ggplot() + 
  geom_contour(data = E_delta_df, mapping = aes(x = x_1, y = x_2, z = delta), 
               breaks = delta, color ="red", size = 1, linetype = "dashed") + # 期待値による分布
  geom_segment(data = E_axis_df, mapping = aes(x = 0, y = 0, xend = xend, yend = yend), 
               color = "red", size = 1, linetype = "dashed", arrow = arrow(length = unit(10, "pt"))) + # 期待値による分布の長軸
  geom_segment(data = res_axis_df, mapping = aes(x = 0, y = 0, xend = xend, yend = yend), 
               color = "orange", alpha = 0.5, arrow = arrow(length = unit(5, "pt"))) + # サンプルによる分布の長軸
  coord_fixed(ratio = 1, xlim = c(min(x_1_vals), max(x_1_vals)), ylim = c(min(x_2_vals), max(x_2_vals))) + # アスペクト比
  labs(title = "Mahalanobis Distance", 
       subtitle = parse(text = paste0("list(nu==", nu, ", W==(list(", paste0(round(w_dd, 2), collapse = ", "), ")), N==", N, ")")), 
       x = expression(x[1]), y = expression(x[2]))

ウィシャート分布の乱数を用いたマハラノビス距離の等高線(楕円)の長軸

 geom_segment()で楕円の軸を描画します。軸の始点の引数x, yに原点0、終点の引数xend, yendaxis_dfxstart, ystart列を指定すると、長さが$\sqrt{\lambda_j^{(n)}}$のベクトルを描画できます。

 サンプルの軸(オレンジ色の矢印)が期待値による軸(赤色の矢印)に近い方向を向いているのを確認できます。

アニメーションによる可視化

 次は、パラメータの値を少しずつ変化させて、期待値とサンプルによるマハラノビス距離の楕円と軸の変化をアニメーションで確認します。

自由度の影響

 まずは、自由度$\nu$の値を変化させ、$\mathbf{W}$を固定します。

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

# 次元数を指定
D <- 2

# 自由度として利用する値を指定
nu_vals <- seq(from = D, to = 50, by = 1)

# 逆スケール行列を指定
w_dd <- matrix(c(6, -0.6, -0.6, 10), nrow = D, ncol = D)

# フレーム数を設定
frame_num <- length(nu_vals)
frame_num
## [1] 49

 値の間隔が一定になるように$\nu$の値をnu_valsとして作成します。パラメータごとにフレームを切り替えるので、nu_valsの要素数がアニメーションのフレーム数になります。

 設定したパラメータに応じて、変数の値$\mathbf{x}$を作成します。

# 作図用の共分散行列を作成
E_sigma_dd <- solve(max(nu_vals) * w_dd)

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

# xの点を作成
x_mat <- tidyr::expand_grid(
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # 格子点を作成
  as.matrix() # マトリクスに変換
head(x_mat)
##            x_1        x_2
## [1,] -0.173727 -0.1345684
## [2,] -0.173727 -0.1318498
## [3,] -0.173727 -0.1291313
## [4,] -0.173727 -0.1264127
## [5,] -0.173727 -0.1236942
## [6,] -0.173727 -0.1209756

 ここで計算するE_sigma_ddは、描画範囲を一定にするための仮の値です。

 パラメータごとに分散共分散行列の期待値$\mathbb{E}[\boldsymbol{\Sigma}]$を計算して、マハラノビス距離を計算します。

# パラメータごとに共分散行列の期待値によるマハラノビス距離(楕円)を計算
anime_E_delta_df <- tidyr::expand_grid(
  nu = nu_vals, # パラメータ
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # パラメータごとに格子点を複製
  dplyr::group_by(nu) |> # 距離の計算用にグループ化
  dplyr::mutate(
    delta = mahalanobis(
      x = x_mat, 
      center = 0, 
      cov = unique(nu) * w_dd, # 精度行列の期待値を指定
      inverted = TRUE
    ), # 確率密度
    parameter = paste0("nu=", unique(nu), ", W=(", paste0(w_dd, collapse = ", "), ")") |> 
      factor(levels = paste0("nu=", nu_vals, ", W=(", paste0(w_dd, collapse = ", "), ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_E_delta_df
## # A tibble: 490,000 × 5
##       nu    x_1    x_2 delta parameter                  
##    <dbl>  <dbl>  <dbl> <dbl> <fct>                      
##  1     2 -0.174 -0.135 0.668 nu=2, W=(6, -0.6, -0.6, 10)
##  2     2 -0.174 -0.132 0.655 nu=2, W=(6, -0.6, -0.6, 10)
##  3     2 -0.174 -0.129 0.642 nu=2, W=(6, -0.6, -0.6, 10)
##  4     2 -0.174 -0.126 0.629 nu=2, W=(6, -0.6, -0.6, 10)
##  5     2 -0.174 -0.124 0.617 nu=2, W=(6, -0.6, -0.6, 10)
##  6     2 -0.174 -0.121 0.604 nu=2, W=(6, -0.6, -0.6, 10)
##  7     2 -0.174 -0.118 0.593 nu=2, W=(6, -0.6, -0.6, 10)
##  8     2 -0.174 -0.116 0.581 nu=2, W=(6, -0.6, -0.6, 10)
##  9     2 -0.174 -0.113 0.570 nu=2, W=(6, -0.6, -0.6, 10)
## 10     2 -0.174 -0.110 0.559 nu=2, W=(6, -0.6, -0.6, 10)
## # … with 489,990 more rows

 パラメータnu_valsと確率変数x_1_vals, x_2_valsの要素の全ての組み合わせをexpand_grid()を作成します。これにより、パラメータごとにx_matを複製できます。
 パラメータ列nuでグループ化することで、x_matごとに距離を計算できます。
 パラメータごとにフレーム切替用のラベルを作成します。文字列型だと文字列の基準で順序が決まるので、因子型にしてパラメータに応じたレベル(順序)を設定します。

 パラメータごとに分散共分散行列の期待値$\mathbb{E}[\boldsymbol{\Sigma}]$の固有値・固有ベクトルを計算して、楕円の軸を計算します。

# パラメータごとに共分散行列の期待値による楕円の軸を計算
anime_E_axis_df <- tidyr::expand_grid(
  nu = nu_vals, # パラメータ
  name = paste0(rep(c("x", "y"), each = 2, times = 2), rep(c("start", "end"), each = 4)) # 列名
) |> # パラメータごとに受け皿を複製
  dplyr::group_by(nu) |> # 軸の計算用にグループ化
  dplyr::mutate(
    axis = rep(c("y_1", "y_2"), times = 4), # pivot_wider用の列
    sign = rep(c(-1, 1), each = 4), # 始点・終点の計算用の符号
    lambda = solve(unique(nu) * w_dd) |> # 共分散行列の期待値を計算
      (\(.){eigen(.)[["values"]]})() |> 
      (\(.){rep(sqrt(.), times = 4)})(), # 固有値
    u = solve(unique(nu) * w_dd) |> # 共分散行列の期待値を計算
      (\(.){t(eigen(.)[["vectors"]])})() |> 
      (\(.){rep(as.vector(.), times = 2)})(), # 固有ベクトル
    value = sign * u * lambda, # 軸の始点・終点を計算
    parameter = paste0("nu=", nu, ", W=(", paste0(w_dd, collapse = ", "), ")") |> 
      factor(levels = paste0("nu=", nu_vals, ", W=(", paste0(w_dd, collapse = ", "), ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() |> # グループ化を解除
  dplyr::select(axis, name, value, parameter) |> # 利用する列を選択
  tidyr::pivot_wider(
    id_cols = c(axis, parameter), 
    names_from = name, 
    values_from = value
  ) # 軸の視点・終点の列を分割
anime_E_axis_df
## # A tibble: 98 × 6
##    axis  parameter                    xstart ystart    xend    yend
##    <chr> <fct>                         <dbl>  <dbl>   <dbl>   <dbl>
##  1 y_1   nu=2, W=(6, -0.6, -0.6, 10)  0.288  0.0422 -0.288  -0.0422
##  2 y_2   nu=2, W=(6, -0.6, -0.6, 10) -0.0323 0.220   0.0323 -0.220 
##  3 y_1   nu=3, W=(6, -0.6, -0.6, 10)  0.235  0.0345 -0.235  -0.0345
##  4 y_2   nu=3, W=(6, -0.6, -0.6, 10) -0.0264 0.180   0.0264 -0.180 
##  5 y_1   nu=4, W=(6, -0.6, -0.6, 10)  0.203  0.0299 -0.203  -0.0299
##  6 y_2   nu=4, W=(6, -0.6, -0.6, 10) -0.0229 0.156   0.0229 -0.156 
##  7 y_1   nu=5, W=(6, -0.6, -0.6, 10)  0.182  0.0267 -0.182  -0.0267
##  8 y_2   nu=5, W=(6, -0.6, -0.6, 10) -0.0204 0.139   0.0204 -0.139 
##  9 y_1   nu=6, W=(6, -0.6, -0.6, 10)  0.166  0.0244 -0.166  -0.0244
## 10 y_2   nu=6, W=(6, -0.6, -0.6, 10) -0.0187 0.127   0.0187 -0.127 
## # … with 88 more rows

 パラメータごとに2つの軸の始点と終点の8つの値を計算します。
 パラメータsigma2_1_valsと8つの値に対応する列名の全ての組み合わせをexpand_grid()で作成します。これにより、パラメータごとに列名を複製できます。
 eigen()で固有値と固有ベクトルをそれぞれ計算し、上手いことして軸の始点と終点を計算します(誰かもっと上手く処理してください…)。
 求めた値をpivot_wider()で列名と対応付けて分割します。

 パラメータごとにウィシャート分布の乱数(精度行列$\boldsymbol{\Lambda}_n$)を生成して、サンプルごとに距離を計算します。

# 楕円の数(サンプルサイズ)を指定
N <- 30

# パラメータごとに楕円のサンプルを計算
anime_sample_delta_df <- tidyr::expand_grid(
  nu = nu_vals, # パラメータ
  n = 1:N, # サンプル番号
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # パラメータとサンプルごとに格子点を複製
  dplyr::group_by(nu, n) |> # 距離の計算用にグループ化
  dplyr::mutate(
    delta = mahalanobis(
      x = x_mat, 
      center = 0, 
      cov = rWishart(n = 1, df = nu, Sigma = w_dd)[, , 1], # 精度行列を生成
      inverted = TRUE
    ), 
    parameter = paste0("nu=", nu, ", W=(", paste0(w_dd, collapse = ", "), ")") |> 
      factor(levels = paste0("nu=", nu_vals, ", W=(", paste0(w_dd, collapse = ", "), ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_sample_delta_df
## # A tibble: 14,700,000 × 6
##       nu     n    x_1    x_2 delta parameter                  
##    <dbl> <int>  <dbl>  <dbl> <dbl> <fct>                      
##  1     2     1 -0.174 -0.135 0.189 nu=2, W=(6, -0.6, -0.6, 10)
##  2     2     1 -0.174 -0.132 0.188 nu=2, W=(6, -0.6, -0.6, 10)
##  3     2     1 -0.174 -0.129 0.188 nu=2, W=(6, -0.6, -0.6, 10)
##  4     2     1 -0.174 -0.126 0.187 nu=2, W=(6, -0.6, -0.6, 10)
##  5     2     1 -0.174 -0.124 0.186 nu=2, W=(6, -0.6, -0.6, 10)
##  6     2     1 -0.174 -0.121 0.186 nu=2, W=(6, -0.6, -0.6, 10)
##  7     2     1 -0.174 -0.118 0.185 nu=2, W=(6, -0.6, -0.6, 10)
##  8     2     1 -0.174 -0.116 0.185 nu=2, W=(6, -0.6, -0.6, 10)
##  9     2     1 -0.174 -0.113 0.184 nu=2, W=(6, -0.6, -0.6, 10)
## 10     2     1 -0.174 -0.110 0.184 nu=2, W=(6, -0.6, -0.6, 10)
## # … with 14,699,990 more rows

 パラメータnu_vals・サンプル番号(1からNの整数)・確率変数x_1_vals, x_2_valsの要素の全ての組み合わせを作成します。これにより、パラメータとサンプルの組み合わせごとにx_matを複製します。
 x_matごとにサンプルを生成してマハラノビス距離を計算します。

 パラメータごとにウィシャート分布の乱数(精度行列$\boldsymbol{\Lambda}_n$)を生成して、サンプルごとに軸を計算します。

# 軸の数(サンプルサイズ)を指定
N <- 100

# パラメータごとに楕円の長軸のサンプルを計算
anime_sample_axis_df <- tidyr::expand_grid(
  nu = nu_vals, # パラメータ
  n = 1:N # サンプル番号
) |> # パラメータごとに受け皿を複製
  dplyr::group_by(nu, n) |> # 軸の計算用にグループ化
  dplyr::mutate(
    eigen_lt = rWishart(n = 1, df = nu, Sigma = w_dd)[, , 1] |> # 精度行列を生成
      solve() |> # 分散共分散行列に変換
      eigen() |> # 固有値・固有ベクトルを計算
      list(), # リストに格納
    lambda_1 = eigen_lt[[1]][["values"]][1], # 固有値
    u_11 = eigen_lt[[1]][["vectors"]][1, 1], # 固有ベクトルの成分
    u_12 = eigen_lt[[1]][["vectors"]][2, 1], # 固有ベクトルの成分
    xend = u_11 * sqrt(lambda_1), # 終点のx軸の値
    yend = u_12 * sqrt(lambda_1), # 終点のy軸の値
    parameter = paste0("nu=", nu, ", W=(", paste0(w_dd, collapse = ", "), ")") |> 
      factor(levels = paste0("nu=", nu_vals, ", W=(", paste0(w_dd, collapse = ", "), ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_sample_axis_df
## # A tibble: 4,900 × 9
##       nu     n eigen_lt lambda_1   u_11   u_12   xend    yend parameter         
##    <dbl> <int> <list>      <dbl>  <dbl>  <dbl>  <dbl>   <dbl> <fct>             
##  1     2     1 <eigen>    0.479  -0.993  0.115 -0.687  0.0794 nu=2, W=(6, -0.6,…
##  2     2     2 <eigen>    0.336  -0.945 -0.326 -0.548 -0.189  nu=2, W=(6, -0.6,…
##  3     2     3 <eigen>    3.90   -0.151  0.989 -0.297  1.95   nu=2, W=(6, -0.6,…
##  4     2     4 <eigen>    0.149  -0.809 -0.588 -0.312 -0.227  nu=2, W=(6, -0.6,…
##  5     2     5 <eigen>    0.326  -0.860 -0.510 -0.491 -0.291  nu=2, W=(6, -0.6,…
##  6     2     6 <eigen>    0.0527 -0.786  0.618 -0.180  0.142  nu=2, W=(6, -0.6,…
##  7     2     7 <eigen>    6.43    0.661  0.751  1.68   1.90   nu=2, W=(6, -0.6,…
##  8     2     8 <eigen>    0.450  -0.567  0.824 -0.381  0.553  nu=2, W=(6, -0.6,…
##  9     2     9 <eigen>    4.06   -0.837  0.547 -1.69   1.10   nu=2, W=(6, -0.6,…
## 10     2    10 <eigen>    0.0929 -0.714  0.700 -0.218  0.213  nu=2, W=(6, -0.6,…
## # … with 4,890 more rows

 パラメータnu_valsとサンプル番号(1からNの整数)の要素の全ての組み合わせを作成します。これにより、パラメータとごとにN個のサンプルを格納するための行を複製します。
 パラメータとサンプル番号の組み合わせごとにサンプルを生成して、固有値と固有ベクトルを計算し、楕円の軸の先端の値を計算します。

 パラメータごとに分散共分散行列の期待値$\mathbb{E}[\boldsymbol{\Sigma}]$を計算して、ラベルを作成します。

# パラメータごとに精度行列の期待値のラベルを作成
anime_label_df <- tibble::tibble(
  nu = nu_vals # パラメータ
) |> 
  dplyr::group_by(nu) |> # 期待値の計算用にグループ化
  dplyr::mutate(
    label = paste0("E[Lambda]=(", paste0(round(unique(nu) * w_dd, 2), collapse = ", "), ")"), # 精度行列ラベル
    parameter = paste0("nu=", nu, ", W=(", paste0(w_dd, collapse = ", "), ")") |> 
      factor(levels = paste0("nu=", nu_vals, ", W=(", paste0(w_dd, collapse = ", "), ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_label_df
## # A tibble: 49 × 3
##       nu label                           parameter                   
##    <dbl> <chr>                           <fct>                       
##  1     2 E[Lambda]=(12, -1.2, -1.2, 20)  nu=2, W=(6, -0.6, -0.6, 10) 
##  2     3 E[Lambda]=(18, -1.8, -1.8, 30)  nu=3, W=(6, -0.6, -0.6, 10) 
##  3     4 E[Lambda]=(24, -2.4, -2.4, 40)  nu=4, W=(6, -0.6, -0.6, 10) 
##  4     5 E[Lambda]=(30, -3, -3, 50)      nu=5, W=(6, -0.6, -0.6, 10) 
##  5     6 E[Lambda]=(36, -3.6, -3.6, 60)  nu=6, W=(6, -0.6, -0.6, 10) 
##  6     7 E[Lambda]=(42, -4.2, -4.2, 70)  nu=7, W=(6, -0.6, -0.6, 10) 
##  7     8 E[Lambda]=(48, -4.8, -4.8, 80)  nu=8, W=(6, -0.6, -0.6, 10) 
##  8     9 E[Lambda]=(54, -5.4, -5.4, 90)  nu=9, W=(6, -0.6, -0.6, 10) 
##  9    10 E[Lambda]=(60, -6, -6, 100)     nu=10, W=(6, -0.6, -0.6, 10)
## 10    11 E[Lambda]=(66, -6.6, -6.6, 110) nu=11, W=(6, -0.6, -0.6, 10)
## # … with 39 more rows

 パラメータごとにラベル用の文字列を作成します。

 等高線図のアニメーション(gif画像)を作成します。

# 等高線を引く距離を指定
delta <- 1

# マハラノビス距離(楕円)のアニメーションを作図
anime_delta_graph <- ggplot() + 
  geom_contour(data = anime_E_delta_df, mapping = aes(x = x_1, y = x_2, z = delta), 
               breaks = delta, color ="red", size = 1, linetype = "dashed") + # 期待値による分布
  geom_contour(data = anime_sample_delta_df, mapping = aes(x = x_1, y = x_2, z = delta, color = factor(n)), 
               breaks = delta, alpha = 0.5, show.legend = FALSE) + # サンプルによる分布
  geom_label(data = anime_label_df, mapping = aes(x = min(x_1_vals), y = max(x_2_vals), label = label), 
             hjust = 0, vjust = 0, alpha = 0.5) + # 分散共分散行列の期待値ラベル
  gganimate::transition_manual(parameter) + # フレーム
  coord_fixed(ratio = 1, xlim = c(min(x_1_vals), max(x_1_vals)), ylim = c(min(x_2_vals), max(x_2_vals))) + # アスペクト比
  labs(title = "Mahalanobis Distance", 
       subtitle = "{current_frame}", 
       x = expression(x[1]), y = expression(x[2]))

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

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

 同様に、マハラノビス距離(楕円)の軸のアニメーションを作成します。

# 楕円の軸のアニメーションを作図
anime_axis_graph <- ggplot() + 
  geom_contour_filled(data = anime_E_delta_df, mapping = aes(x = x_1, y = x_2, z = delta, fill = ..level..), 
                      alpha = 0.8) + # 期待値による分布
  geom_contour(data = anime_E_delta_df, mapping = aes(x = x_1, y = x_2, z = delta), 
               breaks = 1, color = "red", size = 1, linetype = "dotted") + # 期待値による分布の断面
  geom_segment(data = anime_E_axis_df, mapping = aes(x = xstart, y = ystart, xend = xend, yend = yend), 
               color = "red", size = 1, linetype = "dashed", arrow = arrow(length = unit(10, "pt"))) + # 期待値による断面の軸
  geom_segment(data = anime_sample_axis_df, mapping = aes(x = 0, y = 0, xend = xend, yend = yend), 
               color = "orange", alpha = 0.5, size = 1, arrow = arrow(length = unit(5, "pt"))) + # サンプルによる分布の長軸
  geom_label(data = anime_label_df, mapping = aes(x = min(x_1_vals), y = max(x_2_vals), label = label), 
             hjust = 0, vjust = 0, alpha = 0.5) + # 分散共分散行列の期待値ラベル
  gganimate::transition_manual(parameter) + # フレーム
  coord_fixed(ratio = 1, xlim = c(min(x_1_vals), max(x_1_vals)), ylim = c(min(x_2_vals), max(x_2_vals))) + # アスペクト比
  labs(title ="Mahalanobis Distance", 
       subtitle = "{current_frame}", 
       fill = "distance", 
       x = expression(x[1]), y = expression(x[2]))

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


ウィシャート分布の自由度とマハラノビス距離とその軸の関係

 $\nu$が大きくなるほど、サンプルによる楕円と軸が期待値による楕円と軸に近付いています(なぜ?)。

スケール行列(1,1成分)の影響

 続いて、$\mathbf{W}$の$1,1$成分$w_{1,1}$の値を変化させ、$\nu, w_{1,2}, w_{2,2}$を固定します。

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

# 次元数を指定
D <- 2

# 自由度を指定
nu <- D + 18

# 逆スケール行列を指定
w_11_vals <- seq(from = 1, to = 10, by = 0.1) |> 
  round(digits = 2)
w_12 <- 0.6
w_22 <- 5

# フレーム数を設定
frame_num <- length(w_11_vals)
frame_num
## [1] 91

 値の間隔が一定になるように$w_{1,1}$の値をw_11_valsとして作成します。また、$w_{1,2} = w_{2,1}$をw_12、$w_{2,2}$の値をw_22として値を指定します。

 設定したパラメータに応じて、変数の値$\mathbf{x}$を作成します。

# 作図用の共分散行列を作成
E_sigma_dd <- solve(nu * matrix(c(min(w_11_vals), w_12, w_12, w_22), nrow = D, ncol = D))

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

# xの点を作成
x_mat <- tidyr::expand_grid(
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # 格子点を作成
  as.matrix() # マトリクスに変換
head(x_mat)
##             x_1        x_2
## [1,] -0.3481788 -0.3114205
## [2,] -0.3481788 -0.3051292
## [3,] -0.3481788 -0.2988379
## [4,] -0.3481788 -0.2925466
## [5,] -0.3481788 -0.2862553
## [6,] -0.3481788 -0.2799639


 パラメータごとに期待値による距離を計算します。

# パラメータごとに共分散行列の期待値によるマハラノビス距離(楕円)を計算
anime_E_delta_df <- tidyr::expand_grid(
  w_11 = w_11_vals, # パラメータ
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # サンプルごとに格子点を複製
  dplyr::group_by(w_11) |> # 距離の計算用にグループ化
  dplyr::mutate(
    delta = mahalanobis(
      x = x_mat, 
      center = 0, 
      cov = nu * matrix(c(unique(w_11), w_12, w_12, w_22), nrow = D, ncol = D), # 精度行列の期待値を指定
      inverted = TRUE
    ), # 確率密度
    parameter = paste0("nu=", nu, ", W=(", unique(w_11), ", ", w_12, ", ", w_12, ", ", w_22, ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11_vals, ", ", w_12, ", ", w_12, ", ", w_22, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_E_delta_df
## # A tibble: 910,000 × 5
##     w_11    x_1    x_2 delta parameter                
##    <dbl>  <dbl>  <dbl> <dbl> <fct>                    
##  1     1 -0.348 -0.311  14.7 nu=20, W=(1, 0.6, 0.6, 5)
##  2     1 -0.348 -0.305  14.3 nu=20, W=(1, 0.6, 0.6, 5)
##  3     1 -0.348 -0.299  13.9 nu=20, W=(1, 0.6, 0.6, 5)
##  4     1 -0.348 -0.293  13.4 nu=20, W=(1, 0.6, 0.6, 5)
##  5     1 -0.348 -0.286  13.0 nu=20, W=(1, 0.6, 0.6, 5)
##  6     1 -0.348 -0.280  12.6 nu=20, W=(1, 0.6, 0.6, 5)
##  7     1 -0.348 -0.274  12.2 nu=20, W=(1, 0.6, 0.6, 5)
##  8     1 -0.348 -0.267  11.8 nu=20, W=(1, 0.6, 0.6, 5)
##  9     1 -0.348 -0.261  11.4 nu=20, W=(1, 0.6, 0.6, 5)
## 10     1 -0.348 -0.255  11.0 nu=20, W=(1, 0.6, 0.6, 5)
## # … with 909,990 more rows

 「自由度の影響」のnu, nu_valsw_11, w_11_valsを置き換えるように処理します。

 パラメータごとに期待値による軸を計算します。

# パラメータごとに共分散行列の期待値による楕円の軸を計算
anime_E_axis_df <- tidyr::expand_grid(
  w_11 = w_11_vals, # パラメータごとに
  name = paste0(rep(c("x", "y"), each = 2, times = 2), rep(c("start", "end"), each = 4)) # 列名
) |> # パラメータごとに受け皿を複製
  dplyr::group_by(w_11) |> # 軸の計算用にグループ化
  dplyr::mutate(
    axis = rep(c("y_1", "y_2"), times = 4), # pivot_wider用の列
    sign = rep(c(-1, 1), each = 4), # 始点・終点の計算用の符号
    lambda = solve(nu * matrix(c(unique(w_11), w_12, w_12, w_22), nrow = D, ncol = D)) |> # 共分散行列の期待値を計算
      (\(.){eigen(.)[["values"]]})() |> 
      (\(.){rep(sqrt(.), times = 4)})(), # 固有値
    u = solve(nu * matrix(c(unique(w_11), w_12, w_12, w_22), nrow = D, ncol = D)) |> # 共分散行列の期待値を計算
      (\(.){t(eigen(.)[["vectors"]])})() |> 
      (\(.){rep(as.vector(.), times = 2)})(), # 固有ベクトル
    value = sign * u * lambda, # 軸の始点・終点を計算
    parameter = paste0("nu=", nu, ", W=(", unique(w_11), ", ", w_12, ", ", w_12, ", ", w_22, ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11_vals, ", ", w_12, ", ", w_12, ", ", w_22, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() |> # グループ化を解除
  dplyr::select(axis, name, value, parameter) |> # 利用する列を選択
  tidyr::pivot_wider(
    id_cols = c(axis, parameter), 
    names_from = name, 
    values_from = value
  ) # 軸の視点・終点の列を分割
anime_E_axis_df
## # A tibble: 182 × 6
##    axis  parameter                   xstart  ystart    xend    yend
##    <chr> <fct>                        <dbl>   <dbl>   <dbl>   <dbl>
##  1 y_1   nu=20, W=(1, 0.6, 0.6, 5)   0.232  -0.0340 -0.232   0.0340
##  2 y_2   nu=20, W=(1, 0.6, 0.6, 5)   0.0144  0.0981 -0.0144 -0.0981
##  3 y_1   nu=20, W=(1.1, 0.6, 0.6, 5) 0.220  -0.0331 -0.220   0.0331
##  4 y_2   nu=20, W=(1.1, 0.6, 0.6, 5) 0.0147  0.0980 -0.0147 -0.0980
##  5 y_1   nu=20, W=(1.2, 0.6, 0.6, 5) 0.210  -0.0324 -0.210   0.0324
##  6 y_2   nu=20, W=(1.2, 0.6, 0.6, 5) 0.0151  0.0979 -0.0151 -0.0979
##  7 y_1   nu=20, W=(1.3, 0.6, 0.6, 5) 0.201  -0.0318 -0.201   0.0318
##  8 y_2   nu=20, W=(1.3, 0.6, 0.6, 5) 0.0155  0.0978 -0.0155 -0.0978
##  9 y_1   nu=20, W=(1.4, 0.6, 0.6, 5) 0.193  -0.0314 -0.193   0.0314
## 10 y_2   nu=20, W=(1.4, 0.6, 0.6, 5) 0.0159  0.0978 -0.0159 -0.0978
## # … with 172 more rows


 パラメータごとにサンプリングして距離を計算します。

# 楕円の数(サンプルサイズ)を指定
N <- 30

# パラメータごとに楕円のサンプルを計算
anime_sample_delta_df <- tidyr::expand_grid(
  w_11 = w_11_vals, # パラメータ
  n = 1:N, # サンプル番号
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # パラメータとサンプルごとに格子点を複製
  dplyr::group_by(w_11, n) |> # 距離の計算用にグループ化
  dplyr::mutate(
    delta = mahalanobis(
      x = x_mat, 
      center = 0, 
      cov = rWishart(
        n = 1, df = nu, Sigma = matrix(c(unique(w_11), w_12, w_12, w_22), nrow = D, ncol = D)
      )[, , 1], # 精度行列を生成
      inverted = TRUE
    ), 
    parameter = paste0("nu=", nu, ", W=(", unique(w_11), ", ", w_12, ", ", w_12, ", ", w_22, ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11_vals, ", ", w_12, ", ", w_12, ", ", w_22, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_sample_delta_df
## # A tibble: 27,300,000 × 6
##     w_11     n    x_1    x_2 delta parameter                
##    <dbl> <int>  <dbl>  <dbl> <dbl> <fct>                    
##  1     1     1 -0.348 -0.311  24.3 nu=20, W=(1, 0.6, 0.6, 5)
##  2     1     1 -0.348 -0.305  23.6 nu=20, W=(1, 0.6, 0.6, 5)
##  3     1     1 -0.348 -0.299  22.8 nu=20, W=(1, 0.6, 0.6, 5)
##  4     1     1 -0.348 -0.293  22.1 nu=20, W=(1, 0.6, 0.6, 5)
##  5     1     1 -0.348 -0.286  21.4 nu=20, W=(1, 0.6, 0.6, 5)
##  6     1     1 -0.348 -0.280  20.7 nu=20, W=(1, 0.6, 0.6, 5)
##  7     1     1 -0.348 -0.274  20.1 nu=20, W=(1, 0.6, 0.6, 5)
##  8     1     1 -0.348 -0.267  19.4 nu=20, W=(1, 0.6, 0.6, 5)
##  9     1     1 -0.348 -0.261  18.7 nu=20, W=(1, 0.6, 0.6, 5)
## 10     1     1 -0.348 -0.255  18.1 nu=20, W=(1, 0.6, 0.6, 5)
## # … with 27,299,990 more rows


 パラメータごとにサンプリングして軸を計算します。

# 軸の数(サンプルサイズ)を指定
N <- 100

# パラメータごとに楕円の長軸のサンプルを計算
anime_sample_axis_df <- tidyr::expand_grid(
  w_11 = w_11_vals, # パラメータ
  n = 1:N # サンプル番号
) |> # パラメータごとに受け皿を複製
  dplyr::group_by(w_11, n) |> # 軸の計算用にグループ化
  dplyr::mutate(
    eigen_lt = rWishart(
      n = 1, 
      df = nu, 
      Sigma = matrix(c(unique(w_11), w_12, w_12, w_22), nrow = D, ncol = D)
    )[, , 1] |> # 精度行列を生成
      solve() |> # 分散共分散行列に変換
      eigen() |> # 固有値・固有ベクトルを計算
      list(), # リストに格納
    lambda_1 = eigen_lt[[1]][["values"]][1], # 固有値
    u_11 = eigen_lt[[1]][["vectors"]][1, 1], # 固有ベクトルの成分
    u_12 = eigen_lt[[1]][["vectors"]][2, 1], # 固有ベクトルの成分
    xend = u_11 * sqrt(lambda_1), # 終点のx軸の値
    yend = u_12 * sqrt(lambda_1), # 終点のy軸の値
    parameter = paste0("nu=", nu, ", W=(", unique(w_11), ", ", w_12, ", ", w_12, ", ", w_22, ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11_vals, ", ", w_12, ", ", w_12, ", ", w_22, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_sample_axis_df
## # A tibble: 9,100 × 9
##     w_11     n eigen_lt lambda_1   u_11   u_12   xend    yend parameter         
##    <dbl> <int> <list>      <dbl>  <dbl>  <dbl>  <dbl>   <dbl> <fct>             
##  1     1     1 <eigen>    0.0470 -0.994 0.112  -0.215 0.0242  nu=20, W=(1, 0.6,…
##  2     1     2 <eigen>    0.0556 -0.977 0.213  -0.230 0.0502  nu=20, W=(1, 0.6,…
##  3     1     3 <eigen>    0.0510 -0.984 0.178  -0.222 0.0402  nu=20, W=(1, 0.6,…
##  4     1     4 <eigen>    0.0839 -0.991 0.136  -0.287 0.0394  nu=20, W=(1, 0.6,…
##  5     1     5 <eigen>    0.0460 -0.999 0.0435 -0.214 0.00933 nu=20, W=(1, 0.6,…
##  6     1     6 <eigen>    0.0835 -0.972 0.235  -0.281 0.0678  nu=20, W=(1, 0.6,…
##  7     1     7 <eigen>    0.124  -0.992 0.129  -0.350 0.0455  nu=20, W=(1, 0.6,…
##  8     1     8 <eigen>    0.0525 -0.997 0.0759 -0.228 0.0174  nu=20, W=(1, 0.6,…
##  9     1     9 <eigen>    0.0658 -0.974 0.227  -0.250 0.0582  nu=20, W=(1, 0.6,…
## 10     1    10 <eigen>    0.0492 -0.997 0.0827 -0.221 0.0183  nu=20, W=(1, 0.6,…
## # … with 9,090 more rows


 パラメータごとに期待値のラベルを作成します。

# パラメータごとに共分散行列の期待値のラベルを作成
anime_label_df <- tibble::tibble(
  w_11 = w_11_vals # パラメータ
) |> 
  dplyr::group_by(w_11) |> # 期待値の計算用にグループ化
  dplyr::mutate(
    label = paste0("E[Lambda]=(", paste0(round(unique(nu) * matrix(c(unique(w_11), w_12, w_12, w_22), nrow = D, ncol = D), 2), collapse = ", "), ")"), # 分散共分散行列ラベル
    parameter = paste0("nu=", nu, ", W=(", unique(w_11), ", ", w_12, ", ", w_12, ", ", w_22, ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11_vals, ", ", w_12, ", ", w_12, ", ", w_22, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_label_df
## # A tibble: 91 × 3
##     w_11 label                       parameter                  
##    <dbl> <chr>                       <fct>                      
##  1   1   E[Lambda]=(20, 12, 12, 100) nu=20, W=(1, 0.6, 0.6, 5)  
##  2   1.1 E[Lambda]=(22, 12, 12, 100) nu=20, W=(1.1, 0.6, 0.6, 5)
##  3   1.2 E[Lambda]=(24, 12, 12, 100) nu=20, W=(1.2, 0.6, 0.6, 5)
##  4   1.3 E[Lambda]=(26, 12, 12, 100) nu=20, W=(1.3, 0.6, 0.6, 5)
##  5   1.4 E[Lambda]=(28, 12, 12, 100) nu=20, W=(1.4, 0.6, 0.6, 5)
##  6   1.5 E[Lambda]=(30, 12, 12, 100) nu=20, W=(1.5, 0.6, 0.6, 5)
##  7   1.6 E[Lambda]=(32, 12, 12, 100) nu=20, W=(1.6, 0.6, 0.6, 5)
##  8   1.7 E[Lambda]=(34, 12, 12, 100) nu=20, W=(1.7, 0.6, 0.6, 5)
##  9   1.8 E[Lambda]=(36, 12, 12, 100) nu=20, W=(1.8, 0.6, 0.6, 5)
## 10   1.9 E[Lambda]=(38, 12, 12, 100) nu=20, W=(1.9, 0.6, 0.6, 5)
## # … with 81 more rows


 「自由度の影響」のコードで作図できます。


ウィシャート分布の逆スケール行列とマハラノビス距離とその軸の関係

 $w_{1,1} < w_{2,2}$なほど楕円の軸の傾きが小さく、$w_{1,1} < w_{2,2}$なほど楕円の軸の傾きが大きくなっています。また、$w_{1,1}, w_{2,2}$の差が小さいと、サンプルによる楕円・軸と期待値による楕円・軸のばらつきが大きくなっています。

スケール行列(1,2成分)の影響

 $w_{1,2} = w_{2,1}$の値を変化させ、$\nu, w_{1,1}, w_{2,2}$を固定します。

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

# 次元数を指定
D <- 2

# 自由度を指定
nu <- D + 18

# 逆スケール行列を指定
w_11 <- 6
w_12_vals <- seq(from = -4, to = 4, by = 0.1) |> 
  round(digits = 2)
w_22 <- 5

# フレーム数を設定
frame_num <- length(w_12_vals)
frame_num
## [1] 81

 値の間隔が一定になるように$w_{1,2} = w_{2,1}$の値をw_12_valsとして作成します。また、$w_{1,1}$をw_11、$w_{2,2}$の値をw_22として値を指定します。

 設定したパラメータに応じて、変数の値$\mathbf{x}$を作成します。

# 作図用の共分散行列を作成
E_sigma_dd <- solve(nu * matrix(c(w_11, median(w_12_vals), median(w_12_vals), w_22), nrow = D, ncol = D))

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

# xの点を作成
x_mat <- tidyr::expand_grid(
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # 格子点を作成
  as.matrix() # マトリクスに変換
head(x_mat)
##             x_1        x_2
## [1,] -0.2738613 -0.3000000
## [2,] -0.2738613 -0.2939394
## [3,] -0.2738613 -0.2878788
## [4,] -0.2738613 -0.2818182
## [5,] -0.2738613 -0.2757576
## [6,] -0.2738613 -0.2696970


 パラメータごとに期待値による距離を計算します。

# パラメータごとに共分散行列の期待値によるマハラノビス距離(楕円)を計算
anime_E_delta_df <- tidyr::expand_grid(
  w_12 = w_12_vals, # パラメータ
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # パラメータごとに格子点を複製
  dplyr::group_by(w_12) |> # 距離の計算用にグループ化
  dplyr::mutate(
    delta = mahalanobis(
      x = x_mat, 
      center = 0, 
      cov = nu * matrix(c(w_11, unique(w_12), unique(w_12), w_22), nrow = D, ncol = D), # 精度行列の期待値を計算
      inverted = TRUE
    ), # 確率密度
    parameter = paste0("nu=", nu, ", W=(", w_11, ", ", unique(w_12), ", ", unique(w_12), ", ", w_22, ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11, ", ", w_12_vals, ", ", w_12_vals, ", ", w_22, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_E_delta_df
## # A tibble: 810,000 × 5
##     w_12    x_1    x_2 delta parameter              
##    <dbl>  <dbl>  <dbl> <dbl> <fct>                  
##  1    -4 -0.274 -0.3    4.85 nu=20, W=(6, -4, -4, 5)
##  2    -4 -0.274 -0.294  4.76 nu=20, W=(6, -4, -4, 5)
##  3    -4 -0.274 -0.288  4.67 nu=20, W=(6, -4, -4, 5)
##  4    -4 -0.274 -0.282  4.59 nu=20, W=(6, -4, -4, 5)
##  5    -4 -0.274 -0.276  4.52 nu=20, W=(6, -4, -4, 5)
##  6    -4 -0.274 -0.270  4.46 nu=20, W=(6, -4, -4, 5)
##  7    -4 -0.274 -0.264  4.40 nu=20, W=(6, -4, -4, 5)
##  8    -4 -0.274 -0.258  4.35 nu=20, W=(6, -4, -4, 5)
##  9    -4 -0.274 -0.252  4.31 nu=20, W=(6, -4, -4, 5)
## 10    -4 -0.274 -0.245  4.27 nu=20, W=(6, -4, -4, 5)
## # … with 809,990 more rows

 「自由度の影響」のnu, nu_valsw_12, w_12_valsを置き換えるように処理します。

 パラメータごとに期待値による軸を計算します。

# パラメータごとに共分散行列の期待値による楕円の軸を計算
anime_E_axis_df <- tidyr::expand_grid(
  w_12 = w_12_vals, # パラメータ
  name = paste0(rep(c("x", "y"), each = 2, times = 2), rep(c("start", "end"), each = 4)) # 列名
) |> # パラメータごとに受け皿を複製
  dplyr::group_by(w_12) |> # 軸の計算用にグループ化
  dplyr::mutate(
    axis = rep(c("y_1", "y_2"), times = 4), # pivot_wider用の列
    sign = rep(c(-1, 1), each = 4), # 始点・終点の計算用の符号
    lambda = solve(nu * matrix(c(w_11, unique(w_12), unique(w_12), w_22), nrow = D, ncol = D)) |> # 共分散行列の期待値を計算
      (\(.){eigen(.)[["values"]]})() |> 
      (\(.){rep(sqrt(.), times = 4)})(), # 固有値
    u = solve(nu * matrix(c(w_11, unique(w_12), unique(w_12), w_22), nrow = D, ncol = D)) |> # 共分散行列の期待値を計算
      (\(.){t(eigen(.)[["vectors"]])})() |> 
      (\(.){rep(as.vector(.), times = 2)})(), # 固有ベクトル
    value = sign * u * lambda, # 軸の始点・終点を計算
    parameter = paste0("nu=", nu, ", W=(", w_11, ", ", unique(w_12), ", ", unique(w_12), ", ", w_22, ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11, ", ", w_12_vals, ", ", w_12_vals, ", ", w_22, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() |> # グループ化を解除
  dplyr::select(axis, name, value, parameter) |> # 利用する列を選択
  tidyr::pivot_wider(
    id_cols = c(axis, parameter), 
    names_from = name, 
    values_from = value
  ) # 軸の視点・終点の列を分割
anime_E_axis_df
## # A tibble: 162 × 6
##    axis  parameter                    xstart  ystart    xend   yend
##    <chr> <fct>                         <dbl>   <dbl>   <dbl>  <dbl>
##  1 y_1   nu=20, W=(6, -4, -4, 5)     -0.122  -0.138   0.122  0.138 
##  2 y_2   nu=20, W=(6, -4, -4, 5)      0.0543 -0.0479 -0.0543 0.0479
##  3 y_1   nu=20, W=(6, -3.9, -3.9, 5) -0.118  -0.134   0.118  0.134 
##  4 y_2   nu=20, W=(6, -3.9, -3.9, 5)  0.0547 -0.0481 -0.0547 0.0481
##  5 y_1   nu=20, W=(6, -3.8, -3.8, 5) -0.114  -0.130   0.114  0.130 
##  6 y_2   nu=20, W=(6, -3.8, -3.8, 5)  0.0550 -0.0483 -0.0550 0.0483
##  7 y_1   nu=20, W=(6, -3.7, -3.7, 5) -0.111  -0.127   0.111  0.127 
##  8 y_2   nu=20, W=(6, -3.7, -3.7, 5)  0.0554 -0.0484 -0.0554 0.0484
##  9 y_1   nu=20, W=(6, -3.6, -3.6, 5) -0.108  -0.123   0.108  0.123 
## 10 y_2   nu=20, W=(6, -3.6, -3.6, 5)  0.0558 -0.0486 -0.0558 0.0486
## # … with 152 more rows


 パラメータごとにサンプリングして距離を計算します。

# 楕円の数(サンプルサイズ)を指定
N <- 30

# パラメータごとに楕円のサンプルを計算
anime_sample_delta_df <- tidyr::expand_grid(
  w_12 = w_12_vals, # パラメータ
  n = 1:N, # サンプル番号
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # パラメータとサンプルごとに格子点を複製
  dplyr::group_by(w_12, n) |> # 距離の計算用にグループ化
  dplyr::mutate(
    delta = mahalanobis(
      x = x_mat, 
      center = 0, 
      cov = rWishart(
        n = 1, df = nu, Sigma = matrix(c(w_11, unique(w_12), unique(w_12), w_22), nrow = D, ncol = D)
      )[, , 1], # 精度行列を生成
      inverted = TRUE
    ), 
    parameter = paste0("nu=", nu, ", W=(", w_11, ", ", unique(w_12), ", ", unique(w_12), ", ", w_22, ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11, ", ", w_12_vals, ", ", w_12_vals, ", ", w_22, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_sample_delta_df
## # A tibble: 24,300,000 × 6
##     w_12     n    x_1    x_2 delta parameter              
##    <dbl> <int>  <dbl>  <dbl> <dbl> <fct>                  
##  1    -4     1 -0.274 -0.3    1.88 nu=20, W=(6, -4, -4, 5)
##  2    -4     1 -0.274 -0.294  1.84 nu=20, W=(6, -4, -4, 5)
##  3    -4     1 -0.274 -0.288  1.80 nu=20, W=(6, -4, -4, 5)
##  4    -4     1 -0.274 -0.282  1.77 nu=20, W=(6, -4, -4, 5)
##  5    -4     1 -0.274 -0.276  1.75 nu=20, W=(6, -4, -4, 5)
##  6    -4     1 -0.274 -0.270  1.73 nu=20, W=(6, -4, -4, 5)
##  7    -4     1 -0.274 -0.264  1.71 nu=20, W=(6, -4, -4, 5)
##  8    -4     1 -0.274 -0.258  1.70 nu=20, W=(6, -4, -4, 5)
##  9    -4     1 -0.274 -0.252  1.69 nu=20, W=(6, -4, -4, 5)
## 10    -4     1 -0.274 -0.245  1.69 nu=20, W=(6, -4, -4, 5)
## # … with 24,299,990 more rows


 パラメータごとにサンプリングして軸を計算します。

# 軸の数(サンプルサイズ)を指定
N <- 100

# パラメータごとに楕円の長軸のサンプルを計算
anime_sample_axis_df <- tidyr::expand_grid(
  w_12 = w_12_vals, # パラメータ
  n = 1:N # サンプル番号
) |> 
  dplyr::group_by(w_12, n) |> # 軸の計算用にグループ化
  dplyr::mutate(
    eigen_lt = rWishart(
      n = 1, 
      df = nu, 
      Sigma = matrix(c(w_11, unique(w_12), unique(w_12), w_22), nrow = D, ncol = D)
    )[, , 1] |> # 精度行列を生成
      solve() |> # 分散共分散行列に変換
      eigen() |> # 固有値・固有ベクトルを計算
      list(), # リストに格納
    lambda_1 = eigen_lt[[1]][["values"]][1], # 固有値
    u_11 = eigen_lt[[1]][["vectors"]][1, 1], # 固有ベクトルの成分
    u_12 = eigen_lt[[1]][["vectors"]][2, 1], # 固有ベクトルの成分
    xend = u_11 * sqrt(lambda_1), # 終点のx軸の値
    yend = u_12 * sqrt(lambda_1), # 終点のy軸の値
    parameter = paste0("nu=", nu, ", W=(", w_11, ", ", unique(w_12), ", ", unique(w_12), ", ", w_22, ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11, ", ", w_12_vals, ", ", w_12_vals, ", ", w_22, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_sample_axis_df
## # A tibble: 8,100 × 9
##     w_12     n eigen_lt lambda_1   u_11   u_12   xend    yend parameter         
##    <dbl> <int> <list>      <dbl>  <dbl>  <dbl>  <dbl>   <dbl> <fct>             
##  1    -4     1 <eigen>    0.0408 -0.825 -0.565 -0.167 -0.114  nu=20, W=(6, -4, …
##  2    -4     2 <eigen>    0.0418  0.652  0.758  0.133  0.155  nu=20, W=(6, -4, …
##  3    -4     3 <eigen>    0.0713 -0.835 -0.551 -0.223 -0.147  nu=20, W=(6, -4, …
##  4    -4     4 <eigen>    0.0999 -0.725 -0.689 -0.229 -0.218  nu=20, W=(6, -4, …
##  5    -4     5 <eigen>    0.0342 -0.711 -0.704 -0.131 -0.130  nu=20, W=(6, -4, …
##  6    -4     6 <eigen>    0.0394  0.666  0.746  0.132  0.148  nu=20, W=(6, -4, …
##  7    -4     7 <eigen>    0.0214 -0.725 -0.689 -0.106 -0.101  nu=20, W=(6, -4, …
##  8    -4     8 <eigen>    0.0251 -0.790 -0.613 -0.125 -0.0971 nu=20, W=(6, -4, …
##  9    -4     9 <eigen>    0.0416  0.612  0.791  0.125  0.161  nu=20, W=(6, -4, …
## 10    -4    10 <eigen>    0.0598  0.445  0.896  0.109  0.219  nu=20, W=(6, -4, …
## # … with 8,090 more rows


 パラメータごとに期待値のラベルを作成します。

# パラメータごとに精度行列の期待値のラベルを作成
anime_label_df <- tibble::tibble(
  w_12 = w_12_vals # パラメータ
) |> 
  dplyr::group_by(w_12) |> # 期待値の計算用にグループ化
  dplyr::mutate(
    label = paste0("E[Lambda]=(", paste0(round(nu * matrix(c(w_11, unique(w_12), unique(w_12), w_22), nrow = D, ncol = D), 2), collapse = ", "), ")"), # 分散共分散行列ラベル
    parameter = paste0("nu=", nu, ", W=(", w_11, ", ", unique(w_12), ", ", unique(w_12), ", ", w_22, ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11, ", ", w_12_vals, ", ", w_12_vals, ", ", w_22, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_label_df
## # A tibble: 81 × 3
##     w_12 label                          parameter                  
##    <dbl> <chr>                          <fct>                      
##  1  -4   E[Lambda]=(120, -80, -80, 100) nu=20, W=(6, -4, -4, 5)    
##  2  -3.9 E[Lambda]=(120, -78, -78, 100) nu=20, W=(6, -3.9, -3.9, 5)
##  3  -3.8 E[Lambda]=(120, -76, -76, 100) nu=20, W=(6, -3.8, -3.8, 5)
##  4  -3.7 E[Lambda]=(120, -74, -74, 100) nu=20, W=(6, -3.7, -3.7, 5)
##  5  -3.6 E[Lambda]=(120, -72, -72, 100) nu=20, W=(6, -3.6, -3.6, 5)
##  6  -3.5 E[Lambda]=(120, -70, -70, 100) nu=20, W=(6, -3.5, -3.5, 5)
##  7  -3.4 E[Lambda]=(120, -68, -68, 100) nu=20, W=(6, -3.4, -3.4, 5)
##  8  -3.3 E[Lambda]=(120, -66, -66, 100) nu=20, W=(6, -3.3, -3.3, 5)
##  9  -3.2 E[Lambda]=(120, -64, -64, 100) nu=20, W=(6, -3.2, -3.2, 5)
## 10  -3.1 E[Lambda]=(120, -62, -62, 100) nu=20, W=(6, -3.1, -3.1, 5)
## # … with 71 more rows


 「自由度の影響」のコードで作図できます。


ウィシャート分布の逆スケール行列とマハラノビス距離とその軸の関係

 $w_{1,2} < 0$なほど軸の傾きが大きく、$w_{1,2} > 0$なほど軸の傾きが小さくなっています。

スケール行列(2,2成分)の影響

 $w_{2,2}$の値を変化させ、$\nu, w_{1,1}, w_{1,2}$を固定します。

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

# 次元数を指定
D <- 2

# 自由度を指定
nu <- D + 18

# 逆スケール行列を指定
w_11 <- 5
w_12 <- 0.6
w_22_vals <- seq(from = 1, to = 10, by = 0.1) |> 
  round(digits = 2)

# フレーム数を設定
frame_num <- length(w_22_vals)
frame_num
## [1] 91

 値の間隔が一定になるように$w_{2,2}$の値をw_22_valsとして作成します。また、$w_{1,1}, w_{2,2}$をw_11, w_22として値を指定します。

 設定したパラメータに応じて、変数の値$\mathbf{x}$を作成します。

# 作図用の共分散行列を作成
E_sigma_dd <- solve(nu * matrix(c(w_11, w_12, w_12, min(w_22_vals)), nrow = D, ncol = D))

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

# xの点を作成
x_mat <- tidyr::expand_grid(
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # 格子点を作成
  as.matrix() # マトリクスに変換
head(x_mat)
##             x_1        x_2
## [1,] -0.3114205 -0.3481788
## [2,] -0.3114205 -0.3411448
## [3,] -0.3114205 -0.3341109
## [4,] -0.3114205 -0.3270770
## [5,] -0.3114205 -0.3200431
## [6,] -0.3114205 -0.3130092

 「自由度の影響」のnu, nu_valsw_22, w_22_valsを置き換えるように処理します。

 パラメータごとに期待値による距離を計算します。

# パラメータごとに共分散行列の期待値によるマハラノビス距離(楕円)を計算
anime_E_delta_df <- tidyr::expand_grid(
  w_22 = w_22_vals, # パラメータ
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # パラメータごとに格子点を複製
  dplyr::group_by(w_22) |> # 距離の計算用にグループ化
  dplyr::mutate(
    delta = mahalanobis(
      x = x_mat, 
      center = 0, 
      cov = nu * matrix(c(w_11, w_12, w_12, unique(w_22)), nrow = D, ncol = D), # 精度行列の期待値を指定
      inverted = TRUE
    ), # 確率密度
    parameter = paste0("nu=", nu, ", W=(", w_11, ", ", w_12, ", ", w_12, ", ", unique(w_22), ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11, ", ", w_12, ", ", w_12, ", ", w_22_vals, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_E_delta_df
## # A tibble: 910,000 × 5
##     w_22    x_1    x_2 delta parameter                
##    <dbl>  <dbl>  <dbl> <dbl> <fct>                    
##  1     1 -0.311 -0.348  14.7 nu=20, W=(5, 0.6, 0.6, 1)
##  2     1 -0.311 -0.341  14.6 nu=20, W=(5, 0.6, 0.6, 1)
##  3     1 -0.311 -0.334  14.4 nu=20, W=(5, 0.6, 0.6, 1)
##  4     1 -0.311 -0.327  14.3 nu=20, W=(5, 0.6, 0.6, 1)
##  5     1 -0.311 -0.320  14.1 nu=20, W=(5, 0.6, 0.6, 1)
##  6     1 -0.311 -0.313  14.0 nu=20, W=(5, 0.6, 0.6, 1)
##  7     1 -0.311 -0.306  13.9 nu=20, W=(5, 0.6, 0.6, 1)
##  8     1 -0.311 -0.299  13.7 nu=20, W=(5, 0.6, 0.6, 1)
##  9     1 -0.311 -0.292  13.6 nu=20, W=(5, 0.6, 0.6, 1)
## 10     1 -0.311 -0.285  13.5 nu=20, W=(5, 0.6, 0.6, 1)
## # … with 909,990 more rows


 パラメータごとに期待値による軸を計算します。

# パラメータごとに共分散行列の期待値による楕円の軸を計算
anime_E_axis_df <- tidyr::expand_grid(
  w_22 = w_22_vals, # パラメータごとに
  name = paste0(rep(c("x", "y"), each = 2, times = 2), rep(c("start", "end"), each = 4)) # 列名
) |> # パラメータごとに受け皿を複製
  dplyr::group_by(w_22) |> # 軸の計算用にグループ化
  dplyr::mutate(
    axis = rep(c("y_1", "y_2"), times = 4), # pivot_wider用の列
    sign = rep(c(-1, 1), each = 4), # 始点・終点の計算用の符号
    lambda = solve(nu * matrix(c(w_11, w_12, w_12, unique(w_22)), nrow = D, ncol = D)) |> # 共分散行列の期待値を計算
      (\(.){eigen(.)[["values"]]})() |> 
      (\(.){rep(sqrt(.), times = 4)})(), # 固有値
    u = solve(nu * matrix(c(w_11, w_12, w_12, unique(w_22)), nrow = D, ncol = D)) |> # 共分散行列の期待値を計算
      (\(.){t(eigen(.)[["vectors"]])})() |> 
      (\(.){rep(as.vector(.), times = 2)})(), # 固有ベクトル
    value = sign * u * lambda, # 軸の始点・終点を計算
    parameter = paste0("nu=", nu, ", W=(", w_11, ", ", w_12, ", ", w_12, ", ", unique(w_22), ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11, ", ", w_12, ", ", w_12, ", ", w_22_vals, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() |> # グループ化を解除
  dplyr::select(axis, name, value, parameter) |> # 利用する列を選択
  tidyr::pivot_wider(
    id_cols = c(axis, parameter), 
    names_from = name, 
    values_from = value
  ) # 軸の視点・終点の列を分割
anime_E_axis_df
## # A tibble: 182 × 6
##    axis  parameter                   xstart  ystart    xend    yend
##    <chr> <fct>                        <dbl>   <dbl>   <dbl>   <dbl>
##  1 y_1   nu=20, W=(5, 0.6, 0.6, 1)   0.0340 -0.232  -0.0340  0.232 
##  2 y_2   nu=20, W=(5, 0.6, 0.6, 1)   0.0981  0.0144 -0.0981 -0.0144
##  3 y_1   nu=20, W=(5, 0.6, 0.6, 1.1) 0.0331 -0.220  -0.0331  0.220 
##  4 y_2   nu=20, W=(5, 0.6, 0.6, 1.1) 0.0980  0.0147 -0.0980 -0.0147
##  5 y_1   nu=20, W=(5, 0.6, 0.6, 1.2) 0.0324 -0.210  -0.0324  0.210 
##  6 y_2   nu=20, W=(5, 0.6, 0.6, 1.2) 0.0979  0.0151 -0.0979 -0.0151
##  7 y_1   nu=20, W=(5, 0.6, 0.6, 1.3) 0.0318 -0.201  -0.0318  0.201 
##  8 y_2   nu=20, W=(5, 0.6, 0.6, 1.3) 0.0978  0.0155 -0.0978 -0.0155
##  9 y_1   nu=20, W=(5, 0.6, 0.6, 1.4) 0.0314 -0.193  -0.0314  0.193 
## 10 y_2   nu=20, W=(5, 0.6, 0.6, 1.4) 0.0978  0.0159 -0.0978 -0.0159
## # … with 172 more rows


 パラメータごとにサンプリングして距離を計算します。

# 楕円の数(サンプルサイズ)を指定
N <- 30

# パラメータごとに楕円のサンプルを計算
anime_sample_delta_df <- tidyr::expand_grid(
  w_22 = w_22_vals, # パラメータ
  n = 1:N, # サンプル番号
  x_1 = x_1_vals, 
  x_2 = x_2_vals
) |> # パラメータとサンプルごとに格子点を複製
  dplyr::group_by(w_22, n) |> # 距離の計算用にグループ化
  dplyr::mutate(
    delta = mahalanobis(
      x = x_mat, 
      center = 0, 
      cov = rWishart(
        n = 1, df = nu, Sigma = matrix(c(w_11, w_12, w_12, unique(w_22)), nrow = D, ncol = D)
      )[, , 1], # 精度行列を生成
      inverted = TRUE
    ), 
    parameter = paste0("nu=", nu, ", W=(", w_11, ", ", w_12, ", ", w_12, ", ", unique(w_22), ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11, ", ", w_12, ", ", w_12, ", ", w_22_vals, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_sample_delta_df
## # A tibble: 27,300,000 × 6
##     w_22     n    x_1    x_2 delta parameter                
##    <dbl> <int>  <dbl>  <dbl> <dbl> <fct>                    
##  1     1     1 -0.311 -0.348  14.0 nu=20, W=(5, 0.6, 0.6, 1)
##  2     1     1 -0.311 -0.341  13.9 nu=20, W=(5, 0.6, 0.6, 1)
##  3     1     1 -0.311 -0.334  13.7 nu=20, W=(5, 0.6, 0.6, 1)
##  4     1     1 -0.311 -0.327  13.6 nu=20, W=(5, 0.6, 0.6, 1)
##  5     1     1 -0.311 -0.320  13.5 nu=20, W=(5, 0.6, 0.6, 1)
##  6     1     1 -0.311 -0.313  13.3 nu=20, W=(5, 0.6, 0.6, 1)
##  7     1     1 -0.311 -0.306  13.2 nu=20, W=(5, 0.6, 0.6, 1)
##  8     1     1 -0.311 -0.299  13.1 nu=20, W=(5, 0.6, 0.6, 1)
##  9     1     1 -0.311 -0.292  12.9 nu=20, W=(5, 0.6, 0.6, 1)
## 10     1     1 -0.311 -0.285  12.8 nu=20, W=(5, 0.6, 0.6, 1)
## # … with 27,299,990 more rows


 パラメータごとにサンプリングして軸を計算します。

# 軸の数(サンプルサイズ)を指定
N <- 100

# パラメータごとに楕円の長軸のサンプルを計算
anime_sample_axis_df <- tidyr::expand_grid(
  w_22 = w_22_vals, # パラメータ
  n = 1:N # サンプル番号
) |> 
  dplyr::group_by(w_22, n) |> # 軸の計算用にグループ化
  dplyr::mutate(
    eigen_lt = rWishart(
      n = 1, 
      df = nu, 
      Sigma = matrix(c(w_11, w_12, w_12, unique(w_22)), nrow = D, ncol = D)
    )[, , 1] |> # 精度行列を生成
      solve() |> # 分散共分散行列に変換
      eigen() |> # 固有値・固有ベクトルを計算
      list(), # リストに格納
    lambda_1 = eigen_lt[[1]][["values"]][1], # 固有値
    u_11 = eigen_lt[[1]][["vectors"]][1, 1], # 固有ベクトルの成分
    u_12 = eigen_lt[[1]][["vectors"]][2, 1], # 固有ベクトルの成分
    xend = u_11 * sqrt(lambda_1), # 終点のx軸の値
    yend = u_12 * sqrt(lambda_1), # 終点のy軸の値
    parameter = paste0("nu=", nu, ", W=(", w_11, ", ", w_12, ", ", w_12, ", ", unique(w_22), ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11, ", ", w_12, ", ", w_12, ", ", w_22_vals, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_sample_axis_df
## # A tibble: 9,100 × 9
##     w_22     n eigen_lt lambda_1    u_11  u_12    xend  yend parameter          
##    <dbl> <int> <list>      <dbl>   <dbl> <dbl>   <dbl> <dbl> <fct>              
##  1     1     1 <eigen>    0.0420 -0.231  0.973 -0.0474 0.199 nu=20, W=(5, 0.6, …
##  2     1     2 <eigen>    0.0622 -0.243  0.970 -0.0607 0.242 nu=20, W=(5, 0.6, …
##  3     1     3 <eigen>    0.0727  0.0629 0.998  0.0170 0.269 nu=20, W=(5, 0.6, …
##  4     1     4 <eigen>    0.0719 -0.171  0.985 -0.0459 0.264 nu=20, W=(5, 0.6, …
##  5     1     5 <eigen>    0.0258 -0.334  0.942 -0.0537 0.151 nu=20, W=(5, 0.6, …
##  6     1     6 <eigen>    0.0803 -0.293  0.956 -0.0830 0.271 nu=20, W=(5, 0.6, …
##  7     1     7 <eigen>    0.0627 -0.110  0.994 -0.0275 0.249 nu=20, W=(5, 0.6, …
##  8     1     8 <eigen>    0.0688 -0.254  0.967 -0.0666 0.254 nu=20, W=(5, 0.6, …
##  9     1     9 <eigen>    0.0760 -0.187  0.982 -0.0515 0.271 nu=20, W=(5, 0.6, …
## 10     1    10 <eigen>    0.0520 -0.129  0.992 -0.0294 0.226 nu=20, W=(5, 0.6, …
## # … with 9,090 more rows


 パラメータごとに期待値のラベルを作成します。

# パラメータごとに分散共分散行列の期待値のラベルを作成
anime_label_df <- tibble::tibble(
  w_22 = w_22_vals # パラメータ
) |> 
  dplyr::group_by(w_22) |> # 期待値の計算用にグループ化
  dplyr::mutate(
    label = paste0("E[Lambda]=(", paste0(round(nu * matrix(c(w_11, w_12, w_12, unique(w_22)), nrow = D, ncol = D), 2), collapse = ", "), ")"), # 分散共分散行列ラベル
    parameter = paste0("nu=", nu, ", W=(", w_11, ", ", w_12, ", ", w_12, ", ", unique(w_22), ")") |> 
      factor(levels = paste0("nu=", nu, ", W=(", w_11, ", ", w_12, ", ", w_12, ", ", w_22_vals, ")")) # フレーム切替用ラベル
  ) |> 
  dplyr::ungroup() # グループ化を解除
anime_label_df
## # A tibble: 91 × 3
##     w_22 label                       parameter                  
##    <dbl> <chr>                       <fct>                      
##  1   1   E[Lambda]=(100, 12, 12, 20) nu=20, W=(5, 0.6, 0.6, 1)  
##  2   1.1 E[Lambda]=(100, 12, 12, 22) nu=20, W=(5, 0.6, 0.6, 1.1)
##  3   1.2 E[Lambda]=(100, 12, 12, 24) nu=20, W=(5, 0.6, 0.6, 1.2)
##  4   1.3 E[Lambda]=(100, 12, 12, 26) nu=20, W=(5, 0.6, 0.6, 1.3)
##  5   1.4 E[Lambda]=(100, 12, 12, 28) nu=20, W=(5, 0.6, 0.6, 1.4)
##  6   1.5 E[Lambda]=(100, 12, 12, 30) nu=20, W=(5, 0.6, 0.6, 1.5)
##  7   1.6 E[Lambda]=(100, 12, 12, 32) nu=20, W=(5, 0.6, 0.6, 1.6)
##  8   1.7 E[Lambda]=(100, 12, 12, 34) nu=20, W=(5, 0.6, 0.6, 1.7)
##  9   1.8 E[Lambda]=(100, 12, 12, 36) nu=20, W=(5, 0.6, 0.6, 1.8)
## 10   1.9 E[Lambda]=(100, 12, 12, 38) nu=20, W=(5, 0.6, 0.6, 1.9)
## # … with 81 more rows


 「自由度の影響」のコードで作図できます。


ウィシャート分布の逆スケール行列とマハラノビス距離とその軸の関係

 「1,1成分の影響」と同様に、$w_{1,1}, w_{2,2}$の大小関係によって軸の向きが変化しています。

 この記事では、ウィシャート分布をマハラノビス距離とその軸を用いて可視化しました。

参考文献

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

おわりに

 やってることは前回の記事と被ってるのですが、目的が違うので別記事にしました。
 色々面倒ながら頑張ったのですが、分かりやすいとは言えませんね。もう少し勉強して、他の記事との対応関係を整理しつつ構成し直したいです、来年くらいには。

 さて、この分布で終わりと思ってたのですが、まだガウス・ウィシャート分布がありました。忘れてました。というわけでまだ続きます。

 2022年9月11日は、Juice=Juiceのメジャーデビュー9周年記念日です!!!!!!

 気付けば私がハマってからの期間の方が長くなってました。これからも末永く聴いていたい。

【次の内容】

つづく