からっぽのしょこ

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

4.3.2:ポアソン混合モデルのギブスサンプリング【緑ベイズ入門のノート】

はじめに

 『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「Rで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。

 この記事は4.3.2項の内容になります。ポアソン混合モデルにおけるギブスサンプリングを用いた近似事後分布を導出し、Rで実装します。

 省略してある内容等ありますので、本と併せて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

4.3.2 ギブスサンプリング

 ポアソン混合モデルの事後分布からパラメータ$\boldsymbol{\lambda},\ \boldsymbol{\pi}$と潜在変数$\boldsymbol{S}$をサンプルするアルゴリズムを導出する。

・潜在変数の条件付き分布の導出

 初めに、潜在変数$\boldsymbol{S}$をサンプルするための条件付き分布(事後分布)を求めていく。

$$ \begin{align*} p(\boldsymbol{S} | \boldsymbol{X}, \boldsymbol{\lambda}, \boldsymbol{\pi}) &= \frac{ p(\boldsymbol{X}, \boldsymbol{S}, \boldsymbol{\lambda}, \boldsymbol{\pi}) }{ p(\boldsymbol{X}, \boldsymbol{\lambda}, \boldsymbol{\pi}) } \\ &\propto p(\boldsymbol{X}, \boldsymbol{S}, \boldsymbol{\lambda}, \boldsymbol{\pi}) \\ &= p(\boldsymbol{X} | \boldsymbol{S}, \boldsymbol{\lambda}) p(\boldsymbol{S} | \boldsymbol{\pi}) p(\boldsymbol{\lambda}) p(\boldsymbol{\pi}) \\ &\propto p(\boldsymbol{X} | \boldsymbol{S}, \boldsymbol{\lambda}) p(\boldsymbol{S} | \boldsymbol{\pi}) \\ &= \prod_{n=1}^N p(x_n | \boldsymbol{s}_n, \boldsymbol{\lambda}) p(\boldsymbol{s}_n | \boldsymbol{\pi}) \tag{4.33} \end{align*} $$

$\boldsymbol{S}$と関係のない項を適宜省略している。

 この式の具体的な確率分布を計算していく。前の項は式(4.28)から

$$ \begin{align*} \ln p(x_n | \boldsymbol{s}_n, \boldsymbol{\lambda}) &= \ln \prod_{k=1}^K \mathrm{Poi}(x_n | \lambda_k)^{s_{n,k}} \\ &= \sum_{k=1}^K s_{n,k} \ln \mathrm{Poi}(x_n | \lambda_k) \\ &= \sum_{k=1}^K s_{n,k} \ln \frac{\lambda_k^{x_n}}{x_n!} \exp(- \lambda_k) \\ &= \sum_{k=1}^K s_{n,k} ( x_n \ln \lambda_k - \lambda_k ) + \mathrm{const.} \tag{4.34} \end{align*} $$

となる。また後の項は式(4.2)から

$$ \begin{align*} \ln p(\boldsymbol{s}_n | \boldsymbol{\pi}) &= \ln \mathrm{Cat}(\boldsymbol{s}_n | \boldsymbol{\pi}) \\ &= \ln \prod_{k=1}^K \pi_k^{s_{n,k}} \\ &= \sum_{k=1}^K s_{n,k} \ln \pi_k \tag{4.35} \end{align*} $$

となる。

 よって2つの式を合わせると、式(4.33)は

$$ \begin{align*} \ln p(\boldsymbol{S} | \boldsymbol{X}, \boldsymbol{\lambda}, \boldsymbol{\pi}) &= \sum_{n=1}^N \sum_{k=1}^K s_{n,k} ( x_n \ln \lambda_k - \lambda_k ) s_{n,k} \ln \pi_k + \mathrm{const.} \\ &= \sum_{n=1}^N \sum_{k=1}^K s_{n,k} ( x_n \ln \lambda_k - \lambda_k + \ln \pi_k ) + \mathrm{const.} \tag{4.36} \end{align*} $$

となる。

 この式について

