からっぽのしょこ

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

3.2.3:ポアソン分布の学習と予測【緑ベイズ入門のノート】

はじめに

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

 この記事は、3.2.3項の内容です。「尤度関数をポアソン分布」、「事前分布をガンマ分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」を導出します。

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

【実装編】

www.anarchive-beta.com

www.anarchive-beta.com

【他の節の内容】

www.anarchive-beta.com

【この節の内容】

3.2.3 ポアソン分布の学習と予測

 ポアソン分布に従うと仮定するN個の観測データ$\mathbf{X} = \{x_1, x_2, \cdots, x_N\}$を用いて、パラメータ$\lambda$の事後分布と未観測のデータ$x_{*}$の予測分布を求めていく。

 $\mathbf{X}$の各データが独立に発生しているとの仮定の下で、尤度$p(\mathbf{X} | \lambda)$は

$$ \begin{aligned} p(\mathbf{X} | \lambda) &= p(x_1, x_2, \cdots, x_N | \lambda) \\ &= p(x_1 | \lambda) p(x_2 | \lambda) \cdots p(x_N | \lambda) \\ &= \prod_{n=1}^N p(x_n | \lambda) \end{aligned} $$

と分解できる。

・事後分布の計算

 まずは、$\lambda$の事後分布$p(\lambda | \mathbf{X})$を導出する。

 $\mathbf{X}$が与えられた下での$\lambda$の事後分布は、ベイズの定理より

$$ \begin{align} p(\lambda | \mathbf{X}) &= \frac{ p(\mathbf{X} | \lambda) p(\lambda | a, b) }{ p(\mathbf{X}) } \\ &= \frac{ \left\{ \prod_{n=1}^N p(x_n | \lambda) \right\} p(\lambda | a, b) }{ p(\mathbf{X}) } \\ &= \frac{ \left\{ \prod_{n=1}^N \mathrm{Poi}(x_n | \lambda) \right\} \mathrm{Gam}(\lambda | a, b) }{ p(\mathbf{X}) } \\ &\propto \left\{ \prod_{n=1}^N \mathrm{Poi}(x_n | \lambda) \right\} \mathrm{Gam}(\lambda | a, b) \tag{3.35} \end{align} $$

となる。分母の$p(\mathbf{X})$は$\lambda$に影響しないため省略して、比例関係のみに注目する。省略した部分については、最後に正規化することで対応できる。

 この分布の具体的な形状を明らかにしていく。対数をとって指数部分の計算を分かりやすくして、$\lambda$に関して整理すると

$$ \begin{align} \ln p(\lambda | \mathbf{X}) &= \ln \prod_{n=1}^N \mathrm{Poi}(x_n | \lambda) + \ln \mathrm{Gam}(\lambda | a, b) - \ln p(\mathbf{X}) \\ &= \sum_{n=1}^N \ln \mathrm{Poi}(x_n | \lambda) + \ln \mathrm{Gam}(\lambda | a, b) + \mathrm{const.} \\ &= \sum_{n=1}^N \ln \left\{ \frac{\lambda^{x_n}}{x_n!} \exp(-\lambda) \right\} + \ln \Bigl\{ C_G(a, b) \lambda^{a-1} \exp(- b \lambda) \Bigr\} + \mathrm{const.} \\ &= \sum_{n=1}^N \Bigl\{ x_n \ln \lambda - \ln x_n! - \lambda \Bigr\} + \ln C_G(a, b) + (a - 1) \ln \lambda - b \lambda + \mathrm{const.} \\ &= \sum_{n=1}^N x_n \ln \lambda - N \lambda + (a - 1) \ln \lambda - b \lambda + \mathrm{const.} \\ &= \left( \sum_{n=1}^N x_n + a - 1 \right) \ln \lambda - (N + b) \lambda + \mathrm{const.} \tag{3.36} \end{align} $$

【途中式の途中式】

  1. 式(3.36)の下から2行目の式を用いている。
  2. 自然対数の性質より、$\ln x y = \ln x + \ln y$である。また、$\lambda$に影響しない$- \ln p(\mathbf{X})$を$\mathrm{const.}$とおく。
  3. 尤度と事前分布に、それぞれ具体的な確率分布の式を代入する。
  4. それぞれ対数をとる。$\ln x^a = a \ln x$、$\ln \frac{x}{y} = \ln x - \ln y$、$\ln \exp(x) = x$である。
  5. $\lambda$と関係のない$- \sum_{n=1}^N \ln x_n!$と$\ln C_G(a, b)$を$\mathrm{const.}$にまとめる。
  6. $\ln \lambda$と$\lambda$の項をそれぞれまとめて式を整理する。

となる。

 式(3.36)について

$$ \begin{align} \hat{a} &= \sum_{n=1}^N x_n + a \\ \hat{b} &= N + b \tag{3.38} \end{align} $$

とおき

$$ \ln p(\lambda | \mathbf{X}) = \hat{a} \ln \lambda - \hat{b} \lambda + \mathrm{const.} \tag{3.36} $$

さらに$\ln$を外し、$\mathrm{const.}$を正規化項に置き換える(正規化する)と

$$ p(\lambda | \mathbf{X}) = C_G(\hat{a}, \hat{b}) \lambda^{\hat{a}-1} \exp ( - \hat{b} \lambda ) = \mathrm{Gam}(\lambda | \hat{a}, \hat{b}) \tag{3.37} $$

