からっぽのしょこ

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

【R】二項分布の計算

はじめに

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

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

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

二項分布の計算

 二項分布(Binomial Distribution)の確率と統計量を計算します。二項分布については「二項分布の定義式の確認 - からっぽのしょこ」を参照してください。

確率の計算

 二項分布に従う確率を計算する方法をいくつか確認します。

パラメータの設定

 まずは、二項分布のパラメータ$\phi, M$を設定します。

# パラメータを指定
phi <- 0.35

# 試行回数を指定
M <- 10

# 確率変数の値を指定:(x <= M)
x <- 3

 成功確率$0 < \phi < 1$、試行回数$M$、確率変数がとり得る値$x \in \{0, 1, \ldots, M\}$を指定します。設定した値に従う確率を計算します。

スクラッチで計算

 定義式から確率を計算します。

# 定義式により確率を計算
C <- gamma(M + 1) / gamma(M - x + 1) / gamma(x + 1)
prob <- C * phi^x * (1 - phi)^(M - x)
prob
## [1] 0.2522196

 二項分布の定義式

$$ \begin{aligned} C_{\mathrm{Bin}} &= \frac{M!}{(M - x)! x!} \\ \mathrm{Bin}(x | M, \phi) &= C_{\mathrm{Bin}} \phi^x (1 - \phi)^{M-x} \end{aligned} $$

で計算します。$C_{\mathrm{Bin}}$は、二項分布の正規化係数です。
 階乗$x!$の計算は、ガンマ関数$\Gamma(x + 1) = x!$に置き換えて計算します。ガンマ関数は、gamma()で計算できます。

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

# 対数をとった定義式により確率を計算
log_C <- lgamma(M + 1) - lgamma(M - x + 1) - lgamma(x + 1)
log_prob <- log_C + x * log(phi) + (M - x) * log(1 - phi)
prob <- exp(log_prob)
prob; log_prob
## [1] 0.2522196
## [1] -1.377455

 対数をとった定義式

$$ \begin{aligned} \log C_{\mathrm{Bin}} &= \log M! - \log (M - x)! - \log x! \\ \log \mathrm{Bin}(x | M, \phi) &= \log C_{\mathrm{Bin}} + x \log \phi + (N - x) \log (1 - \phi) \end{aligned} $$

を計算します。
 対数をとったガンマ関数$\log \Gamma(x)$は、lgamma()で計算できます。
 計算結果の指数をとると確率が得られます。

$$ \mathrm{Bin}(x | M, \phi) = \exp \Bigr( \log \mathrm{Bin}(x | M, \phi) \Bigr) $$

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

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

二項分布の関数による計算

 二項分布の確率関数dbinom()を使って計算します。

# 二項分布の関数により確率を計算
prob <- dbinom(x = x, size = M, prob = phi)
prob
## [1] 0.2522196

 成功回数の引数xx、試行回数の引数sizeM、成功確率の引数probphiを指定します。

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

# 二項分布の対数をとった関数により確率を計算
log_prob <- dbinom(x = x, size = M, prob = phi, log = TRUE)
prob <- exp(log_prob)
prob; log_prob
## [1] 0.2522196
## [1] -1.377455

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

多項分布の関数による計算

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

# ベクトルに変換
phi_v <- c(1 - phi, phi)
x_v <- c(M - x, x)
phi_v; x_v
## [1] 0.65 0.35
## [1] 7 3

 失敗確率(クラス0の出現確率)$1 - \phi$と成功確率(クラス1の出現確率)$\phi$を持つベクトルを作成してphi_vとします。
 失敗回数(クラス0の出現頻度)$M - x$と成功回数(クラス1の出現頻度)$x$を持つベクトルを作成してx_vとします。

 多項分布の確率関数dmultinom()を使って計算します。

# 多項分布の関数により確率を計算
prob <- dmultinom(x = x_v, size = M, prob = phi_v)
prob
## [1] 0.2522196

 出現頻度の引数xx_v、試行回数の引数sizeM、出現確率の引数probphi_vを指定します。

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

# 多項分布の対数をとった関数により確率を計算
log_prob <- dmultinom(x = x_v, size = M, prob = phi_v, log = TRUE)
prob <- exp(log_prob)
prob; log_prob
## [1] 0.2522196
## [1] -1.377455

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

統計量の計算

 二項分布の平均と分散を計算します。詳しくは「二項分布の平均と分散の導出:定義式を利用 - からっぽのしょこ」を参照してください。

 平均を計算します。

# 平均を計算
E_x <- M * phi
E_x
## [1] 3.5

 二項分布の平均は、次の式で計算できます。

$$ \mathbb{E}[x] = M \phi $$

 $M, \phi$が大きいほど平均が大きくなります。

 分散を計算します。

# 分散を計算
V_x <- M * phi * (1 - phi)
V_x
## [1] 2.275

 二項分布の分散は、次の式で計算できます。

$$ \mathbb{V}[x] = M \phi (1 - \phi) $$

 $M$が大きいまたは$\phi$が0.5に近いほど分散が大きくなります。

 最頻値を計算します。

# 最頻値を計算
mode_x <- floor(phi * (M + 1))
mode_x
## [1] 3

 二項分布の最頻値は、次の式で計算できます。

$$ \phi (M + 1) - 1 \leq \mathrm{mode}[x] \leq \phi (M + 1) $$

 $\phi (M + 1)$以下の最大の整数で、$\phi (M + 1)$が整数のときは$\phi (M + 1) - 1$も最頻値になります。
 floor()で小数点以下を切り捨てると最大の整数が得られます。ただし、$\phi (M + 1)$が整数のときに$\phi (M + 1) - 1$は得られません。

歪度と尖度の計算

 二項分布の歪度と尖度を計算します。詳しくはいつか書きます。

 歪度を計算します。

# 歪度を計算
skewness <- (1 - 2 * phi) / sqrt(M * phi * (1 - phi))
skewness
## [1] 0.1988981

 二項分布の歪度は、次の式で計算できます。

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

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

 尖度を計算します。

# 尖度を計算
kurtosis <- (1 - 6 * phi * (1 - phi)) / (M * phi * (1 - phi))
kurtosis
## [1] -0.1604396

 二項分布の尖度は、次の式で計算できます。

$$ \mathrm{Kurtosis} = \frac{\mathbb{E}[(x - \mu)^4]}{\sigma^4} - 3 = \frac{ 1 - 6 \phi (1 - \phi) }{ M \phi (1 - \phi) } $$


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

参考文献

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

おわりに

 記事を分割しました。

  • 2022.06.23:追記しました。


【次の内容】

www.anarchive-beta.com