はじめに
機械学習プロフェッショナルシリーズの『トピックモデル』の勉強時に自分の理解の助けになったことや勉強会資料のまとめです。トピックモデルの各種アルゴリズムを「数式」と「プログラム」から理解することを目指します。
この記事は、1.2.3項「ベータ分布」の内容です。ベータ分布の定義を説明して、正規化項と基本的な統計量を導出します。そして最後にR言語で可視化します。
【前節の内容】
【他の節一覧】
【この節の内容】
1.2.3 ベータ分布
・定義
ベータ分布とは、ベルヌーイ分布・二項分布のパラメータ$0 \leq \phi \leq 1$を生成するための確率分布である。このような分布を共役事前分布と呼ぶ。
ベータ分布は次の式で定義される。
ここで$\alpha,\ \beta$はベータ分布のパラメータであり、次の条件を満たす必要がある。
共役事前分布のパラメータのことをハイパーパラメータとも呼ぶ。
この式に従うと、ベルヌーイ分布・二項分布のパラメータ$\phi,\ 1 - \phi$を求めることができる。
まずは、$\phi$の関数である$\phi^{\alpha-1} (1 - \phi)^{\beta-1}$について考える。
この関数を確率分布(確率密度関数)として扱うためには、全ての事象を考慮した際の値が1となる必要がある。しかしこの式のままではそうはならない。そこで、正規化項(この式の逆数となる値)$\frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)}$を掛けることで対応する。
・ガンマ関数
正規化項の導出の前に、ガンマ関数
について確認する。($\exp(-u)$は指数関数で$e^{-u}$のこと)
ガンマ関数には次の3つの性質がある。
・性質1
【途中式の途中式】
ガンマ関数の定義式(1.9)より$\Gamma(x)$のとき$u^{x-1}$なので、$\Gamma(x + 1)$のとき$u^{(x+1)-1} = u^x$になる。
- 部分積分$\int_a^b f'(x) g(x) dx = [f(x) g(x)]_a^b - \int_a^b f(x) g'(x) dx$を行うために、式を変形する。
$e^{-u}$について、指数関数の積分$\int e^{ax} dx = \frac{1}{a} e^{ax}$を行うと、$\int e^{-u} du = \frac{1}{-1} e^{-u} = -e^{-u}$になる。
この$-e^{-u}$について、指数関数の微分$(e^x)' = e^x$と合成関数の微分$\{ f( g(x) )\}' = f'( g(x) ) g'(x)$より$(e^t)' = e^t * t'$を行うと、$(-e^{-u})' = (-e^{-u})' (-u)' = -e^{-u} (-1) = e^{-u}$になる。
つまり$e^{-u} = (-e^{-u})'$である(積分したものを微分したら元に戻るということ)。
- $-\exp(-u)$を$f(x)$、$(-\exp(-u))'$を$f'(x)$、$u^x$を$g(x)$、$x u^{x-1}$を$g'(x)$として部分積分する。
- 式を整理する。
- 前の因子は、$u \rightarrow \infty$のとき$\exp(-u) = 0$、また$u = 0$のとき$u^x = 0$なので、$[u^x ( -\exp(-u) )]_{0}^{\infty} = u^x * 0 - 0 * ( - \exp(-u) ) = 0$となる。
- 後の因子について、$u$と無関係な$-x$を$\int$の外に出す。
- ガンマ関数の定義より、$\Gamma(x) = \int_{0}^{\infty}u^{x-1} \exp(-u) du$で置き換える。
・性質2
【途中式の途中式】
ガンマ関数の定義式(1.9)より、$\Gamma(1)$の式を立てる。
- 部分積分を行えるように式を変形する。また$u0 = 1$である。
- 1を$g(x)$として部分積分する。
- $u \rightarrow \infty$のとき$e^{-u} = e^{-\infty} = 0$、また$u = 0$のとき$e^{-u} = e^0 = 1$となる。
・性質3
$x$が自然数の場合。
【途中式の途中式】
ガンマ関数の定義式(1.9)より、$\Gamma(x + 1)$の式を立てる。
- ガンマ関数の性質(1.10)より$\Gamma(x + 1) = x \Gamma(x)$なので、$\Gamma(x) = (x - 1) \Gamma(x - 1)$になる。
- 同様に$\Gamma(x - 1)$も$(x - 2) \Gamma(x - 2)$と置き換えられる。これを繰り返す。
- ガンマ関数の性質(1.11)より$\Gamma(1) = 1$なので、式全体が$x$から1までの自然数の積$x!$になる。
・ベータ関数
・正規化項の導出
続いて、ベータ関数
を用いて、ベータ分布の正規化項を導出する。
ベータ関数について$m = \alpha - 1,\ n = \beta - 1$とおき、$I(m, n) = \int_0^1 x^m (1 - x)^n dx$を解く(表記を分かりやすくするため、$B(\alpha,\beta) = B(m + 1, n + 1) = I(m, n)$としておく)。
【途中式の途中式】
- 部分積分$\int_a^b f'(x) g(x) dx = [f(x) g(x)]_a^b - \int_a^b f'(x) g'(x) dx$を行うために、式を変形する。
- $x^m$を$f'(x)$とすると、$x^m$を積分した$\frac{1}{m + 1} x^{m+1}$が$f(x)$になる。
- $(1 - x)^n$を$g(x)$とすると、合成関数の微分$\{g(h(x))\}' = g'( h(x) ) h'(x)$をした$((1 - x)^n)' = n (1 - x)^{n-1} (- x^{1-1}) = -n (1 - x)^{n-1}$が$g'(x)$となる。
- 式を整理する。
- $x = 0$のとき$\frac{1}{m + 1} x^{m + 1} = 0$、$x = 1$のとき$(1 - x)^n = 0$なので、$[\frac{1}{m + 1} x^{m+1}]_{0}^{1} = 0 - 0 = 0$となる。
- $-n$と$\frac{1}{m + 1}$は定数なので$\int$の外に出す。
- $I(m, n) = \int_0^1 x^m (1 - x)^n dx$なので、$I(m + 1, n - 1) = \int_0^1 x^{m+1} n (1 - x)^{n-1} dx$になる。
- $I(m, n) = \frac{n}{m + 1} I(m + 1, n - 1)$であることから、$I(m + 1, n - 1) = \frac{n - 1}{m + 2} I(m + 2, n - 2)$になる。
- 同様に$I(m + 2, n - 2) = \frac{n - 2}{m + 3} I(m + 3, n - 3)$とできる。これを$\frac{1}{m + n} I(m + n, 0)$まで繰り返す。
- それぞれ置き換える。
- 分子の$n (n - 1) \cdots 1$は、$n$から1までの自然数の積なので$n!$である。
- 定義式より、$I(m + n, 0) = \int_0^1 x^{m+n} (1 - x)^0 dx = \int_0^1 x^{m+n} dx$である。
- 後の因子を積分すると、$\int_0^1 x^{m+n} dx = \frac{1}{m + n + 1} 1^{m+n} - \frac{1}{m + n + 1} 0^{m+n} = \frac{1}{m + n + 1}$になる。
- 分母は$(m + 1)$から$(m + n + 1)$までの積なので、(分子と共に)$m!$を掛けると$(m + n + 1)!$になる。
- $m = \alpha - 1,\ n = \beta - 1$に戻す。
- ガンマ関数の性質$(x + 1)! = \Gamma(x)$より、各項を置き換える。
ベータ関数$B(\alpha,\beta) = \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha + \beta)}$が求まった。この式の逆数$\frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)}$を正規化項(積分して1とするための項)として用いることで、ベータ分布(1.8)が定義される。
・統計量
ここからはベータ分布の平均、分散、最頻値を求めていく。
・平均の導出
分布が取り得る値$\phi$とその確率密度$p(\phi)$を掛けて、0から1の範囲で積分した値が平均となる。
【途中式の途中式】
平均の定義式(1.5')より、式を立てる。
- $p(\phi) = \mathrm{Beta}(\phi | \alpha, \beta)$で置き換える。
- 正規化項はベータ関数の逆数であるため、$\frac{1}{B(\alpha, \beta)} = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)}$で置き換える。また$\phi \phi^{\alpha-1} = \phi^\alpha$で$\phi$の項をまとめる。
- $\alpha = \alpha' - 1$とおく。
- ベータ関数の定義式$B(\alpha, \beta) = \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha + \beta)} = \frac{(\alpha - 1)! (\beta - 1)!}{(\alpha + \beta - 1)!}$用いて、項を変形する。
ベータ関数の定義式
より、$B(\alpha - 1, \beta)$は
となる。この式の両辺に$\frac{\alpha - 1}{\alpha + \beta - 1}$を掛けると
$B(\alpha, \beta)$となる。この式の両辺の逆数をとり、式を整理すると
の関係が求まる。
- 定数の項を$\int$の外に出す。また$\frac{1}{B(\alpha' - 1, \beta)}$を正規化項の形に戻す。
- パラメータが$\alpha',\ \beta$のベータ分布について、 $\phi$が取り得る全ての値の積分であるため1になる。
- $\alpha' = \alpha + 1$に戻す。
・分散の導出
「$\phi$の2乗の平均」と「$\phi$の平均の2乗」との差が分散となる。そこでまずは、$\phi$の2乗の平均を求める。
【途中式の途中式】
平均の定義式(1.5')より、式を立てる。
- $p(\phi) = \mathrm{Beta}(\phi | \alpha, \beta)$で置き換える。
- 正規化項はベータ関数の逆数なので、$\frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)} = \frac{1}{B(\alpha, \beta)}$で置き換える。また$\phi^2 \phi^{\alpha-1} = \phi^{\alpha+1}$である。
- $\alpha + 1 = \alpha' - 1$とおく。
- ベータ関数の定義式$B(\alpha, \beta) = \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha + \beta)} = \frac{(\alpha - 1)! (\beta - 1)!}{(\alpha + \beta - 1)!}$を用いて、項を置き換える。
先ほどと同様に、ベータ関数の定義式から
となることが分かる。この式の両辺の逆数をとり、式を整理すると
の関係が求まる。
- 定数の項を$\int$の外に出す。また$\frac{1}{B(\alpha', \beta)}$を正規化項の形に戻す。
- パラメータが$\alpha',\ \beta$のベータ分布について、$0 \leq \phi \leq 1$の範囲での積分なので1になる。
- $\alpha' = \alpha + 2$に戻す。
「$\phi$の2乗の平均」と「$\phi$の平均の2乗」との差を求める。
・最頻値の導出
最頻値は曲線(グラフ)の頂点の値であることから、$p(\phi)$を微分して、$\frac{d p(\phi)}{d \phi} = 0$となる$\phi$である。
ベータ分布の確率密度関数
を$\phi$で微分する。
【途中式の途中式】
そのためにまずは、$\phi^{\alpha-1} (1 - \phi)^{\beta-1}$の微分について考える。
$\phi^{\alpha-1}$を$f(\phi)$、$(1 - \phi)^{\beta-1}$を$g(\phi)$とするとこの式の微分は、積の微分$\{f(\phi) g(\phi)\}' = f'(\phi) g(\phi) + f(\phi) g'(\phi)$である。
$\phi^{\alpha-1}$を微分すると$(\alpha - 1) \phi^{\alpha-2}$になる。
$(1 - \phi)^{\beta-1}$の微分は、合成関数の微分$\{ g( h(\phi) ) \} ' = g'( h(\phi) ) h'(\phi)$である。$(1 - \phi)^{\beta-1}$を$g( h(\phi) )$、$1 - \phi$を$h(\phi)$とすると
となる。
【途中式の途中式】
- $(1 - \phi)^{\beta-1} = (1 - \phi) (1 - \phi)^{\beta-2}$、$\phi^{\alpha-1} = \phi \phi^{\alpha-2}$に項を分割する。
- $\phi^{\alpha-2}$と$(1 - \phi)^{\beta-2}$を括り出す。
$\frac{d p(\phi)}{d \phi} = 0$となる$\phi$を求める。
ここで式よりパラメータ$\alpha,\ \beta$は、$\alpha > 1,\ \beta > 1$を満たす必要がある。
・可視化
$\alpha = 0.5, \beta = 0.5$のときの分布を確認する。
# 利用パッケージ library(tidyverse) # パラメーターを指定 alpha <- 0.5 beta <- 0.5 # 作図 tibble( phi = seq(0.01, 0.99, 0.01), # phiの値 B = gamma(alpha + beta) / gamma(alpha) / gamma(beta), # 正規化項 density = B * phi^(alpha - 1) * (1 - phi)^(beta - 1) # 確率密度 ) %>% ggplot(mapping = aes(x = phi, y = density)) + # データ geom_line(color = "#00A968") + # 折れ線グラフ labs(title = "Beta Distribution:", subtitle = paste0("alpha=", alpha, ", beta=", beta), x = expression(phi)) # ラベル
確率密度の計算は、ベータ関数beta()
やベータ分布に従う確率密度dbeta()
を使っても求められる。
・パラメータの違いによる分布の変化
・コード(クリックで展開)
# パラメーターの指定 alpha_vec <- c(1, 2, 2, 0.9, 0.8) beta_vec <- c(1, 2, 4, 0.7, 1.2) # 作図用のデータフレームを作成 res_df <- tibble() for(i in seq_along(alpha_vec)) { # 分布を計算 tmp_df <- tibble( phi = seq(0.01, 0.99, 0.01), # phiの値 B = 1 / beta(alpha_vec[i], beta_vec[i]), # 正規化項 density = B * phi^(alpha_vec[i] - 1) * (1 - phi)^(beta_vec[i] - 1), # 確率密度 parameter = paste0("alpha=", alpha_vec[i], ", beta=", beta_vec[i]) # 作図用のラベル ) # 計算結果を結合 res_df <- rbind(res_df, tmp_df) } # 作図 ggplot(data = res_df, mapping = aes(x = phi, y = density, color = parameter)) + # データの指定 geom_line() + # 折れ線グラフ scale_color_manual(values = c("#FFC0CB", "#FF0000", "#FFFF00", "#EE82EE", "#7FFF00")) + # 線の色(不必要) theme_dark() + # 背景色(不必要) labs(title = "Beta Distribution", x = expression(phi)) # ラベル
(色の指定はただの趣味です。指定するパラメータの個数と色の数を合わせる必要があるため、そもそも使わない方がいいです。)
パラメータ$\phi$の変化による分布の変化をgif画像で確認する。
・コード(クリックで展開)
# 利用するライブラリ library(tidyverse) library(gganimate) # パラメータの組み合わせ parameter_df = tibble( alpha = rep(seq(0, 10, by = 0.25), times = 41), beta = rep(seq(0, 10, by = 0.25), each = 41) ) # 作図用のデータフレームを作成 res_df <- tibble() for(i in 1:nrow(parameter_df)) { # 分布を計算 tmp_df <- tibble( phi = seq(0.01, 0.99, 0.01), # phiの値 density = dbeta( x = phi, shape1 = parameter_df[["alpha"]][i], shape2 = parameter_df[["beta"]][i] ), # 確率密度 parameter = paste0( "alpha=", parameter_df[["alpha"]][i], ", beta=", parameter_df[["beta"]][i] ) %>% as.factor() # 作図用のラベル ) # 計算結果を結合 res_df <- rbind(res_df, tmp_df) } # 作図 graph_data <- ggplot(data = res_df, mapping = aes(x = phi, y = density)) + # データの指定 geom_line(color = "#00A968") + # 折れ線グラフ transition_manual(parameter) + # フレーム labs(title = "Beta Distribution", subtitle = "{current_frame}", x = expression(phi)) # ラベル # gif画像を作成 animate(graph_data, nframes = nrow(parameter_df), fps = 25)
gganimate::transition_manual()
に時系列となる列を指定する。図を構成できたらgganimate::animate()
で出力する。
Hello World!から一歩進んだ pic.twitter.com/y4Xou5Wxzz
— しょこ📚 (@anemptyarchive) 2020年6月23日
重くて貼れないようなのでこれで代用…
参考書籍
- 岩田具治(2015)『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社
- 須山敦志(2017)『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)講談社
おわりに
2019/08/19:加筆修正しました。
部分積分ちんぷんかんぷん。ベータ分布のy軸はx軸の値の出やすさなんでしょうけど、具体的な数値って何を意味しているのでしょうか?(確率密度ですよ。)
2020/06/25:加筆修正しました。
【次節の内容】