$$ \eta_{n,k} \propto \exp\{ x_n \ln \lambda_k - \lambda_k + \ln \pi_k \} \tag{4.38} $$

とおくと、式の形状から

$$ p(\boldsymbol{S} | \boldsymbol{X}, \boldsymbol{\lambda}, \boldsymbol{\pi}) = \prod_{n=1}^N \prod_{k=1}^K \eta_{n,k}^{s_{n,k}} = \prod_{n=1}^N \mathrm{Cat}(\boldsymbol{s}_n | \boldsymbol{\eta}_n) $$

パラメータ$\boldsymbol{\eta}$を持つカテゴリ分布となることが分かる。ただし

$$ \boldsymbol{\eta}_n = (\eta_{n,1},\eta_{n,2} \cdots, \eta_{n,K}),\ \sum_{k=1}^K \eta_{n,k} = 1 $$

である。

 つまりハイパーパラメータ$\boldsymbol{\eta}_n$を計算することで、カテゴリ分布から潜在変数$\boldsymbol{s}_n$をサンプリングできる。

$$ \boldsymbol{s}_n \sim \mathrm{Cat}(\boldsymbol{s}_n | \boldsymbol{\eta}_n) \tag{4.37} $$


・パラメータの条件付き分布の導出

 続いて、式(4.32)からパラメータ$\boldsymbol{\lambda},\ \boldsymbol{\pi}$をサンプルするための条件付き分布を求めていく。

$$ \begin{align*} p(\boldsymbol{\lambda}, \boldsymbol{\pi} | \boldsymbol{X}, \boldsymbol{S}) &= \frac{ p(\boldsymbol{X}, \boldsymbol{S}, \boldsymbol{\lambda}, \boldsymbol{\pi}) }{ p(\boldsymbol{X}, \boldsymbol{S}) } \\ &\propto p(\boldsymbol{X}, \boldsymbol{S}, \boldsymbol{\lambda}, \boldsymbol{\pi}) \\ &= p(\boldsymbol{X} | \boldsymbol{S}, \boldsymbol{\lambda}) p(\boldsymbol{S} | \boldsymbol{\pi}) p(\boldsymbol{\lambda}) p(\boldsymbol{\pi}) \tag{4.39} \end{align*} $$

 この式を用いて、パラメータ$\boldsymbol{\lambda},\ \boldsymbol{\pi}$の具体的な確率分布を計算していく。

 まずは$\boldsymbol{\lambda}$に関する項を取り出すと

$$ \begin{align*} \ln p(\boldsymbol{X} | \boldsymbol{S}, \boldsymbol{\lambda}) p(\boldsymbol{\lambda}) &= \ln \prod_{n=1}^N p(x_n | \boldsymbol{s}_n, \boldsymbol{\lambda}) + \ln \prod_{k=1}^K p(\lambda_k) \\ &= \ln \prod_{n=1}^N \prod_{k=1}^K \mathrm{Poi}(x_n | \lambda_k)^{s_{n,k}} + \ln \prod_{k=1}^K \mathrm{Gam}(\lambda_k | a, b) \\ &= \sum_{k=1}^K \left\{ \sum_{n=1}^N s_{n,k} \ln \frac{\lambda_k^{x_n}}{x_n!} \exp(- \lambda_k) + \ln C_G(a, b) \lambda_k^{a-1} \exp(- b \lambda_k) \right\} \\ &= \sum_{k=1}^K \left\{ \sum_{n=1}^N s_{n,k} ( x_n \ln \lambda_k - \lambda_k ) + (a - 1) \ln \lambda_k - b \lambda_k \right\} + \mathrm{const.} \\ &= \sum_{k=1}^K \left\{ \Bigl(\sum_{n=1}^N s_{n,k} x_n + a - 1\Bigr) \ln \lambda_k - \Bigl(\sum_{n=1}^N s_{n,k} + b\Bigr) \lambda_k \right\} + \mathrm{const.} \tag{4.40} \end{align*} $$

となる。

 この式について

