からっぽのしょこ

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

【R】ガウス-ガンマ分布の作図

はじめに

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

 ガウス-ガンマ分布(正規ガンマ分布)の計算とグラフの作成をR言語で行います。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

ガウス-ガンマ分布の作図

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

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

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

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

定義式の確認

 まずは、ガウス-ガンマ分布の定義式を確認します。1次元ガウス分布については「【R】1次元ガウス分布の作図 - からっぽのしょこ」、ガンマ分布については「【R】ガンマ分布の作図 - からっぽのしょこ」を参照してください。

 ガウス-ガンマ分布は、1次元ガウス分布とガンマ分布の積で定義される$\mu, \lambda$の同時分布です。

$$ \begin{aligned} \mathrm{NG}(\mu, \lambda | m, \beta, a, b) &= \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) \mathrm{Gam}(\lambda | a, b) \\ &= \frac{1}{\sqrt{2 \pi (\beta \lambda)^{-1}}} \exp \left( - \frac{\beta \lambda}{2} (\mu - m)^2 \right) \frac{b^a}{\Gamma(a)} \lambda^{a-1} \exp(- b \lambda) \end{aligned} $$

 ここで、$\mathcal{N}(\cdot)$はガウス分布、$\mathrm{Gam}(\cdot)$はガンマ分布を表します。$m$は1次元ガウス分布の平均パラメータ、$\beta \lambda$は精度パラメータです。精度の逆数$\frac{1}{\beta \lambda}$が分散パラメータで、その平方根$\sqrt{\frac{1}{\beta \lambda}}$が標準偏差になります。また、$\lambda$はガンマ分布の変数でもあり、$a, b$はガンマ分布のパラメータです。
 確率変数の値$\mu$は実数、$\lambda$は0より大きい値$\lambda > 0$となります。パラメータ$\lambda, \beta, a, b$は、0より大きい値$\lambda > 0, \beta > 0, a > 0, b > 0$を満たす必要があります。

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

$$ \begin{aligned} \log \mathrm{NG}(\mu, \lambda | m, \beta, a, b) &= \log \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) + \log \mathrm{Gam}(\lambda | a, b) \\ &= - \frac{1}{2} \Bigl\{ \beta \lambda (x - m)^2 - \log(\beta \lambda) + \log(2 \pi) \Bigr\} \\ &\quad + (a - 1) \log \lambda - b \lambda + a \log b - \log \Gamma(a) \end{aligned} $$

 ガウス-ガンマ分布の$\mu$の平均と分散は、定義より次となります。

$$ \begin{aligned} \mathbb{E}[\mu] &= m \\ \mathbb{V}[\mu] &= (\beta \lambda)^{-1} = \frac{1}{\beta \lambda} \end{aligned} $$

 また、$\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} $$


 ガウス-ガンマ分布は、1次元ガウス分布の事前分布として利用されます。

確率密度の計算

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

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

# 1次元ガウス分布のパラメータを指定
m <- 0

# 1次元ガウス分布の精度パラメータの係数を指定
beta <- 2

# ガンマ分布のパラメータを指定
a <- 5
b <- 6

# 確率変数の値を指定
mu     <- 1.5
lambda <- 2.5

 1次元ガウス分布の平均パラメータ$m$、精度パラメータの係数$\beta > 0$、ガンマ分布のパラメータ$a > 0, b > 0$、確率変数がとり得る値$\mu$、$\lambda > 0$を指定します。設定した値に対する確率密度を計算します。

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

# 定義式により確率密度を計算
C_N      <- 1 / sqrt(2 * pi / beta / lambda)
dens_N   <- C_N * exp(-0.5 * beta * lambda * (mu - m)^2)
C_Gam    <- b^a / gamma(a)
dens_Gam <- C_Gam * lambda^(a - 1) * exp(-b * lambda)
dens     <- dens_N * dens_Gam
dens
## [1] 1.245594e-05

 1次元ガウス分布

$$ \begin{aligned} C_{\mathcal{N}} &= \frac{1}{\sqrt{2 \pi (\beta \lambda)^{-1}}} \\ \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) &= C_{\mathcal{N}} \exp \left( - \frac{\beta \lambda}{2} (\mu - m)^2 \right) \end{aligned} $$

と、ガンマ分布

$$ \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} $$

をそれぞれ計算して、2つの積

$$ \begin{aligned} \mathrm{NG}(\mu, \lambda | m, \beta, a, b) &= \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) \mathrm{Gam}(\lambda | a, b) \end{aligned} $$

を求めます。$C_{\mathcal{N}}$は1次元ガウス分布の正規化係数、$C_{\mathrm{Gam}}$は、ガンマ分布の正規化係数です。また、$\pi$は円周率です。
 $\Gamma(x)$はガンマ関数で、gamma()で計算できます。

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

# 対数をとった定義式により確率密度を計算
log_C_N      <- -0.5 * (log(2 * pi) - log(beta * lambda))
log_dens_N   <- log_C_N - 0.5 * beta * lambda * (mu - m)^2
log_C_Gam    <- a * log(b) - lgamma(a)
log_dens_Gam <- log_C_Gam + (a - 1) * log(lambda) - b * lambda
log_dens     <- log_dens_N + log_dens_Gam
dens         <- exp(log_dens)
dens; log_dens
## [1] 1.245594e-05
## [1] -11.29331

 対数をとった1次元ガウス分布

