はじめに

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

　この記事では、R言語でガンマ分布に関する計算をします。

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

はじめに
ガンマ分布の計算
参考文献
おわりに

ガンマ分布の計算

　ガンマ分布(Gamma Distribution)の確率密度と統計量を計算します。ガンマ分布については「【R】ガンマ分布の作図 - からっぽのしょこ」を参照してください。

確率密度の計算

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

パラメータの設定

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

# パラメータを指定
a <- 5
b <- 2

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

　2つのパラメータ$a > 0, b > 0$、正の実数$\lambda > 0$を指定します。設定した値に従う確率密度を計算します。

スクラッチで計算

　定義式から計算します。

# 定義式により確率密度を計算
C    <- b^a / gamma(a)
dens <- C * lambda^(a - 1) * exp(-b * lambda)
dens

## [1] 0.3509347

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

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

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

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

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

## [1] 0.3509347
## [1] -1.047155

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

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

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

$$ \mathrm{Gam}(\lambda | a, b) = \exp \Bigr( \log \mathrm{Gam}(\lambda | a, b) \Bigr) $$

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

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

ガンマ分布の関数による計算

　ガンマ分布の確率密度関数dgamma()で計算します。

# ガンマ分布の関数により確率密度を計算
dens <- dgamma(x = lambda, shape = a, rate = b)
dens

## [1] 0.3509347

　確率変数の引数xにlambda、形状の引数shapeにa、尺度の引数rateにbを指定します。

　rate引数の代わりに、scale引数でも計算できます。

# ガンマ分布の関数により確率密度を計算
dens <- dgamma(x = lambda, shape = a, scale = 1/b)
dens

## [1] 0.3509347

　尺度の引数scaleに$b$の逆数1/bを指定します。

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

# ガンマ分布の対数をとった関数により確率密度を計算
log_dens <- dgamma(x = lambda, shape = a, rate = b, log = TRUE)
dens     <- exp(log_dens)
dens; log_dens

## [1] 0.3509347
## [1] -1.047155

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

　scaleの場合も同じです。

# ガンマ分布の対数をとった関数により確率密度を計算
log_dens <- dgamma(x = lambda, shape = a, scale = 1 / b, log = TRUE)
dens     <- exp(log_dens)
dens; log_dens

## [1] 0.3509347
## [1] -1.047155

統計量の計算

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

　期待値を計算します。

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

## [1] 2.5

　ガンマ分布の期待値は、次の式で計算できます。

$$ \mathbb{E}[\lambda] = \frac{a}{b} $$

　分散を計算します。

# 分散を計算
V_lambda <- a / b^2
V_lambda

## [1] 1.25

　ガンマ分布の分散は、次の式で計算できます。

$$ \mathbb{V}[\lambda] = \frac{a}{b^2} $$

　最頻値を計算します。

# 最頻値を計算
mode_lambda <- (a - 1) / b
mode_lambda

## [1] 2

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

$$ \mathrm{mode}[\lambda] = \frac{a - 1}{b} $$

歪度と尖度の計算

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

　歪度を計算します。

# 歪度を計算
skewness <- 2 / sqrt(a)
skewness

## [1] 0.8944272

　ガンマ分布の歪度は、次の式で計算できます。

$$ \mathrm{Skewness} = \frac{\mathbb{E}[(\lambda - \mu)^3]}{\sigma^3} = \frac{2}{\sqrt{a}} $$

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

　尖度を計算します。

# 尖度を計算
kurtosis <- 6 / a
kurtosis

## [1] 1.2

　ガンマ分布の尖度は、次の式で計算できます。

$$ \mathrm{Kurtosis} = \frac{\mathbb{E}[(\lambda - \mu)^4]}{\sigma^4} - 3 = \frac{6}{a} $$

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

参考文献

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

おわりに

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

【次の内容】

からっぽのしょこ

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

【R】ガンマ分布の計算