$$ \begin{align*} \hat{a}_k &= \sum_{n=1}^N s_{n,k} x_n + a \\ \hat{b}_k &= \sum_{n=1}^N s_{n,k} + b \tag{4.42} \end{align*} $$

とおくと、式の形状から

$$ \ln p(\boldsymbol{X} | \boldsymbol{S}, \boldsymbol{\lambda}) p(\boldsymbol{\lambda}) = \sum_{k=1}^K \hat{a}_k \ln \lambda_k - \hat{b}_k \lambda_k + \mathrm{const.} = \sum_{k=1}^K \ln \mathrm{Gam}(\lambda_k | \hat{a}_k, \hat{b}_k) $$

パラメータ$\hat{a}_k,\ \hat{b}_k$を持つガンマ分布であることが分かる。

 つまりハイパーパラメータ$\hat{a}_k,\ \hat{b}_k$を計算することで、ガンマ分布からパラメータ$\lambda_k$をサンプリングできる。

$$ \lambda_k \sim \mathrm{Gam}(\lambda_k | \hat{a}_k, \hat{b}_k) \tag{4.41} $$


 $\boldsymbol{\pi}$についても同様に、式(4.39)から$\boldsymbol{\pi}$に関係する項を取り出すと

$$ \begin{align*} \ln p(\boldsymbol{S} | \boldsymbol{\pi}) p(\boldsymbol{\pi}) &= \ln \prod_{n=1}^N p(\boldsymbol{s}_n | \boldsymbol{\pi}) + \ln p(\boldsymbol{\pi}) \\ &= \ln \prod_{n=1}^N \mathrm{Cat}(\boldsymbol{s}_n | \boldsymbol{\pi}) + \ln \mathrm{Dir}(\boldsymbol{\pi} | \boldsymbol{\alpha}) \\ &= \ln \prod_{n=1}^N \prod_{k=1}^K \pi_k^{s_{n,k}} + \ln C_D(\boldsymbol{\alpha}) \prod_{k=1}^K \pi_k^{\alpha_k-1} \\ &= \sum_{k=1}^K \left\{ \sum_{n=1}^N s_{n,k} \ln \pi_k^{} + (\alpha_k - 1) \ln \pi_k \right\} + \mathrm{const.} \\ &= \sum_{k=1}^K \left(\sum_{n=1}^N s_{n,k} + \alpha_k - 1\right) \ln \pi_k + \mathrm{const.} \tag{4.43} \end{align*} $$

となる。

 この式について

$$ \hat{\alpha}_k = \sum_{n=1}^N s_{n,k} + \alpha_k \tag{4.45} $$

とおくと、式の形状から

$$ \ln p(\boldsymbol{S} | \boldsymbol{\pi}) p(\boldsymbol{\pi}) = \sum_{k=1}^K (\hat{\alpha}_k - 1) \ln \pi_k + \mathrm{const.} = \sum_{k=1}^K \ln \mathrm{Dir}(\boldsymbol{\pi} | \hat{\boldsymbol{\alpha}}) $$

パラメータ$\hat{\boldsymbol{\alpha}}$を持つディリクレ分布となることが分かる。ただし

$$ \hat{\boldsymbol{\alpha}} = (\hat{\alpha}_1, \hat{\alpha}_2, \cdots, \hat{\alpha}_K) $$

である。

 つまりハイパーパラメータ$\hat{\boldsymbol{\alpha}}$を計算することで、ディリクレ分布からパラメータ$\boldsymbol{\pi}$をサンプリングできる。

$$ \boldsymbol{\pi} \sim \mathrm{Dir}(\boldsymbol{\pi} | \hat{\boldsymbol{\alpha}}) \tag{4.44} $$


・Rでやってみる

 では、手を動かしてプログラムからアルゴリズムの理解を深めましょう。アルゴリズム4.2を参考に実装します。

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

 必要なパッケージを読み込みます。

・真の観測モデルの設定
# (観測)データ数を指定
N <- 100