$$ \begin{aligned} \log C_{\mathcal{N}} &= - \frac{1}{2} \Bigl\{ \log (2 \pi) - \log (\beta \lambda) \Bigr\} \\ \log \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) &= \log C_{\mathcal{N}} - \frac{\beta \lambda}{2} (\mu - m)^2 \end{aligned} $$

と、対数をとったガンマ分布

$$ \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} $$

の和

$$ \log \mathrm{NG}(\mu, \lambda | m, \beta, a, b) = \log \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) + \log \mathrm{Gam}(\lambda | a, b) $$

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

$$ \mathrm{NG}(\mu, \lambda | m, \beta, a, b) = \exp \Bigr( \log \mathrm{NG}(\mu, \lambda | m, \beta, a, b) \Bigr) $$

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

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

# 関数により確率密度を計算
dens_N   <- dnorm(x = mu, mean = m, sd = sqrt(1 / beta / lambda))
dens_Gam <- dgamma(x = lambda, shape = a, rate = b)
dens     <- dens_N * dens_Gam
dens
## [1] 1.245594e-05

 dnorm()の変数の引数xmu、平均の引数meanm、標準偏差の引数sdに精度beta * lambdaの逆数(分散)1 / beta / lambdaの平方根(標準偏差)を指定します。平方根はsqrt()で計算できます。
 dgamma()の変数の引数xlambda、形状の引数shapea、尺度の引数ratebを指定します。

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

# 対数をとった関数により確率密度を計算
log_dens_N   <- dnorm(x = mu, mean = m, sd = sqrt(1 / beta / lambda), log = TRUE)
log_dens_Gam <- dgamma(x = lambda, shape = a, rate = b, log = TRUE)
log_dens     <- log_dens_N + log_dens_Gam
dens         <- exp(log_dens)
dens; log_dens
## [1] 1.245594e-05
## [1] -11.29331

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

統計量の計算

 ガウス-ガンマ分布の平均と分散を計算します。

 平均を計算します。

# 平均を計算
E_mu <- m
E_lambda <- a / b
E_mu; E_lambda
## [1] 0
## [1] 0.8333333

 1次元ガウス分布($\mu$)とガンマ分布($\lambda$)の平均は、次の式で計算できます。

$$ \begin{aligned} \mathbb{E}[\mu] &= m \\ \mathbb{E}[\lambda] &= \frac{a}{b} \end{aligned} $$

 分散を計算します。

# 分散を計算
V_mu <- 1 / beta / E_lambda
V_lambda <- a / b^2
V_mu; V_lambda
## [1] 0.6
## [1] 0.1388889

 1次元ガウス分布($\mu$)とガンマ分布($\lambda$)の分散は、次の式で計算できます。

$$ \begin{aligned} \mathbb{V}[\mu] &= \frac{1}{\beta \lambda} \\ \mathbb{V}[\lambda] &= \frac{a}{b^2} \end{aligned} $$

 最頻値を計算します。

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

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

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


分布の可視化

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

 作図に利用する$\mu, \lambda$の値を設定します。

# 1次元ガウス分布の平均パラメータを指定
m <- 0

# 1次元ガウス分布の精度パラメータの係数を指定
beta <- 2

# ガンマ分布のパラメータを指定
a <- 5
b <- 6

# 作図用のmuの値を作成
mu_vals <- seq(from = -3, to = 3, length.out = 201)

# 作図用のlambdaの値を作成
lambda_vals <- seq(from = 0.01, to = 2, length.out = 200)
mu_vals[1:5]; lambda_vals[1:5]
## [1] -3.00 -2.97 -2.94 -2.91 -2.88
## [1] 0.01 0.02 0.03 0.04 0.05

 $\mu$がとり得る値を作成してmu_valsとします。この例では、-3から3を範囲とします。
 $\lambda$がとり得る値を作成してlambda_valsとします。この例では、0.01から2を範囲とします。

 作図用の点$(\mu, \lambda)$を作成します。

# 作図用のmuとlambdaの点を作成
mu_lambda_points <- expand.grid(mu = mu_vals, lambda = lambda_vals) %>% 
  as.matrix()
mu_lambda_points <- tidyr::tibble(
  mu = rep(mu_vals, times = length(lambda_vals)), # muの値
  lambda = rep(lambda_vals, each = length(mu_vals)) # lambdaの値
) %>% 
  as.matrix()
mu_lambda_points[1:5, ]
##         mu lambda
## [1,] -3.00   0.01
## [2,] -2.97   0.01
## [3,] -2.94   0.01
## [4,] -2.91   0.01
## [5,] -2.88   0.01

 expand.grid()またはrep()を使って、mu_vals, lambda_valsそれぞれの値に対して、全ての組み合わせを持つ(格子状の点となる)ように値を複製します。
 各行が1つの点$(\mu, \lambda)$に対応します。

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

# ガウス-ガンマ分布を計算
dens_df <- tidyr::tibble(
  mu = mu_lambda_points[, 1], # 確率変数mu
  lambda = mu_lambda_points[, 2], # 確率変数lambda
  N_dens = dnorm(x = mu, mean = m, sd = sqrt(1 / beta / lambda)), # muの確率密度
  Gam_dens = dgamma(x = lambda, shape = a, rate = b), # lambdaの確率密度
  density = N_dens * Gam_dens # 確率密度
)
head(dens_df)
## # A tibble: 6 x 5
##      mu lambda N_dens   Gam_dens     density
##   <dbl>  <dbl>  <dbl>      <dbl>       <dbl>
## 1 -3      0.01 0.0516 0.00000305 0.000000157
## 2 -2.97   0.01 0.0517 0.00000305 0.000000158
## 3 -2.94   0.01 0.0517 0.00000305 0.000000158
## 4 -2.91   0.01 0.0518 0.00000305 0.000000158
## 5 -2.88   0.01 0.0519 0.00000305 0.000000158
## 6 -2.85   0.01 0.0520 0.00000305 0.000000159

 mu_lambda_pointsの各行に対応する確率密度を求めてデータフレームに格納します。

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

