からっぽのしょこ

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

トップ > 確率分布 > 【R】ウィシャート分布の計算

【R】ウィシャート分布の計算

確率分布 確率分布-ウィシャート分布 R

はじめに

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

　この記事では、R言語でウィシャート分布に関する計算をします。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

ウィシャート分布の計算

　ウィシャート分布(Wishart Distribution)の確率密度と統計量を計算します。ウィシャート分布については「ウィシャート分布の定義式 - からっぽのしょこ」を参照してください。

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

# 利用パッケージ
library(MCMCpack)

　この記事では、パッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。
　magrittrパッケージのパイプ演算子%>%ではなく、ベースパイプ(ネイティブパイプ)演算子|>を使っています。%>%に置き換えても処理できます。

確率密度の計算

　ウィシャート分布に従う確率密度を計算する方法をいくつか確認します。

パラメータの設定

　ウィシャート分布の次元数$D$とパラメータ$\nu, \mathbf{W}$、確率変数の実現値$\boldsymbol{\Lambda}$を設定します。

# 次元数を指定
D <- 3

# 自由度を指定
nu <- D

# 逆スケール行列を指定
w_dd <- diag(D) # 単位行列
w_dd <- diag(c(10, 6)) # 対角行列
w_dd <- c(
  1.2, 0.6, -1.5, 
  0.6, 2.4, 0.4, 
  -1.5, 0.4, 3.6
) |> # 値を指定
  matrix(nrow = D, ncol = D) # マトリクスに変換

# 確率変数の値を指定
lambda_dd <- c(
  2, 0, -1, 
  0, 2.4, 0.4, 
  -1, 0.4, 4
) |> # 値を指定
  matrix(nrow = D, ncol = D) # マトリクスに変換

　自由度$\nu$と逆スケール行列$\mathbf{W}$を指定します。
　自由度は$\nu > D - 1$、逆スケール行列は$D \times D$の正定値行列を満たす必要があります。設定した値に従う確率密度を計算します。

スクラッチで計算

　定義式から確率密度を計算します。

# 定義式により確率密度を計算
C_w1 <- 1 / sqrt(det(w_dd))^nu / sqrt(2)^(nu*D) / sqrt(pi)^(0.5*D*(D-1))
C_w2 <- 1 / prod(gamma(0.5 * (nu + 1 - 1:D)))
dens <- C_w1 * C_w2 * sqrt(det(lambda_dd))^(nu-D-1) * exp(-0.5 * sum(diag(solve(w_dd) %*% lambda_dd)))
dens

## [1] 4.04105e-06

　ウィシャート分布は、次の式で定義されます。

$$ \begin{aligned} C_{\mathcal{W}} &= \left( |\mathbf{W}|^{\frac{\nu}{2}} 2^{\frac{\nu D}{2}} \pi^{\frac{D(D-1)}{4}} \prod_{d=1}^D \Gamma \Bigl( \frac{\nu + 1 - d}{2} \Bigr) \right)^{-1} \\ \mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W}) &= C_{\mathcal{W}} |\boldsymbol{\Lambda}|^{\frac{\nu-D-1}{2}} \exp \Bigl( - \frac{1}{2} \mathrm{Tr}(\mathbf{W}^{-1} \boldsymbol{\Lambda}) \Bigr) \end{aligned} $$

　ここで、$C_{\mathcal{W}}$はウィシャート分布の正規化係数、$\pi$は円周率、$|\mathbf{A}|$は行列式、$\mathrm{Tr}(\mathbf{A})$はトレース、$\mathbf{A}^{-1}$は逆行列です。
　円周率はpi、行列式はdet()、トレースはsum(diag())、逆行列solve()、行列の積は%*%演算子で計算できます。

　対数をとった定義式から計算します。

# 対数をとった定義式により確率密度を計算
log_C_w1 <- -0.5 * (nu*log(det(w_dd)) + nu*D*log(2) + 0.5*D*(D-1)*log(pi))
log_C_w2 <- -sum(lgamma(0.5 * (nu + 1 - 1:D)))
log_dens <- log_C_w1 + log_C_w2 + 0.5 * ((nu - D - 1) * log(det(lambda_dd)) - sum(diag(solve(w_dd) %*% lambda_dd)))
dens <- exp(log_dens)
dens; log_dens

## [1] 4.04105e-06
## [1] -12.41901

　対数をとった定義式を計算します。

