はじめに

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

　この記事では、R言語で多次元(多変量)スチューデントのt分布に関する計算をします。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

はじめに
多次元スチューデントのt分布の計算
- 確率密度の計算
- 統計量の計算
参考文献
おわりに

多次元スチューデントのt分布の計算

　多次元スチューデントのt分布(Multivariate Student's t Distribution)の確率密度を計算します。多次元t分布については「多次元スチューデントのt分布の定義式 - からっぽのしょこ」を参照してください。

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

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

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

確率密度の計算

　多次元スチューデントのt分布に従う確率密度を計算する方法をいくつか確認します。

スケール行列を使用

　まずは、スケール行列を用いてt分布を計算します。

パラメータの設定

　t分布の次元数$D$とパラメータ$\boldsymbol{\mu}, \boldsymbol{\Sigma}$、確率変数の値$\mathbf{x}$を設定します。

# 次元数を指定
D <- 3

# 自由度を指定
nu <- 5

# 位置ベクトルを指定
mu_d <- c(10, -6, 1.5)

# スケール行列を指定
sigma_dd <- c(
  4, 1.8, -0.1, 
  1.8, 9, 2.4, 
  -0.1, 2.4, 1
) |> # 値を指定
  matrix(nrow = D, ncol = D, byrow = TRUE) # マトリクスに変換

# 確率変数の値を指定
x_d <- c(11.5, -5, 0)

　自由度(形状パラメータ)$\nu$、位置ベクトル$\boldsymbol{\mu}$、スケール行列$\boldsymbol{\Sigma}$、確率変数の値$\mathbf{x}$を指定します。$\boldsymbol{\mu}, \mathbf{x}$は$D$次元ベクトル、$\boldsymbol{\Sigma}$は$D \times D$の行列です。

$$ \boldsymbol{\mu} = \begin{pmatrix} \mu_1 \\ \mu_2 \\ \vdots \\ \mu_D \end{pmatrix} ,\ \boldsymbol{\Sigma} = \begin{pmatrix} \sigma_{1,1} & \sigma_{1,2} & \cdots & \sigma_{1,D} \\ \sigma_{2,1} & \sigma_{2,2} & \cdots & \sigma_{2,D} \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{D,1} & \sigma_{D,2} & \cdots & \sigma_{D,D} \end{pmatrix} ,\ \mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_D \end{pmatrix} $$

　$x_d$は実数をとり、$\nu$は正の実数、$\mu_d$は実数、$\sigma_{d,d}$は正の実数、$\sigma_{i,j} = \sigma_{j,i}\ (i \neq j)$は実数、また$\boldsymbol{\Sigma}$は正定値行列を満たす必要があります。設定した値に従う確率密度を計算します。

スクラッチで計算

　定義式から計算します。

# 定義式により確率密度を計算
C_St1    <- gamma(0.5 * (nu + D)) / gamma(0.5 * nu)
C_St2    <- 1 / sqrt(pi * nu)^D / sqrt(det(sigma_dd))
tmp_term <- (t(x_d - mu_d) %*% solve(sigma_dd) %*% (x_d - mu_d)) |> 
  as.numeric()
term     <- 1 / sqrt(1 + tmp_term / nu)^(nu + D)
dens     <- C_St1 * C_St2 * term
dens

## [1] 0.0003309269

　t分布は、$\boldsymbol{\Sigma}$を用いて次の式で定義されます。

$$ \begin{aligned} C_{\mathrm{St}}(\nu, \boldsymbol{\Sigma}) &= \frac{ \Gamma(\frac{\nu + D}{2}) }{ \Gamma(\frac{\nu}{2}) } \frac{ 1 }{ (\pi \nu)^{\frac{D}{2}} |\boldsymbol{\Sigma}|^{\frac{1}{2}} } \\ \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Sigma}) &= C_{\mathrm{St}}(\nu, \boldsymbol{\Sigma}) \left\{ 1 + \frac{1}{\nu} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) \right\}^{-\frac{\nu+D}{2}} \end{aligned} $$

　ここで、$C_{\mathrm{St}}$はt分布の正規化係数、$\pi$は円周率、$\mathbf{A}^{\top}$は転置行列、$\mathbf{A}^{-1}$は逆行列、$|\mathbf{A}|$は行列式、2分の1乗は平方根$\sqrt{a} = a^{\frac{1}{2}}$です。
　円周率はpi、転置はt()、逆行列solve()、行列式はdet()、行列の積は%*%演算子で計算できます。

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