# ガウス-ガンマ分布を作図
ggplot() + 
  geom_contour(data = dens_df, aes(x = mu, y = lambda, z = density, color = ..level..)) + # 等高線図
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = paste0("m=", m, ", beta=", beta, ", a=", a, ", b=", b), 
       x = expression(mu), y = expression(lambda), 
       color = "density") # ラベル

f:id:anemptyarchive:20220203222544p:plain
ガウス-ガンマ分布のグラフ


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

# 統計量を計算
E_lambda    <- a / b
s_lambda    <- sqrt(a / b^2)
mode_lambda <- (a - 1) / b
E_mu    <- m
E_s_mu  <- sqrt(1 / beta / E_lambda)
s_mu_df <- tidyr::tibble(
  lambda = lambda_vals, 
  s_minus = E_mu - sqrt(1 / beta / lambda_vals), # 平均 - 標準偏差
  s_plus = E_mu + sqrt(1 / beta / lambda_vals)   # 平均 + 標準偏差
)

# 統計量を重ねたガウス-ガンマ分布を作図
ggplot() + 
  geom_contour(data = dens_df, aes(x = mu, y = lambda, z = density, color = ..level..)) + # 分布
  geom_vline(xintercept = m, color = "#00A968", size = 1.5, linetype = "dashed") + # 平均
  geom_line(data = s_mu_df, mapping = aes(x = s_minus, y = lambda), color = "#00A968", size = 1.5, linetype = "dotted") + # 平均 - 標準偏差
  geom_line(data = s_mu_df, mapping = aes(x = s_plus, y = lambda), color = "#00A968", size = 1.5, linetype = "dotted") + # 平均 + 標準偏差
  geom_hline(yintercept = E_lambda, color = "orange", size = 1.5, linetype = "dashed") + # 平均
  geom_hline(yintercept = E_lambda - s_lambda, color = "orange", size = 1.5, linetype = "dotted") + # 平均 - 標準偏差
  geom_hline(yintercept = E_lambda + s_lambda, color = "orange", size = 1.5, linetype = "dotted") + # 平均 + 標準偏差
  geom_hline(yintercept = mode_lambda, color = "chocolate", size = 1.5, linetype = "dashed") + # 最頻値
  xlim(c(min(mu_vals), max(mu_vals))) + # x軸の表示範囲
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = paste0("m=", m, ", beta=", beta, ", a=", a, ", b=", b), 
       x = expression(mu), y = expression(lambda), 
       color = "density")

f:id:anemptyarchive:20220203222602p:plain
ガウス-ガンマ分布のグラフ

 ガンマ分布は、$\mu$の影響を受けない(定義式に$\mu$を含まない)ので、平均・標準偏差・最頻値は一定(水平線)になります。また、ガンマ分布は対称な形ではないので、平均(オレンジ色の破線)と最頻値(茶色の破線)が一致しません。オレンジ色の点線は、平均から標準偏差を足し引きした値です。
 1次元ガウス分布の標準偏差は、$\lambda$の影響を受けます(計算式に$\lambda$を含みます)。よって、平均から標準偏差を足し引きした値(青緑色の点線)は$\lambda$(y軸)の値によって変化します。この曲線の影響を受けて、ガウス-ガンマ分布がおにぎり型になるのを確認できます。

 ガウス-ガンマ分布のグラフを描画できました。続いて、1つの変数(各軸)$\mu, \lambda$に注目します。

 $\lambda = \mathbb{E}[\lambda]$における1次元ガウス分布のグラフを作成します。

# 精度の期待値による1次元ガウス分布を計算
dens_N_df <- tidyr::tibble(
  mu = mu_vals, # 確率変数
  density = dnorm(x = mu_vals, mean = m, sd = sqrt(1 / beta / E_lambda)) # 確率密度
)

# 統計量を重ねた精度の期待値による1次元ガウス分布を作図
ggplot() + 
  geom_line(data = dens_N_df, mapping = aes(x = mu, y = density), color = "#00A968") + # 分布
  geom_vline(xintercept = E_mu, color = "#00A968", size = 1.5, linetype = "dashed") + # 平均
  geom_vline(xintercept = E_mu - E_s_mu, color = "#00A968", size = 1.5, linetype = "dotted") + # 平均 - 標準偏差
  geom_vline(xintercept = E_mu + E_s_mu, color = "#00A968", size = 1.5, linetype = "dotted") + # 平均 + 標準偏差
  labs(title = "Gaussian Distribution", 
       subtitle = paste0("m=", m, ", beta=", beta, ", E[lambda]=", round(E_lambda, 2)), 
       x = expression(mu), y = expression(p(mu)))

f:id:anemptyarchive:20220203222622p:plain
x軸側から見たグラフ

 この図は、ガウス-ガンマ分布を$\lambda$軸(y軸)が$\mathbb{E}[\lambda]$の点で切断した、断面図に対応します。

 同様に、全ての$\lambda$における分布をアニメーションで確認します。

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

