からっぽのしょこ

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

【R】二項分布の作図

はじめに

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

 二項分布の計算とグラフの作成をR言語で行います。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【目次】

二項分布

 二項分布(Binomial Distribution)の計算と作図を行います。

 利用するパッケージを読み込みます。

# 利用するパッケージ
library(tidyverse)
library(gganimate)

 分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。

定義式の確認

 表(当たり)の出る確率が$\phi$であるコイン(くじ)を$M$回試行するとき、表(当たり)の出る回数$x$の確率分布を二項分布と言います。
 $x$は、コインを$M$回投げて表が出た回数なので、全て裏のときの0から全て表のときの$M$の整数になります。$x$が0から$M$の値をとることを$x = \{0, 1, \cdots, M\}$で表します。
 ベルヌーイ分布と同様に、成功確率(表・当たりとなる確率)を$0 \leq \phi \leq 1$、失敗確率(裏・外れとなる確率)を$1 - \phi$とします。

 二項分布は、パラメータ$\phi$を用いて次の式で定義されます。

$$ \mathrm{Bin}(x | M, \phi) = \frac{M!}{(M - x)! x!} \phi^x (1 - \phi)^{M-x} $$

 ここで、正規化項

$$ {}_M\mathrm{C}_x = \frac{M!}{(M - x)! x!} $$

は、試行回数が$M$回のとき表が$x$回となる組合せの数に対応します。

 例えば、コインを3回投げて($M = 3$で)表が2回出る($x = 2$となる)とき、二項分布の式は

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

表が2回出る確率$\phi^2$と裏が1回出る確率$1 - \phi$と正規化項の積となります。このとき正規化項は、3回中2回表となる組み合わせを表し、「表・表・裏」「表・裏・表」「裏・表・表」の3通りであることを求めています。
 このように、$x$の値に対応した確率となるように式が定義されています。

 また試行回数が1回($M = 1$)のとき、二項分布は

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

ベルヌーイ分布と等しくなります。階乗の定義より$0! = 1$なので、$x$が0か1のどちらであっても正規化項は1になります。
 ベルヌーイ分布は、コインを1回投げて表が1回出ることを$x = 1$、表が1回も出ない(裏が出る)ことを$x = 0$で表していると言えます。

 二項分布の対数をとると

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

となります。対数の性質より$\log x^a = a \log x$です。

 二項分布の平均と分散は、次の式で計算できます。詳しくは「1.2.1:二項分布【『トピックモデル』の勉強ノート】 - からっぽのしょこ」を参照してください。

$$ \begin{aligned} \mathbb{E}[x] &= M \phi \\ \mathbb{V}[x] &= M \phi (1 - \phi) \end{aligned} $$


 これらの計算を行いグラフを作成します。