# 真のパラメータを指定
lambda_truth <- c(5, 25)
pi_truth <- c(0.3, 0.7)

 生成するデータ数$N$を指定します。

 パラメータの真の値を指定します。
 データ生成に直接関わるパラメータ$\boldsymbol{\lambda}$の真の値をlambda_truth、その混合率でもあるパラメータ$\boldsymbol{\pi}$をpi_truthとします。ここで設定した値を推定するのが目的です。

# クラスタ数
K <- length(lambda_truth)

# クラスタ(潜在変数)
s_nk <- rmultinom(n =  N, size = 1, prob = pi_truth) %>% 
  t()

# (観測)データXを生成
x_n <- rpois(n = N, lambda = apply(lambda_truth ^ t(s_nk), 2, prod))

 lambda_truthpi_truthに設定した要素数がクラスタ数$K$となります。この例では2とします。

 データを生成するにあたり、まずは各データのクラスタの集合$\boldsymbol{S}$をカテゴリ分布に従い生成します。ただしカテゴリ分布に従う乱数生成関数は用意されていないため、多項分布rmultinom()の引数をsize = 1として用います。結果はマトリクスで返ってくるため、ベクトルに変換しておきます。

 s_nklambda_truthを用いてN個のデータを生成します。
 ポアソン分布に従う乱数にはrpois()を使います。引数lambdaには、apply(lambda_truth ^ t(s_nk), 2, prod)(式(4.28))を指定します。

 観測データを確認しましょう。

# (観測)データを確認
tibble(
  x = x_n
) %>% 
  ggplot(aes(x = x)) + 
    geom_bar(fill = "#56256E")

f:id:anemptyarchive:20200417214512p:plain
観測データ$\boldsymbol{X}$のヒストグラム

このデータを用いて推定を行います。当然データ数が多いほど推定精度が良くなります。

・パラメータの設定
# 試行回数を指定
Iter <- 50

 繰り返し学習する回数を指定します。

# ハイパーパラメータa,bを指定
a <- 1
b <- 1

# パラメータlambdaを生成
lambda_k <- rgamma(n = K, shape = a, rate = b)

 $\boldsymbol{\lambda}$の事前分布のパラメータ$a,\ b$をそれぞれ指定します。

 指定したa, bを使い、$\boldsymbol{\lambda}$を生成します。

# ハイパーパラメータalphaを指定
alpha_k <- rep(2, K)

# パラメータpiを生成
pi_k <- MCMCpack::rdirichlet(n = 1, alpha = alpha_k) %>% 
  as.vector()

 同様にハイパーパラメータ$\boldsymbol{\alpha}$を指定し、$\boldsymbol{\pi}$を生成します。
 ディリクレ分布に従う乱数生成は、MCMCpackパッケージのrdirichletを用います。この関数も結果がマトリクスで返ってくるため、ベクトルに変換します。

・ギブスサンプリング

・コード全体(クリックで展開)

# 受け皿を用意
eta_nk <- matrix(0, nrow = N, ncol = K)
a_hat_k <- rep(0, K)
b_hat_k <- rep(0, K)

# 推移の確認用データフレームを作成
trace_parameter_df <- tibble(
  a_hat = rep(a, K), 
  b_hat = rep(b, K), 
  alpha_hat = alpha_k, 
  cluster = as.factor(1:K), 
  Iteration = 0 # 初期値
)