$$ \begin{aligned} \log C_{\mathcal{W}} &= - \frac{\nu}{2} \log |\mathbf{W}| - \frac{\nu D}{2} \log 2 - \frac{D(D-1)}{4} \log \pi - \sum_{d=1}^D \log \Gamma \Bigl( \frac{\nu + 1 - d}{2} \Bigr) \\ \log \mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W}) &= \log C_{\mathcal{W}} + \frac{\nu - D - 1}{2} \log |\boldsymbol{\Lambda}| - \frac{1}{2} \mathrm{Tr}(\mathbf{W}^{-1} \boldsymbol{\Lambda}) \end{aligned} $$

　計算結果の指数をとると確率密度が得られます。

$$ \mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W}) = \exp \Bigr( \log \mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W}) \Bigr) $$

　指数と対数の性質より$\exp(\log x) = x$です。

　次は、関数を使って確率密度を計算します。

関数で計算

　MCMCpackパッケージのウィシャート分布の確率密度関数dwish()で計算します。

# ウィシャート分布の関数により確率密度を計算
dens <- MCMCpack::dwish(W = lambda_dd, v = nu, S = w_dd)
dens

## [1] 4.04105e-06

　確率変数の引数Wにlambda_dd、自由度の引数vにnu、逆スケール行列の引数Sにw_ddを指定します。

統計量の計算

　ウィシャート分布の統計量を計算します。

　期待値を計算します。

# 期待値を計算
E_lambda_dd <- nu * w_dd
E_lambda_dd

##      [,1] [,2] [,3]
## [1,]  3.6  1.8 -4.5
## [2,]  1.8  7.2  1.2
## [3,] -4.5  1.2 10.8

　ウィシャート分布の期待値は、次の式で計算できます。

$$ \mathbb{E}[\boldsymbol{\Lambda}] = \nu \mathbf{W} $$

　1つの成分の分散を計算します。

# 成分を指定
i <- 1
j <- 2

# 指定した成分の分散を計算
V_lambda_ij <- nu * (w_dd[i, j]^2 + w_dd[i, i] * w_dd[j, j])
V_lambda_ij

## [1] 9.72

　ウィシャート分布の各成分の分散は、次の式で計算できます。

$$ \mathbb{V}[\lambda_{i,j}] = \nu (w_{i,j}^2 + w_{i,i} w_{j,j}) $$

　全ての成分の分散を計算します。

# 対角成分を行方向に複製した行列を作成
w_ii <- w_dd |> 
  diag() |> # 対角成分を抽出
  rep(each = D) |> # 対角成分を複製
  matrix(nrow = D, ncol = D, byrow = TRUE) # マトリクスに変換

# 全ての成分の分散を計算
V_lambda_dd <- nu * (w_dd^2 + w_ii * t(w_ii))
V_lambda_dd; w_ii

##       [,1]  [,2]  [,3]
## [1,]  8.64  9.72 19.71
## [2,]  9.72 34.56 26.40
## [3,] 19.71 26.40 77.76
##      [,1] [,2] [,3]
## [1,]  1.2  1.2  1.2
## [2,]  2.4  2.4  2.4
## [3,]  3.6  3.6  3.6

　対角成分の積$w_{i,i} w_{j,j}$を計算するために、対角要素を各行(または列)の全ての要素に複製したマトリクスw_iiを作成して計算します。

　最頻値を計算します。

# 最頻値を計算:(ν > D+1)
mode_lambda_dd <- (nu - D - 1) * w_dd
mode_lambda_dd

##      [,1] [,2] [,3]
## [1,] -1.2 -0.6  1.5
## [2,] -0.6 -2.4 -0.4
## [3,]  1.5 -0.4 -3.6

　ウィシャート分布の最頻値は、次の式で計算できます。

$$ \mathrm{mode}[\boldsymbol{\Lambda}] = (\nu - D - 1) \mathbf{W} \quad (\nu > D + 1) $$


　この記事では、ウィシャート分布の計算を確認しました。次は、乱数を生成します。

参考文献

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

おわりに

　これまでは「計算」→「作図」→「乱数生成」の流れで進めてきましたが、ウィシャート分布のグラフ化は難しいので、「計算」→「乱数生成」→「可視化」の順番に進めます。「可視化」では 分散共分散行列(確率変数)を可視化することでウィシャート分布の理解を目指します。

【次の内容】

www.anarchive-beta.com