からっぽのしょこ

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

トップ > 確率分布 > 【R】ディリクレ分布の計算

【R】ディリクレ分布の計算

確率分布 確率分布-ディリクレ分布 R

はじめに

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

　この記事では、R言語でディリクレ分布に関する計算をします。

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

ディリクレ分布の計算

　ディリクレ分布(Dirichlet Distribution)の確率密度と統計量を計算します。ディリクレ分布については「ディリクレ分布の定義式 - からっぽのしょこ」を参照してください。

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

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

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

確率密度の計算

　ディリクレ分布に従う確率密度を計算する方法をいくつか確認します。

パラメータの設定

　ディリクレ分布のパラメータ$\boldsymbol{\beta}$、確率変数の実現値$\boldsymbol{\phi}$を設定します。

# パラメータを指定
beta_v <- c(4, 2, 3)

# 確率変数の値を指定
phi_v <- c(0.5, 0.3, 0.2)

　$V$次元ベクトル$\boldsymbol{\beta} = (\beta_1, \beta_2, \cdots, \beta_V)$、$\boldsymbol{\phi} = (\phi_1, \phi_2, \cdots, \phi_V)$を指定します。
　パラメータの各要素は$\beta_v > 0$を満たす必要があり、確率変数の各要素は$0 < \phi_v < 1$、$\sum_{v=1}^V \phi_v = 1$をとります。

スクラッチで計算

　定義式から計算します。

# 定義式により確率密度を計算
C    <- gamma(sum(beta_v)) / prod(gamma(beta_v))
dens <- C * prod(phi_v^(beta_v - 1))
dens

## [1] 5.04

　ディリクレ分布は、次の式で定義されます。

$$ \begin{aligned} C_{\mathrm{Dir}} &= \frac{ \Gamma(\sum_{v=1}^V \beta_v) }{ \prod_{v=1}^V \Gamma(\beta_v) } \\ \mathrm{Dir}(\boldsymbol{\phi} | \boldsymbol{\beta}) &= C_{\mathrm{Dir}} \prod_{v=1}^V \phi_v^{\beta_v-1} \end{aligned} $$

　ここで、$C_{\mathrm{Dir}}$はディリクレ分布の正規化係数、$\Gamma(x)$はガンマ関数です。
　ガンマ関数はgamma()で計算できます。

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

# 対数をとった定義式により確率密度を計算
log_C    <- lgamma(sum(beta_v)) - sum(lgamma(beta_v))
log_dens <- log_C + sum((beta_v - 1) * log(phi_v))
dens     <- exp(log_dens)
dens; log_dens

## [1] 5.04
## [1] 1.617406

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

$$ \begin{aligned} \log C_{\mathrm{Dir}} &= \log \Gamma \Bigl( \sum_{v=1}^V \beta_v \Bigr) - \sum_{v=1}^V \log \Gamma(\beta_v) \\ \log \mathrm{Dir}(\boldsymbol{\phi} | \boldsymbol{\beta}) &= \log C_{\mathrm{Dir}} + \sum_{v=1}^V (\beta_v - 1 )\log \phi_v \end{aligned} $$

　対数をとったガンマ関数は、lgamma()で計算できます。引数の値が大きいとgamma()の計算結果が発散してしまいます。その場合でも、lgamma()で計算できます。
　計算結果の指数をとると確率密度が得られます。

$$ \mathrm{Dir}(\boldsymbol{\phi} | \boldsymbol{\beta}) = \exp \Bigr( \log \mathrm{Dir}(\boldsymbol{\phi} | \boldsymbol{\beta}) \Bigr) $$

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

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

関数で計算

　MCMCpackパッケージのディリクレ分布の確率密度関数ddirichlet()で計算します。

# 関数により確率密度を計算
dens <- MCMCpack::ddirichlet(x = phi_v, alpha = beta_v)
dens

## [1] 5.04

　確率変数の引数xにphi_v、パラメータの引数alphaにbeta_vを指定します。

統計量の計算

