はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
二項分布の乱数の生成をR言語で行います。
【前の内容】
【他の記事一覧】
【この記事の内容】
Rで二項分布の乱数生成
二項分布(Binomial Distribution)の乱数を生成します。二項分布については「二項分布の定義式の確認 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用するパッケージ library(tidyverse) library(gganimate)
分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。
サンプリング
二項分布の乱数を生成します。
二項分布のパラメータ$\phi, M$とデータ数$N$を設定します。
# パラメータを指定 phi <- 0.3 # 試行回数を指定 M <- 10 # データ数を指定 N <- 1000
成功確率$0 < \phi < 1$、試行回数$M$とデータ数$N$を指定します。
二項分布に従う乱数を生成します。
# 二項分布に従う乱数を生成 x_n <- rbinom(n = N, size = M, prob = phi) x_n[1:5]
## [1] 0 2 3 3 3
二項分布の乱数はrbinom()で生成できます。データ数(サンプルサイズ)の引数nにN、試行回数の引数sizeにM、成功確率の引数probにphiを指定します。
乱数の可視化
サンプルのヒストグラムを作成します。
サンプルの値を集計します。
# 乱数を集計して格納 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> <int> <dbl> ## 1 0 27 0.027 ## 2 1 112 0.112 ## 3 2 241 0.241 ## 4 3 255 0.255 ## 5 4 205 0.205 ## 6 5 105 0.105
サンプリングした値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 = ", "), ")")) # ラベル
相対度数を分布と重ねて描画します。
# 作図用のxの点を作成 x_vals <- 0:M # 二項分布の情報を格納 prob_df <- tidyr::tibble( x = x_vals, # 確率変数 probability = dbinom(x = x_vals, size = M, prob = phi) # 確率 ) # サンプルの相対度数を作図 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 = ", "), ")")) # ラベル
二項分布のグラフの作成については「【R】二項分布の作図 - からっぽのしょこ」を参照してください。
データ数が十分に増えると、ヒストグラムの形が分布の形に近付きます。
推移の可視化
サンプルサイズとヒストグラムの変化をアニメーションで確認します。
1データずつアニメーションで可視化
1データずつサンプルを増やしてヒストグラムの変化を確認します。
・作図コード(クリックで展開)
# データ数(フレーム数)を指定 N_frame <- 100 # 乱数を1つずつ生成 #x_n <- rep(NA, times = N_frame) anime_freq_df <- tidyr::tibble() anime_data_df <- tidyr::tibble() anime_prob_df <- tidyr::tibble() for(n in 1:N_frame) { # 二項分布に従う乱数を生成 #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> <int> <dbl> <fct> ## 1 0 1 1 phi=0.3, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) ## 2 1 0 0 phi=0.3, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) ## 3 2 0 0 phi=0.3, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) ## 4 3 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)
for()ループを使って、「x_nから1データずつ取り出す」または「1データずつ生成」して、結果をanime_***_dfに追加していきます。
フレームごとに分布の情報prob_dfを複製し、対応するparameter列を追加する必要があります。
ヒストグラムのアニメーションを作成します。
# アニメーション用のサンプルのヒストグラムを作成 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)
サンプルが増えるに従って、ヒストグラムが元の分布に近付くのを確認できます。
全データをアニメーションで可視化
続いて、フレーム数を指定して、フレームごとに均等に分けたサンプルを増やしてヒストグラムの変化をアニメーションで確認します。
・作図コード(クリックで展開)
# フレーム数を指定 frame_num <- 100 # 1フレーム当たりのデータ数を計算 n_per_frame <- N %/% frame_num # 乱数を1つずつ生成 anime_freq_df <- tidyr::tibble() anime_prob_df <- tidyr::tibble() for(n in 1:frame_num) { # n個の乱数を集計して格納 tmp_freq_df <- tidyr::tibble(x = x_n[1:(n*n_per_frame)]) %>% # 乱数を格納 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*n_per_frame)) # 相対度数を計算 # ラベル用のテキストを作成 label_text <- paste0( "phi=", phi, ", N=", n*n_per_frame, "=(", paste0(tmp_freq_df[["frequency"]], collapse = ", "), ")" ) # フレーム切替用のラベルを付与 tmp_freq_df <- tmp_freq_df %>% dplyr::mutate(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_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 1 0.1 phi=0.3, N=10=(1, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0) ## 2 1 1 0.1 phi=0.3, N=10=(1, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0) ## 3 2 3 0.3 phi=0.3, N=10=(1, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0) ## 4 3 5 0.5 phi=0.3, N=10=(1, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0) ## 5 4 0 0 phi=0.3, N=10=(1, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0) ## 6 5 0 0 phi=0.3, N=10=(1, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0)
フレーム数frame_numを指定して、1フレーム当たりのデータ数n_per_frameを計算します。%/%は、整数商(割り算の整数部分)を計算します。
n*n_per_frame個ずつサンプルを取り出して、集計結果をanime_freq_dfに追加していきます。
ヒストグラムのアニメーションを作成します。
# アニメーション用のサンプルのヒストグラムを作成 anime_freq_graph <- ggplot() + geom_bar(data = anime_freq_df, mapping = aes(x = x, y = frequency), 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_freq_graph, nframes = frame_num, 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") + # 元の分布 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 = frame_num, fps = 100)
以上で、二項分布を確認できました。
参考文献
- 岩田具治『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社,2015年.
おわりに
記事を分割しました。
【次の内容】