はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
この記事では、R言語でベルヌーイ分布のグラフを作成します。
【前の内容】
【他の記事一覧】
【目次】
ベルヌーイ分布の作図
ベルヌーイ分布(Bernoulli Distribution)のグラフを作成します。ベルヌーイ分布については「ベルヌーイ分布の定義式 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用するパッケージ library(tidyverse) library(gganimate)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。
magrittrパッケージのパイプ演算子%>%ではなく、ベースパイプ(ネイティブパイプ)演算子|>を使っています。%>%に置き換えても処理できます。
分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。
定義式の確認
まずは、ベルヌーイ分布の定義式を確認します。
ベルヌーイ分布は、次の式で定義されます。
ここで、$x$は成功・失敗を表す値、$\phi$は成功確率($x = 1$となる確率)です。
確率変数の値$x$は、$x \in \{0, 1\}$となり、$x = 1$が成功・$x = 0$が失敗を表します。パラメータ$\phi$は、$\phi \in (0, 1)$を満たす必要があります。また、失敗確率($x = 0$となる確率)は$1 - \phi$で表せます。
ベルヌーイ分布の平均と分散は、次の式で計算できます。詳しくは「ベルヌーイ分布の平均と分散の導出 - からっぽのしょこ」を参照してください。
最頻値は、次の条件で計算できます。
ベルヌーイ分布の歪度と尖度は、次の式で計算できます。
これらの計算を行いグラフを作成します。
グラフの作成
ggplot2パッケージを利用して、ベルヌーイ分布のグラフを作成します。ベルヌーイ分布の確率や統計量の計算については「【R】ベルヌーイ分布の計算 - からっぽのしょこ」を参照してください。
ベルヌーイ分布のパラメータ$\phi$を設定します。
# パラメータを指定 phi <- 0.35
成功確率$0 < \phi < 1$を指定します。
ベルヌーイ分布の確率変数がとり得る値を作成します。
# xがとり得る値を作成 x_vals <- 0:1
$x$がとり得る値$x \in \{0, 1\}$を作成して、x_valsとします。
$x$の値ごとに確率を計算します。
# ベルヌーイ分布を計算 prob_df <- tidyr::tibble( x = x_vals, # 確率変数 probability = c(1 - phi, phi) # 確率 ) prob_df
## # A tibble: 2 × 2 ## x probability ## <int> <dbl> ## 1 0 0.65 ## 2 1 0.35
$x$がとり得る値0, 1と、それぞれに対応する確率1 - phi, phiをデータフレームに格納します。
ベルヌーイ分布のグラフを作成します。
# ベルヌーイ分布を作図 ggplot(data = prob_df, mapping = aes(x = x, y = probability)) + # データ geom_bar(stat = "identity", fill = "#00A968") + # 棒グラフ scale_x_continuous(breaks = x_vals, labels = x_vals) + # x軸目盛 ylim(c(0, 1)) + # y軸の表示範囲 labs(title = "Bernoulli Distribution", subtitle = paste0("phi=", phi), # (文字列表記用) #subtitle = parse(text = paste0("phi==", phi)), # (数式表記用) x = "x", y = "probability") # ラベル
ギリシャ文字などの記号を使った数式を表示する場合は、expression()の記法を使います。等号は"=="、複数の(数式上の)変数を並べるには"list(変数1, 変数2)"とします。(プログラム上の)変数の値を使う場合は、parse()のtext引数に指定します。
この分布に統計量の情報を重ねて表示します。
# 補助線用の統計量を計算 E_x <- phi s_x <- sqrt(phi * (1 - phi)) mode_x <- which.max(c(1 - phi, phi)) - 1 # ベルヌーイ分布を作図:線のみ ggplot(data = prob_df, mapping = aes(x = x, y = probability)) + # データ geom_bar(stat = "identity", fill = "#00A968") + # 分布 geom_vline(xintercept = E_x, color = "blue", size = 1, linetype = "dashed") + # 期待値 geom_vline(xintercept = E_x-s_x, color = "orange", size = 1, linetype = "dotted") + # 期待値 - 標準偏差 geom_vline(xintercept = E_x+s_x, color = "orange", size = 1, linetype = "dotted") + # 期待値 + 標準偏差 geom_vline(xintercept = mode_x, color = "chocolate", size = 1, 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("phi==", phi)), x = "x", y = "probability") # ラベル
期待値・標準偏差(分散の平方根)・最頻値を計算して、それぞれgeom_vline()で垂直線を引きます。標準偏差については期待値から足し引きした値を使います。
期待値(青色の破線)に近い整数$x$となる確率が最大であり、最頻値(茶色の破線)と一致するのを確認できます。また、各線の位置関係から分布の非対称性が分かります。
凡例を表示する場合は、垂線を引く値をデータフレームにまとめます。
# 統計量を格納 stat_df <- tibble::tibble( statistic = c(E_x, E_x-s_x, E_x+s_x, mode_x), # 統計量 type = c("mean", "sd", "sd", "mode") # 色分け用ラベル ) stat_df
## # A tibble: 4 × 2 ## statistic type ## <dbl> <chr> ## 1 0.35 mean ## 2 -0.127 sd ## 3 0.827 sd ## 4 0 mode
各統計量をそれぞれを区別するための文字列と共にデータフレームに格納します。期待値と標準偏差の差・和は同じタイプとします。
折れ線の色や種類、凡例テキストを指定する場合は、設定用の名前付きベクトルを作成します。
# 凡例用の設定を作成:(数式表示用) color_vec <- c(mean = "blue", sd = "orange", mode = "chocolate") linetype_vec <- c(mean = "dashed", sd = "dotted", mode = "dashed") label_vec <- c(mean = expression(E(x)), sd = expression(E(x) %+-% sqrt(V(x))), mode = expression(mode(x))) color_vec; linetype_vec; label_vec
## mean sd mode ## "blue" "orange" "chocolate" ## mean sd mode ## "dashed" "dotted" "dashed" ## expression(mean = E(x), sd = E(x) %+-% sqrt(V(x)), mode = mode(x))
各要素(設定)の名前をtype列の文字列と対応させて、設定を指定します。
凡例付きのベルヌーイ分布のグラフを作成します。
# 統計量を重ねたベルヌーイ分布のグラフを作成:凡例付き ggplot() + # データ geom_bar(data = prob_df, mapping = aes(x = x, y = probability), stat = "identity", fill = "#00A968") + # 分布 geom_vline(data = stat_df, mapping = aes(xintercept = statistic, color = type, linetype = type), size = 1) + # 統計量 scale_x_continuous(breaks = x_vals, labels = x_vals) + # x軸目盛 scale_color_manual(values = color_vec, labels = label_vec, name = "statistic") + # 線の色:(色指定と数式表示用) scale_linetype_manual(values = linetype_vec, labels = label_vec, name = "statistic") + # 線の種類:(線指定と数式表示用) theme(legend.text.align = 0) + # 図の体裁:凡例 ylim(c(0, 1)) + # y軸の表示範囲 labs(title = "Bernoulii Distribution", subtitle = parse(text = paste0("phi==", phi)), x = "x", y = "probability") # ラベル
scale_color_manual()のvalues引数に線の色、scale_linetype_manual()のvalues引数に線の種類を指定します。凡例テキストを指定する場合は、それぞれのlabels引数に指定します。names引数は凡例ラベルです。
凡例テキストは、theme()のlegend.text.align引数に0を指定すると左寄せ、デフォルト(1)だと右寄せになります。
ベルヌーイ分布のグラフを描画できました。
パラメータと分布の関係をアニメーションで可視化
次は、パラメータの値を少しずつ変化させて、分布の形状の変化をアニメーションで確認します。
・作図コード(クリックで展開)
パラメータ$\phi$の値を作成して、それぞれ分布を計算します。
# パラメータとして利用する値を指定 phi_vals <- seq(from = 0, to = 1, by = 0.01) length(phi_vals) # フレーム数 # xがとり得る値を作成 x_vals <- 0:1 # パラメータごとにベルヌーイ分布を計算 anime_prob_df <- tidyr::expand_grid( x = x_vals, phi = phi_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(phi, x) |> # パラメータごとに並べ替え dplyr::mutate( probability = dbinom(x = x, size = 1, prob = phi), parameter = paste0("phi=", phi) |> factor(levels = paste0("phi=", phi_vals)) # フレーム切替用ラベル ) # 確率を計算 anime_prob_df
## # A tibble: 202 × 4 ## x phi probability parameter ## <int> <dbl> <dbl> <fct> ## 1 0 0 1 phi=0 ## 2 1 0 0 phi=0 ## 3 0 0.01 0.99 phi=0.01 ## 4 1 0.01 0.01 phi=0.01 ## 5 0 0.02 0.98 phi=0.02 ## 6 1 0.02 0.02 phi=0.02 ## 7 0 0.03 0.97 phi=0.03 ## 8 1 0.03 0.03 phi=0.03 ## 9 0 0.04 0.96 phi=0.04 ## 10 1 0.04 0.04 phi=0.04 ## # … with 192 more rows
値の間隔が一定になるようにphi_valsを作成します。パラメータごとにフレームを切り替えるので、phi_valsの要素数がアニメーションのフレーム数になります。
$x$の値x_valsとパラメータphi_valsの各要素の全ての組み合わせをexpand_grid()で作成します。
組み合わせごとに確率を計算して、パラメータごとにラベルを作成します。ラベルは、フレームの制御に使います。文字列型だと文字列の基準で順序が決まるので、因子型にしてパラメータに応じたレベル(順序)を設定します。
ベルヌーイ分布のアニメーション(gif画像)を作成します。
# ベルヌーイ分布のアニメーションを作図 anime_prob_graph <- ggplot(data = anime_prob_df, mapping = aes(x = x, y = probability)) + # データ geom_bar(stat = "identity", fill = "#00A968") + # 分布 gganimate::transition_manual(parameter) + # フレーム scale_x_continuous(breaks = x_vals, labels = x_vals) + # x軸目盛 labs(title = "Bernoulli Distribution", subtitle = "{current_frame}") # ラベル # gif画像を作成 gganimate::animate(anime_prob_graph, nframes = length(phi_vals), fps = 100, width = 600, height = 600)
transition_manual()にフレームの順序を表す列を指定します。
animate()のフレーム数の引数にnframesにパラメータ数、フレームレートの引数fpsに1秒当たりのフレーム数を指定します。fps引数の値が大きいほどフレームが早く切り替わります。
(定義のままですが、)パラメータ$\phi$の値が大きくなるほど、失敗確率($x = 0$となる確率)が下がり、成功確率($x = 1$となる確率)が上がるのを確認できます。
パラメータと統計量の関係をアニメーションで可視化
前節では、パラメータと分布の関係を確認しました。次は、パラメータと統計量と歪度・尖度の関係をアニメーションで確認します。
・作図コード(クリックで展開)
パラメータ$\phi$の値を作成して、それぞれ歪度と尖度を計算し、フレーム切替用のラベルを作成します。
# パラメータとして利用する値を作成 phi_vals <- seq(from = 0, to = 1, by = 0.01) length(phi_vals) # フレーム数 # xがとり得る値を作成 x_vals <- 0:1 # 歪度を計算 skewness_vec <- (1 - 2 * phi_vals) / sqrt(phi_vals * (1 - phi_vals)) # 尖度を計算 kurtosis_vec <- (1 - 6 * phi_vals * (1 - phi_vals)) / (phi_vals * (1 - phi_vals)) # ラベル用のテキストを作成 label_vec <- paste0( "phi=", phi_vals, ", skewness=", round(skewness_vec, 3), ", kurtosis=", round(kurtosis_vec, 3) ) head(label_vec)
## [1] "phi=0, skewness=Inf, kurtosis=Inf" ## [2] "phi=0.01, skewness=9.849, kurtosis=95.01" ## [3] "phi=0.02, skewness=6.857, kurtosis=45.02" ## [4] "phi=0.03, skewness=5.51, kurtosis=28.364" ## [5] "phi=0.04, skewness=4.695, kurtosis=20.042" ## [6] "phi=0.05, skewness=4.129, kurtosis=15.053"
パラメータごとに歪度と尖度を計算して、フレーム切替用のラベルとして利用します。
パラメータごとに分布を計算します。
# パラメータごとにベルヌーイ分布を計算 anime_prob_df <- tidyr::expand_grid( x = x_vals, phi = phi_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(phi, x) |> # パラメータごとに並べ替え dplyr::mutate( probability = dbinom(x = x, size = 1, prob = phi), parameter = rep(label_vec, each = length(x_vals)) |> factor(levels = label_vec) # フレーム切替用ラベル ) # 確率を計算 anime_prob_df
## # A tibble: 202 × 4 ## x phi probability parameter ## <int> <dbl> <dbl> <fct> ## 1 0 0 1 phi=0, skewness=Inf, kurtosis=Inf ## 2 1 0 0 phi=0, skewness=Inf, kurtosis=Inf ## 3 0 0.01 0.99 phi=0.01, skewness=9.849, kurtosis=95.01 ## 4 1 0.01 0.01 phi=0.01, skewness=9.849, kurtosis=95.01 ## 5 0 0.02 0.98 phi=0.02, skewness=6.857, kurtosis=45.02 ## 6 1 0.02 0.02 phi=0.02, skewness=6.857, kurtosis=45.02 ## 7 0 0.03 0.97 phi=0.03, skewness=5.51, kurtosis=28.364 ## 8 1 0.03 0.03 phi=0.03, skewness=5.51, kurtosis=28.364 ## 9 0 0.04 0.96 phi=0.04, skewness=4.695, kurtosis=20.042 ## 10 1 0.04 0.04 phi=0.04, skewness=4.695, kurtosis=20.042 ## # … with 192 more rows
「パラメータと分布の関係」と同様に処理します。フレーム切替用のラベル(parameter列)の値にはlabel_vecを使います。
パラメータごとに統計量を計算します。
# パラメータごとに統計量を計算 anime_stat_df <- tibble::tibble( mean = phi_vals, # 期待値 sd = sqrt(phi_vals * (1 - phi_vals)), # 標準偏差 mode = max.col(cbind(1 - phi_vals, phi_vals)) - 1, # 最頻値 parameter = factor(label_vec, levels = label_vec) # フレーム切替用のラベル ) |> # 統計量を計算 dplyr::mutate( sd_m = mean - sd, sd_p = mean + sd ) |> # 期待値±標準偏差を計算 dplyr::select(!sd) |> # 不要な列を削除 tidyr::pivot_longer( cols = !parameter, names_to = "type", values_to = "statistic" ) |> # 統計量の列をまとめる dplyr::mutate( type = stringr::str_replace(type, pattern = "sd_.", replacement = "sd")) # 期待値±標準偏差のカテゴリを統一 anime_stat_df
## # A tibble: 404 × 3 ## parameter type statistic ## <fct> <chr> <dbl> ## 1 phi=0, skewness=Inf, kurtosis=Inf mean 0 ## 2 phi=0, skewness=Inf, kurtosis=Inf mode 0 ## 3 phi=0, skewness=Inf, kurtosis=Inf sd 0 ## 4 phi=0, skewness=Inf, kurtosis=Inf sd 0 ## 5 phi=0.01, skewness=9.849, kurtosis=95.01 mean 0.01 ## 6 phi=0.01, skewness=9.849, kurtosis=95.01 mode 0 ## 7 phi=0.01, skewness=9.849, kurtosis=95.01 sd -0.0895 ## 8 phi=0.01, skewness=9.849, kurtosis=95.01 sd 0.109 ## 9 phi=0.02, skewness=6.857, kurtosis=45.02 mean 0.02 ## 10 phi=0.02, skewness=6.857, kurtosis=45.02 mode 0 ## # … with 394 more rows
期待値・標準偏差・最頻値を計算して、さらに期待値から標準偏差を引いた値と足した値を求めます。
期待値・標準偏差の和と差・最頻値の4つの列をpivot_longer()でまとめます。列名がラベルになるので、標準偏差の和と差のラベルを統一します。
折れ線と凡例テキストの設定用の名前付きベクトルを作成します。
# 凡例用の設定を作成:(数式表示用) color_vec <- c(mean = "blue", sd = "orange", mode = "chocolate") linetype_vec <- c(mean = "dashed", sd = "dotted", mode = "dashed") label_vec <- c(mean = expression(E(x)), sd = expression(E(x) %+-% sqrt(V(x))), mode = expression(mode(x)))
「グラフの作成」のときと同じ処理です。
統計量の情報を重ねた二項分布のアニメーション(gif画像)を作成します。
# 統計量を重ねたベルヌーイ分布のアニメーションを作図 anime_prob_graph <- ggplot(data = anime_stat_df, ) + geom_bar(data = anime_prob_df, mapping = aes(x = x, y = probability), stat = "identity", fill = "#00A968", color = "#00A968") + # 分布 geom_vline(data = anime_stat_df, mapping = aes(xintercept = statistic, color = type, linetype = type), size = 1) + # 統計量 gganimate::transition_manual(parameter) + # フレーム scale_x_continuous(breaks = x_vals, labels = x_vals) + # x軸目盛 scale_linetype_manual(values = linetype_vec, labels = label_vec, name = "statistic") + # 線の種類:(線指定と数式表示用) scale_color_manual(values = color_vec, labels = label_vec, name = "statistic") + # 線の色:(色指定と数式表示用) theme(legend.text.align = 0) + # 図の体裁:凡例 labs(title = "Bernoulii Distribution", subtitle = "{current_frame}", x = "x", y = "probability") # ラベル # gif画像を作成 gganimate::animate(anime_prob_graph, nframes = length(phi_vals), fps = 100, width = 700, height = 600)
この記事では、ベルヌーイ分布の作図を確認しました。次は、乱数を生成します。
参考文献
- 岩田具治『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社,2015年.
おわりに
加筆修正の際に青トピシリーズから独立させました。
元記事的にも現記事的にも確率分布ネタシリーズ1つ目の記事です。これの元を書いたのは2年半前です。懐かしい、あの頃の自分に色々教えてあげたい。
- 2022.06.29：加筆修正しました。その際に「計算」と「乱数生成」の内容を分割しました。
【次の内容】
