からっぽのしょこ

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

【R】ガンマ分布の作図

はじめに

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

 ガンマ分布の計算とグラフの作成をR言語で行います。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

ガンマ分布の作図

 ガンマ分布(Gamma Distribution)の計算と作図を行います。

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

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

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

定義式の確認

 まずは、ガンマ分布の定義式を確認します。

 ガンマ分布は、次の式で定義されます。

$$ \mathrm{Gam}(\lambda | a, b) = \frac{b^a}{\Gamma(a)} \lambda^{a-1} e^{-b\lambda} $$

 ここで、$a$は形状に関するパラメータ、$b$は尺度に関するパラメータです。
 確率変数の値$\lambda$は、$\lambda > 0$となります。パラメータ$a, b$は、$a > 0, b > 0$を満たす必要があります。

 この式の対数をとると、次の式になります。

$$ \log \mathrm{Gam}(\lambda | a, b) = a \log b - \log \Gamma(a) + (a - 1) \log \lambda - b \lambda $$

 ガンマ分布の平均と・分散・最頻値は、それぞれ次の式で計算できます。

$$ \begin{aligned} \mathbb{E}[\lambda] &= \frac{a}{b} \\ \mathbb{V}[\lambda] &= \frac{a}{b^2} \\ \mathrm{mode}[\lambda] &= \frac{a - 1}{b} \end{aligned} $$


 ガンマ分布は、0以上の値$\lambda > 0$を生成することから、1次元ガウス分布の精度パラメータやポアソン分布のパラメータの事前分布として利用されます。

確率密度の計算

 ガンマ分布に従う確率密度を計算する方法をいくつか確認します。

 パラメータを設定します。

# パラメータを指定
a <- 2
b <- 2

# 確率変数の値を指定
lambda <- 2

 ガンマ分布のパラメータ$a > 0, b > 0$、確率変数がとり得る値$\lambda > 0$を指定します。設定した値に従う確率密度を計算します。

 まずは、定義式から確率密度を計算します。

# 定義式により確率密度を計算
C <- b^a / gamma(a)
dens <- C * lambda^(a - 1) * exp(-b * lambda)
dens
## [1] 0.1465251

 ガンマ分布の定義式

$$ \begin{aligned} C_{\mathrm{Gam}} &= \frac{b^a}{\Gamma(a)} \\ \mathrm{Gam}(\lambda | a, b) &= C_{\mathrm{Gam}} \lambda^{a-1} \exp(- b \lambda) \end{aligned} $$

で計算します。$C_{\mathrm{Gam}}$は、ガンマ分布の正規化係数です。
 $\Gamma(x)$はガンマ関数で、gamma()で計算できます。

 対数をとった定義式から計算します。

# 対数をとった定義式により確率密度を計算
log_C <- a * log(b) - lgamma(a)
log_dens <- log_C + (a - 1) * log(lambda) - b * lambda
dens <- exp(log_dens)
dens; log_dens
## [1] 0.1465251
## [1] -1.920558

 対数をとった定義式

$$ \begin{aligned} \log C_{\mathrm{Gam}} &= a \log b - \log \Gamma(a) \\ \log \mathrm{Gam}(\lambda | a, b) &= \log C_{\mathrm{Gam}} + (a - 1) \log \lambda - b \lambda \end{aligned} $$

を計算します。計算結果の指数をとると確率が得られます。

$$ \mathrm{Gam}(\lambda | a, b) = \exp \Bigr( \log \mathrm{Gam}(\lambda | a, b) \Bigr) $$

 指数と対数の性質より$\exp(\log x) = x$です。

 次は、関数を使って確率を計算します。
 ガンマ分布の確率密度関数dgamma()を使って計算します。

# ガンマ分布の関数により確率密度を計算
dens <- dgamma(x = lambda, shape = a, rate = b)
dens
## [1] 0.1465251

 変数の引数xlambda、形状の引数shapea、尺度の引数ratebを指定します。

 rate引数の代わりに、scale引数でも計算できます。

# ガンマ分布の関数により確率密度を計算
dens <- dgamma(x = lambda, shape = a, scale = 1 / b)
dens
## [1] 0.1465251

 尺度の引数scaleに$b$の逆数1 / bを指定します。

 log = TRUEを指定すると対数をとった確率密度を返します。

