からっぽのしょこ

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

トップ > 確率分布 > 【R】ガウス-ガンマ分布の計算

【R】ガウス-ガンマ分布の計算

確率分布 確率分布-ガウス・ガンマ分布 R

はじめに

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

　この記事では、ガウス-ガンマ分布(正規-ガンマ分布)に関する計算をします。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

ガウス-ガンマ分布の計算

　ガウス-ガンマ分布(Gaussian-Gamma Distribution)または正規-ガンマ分布(Normal-Gamma Distribution)の確率密度と統計量を計算します。ガウス-ガンマ分布については「確率分布ネタの記事まとめ - からっぽのしょこ」を参照してください。

確率密度の計算

　ガウス-ガンマ分布に従う確率密度を計算する方法をいくつか確認します。1次元ガウス分布の計算については「【R】1次元ガウス分布の計算 - からっぽのしょこ」、ガンマ分布の計算については「【R】ガンマ分布の計算 - からっぽのしょこ」を参照してください。

パラメータの設定

　まずは、ガウス分布のパラメータ$m, \beta$と確率変数の実現値$\mu$、ガンマ分布のパラメータ$a, b$と確率変数の実現値$\lambda$を設定します。

# 1次元ガウス分布の平均パラメータを指定
m <- 0

# 1次元ガウス分布の精度パラメータの係数を指定
beta <- 2

# ガンマ分布のパラメータを指定
a <- 5
b <- 6

# 確率変数の値を指定
mu     <- 1.5
lambda <- 2.5

　ガウス分布の平均パラメータ$m$と精度パラメータの係数$\beta > 0$、ガンマ分布の形状パラメータ$a > 0$と尺度パラメータ$b > 0$、確率変数がとり得る値$\mu, \lambda > 0$を指定します。設定した値に従う確率密度を計算します。

スクラッチで計算

　定義式から計算します。

# 定義式により確率密度を計算
C_N      <- sqrt(beta * lambda / 2 / pi)
dens_N   <- C_N * exp(-0.5 * beta * lambda * (mu - m)^2)
C_Gam    <- b^a / gamma(a)
dens_Gam <- C_Gam * lambda^(a - 1) * exp(-b * lambda)
dens     <- dens_N * dens_Gam
dens; dens_N; dens_Gam

## [1] 1.245594e-05
## [1] 0.003217278
## [1] 0.003871576

　1次元ガウス分布

$$ \begin{aligned} C_{\mathcal{N}} &= \Bigl( \frac{\beta \lambda}{2 \pi} \Bigr)^{\frac{1}{2}} \\ \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) &= C_{\mathcal{N}} \exp \left( - \frac{\beta \lambda}{2} (\mu - m)^2 \right) \end{aligned} $$

と、ガンマ分布

$$ \begin{aligned} C_{\mathrm{Gam}} &= \frac{b^a}{\Gamma(a)} \\ \mathrm{Gam}(\lambda | a, b) &= C_{\mathrm{Gam}} \lambda^{a-1} \exp(- b \lambda) \end{aligned} $$

の積により、ガウス-ガンマ分布が定義されます。

$$ \mathrm{NG}(\mu, \lambda | m, \beta, a, b) = \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) \mathrm{Gam}(\lambda | a, b) $$

　ここで、$C_{\mathcal{N}}$はガウス分布の正規化係数、$C_{\mathrm{Gam}}$はガンマ分布の正規化係数、$\pi$は円周率、$\Gamma(x)$はガンマ関数です。
　円周率はpi、ガンマ関数はgamma()で計算できます。

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

# 対数をとった定義式により確率密度を計算
log_C_N      <- 0.5 * (log(beta * lambda) - log(2 * pi))
log_dens_N   <- log_C_N - 0.5 * beta * lambda * (mu - m)^2
log_C_Gam    <- a * log(b) - lgamma(a)
log_dens_Gam <- log_C_Gam + (a - 1) * log(lambda) - b * lambda
log_dens     <- log_dens_N + log_dens_Gam
dens         <- exp(log_dens)
dens; log_dens

## [1] 1.245594e-05
## [1] -11.29331

　対数をとった1次元ガウス分布

