【R】ベータ分布の計算

はじめに

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

　この記事では、R言語でベータ分布の計算を行います。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

ベータ分布の計算

　ベータ分布(Beta Distribution)の確率密度と統計量を計算します。

確率密度の計算

　ベータ分布に従う確率密度を計算する方法をいくつか確認します。

パラメータの設定

　まずは、ベータ分布のパラメータ$\alpha, \beta$と確率変数の実現値$\phi$を設定します。

# パラメータを指定
alpha <- 5
beta  <- 2

# 確率変数の値を指定
phi <- 0.6

　2つのパラメータ$\alpha > 0, \beta > 0$、0から1の値$0 < \phi < 1$を指定します。設定した値に従う確率密度を計算します。$\phi$がとり得る範囲外の値の場合は、確率密度が0になります。

スクラッチで計算

　定義式から計算します。

# 定義式により確率密度を計算
C    <- gamma(alpha + beta) / gamma(alpha) / gamma(beta)
dens <- C * phi^(alpha - 1) * (1 - phi)^(beta - 1)
dens

## [1] 1.5552

　ベータ分布は、次の式で定義されます。

$$ \begin{aligned} C_{\mathrm{Beta}} &= \frac{ \Gamma(\alpha + \beta) }{ \Gamma(\alpha) \Gamma(\beta) } \\ \mathrm{Beta}(\phi | \alpha, \beta) &= C_{\mathrm{Beta}} \phi^{\alpha-1} (1 - \phi)^{\beta-1} \end{aligned} $$

　ここで、$C_{\mathrm{Beta}}$はベータ分布の正規化係数、$\Gamma(x)$はガンマ関数です。
　ガンマ関数は、gamma()で計算できます。

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

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

## [1] 1.5552
## [1] 0.4416042

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

$$ \begin{aligned} \log C_{\mathrm{Beta}} &= \log \Gamma(\alpha + \beta) - \log \Gamma(\alpha) - \log \Gamma(\beta) \\ \log \mathrm{Beta}(\phi | \alpha, \beta) &= \log C_{\mathrm{Beta}} + (\alpha - 1) \log \phi + (\beta - 1) \log (1 - \phi) \end{aligned} $$

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

$$ \mathrm{Beta}(\phi | \alpha, \beta) = \exp \Bigr( \log \mathrm{Beta}(\phi | \alpha, \beta) \Bigr) $$

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

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

ベータ分布の関数による計算

　ベータ分布の確率密度関数dbeta()で計算します。

# ベータ分布の関数により確率密度を計算
dens <- dbeta(x = phi, shape1 = alpha, shape2 = beta)
dens

## [1] 1.5552

　確率変数の引数xにphi、パラメータの引数shape1, shape2にalpha, betaを指定します。

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

# ベータ分布の対数をとった関数により確率密度を計算
log_dens <- dbeta(x = phi, shape1 = alpha, shape2 = beta, log = TRUE)
dens     <- exp(log_dens)
dens; log_dens

## [1] 1.5552
## [1] 0.4416042

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

ディリクレ分布の関数による計算

　以降の計算は、変数とパラメータをベクトル形式にしておく必要があります。

# ベクトルに変換
alpha_v <- c(alpha, beta)
phi_v   <- c(phi, 1 - phi)
alpha_v; phi_v

## [1] 5 2
## [1] 0.6 0.4

　2つのパラメータをベクトルにまとめてalpha_vとします。
　成功確率(クラス1の出現確率)$\phi$と失敗確率(クラス0の出現確率)$1 - \phi$をベクトルにまとめてphi_vとします。
　「alphaとphi」「betaと1-phi」のインデックスが対応していれば、順番は計算結果に影響しません。

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

# ディリクレ分布の関数により確率密度を計算
dens <- MCMCpack::ddirichlet(x = phi_v, alpha = alpha_v)
dens

## [1] 1.5552

　確率変数の引数xにphi_v、パラメータの引数alphaにalpha_vを指定します。$V = 2$のときベータ分布の確率密度になります。

統計量の計算

　次は、ベータ分布の期待値・分散・最頻値を計算します。詳しくは「1.2.3：ベータ分布【『トピックモデル』の勉強ノート】 - からっぽのしょこ」を参照してください。

　期待値を計算します。

# 期待値を計算
E_phi <- alpha / (alpha + beta)
E_phi

## [1] 0.7142857

　ベータ分布の期待値は、次の式で計算できます。

$$ \mathbb{E}[\phi] = \frac{\alpha}{\alpha + \beta} $$

　分散を計算します。

# 分散を計算
V_phi <- alpha * beta / (alpha + beta)^2 / (alpha + beta + 1)
V_phi

## [1] 0.0255102

　ベータ分布の分散は、次の式で計算できます。

$$ \mathbb{V}[\phi] = \frac{ \alpha \beta }{ (\alpha + \beta)^2 (\alpha + \beta + 1) } $$

　最頻値を計算します。

# 最頻値を計算
mode_phi <- (alpha - 1) / (alpha + beta - 2)
mode_phi

## [1] 0.8

　ベータ分布の最頻値(モード)は、次の式で計算できます。

$$ \mathrm{mode}[\phi] = \frac{\alpha - 1}{\alpha + \beta - 2} $$

歪度と尖度の計算

　最後に、ベータ分布の歪度と尖度を計算します。詳しくはいつか書きます。

　歪度を計算します。

# 歪度を計算
numer <- 2 * (beta - alpha) * sqrt(alpha + beta + 1)
denom <- (alpha + beta + 2) * sqrt(alpha * beta)
skewness <- numer / denom
skewness

## [1] -0.5962848

　ベータ分布の歪度は、次の式で計算できます。

$$ \mathrm{Skewness} = \frac{\mathbb{E}[(\phi - \mu)^3]}{\sigma^3} = \frac{ 2 (\beta - \alpha) \sqrt{\alpha + \beta + 1} }{ (\alpha + \beta + 2) \sqrt{\alpha \beta} } $$

　ここで、$\phi$の平均$\mu = \mathbb{E}[\phi]$、標準偏差$\sigma = \sqrt{\mathbb{E}[(\phi - \mu)^2]}$です。

　尖度を計算します。

# 尖度を計算
numer <- (6 * (alpha - beta)^2 * (alpha + beta + 1) - alpha * beta * (alpha + beta + 2))
denom <- alpha * beta * (alpha + beta + 2) * (alpha + beta + 3)
kurtosis <- numer / denom
kurtosis

## [1] 0.38

　ベータ分布の尖度は、次の式で計算できます。

$$ \mathrm{Kurtosis} = \frac{\mathbb{E}[(\phi - \mu)^4]}{\sigma^4} - 3 = \frac{ 6 \Bigl\{ (\alpha - \beta)^2 (\alpha + \beta + 1) - \alpha \beta (\alpha + \beta + 2) \Bigr\} }{ (\alpha + \beta + 2) (\alpha + \beta + 3) } $$

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

参考文献

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

おわりに

　途中まで書いて放置してたのを続きも含めて書き上げられました。このまま他の分布も進めたい。

【次の内容】

www.anarchive-beta.com