# lambdaの値ごとに分布を計算
anime_dens_N_df <- tidyr::tibble()
anime_s_mu_df    <- tidyr::tibble()
for(lambda in lambda_vals) {
  # ラベル用のテキストを作成
  label_text <- paste0("m=", m, ", beta=", beta, ", lambda=", round(lambda, 2))
  
  # 1次元ガウス分布を計算
  tmp_dens_N_df <- tidyr::tibble(
    mu = mu_vals, # 確率変数
    density = dnorm(x = mu_vals, mean = m, sd = sqrt(1 / beta / lambda)), # 確率密度
    parameter = as.factor(label_text) # フレーム切替用のラベル
  )
  
  # 標準偏差を格納
  tmp_s_mu_df <- tidyr::tibble(
    s_minus = m - sqrt(1 / beta / lambda), # 平均 - 標準偏差
    s_plus = m + sqrt(1 / beta / lambda),  # 平均 + 標準偏差
    parameter = as.factor(label_text)      # フレーム切替用のラベル
  )
  
  # 結果を結合
  anime_dens_N_df <- rbind(anime_dens_N_df, tmp_dens_N_df)
  anime_s_mu_df   <- rbind(anime_s_mu_df, tmp_s_mu_df)
}
head(anime_dens_N_df); head(anime_s_mu_df)
## # A tibble: 6 x 3
##      mu density parameter               
##   <dbl>   <dbl> <fct>                   
## 1 -3     0.0516 m=0, beta=2, lambda=0.01
## 2 -2.97  0.0517 m=0, beta=2, lambda=0.01
## 3 -2.94  0.0517 m=0, beta=2, lambda=0.01
## 4 -2.91  0.0518 m=0, beta=2, lambda=0.01
## 5 -2.88  0.0519 m=0, beta=2, lambda=0.01
## 6 -2.85  0.0520 m=0, beta=2, lambda=0.01
## # A tibble: 6 x 3
##   s_minus s_plus parameter               
##     <dbl>  <dbl> <fct>                   
## 1   -7.07   7.07 m=0, beta=2, lambda=0.01
## 2   -5      5    m=0, beta=2, lambda=0.02
## 3   -4.08   4.08 m=0, beta=2, lambda=0.03
## 4   -3.54   3.54 m=0, beta=2, lambda=0.04
## 5   -3.16   3.16 m=0, beta=2, lambda=0.05
## 6   -2.89   2.89 m=0, beta=2, lambda=0.06

 for()ループを使ってlambda_valsの値ごとに確率密度を計算して、anime_dens_dfに追加していきます。

 gif画像を作成します。

# アニメーション用の統計量を重ねた1次元ガウス分布を作図
anime_dens_N_graph <- ggplot() + 
  geom_line(data = anime_dens_N_df, mapping = aes(x = mu, y = density), color = "#00A968") + # 分布
  geom_vline(xintercept = E_mu, color = "#00A968", size = 1.5, linetype = "dashed") + # 平均
  geom_vline(data = anime_s_mu_df, mapping = aes(xintercept = s_minus), color = "#00A968", size = 1.5, linetype = "dotted") + # 平均 - 標準偏差
  geom_vline(data = anime_s_mu_df, mapping = aes(xintercept = s_plus), color = "#00A968", size = 1.5, linetype = "dotted") + # 平均 + 標準偏差
  gganimate::transition_manual(parameter) + # フレーム
  xlim(c(min(mu_vals), max(mu_vals))) + # x軸の表示範囲
  labs(title = "Gaussian Distribution", 
       subtitle = "{current_frame}", 
       x = expression(mu), y = expression(p(mu)))

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


f:id:anemptyarchive:20220203222636g:plain
x軸側から見たグラフ

 $\lambda$が大きくなるに従って、1次元ガウス分布の裾が狭くなり(山が細く高くなり)ます。これが、ガウス-ガンマ分布がおにぎり型になるのに対応します。
 $\mu \pm \sigma$の点線が移動するのが、ガウス-ガンマ分布の$\mu \pm \sigma$の点線が曲線になるのに対応しています。

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

# ガンマ分布を計算
dens_Gam_df <- tidyr::tibble(
  lambda = lambda_vals, # 確率変数
  density = dgamma(x = lambda_vals, shape = a, rate = b) # 確率密度
)

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

f:id:anemptyarchive:20220203222732p:plain
y軸側から見たグラフ

 ガウス-ガンマ分布をy軸側から見た図に対応します。
 ガンマ分布は対称な形ではないので、平均(オレンジ色の破線)と最頻値(茶色の破線)が一致しないのを確認できます。

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

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

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

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

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

# 固定するパラメータを指定
m    <- 0
beta <- 2
a    <- 5
b    <- 6

# 作図用の変数の値を作成
mu_vals     <- seq(from = -3, to = 3, length.out = 201)
lambda_vals <- seq(from = 0.01, to = 2, length.out = 200)

# フレーム数を設定
frame_num <- length(m_vals)
#frame_num <- length(beta_vals)
#frame_num <- length(a_vals)
#frame_num <- length(b_vals)