# ガンマ分布の対数をとった関数により確率密度を計算
log_dens <- dgamma(x = lambda, shape = a, rate = b, log = TRUE)
dens <- exp(log_dens)
dens; log_dens
## [1] 0.1465251
## [1] -1.920558

 計算結果の指数をとると確率密度が得られます。

 scaleの場合も同じです。

# ガンマ分布の対数をとった関数により確率密度を計算
log_dens <- dgamma(x = lambda, shape = a, scale = 1 / b, log = TRUE)
dens <- exp(log_dens)
dens; log_dens
## [1] 0.1465251
## [1] -1.920558


統計量の計算

 ガンマ分布の平均と分散、最頻値を計算します。

 平均を計算します。

# 平均を計算
E_lambda <- a / b
E_lambda
## [1] 1

 ガンマ分布の平均は、次の式で計算できます。

$$ \mathbb{E}[x] = \frac{a}{b} $$

 分散を計算します。

# 分散を計算
V_lambda <- a / b^2
V_lambda
## [1] 0.5

 ガンマ分布の分散は、次の式で計算できます。

$$ \mathbb{V}[x] = \frac{a}{b^2} $$

 最頻値を計算します。

# 最頻値を計算
mode_lambda <- (a - 1) / b
mode_lambda
## [1] 0.5

 ガンマ分布の最頻値は、次の式で計算できます。

$$ \mathrm{mode}[\lambda] = \frac{a - 1}{b} $$


分布の可視化

 ggplot2パッケージを利用してガンマ分布のグラフを作成します。

 ガンマ分布の確率変数がとり得る値$\lambda$ごとの確率密度を計算します。

# パラメータを指定
a <- 2
b <- 2

# 作図用のlambdaの点を作成
lambda_vals <- seq(from = 0, to = 5, length.out = 250)

# ガンマ分布を計算
dens_df <- tidyr::tibble(
  lambda = lambda_vals, # 確率変数
  density = dgamma(x = lambda_vals, shape = a, rate = b) # 確率密度
)
head(dens_df)
## # A tibble: 6 x 2
##   lambda density
##    <dbl>   <dbl>
## 1 0       0     
## 2 0.0201  0.0772
## 3 0.0402  0.148 
## 4 0.0602  0.214 
## 5 0.0803  0.274 
## 6 0.100   0.329

 $\lambda$がとり得る値を作成してlambda_valsとします。この例では、0から5を範囲とします。
 lambda_valsの各要素に対応する確率密度を求めてデータフレームに格納します。

 $\lambda > 0$なので、lambda0の行はdensity0になっています。

 ガンマ分布のグラフを作成します。

# ガンマ分布を作図
ggplot(data = dens_df, mapping = aes(x = lambda, y = density)) + # データ
  geom_line(color = "#00A968") + # 折れ線グラフ
  labs(title = "Gamma Distribution", 
       subtitle = paste0("a=", a, ", b=", b), 
       x = expression(lambda)) # ラベル

f:id:anemptyarchive:20220128155850p:plain
ガンマ分布のグラフ


 この分布に平均と最頻値、標準偏差の情報を重ねて表示します。

# 統計量の計算
E_lambda <- a / b
s_lambda <- sqrt(a / b^2)
mode_lambda <- (a - 1) / b

# 統計量を重ねたガンマ分布を作図
ggplot(data = dens_df, mapping = aes(x = lambda, y = density)) + # データ
  geom_line(color = "#00A968") + # 分布
  geom_vline(xintercept = E_lambda, color = "orange", size = 1, linetype = "dashed") + # 平均
  geom_vline(xintercept = E_lambda - s_lambda, color = "orange", size = 1, linetype = "dotted") + # 平均 - 標準偏差
  geom_vline(xintercept = E_lambda + s_lambda, color = "orange", size = 1, linetype = "dotted") + # 平均 + 標準偏差
  geom_vline(xintercept = mode_lambda, color = "chocolate", size = 1, linetype = "dashed") + # 最頻値
  labs(title = "Gamma Distribution", 
       subtitle = paste0("a=", a, ", b=", b), 
       x = expression(lambda)) # ラベル

f:id:anemptyarchive:20220128155903p:plain
ガンマ分布のグラフ

 ガンマ分布は対称な形ではないので、平均(オレンジ色の破線)と最頻値(茶色の破線)が一致しないのを確認できます。

 ガンマ分布のグラフを描画できました。