# 対数をとった定義式により確率密度を計算
log_C_St1 <- lgamma(0.5 * (nu + D)) - lgamma(0.5 * nu)
log_C_St2 <- -D * 0.5 * log(pi * nu) - 0.5 * log(det(sigma_dd))
tmp_term  <- (t(x_d - mu_d) %*% solve(sigma_dd) %*% (x_d - mu_d)) |> 
  as.numeric()
log_term  <- -(nu + D) * 0.5 * log(1 + tmp_term / nu)
log_dens  <- log_C_St1 + log_C_St2 + log_term
dens      <- exp(log_dens)
dens; log_dens

## [1] 0.0003309269
## [1] -8.013613

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

$$ \begin{aligned} \log C_{\mathrm{St}}(\nu, \boldsymbol{\Sigma}) &= \log \Gamma \Bigl(\frac{\nu + D}{2}\Bigr) - \log \Gamma \Bigl(\frac{\nu}{2}\Bigr) - \frac{D}{2} \log (\pi \nu) - \frac{1}{2} \log |\boldsymbol{\Sigma}| \\ \log \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Sigma}) &= \log C_{\mathrm{St}}(\nu, \boldsymbol{\Sigma}) - \frac{\nu+D}{2} \log \Bigl\{ 1 + \frac{1}{\nu} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) \Bigr\} \end{aligned} $$

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

$$ \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \exp \Bigr( \log \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Sigma}) \Bigr) $$

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

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

関数で計算

　mvnfastパッケージの多次元t分布の確率密度関数dmvt()で計算します。

# 関数により確率密度を計算
dens <- mvnfast::dmvt(X = x_d, mu = mu_d, sigma = sigma_dd, df = nu)
dens

## [1] 0.0003309269

　確率変数の引数Xにx_d、位置ベクトルの引数muにmu_d、スケール行列の引数sigmaにsigma_dd、自由度の引数dfにnuを指定します。

　log = TRUEを指定すると対数をとった確率密度を返します。

# 関数により確率密度を計算
log_dens <- mvnfast::dmvt(X = x_d, mu = mu_d, sigma = sigma_dd, df = nu, log = TRUE)
dens     <- exp(log_dens)
dens; log_dens

## [1] 0.0003309269
## [1] -8.013613

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

　LaplacesDemonパッケージのdmvt()でも計算できます。

# 関数により確率密度を計算
dens <- LaplacesDemon::dmvt(x = x_d, mu = mu_d, S = sigma_dd, df = nu)
dens

# 関数により確率密度を計算
log_dens <- LaplacesDemon::dmvt(x = x_d, mu = mu_d, S = sigma_dd, df = nu)
dens     <- exp(log_dens)
dens; log_dens

## [1] 0.0003309269
## [1] 1.000331
## [1] 0.0003309269

　確率変数の引数xにx_d、位置ベクトルの引数muにmu_d、スケール行列の引数Sにsigma_dd、自由度の引数dfにnuを指定します。

逆スケール行列を使用

　次は、逆スケール行列を用いてt分布を計算します。

パラメータの設定

　スケール行列$\boldsymbol{\Sigma}$の代わりに、逆スケール行列$\boldsymbol{\Lambda}$を指定します。

# 逆スケール行列を指定
lambda_dd <- c(
  4, 1.8, -0.1, 
  1.8, 9, 2.4, 
  -0.1, 2.4, 1
) |> # 値を指定
  matrix(nrow = D, ncol = D, byrow = TRUE) # マトリクスに変換

　あるいは、$\boldsymbol{\Sigma}$を指定して$\boldsymbol{\Lambda}$を計算します。

# 逆スケール行列を計算
lambda_dd <- solve(sigma_dd)
lambda_dd

##            [,1]       [,2]       [,3]
## [1,]  0.3696099 -0.2327173  0.5954825
## [2,] -0.2327173  0.4551677 -1.1156742
## [3,]  0.5954825 -1.1156742  3.7371663

　逆スケール行列は、スケール行列の逆行列$\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{-1}$です。

スクラッチで計算

　定義式から計算します。

# 定義式により確率密度を計算
C_St1    <- gamma(0.5 * (nu + D)) / gamma(0.5 * nu)
C_St2    <- sqrt(det(lambda_dd)) / sqrt(pi * nu)^D
tmp_term <- (t(x_d - mu_d) %*% lambda_dd %*% (x_d - mu_d)) |> 
  as.numeric()
term     <- 1 / sqrt(1 + tmp_term / nu)^(nu + D)
dens     <- C_St1 * C_St2 * term
dens

## [1] 0.0003309269

　t分布は、$\boldsymbol{\Lambda}$を用いて次の式で定義されます。