事後分布は式の形状から、パラメータ$\hat{a},\ \hat{b}$を持つガンマ分布であることが分かる。
 また、式(3.38)が超パラメータ$\hat{a},\ \hat{b}$の計算式(更新式)である。

・予測分布の計算

 続いて、ポアソン分布に従う未観測のデータ$x_{*}$に対する予測分布を導出する。

 事前分布(観測データによる学習を行っていない分布)$p(\lambda)$を用いて、パラメータ$\lambda$を周辺化することで予測分布$p(x_{*})$となる。

$$ \begin{align} p(x_{*}) &= \int p(x_{*} | \lambda) p(\lambda) d\lambda \\ &= \int \mathrm{Poi}(x_{*} | \lambda) \mathrm{Gam}(\lambda | a, b) d\lambda \tag{3.39}\\ &= \int \frac{\lambda^{x_{*}}}{x_{*}!} \exp(- \lambda) C_G(a, b) \lambda^{a-1} \exp(- b \lambda) d\lambda \\ &= \frac{C_G(a, b)}{x_{*}!} \int \lambda^{x_{*}+a-1} \exp \Bigl( - (1 + b) \lambda \Bigr) d\lambda \tag{3.40} \end{align} $$

 積分部分に注目すると、パラメータ$x_{*} + a,\ 1 + b$を持つ正規化項のないガンマ分布の形をしている。
 これはガンマ分布の正規化項の逆数に変形できる(そもそもこの部分の逆数が正規化項である)。

$$ \int \lambda^{x_{*}+a-1} \exp \Bigl( - (1 + b) \lambda \Bigr) d\lambda = \frac{1}{C_G(x_{*} + a, 1 + b)} $$

 これを式(3.40)に代入すると、予測分布は

$$ p(x_{*}) = \frac{ C_G(a, b) }{ x_{*}! C_G(x_{*} + a, 1 + b) } \tag{3.42} $$

となる。
 さらにガンマ分布の正規化項(2.57)を用いると

$$ \begin{aligned} p(x_{*}) &= \frac{C_G(a, b)}{x_{*}!} \frac{ 1 }{ C_G(x_{*} + a, 1 + b) } \\ &= \frac{ b^a }{ x_{*}! \Gamma(a) } \frac{ \Gamma(x_{*} + a) }{ (1 + b)^{x_{*} + a} } \\ &= \frac{ \Gamma(x_{*} + a) }{ x_{*}! \Gamma(a) } \frac{ b^a }{ (1 + b)^{a} } \frac{1}{(1 + b)^{x_{*}}} \\ &= \frac{ \Gamma(x_{*} + a) }{ x_{*}! \Gamma(a) } \Bigl( \frac{ b }{ 1 + b } \Bigr)^{a} \Bigl( \frac{1}{1 + b} \Bigr)^{x_{*}} \\ &= \frac{ \Gamma(x_{*} + a) }{ x_{*}! \Gamma(a) } \Bigl( 1 - \frac{ 1 }{ 1 + b } \Bigr)^{a} \Bigl( \frac{1}{1 + b} \Bigr)^{x_{*}} \end{aligned} $$

となる。この式の形状の確率分布は負の二項分布である。

 この式について

$$ \begin{align} r &= a \\ p &= \frac{1}{1 + b} \tag{3.44} \end{align} $$

とおくと

$$ p(x_{*}) = \frac{ \Gamma(x_{*} + r) }{ x_{*}! \Gamma(r) } (1 - p)^r p^{x_{*}} = \mathrm{NB}(x_{*} | r, p) \tag{3.43} $$

予測分布は式の形状から、パラメータ$r,\ p$を持つ負の二項分布であることが分かる。

 予測分布の計算に事前分布$p(\lambda)$を用いることで、観測データによる学習を行っていない予測分布$p(x_{*})$(のパラメータ$r,\ p$)を求めた。事後分布$p(\lambda | \mathbf{X})$を用いると、同様の手順で観測データ$\mathbf{X}$によって学習した予測分布$p(x_{*} | \mathbf{X})$(のパラメータ)を求められる。

 そこで、$r,\ p$を構成する事前分布のパラメータ$a,\ b$について、事後分布のパラメータ(3.38)に置き換えたものをそれぞれ$\hat{r},\ \hat{p}$とおくと

$$ \begin{aligned} \hat{r} &= \hat{a} \\ &= \sum_{n=1}^N x_n + a \end{aligned} $$

また

$$ \begin{aligned} \hat{p} &= \frac{1}{1 + \hat{b}} \\ &= \frac{1}{N + 1 + b} \end{aligned} $$

となり、予測分布

$$ p(x_{*} | \mathbf{X}) = \frac{ \Gamma(x_{*} + \hat{r}) }{ x_{*}! \Gamma(\hat{r}) } (1 - \hat{p})^{\hat{r}} \hat{p}^{x_{*}} = \mathrm{NB}(x_{*} | \hat{r}, \hat{p}) $$

が得られる。
 また、上の式が予測分布のパラメータ$\hat{r},\ \hat{p}$の計算式(更新式)である。

参考文献

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

おわりに

 順調である。

【次節の内容】

www.anarchive-beta.com

2020/03/05:加筆修正しました。
2021/04/04:加筆修正しました。その際にRで実装編と記事を分割しました。と言ってもこの記事は表記を統一したくらいですかね。