はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
この記事では、ガンマ分布からポアソン分布と1次元ガウス分布(1変量正規分布)を生成します。
【前の内容】
【他の記事一覧】
【この記事の内容】
ガンマ分布から確率分布の生成
ガンマ分布(Gamma Distribution)からポアソン分布(Poisson Distribution)と1次元ガウス分布(gaussian Distribution)を生成します。ガンマ分布については「【R】ガンマ分布の作図 - からっぽのしょこ」、ポアソン分布については「ポアソン分布の定義式の確認 - からっぽのしょこ」、1次元ガウス分布については「【R】1次元ガウス分布の作図 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用パッケージ library(tidyverse) library(patchwork)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。
magrittrパッケージのパイプ演算子%>%ではなく、ベースパイプ(ネイティブパイプ)演算子|>を使っています。%>%に置き換えても処理できます。
生成分布の設定
まずは、パラメータの生成分布としてガンマ分布を設定して、ポアソン分布と1次元ガウス分布のパラメータを生成(サンプリング)します。
ガンマ分布のパラメータ$a, b$とサンプルサイズ$N$を設定します。
# パラメータを指定 a <- 5 b <- 2 # サンプルサイズを指定 N <- 10
$a > 0, b > 0$の値とサンプルサイズ(パラメータ数)$N$を指定します。
ガンマ分布に従う乱数を生成します。
# ポアソン分布・ガウス分布のパラメータを生成 lambda_n <- rgamma(n = N, shape = a, rate = b) |> # ガンマ分布の乱数を生成 sort() # 昇順に並べ替え # パラメータを格納 param_df <- tibble::tibble( lambda = lambda_n, parameter = round(lambda_n, 3) |> factor(), # 色分け用ラベル param_rev = round(lambda_n, 3) |> factor(levels = round(rev(lambda_n), 3)) # 色分け用ラベル:逆順 ) param_df
## # A tibble: 10 × 3 ## lambda parameter param_rev ## <dbl> <fct> <fct> ## 1 1.09 1.093 1.093 ## 2 1.28 1.277 1.277 ## 3 1.48 1.484 1.484 ## 4 1.68 1.678 1.678 ## 5 2.37 2.369 2.369 ## 6 2.90 2.901 2.901 ## 7 3.20 3.201 3.201 ## 8 3.40 3.399 3.399 ## 9 4.61 4.608 4.608 ## 10 5.50 5.497 5.497
ガンマ分布の乱数は、rgamma()で生成できます。データ数(サンプルサイズ)の引数nにN、形状の引数shapeにa、尺度の引数rateにbを指定します。
作図の都合上、降順にレベル付けしたラベル列も作成しておきます。
生成した値をポアソン分布とガウス分布のパラメータ$\lambda$として使います。
ガンマ分布を計算します。
# lambdaの値を作成 lambda_vals <- seq(from = 0, to = a/b * 4, length.out = 501) # ガンマ分布を計算 gamma_df <- tidyr::tibble( lambda = lambda_vals, # 確率変数 density = dgamma(x = lambda_vals, shape = a, rate = b) # 確率密度 ) gamma_df
## # A tibble: 501 × 2 ## lambda density ## <dbl> <dbl> ## 1 0 0 ## 2 0.02 0.000000205 ## 3 0.04 0.00000315 ## 4 0.06 0.0000153 ## 5 0.08 0.0000465 ## 6 0.1 0.000109 ## 7 0.12 0.000217 ## 8 0.14 0.000387 ## 9 0.16 0.000635 ## 10 0.18 0.000977 ## # … with 491 more rows
ガンマ分布の確率変数がとり得る値$\lambda > 0$を作成して、確率密度を計算します。
ガンマ分布の確率密度は、dgamma()で計算できます。確率変数の引数xにlambda_vals、形状の引数shapeにa、尺度の引数rateにbを指定します。
lambda_valsと、lambda_valsの各要素に対応する確率密度をデータフレームに格納します。
ガンマ分布の期待値$\mathbb{E}[\lambda]$を計算します。
# ガンマ分布の期待値を計算 E_lambda <- a / b E_lambda
## [1] 2.5
サンプリングされたパラメータとその分布の基準を示すのに使います。
生成分布(ガンマ分布)とパラメータのサンプルをプロットします。
# 生成分布とパラメータのサンプルを作図:通常 gamma_graph <- ggplot() + geom_vline(xintercept = E_lambda, color = "red", size = 1, linetype = "dashed") + # パラメータの期待値 geom_line(data = gamma_df, mapping = aes(x = lambda, y = density), color = "#00A968", size = 1) + # パラメータの生成分布 geom_point(data = param_df, mapping = aes(x = lambda, y = 0, color = parameter), size = 6, alpha = 0.8, show.legend = FALSE) + # パラメータのサンプル guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))) + # 凡例の体裁 labs(title = "Gamma Distribution", subtitle = paste0("a=", a, ", b=", b), color = expression(lambda), x = expression(lambda), y = "density") # ラベル gamma_graph
期待値を破線で示します。
x軸とy軸を入れ替えたグラフも作成しておきます。
# 生成分布とパラメータのサンプルを作図:軸の入替 gamma_graph_flip <- ggplot() + geom_vline(xintercept = E_lambda, color = "red", size = 1, linetype = "dashed") + # パラメータの期待値 geom_line(data = gamma_df, mapping = aes(x = lambda, y = density), color = "#00A968", size = 1) + # パラメータの生成分布 geom_point(data = param_df, mapping = aes(x = lambda, y = 0, color = param_rev), alpha = 0.8, size = 6, show.legend = FALSE) + # パラメータのサンプル coord_flip() + # 軸の入れ替え labs(title = "Gamma Distribution", subtitle = paste0("a=", a, ", b=", b), color = expression(lambda), x = expression(lambda), y = "density") # ラベル gamma_graph_flip
coord_flip()で軸を入れ替えられます。
以上で、生成分布を設定して、パラメータを生成しました。次は、パラメータのサンプルを用いて、ポアソン分布とガウス分布を作図します。
分布の作図：ポアソン分布
次に、生成した値をポアソン分布のパラメータとして利用します。グラフ作成については「【R】ポアソン分布の作図 - からっぽのしょこ」を参照してください。
パラメータのサンプルごとにポアソン分布を計算します。
# xの値を作成 x_vals <- seq(from = 0, to = ceiling(E_lambda) * 4) x_vals <- seq(from = 0, to = ceiling(E_lambda*4)) # パラメータのサンプルごとにガウス分布を計算 res_poisson_df <- tidyr::expand_grid( x = x_vals, lambda = lambda_n ) |> dplyr::arrange(lambda, x) |> # パラメータごとに並べ替え dplyr::mutate( probability = dpois(x = x, lambda = lambda), # 確率密度 parameter = round(lambda, 3) |> factor(levels = round(lambda_n, 3)) # 色分け用ラベル ) res_poisson_df
## # A tibble: 110 × 4 ## x lambda probability parameter ## <int> <dbl> <dbl> <fct> ## 1 0 1.09 0.335 1.093 ## 2 1 1.09 0.366 1.093 ## 3 2 1.09 0.200 1.093 ## 4 3 1.09 0.0730 1.093 ## 5 4 1.09 0.0199 1.093 ## 6 5 1.09 0.00436 1.093 ## 7 6 1.09 0.000795 1.093 ## 8 7 1.09 0.000124 1.093 ## 9 8 1.09 0.0000170 1.093 ## 10 9 1.09 0.00000206 1.093 ## # … with 100 more rows
ポアソン分布の確率変数がとり得る値(非負の整数)$x$を作成してx_valsとします。
x_valsとlambda_nそれぞれの要素の全ての組み合わせをexpand_grid()で作成して、組み合わせごとに確率を計算します。
ポアソン分布の確率は、dpois()で計算できます。確率変数の引数xにx_valsの値、パラメータの引数lambdaにlambda_nの値を指定します。
ガンマ分布の期待値(ポアソン分布のパラメータの期待値)を用いて、ガンマ分布を計算します。
# パラメータの期待値によるガンマ分布を計算 E_poisson_df <- tidyr::tibble( x = x_vals, # 確率変数 probability = dpois(x = x_vals, lambda = E_lambda) # 確率 ) E_poisson_df
## # A tibble: 11 × 2 ## x probability ## <int> <dbl> ## 1 0 0.0821 ## 2 1 0.205 ## 3 2 0.257 ## 4 3 0.214 ## 5 4 0.134 ## 6 5 0.0668 ## 7 6 0.0278 ## 8 7 0.00994 ## 9 8 0.00311 ## 10 9 0.000863 ## 11 10 0.000216
lambda引数にE_lambdaを指定します。
$N + 1$個のポアソン分布を作図します。
# サンプルと期待値によるポアソン分布を作図 poisson_graph <- ggplot() + geom_line(data = E_poisson_df, mapping = aes(x = x, y = probability), color = "red", size = 1, linetype ="dashed") + # 期待値による分布 geom_point(data = E_poisson_df, mapping = aes(x = x, y = probability), color = "red", size = 3) + # 期待値による分布 geom_line(data = res_poisson_df, mapping = aes(x = x, y = probability, color = parameter), alpha = 0.8, size = 1) + # サンプルによる分布 geom_point(data = res_poisson_df, mapping = aes(x = x, y = probability, color = parameter), alpha = 0.8, size = 3) + # サンプルによる分布 scale_x_continuous(breaks = x_vals, minor_breaks = FALSE) + # x軸目盛 labs(title = "Poisson Distribution", subtitle = parse(text = paste0("E(lambda)==", round(E_lambda, 3))), color = expression(lambda), x = "x") # ラベル poisson_graph
(棒グラフだと分かりにくいので)散布図と折れ線グラフで分布を描画します。期待値による分布を破線で示します。
パラメータの生成分布(ガンマ分布)と生成された分布(ポアソン分布)を並べて描画します。
# グラフを並べて描画 gamma_graph / poisson_graph + patchwork::plot_layout(guides = "collect")
patchworkパッケージの/演算子を使うと上下に並べて描画できます。
ポアソン分布の期待値はパラメータ$\lambda$、最頻値は$\lambda$以下の最大の整数になります。パラメータのサンプルの点と、分布のピークが対応しているのが分かります。
分布の作図：1次元ガウス分布
続いて、生成した値を1次元ガウス分布の精度パラメータ(分散パラメータの逆数)として利用します。グラフ作成については「【R】1次元ガウス分布の作図 - からっぽのしょこ」を参照してください。
平均パラメータ$\mu$を指定します。また、標準偏差パラメータの期待値を計算します。
# 平均パラメータを指定 mu = 0 # 標準偏差の期待値を計算 E_sigma <- sqrt(1 / E_lambda) E_sigma
## [1] 0.6324555
精度$\lambda$の逆数$\frac{1}{\lambda}$が分散$\sigma^2$であり、分散$\sigma^2$の平方根$\sqrt{\sigma^2}$が標準偏差$\sigma$なので、$\sigma = \sqrt{\frac{1}{\lambda}}$です。よって、精度パラメータの期待値$\mathbb{E}[\lambda]$を用いて、標準偏差の期待値$\mathbb{E}[\sigma] = \sqrt{\frac{1}{\mathbb{E}[\lambda]}}$を計算します。
パラメータのサンプルごとにガウス分布を計算します。
# xの値を作成 x_vals <- seq(from = mu-E_sigma * 4, to = mu+E_sigma * 4, length.out = 500) # パラメータのサンプルごとにガウス分布を計算 res_gaussian_df <- tidyr::expand_grid( x = x_vals, lambda = lambda_n ) |> # 全ての組み合わせを作成 dplyr::arrange(lambda, x) |> # パラメータごとに並べ替え dplyr::mutate( density = dnorm(x = x, mean = mu, sd = sqrt(1/lambda)), # 確率密度 param_rev = round(lambda, 3) |> factor(levels = round(rev(lambda_n), 3)) # 色分け用ラベル:逆順 ) res_gaussian_df
## # A tibble: 5,000 × 4 ## x lambda density param_rev ## <dbl> <dbl> <dbl> <fct> ## 1 -2.53 1.09 0.0126 1.093 ## 2 -2.52 1.09 0.0130 1.093 ## 3 -2.51 1.09 0.0133 1.093 ## 4 -2.50 1.09 0.0137 1.093 ## 5 -2.49 1.09 0.0141 1.093 ## 6 -2.48 1.09 0.0145 1.093 ## 7 -2.47 1.09 0.0149 1.093 ## 8 -2.46 1.09 0.0153 1.093 ## 9 -2.45 1.09 0.0157 1.093 ## 10 -2.44 1.09 0.0162 1.093 ## # … with 4,990 more rows
ガウス分布の確率変数がとり得る値$x$を作成してx_valsとします。
x_valsとlambda_nそれぞれの要素の全ての組み合わせをexpand_grid()で作成して、組み合わせごとに確率を計算します。
ガウス分布の確率密度は、dnorm()で計算できます。確率変数の引数xにx_valsの値、平均の引数meanにmu、標準偏差の引数sdにlambda_nの値の逆数の平方根を指定します。
ガンマ分布の期待値(ガウス分布の精度パラメータの期待値)を用いて、ガウス分布を計算します。
# パラメータの期待値によるガウス分布を計算 E_gaussian_df <- tidyr::tibble( x = x_vals, # 確率変数 density = dnorm(x = x_vals, mean = mu, sd = E_sigma) # 確率密度 ) E_gaussian_df
## # A tibble: 500 × 2 ## x density ## <dbl> <dbl> ## 1 -2.53 0.000212 ## 2 -2.52 0.000226 ## 3 -2.51 0.000240 ## 4 -2.50 0.000256 ## 5 -2.49 0.000273 ## 6 -2.48 0.000291 ## 7 -2.47 0.000309 ## 8 -2.46 0.000329 ## 9 -2.45 0.000351 ## 10 -2.44 0.000373 ## # … with 490 more rows
sd引数にE_sigmaを指定します。
$N + 1$個のガウス分布を作図します。
# サンプルによるガウス分布を作図 gaussian_graph <- ggplot() + geom_line(data = E_gaussian_df, mapping = aes(x = x, y = density,), color = "red", size = 1, linetype = "dashed") + # 期待値による分布 geom_line(data = res_gaussian_df, mapping = aes(x = x, y = density, color = param_rev), alpha = 0.8, size = 1) + # サンプルによる分布 theme(legend.text.align = 0) + # 図の体裁:凡例 labs(title = "Gaussian Distribution", subtitle = parse(text = paste0("list(E(lambda)==", round(E_lambda, 3), ", mu==", mu, ")")), color = expression(lambda), x = "x", y = "density") # ラベル gaussian_graph
期待値による分布を破線で示します。
パラメータの生成分布(ガンマ分布)と生成された分布(ガウス分布)を並べて描画します。
# グラフを並べて描画 gamma_graph_flip + gaussian_graph + patchwork::plot_layout(widths = c(1, 2), guides = "collect")
patchworkパッケージの+演算子を使うと左右に並べて描画できます。
精度$\lambda$が大きいほど、分散が小さくなるので、先の細い形状になります。(軸を入れ替えた)パラメータのサンプルの点と、分布のピークが対応しているのが分かります(？)。
この記事では、ベータ分布からの分布生成を確認しました。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
加筆修正の際に「ガンマ分布の作図」から記事を分割しました。
【次の内容】
