はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
この記事では、R言語でベルヌーイ分布の乱数を生成します。
【前の内容】
【他の記事一覧】
【この記事の内容】
ベルヌーイ分布の乱数生成
ベルヌーイ分布(Bernoulli Distribution)の乱数を生成します。ベルヌーイ分布については「ベルヌーイ分布の定義式 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用するパッケージ library(tidyverse) library(gganimate)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。また、基本的にベースパイプ(ネイティブパイプ)演算子|>を使いますが、パイプ演算子%>%ででないと処理できない部分があるため、magrittrも読み込む必要があります。
分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。
サンプリング
ベルヌーイ分布の乱数を生成します。
ベルヌーイ分布のパラメータ$\phi$とデータ数$N$を設定します。
# パラメータを指定 phi <- 0.35 # データ数を指定 N <- 1000
成功確率$0 < \phi < 1$とデータ数(サンプルサイズ)$N$を指定します。
ベルヌーイ分布に従う乱数を生成します。
# ベルヌーイ分布に従う乱数を生成 x_n <- rbinom(n = N, size = 1, prob = phi) head(x_n)
## [1] 0 0 1 1 0 0
二項分布の乱数生成関数rbinom()の試行回数の引数sizeに1を指定することで、ベルヌーイ分布に従う乱数を生成できます。データ数(サンプルサイズ)の引数nにN、成功確率の引数probにphiを指定します。
乱数の可視化
サンプルのヒストグラムを作成します。
サンプルの値を集計します。
# xがとり得る値を作成 x_vals <- 0:1 # サンプルを集計 freq_df <- tidyr::tibble( x = x_vals, # 確率変数 frequency = c(N - sum(x_n), sum(x_n)) # 度数 ) freq_df
## # A tibble: 2 × 2 ## x frequency ## <int> <int> ## 1 0 638 ## 2 1 362
ベルヌーイ分布の乱数(確率変数)は0か1の値をとるので、1の要素数はx_nの総和で得られます。また、データ数はNなので、0の要素数はNと1の要素数の差で得られます。
確率変数がとり得る値と対応する度数をデータフレームに格納します。
サンプルのヒストグラムを作成します。
# サンプルのヒストグラムを作成:度数 ggplot(data = freq_df, mapping = aes(x = x, y = frequency)) + # データ geom_bar(stat = "identity", fill = "#00A968") + # 度数 scale_x_continuous(breaks = x_vals, labels = x_vals) + # x軸目盛 labs(title = "Bernoulli Distribution", subtitle = paste0("phi=", phi, ", N=", N), x = "x", y = "frequency") # ラベル
度数を高さとする棒グラフを作成します。
相対度数を分布と重ねて描画します。
# ベルヌーイ分布を計算 prob_df <- tidyr::tibble( x = x_vals, # 確率変数 probability = c(1 - phi, phi) # 確率 ) # サンプルのヒストグラムを作成:相対度数 ggplot() + geom_bar(data = freq_df, mapping = aes(x = x, y = frequency/N), stat = "identity", fill = "#00A968") + # 構成比 geom_bar(data = prob_df, mapping = aes(x = x, y = probability), stat = "identity", alpha = 0, color = "darkgreen", linetype = "dashed") + # 元の分布 scale_x_continuous(breaks = x_vals, labels = x_vals) + # x軸目盛 ylim(c(0, 1)) + # y軸の表示範囲 labs( title = "Bernoulli Distribution", subtitle = parse( text = paste0( "list(phi==", phi, ", N==(list(", paste0(freq_df[["frequency"]], collapse = ", "), ")))") ), x = "x", y = "relative frequency" ) # ラベル
度数frequencyをデータ数Nで割った相対度数を高さとする棒グラフを作成します。
ギリシャ文字などの記号を使った数式を表示する場合は、expression()の記法を使います。等号は"=="、複数の(数式上の)変数を並べるには"list(変数1, 変数2)"とします。(プログラム上の)変数の値を使う場合は、parse()のtext引数に指定します。
データ数(サンプルサイズ)が十分に増えると、ヒストグラムの形が分布の形に近付きます。
乱数と分布の関係をアニメーションで可視化
サンプルサイズとヒストグラムの形状の関係をアニメーションで確認します。
・作図コード(クリックで展開)
パラメータ$\phi$とデータ数$N$を指定して、ベルヌーイ分布の乱数を生成します。
# パラメータを指定 phi <- 0.35 # # データ数(フレーム数)を指定 N <- 300 # ベルヌーイ分布に従う乱数を生成 x_n <- rbinom(n = N, size = 1, prob = phi) head(x_n)
## [1] 0 1 0 1 1 0
x_nのn番目の要素を、n回目にサンプリングされた値とみなします。アニメーションのn番目のフレームでは、n個のサンプルx_n[1:n]のヒストグラムを描画します。
サンプリング回数ごとに、$x$がとり得る値x_valsごとの度数を求めます。
# xがとり得る値を作成 x_vals <- 0:1 # サンプルを集計 freq_df <- tibble::tibble( x = x_n, # サンプル n = 1:N, # データ番号 frequency = 1 # 集計用の値 ) |> dplyr::right_join(tidyr::expand_grid(x = x_vals, n = 1:N), by = c("x", "n")) |> # 全てのパターンに結合 dplyr::mutate(frequency = tidyr::replace_na(frequency, replace = 0)) |> # サンプルにない場合の欠損値を0に置換 dplyr::arrange(n, x) |> # 集計用に昇順に並べ替え dplyr::group_by(x) |> # 集計用にグループ化 dplyr::mutate(frequency = cumsum(frequency)) |> # 累積和を計算 dplyr::ungroup() # グループ化を解除 freq_df
## # A tibble: 600 × 3 ## x n frequency ## <int> <int> <dbl> ## 1 0 1 1 ## 2 1 1 0 ## 3 0 2 1 ## 4 1 2 1 ## 5 0 3 2 ## 6 1 3 1 ## 7 0 4 2 ## 8 1 4 2 ## 9 0 5 2 ## 10 1 5 3 ## # … with 590 more rows
サンプルx_n、データ番号(サンプリング回数)1:N、度数計算用の値1をデータフレームに格納します。各データ番号において、実際にサンプリングされた値がデータフレームに含まれます。
$x$の値とデータ番号の全ての組み合わせをexpand_grid()で作成して、そこにサンプルをright_join()で結合します。サンプルにない場合は、frequency列が欠損値NAになるので、replace_na()で0に置換します。
x列($x$の値)でグループ化して、frequency列(サンプルは1でそれ以外は0)の累積和をcumsum()で計算して、各試行回数までの度数を求めます。
ヒストグラムの作図に利用する集計結果が得られました。以降は、グラフの装飾のための処理です。
$x$がとり得る値(x_valseの各要素)ごと度数を使って、フレーム切替用のラベルを作成します。
# フレーム切替用のラベルを作成 label_vec <- freq_df |> tidyr::pivot_wider( id_cols = n, names_from = x, names_prefix = "x", values_from = frequency ) |> # 度数列を展開 tidyr::unite(col = "label", dplyr::starts_with("x"), sep = ", ") |> # 度数情報をまとめて文字列化 dplyr::mutate( label = paste0("phi=", phi, ", N=", n, "=(", label, ")") %>% factor(., levels = .) ) |> # パラメータ情報をまとめて因子化 dplyr::pull(label) # ベクトルとして取得 head(label_vec)
## [1] phi=0.35, N=1=(1, 0) phi=0.35, N=2=(1, 1) phi=0.35, N=3=(2, 1) ## [4] phi=0.35, N=4=(2, 2) phi=0.35, N=5=(2, 3) phi=0.35, N=6=(3, 3) ## 300 Levels: phi=0.35, N=1=(1, 0) phi=0.35, N=2=(1, 1) ... phi=0.35, N=300=(199, 101)
pivot_wider()で、x列(x_vals)の値ごとの度数列に展開します。データ番号(サンプリング回数)の列nとx_valsに対応したM+1個の列になります。
度数列の値をunite()で文字列としてまとめて、さらにpaste0()でパラメータの値と結合し、因子型に変換します。
作成したラベル列をpull()でベクトルとして取り出します。
集計結果のデータフレームにフレーム切替用のラベル列を追加します。
# フレーム切替用のラベルを追加 anime_freq_df <- freq_df |> tibble::add_column(parameter = rep(label_vec, each = length(x_vals))) anime_freq_df
## # A tibble: 600 × 4 ## x n frequency parameter ## <int> <int> <dbl> <fct> ## 1 0 1 1 phi=0.35, N=1=(1, 0) ## 2 1 1 0 phi=0.35, N=1=(1, 0) ## 3 0 2 1 phi=0.35, N=2=(1, 1) ## 4 1 2 1 phi=0.35, N=2=(1, 1) ## 5 0 3 2 phi=0.35, N=3=(2, 1) ## 6 1 3 1 phi=0.35, N=3=(2, 1) ## 7 0 4 2 phi=0.35, N=4=(2, 2) ## 8 1 4 2 phi=0.35, N=4=(2, 2) ## 9 0 5 2 phi=0.35, N=5=(2, 3) ## 10 1 5 3 phi=0.35, N=5=(2, 3) ## # … with 590 more rows
label_vecの各要素をM+1個(x_valsの要素数)に複製して、add_column()で列として追加します。
(M+1)N行×のデータフレームになります。これをヒストグラムの描画に使います。
サンプルと対応するラベルをデータフレームに格納します。
# サンプルを格納 anime_data_df <- tibble::tibble( x = x_n, # サンプル parameter = label_vec # フレーム切替用ラベルを追加 ) anime_data_df
## # A tibble: 300 × 2 ## x parameter ## <int> <fct> ## 1 0 phi=0.35, N=1=(1, 0) ## 2 1 phi=0.35, N=2=(1, 1) ## 3 0 phi=0.35, N=3=(2, 1) ## 4 1 phi=0.35, N=4=(2, 2) ## 5 1 phi=0.35, N=5=(2, 3) ## 6 0 phi=0.35, N=6=(3, 3) ## 7 1 phi=0.35, N=7=(3, 4) ## 8 0 phi=0.35, N=8=(4, 4) ## 9 0 phi=0.35, N=9=(5, 4) ## 10 1 phi=0.35, N=10=(5, 5) ## # … with 290 more rows
N行のデータフレームになります。各サンプリングにおけるサンプルの値を描画するのに使います。
ベルヌーイ分布を計算して、データ数分に複製し、対応するラベルを追加します。
# ベルヌーイ分布の情報を複製 anime_prob_df <- tibble::tibble( x = x_vals, # 確率変数 probability = c(1 - phi, phi), # 確率 num = N # 複製数 ) |> tidyr::uncount(num) |> # データ数分に複製 tibble::add_column(parameter = rep(label_vec, times = length(x_vals))) |> # フレーム切替用ラベルを追加 dplyr::arrange(parameter) # サンプリング回数ごとに並べ替え anime_prob_df
## # A tibble: 600 × 3 ## x probability parameter ## <int> <dbl> <fct> ## 1 0 0.65 phi=0.35, N=1=(1, 0) ## 2 1 0.35 phi=0.35, N=1=(1, 0) ## 3 0 0.65 phi=0.35, N=2=(1, 1) ## 4 1 0.35 phi=0.35, N=2=(1, 1) ## 5 0 0.65 phi=0.35, N=3=(2, 1) ## 6 1 0.35 phi=0.35, N=3=(2, 1) ## 7 0 0.65 phi=0.35, N=4=(2, 2) ## 8 1 0.35 phi=0.35, N=4=(2, 2) ## 9 0 0.65 phi=0.35, N=5=(2, 3) ## 10 1 0.35 phi=0.35, N=5=(2, 3) ## # … with 590 more rows
確率分布の情報とデータ数Nをデータフレームに格納して、uncount()で複製します。
(M+1)×N行のデータフレームになります。これをベルヌーイ分布の棒グラフを描画するのに使います。
サンプルのヒストグラムのアニメーションを作成します。
# ベルヌーイ乱数のヒストグラムのアニメーションを作図:度数 anime_hist_graph <- ggplot() + geom_bar(data = anime_freq_df, mapping = aes(x = x, y = frequency), stat = "identity", fill = "#00A968") + # 度数 geom_point(data = anime_data_df, mapping = aes(x = x, y = 0), color = "orange", size= 5) + # サンプル gganimate::transition_manual(parameter) + # フレーム scale_x_continuous(breaks = x_vals, labels = x_vals) + # x軸目盛 labs(title = "Bernoulli Distribution", subtitle = "{current_frame}", x = "x", y = "frequency") # ラベル # gif画像を作成 gganimate::animate(anime_hist_graph, nframes = N, fps = 100, width = 600, height = 600)
transition_manual()にフレームの順序を表す列を指定します。
animate()のフレーム数の引数にnframesにデータ数(サンプルサイズ)、フレームレートの引数fpsに1秒当たりのフレーム数を指定します。fps引数の値が大きいほどフレームが早く切り替わります。
相対度数を分布と重ねたアニメーションを作成します。
# ベルヌーイ乱数のヒストグラムのアニメーションを作図:相対度数 anime_prop_graph <- ggplot() + geom_bar(data = anime_freq_df, mapping = aes(x = x, y = frequency/n), stat = "identity", fill = "#00A968") + # 相対度数 geom_bar(data = anime_prob_df, mapping = aes(x = x, y = probability), stat = "identity", alpha = 0, color = "darkgreen", linetype = "dashed") + # 元の分布 geom_point(data = anime_data_df, mapping = aes(x = x, y = 0), color = "orange", size= 5) + # サンプル gganimate::transition_manual(parameter) + # フレーム scale_x_continuous(breaks = x_vals, labels = x_vals) + # x軸目盛 ylim(c(0, 1)) + # y軸の表示範囲 labs(title = "Bernoulli Distribution", subtitle = "{current_frame}", x = "x", y = "relative frequency") # ラベル # gif画像を作成 gganimate::animate(anime_prop_graph, nframes = N, fps = 100, width = 600, height = 600)
サンプルが増えるに従って、ヒストグラムが元の分布に近付くのを確認できます。
この記事では、ベルヌーイ分布の乱数を確認しました。
参考文献
- 岩田具治『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社,2015年.
おわりに
加筆修正の際に「ベルヌーイ分布の作図」から記事を分割しました。
【次の内容】
