からっぽのしょこ

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

トップ > 確率分布 > 【R】多次元ガウス分布の計算

【R】多次元ガウス分布の計算

確率分布 確率分布-ガウス分布(正規分布) R

はじめに

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

　この記事では、R言語で多次元ガウス分布(多変量正規分布)に関する計算をします。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

多次元ガウス分布の計算

　多次元ガウス分布(Multivariate Gaussian Distribution)・多変量正規分布(Multivariate Normal Distribution)の確率密度を計算します。多次元ガウス分布については「多次元ガウス分布の定義式 - からっぽのしょこ」を参照してください。

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

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

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

確率密度の計算

　多次元ガウス分布に従う確率密度を計算する方法をいくつか確認します。

分散共分散行列を使用

　まずは、分散共分散行列を用いてガウス分布を計算します。

パラメータの設定

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

# 次元数を指定
D <- 3

# 平均ベクトルを指定
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)

　平均ベクトル$\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^2 & \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^2 \end{pmatrix} ,\ \mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_D \end{pmatrix} $$

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

スクラッチで計算

　定義式から計算します。

# 定義式により確率密度を計算
C_N  <- 1 / sqrt((2 * pi)^D * det(sigma_dd))
dens <- C_N * exp(-0.5 * t(x_d - mu_d) %*% solve(sigma_dd) %*% (x_d - mu_d)) |> 
  as.numeric()
dens

## [1] 0.0001708777

　ガウス分布は、次の式で定義されます。

$$ \begin{aligned} C_{\mathcal{N}} &= \frac{1}{\sqrt{(2 \pi)^D |\boldsymbol{\Sigma}|}} \\ \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}, \boldsymbol{\Sigma}) &= C_{\mathcal{N}} \exp \left\{ - \frac{1}{2} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) \right\} \end{aligned} $$

　ここで、$C_{\mathcal{N}}$はガウス分布の正規化係数、$\pi$は円周率、$\mathbf{A}^{\top}$は転置行列、$|\mathbf{A}|$は行列式、$\mathbf{A}^{-1}$は逆行列です。
　円周率はpi、転置はt()、行列式はdet()、逆行列solve()、行列の積は%*%演算子で計算できます。

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

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

## [1] 0.0001708777
## [1] -8.674563

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

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

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

$$ \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \exp \Bigr( \log \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}, \boldsymbol{\Sigma}) \Bigr) $$

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

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

関数で計算

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

# 多次元ガウス分布の関数により確率密度を計算
dens <- mvnfast::dmvn(X = x_d, mu = mu_d, sigma = sigma_dd)
dens

## [1] 0.0001708777

　確率変数の引数Xにx_d、平均の引数muにmu_d、分散共分散行列の引数sigmaにsigma_ddを指定します。

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

# 多次元ガウス分布の対数をとった関数により確率密度を計算
log_dens <- mvnfast::dmvn(X = x_d, mu = mu_d, sigma = sigma_dd, log = TRUE)
dens <- exp(log_dens)
dens; log_dens

## [1] 0.0001708777
## [1] -8.674563

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

精度行列を使用

　続いて、精度行列を用いてガウス分布を計算します。

パラメータの設定

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

# 次元数を指定
D <- 3

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

# 精度行列を指定
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) # マトリクスに変換

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

　分散共分散行列$\boldsymbol{\Sigma}$の代わりに、精度行列$\boldsymbol{\Lambda}$を指定します。

　あるいは、$\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_N  <- sqrt(det(lambda_dd) / (2 * pi)^D)
dens <- C_N * exp(-0.5 * t(x_d - mu_d) %*% lambda_dd %*% (x_d - mu_d)) |> 
  as.numeric()
dens

## [1] 0.0001708777

　ガウス分布は、次の式で定義されます。

$$ \begin{aligned} C_{\mathcal{N}} &= \left\{ \frac{|\boldsymbol{\Lambda}|}{(2 \pi)^D} \right\}^{\frac{1}{2}} \\ \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1}) &= C_{\mathcal{N}} \exp \left\{ - \frac{1}{2} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Lambda} (\mathbf{x} - \boldsymbol{\mu}) \right\} \end{aligned} $$

　平方根の性質より$\sqrt{x} = x^{\frac{1}{2}}$、$\sqrt{\frac{x}{y}} = \frac{\sqrt{x}}{\sqrt{y}}$です。

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

# 対数をとった定義式により確率密度を計算
log_C_N  <- -0.5 * (D * log(2 * pi) - log(det(lambda_dd)))
log_dens <- log_C_N - 0.5 * t(x_d - mu_d) %*% lambda_dd %*% (x_d - mu_d) |> 
  as.numeric()
dens <- exp(log_dens)
dens; log_dens

## [1] 0.0001708777
## [1] -8.674563

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

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

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

関数で計算

　ガウス分布の確率密度関数dmvn()で計算します。

# 多次元ガウス分布の関数により確率密度を計算
dens <- mvnfast::dmvn(X = x_d, mu = mu_d, sigma = solve(lambda_dd))
dens

# 多次元ガウス分布の対数をとった関数により確率密度を計算
log_dens <- mvnfast::dmvn(X = x_d, mu = mu_d, sigma = solve(lambda_dd), log = TRUE)
dens <- exp(log_dens)
dens; log_dens

## [1] 0.0001708777
## [1] 0.0001708777
## [1] -8.674563

　分散共分散行列の引数sigmaにlambda_ddの逆行列を指定します。

統計量の計算

　ガウス分布の分散共分散行列に関する統計量の計算については

www.anarchive-beta.com

を参照してください。

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

参考文献

  • 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.

おわりに

　書けば書くほど書きたいことが増えていく。

【次の内容】

www.anarchive-beta.com