# パラメータの値ごとに分布を計算
anime_dens_df <- tidyr::tibble()
for(i in 1:frame_num) {
  # i番目の値を取得
  m    <- m_vals[i]
  #beta <- beta_vals[i]
  #a    <- a_vals[i]
  #b    <- b_vals[i]
  
  # ガウス-ガンマ分布を計算
  tmp_dens_df <- tidyr::tibble(
    mu = rep(mu_vals, times = length(lambda_vals)), # 確率変数mu
    lambda = rep(lambda_vals, each = length(mu_vals)), # 確率変数lambda
    N_dens = dnorm(x = mu, mean = m, sd = sqrt(1 / beta / lambda)), # muの確率密度
    Gam_dens = dgamma(x = lambda, shape = a, rate = b), # lambdaの確率密度
    density = N_dens * Gam_dens, # 確率密度
    parameter = paste0(
      "m=", round(m, 1), ", beta=", round(beta, 1), ", a=", round(a, 1), ", b=", round(b, 1)
    ) %>% 
      as.factor() # フレーム切替用のラベル
  )
  
  # 結果を結合
  anime_dens_df <- rbind(anime_dens_df, tmp_dens_df)
  
  # 途中経過を表示
  #message("\r", appendLF = FALSE) # 前回の表示を消去
  #message("\r", "i=", i, " (", round(i / frame_num * 100, 2), "%)", appendLF = FALSE)
}
head(anime_dens_df)
## # A tibble: 6 x 6
##      mu lambda N_dens   Gam_dens     density parameter             
##   <dbl>  <dbl>  <dbl>      <dbl>       <dbl> <fct>                 
## 1 -3      0.01 0.0559 0.00000305 0.000000170 m=-2, beta=2, a=5, b=6
## 2 -2.97   0.01 0.0559 0.00000305 0.000000171 m=-2, beta=2, a=5, b=6
## 3 -2.94   0.01 0.0559 0.00000305 0.000000171 m=-2, beta=2, a=5, b=6
## 4 -2.91   0.01 0.0560 0.00000305 0.000000171 m=-2, beta=2, a=5, b=6
## 5 -2.88   0.01 0.0560 0.00000305 0.000000171 m=-2, beta=2, a=5, b=6
## 6 -2.85   0.01 0.0560 0.00000305 0.000000171 m=-2, beta=2, a=5, b=6

 各パラメータがとり得る値を作成して*_valsとします。
 for()ループを使って、注目するパラメータ*_valsから順番に値を取り出して確率密度を計算して、anime_dens_dfに追加していきます。

 gif画像を作成します。

# アニメーション用のガウス-ガンマ分布を作図
anime_dens_graph <- ggplot() + 
  geom_contour(data = anime_dens_df, aes(x = mu, y = lambda, z = density, color = ..level..)) + # 分布
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = "{current_frame}", 
       x = expression(mu), y = expression(lambda), 
       color = "density")

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


f:id:anemptyarchive:20220203222800g:plain
$m$とガウス-ガンマ分布の形状の関係

 $m$は$\mu$の平均なので、$m$が大きくなるに従って、$\mu$(x軸)の値が大きいほど確率密度が大きくなり(山が右に移動し)ます。

f:id:anemptyarchive:20220203222936g:plain
$\beta$とガウス-ガンマ分布の形状の関係

 $\beta$は$\mu$の精度(分散)に影響するので、$\beta$が大きくなるに従って、$\mu$(x軸)における確率密度が大きい範囲が狭くなり(山が尖り)ます。

f:id:anemptyarchive:20220203223014g:plain
$a$とガウス-ガンマ分布の形状の関係

 $a, b$は$\lambda$に対して影響します。$a$が大きくなるに従って、$\lambda$(y軸)が大きいほど確率密度が大きくなります。

f:id:anemptyarchive:20220203223045g:plain
$b$とガウス-ガンマ分布の形状の関係

 $b$が大きくなるに従って、$\lambda$(y軸)が小さい(0に近い)ほど確率密度が大きくなります。

乱数の生成

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

サンプリング

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

# 1次元ガウス分布の平均パラメータを指定
m <- 0

# 1次元ガウス分布の精度パラメータの係数を指定
beta <- 2

# ガンマ分布のパラメータを指定
a <- 5
b <- 6

# データ数を指定
N <- 10000

# ガンマ分布に従う乱数を生成
lambda_n <- rgamma(n = N, shape = a, rate = b)

# 1次元ガウス分布に従う乱数を生成
mu_n <- rnorm(n = N, mean = m, sd = sqrt(1 / beta / lambda_n))
mu_n[1:5]; lambda_n[1:5]
## [1]  3.55747925  1.04508782 -0.08643861 -0.16471956 -0.52733525
## [1] 0.1345578 0.7027904 1.3111707 0.9523298 1.3327713

 ガンマ分布の乱数生成関数rgamma()のデータ数(サンプルサイズ)の引数nN、パラメータの引数shape, ratea, bを指定して、精度パラメータ$\lambda$を生成します。
 $\lambda$のサンプルlambda_nを用いて、平均パラメータ$\mu$をサンプリングします。1次元ガウス分布の乱数生成関数rnorm()のデータ数(サンプルサイズ)の引数nN、平均の引数meanm、標準偏差の引数sdに精度beta * lambda_nの逆数1 / beta / lambda_nの平方根を指定します。N個の乱数ごとに標準偏差が異なります。

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

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

# 作図用のmuの値を作成
mu_vals <- seq(from = -3, to = 3, length.out = 101)

# 作図用のlambdaの値を作成
lambda_vals <- seq(from = 0.02, to = 2, length.out = 100)

