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

はじめに

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

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

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

ディリクレ分布の計算

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

　利用するライブラリを読み込みます。

# 利用ライブラリ
import numpy as np
from scipy.special import gamma, loggamma, digamma
from scipy.stats import dirichlet

確率密度の計算

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

パラメータの設定

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

# パラメータを指定
beta_v = np.array([0.4, 0.2, 0.3])

# 確率変数の値を指定
phi_v = np.array([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 \leq \phi_v \leq 1$、$\sum_{v=1}^V \phi_v = 1$をとります。

スクラッチで計算

　定義式から計算します。

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

0.4297763573774013

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

$$ \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)$はガンマ関数です。
　ガンマ関数はSciPyライブラリのspecialモジュールのgamma()で計算できます。

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

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

0.42977635737740116
-0.8444903047153147

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

$$ \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} $$

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

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

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

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

モジュールで計算

　SciPyライブラリのstatsモジュールのdirichletの確率密度関数pdf()で計算します。

# 関数により確率密度を計算
dens = dirichlet.pdf(x=phi_v, alpha=beta_v)
print(dens)

0.42977635737740116

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

　logpdf()だと対数をとった確率密度を計算します。

# 対数をとった関数により確率密度を計算
log_dens = dirichlet.logpdf(x=phi_v, alpha=beta_v)
dens = np.exp(log_dens)
print(dens)
print(log_dens)

0.42977635737740116
-0.8444903047153147

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

　全てのパラメータが1未満で確率変数に0を含むとエラーになるようです。

# パラメータを指定
beta_v = np.array([0.4, 0.2, 0.3])

# 確率変数の値を指定
phi_v = np.array([0.6, 0.4, 0.0])

# 関数により確率密度を計算
dens = dirichlet.pdf(x=phi_v, alpha=beta_v)
print(dens)

---------------------------------------------------------------------------

ValueError                                Traceback (most recent call last)

Input In [8], in <cell line: 8>()
      5 phi_v = np.array([0.6, 0.4, 0.0])
      7 # 関数により確率密度を計算
----> 8 dens = dirichlet.pdf(x=phi_v, alpha=beta_v)
      9 print(dens)
(省略)

ValueError: Each entry in 'x' must be greater than zero if its alpha is less than one.

　1以上の値だとエラーになりません。

# パラメータを指定
beta_v = np.array([0.4, 0.2, 3.0])

# 確率変数の値を指定
phi_v = np.array([0.6, 0.4, 0.0])

# 関数により確率密度を計算
dens = dirichlet.pdf(x=phi_v, alpha=beta_v)
print(dens)

0.0

　なぜこんな挙動を？

統計量の計算

　次は、ディリクレ分布の統計量を計算します。詳しくは「1.2.4：ディリクレ分布【『トピックモデル』の勉強ノート】 - からっぽのしょこ」を参照してください。

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

# パラメータを指定
beta_v = np.array([4.0, 2.0, 3.0])

# 次元数を設定
V = len(beta_v)

　期待値を計算します。

# 期待値を計算
E_phi_v = beta_v / np.sum(beta_v)
print(E_phi_v)

[0.44444444 0.22222222 0.33333333]

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

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

　分散を計算します。

# 分散を計算
V_phi_v = beta_v * (np.sum(beta_v) - beta_v)
V_phi_v /= np.sum(beta_v)**2 * (np.sum(beta_v) + 1.0)
print(V_phi_v)

[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 = 0
j = 1

# 共分散を計算
Cov_phi_ij = -beta_v[i] * beta_v[j]
Cov_phi_ij /= np.sum(beta_v)**2 * (np.sum(beta_v) + 1.0)
print(Cov_phi_ij)

-0.009876543209876543

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

$$ \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 = -np.dot(beta_v.reshape(V, 1), beta_v.reshape(1, V))
Cov_phi_vv /= np.sum(beta_v)**2 * (np.sum(beta_v) + 1.0)
Cov_phi_vv[np.arange(V), np.arange(V)] = np.nan
print(Cov_phi_vv)

[[        nan -0.00987654 -0.01481481]
 [-0.00987654         nan -0.00740741]
 [-0.01481481 -0.00740741         nan]]

　ただし、$i = j$の場合は計算できないので、対角要素を欠損値np.nanに置き換えます。

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

# 最頻値を計算:(β_v > 1)
mode_phi_v = (beta_v - 1.0) / (np.sum(beta_v) - V)
print(mode_phi_v)

[0.5        0.16666667 0.33333333]

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

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

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

# 対数の期待値を計算
E_log_phi_v = digamma(beta_v) - digamma(np.sum(beta_v))
print(E_log_phi_v)

[-0.88452381 -1.71785714 -1.21785714]

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

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

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

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

参考文献

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

おわりに

　謎挙動がありましたが深入りしないことにします。

【次の内容】

www.anarchive-beta.com