はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
この記事では、R言語でベルヌーイ分布の計算を行います。
【前の内容】
【他の記事一覧】
【この記事の内容】
ベルヌーイ分布の計算
ベルヌーイ分布(Bernoulli Distribution)の確率と統計量を計算します。ベルヌーイ分布については「ベルヌーイ分布の定義式 - からっぽのしょこ」を参照してください。
確率の計算
ベルヌーイ分布に従う確率を計算する方法をいくつか確認します。
パラメータの設定
まずは、ベルヌーイ分布のパラメータ$\phi$を設定します。
# パラメータを指定 phi <- 0.35 # 確率変数の値を指定:(0, 1) x <- 1
ベルヌーイ分布のパラメータ$0 \leq \phi \leq 1$と確率変数がとり得る値$x \in \{0, 1\}$を指定します。設定した値に従う確率を計算します。
スクラッチで計算
定義式から確率を計算します。
# 定義式により確率を計算 prob <- phi^x * (1 - phi)^(1 - x) prob
## [1] 0.35
ベルヌーイ分布の定義式
で計算します。
対数をとった定義式から計算します。
# 対数をとった定義式により確率を計算 log_prob <- x * log(phi) + (1 - x) * log(1 - phi) prob <- exp(log_prob) prob; log_prob
## [1] 0.35 ## [1] -1.049822
対数をとった定義式
を計算します。計算結果の指数をとると確率が得られます。
指数と対数の性質より$\exp(\log x) = x$です。
次は、関数を使って確率を計算します。
二項分布の関数による計算
二項分布の確率関数dbinom()を使って計算します。
# 二項分布の関数により確率を計算 prob <- dbinom(x = x, size = 1, prob = phi) prob
## [1] 0.35
試行回数の引数sizeに1を指定することで、ベルヌーイ分布の確率を計算できます。成功回数の引数xにx、成功確率の引数probにphiを指定します。
log = TRUEを指定すると対数をとった確率を返します。
# 二項分布の対数をとった関数により確率を計算 log_prob <- dbinom(x = x, size = 1, prob = phi, log = TRUE) prob <- exp(log_prob) prob; log_prob
## [1] 0.35 ## [1] -1.049822
計算結果の指数をとると確率が得られます。
多項式分布の関数による計算
以降の計算は、変数とパラメータをベクトルに変換しておく必要があります。
# ベクトルに変換 phi_v <- c(1 - phi, phi) x_v <- c(1 - x, x) phi_v; x_v
## [1] 0.65 0.35 ## [1] 0 1
失敗確率($x = 0$となる確率)$1 - \phi$と成功確率($x = 1$となる確率)$\phi$をベクトルに格納してphi_vとします。これは、定義式にそれぞれ代入すると、0乗は1になるので($a^0 = 1$なので)、$x = 0$のとき$\mu^0 (1 - \mu)^{1-0} = 1 - \mu$、$x = 1$のとき$\mu^1 (1 - \mu)^{1-1} = \mu$になることに対応しています。
$x = 0$のとき$1 - x = 1$、$x = 1$のとき$1 - x = 0$となるのを利用して、one-hotベクトル(1-of-K符号化法)の$x$を作成してx_vとします。xが0のときx_vはc(1, 0)、xが1のときx_vはc(0, 1)になります。x_vの1のインデックスによってxの値を表します。
多項分布の確率関数dmultinom()を使って計算します。
# 多項分布の関数により確率を計算 prob <- dmultinom(x = x_v, size = 1, prob = phi_v) prob
## [1] 0.35
二項分布のときと同様に試行回数の引数をsize = 1として、出現頻度の引数xにx_v、出現確率の引数probにphi_vを指定します。
x_vの値が1の要素(インデックス)に対応したphi_vの要素(値)が出力されます。xが0のときx_vの1番目の要素が1なので1-phi、xが1のときx_vの2番目の要素が1なのでphiになります。
log = TRUEを指定すると対数をとった確率を返します。
# 多項分布の対数をとった関数により確率を計算 log_prob <- dmultinom(x = x_v, size = 1, prob = phi_v, log = TRUE) prob <- exp(log_prob) prob; log_prob
## [1] 0.35 ## [1] -1.049822
計算結果の指数をとると確率が得られます。
インデックスで抽出
最後に、添字による抽出で確率を取り出します。
# インデックスにより確率を抽出 prob <- phi_v[x+1] prob
## [1] 0.35
xが0となる確率はphi_vの1番目の要素で、1となる確率はphi_vの2番目の要素なので、xの値に対応するパラメータ(確率)のインデックスはx+1になります。
統計量の計算
ベルヌーイ分布の平均と分散を計算します。
平均を計算します。
# 平均を計算 E_x <- phi E_x
## [1] 0.35
ベルヌーイ分布の平均は、パラメータ$\phi$です。
分散を計算します。
# 分散を計算 V_x <- phi * (1 - phi) V_x
## [1] 0.2275
ベルヌーイ分布の分散は、次の式で計算できます。
$\phi$が0.5に近いほど分散が大きくなります。
最頻値を計算します。
# 最頻値を計算:(注:複数の場合も1つしか返さない) mode_x <- which.max(c(1 - phi, phi)) - 1 mode_x
## [1] 0
ベルヌーイ分布の最頻値は、次の条件で計算できます。
which.max()にベクトルを指定すると、最大値のインデックスを返します。そこで、1-phiとphiをベクトルに格納して、大きい方のインデックスから1を引くと、最頻値が得られます。ただし、$\phi = 0.5$のとき$x = 0$しか得られません。
歪度と尖度の計算
二項分布の歪度と尖度を計算します。詳しくはいつか書きます。
歪度を計算します。
# 歪度を計算 skewness <- (1 - 2 * phi) / sqrt(phi * (1 - phi)) skewness
## [1] 0.6289709
ベルヌーイ分布の歪度は、次の式で計算できます。
ここで、$x$の平均$\mu = \mathbb{E}[x]$、標準偏差$\sigma = \sqrt{\mathbb{E}[(x - \mu)^2]}$です。
尖度を計算します。
# 尖度を計算 kurtosis <- (1 - 6 * phi * (1 - phi)) / (phi * (1 - phi)) kurtosis
## [1] -1.604396
ベルヌーイ分布の尖度は、次の式で計算できます。
この記事では、ベルヌーイ分布の計算を確認しました。次は、グラフを作成します。
参考文献
- 岩田具治『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社,2015年.
おわりに
追記の際に「ベルヌーイ分布の作図」から記事を分割しました。
【次の内容】