# ガウス-ガンマ分布を計算
dens_df <- tidyr::tibble(
  mu = rep(mu_vals, times = length(lambda_vals)), # 確率変数mu
  lambda = rep(lambda_vals, each = length(mu_vals)), # 確率変数lambda
  N_dens = dnorm(x = mu, mean = m, sd = sqrt(1 / beta / lambda)), # muの確率密度
  Gam_dens = dgamma(x = lambda, shape = a, rate = b), # lambdaの確率密度
  density = N_dens * Gam_dens # 確率密度
)


乱数の可視化

 サンプルの散布図を作成します。

# サンプルの散布図を作成
ggplot() + 
  geom_point(data = data_df, mapping = aes(x = mu, y = lambda), color = "orange") + # サンプル
  geom_contour(data = dens_df, mapping = aes(x = mu, y = lambda, z = density, color = ..level..)) + # 元の分布
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = paste0("m=", m, ", beta=", beta, ", a=", a, ", b=", b, ", N=", N), 
       x = expression(mu), y = expression(lambda), 
       color = "density") # ラベル

f:id:anemptyarchive:20220203223940p:plain
ガウス-ガンマ分布の乱数の散布図


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

# サンプルのヒストグラムを作成
ggplot() + 
  geom_bin2d(data = data_df, mapping = aes(x = mu, y = lambda), alpha = 0.8) + # 頻度
  geom_contour(data = dens_df, mapping = aes(x = mu, y = lambda, z = density, color = ..level..)) + # 元の分布
  scale_fill_distiller(palette = "Spectral") + # 塗りつぶしの色
  scale_color_distiller(palette = "Spectral") + # 等高線の色
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = paste0("m=", m, ", beta=", beta, ", a=", a, ", b=", b, ", N=", N), 
       x = expression(mu), y = expression(lambda), 
       fill = "frequency", color = "density") # ラベル

f:id:anemptyarchive:20220203223955p:plain
ガウス-ガンマ分布の乱数のヒストグラム

 geom_bin2d()で2次元の変数に対するヒストグラムを描画できます。

 頻度を密度に変換して作図します。

# サンプルの密度を作図
ggplot() + 
  geom_density2d_filled(data = data_df, mapping = aes(x = mu, y = lambda, fill = ..level..), alpha = 0.8) + # 密度
  geom_contour(data = dens_df, mapping = aes(x = mu, y = lambda, z = density, color = ..level..)) + # 元の分布
  scale_color_viridis_c(option = "D") + # 等高線の色
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = paste0("m=", m, ", beta=", beta, ", a=", a, ", b=", b, ", N=", N), 
       x = expression(mu), y = expression(lambda), 
       fill = "frequency", color = "density") # ラベル

f:id:anemptyarchive:20220203224008p:plain
ガウス-ガンマ分布の乱数の密度

 geom_density2d_filled()で頻度を密度に変換したグラフを作成できます。

1データずつアニメーションで可視化

 データ数を指定して、サンプルを1つずつ増やしてヒストグラムの変化をアニメーションで確認します。

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

# フレーム数(データ数)を指定
N_frame <- 100

# 乱数を1つずつ生成
#lambda_n <- rep(NA, times = N_frame)
#mu_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)
  
  # 1次元ガウス分布に従う乱数を生成
  #mu_n[n] <- rnorm(n = 1, mean = m, sd = sqrt(1 / beta / lambda_n[n]))
  
  # ラベル用のテキストを作成
  label_text <- paste0("m=", m, ", beta=", beta, ", a=", a, ", b=", b, ", N=", n)
  
  # n個のサンプルを格納
  tmp_freq_df <- tidyr::tibble(
    mu = mu_n[1:n], # muのサンプル
    lambda = lambda_n[1:n], # lambdaのサンプル
    n = n, # データ数
    parameter = as.factor(label_text) # フレーム切替用のラベル
  )
  
  # n番目のサンプルを格納
  tmp_data_df <- tidyr::tibble(
    mu = mu_n[n], # muのサンプル
    lambda = lambda_n[n], # lambdaのサンプル
    n = n, # データ数
    parameter = as.factor(label_text) # フレーム切替用のラベル
  )
  
  # n回目のラベルを付与
  tmp_dens_df <- dens_df %>% 
    dplyr::mutate(
      n = n, 
      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)
  
  # 途中経過を表示
  #message("\r", appendLF = FALSE) # 前回の表示を消去
  #message("\r", "n=", n, " (", round(n / N_frame * 100, 2), "%)", appendLF = FALSE)
}
head(anime_freq_df)
## # A tibble: 6 x 4
##        mu lambda     n parameter                 
##     <dbl>  <dbl> <int> <fct>                     
## 1  3.56    0.135     1 m=0, beta=2, a=5, b=6, N=1
## 2  3.56    0.135     2 m=0, beta=2, a=5, b=6, N=2
## 3  1.05    0.703     2 m=0, beta=2, a=5, b=6, N=2
## 4  3.56    0.135     3 m=0, beta=2, a=5, b=6, N=3
## 5  1.05    0.703     3 m=0, beta=2, a=5, b=6, N=3
## 6 -0.0864  1.31      3 m=0, beta=2, a=5, b=6, N=3

 for()ループを使って、乱数を1つずつ取り出すまたは生成して、結果をanime_***_dfに追加していきます。

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

# メッシュ用の値を設定
mu_breaks     <- seq(from = min(mu_vals), to = max(mu_vals), length.out = 30)
lambda_breaks <- seq(from = min(lambda_vals), to = max(lambda_vals), length.out = 30)