パラメータと分布の形状の関係

 パラメータが及ぼす分布への影響をアニメーション(gif画像)で可視化します。

 パラメータ$a, b$の値を少しずつ変更して、分布の変化をアニメーションで確認します。

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

# パラメータとして利用する値を指定
a_vals <- seq(from = 0.1, to = 10, by = 0.1)
b_vals <- seq(from = 0.1, to = 10, by = 0.1)

# 固定するパラメータを指定
a <- 2
b <- 2

# 作図用のlambdaの点を作成
lambda_vals <- seq(from = 0, to = 5, length.out = 250)

# a, bの値ごとに分布を計算
anime_dens_df <- tidyr::tibble()
for(a in a_vals) {
#for(b in b_vals) {
  # ガンマ分布を計算
  tmp_dens_df <- tidyr::tibble(
    lambda = lambda_vals, # 確率変数
    density = dgamma(x = lambda_vals, shape = a, rate = b), # 確率密度
    parameter = paste0("a=", a, ", b=", b) %>% 
      as.factor() # フレーム切替用のラベル
  )
  
  # 結果を結合
  anime_dens_df <- rbind(anime_dens_df, tmp_dens_df)
}
head(anime_dens_df)
## # A tibble: 6 x 3
##   lambda density parameter 
##    <dbl>   <dbl> <fct>     
## 1 0      Inf     a=0.1, b=2
## 2 0.0201   3.65  a=0.1, b=2
## 3 0.0402   1.88  a=0.1, b=2
## 4 0.0602   1.25  a=0.1, b=2
## 5 0.0803   0.928 a=0.1, b=2
## 6 0.100    0.729 a=0.1, b=2

 $a, b$がとり得る値を作成してa_vals, b_valsとします。
 for()ループを使ってa_valsまたはb_valsの値ごとに確率密度を計算して、anime_dens_dfに追加していきます。

 gif画像を作成します。

# アニメーション用のガウス分布を作図
anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = lambda, y = density)) + # データ
  geom_line(color = "#00A968") + # 折れ線グラフ
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Gamma Distribution", 
       subtitle = "{current_frame}", 
       x = expression(lambda)) # ラベル

# gif画像を作成
gganimate::animate(anime_dens_graph, nframes = length(a_vals), fps = 100)


f:id:anemptyarchive:20220128155919g:plainf:id:anemptyarchive:20220128155922g:plain
パラメータとガンマ分布の形状の関係

 平均$\mathbb{E}[\lambda] = \frac{a}{b}$、最頻値$\mathrm{mode}[\lambda] = \frac{a - 1}{b}$なので、$a$が大きくなるに従って$\lambda$が大きいほど確率密度が高くなり(山が右に移動し)ます。
 分散$\mathbb{V}[\lambda] = \frac{a}{b^2}$なので、$b$が大きくなるに従って分布が広がり(山が低くなり)ます。

乱数の生成

 ガンマ分布の乱数を生成してヒストグラムを確認します。

 パラメータを指定して、ガンマ分布に従う乱数を生成します。

# パラメータを指定
a <- 2
b <- 2

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

# ガンマ分布に従う乱数を生成
lambda_n <- rgamma(n = N, shape = a, rate = b)
lambda_n[1:5]
## [1] 0.8773600 0.6751416 2.9782231 0.2676596 0.1757919

 ガンマ分布の乱数生成関数rgamma()のデータ数(サンプルサイズ)の引数nN、パラメータの引数shape, ratea, bを指定します。

 作図に利用するため、サンプルをデータフレームに格納し、また$\lambda$の値と分布を作成しておきます。

# サンプルを格納
data_df <- tidyr::tibble(lambda = lambda_n)

# 作図用の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), 
                 breaks = seq(from = min(lambda_vals), to = max(lambda_vals), length.out = 30), 
                 fill = "#00A968", color = "#00A968") + # ヒストグラム
  labs(title = "Gamma Distribution", 
       subtitle = paste0("a=", a, ", b=", b), 
       x = expression(lambda)) # ラベル

f:id:anemptyarchive:20220128160008p:plain
ガンマ分布の乱数のヒストグラム:頻度

 geom_histogram()でヒストグラムを作成します。breaks引数に区切り位置を指定します。この例では、seq()を使ってlambda_valsの最小値から最大値の範囲をlength.out引数に指定した数に区切ります。

 サンプルの密度を分布と重ねて描画します。