確率の計算

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

 パラメータを設定します。

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

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

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

 二項分布のパラメータ$0 \leq \phi \leq 1$、試行回数$M$、確率変数がとり得る値$x \in \{0, 1, \cdots, 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!} \\ p(x | \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 p(x | \phi) &= \log C_{\mathrm{Bin}} + x \log \phi + (N - x) \log (1 - \phi) \end{aligned} $$

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

$$ p(x | \phi) = \exp \Bigr( \log p(x | \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 $$

 分散を計算します。

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

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

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


グラフの作成

 ggplot2パッケージを利用して二項分布のグラフを作成します。

 二項分布の確率変数$x$がとり得る値ごとの確率を計算します。

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

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

# 作図用のxの値を作成
x_vals <- 0:M

# 二項分布の情報を格納
prob_df <- tidyr::tibble(
  x = x_vals, # 確率変数
  probability = dbinom(x = x_vals, size = M, prob = phi) # 確率
)
head(prob_df)
## # A tibble: 6 x 2
##       x probability
##   <int>       <dbl>
## 1     0      0.0282
## 2     1      0.121 
## 3     2      0.233 
## 4     3      0.267 
## 5     4      0.200 
## 6     5      0.103

 $x$がとり得る値0:Mと、それに対応する確率をデータフレームに格納します。

 二項分布のグラフを作成します。

# 二項分布のグラフを作成
ggplot(data = prob_df, mapping = aes(x = x, y = probability)) + # データ
  geom_bar(stat = "identity", position = "dodge", fill = "#00A968") + # 分布
#  geom_vline(xintercept = E_x, color = "orange", size = 1, linetype = "dashed") + # 平均
#  geom_vline(xintercept = E_x - sqrt(V_x), color = "orange", size = 1, linetype = "dotted") + # 平均 - 標準偏差
#  geom_vline(xintercept = E_x + sqrt(V_x), color = "orange", size = 1, linetype = "dotted") + # 平均 + 標準偏差
  scale_x_continuous(breaks = 0:M, labels = 0:M) + # x軸目盛
  labs(title = "Binomial Distribution", 
       subtitle = paste0("phi=", phi, ", M=", M)) # ラベル

f:id:anemptyarchive:20220113043946p:plain
二項分布のグラフ

 二項分布のグラフを描画できました。

パラメータと分布の形状の関係

 続いて、パラメータ$\phi$の値を少しずつ変更して、分布の変化をアニメーションで確認します。

・作図コード(クリックで展開)

# 作図用のphiの値を作成
phi_vals <- seq(from = 0, to = 1, by = 0.01)

# phiの値ごとに分布を計算
anime_prob_df <- tidyr::tibble()
for(phi in phi_vals) {
  # 二項分布の情報を格納
  tmp_prob_df <- tidyr::tibble(
    x = x_vals, # 確率変数
    probability = dbinom(x = x_vals, size = M, prob = phi), # 確率
    parameter = paste0("phi=", phi, ", M=", M) %>% 
      as.factor() # フレーム切替用のラベル
  )
  
  # 結果を結合
  anime_prob_df <- rbind(anime_prob_df, tmp_prob_df)
}
head(anime_prob_df)
##       x probability parameter  
##   <int>       <dbl> <fct>      
## 1     0           1 phi=0, M=10
## 2     1           0 phi=0, M=10
## 3     2           0 phi=0, M=10
## 4     3           0 phi=0, M=10
## 5     4           0 phi=0, M=10
## 6     5           0 phi=0, M=10

 $\phi$がとり得る値を作成してphi_valsとします。
 for()ループを使ってphi_valsの値ごとに確率を計算して、anime_prob_dfに追加していきます。

 gif画像を作成します。

# アニメーション用の二項分布を作図
anime_prob_graph <- ggplot(data = anime_prob_df, mapping = aes(x = x, y = probability)) + # データ
  geom_bar(stat = "identity", position = "dodge", fill = "#00A968") + # 分布
  scale_x_continuous(breaks = 0:M, labels = 0:M) + # x軸目盛
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Binomial Distribution", 
       subtitle = "{current_frame}") # ラベル

# gif画像を作成
gganimate::animate(anime_prob_graph, nframes = length(phi_vals), fps = 100)


f:id:anemptyarchive:20220113044037g:plain
二項分布とパラメータの関係

 パラメータ$\phi$の値が大きくなるのに従って、成功回数$x$が大きいほど確率が高くなる(山が右に移動する)のを確認できます。

 続いて、$\phi$を固定して$M$を変更してみます。

・作図コード(クリックで展開)

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

# 試行回数の最大値を指定
M_max <- 100

# Mの値ごとに確率を計算
anime_prob_df <- tidyr::tibble()
for(M in 1:M_max) {
  # 作図用のxの値を作成
  x_vals <- 0:M
  
  # 二項分布の情報を格納
  tmp_prob_df <- tidyr::tibble(
    x = x_vals, # 表の回数
    prob = dbinom(x = x_vals, size = M, prob = phi), # 確率
    parameter = paste0("phi=", phi, ", M=", M) %>% 
      as.factor() # フレーム切替用ラベル
  )
  
  # 計算結果を結合
  anime_prob_df <- rbind(anime_prob_df, tmp_prob_df)
}
head(anime_prob_df)
## # A tibble: 6 x 3
##       x   prob parameter   
##   <int>  <dbl> <fct>       
## 1     0 0.7    phi=0.3, M=1
## 2     1 0.3    phi=0.3, M=1
## 3     0 0.490  phi=0.3, M=2
## 4     1 0.42   phi=0.3, M=2
## 5     2 0.0900 phi=0.3, M=2
## 6     0 0.343  phi=0.3, M=3

 for()ループを使って1:M_maxの値ごとに確率を計算して、anime_prob_dfに追加していきます。

 gif画像を作成します。

# アニメーション用の二項分布のグラフを作成
anime_prob_graph <- ggplot(data = anime_prob_df, mapping = aes(x = x, y = prob)) + # データ
  geom_bar(stat = "identity", position = "dodge", fill = "#00A968", color = "#00A968") + # 棒グラフ
#  scale_x_continuous(breaks = 0:M_max, labels = 0:M_max) + # x軸目盛
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Binomial Distribution", 
       subtitle = "{current_frame}") # ラベル

# gif画像を作成
gganimate::animate(anime_prob_graph, nframes = M_max, fps = 100)


f:id:anemptyarchive:20220113044134g:plain
二項分布と試行回数の関係

 試行回数$M$が増えるのに従って、成功回数$x$が大きいほど確率が高くなる(山が右に移動する)のを確認できます。ただし、$x$がとり得る値の範囲x_valsも広がっていくため、x_vals全体における山の位置は変わりません。

乱数の生成

 二項分布の乱数を生成してヒストグラムを確認します。

 パラメータを指定して、二項分布に従う乱数を生成します。

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

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

# データ数を指定
N <- 1000

# 二項分布に従う乱数を生成
x_n <- rbinom(n = N, size = M, prob = phi)
x_n[1:5]
## [1] 4 4 2 4 2

 二項分布の乱数生成関数rbinom()のデータ数(サンプルサイズ)の引数nN、試行回数の引数size1、成功確率の引数probphiを指定します。

 サンプルの値を集計します。

# 乱数を集計して格納
freq_df <- tidyr::tibble(x = x_n) %>% # 乱数を格納
  dplyr::count(x, name = "frequency") %>% # 頻度を測定
  dplyr::right_join(tidyr::tibble(x = 0:M), by = "x") %>% # 確率変数列を追加
  dplyr::mutate(frequency = tidyr::replace_na(frequency, 0)) %>% # サンプルにない場合のNAを0に置換
  dplyr::mutate(proportion = frequency / N) # 構成比を計算
head(freq_df)
## # A tibble: 6 x 3
##       x frequency proportion
##   <int>     <dbl>      <dbl>
## 1     0        34      0.034
## 2     1       122      0.122
## 3     2       247      0.247
## 4     3       284      0.284
## 5     4       176      0.176
## 6     5        97      0.097

 サンプリングした値x_nをデータフレームに格納します。
 count()で重複する値の数をカウントします。
 サンプルx_nにない値はデータフレームに保存されないため、0からMの値を持つデータフレームを作成し、right_join()で結合します。
 サンプルにない場合は頻度列frequencyが欠損値NAになるので、replace_na()0に置換します。
 得られた頻度frequencyをデータ数Nで割り、各値の構成比を計算してproportion列とします。

 ヒストグラムを作成します。

# サンプルのヒストグラムを作成
ggplot(data = freq_df, mapping = aes(x = x, y = frequency)) + # データ
  geom_bar(stat = "identity", position = "dodge", fill = "#00A968") + # ヒストグラム
  scale_x_continuous(breaks = 0:M, labels = 0:M) + # x軸目盛
  labs(title = "Binomial Distribution", 
       subtitle = paste0("phi=", phi, ", M=", M, 
                         ", N=", N, "=(", paste0(freq_df[["frequency"]], collapse = ", "), ")")) # ラベル

f:id:anemptyarchive:20220113044225p:plain
二項分布の乱数のヒストグラム


 構成比を分布と重ねて描画します。

# サンプルの構成比を作図
ggplot() + 
  geom_bar(data = freq_df, mapping = aes(x = x, y = proportion), 
           stat = "identity", position = "dodge", fill = "#00A968") + # 構成比
  geom_bar(data = prob_df, mapping = aes(x = x, y = probability), 
           stat = "identity", position = "dodge", alpha = 0, color = "darkgreen", linetype = "dashed") + # 真の分布
  scale_x_continuous(breaks = 0:M, labels = 0:M) + # x軸目盛
  labs(title = "Binomial Distribution", 
       subtitle = paste0("phi=", phi, ", M=", M, 
                         ", N=", N, "=(", paste0(freq_df[["frequency"]], collapse = ", "), ")")) # ラベル

f:id:anemptyarchive:20220113044245p:plain
二項分布の乱数の構成比

 データ数が十分に増えると分布のグラフに形が近づきます。

 サンプルサイズとヒストグラムの変化をアニメーションで確認します。

・作図コード(クリックで展開)

# データ数を指定
N <- 100

# 乱数を1つずつ生成
x_n <- rep(NA, times = N)
anime_freq_df <- tidyr::tibble()
anime_data_df <- tidyr::tibble()
anime_prob_df <- tidyr::tibble()
for(n in 1:N) {
  # 二項分布に従う乱数を生成
  x_n[n] <- rbinom(n = 1, size = M, prob = phi)
  
  # n個の乱数を集計して格納
  tmp_freq_df <- tidyr::tibble(x = x_n[1:n]) %>% # 乱数を格納
    dplyr::count(x, name = "frequency") %>% # 頻度を測定
    dplyr::right_join(tidyr::tibble(x = 0:M), by = "x") %>% # 確率変数列を追加
    dplyr::mutate(frequency = tidyr::replace_na(frequency, 0)) %>% # サンプルにない場合のNAを0に置換
    dplyr::mutate(proportion = frequency / n) # 構成比を計算
  
  # ラベル用のテキストを作成
  label_text <- paste0(
    "phi=", phi, ", N=", n, "=(", paste0(tmp_freq_df[["frequency"]], collapse = ", "), ")"
  )
  
  # フレーム切替用のラベルを付与
  tmp_freq_df <- tmp_freq_df %>% 
    dplyr::mutate(parameter = as.factor(label_text))
  
  # n番目の乱数を格納
  tmp_data_df <- tidyr::tibble(
    x = x_n[n], 
    parameter = as.factor(label_text) # フレーム切替用のラベル
  )
  
  # n回目のラベルを付与
  tmp_prob_df <- prob_df %>% 
    dplyr::mutate(parameter = as.factor(label_text))
  
  # 結果を結合
  anime_freq_df <- rbind(anime_freq_df, tmp_freq_df)
  anime_data_df <- rbind(anime_data_df, tmp_data_df)
  anime_prob_df <- rbind(anime_prob_df, tmp_prob_df)
}
head(anime_freq_df)
## # A tibble: 6 x 4
##       x frequency proportion parameter                                     
##   <int>     <dbl>      <dbl> <fct>                                         
## 1     3         1          1 phi=0.3, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
## 2     0         0          0 phi=0.3, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
## 3     1         0          0 phi=0.3, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
## 4     2         0          0 phi=0.3, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
## 5     4         0          0 phi=0.3, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
## 6     5         0          0 phi=0.3, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

 乱数を1つず生成して、結果をanime_***_dfに追加していきます。

 ヒストグラムのアニメーションを作成します。

# アニメーション用のサンプルのヒストグラムを作成
anime_freq_graph <- ggplot() + 
  geom_bar(data = anime_freq_df, mapping = aes(x = x, y = frequency), 
           stat = "identity", position = "dodge", fill = "#00A968") + # ヒストグラム
  geom_point(data = anime_data_df, mapping = aes(x = x, y = 0), 
             color = "orange", size= 5) + # サンプル
  scale_x_continuous(breaks = 0:M, labels = 0:M) + # x軸目盛
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Binomial Distribution", 
       subtitle = "{current_frame}") # ラベル