# アニメーション用のサンプルのヒストグラムを作成
anime_freq_graph <- ggplot() + 
  geom_bin2d(data = anime_freq_df, mapping = aes(x = mu, y = lambda), 
             breaks = list(x = mu_breaks, y = lambda_breaks), alpha = 0.8) + # 頻度
  geom_point(data = anime_data_df, mapping = aes(x = mu, y = lambda), color = "orange", size = 5) + # サンプル
  geom_contour(data = anime_dens_df, mapping = aes(x = mu, y = lambda, z = density, color = ..level..)) + # 元の分布
  scale_fill_distiller(palette = "Spectral") + # 塗りつぶしの色
  scale_color_distiller(palette = "Spectral") + # 等高線の色
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = "{current_frame}", 
       x = expression(mu), y = expression(lambda), 
       fill = "frequency", color = "density") # ラベル

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

 geom_bin2d()breaks引数で頻度をカウントする位置を指定できます。x軸とy軸を分けて指定する場合は、リストにして渡します。カウントする位置を指定することで、サンプルの範囲によってカウントの仕方が変わらないようにします。

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

# (データが少ないと密度が計算できないため)最初のフレームを除去
n_min <- 10
tmp_freq_df <- anime_freq_df %>% 
  dplyr::filter(n > n_min)
tmp_data_df <- anime_data_df %>% 
  dplyr::filter(n > n_min)
tmp_dens_df <- anime_dens_df %>% 
  dplyr::filter(n > n_min)

# (表示が上手くいかなくなるため)因子のレベルを再設定
tmp_freq_df[["parameter"]] <- factor(as.character(tmp_freq_df[["parameter"]]), levels = unique(tmp_freq_df[["parameter"]]))
tmp_data_df[["parameter"]] <- factor(as.character(tmp_data_df[["parameter"]]), levels = unique(tmp_data_df[["parameter"]]))
tmp_dens_df[["parameter"]] <- factor(as.character(tmp_dens_df[["parameter"]]), levels = unique(tmp_dens_df[["parameter"]]))

# アニメーション用のサンプルの密度を作図
anime_prop_graph <- ggplot() + 
  geom_density2d_filled(data = tmp_freq_df, mapping = aes(x = mu, y = lambda), alpha = 0.8) + # 密度
  geom_point(data = tmp_data_df, mapping = aes(x = mu, y = lambda), color = "orange", size = 5) + # サンプル
  geom_contour(data = tmp_dens_df, mapping = aes(x = mu, y = lambda, z = density, color = ..level..)) + # 元の分布
  scale_color_viridis_c(option = "D") + # 等高線の色
  gganimate::transition_manual(parameter) + # フレーム
  xlim(c(min(mu_vals), max(mu_vals))) + # x軸の表示範囲
  ylim(c(min(lambda_vals), max(lambda_vals))) + # y軸の表示範囲
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = "{current_frame}", 
       x = expression(mu), y = expression(lambda), 
       fill = "density", color = "density") # ラベル

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


f:id:anemptyarchive:20220203224028g:plainf:id:anemptyarchive:20220203224032g:plain
ガウス-ガンマ分布の乱数の推移

 データ数が少ないと密度を計算できない(?)ので、始めの頃のデータ(行)取り除いておきます。
 データフレームからデータを除いても因子のレベル情報に含まれるとアニメーションに含まれてしまう(?)ので、因子の情報を再設定します。

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

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

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

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

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

# 乱数を1つずつ生成
anime_freq_df <- tidyr::tibble()
anime_dens_df <- tidyr::tibble()
for(i in 1:frame_num) {
  # ラベル用のテキストを作成
  label_text <- paste0(
    "m=", m, ", beta=", beta, ", a=", a, ", b=", b, ", N=", i * n_per_frame
  )
  
  # n個のサンプルを格納
  tmp_freq_df <- tidyr::tibble(
    mu = mu_n[1:(i*n_per_frame)], # muのサンプル
    lambda = lambda_n[1:(i*n_per_frame)], # lambdaのサンプル
    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_dens_df <- rbind(anime_dens_df, tmp_dens_df)
  
  # 途中経過を表示
  #message("\r", appendLF = FALSE) # 前回の表示を消去
  #message("\r", "i=", i, " (", round(i / frame_num * 100, 2), "%)", appendLF = FALSE)
}
head(anime_freq_df)
## # A tibble: 6 x 3
##        mu lambda parameter                   
##     <dbl>  <dbl> <fct>                       
## 1  3.56    0.135 m=0, beta=2, a=5, b=6, N=100
## 2  1.05    0.703 m=0, beta=2, a=5, b=6, N=100
## 3 -0.0864  1.31  m=0, beta=2, a=5, b=6, N=100
## 4 -0.165   0.952 m=0, beta=2, a=5, b=6, N=100
## 5 -0.527   1.33  m=0, beta=2, a=5, b=6, N=100
## 6  0.820   1.35  m=0, beta=2, a=5, b=6, N=100

 1フレーム当たりのデータ数を計算してn_per_frameとします。
 乱数をn_per_frame個ずつ取り出して、結果をanime_freq_dfに追加していきます。

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

# メッシュ用の値を設定
mu_breaks     <- seq(from = min(mu_vals), to = max(mu_vals), length.out = 30)
lambda_breaks <- seq(from = min(lambda_vals), to = max(lambda_vals), length.out = 30)

# アニメーション用のサンプルのヒストグラムを作成
anime_freq_graph <- ggplot() + 
  geom_bin2d(data = anime_freq_df, mapping = aes(x = mu, y = lambda), 
             breaks = list(x = mu_breaks, y = lambda_breaks), alpha = 0.8) + # 頻度
  geom_contour(data = anime_dens_df, mapping = aes(x = mu, y = lambda, z = density, color = ..level..)) + # 元の分布
  scale_fill_distiller(palette = "Spectral") + # 塗りつぶしの色
  scale_color_distiller(palette = "Spectral") + # 等高線の色
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = "{current_frame}", 
       x = expression(mu), y = expression(lambda), 
       fill = "frequency", color = "density") # ラベル

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


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