for(i in 1:Iter) {
  
  for(n in 1:N) {
    
    # ハイパーパラメータeta_nを計算:式(4.38)
    tmp_eta_k <- exp(x_n[n] * log(lambda_k) - lambda_k + log(pi_k))
    eta_nk[n, ] <- tmp_eta_k / sum(tmp_eta_k) # 正規化
    
    # 潜在変数s_nをサンプリング:式(4.37)
    s_nk[n, ] <- rmultinom(n =  1, size = 1, prob = eta_nk[n, ]) %>% 
      as.vector()
  }
  
  for(k in 1:K) {
    
    # ハイパーパラメータhat{a}_k,hat{b}_kを計算:式(4.42)
    a_hat_k[k] <- sum(s_nk[, k] * x_n) + a
    b_hat_k[k] <- sum(s_nk[, k]) + b
    
    # パラメータlambda_kをサンプリング:式(4.41)
    lambda_k <- rgamma(n = K, shape = a_hat_k, rate = b_hat_k)
  }
  
  # ハイパーパラメータをhat{alpha}を計算:式(4.45)
  alpha_hat_k <- apply(s_nk, 2, sum) + alpha_k
  
  # パラメータpiをサンプリング:式(4.44)
  pi_k <- MCMCpack::rdirichlet(n = 1, alpha = alpha_hat_k) %>% 
    as.vector()
  
  # 推移の確認用データフレームを作成
  tmp_parameter_df <- data.frame(
    a_hat = a_hat_k, 
    b_hat = b_hat_k, 
    alpha_hat = alpha_hat_k, 
    cluster = as.factor(1:K), 
    Iteration = i # 試行回数
  )
  # 結合
  trace_parameter_df <- rbind(trace_parameter_df, tmp_parameter_df)
}

 添字を使って代入するために、予め変数を用意しておく必要があります。
 また各パラメータの更新値の推移を確認するため、学習する度に値を保存しておきます。


・潜在変数$\boldsymbol{S}$のサンプリング

for(n in 1:N) {
  
  # ハイパーパラメータeta_nを計算:式(4.38)
  tmp_eta_k <- exp(x_n[n] * log(lambda_k) - lambda_k + log(pi_k))
  eta_nk[n, ] <- tmp_eta_k / sum(tmp_eta_k) # 正規化
  
  # 潜在変数s_nをサンプリング:式(4.37)
  s_nk[n, ] <- rmultinom(n =  1, size = 1, prob = eta_nk[n, ]) %>% 
    as.vector()
}

 式(4.38)の計算を行い、ハイパーパラメータ$\boldsymbol{\eta}_n$を更新します。

 更新したパラメータeta_nk[n, ]を用いて、潜在変数$\boldsymbol{s}_n$をサンプリングします(式(4.37))。

 これで学習したパラメータを用いた(事後分布からの)$\boldsymbol{s}_n$のサンプリングを行えます。この処理を1から$N$まで繰り返し行うことで、全ての潜在変数のサンプルが得られます。

・パラメータ$\boldsymbol{\lambda}$のサンプリング

for(k in 1:K) {
  
  # ハイパーパラメータhat{a}_k,hat{b}_kを計算:式(4.42)
  a_hat_k[k] <- sum(s_nk[, k] * x_n) + a
  b_hat_k[k] <- sum(s_nk[, k]) + b
  
  # パラメータlambda_kをサンプリング:式(4.41)
  lambda_k <- rgamma(n = K, shape = a_hat_k, rate = b_hat_k)
}

 式(4.42)の計算を行い、ハイパーパラメータ$\hat{a}_k,\ \hat{b}_k$を更新します。

 更新したパラメータa_hat_k[k], b_hat_k[k]を用いて、パラメータ$\lambda_k$をサンプリングします(式(4.41))。

 これで学習したパラメータを用いた(事後分布からの)$\lambda_k$のサンプリングを行えます。この処理を1から$K$まで繰り返し行うことで、全てのパラメータのサンプルが得られます。

・パラメータ$\boldsymbol{\pi}$のサンプリング

# ハイパーパラメータをhat{alpha}を計算:式(4.45)
alpha_hat_k <- apply(s_nk, 2, sum) + alpha_k

# パラメータpiをサンプリング:式(4.44)
pi_k <- MCMCpack::rdirichlet(n = 1, alpha = alpha_hat_k) %>% 
  as.vector()

 式(4.45)の計算を行い、ハイパーパラメータ$\hat{\alpha}_k$を更新します。

 更新したパラメータalpha_hat_kを用いて、パラメータ$\boldsymbol{\pi}$をサンプリングします(式(4.44))。

 $\boldsymbol{\pi}$については、全て要素を1回でサンプリングするためループ処理は行いません。

