からっぽのしょこ

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

【R】ポアソン分布の乱数生成

はじめに

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

 ポアソン分布の乱数の生成をR言語で行います。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

Rでポアソン分布の乱数生成

 ポアソン分布(Poisson Distribution)の乱数を生成します。ポアソン分布については「ポアソン分布の定義式の確認 - からっぽのしょこ」を参照してください。

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

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

 分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。

サンプリング

 ポアソン分布の乱数を生成します。

 ポアソン分布のパラメータ$\lambda$を指定しします。

# パラメータを指定
lambda <- 4

# データ数(サンプルサイズ)を指定
N <- 1000

 発生回数の期待値$\lambda > 0$とデータ数$N$を指定します。

 ポアソン分布に従う乱数を生成します。

# ポアソン分布に従う乱数を生成
x_n <- rpois(n = N, lambda = lambda)
x_n[1:10]
##  [1] 4 6 3 2 7 3 2 6 2 2

 ポアソン分布の乱数はrpois()で生成できます。データ数(サンプルサイズ)の引数nN、パラメータの引数lambdalmandaを指定します。

乱数の可視化

 サンプルのヒストグラムを作成します。

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

# 乱数を集計して格納
freq_df <- tidyr::tibble(x = x_n) %>% # 乱数を格納
  dplyr::count(x, name = "frequency") %>% # 頻度を測定
  dplyr::mutate(proportion = frequency / N) # 相対度数を計算
head(freq_df)
## # A tibble: 6 x 3
##       x frequency proportion
##   <int>     <int>      <dbl>
## 1     0        20      0.02 
## 2     1        77      0.077
## 3     2       163      0.163
## 4     3       213      0.213
## 5     4       202      0.202
## 6     5       142      0.142

 サンプリングした値x_nをデータフレームに格納します。
 count()で重複する値をカウントします。
 得られた頻度frequencyをデータ数Nで割り、各値の相対度数を計算してproportion列とします。

 作図に利用するために、ポアソン分布の情報を作成します。

# 作図用のxの点を作成
x_vals <- seq(from = 0, to = max(x_n) + 3)

# ポアソン分布を計算
prob_df <- tidyr::tibble(
  x = x_vals, # 確率変数
  probability = dpois(x = x_vals, lambda = lambda) # 確率
)

 ポアソン分布のグラフの作成については「作図」を参照してください。

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

# サンプルのヒストグラムを作成
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 = "Poisson Distribution", 
       subtitle = paste0("lambda=", lambda, ", N=", N)) # ラベル

f:id:anemptyarchive:20220221234326p:plain
ポアソン分布の乱数のヒストグラム:頻度


 相対度数を分布と重ねて描画します。

# サンプルの相対度数を作図
ggplot() + 
  geom_bar(data = freq_df, mapping = aes(x = x, y = proportion), 
           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軸目盛
  labs(title = "Poisson Distribution", 
       subtitle = paste0("lambda=", lambda, ", N=", N)) # ラベル

f:id:anemptyarchive:20220221234346p:plain
ポアソン分布の乱数のヒストグラム:相対頻度

 データ数が十分に増えると、ヒストグラムが分布の形に近付きます。

推移の可視化

 サンプルサイズとヒストグラムの変化をアニメーションで確認します。

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] <- rpois(n = 1, lambda = lambda)
  
  # n個の乱数を集計して格納
  tmp_freq_df <- tidyr::tibble(x = x_n[1:n]) %>% # 乱数を格納
    dplyr::count(x, name = "frequency") %>% # 頻度を測定
    dplyr::full_join(tidyr::tibble(x = x_vals), by = "x") %>% # サンプルにない値を補完
    dplyr::mutate(frequency = tidyr::replace_na(frequency, 0)) %>% # サンプルにない場合のNAを0に置換
    dplyr::mutate(proportion = frequency / n) # 相対度数を計算

  # ラベル用のテキストを作成
  label_text <- paste0(
    "lambda=", lambda, ", 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     4         1          1 lambda=4, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
## 2     0         0          0 lambda=4, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
## 3     1         0          0 lambda=4, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
## 4     2         0          0 lambda=4, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
## 5     3         0          0 lambda=4, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
## 6     5         0          0 lambda=4, N=1=(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~

 for()ループを使って、「x_nから1データずつ取り出す」または「1データずつ生成」して、結果をanime_***_dfに追加していきます。

 ただし、先ほどの集計処理に、サンプルx_nにない値のカウント(行)を追加する処理を行っています。
 x_valsの値を持つデータフレームを作成し、full_join()で結合します。
 サンプルにない場合は頻度列frequencyが欠損値NAになるので、replace_na()0に置換します。

 また、フレームごとに分布の情報prob_dfを複製し、対応するparameter列を追加する必要があります。

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

# アニメーション用のサンプルのヒストグラムを作成
anime_freq_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 = "Poisson Distribution", 
       subtitle = "{current_frame}") # ラベル

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


 相対度数のアニメーションを作成します。