$$ \begin{aligned} C_{\mathrm{St}}(\nu, \boldsymbol{\Lambda}) &= \frac{ \Gamma(\frac{\nu + D}{2}) }{ \Gamma(\frac{\nu}{2}) } \frac{ |\boldsymbol{\Lambda}|^{\frac{1}{2}} }{ (\pi \nu)^{\frac{D}{2}} } \\ \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Lambda}) &= C_{\mathrm{St}}(\nu, \boldsymbol{\Lambda}) \left\{ 1 + \frac{1}{\nu} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Lambda} (\mathbf{x} - \boldsymbol{\mu}) \right\}^{-\frac{\nu+D}{2}} \end{aligned} $$

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

# 対数をとった定義式により確率密度を計算
log_C_St1 <- lgamma(0.5 * (nu + D)) - lgamma(0.5 * nu)
log_C_St2 <- 0.5 * log(det(lambda_dd)) - D * 0.5 * log(pi * nu)
tmp_term  <- (t(x_d - mu_d) %*% lambda_dd %*% (x_d - mu_d)) |> 
  as.numeric()
term      <- -(nu + D) * 0.5 * log(1 + tmp_term / nu)
log_dens  <- log_C_St1 + log_C_St2 + log_term
dens      <- exp(log_dens)
dens; log_dens

## [1] 0.0003309269
## [1] -8.013613

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

$$ \begin{aligned} \log C_{\mathrm{St}}(\nu, \boldsymbol{\Lambda}) &= \log \Gamma \Bigl(\frac{\nu + D}{2}\Bigr) - \log \Gamma \Bigl(\frac{\nu}{2}\Bigr) - \frac{D}{2} \log (\pi \nu) + \frac{1}{2} \log |\boldsymbol{\Lambda}| \\ \log \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Lambda}) &= \log C_{\mathrm{St}}(\nu, \boldsymbol{\Lambda}) - \frac{\nu+D}{2} \log \Bigl\{ 1 + \frac{1}{\nu} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Lambda} (\mathbf{x} - \boldsymbol{\mu}) \Bigr\} \end{aligned} $$

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

関数で計算

　mvnfastパッケージのdmvt()で計算します。

# 関数により確率密度を計算
dens <- mvnfast::dmvt(X = x_d, mu = mu_d, sigma = solve(lambda_dd), df = nu)
dens

## [1] 0.0003309269

　スケール行列の引数sigmaにlambda_ddの逆行列を指定します。

　LaplacesDemonパッケージのdmvt()で計算します。

# 関数により確率密度を計算
dens <- LaplacesDemon::dmvt(x = x_d, mu = mu_d, S = sigma_dd, df = nu)
dens

## [1] 0.0003309269

　スケール行列の引数Sにlambda_ddの逆行列を指定します。

標準化t分布による計算

　続いて、確率変数の値を変換しておき、標準化t分布により確率密度を計算します。スケール行列と固有値・固有ベクトルの関係については「【R】2.3.0：分散共分散行列と固有値・固有ベクトルの関係の計算【PRMLのノート】 - からっぽのしょこ」を参照してください。

　スケール行列の固有値と固有ベクトルを計算します。

# 固有値・固有ベクトルを計算
res_eigen <- eigen(sigma_dd)
lambda_d <- res_eigen[["values"]]
u_dd <- res_eigen[["vectors"]] |> 
  t()
lambda_d; u_dd

## [1] 10.1423395  3.6188268  0.2388337
##           [,1]       [,2]       [,3]
## [1,] 0.2692509  0.9322270  0.2417783
## [2,] 0.9488715 -0.2138358 -0.2322006
## [3,] 0.1647628 -0.2919368  0.9421391

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

　確率変数の値を変換します。

# 変数を変換
y_d <- (u_dd %*% (x_d - mu_d) / sqrt(lambda_d)) |> 
  as.vector()
y_d

## [1]  0.3056598  0.8188797 -2.9833909

　$D$個の固有値$\lambda_d$を対角成分とする対角行列を$\mathbf{V}$(固有値の逆数$\frac{1}{\lambda_d}$を対角成分とする対角行列を$\mathbf{V}^{-1}$)、固有ベクトル$\mathbf{u}_d^{\top}$を行方向に並べた行列を$\mathbf{U}$として、次の式を計算します。

$$ \begin{aligned} \mathbf{y} &= \mathbf{V}^{-\frac{1}{2}} \mathbf{U} (\mathbf{x} - \boldsymbol{\mu}) \\ y_d &= \frac{1}{\sqrt{\lambda_d}} \mathbf{u}_d^{\top} (\mathbf{x} - \boldsymbol{\mu}) \end{aligned} $$