・推定結果の確認

 ggplot2パッケージを利用して、パラメータ$\boldsymbol{\lambda},\ \boldsymbol{\pi}$の事後分布を可視化します。

## lambdaの事後分布
# 作図用のデータフレームを作成
posterior_lambda_df <- tibble()
for(k in 1:K) {
  tmp_lambda_df <- tibble(
    x = seq(0, max(x_n), by = 0.01), 
    density = dgamma(x, shape = a_hat_k[k], rate = b_hat_k[k]), 
    cluster = as.factor(k)
  )
  posterior_lambda_df <- rbind(posterior_lambda_df, tmp_lambda_df)
}

# 作図
ggplot(posterior_lambda_df, aes(x, density, color = cluster)) + 
  geom_line() + # 折れ線グラフ
  scale_color_manual(values = c("#00A968", "orange")) + # グラフの色(不必要)
  geom_vline(xintercept = lambda_truth, color = "pink", linetype = "dashed") + # 垂直線
  labs(title = "Poisson mixture model:Gibbs sampling", 
       subtitle = paste0("a_hat=(", paste0(a_hat_k, collapse = ", "), 
                         "), b_hat=(", paste0(b_hat_k, collapse = ", "), ")")) # ラベル

f:id:anemptyarchive:20200417215550p:plain
$\boldsymbol{\lambda}$の近似事後分布

## piの事後分布(K=2の場合)
# 作図用のデータフレームを作成
posterior_pi_df <- tibble()
for(k in 1:K) {
  tmp_pi_df <- tibble(
    x = seq(0, 1, by = 0.01), 
    density = dbeta(x, shape1 = alpha_hat_k[k], shape2 = alpha_hat_k[2 / k]), 
    cluster = as.factor(k)
  )
  posterior_pi_df <- rbind(posterior_pi_df, tmp_pi_df)
}

# 作図
ggplot(posterior_pi_df, aes(x, density, color = cluster)) + 
  geom_line() + # 折れ線グラフ
  scale_color_manual(values = c("#00A968", "orange")) + # グラフの色(不必要)
  geom_vline(xintercept = pi_truth, color = "pink", linetype = "dashed") + # 垂直線
  labs(title = "Poisson mixture model:Gibbs sampling", 
       subtitle = paste0("alpha_hat=(", paste0(alpha_hat_k, collapse = ", "), ")")) # ラベル

f:id:anemptyarchive:20200417215735p:plain
$\boldsymbol{\pi}$の近似事後分布

 $\boldsymbol{\pi}$については、このコードでグラフ化できるのは$N = 2$のときだけになります。(詳しくは 3.2.2:カテゴリ分布の学習と予測【緑ベイズ入門のノート】 - からっぽのしょこをご参照ください。)

 またscale_color_manual()は使わない方が無難です(私の趣味の色です)。色を指定する場合は、クラスタ数$K$と同じ数の色を指定する必要があります。

・推移の確認

・パラメータの事後分布の推移

 gganimateパッケージを利用して、パラメータ$\boldsymbol{\lambda}$の事後分布の推移をgif画像として出力するためのコードです。

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

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

## lambdaの事後分布の推移
# 作図用のデータフレームを作成
trace_lambda_df <- tibble()
for(i in 1:(Iter + 1)) {
  for(k in 1:K) {
    tmp_parameter <- trace_parameter_df %>% 
      filter(Iteration == i - 1, cluster == k)
    tmp_lambda_df <- tibble(
      x = seq(0, max(x_n), by = 0.001), 
      density = dgamma(x, shape = tmp_parameter[["a_hat"]], rate = tmp_parameter[["b_hat"]]), 
      cluster = tmp_parameter[["cluster"]], 
      Iteration = tmp_parameter[["Iteration"]]
    )
    trace_lambda_df <- rbind(trace_lambda_df, tmp_lambda_df)
  }
}