# アニメーション用のサンプルの相対度数を作図
anime_prop_graph <- ggplot() + 
  geom_bar(data = anime_freq_df, mapping = aes(x = x, y = proportion), 
           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.01, 0.5)) + # y軸の表示範囲
   labs(title = "Poisson Distribution", 
       subtitle = "{current_frame}") # ラベル

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


f:id:anemptyarchive:20220221234410g:plainf:id:anemptyarchive:20220221234412g:plain
ポアソン分布の乱数の推移

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

全データをアニメーションで可視化

 続いて、フレーム数を指定して、フレームごとに均等に分けたサンプルを増やしてヒストグラムの変化をアニメーションで確認します。

・作図コード(クリックで展開)

# フレーム数を指定
frame_num <- 100

# 1フレーム当たりのデータ数を計算
n_per_frame <- N %/% frame_num

# フレームごとに乱数を抽出
anime_freq_df <- tidyr::tibble()
anime_prob_df <- tidyr::tibble()
for(n in 1:N_frame) {
  # 乱数を集計して格納
  tmp_freq_df <- tidyr::tibble(x = x_n[1:(n*n_per_frame)]) %>% # 乱数を格納
    dplyr::count(x, name = "frequency") %>% # 頻度を測定
    dplyr::full_join(tidyr::tibble(x = x_vals), by = "x") %>% # サンプルにない値を補完
    dplyr::mutate(frequency = tidyr::replace_na(frequency, 0)) %>% # サンプルにない場合のNAを0に置換
    dplyr::mutate(proportion = frequency / (n*n_per_frame)) # 相対度数を計算
  
  # ラベル用のテキストを作成
  label_text <- paste0(
    "lambda=", lambda, ", 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     2         4        0.4 lambda=4, N=10=(4, 2, 1, 2, 1, 0, 0, 0, 0, 0, 0, 0~
## 2     3         2        0.2 lambda=4, N=10=(4, 2, 1, 2, 1, 0, 0, 0, 0, 0, 0, 0~
## 3     4         1        0.1 lambda=4, N=10=(4, 2, 1, 2, 1, 0, 0, 0, 0, 0, 0, 0~
## 4     6         2        0.2 lambda=4, N=10=(4, 2, 1, 2, 1, 0, 0, 0, 0, 0, 0, 0~
## 5     7         1        0.1 lambda=4, N=10=(4, 2, 1, 2, 1, 0, 0, 0, 0, 0, 0, 0~
## 6     0         0        0   lambda=4, N=10=(4, 2, 1, 2, 1, 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", 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_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", fill = "#00A968") + # 相対度数
  geom_bar(data = anime_prob_df, mapping = aes(x = x, y = probability), 
           stat = "identity", alpha = 0, color = "darkgreen", linetype = "dashed") + # 元の分布
  gganimate::transition_manual(parameter) + # フレーム
  scale_x_continuous(breaks = x_vals, labels = x_vals) + # x軸目盛
  ylim(c(-0.01, 0.5)) + # y軸の表示範囲
  labs(title = "Poisson Distribution", 
       subtitle = "{current_frame}") # ラベル

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


f:id:anemptyarchive:20220221234452g:plainf:id:anemptyarchive:20220221234455g:plain
ポアソン分布の乱数の推移


 以上で、ポアソン分布を確認できました。

参考文献

  • 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.

おわりに

 記事を分割しました。

【次の内容】

つづく