# サンプルの密度を作図
ggplot() + 
  geom_histogram(data = data_df, mapping = aes(x = lambda, y = ..density..), 
                 breaks = seq(from = min(lambda_vals), to = max(lambda_vals), length.out = 30), 
                 fill = "#00A968", color = "#00A968") + # ヒストグラム
  geom_line(data = dens_df, mapping = aes(x = lambda, y = density), 
            color = "darkgreen", linetype = "dashed") + # 元の分布
  labs(title = "Gamma Distribution", 
       subtitle = paste0("a=", a, ", b=", b), 
       x = expression(lambda)) # ラベル

f:id:anemptyarchive:20220128160025p:plain
ガンマ分布の乱数のヒストグラム:密度

 geom_histogram()y = ..density..を指定すると、頻度を密度に変換して描画します。

 データ数が十分に増えると元の分布(破線のグラフ)に形が近付きます。

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

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

# フレーム数を指定
N_frame <- 300

# 乱数を1つずつ生成
lambda_n <- rep(NA, times = N_frame)
anime_freq_df <- tidyr::tibble()
anime_data_df <- tidyr::tibble()
anime_dens_df <- tidyr::tibble()
for(n in 1:N_frame) {
  # ガンマ分布に従う乱数を生成
  lambda_n[n] <- rgamma(n = 1, shape = a, rate = b)
  
  # ラベル用のテキストを作成
  label_text <- paste0("a=", a, ", b=", b, ", N=", n)
  
  # n個の乱数を格納
  tmp_freq_df <- tidyr::tibble(
    lambda = lambda_n[1:n], # サンプル
    parameter = as.factor(label_text) # フレーム切替用のラベル
  )
  
  # n番目の乱数を格納
  tmp_data_df <- tidyr::tibble(
    lambda = lambda_n[n], # サンプル
    parameter = as.factor(label_text) # フレーム切替用のラベル
  )
  
  # n回目のラベルを付与
  tmp_dens_df <- dens_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_dens_df <- rbind(anime_dens_df, tmp_dens_df)
}
head(anime_freq_df)
## # A tibble: 6 x 2
##   lambda parameter    
##    <dbl> <fct>        
## 1 0.410  a=2, b=2, N=1
## 2 0.410  a=2, b=2, N=2
## 3 0.0549 a=2, b=2, N=2
## 4 0.410  a=2, b=2, N=3
## 5 0.0549 a=2, b=2, N=3
## 6 2.03   a=2, b=2, N=3

 乱数を1つずつ生成して、結果をanime_***_dfに追加していきます。

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

# アニメーション用のサンプルの頻度を作図
anime_freq_graph <- ggplot() + 
  geom_histogram(data = anime_freq_df, mapping = aes(x = lambda), 
                 breaks = seq(from = min(lambda_vals), to = max(lambda_vals), length.out = 30), 
                 fill = "#00A968", color = "#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}") # ラベル

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


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

# アニメーション用のサンプルの密度を作図
anime_prop_graph <- ggplot() + 
  geom_histogram(data = anime_freq_df, mapping = aes(x = lambda, y = ..density..), 
                 breaks = seq(from = min(lambda_vals), to = max(lambda_vals), length.out = 30), 
                 fill = "#00A968", color = "#00A968") + # ヒストグラム
  geom_line(data = anime_dens_df, mapping = aes(x = lambda, y = density), 
            color = "darkgreen", linetype = "dashed") + # 元の分布
  geom_point(data = anime_data_df, mapping = aes(x = lambda, y = 0), 
             color = "orange", size = 5) + # サンプル
  gganimate::transition_manual(parameter) + # フレーム
  ylim(c(-0.01, max(dens_df[["density"]]) * 2)) + # y軸の表示範囲
  labs(title = "Gamma Distribution", 
       subtitle = "{current_frame}") # ラベル

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


f:id:anemptyarchive:20220128161350g:plainf:id:anemptyarchive:20220128161354g:plain
ガンマ分布の乱数の推移

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

分布の生成

 ガンマ分布が共役事前分布となる1次元ガウス分布とポアソン分布のパラメータを生成して分布を作図します。

 パラメータ$\lambda$を生成します。

# パラメータを指定
a <- 5
b <- 2

# サンプルサイズを指定
N <- 10

# ガウス分布・ポアソン分布のパラメータを生成
lambda_n <- rgamma(n = N, shape = a, rate = b)
lambda_n[1:5]
## [1] 1.386198 2.732954 2.000720 2.198415 1.714121

 ガンマ分布に従う乱数を生成して、パラメータ$\lambda$として利用します。