# 作図
trace_lambda_graph <- ggplot(trace_lambda_df, aes(x, density, color = cluster)) + 
  geom_line() + # 折れ線グラフ
  scale_color_manual(values = c("#00A968", "orange")) + # グラフの色(不必要)
  geom_vline(xintercept = lambda_truth, 
             color = "pink", linetype = "dashed") + # 垂直線
  transition_manual(Iteration) + # フレーム
  labs(title = "Poisson mixture model:Gibbs sampling", 
       subtitle = "i={current_frame}") # ラベル

# 描画
animate(trace_lambda_graph, nframes = Iter + 1, fps = 10)

f:id:anemptyarchive:20200417220059g:plain
$\boldsymbol{\lambda}$の近似事後分布の推移

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

## piの事後分布の推移
# 作図用のデータフレームを作成
trace_pi_df <- tibble()
for(i in 1:(Iter + 1)) {
  for(k in 1:K) {
    tmp_parameter <- trace_parameter_df %>% 
      filter(Iteration == i - 1)
    tmp_pi_df <- tibble(
      x = seq(0, 1, by = 0.01), 
      density = dbeta(
        x, 
        shape1 = tmp_parameter[["alpha_hat"]][k], 
        shape2 = tmp_parameter[["alpha_hat"]][2 / k]
      ), 
      cluster = as.factor(k), 
      Iteration = i - 1
    )
    trace_pi_df <- rbind(trace_pi_df, tmp_pi_df)
  }
}

# 作図
trace_pi_graph <- ggplot(trace_pi_df, aes(x, density, color = cluster)) + 
  geom_line() + # 折れ線グラフ
  scale_color_manual(values = c("#00A968", "orange")) + # グラフの色(不必要)
  geom_vline(xintercept = pi_truth, 
             color = "pink", linetype = "dashed") + # 垂直線
  transition_manual(Iteration) + # フレーム
  labs(title = "Poisson mixture model:Gibbs sampling", 
       subtitle = "i={current_frame}") # ラベル

# 描画
animate(trace_pi_graph, nframes = Iter + 1, fps = 10)

f:id:anemptyarchive:20200417220217g:plain
$\boldsymbol{\pi}$の近似事後分布の推移

 (作図用のデータフレーム作成はもっと上手いやり方があるはず。取り急ぎfor()でやっつけた。)

 gganimate::transition_manual()に時系列を指定して、gganimate::animate()でgifファイルに変換します。フレーム数は初期値の0を含むためイタレーション数+1となります。

・ハイパーパラメータの推移

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

# 作図
trace_parameter_df %>% 
  pivot_longer(
    cols = -c(cluster, Iteration), 
    names_to = "parameter", 
    names_ptypes = list(parameter = factor()), 
    values_to = "value"
  ) %>% 
  mutate(parameters = paste(parameter, cluster, sep = "_")) %>% 
  ggplot(aes(Iteration, value, color = parameters)) + 
    geom_line() + # 垂直線
    labs(title = "Poisson mixture model:Gibbs sampling") # ラベル

f:id:anemptyarchive:20200417220905p:plain
ハイパーパラメータの推移

参考文献

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

おわりに

 白トピ本3章の実装を進めていました(ほぼほぼ完成しました(でも3月中にやりきるはずだった))が、ギブスサンプリングの理解不足を感じて色々飛ばしてこれをやることにしました。

 世界中で制限された日々が続いていますが、そんな中で立ち上がったOsaka.Rのオンライン朝モク会に参加しています。この記事の半分はそのもくもく会中に書きました。
 普段ならほげっとしてるないし寝ている時間に勉強すると、1日を長く感じてとても良いです。それだけでなく、朝に頭を使うと他の時間の頭の働きも良くなったのも嬉しいです。「やる気を出すにはまずやること」の最初の1歩のきっかけをもらえるのは、(意思の弱い)自分にとってとてもありがたいことです。折角なのでこのまま習慣にしたいなぁ。

 というわけで意思の弱い人もそうでない人も、大阪関係なく参加OKとのことなので是非ご一緒しましょう!

osaka-r.connpass.com


【次節の内容】

www.anarchive-beta.com