からっぽのしょこ

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

【R】二項分布の計算

はじめに

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

 二項分布の計算をR言語で行います。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

Rで二項分布の計算

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

確率の計算

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

パラメータの設定

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

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

# 試行回数を指定
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.2668279

 二項分布の定義式

$$ \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()で計算できます。

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

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

 対数をとった定義式

$$ \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.2668279

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

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

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

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

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

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

# ベクトルに変換
phi_v <- c(1 - phi, phi)
x_v <- c(M - x, x)
phi_v; x_v
## [1] 0.7 0.3
## [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.2668279

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

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

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

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

統計量の計算

 二項分布の平均と分散を計算します。二項分布の統計量については「二項分布の平均と分散の導出 - からっぽのしょこ」を参照してください。

 平均を計算します。

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

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

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

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

 分散を計算します。

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

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

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

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

 以上で、二項分布の確率と統計量の計算を確認できました。次は、グラフを作成します。

参考文献

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

おわりに

 記事を分割しました。

【次の内容】

www.anarchive-beta.com