からっぽのしょこ

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

【R】ベルヌーイ分布の作図

はじめに

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

 この記事では、ベルヌーイ分布の計算とグラフの作成をR言語で行います。

【他の記事一覧】

www.anarchive-beta.com

【目次】

ベルヌーイ分布

 ベルヌーイ分布(Bernoulli Distribution)の計算と作図を行います。

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

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

 分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。

定義式の確認

 コインの裏表やくじの当たり外れのように、2値をとる変数の確率分布をベルヌーイ分布と言います。
 例えば、コインを投げて表なら1、裏なら0である変数を$x$とします。$x$が0か1の値をとることを$x \in {0,1}$で表します。ちなみに、コインを1回投げて表が1回出ることを$x = 1$、表が1回も出ない(裏が出る)ことを$x = 0$で表していると理解しておくと、次の二項分布との関係が分かりやすくなります。
 $x = 1$となる(表が出る)確率を$\phi$を使って表すことにします。$0 \leq \phi \leq 1$であり、$x = 0$となる(裏が出る)確率は$1 - \phi$になります。

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

$$ \mathrm{Bern}(x | \phi) = \phi^x (1 - \phi)^{1 - x} $$


 この式は、コインが表つまり$x = 1$のとき

$$ \begin{aligned} \mathrm{Bern}(x = 1 | \phi) &= \phi^1 (1 - \phi)^{1-1} \\ &= \phi * 1 \\ &= \phi \end{aligned} $$

となり、またコインが裏つまり$x = 0$のときは

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

となります。$x^0 = 1$です。
 このように、$x$の値に対応した確率となるように式が定義されています。

 ベルヌーイ分布の対数をとると

$$ \log \mathrm{Bern}(x | \phi) = x \log \phi + (1 - x) \log (1 - \phi) $$

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

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

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


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

確率の計算

 ベルヌーイ分布に従う確率を計算する方法をいくつか確認します。

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

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

# 確率変数の値を指定
x <- 0

 ベルヌーイ分布のパラメータ$0 \leq \phi \leq 1$と確率変数がとり得る値$x \in {0, 1}$を指定します。設定した値に従う確率を計算します。

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

# 定義式により確率を計算
prob <- phi^x * (1 - phi)^(1 - x)
prob
## [1] 0.7

 ベルヌーイ分布の定義式

$$ p(x | \phi) = \phi^x (1 - \phi)^{1 - x} $$

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

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

 対数をとった定義式

$$ \log p(x | \phi) = x \log \phi + (1 - x) \log (1 - \phi) $$

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

$$ p(x | \phi) = \exp \Bigr( \log p(x | \phi) \Bigr) $$

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

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

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

 試行回数の引数size1を指定することで、ベルヌーイ分布の確率を計算できます。
 成功回数の引数xx、成功確率の引数probphiを指定します。

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

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

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

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

# ベクトルに変換
phi_v <- c(1 - phi, phi)
x_v <- c(1 - x, x)
phi_v; x_v
## [1] 0.7 0.3
## [1] 1 0

 $x = 0$となる確率$1 - \phi$と、$x = 1$となる確率$\phi$を持つベクトルを作成してphi_vとします。
 $x = 0$のとき$1 - x = 1$、$x = 1$のとき$1 - x = 0$となるのを利用して、one-hotベクトル(1-of-K符号化法)の$x$を作成してx_vとします。

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

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

 二項分布と同様に、試行回数の引数をsize = 1として、出現頻度の引数xx、出現確率の引数probphi_vを指定します。

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

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

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

 最後に、スライス機能を使って確率を取り出します。

# インデックスにより確率を抽出
prob <- phi_v[x+1]
prob
## [1] 0.7

 x0となる確率はphi_vの1番目の要素で、1となる確率はphi_vの2番目の要素なので、xの値に対応するパラメータのインデックスはx + 1になります。

統計量の計算

 ベルヌーイ分布の平均と分散を計算します。

 平均を計算します。

# 平均を計算
E_x <- phi
E_x
## [1] 0.3

 ベルヌーイ分布の平均は、次の式で計算できます。

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

 分散を計算します。

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

 ベルヌーイ分布の分散は、次の式で計算できます。

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


グラフの作成

 ggplot2パッケージを利用してベルヌーイ分布のグラフを作成します。

 ベルヌーイ分布の確率変数$x$がとり得る値ごとの確率を計算します。

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

# ベルヌーイ分布の情報を格納
prob_df <- tidyr::tibble(
  x = 0:1, # 確率変数
  probability = c(1 - phi, phi) # 確率
)
prob_df
## # A tibble: 2 x 2
##       x probability
##   <int>       <dbl>
## 1     0         0.7
## 2     1         0.3

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

 ベルヌーイ分布のグラフを作成します。

# ベルヌーイ分布を作図
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:1, labels = 0:1) + # x軸目盛
  labs(title = "Bernoulli Distribution", 
       subtitle = paste0("phi=", phi)) # ラベル