# アニメーション用のサンプルの密度を作図
anime_prop_graph <- ggplot() + 
  geom_density2d_filled(data = anime_freq_df, mapping = aes(x = mu, y = lambda), alpha = 0.8) + # 密度
  geom_contour(data = anime_dens_df, mapping = aes(x = mu, y = lambda, z = density, color = ..level..)) + # 元の分布
  scale_color_viridis_c(option = "D") + # 等高線の色
  xlim(c(min(mu_vals), max(mu_vals))) + # x軸の表示範囲
  ylim(c(min(lambda_vals), max(lambda_vals))) + # y軸の表示範囲
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Gaussian-Gamma Distribution", 
       subtitle = "{current_frame}", 
       x = expression(mu), y = expression(lambda), 
       fill = "density", color = "density") # ラベル

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


f:id:anemptyarchive:20220203224117g:plainf:id:anemptyarchive:20220203224119g:plain
ガウス-ガンマ分布の乱数の推移

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

分布の生成

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

 パラメータを$m, \beta, a, b$とするガウス-ガンマ分布

$$ \mathrm{NG}(\mu, \lambda | m, \beta, a, b) = \mathcal{N}(\mu | m, (\beta \lambda)^{-1}) \mathrm{Gam}(\lambda | a, b) $$

を事前分布として、$\mu, \lambda$を生成します。
 生成した$\mu, \lambda$を用いて、平均$\mu$・精度$\lambda$の1次元ガウス分布

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

を描画します。

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

# 1次元ガウス分布の平均パラメータを指定
m <- 0

# 1次元ガウス分布の精度パラメータの係数を指定
beta <- 2

# ガンマ分布のパラメータを指定
a <- 5
b <- 6

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

# 精度パラメータを生成
lambda_n <- rgamma(n = N, shape = a, rate = b)

# 平均パラメータを生成
mu_n <- rnorm(n = N, mean = m, sd = sqrt(1 / beta / lambda_n))
mu_n[1:5]; lambda_n[1:5]
## [1] -0.13619847 -0.64993751  1.89778694  0.79559072 -0.02703327
## [1] 0.2592207 0.4917447 0.4584893 0.8913168 1.4162055

 ガンマ分布に従う乱数を生成して、精度パラメータ$\lambda$として利用します。
 $\lambda$のサンプルlambda_nを用いて、平均パラメータ$\mu$をサンプリングします。

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

# 平均パラメータの期待値を計算
E_mu <- m

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

# 標準偏差パラメータの期待値を計算
E_sigma <- sqrt(1 / E_lambda)

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

# パラメータの期待値による1次元ガウス分布を計算
res_dens_df <- tidyr::tibble(
  x = x_vals, # 確率変数
  density = dnorm(x = x_vals, mean = E_mu, sd = E_sigma), # 確率密度
  parameter = paste0("E[mu]=", E_mu, ", E[lambda]=", round(E_lambda, 2)) %>% 
    as.factor() # 色分け用のラベル
)
head(res_dens_df)
## # A tibble: 6 x 3
##       x    density parameter              
##   <dbl>      <dbl> <fct>                  
## 1 -5.48 0.00000136 E[mu]=0, E[lambda]=0.83
## 2 -5.43 0.00000166 E[mu]=0, E[lambda]=0.83
## 3 -5.39 0.00000202 E[mu]=0, E[lambda]=0.83
## 4 -5.35 0.00000246 E[mu]=0, E[lambda]=0.83
## 5 -5.30 0.00000299 E[mu]=0, E[lambda]=0.83
## 6 -5.26 0.00000363 E[mu]=0, E[lambda]=0.83

 平均パラメータの期待値$\mathbb{E}[\mu] = m$をE_mu、精度パラメータの期待値$\mathbb{E}[\lambda] = \frac{a}{b}$をE_lambdaとします。
 また、標準偏差パラメータの期待値$\mathbb{E}[\sigma] = \sqrt{\frac{1}{\mathbb{E}[\lambda]}}$を計算してE_sigmaとします。

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

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

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

 for()ループを使って、mu_n, lambda_nから順番に値を取り出して確率密度を計算して、res_dens_dfに追加していきます。

 $N + 1$個の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("m=", m, ", beta=", beta, ", a=", a, ", b=", b, ", N=", N), 
       x = expression(x)) # ラベル

f:id:anemptyarchive:20220203224214p:plain
サンプルしたパラメータによる1次元ガウス分布のグラフ

 期待値による分布(破線)を中心に分布しているのを確認できます。

参考文献

  • C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上』,丸善出版,2012年.

おわりに

 書いておきたかった2つ目の分布まで辿り着いた。

 2022年2月2日は、モーニング娘。'22の牧野真莉愛さんの21歳のお誕生日です!!

 この圧倒的美少女感とキャラのギャップが凄い。

【次の内容】

 これの次はガウス-ウィシャート分布だと思うけどグラフにできないしどうしたものか。