1次元ガウス分布

 生成した$\lambda$を1次元ガウス分布の精度パラメータとして利用します。1次元ガウス分布の計算と可視化については「(2日後に更新予定の)1次元ガウス分布の作図」を参照してください。

 1次元ガウス分布は、平均パラメータ$\mu$、精度パラメータ$\lambda$を用いて次の式で定義されます。

$$ \mathcal{N}(x | \mu, \lambda^{-1}) = \frac{1}{\sqrt{2 \pi \lambda^{-1}}} \exp \left( - \frac{\lambda}{2} (x - \mu)^2 \right) $$

 ここで、精度$\lambda$は分散$\sigma^2$の逆数$\lambda = \frac{1}{\sigma^2}$です。

 まずは、目安となるように$\lambda$の期待値$\mathbb{E}[\lambda]$による分布を求めます。

# 平均パラメータを指定
mu = 0.0

# 精度パラメータの期待値を計算
E_lambda <- a / b

# 標準偏差の期待値を計算
E_sigma <- sqrt(1.0 / E_lambda)

# 作図用のxの点を作成
x_vals <- seq(from = mu - E_sigma*4, to = mu + E_sigma*4, length.out = 250)

# 精度パラメータの期待値による1次元ガウス分布を計算
res_dens_df <- tidyr::tibble(
  x = x_vals, # 確率変数
  density = dnorm(x = x_vals, mean = mu, sd = E_sigma), # 確率密度
  parameter = paste0("E[lambda]=", round(E_lambda, 2)) %>% 
    as.factor() # 色分け用のラベル
)
head(res_dens_df)
## # A tibble: 6 x 3
##       x  density parameter    
##   <dbl>    <dbl> <fct>        
## 1 -2.53 0.000212 E[lambda]=2.5
## 2 -2.51 0.000240 E[lambda]=2.5
## 3 -2.49 0.000273 E[lambda]=2.5
## 4 -2.47 0.000310 E[lambda]=2.5
## 5 -2.45 0.000351 E[lambda]=2.5
## 6 -2.43 0.000397 E[lambda]=2.5

 1次元ガウス分布の平均パラメータ$\mu$を指定します。

 精度パラメータの期待値$\mathbb{E}[\lambda] = \frac{a}{b}$を計算して、E_lambdaとします。精度$\lambda$の逆数$\frac{1}{\lambda}$が分散$\sigma^2$であり、分散$\sigma^2$の平方根$\sqrt{\sigma^2}$が標準偏差$\sigma$になります。
 標準偏差の期待値$\mathbb{E}[\sigma] = \sqrt{\frac{1}{\mathbb{E}[\lambda]}}$を計算して、E_sigmaとします。

 $x$として利用する値を作成してx_valsとします。この例では、平均muを中心に標準偏差のE_sigma4倍を範囲とします。
 dnorm()で1次元ガウス分布の確率密度を計算します。平均の引数meanmu、標準偏差の引数sdE_sigmaを指定します。

 次に、パラメータのサンプルlambda_nの値ごとに分布を計算します。

# lambdaの値ごとに分布を計算
for(lambda in lambda_n) {
  # 1次元ガウス分布を計算
  tmp_dens_df <- tidyr::tibble(
    x = x_vals, # 確率変数
    density = dnorm(x = x_vals, mean = mu, sd = sqrt(1 / lambda)), # 確率密度
    parameter = paste0("lambda=", round(lambda, 2)) %>% 
      as.factor() # 色分け用のラベル
  )
  
  # 結果を結合
  res_dens_df <- rbind(res_dens_df, tmp_dens_df)
}

 各分布の情報をres_dens_dfに追加します。

 $N + 1$個の1次元ガウス分布を作図します。

# サンプルによる分布を作図
ggplot() + 
  geom_line(data = res_dens_df, 
            mapping = aes(x = x, y = density, color = parameter, 
                          alpha = parameter, linetype = parameter), 
            size = 0.8) + # 分布
  scale_linetype_manual(values = c("dashed", rep("solid", times = N))) + # 線の種類
  scale_alpha_manual(values = c(1, rep(0.5, times = N))) + # 透過度
  labs(title = "Gaussian Distribution", 
       subtitle = paste0("a=", a, ", b=", b), 
       x = expression(lambda)) # ラベル