　$\mathbf{x}$の確率密度を計算します。

# Σを用いて標準化t分布により確率密度を計算
tmp_dens <- mvnfast::dmvt(X = y_d, mu = rep(0, times = D), sigma = diag(D), df = nu)
dens <- tmp_dens / sqrt(det(sigma_dd))
dens

## [1] 0.0003309269

　標準化t分布に従う$\mathbf{y}$の確率密度を計算して、$\boldsymbol{\Sigma}$の行列式の平方根の逆数を掛けます。

$$ \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \frac{ 1 }{ |\boldsymbol{\Sigma}|^{\frac{1}{2}} } \mathrm{St}(\mathbf{y} | \nu) $$

　標準化t分布に従う$\mathbf{y}$の確率密度を計算して、$\boldsymbol{\Sigma}$の行列式の平方根の逆数を掛けます。
　標準化t分布の計算は、dmvt()の位置ベクトルの引数muに0ベクトル、スケール行列の引数sigmaに単位行列を指定して計算できます。

　逆スケール行列を用いる場合は、次のように計算します。

# Λを用いて標準化t分布により確率密度を計算
tmp_dens <- mvnfast::dmvt(X = y_d, mu = rep(0, times = D), sigma = diag(D), df = nu)
dens <- sqrt(det(lambda_dd)) * tmp_dens
dens

## [1] 0.0003309269

　標準化t分布に従う$\mathbf{y}$の確率密度を計算して、$\boldsymbol{\Lambda}$の行列式の平方根を掛けます。

$$ \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Lambda}) = |\boldsymbol{\Lambda}|^{\frac{1}{2}} \mathrm{St}(\mathbf{y} | \nu) $$

統計量の計算

　最後に、多次元スチューデントのt分布の期待値・共分散・最頻値を計算します。

　期待値を計算します。

# 期待値を計算:(ν > 1)
E_x_d <- mu_d
E_x_d

## [1] 10.0 -6.0  1.5

　t分布の期待値は$\boldsymbol{\mu}$です。ただし、$\nu > 1$の場合に定義されます。

$$ \mathbb{E}[\mathbf{x}] = \boldsymbol{\mu} \quad (\nu > 1) $$

　共分散を計算します。

# Σを使って共分散を計算:(ν > 2)
cov_x_dd <- nu / (nu - 2) * sigma_dd
cov_x_dd

# Λを使って共分散を計算:(ν > 2)
cov_x_dd <- nu / (nu - 2) * solve(lambda_dd)
cov_x_dd

##            [,1] [,2]       [,3]
## [1,]  6.6666667    3 -0.1666667
## [2,]  3.0000000   15  4.0000000
## [3,] -0.1666667    4  1.6666667
##            [,1] [,2]       [,3]
## [1,]  6.6666667    3 -0.1666667
## [2,]  3.0000000   15  4.0000000
## [3,] -0.1666667    4  1.6666667

　t分布の分散は、$\boldsymbol{\Sigma}$または$\boldsymbol{\Lambda}$を使って、次の式で計算できます。ただし、$\nu > 2$の場合に定義されます。

$$ \mathrm{Cov}[\mathbf{x}] = \frac{\nu}{\nu - 2} \boldsymbol{\Sigma} = \frac{\nu}{\nu - 2} \boldsymbol{\Lambda}^{-1} \quad (\nu > 2) $$

　$\boldsymbol{\Sigma} = \boldsymbol{\Lambda}^{-1}$が成り立つのが分かります。

　最頻値を計算します。

# 最頻値を計算
mode_x_d <- mu_d
mode_x_d

## [1] 10.0 -6.0  1.5

　t分布の最頻値は$\boldsymbol{\mu}$です。

$$ \mathrm{mode}[\mathbf{x}] = \boldsymbol{\mu} $$

　この記事では、多次元t分布の計算を確認しました。次は、グラフを作成します。

参考文献

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

おわりに

　標準化t分布による計算の件はややこしいだけのような気もしますが、一応そうなるのを確認するためにやっておきました。1次元のときは関数の仕様で便利な面もありましたが、多次元でも逆行列の計算を回避できて嬉しいみたいなことがありますかね。

【次の内容】

からっぽのしょこ

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

【R】多次元スチューデントのt分布の計算

はじめに

多次元スチューデントのt分布の計算

確率密度の計算

スケール行列を使用

パラメータの設定

スクラッチで計算

関数で計算

逆スケール行列を使用

パラメータの設定

スクラッチで計算

関数で計算

標準化t分布による計算

統計量の計算

参考文献

おわりに