f:id:anemptyarchive:20220111050841p: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 = 0:1, # 確率変数
    probability = c(1 - phi, phi), # 確率
    parameter = paste0("phi=", phi) %>% 
      as.factor() # フレーム切替用のラベル
  )
  
  # 結果を結合
  anime_prob_df <- rbind(anime_prob_df, tmp_prob_df)
}
head(anime_prob_df)
## # A tibble: 6 x 3
##       x probability parameter
##   <int>       <dbl> <fct>    
## 1     0        1    phi=0    
## 2     1        0    phi=0    
## 3     0        0.99 phi=0.01 
## 4     1        0.01 phi=0.01 
## 5     0        0.98 phi=0.02 
## 6     1        0.02 phi=0.02

 $\phi$がとり得る値を作成してphi_valsとします。
 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") + # 分布
  gganimate::transition_manual(parameter) + # フレーム
  scale_x_continuous(breaks = 0:1, labels = 0:1) + # x軸目盛
  labs(title = "Bernoulli Distribution", 
       subtitle = "{current_frame}") # ラベル

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


f:id:anemptyarchive:20220111050901g:plain
ベルヌーイ分布のグラフ

 (定義のままですが、)パラメータ$\phi$の値が大きくなるほど、$x = 0$となる確率が下がり、$x = 1$となる確率が上がるのを確認できます。

乱数の生成

 ベルヌーイ分布の乱数を生成してヒストグラムを確認します。

 パラメータを指定して、ベルヌーイ分布に従う乱数を生成します。

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

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

# ベルヌーイ分布に従う乱数を生成
x_n <- rbinom(n = N, size = 1, prob = phi)
x_n[1:5]
## [1] 1 1 0 0 0

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

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

# 乱数を集計して格納
freq_df <- tidyr::tibble(
  x = 0:1, # 確率変数
  frequency = c(sum(x_n == 0), sum(x_n == 1)), # 度数
  proportion = frequency / N # 構成比
)
freq_df
## # A tibble: 2 x 3
##       x frequency proportion
##   <int>     <int>      <dbl>
## 1     0       704      0.704
## 2     1       296      0.296

 x_nに含まれる0の要素数は、sum(x_n == 0)で得られます。値が1の要素数も同様に求めて、データフレームに格納します。
 また、得られた頻度frequencyをデータ数Nで割り、01の構成比を計算します。

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

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

f:id:anemptyarchive:20220111051000p: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:1, labels = 0:1) + # x軸目盛
  labs(title = "Bernoulli Distribution", 
       subtitle = paste0("phi=", phi, ", N=", N, "=(", sum(x_n == 0), ", ", sum(x_n == 1), ")")) # ラベル

f:id:anemptyarchive:20220111051012p: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 = 1, prob = phi)
  
  # ラベル用のテキストを作成
  label_text <- paste0(
    "phi=", phi, ", N=", n, "=(", sum(x_n[1:n] == 0), ", ", sum(x_n[1:n] == 1), ")"
  )
  
  # n個の乱数を集計して格納
  tmp_freq_df <- tidyr::tibble(
    x = 0:1, # 確率変数
    frequency = c(sum(x_n[1:n] == 0), sum(x_n[1:n] == 1)), # 度数
    proportion = frequency / n, # 構成比
    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 %>% 
    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>     <int>      <dbl> <fct>              
## 1     0         0      0     phi=0.3, N=1=(0, 1)
## 2     1         1      1     phi=0.3, N=1=(0, 1)
## 3     0         0      0     phi=0.3, N=2=(0, 2)
## 4     1         2      1     phi=0.3, N=2=(0, 2)
## 5     0         1      0.333 phi=0.3, N=3=(1, 2)
## 6     1         2      0.667 phi=0.3, N=3=(1, 2)

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

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

# アニメーション用のサンプルのヒストグラムを作成
anime_hist_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:1, labels = 0:1) + # x軸目盛
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Bernoulli Distribution", 
       subtitle = "{current_frame}") # ラベル

# gif画像を作成
gganimate::animate(anime_hist_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:1, labels = 0:1) + # x軸目盛
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Bernoulli Distribution", 
       subtitle = "{current_frame}") # ラベル

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


f:id:anemptyarchive:20220111125937g:plainf:id:anemptyarchive:20220111130004g:plain
乱数の推移

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

 以上で、ベルヌーイ分布を確認できました。

関連する記事

 ベルヌーイ分布を扱う機械学習関連の記事をいくつか貼っておきます。

 Python版です。

www.anarchive-beta.com

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

www.anarchive-beta.com

 ベイズ推論を導出します。

www.anarchive-beta.com

 ベイズ推論を実装します。

www.anarchive-beta.com

www.anarchive-beta.com

 その他、混合ベルヌーイ分布やロジスティック回帰も一応ベルヌーイ分布ネタかな?がこちらにあります。

www.anarchive-beta.com


おわりに

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

 元記事的にも現記事的にも確率分布ネタシリーズ1つ目の記事です。これの元を書いたのは2年半前です。懐かしい、あの頃の自分に色々教えてあげたい。

【次の内容】

www.anarchive-beta.com