f:id:anemptyarchive:20220128160156p:plain
ガンマ分布からサンプリングした1次元ガウス分布のグラフ

 平均の分布(破線)を中心に分布しています。

ポアソン分布

 続いて、生成した$\lambda$をポアソン分布のパラメータとして利用します。ポアソン分布の計算と可視化については「【R】ポアソン分布の作図 - からっぽのしょこ」を参照してください。

 ポアソン分布は、パラメータ$\lambda$を用いて次の式で定義されます。

$$ \mathrm{Poi}(x | \lambda) = \frac{\lambda^x}{x!} e^{-\lambda} $$

 まずは、目安となるように$\lambda$の期待値$\mathbb{E}[\lambda]$による分布を求めます。

# パラメータの期待値を計算
E_lambda <- a / b

# 作図用のxの点を作成
x_vals <- seq(from = 0, to = ceiling(E_lambda) * 4)

# パラメータの期待値によるガンマ分布を計算
res_prob_df <- tidyr::tibble(
  x = x_vals, # 確率変数
  probability = dpois(x = x_vals, lambda = E_lambda), # 確率
  parameter = paste0("E[lambda]=", round(E_lambda, 2)) %>% 
    as.factor() # 色分け用のラベル
)
head(res_prob_df)
## # A tibble: 6 x 3
##       x probability parameter    
##   <int>       <dbl> <fct>        
## 1     0      0.0821 E[lambda]=2.5
## 2     1      0.205  E[lambda]=2.5
## 3     2      0.257  E[lambda]=2.5
## 4     3      0.214  E[lambda]=2.5
## 5     4      0.134  E[lambda]=2.5
## 6     5      0.0668 E[lambda]=2.5

 こちらも、パラメータの期待値$\mathbb{E}[\lambda] = \frac{a}{b}$を計算して、E_lambdaとします。

 $x$として利用する値を作成してx_valsとします。この例では、0からパラメータE_lambdaの最大値の4倍を範囲とする整数とします。整数となるように、E_lambdaceiling()で切り上げて利用します。
 dpis()でポアソン分布の確率を計算します。変数の引数xx_vals、パラメータの引数lambdaE_lambdaを指定します。

 次に、パラメータのサンプルlambda_nの値ごとにポアソン分布を計算します。

# lambdaの値ごとに分布を計算
for(lambda in lambda_n) {
  # ガンマ分布を計算
  tmp_prob_df <- tidyr::tibble(
    x = x_vals, # 確率変数
    probability = dpois(x = x_vals, lambda = lambda), # 確率
    parameter = paste0("lambda=", round(lambda, 2)) %>% 
      as.factor() # 色分け用のラベル
  )
  
  # 結果を結合
  res_prob_df <- rbind(res_prob_df, tmp_prob_df)
}

 各分布の情報をres_prob_dfに追加します。

 $N + 1$個のポアソン分布を作図します。

# サンプルによる分布を作図
ggplot() + 
  geom_step(data = res_prob_df, 
            mapping = aes(x = x, y = probability, color = parameter, 
                          alpha = parameter, linetype = parameter), 
            direction = "mid", size = 0.8) + # 分布
  #geom_bar(data = res_prob_df, mapping = aes(x = x, y = probability, color = parameter, linetype = parameter), 
  #         stat = "identity", alpha = 0) + # 分布
  scale_linetype_manual(values = c("dashed", rep("solid", times = N))) + # 線の種類
  scale_alpha_manual(values = c(1, rep(0.5, times = N))) + # 透過度
  scale_x_continuous(breaks = x_vals) + # x軸目盛
  labs(title = "Poisson Distribution", 
       subtitle = paste0("a=", a, ", b=", b), 
       x = expression(lambda)) # ラベル

f:id:anemptyarchive:20220128160225p:plain
ガンマ分布からサンプリングしたポアソン分布のグラフ

 複数の棒グラフを重ねると分かりにくいので、geom_step()で階段状のグラフとして描画します。

 平均の分布(破線)を中心に分布しています。

参考文献

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

おわりに

 統計量の導出もやりたいけどしばらくお預けかな。

 1月28日はモーニング娘。のメジャーデビュー日です!

 24周年って凄すぎる。私がハマってからは4年半くらい、これからも楽しい日々が続きますよーに。

【次の内容】

www.anarchive-beta.com