　次は、ディリクレ分布の統計量を計算します。詳しくは「統計量の導出」を参照してください。

　ディリクレ分布のパラメータ$\boldsymbol{\beta}$と次元数$V$を設定します。

# パラメータを指定
beta_v <- c(4, 2, 3)

# 次元数を設定
V <- length(beta_v)


　期待値を計算します。

# 期待値を計算
E_phi_v <- beta_v / sum(beta_v)
E_phi_v

## [1] 0.4444444 0.2222222 0.3333333

　ディリクレ分布の期待値は、次の式で計算できます。

$$ \mathbb{E}[\phi_v] = \frac{\beta_v}{\sum_{v'=1}^V \beta_{v'}} $$

　対数の期待値を計算します。

# 対数の期待値を計算
E_log_phi_v <- digamma(beta_v) - digamma(sum(beta_v))
E_log_phi_v

## [1] -0.8845238 -1.7178571 -1.2178571

　対数をとった変数の期待値は、次の式で計算できます。

$$ \mathbb{E}[\log \phi_v] = \Psi(\beta_v) - \Psi \Bigl( \sum_{v=1}^V \beta_v \Bigr) $$

　ここで、$\Psi(x)$はディガンマ関数で、digamma()で計算できます。

　分散を計算します。

# 分散を計算
V_phi_v <- beta_v * (sum(beta_v) - beta_v) / sum(beta_v)^2 / (sum(beta_v) + 1)
V_phi_v

## [1] 0.02469136 0.01728395 0.02222222

　ディリクレ分布の分散は、次の式で計算できます。

$$ \mathbb{V}[\phi_v] = \frac{ \beta_v \left\{ \left( \sum_{v=1}^V \beta_v \right) - \beta_v \right\} }{ \left( \sum_{v=1}^V \beta_v \right)^2 \left\{ \left( \sum_{v=1}^V \beta_v \right) + 1 \right\} } $$

　要素を指定して、共分散を計算します。

# インデックスを指定:(i ≠ j)
i <- 1
j <- 2

# 共分散を計算
Cov_phi_ij <- -beta_v[i] * beta_v[j] / sum(beta_v)^2 / (sum(beta_v) + 1)
Cov_phi_ij

## [1] -0.009876543

　ディリクレ分布の共分散は、次の式で計算できます。

$$ \mathrm{Cov}[\phi_i, \phi_j] = \frac{ \beta_i \beta_j }{ \left( \sum_{v=1}^V \beta_v \right)^2 \left\{ \left( \sum_{v=1}^V \beta_v \right) + 1 \right\} } \qquad (i \neq j) $$

　全ての要素の組み合わせで、共分散を計算します。

# 共分散を計算
Cov_phi_vv <- -beta_v %*% t(beta_v) / sum(beta_v)^2 / (sum(beta_v) + 1)
diag(Cov_phi_vv) <- as.numeric(NA)
Cov_phi_vv

##              [,1]         [,2]         [,3]
## [1,]           NA -0.009876543 -0.014814815
## [2,] -0.009876543           NA -0.007407407
## [3,] -0.014814815 -0.007407407           NA

　ただし、$i = j$の場合は計算できないので、diag()を使って欠損値NAに置き換えます。

　最頻値(モード)を計算します。

# 最頻値を計算:(β_v > 1)
mode_phi_v <- (beta_v - 1) / (sum(beta_v) - V)
mode_phi_v

## [1] 0.5000000 0.1666667 0.3333333

　ディリクレ分布の最頻値は、次の式で計算できます。

$$ \mathrm{mode}[\phi_v] = \frac{ \beta_v - 1 }{ \left( \sum_{v=1}^V \beta_v \right) - V } \qquad (\beta_v > 1) $$


　この記事では、ディリクレ分布の計算を確認しました。次は、グラフを作成します。

参考文献

  • 岩田具治『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社,2015年.

おわりに

　飛ばしていたディリクレ分布をやっていきます。

【次の内容】

www.anarchive-beta.com