はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
この記事では、R言語でガンマ分布の乱数を生成します。
【前の内容】
【他の記事一覧】
【この記事の内容】
ガンマ分布の乱数生成
ガンマ分布(Binomial Distribution)の乱数を生成します。ガンマ分布については「【R】ガンマ分布の作図 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用パッケージ library(tidyverse) library(gganimate)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。
magrittrパッケージのパイプ演算子%>%ではなく、ベースパイプ(ネイティブパイプ)演算子|>を使っています。%>%に置き換えても処理できます。
分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。
サンプリング
ガンマ分布の乱数を生成します。
ガンマ分布のパラメータ$\lambda$とデータ数$N$を設定します。
# パラメータを指定 a <- 5 b <- 2 # データ数(サンプルサイズ)を指定 N <- 1000
$\lambda > 0$とデータ数(サンプルサイズ)$N$を指定します。
ガンマ分布に従う乱数を生成します。
# ガンマ分布に従う乱数を生成 lambda_n <- rgamma(n = N, shape = a, rate = b) head(lambda_n)
## [1] 1.377111 1.157740 5.666513 2.121982 4.624966 3.183849
ガンマ分布の乱数は、rgamma()で生成できます。データ数(サンプルサイズ)の引数nにN、形状の引数shapeにa、尺度の引数rateにbを指定します。
乱数の可視化
サンプルのヒストグラムを作成します。
# サンプルを格納 data_df <- tidyr::tibble(lambda = lambda_n) # サンプルのヒストグラムを作成:度数 ggplot(data = data_df, mapping = aes(x = lambda, y = ..count..)) + # データ geom_histogram(fill = "#00A968", bins = 30) + # 度数 labs(title = "Gamma Distribution", subtitle = paste0("a=", a, ", b=", b, ", N=", N), x = expression(lambda), y = "frequency") # ラベル
サンプルをデータフレームに格納して、geom_histgram()でヒストグラムを描画します。デフォルト(y = ..count..)では、度数のヒストグラムを作成します。集計の範囲については、バーの数の引数binsまたはバーのサイズの引数binwidthを指定します。
サンプルの密度を確率分布と重ねて描画します。
# lambdaの値を作成 lambda_vals <- seq(from = 0, to = max(lambda_n)+1, length.out = 250) # ガンマ分布を計算 dens_df <- tidyr::tibble( lambda = lambda_vals, # 確率変数 density = dgamma(x = lambda_vals, shape = a, rate = b) # 確率密度 ) # サンプルのヒストグラムを作成:密度 ggplot() + geom_histogram(data = data_df, mapping = aes(x = lambda, y = ..density..), fill = "#00A968", bins = 30) + # 密度 geom_line(data = dens_df, mapping = aes(x = lambda, y = density), color = "darkgreen", size = 1, linetype = "dashed") + # 元の分布 labs(title = "Gamma Distribution", subtitle = paste0("a=", a, ", b=", b, ", N=", N), x = expression(lambda), y = "density") # ラベル
geom_histogram()のy軸の引数yに..density..を指定すると、密度に変換したヒストグラムを描画します。
データ数(サンプルサイズ)が十分に増えると、ヒストグラムの形が分布の形に近付きます。
乱数と分布の関係をアニメーションで可視化
次は、サンプルサイズとヒストグラムの形状の関係をアニメーションで確認します。
・作図コード(クリックで展開)
パラメータ$\lambda$とデータ数$N$を指定して、ガンマ分布の乱数を生成します。
# パラメータを指定 a <- 5 b <- 2 # データ数(フレーム数)を指定 N <- 100 # ガンマ分布に従う乱数を生成 lambda_n <- rgamma(n = N, shape = a, rate = b) head(lambda_n)
## [1] 1.838699 2.718451 2.650305 1.599120 1.664898 2.952641
lambda_nのn番目の要素を、n回目にサンプリングされた値とみなします。アニメーションのn番目のフレームでは、n個のサンプルlambda_n[1:n]のヒストグラムを描画します。
サンプリング回数ごとに、それまでのサンプルを持つデータフレームを作成します。
# サンプルを複製して格納 anime_freq_df <- tibble::tibble( lambda = rep(lambda_n, times = N), # サンプル n = rep(1:N, times = N), # データ番号 frame = rep(1:N, each = N) # フレーム番号 ) |> # 全ての組み合わせを作成 dplyr::filter(n <= frame) |> # サンプリング回数以前のサンプルを抽出 dplyr::mutate( parameter = paste0("a=", a, ", b=", b, ", N=", frame) |> factor(levels = paste0("a=", a, ", b=", b, ", N=", 1:N)) ) # フレーム切替用ラベルを追加 anime_freq_df
## # A tibble: 5,050 × 4 ## lambda n frame parameter ## <dbl> <int> <int> <fct> ## 1 1.84 1 1 a=5, b=2, N=1 ## 2 1.84 1 2 a=5, b=2, N=2 ## 3 2.72 2 2 a=5, b=2, N=2 ## 4 1.84 1 3 a=5, b=2, N=3 ## 5 2.72 2 3 a=5, b=2, N=3 ## 6 2.65 3 3 a=5, b=2, N=3 ## 7 1.84 1 4 a=5, b=2, N=4 ## 8 2.72 2 4 a=5, b=2, N=4 ## 9 2.65 3 4 a=5, b=2, N=4 ## 10 1.60 4 4 a=5, b=2, N=4 ## # … with 5,040 more rows
データ番号(n列)は1からNの整数がN回繰り返すように、フレーム番号(frame列)は1からNの整数がN個ずつ並ぶように作成します。サンプルはデータ番号に対応するように複製します。
フレーム番号ごとに、それ以下の番号のサンプルをfilter()で抽出します。
フレーム番号をデータ番号として、フレーム切替用のラベル列を作成します。これにより、各フレーム(サンプリング回数)ごとにそれ以前のサンプルを使ってグラフを作成できます。ラベルが文字列型だと文字列の基準で順序が決まるので、因子型にしてサンプリング回数に応じたレベル(順序)を設定します。
このデータフレームは、ヒストグラムを描画するのに使います。
サンプルと対応するラベルをデータフレームに格納します。
# サンプルを格納 anime_data_df <- tidyr::tibble( lambda = lambda_n, # サンプル parameter = paste0("a=", a, ", b=", b, ", N=", 1:N) |> factor(levels = paste0("a=", a, ", b=", b, ", N=", 1:N)) # フレーム切替用ラベル ) anime_data_df
## # A tibble: 100 × 2 ## lambda parameter ## <dbl> <fct> ## 1 1.84 a=5, b=2, N=1 ## 2 2.72 a=5, b=2, N=2 ## 3 2.65 a=5, b=2, N=3 ## 4 1.60 a=5, b=2, N=4 ## 5 1.66 a=5, b=2, N=5 ## 6 2.95 a=5, b=2, N=6 ## 7 5.32 a=5, b=2, N=7 ## 8 1.56 a=5, b=2, N=8 ## 9 0.765 a=5, b=2, N=9 ## 10 2.58 a=5, b=2, N=10 ## # … with 90 more rows
このデータフレームは、各サンプリングにおけるサンプルの値を描画するのに使います。
サンプルのヒストグラムのアニメーション(gif画像)を作成します。
# サンプルのヒストグラムを作成:度数 anime_freq_graph <- ggplot() + geom_histogram(data = anime_freq_df, mapping = aes(x = lambda, y = ..count..), breaks = seq(from = 0, to = max(lambda_vals), length.out = 30), fill = "#00A968") + # 度数 geom_point(data = anime_data_df, mapping = aes(x = lambda, y = 0), color = "orange", size = 5) + # サンプル gganimate::transition_manual(parameter) + # フレーム labs(title = "Gamma Distribution", subtitle = "{current_frame}", x = expression(lambda), y = "frequency") # ラベル # gif画像を作成 gganimate::animate(anime_freq_graph, nframes = N, fps = 10, width = 800, height = 600)
集計範囲を固定するために、geom_histogram()のbreaks引数に区切り位置を指定します。この例では、seq()を使って0からlambda_valsの最大値の範囲を30等分しています。
transition_manual()にフレームの順序を表す列を指定します。この例では、因子型のラベルのレベルの順に描画されます。
animate()のフレーム数の引数nframesにデータ数(サンプルサイズ)、フレームレートの引数fpsに1秒当たりのフレーム数を指定します。fps引数の値が大きいほどフレームが早く切り替わります。
サンプルの密度を確率分布と重ねたアニメーションを作成します。
# lambdaの値を作成 lambda_vals <- seq(from = 0, to = max(lambda_n)+1, length.out = 500) # ガンマ分布を計算 dens_df <- tidyr::tibble( lambda = lambda_vals, # 確率変数 density = dgamma(x = lambda_vals, shape = a, rate = b) # 確率密度 ) # サンプルのヒストグラムを作成:密度 anime_freq_graph <- ggplot() + geom_histogram(data = anime_freq_df, mapping = aes(x = lambda, y = ..density..), breaks = seq(from = 0, to = max(lambda_vals), length.out = 30), fill = "#00A968") + # 密度 geom_line(data = dens_df, mapping = aes(x = lambda, y = density), color = "darkgreen", size = 1, linetype = "dashed") + # 元の分布 geom_point(data = anime_data_df, mapping = aes(x = lambda, y = 0), color = "orange", size = 5) + # サンプル gganimate::transition_manual(parameter) + # フレーム coord_cartesian(ylim = c(NA, max(dens_df[["density"]]*2))) + # 軸の表示範囲 labs(title = "Gamma Distribution", subtitle = "{current_frame}", x = expression(lambda), y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_freq_graph, nframes = N, fps = 10, width = 800, height = 600)
サンプルが増えるに従って、ヒストグラムが元の分布に近付くのを確認できます。
この記事では、ガンマ分布の乱数生成を確認しました。次は、ガンマ分布から確率分布を生成します。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
加筆修正の際に「ガンマ分布の作図」から記事を分割しました。
【次の内容】
