はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
ポアソン分布のグラフをR言語で作成します。
【前の内容】
【他の記事一覧】
【この記事の内容】
Rでポアソン分布の作図
ポアソン分布(Poisson Distribution)のグラフを作成します。ポアソン分布については「ポアソン分布の定義式の確認 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用するパッケージ library(tidyverse) library(gganimate)
分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。
定義式の確認
まずは、ポアソン分布の定義式を確認します。
ポアソン分布は、次の式で定義されます。
ここで、$x$は単位時間における事象の発生回数、$\lambda$は発生回数の期待値です。
確率変数の値$x$は0以上の整数となります。パラメータ$\lambda$は、$\lambda > 0$を満たす必要があります。
ポアソン分布の平均と分散は、どちらもパラメータ$\lambda$になります。
これらの計算を行いグラフを作成します。
グラフの作成
ggplot2パッケージを利用してポアソン分布のグラフを作成します。
ポアソン分布のパラメータを設定します。
# パラメータを指定 lambda <- 4
発生回数の期待値$\lambda > 0$を指定します。
ポアソン分布の確率変数がとり得る値を作成します。
# 作図用のxの点を作成 x_vals <- seq(from = 0, to = ceiling(lambda) * 4) x_vals[1:5]
## [1] 0 1 2 3 4
$x$がとり得る値を作成して、x_valsとします。この例では、0からlambdaの4倍を範囲とします。ただし、$x$は非負の整数なので、ceiling()で値を切り上げて使います。
$x$の値ごとに確率を計算します。
# ポアソン分布を計算 prob_df <- tidyr::tibble( x = x_vals, # 確率変数 probability = dpois(x = x_vals, lambda = lambda) # 確率 ) head(prob_df)
## # A tibble: 6 x 2 ## x probability ## <int> <dbl> ## 1 0 0.0183 ## 2 1 0.0733 ## 3 2 0.147 ## 4 3 0.195 ## 5 4 0.195 ## 6 5 0.156
ポアソン分布の確率は、dpois()で計算できます。変数の引数xにx、パラメータの引数lambdaにlambdaを指定します。
x_valsと、x_valsの各要素に対応する確率をデータフレームに格納します。
ポアソン分布のグラフを作成します。
# ポアソン分布を作図 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軸目盛 labs(title = "Poisson Distribution", subtitle = paste0("lambda=", lambda)) # ラベル
ポアソン分布のグラフを描画できました。
この分布に平均と標準偏差の情報を重ねて表示します。
# 補助線用の統計量の計算 E_x <- lambda s_x <- sqrt(lambda) # 統計量を重ねたポアソン分布を作図 ggplot(data = prob_df, mapping = aes(x = x, y = probability)) + # データ geom_bar(stat = "identity", fill = "#00A968") + # 分布 geom_vline(xintercept = E_x, color = "orange", 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") + # 平均 + 標準偏差 scale_x_continuous(breaks = x_vals, labels = x_vals) + # x軸目盛 labs(title = "Poisson Distribution", subtitle = paste0("lambda=", lambda)) # ラベル
平均(破線)となる確率が最大なのを確認できます。また、平均を中心に標準偏差の範囲を点線で示しています。
パラメータと分布の関係
パラメータと分布の形状の関係を可視化します。
並べて可視化
複数のパラメータを指定してグラフを比較します。
・作図コード(クリックで展開)
複数個の$\lambda$の値を指定して、分布を計算します。
# パラメータを指定 lambda_vals <- c(1, 3.5, 8, 15) # 作図用のxの点を作成 x_vals <- seq(from = 0, to = ceiling(max(lambda_vals)) * 2) # lambdaの値ごとに分布を計算 res_prob_df <- tidyr::tibble() for(lambda in lambda_vals) { # ポアソン分布の情報を格納 tmp_prob_df <- tidyr::tibble( x = x_vals, # 確率変数 probability = dpois(x = x_vals, lambda = lambda), # 確率 parameter = paste0("lambda=", lambda) %>% as.factor() # 作図用のラベル ) # 結果を結合 res_prob_df <- rbind(res_prob_df, tmp_prob_df) } head(res_prob_df)
## # A tibble: 6 x 3 ## x probability parameter ## <int> <dbl> <fct> ## 1 0 0.368 lambda=1 ## 2 1 0.368 lambda=1 ## 3 2 0.184 lambda=1 ## 4 3 0.0613 lambda=1 ## 5 4 0.0153 lambda=1 ## 6 5 0.00307 lambda=1
for()ループを使って、lambda_valsの値を順番に取り出して確率を計算して、res_prob_dfに追加していきます。
ポアソン分布のグラフを作成します。
# ポアソン分布のグラフを作成 ggplot(data = res_prob_df, mapping = aes(x = x, y = probability, fill = parameter, color = parameter)) + # データ geom_bar(stat = "identity", position = "dodge") + # 棒グラフ #scale_x_continuous(breaks = x_vals, labels = x_vals) + # x軸目盛 labs(title = "Poisson Distribution") # ラベル
棒グラフを並べると分かりにくいので、散布図と折れ線グラフを重ねて可視化します。
# ポアソン分布のグラフを作成 ggplot(data = res_prob_df, mapping = aes(x = x, y = probability, fill = parameter, color = parameter)) + # データ geom_point(size = 2) + # 散布図 geom_line(size = 1) + # 折れ線グラフ #scale_x_continuous(breaks = x_vals, labels = x_vals) + # x軸目盛 labs(title = "Poisson Distribution") # ラベル
発生回数の期待値$\lambda$が大きいと、発生回数$x$が大きいほど確率が高くなる(山が右に位置する)のが分かります。このことは、平均の計算式からも分かります。
また分散の計算式からも分かるように、$\lambda$が大きいほど分布の裾が広く確率の最大値が小さく(なだらかな山に)なります。
アニメーションによる可視化
続いて、$\lambda$の値を少しずつ変更して、分布の変化をアニメーションで確認します。
・作図コード(クリックで展開)
# lambdaとして利用する値を指定 lambda_vals <- seq(from = 0, to = 10, by = 0.1) #length(lambda_vals) # フレーム数 # 作図用のxの点を作成 x_vals <- seq(from = 0, to = ceiling(max(lambda_vals)) * 2) # lambdaの値ごとに分布を計算 anime_prob_df <- tidyr::tibble() for(lambda in lambda_vals) { # ポアソン分布を計算 tmp_prob_df <- tidyr::tibble( x = x_vals, # 確率変数 probability = dpois(x = x_vals, lambda = lambda), # 確率 parameter = paste0("lambda=", lambda) %>% 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 lambda=0 ## 2 1 0 lambda=0 ## 3 2 0 lambda=0 ## 4 3 0 lambda=0 ## 5 4 0 lambda=0 ## 6 5 0 lambda=0
$\lambda$がとり得る値を作成してlambda_valsとします。
for()ループを使って、lambda_valsの値を順番に取り出して確率を計算して、anime_prob_dfに追加していきます。
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 = "Poisson Distribution", subtitle = "{current_frame}") # ラベル # gif画像を作成 gganimate::animate(anime_prob_graph, nframes = length(lambda_vals), fps = 100)
$\lambda$が大きくなるに従って、$x$が大きいほど確率が高くなる(山が右に移動する)のを確認できます。
以上で、ポアソン分布のグラフ作成を確認できました。次は、乱数を生成します。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
青トピに続けて緑と黄色のベイズ本で登場する分布をやっていきます。
- 2022.02.21：追記して、記事を分割しました。
【次の内容】