$$ \begin{aligned} \log C_{\mathcal{N}} &= \frac{1}{2} \Bigl( \log (\beta \lambda) - \log (2 \pi) \Bigr) \\ \log \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) &= \log C_{\mathcal{N}} - \frac{\beta \lambda}{2} (\mu - m)^2 \end{aligned} $$

と、対数をとったガンマ分布

$$ \begin{aligned} \log C_{\mathrm{Gam}} &= a \log b - \log \Gamma(a) \\ \log \mathrm{Gam}(\lambda | a, b) &= \log C_{\mathrm{Gam}} + (a - 1) \log \lambda - b \lambda \end{aligned} $$

の和を計算します。

$$ \log \mathrm{NG}(\mu, \lambda | m, \beta, a, b) = \log \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) + \log \mathrm{Gam}(\lambda | a, b) $$

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

$$ \mathrm{NG}(\mu, \lambda | m, \beta, a, b) = \exp \Bigr( \log \mathrm{NG}(\mu, \lambda | m, \beta, a, b) \Bigr) $$

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

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

関数による計算

　ガウス分布の確率密度関数dnorm()とガンマ関数の確率密度関数gamma()で計算します。

# 関数により確率密度を計算
dens_N   <- dnorm(x = mu, mean = m, sd = 1/sqrt(beta*lambda))
dens_Gam <- dgamma(x = lambda, shape = a, rate = b)
dens     <- dens_N * dens_Gam
dens

## [1] 1.245594e-05

　dnorm()の確率変数の引数xにmu、平均の引数meanにm、標準偏差の引数sdにbetaとlambdaの積の平方根の逆数を指定して、ガウス分布の確率密度を計算します。平方根はsqrt()で計算できます。
　dgamma()の確率変数の引数xにlambda、形状の引数shapeにa、尺度の引数rateにbを指定して、ガンマ分布の確率密度を計算します。
　2つの確率密度の積を計算します。

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

# 対数をとった関数により確率密度を計算
log_dens_N   <- dnorm(x = mu, mean = m, sd = 1/sqrt(beta*lambda), log = TRUE)
log_dens_Gam <- dgamma(x = lambda, shape = a, rate = b, log = TRUE)
log_dens     <- log_dens_N + log_dens_Gam
dens         <- exp(log_dens)
dens; log_dens

## [1] 1.245594e-05
## [1] -11.29331

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

統計量の計算

　次は、ガウス-ガンマ分布の期待値・分散・最頻値を計算します。詳しくはいつか書きます。

　$\mu, \lambda$それぞれの期待値を計算します。

# 期待値を計算
E_mu 　　<- m
E_lambda <- a / b
E_mu; E_lambda

## [1] 0
## [1] 0.8333333

　$\mu$の期待値は$m$で、$\lambda$の期待値は次の式で計算できます。

$$ \begin{aligned} \mathbb{E}[\mu] &= m \\ \mathbb{E}[\lambda] &= \frac{a}{b} \end{aligned} $$

　分散を計算します。

# 分散を計算
V_mu 　　<- 1 / beta / E_lambda
V_lambda <- a / b^2
V_mu; V_lambda

## [1] 0.6
## [1] 0.1388889

　$\mu$の分散と$\lambda$の分散は、それぞれ次の式で計算できます。

$$ \begin{aligned} \mathbb{V}[\mu] &= \frac{1}{\beta \lambda} \\ \mathbb{V}[\lambda] &= \frac{a}{b^2} \end{aligned} $$

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

# 最頻値を計算:(a >= 1)
mode_mu     <- m
mode_lambda <- (a - 1) / b
mode_mu; mode_lambda

## [1] 0
## [1] 0.6666667

　$\mu$の最頻値は$m$で、$\lambda$の最頻値は次の式で計算できます。

$$ \begin{aligned} \mathrm{mode}[\mu] &= m \\ \mathrm{mode}[\lambda] &= \frac{a - 1}{b} \quad (a \geq 1) \end{aligned} $$

　ガンマ分布の確率変数は$\lambda > 0$の値をとるので、$a \geq 1$のとき最頻値が定義されます($a < 1$だと計算結果が負の値になってしまいます)。

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

参考文献

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

おわりに

　加筆修正の際に「ガウス-ガンマの作図」から記事を分割しました。

【次の内容】

www.anarchive-beta.com