# gif画像を作成
gganimate::animate(anime_freq_graph, nframes = N, fps = 100)


 構成比のアニメーションを作成します。

# アニメーション用のサンプルの構成比を作図
anime_prop_graph <- ggplot() + 
  geom_bar(data = anime_freq_df, mapping = aes(x = x, y = proportion), 
           stat = "identity", position = "dodge", fill = "#00A968") + # 構成比
  geom_bar(data = anime_prob_df, mapping = aes(x = x, y = probability), 
           stat = "identity", position = "dodge", alpha = 0, color = "darkgreen", linetype = "dashed") + # 真の分布
  geom_point(data = anime_data_df, mapping = aes(x = x, y = 0), 
             color = "orange", size= 5) + # サンプル
  scale_x_continuous(breaks = 0:M, labels = 0:M) + # x軸目盛
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Binomial Distribution", 
       subtitle = "{current_frame}") # ラベル

# gif画像を作成
gganimate::animate(anime_prop_graph, nframes = N, fps = 100)


f:id:anemptyarchive:20220113044317g:plainf:id:anemptyarchive:20220113044320g:plain
二項分布の乱数の推移

 サンプルが増えるに従って、真の分布に近付いていくのを確認できます。

 以上で、二項分布を確認できました。

関連する記事

 Python版です。

www.anarchive-beta.com

 平均と分散を導出します。

www.anarchive-beta.com


おわりに

 加筆修正の際に青トピシリーズから独立させました。

 もっとサクサク進むと思ってましたが、意外と大変でした。

【次の内容】

www.anarchive-beta.com