正規分布の共役事前分布(正規ガンマ分布)

黒木玄

2017-11-28~2018-01-26, 2018-04-02, 2018-10-11

解説

正規分布とガンマ分布

正規分布の確率密度函数を

$$ p_N(x|\mu,\lambda) = \frac{1}{\sqrt{2\pi}}e^{-\lambda(x-\mu)^2/2}\lambda^{1/2} $$

と表わすことにする. この正規分布の平均と標準偏差はそれぞれ $\mu$, $\sigma=1/\sqrt{\lambda}$ である.

ガンマ分布の確率密度函数を

$$ p_\Gamma(\lambda|\alpha,\theta) = \frac{1}{\Gamma(\alpha)\theta^\alpha}e^{-\lambda/\theta}\lambda^{\alpha-1} $$

と表わす. このガンマ分布の平均と標準偏差はそれぞれ $\mu_\lambda=\alpha\theta$, $\sigma_\lambda=\sqrt{\alpha\theta^2}$ である.

共役事前分布(正規ガンマ分布)

正規分布の共役事前分布 $\varphi$ が次のように定義される:

$$ \begin{aligned} \varphi(\mu,\lambda|\mu_0,\lambda_0,\alpha,\theta) &= p_N(\mu|\mu_0, \lambda\lambda_0)p_\Gamma(\lambda|\alpha,\theta) \\ & = \sqrt{\frac{\lambda_0}{2\pi}} \frac{1}{\Gamma(\alpha)\theta^\alpha} \exp\left[-\frac{\lambda}{2}\left( \lambda_0(\mu-\mu_0)^2 + 2\theta^{-1} \right)\right] \lambda^{\alpha+1/2-1}. \end{aligned} $$

この確率密度函数で定義される確率分布を正規ガンマ分布と呼ぶことにする.

正規ガンマ分布 $\varphi(\mu,\lambda)=\varphi(\mu,\lambda|\mu_0,\lambda_0,\alpha,\theta)$ に関する $f(\mu,\lambda)$ の平均を

$$ E_{\varphi(\mu,\lambda)}[f(\mu,\lambda)] =\int_0^\infty d\lambda\int_{-\infty}^\infty d\mu\;f(\mu,\lambda)\varphi(\mu,\lambda) $$

と書くことにする.

$\mu_\lambda = \alpha\theta$ とおく. このとき

$$ \begin{aligned} & E_{\varphi(\mu,\lambda)}[\mu] = \mu_0, & & E_{\varphi(\mu,\lambda)}[(\mu-\mu_0)^2] = \begin{cases} \dfrac{1}{\lambda_0\mu_\lambda(1-\alpha^{-1})} & (\alpha > 1),\\ \infty & (0<\alpha\leqq 1), \end{cases} \\ & E_{\varphi(\mu,\lambda)}[\lambda] = \alpha\theta = \mu_\lambda, & & E_{\varphi(\mu,\lambda)}[(\lambda-\mu_\lambda)^2] = \alpha\theta^2 = \frac{\mu_\lambda^2}{\alpha}. \end{aligned} $$

$\mu_0$ と $\mu_\lambda$ を一定とするとき, $\mu,\lambda$ の分散を大きくするためには $\alpha$ と $\lambda_0$ を小さくすればよい.

共役事前分布のベイズ更新

事前分布 $\varphi(\mu,\lambda)$ と サンプル $X_1,\ldots,X_n$ に関する逆温度 $\beta$ の事後分布 $\tilde\varphi(\mu,\lambda)$ は次の式で定義される:

$$ \tilde\varphi(\mu,\lambda) = \frac {\left(\prod_{i=1}^n p_N(X_i|\mu,\lambda)\right)^\beta\varphi(\mu,\lambda)} {\int_0^\infty d\lambda\int_{-\infty}^\infty d\mu\,\left(\prod_{i=1}^n p_N(X_i|\mu,\lambda)\right)^\beta\varphi(\mu,\lambda)}. $$

$\varphi(\mu,\lambda) = \varphi(\mu,\lambda|\mu_0,\lambda_0,\alpha,\theta)$ のとき

$$ \tilde\varphi(\mu,\lambda) = \varphi(\mu,\lambda|\tilde\mu_0, \tilde\lambda_0, \tilde\alpha, \tilde\theta) $$

となる. ただし, $\tilde\mu_0, \tilde\lambda_0, \tilde\alpha, \tilde\theta$ は, サンプルの平均と分散を

$$ \bar X = \frac{1}{n}\sum_{i=1}^n X_i, \quad V(X) = \frac{1}{n}\sum_{i=1}^n (X_i - \bar X)^2. $$

書くことにすると, 以下のように定義される:

$$ \begin{aligned} & \tilde\mu_0 = \frac{\lambda_0\mu_0 + \beta n \bar X}{\lambda_0 + \beta n}, %\quad & & \tilde\lambda_0 = \lambda_0 + \beta n, \quad \\ & \tilde\alpha = \alpha + \frac{1}{2}\beta n, %\quad & & {\tilde\theta}^{-1} = \theta^{-1}\left[ 1 + \frac{\theta}{2}\left( \frac{\beta n\lambda_0}{\lambda_0 + \beta n}\left(\bar X - \mu_0\right)^2 + \beta n V(X) \right) \right]. \end{aligned} $$

これを共役事前分布のベイズ更新と呼ぶ. 共役分布のベイズ更新はもとのパラメーターと $\beta n$, $\bar X$, $V(X)$ で表わされる. サンプルのサイズと逆温度の積, サンプルの平均, サンプルの分散の情報があればベイズ更新可能である.

予測分布

$\varphi(\mu,\lambda) = \varphi(\mu,\lambda|\mu_0,\lambda_0,\alpha,\theta)$ に関する予測分布 $p^*(x)$ が

$$ p^*(x) = E_{\varphi(\mu,\lambda)}[p_N(x|\mu,\lambda)] = \int_0^\infty d\lambda\int_{-\infty}^\infty p_N(x|\mu,\lambda)\varphi(\mu,\lambda) $$

と定義される. 予測分布は次の形になる:

$$ p^*(x) = \sqrt{\frac{1}{2\pi}\frac{\theta\lambda_0}{\lambda_0+1}}\; \frac{\Gamma(\alpha+1/2)}{\Gamma(\alpha)} \left[ 1 + \frac{1}{2}\frac{\theta\lambda_0}{\lambda_0+1}(x-\mu_0)^2 \right]^{-(\alpha+1/2)}. $$

この分布は自由度 $\nu = 2\alpha$ のt分布を $\mu_0$ が中心になるように平行移動し, $\displaystyle \rho = \sqrt{\frac{\lambda_0+1}{\theta\lambda_0}}$ 倍でスケール変換して得らえる確率分布である. この予測分布に従う確率変数 $X^*$ は自由度 $\nu$ のt分布に従う確率変数 $T(\nu)$ によって

$$ X^* = \mu_0 + \rho\,T(\nu), \quad \rho = \sqrt{\frac{\lambda_0+1}{\theta\lambda_0}}, \quad \nu=2\alpha $$

と表わされる. 自由度 $\nu$ を大きくするとt分布は正規分布に近付く.

サンプルサイズ $n$ を大きくすると, 予測分布は平均 $\bar X$, 分散 $V(X)$ の正規分布に漸近して行く.

分配函数

$\varphi(\mu,\lambda) = \varphi(\mu,\lambda|\mu_0,\lambda_0,\alpha,\theta)$ に関する逆温度 $\beta$ の分配函数が次のように定義される:

$$ Z^\beta(X_1,\ldots,X_n|\mu_0,\lambda_0,\alpha,\theta) = \int_0^\infty d\lambda\int_{-\infty}^\infty d\mu\; \left(\prod_{i=1}^n p_N(X_i|\mu,\lambda)\right)^\beta\varphi(\mu,\lambda). $$

これは次のように表わされる:

$$ Z^\beta(X_1,\ldots,X_n|\mu_0,\lambda_0,\alpha,\theta) = \frac {z(\tilde\mu_0, \tilde\lambda_0, \tilde\alpha, \tilde\theta)} {(2\pi)^{\beta n/2}\; z(\mu_0,\lambda_0,\alpha,\theta)}. $$

ここで,

$$ z(\mu_0,\lambda_0,\alpha,\theta) = \sqrt{\frac{2\pi}{\lambda_0}} \Gamma(\alpha)\theta^\alpha. $$

この公式を用いれば対数分配函数 $\log Z^\beta(\mu_0,\lambda_0,\alpha,\theta)$ に関する公式も得られる:

$$ \log Z^\beta(X_1,\ldots,X_n|\mu_0,\lambda_0,\alpha,\theta) = \log z(\tilde\mu_0, \tilde\lambda_0, \tilde\alpha, \tilde\theta) - \log z(\mu_0,\lambda_0,\alpha,\theta) -\frac{\beta n}{2}\log(2\pi). $$

ここで

$$ \log z(\mu_0,\lambda_0,\alpha,\theta) = \frac{1}{2}\left(\log(2\pi) - \log\lambda_0\right) + \log\Gamma(\alpha) + \alpha\log\theta. $$

LOOCV (1個抜き出し交差検証, leave-one-out cross-validation)

このノートでは情報量規準のスケールとして, AICやBICで採用されている伝統的なスケールを採用する. そのスケールのもとで LOOCV は次のように定義される:

$$ \mathrm{LOOCV} = - 2\sum_{i=1}^n \log \left(\text{逆温度 $\beta$ におけるサンプル $X_1,\ldots,\widehat X_i,\ldots,X_n$ の事後分布に関する $p_N(X_i|\mu,\lambda)$} の期待値\right). $$

ここで $\widehat X_i$ は $X_i$ を除くという意味である. $X_i$ を除いたサンプルに関する事後分布に関する期待値は除かない場合の期待値で書ける. 除かない場合の通常の逆温度は $\beta$ の事後分布に関する期待値を $E_{\tilde\varphi(\mu,\lambda)}[f(\mu,\lambda)]$ と書くことにすると,

$$ \mathrm{LOOCV} = - 2\sum_{i=1}^n \log\frac {E_{\tilde\varphi(\mu,\lambda)}\left[p_N(X_i|\mu,\lambda)^{1-\beta}\right]} {E_{\tilde\varphi(\mu,\lambda)}\left[p_N(X_i|\mu,\lambda)^{-\beta}\right]}. $$

分配函数の定義より,

$$ Z^\kappa(x|\tilde\mu_0, \tilde\lambda_0, \tilde\alpha, \tilde\theta) = \int_0^\infty d\lambda\int_{-\infty}^\infty d\mu\; p_N(x|\mu,\lambda)^\kappa\;\varphi(\mu,\lambda|\tilde\mu_0, \tilde\lambda_0, \tilde\alpha, \tilde\theta). $$

であることに注意すれば, LOOCVは事後分布を事前分布とする分配函数によって次のように表わされることがわかる:

$$ \mathrm{LOOCV} = - 2\sum_{i=1}^n \log\frac {Z^{1-\beta}(X_i|\tilde\mu_0, \tilde\lambda_0, \tilde\alpha, \tilde\theta)} {Z^{-\beta}(X_i|\tilde\mu_0, \tilde\lambda_0, \tilde\alpha, \tilde\theta)}. $$

この公式に前節の公式を適用すればLOOCVを計算できる.

ベイズ自由エネルギー

AICやBICの伝統的なスケールのもとでの逆温度 $\beta$ のベイズ自由エネルギーは次のように定義される:

$$ \mathrm{BayesianFreeEnergy} = -\frac{2}{\beta}\log Z^\beta(X_1,\ldots,X_n|\tilde\mu_0, \tilde\lambda_0, \tilde\alpha, \tilde\theta). $$

対数尤度函数およびその二乗の期待値

WBICやWAICの公式を作るためには $\varphi(\mu,\lambda) = \varphi(\mu,\lambda|\mu_0,\lambda_0,\alpha,\theta)$ に関する $\log p_N(x|\mu,\lambda)$ およびその二乗の平均の公式を作っておく必要がある. その結果は以下の通り. 函数 $f(\mu,\lambda)$ に対して

$$ E_{\varphi(\mu,\lambda)}[f(\mu,\lambda)] = \int_0^\infty d\lambda \int_{-\infty}^\infty d\mu\; f(\mu,\lambda)\;\varphi(\mu,\lambda|\mu_0,\lambda_0,\alpha,\theta) $$

と書くことにする. このとき,

$$ -2 E_{\varphi(\mu,\lambda)}[\log p_N(x|\mu,\lambda)] = \log(2\pi) + \alpha\theta(x-\mu_0)^2 + \frac{1}{\lambda_0} - (\psi(\alpha)+\log\theta). $$

ここで $\psi(\alpha)$ はガンマ函数の対数微分であり, digamma函数と呼ばれる. さらに

$$ \begin{aligned} 4 E_{\varphi(\mu,\lambda)}\left[(\log p_N(x|\mu,\lambda))^2\right] &= (\log(2\pi))^2 \\ & + (\alpha+1)\alpha\theta^2(x-\mu_0)^4 + \frac{6\alpha\theta}{\lambda_0}(x-\mu_0)^2 + \frac{3}{\lambda_0^2} \\ & + \psi'(\alpha) + (\psi(\alpha)+\log\theta)^2 \\ & + 2\log(2\pi)\left(\alpha\theta(x-\mu_0)^2+\frac{1}{\lambda_0}\right) \\ & - 2\log(2\pi)(\psi(\alpha)+\log\theta) \\ & - 2\left( (x-\mu_0)^2 (\theta+\alpha\theta\;(\psi(\alpha)+\log\theta)) + \frac{\psi(\alpha)+\log\theta}{\lambda_0} \right). \end{aligned} $$

これらの公式を使えば, WBICとWAICを正確に計算するための公式が得られる.

このノートブックでは実際にその公式を使ってWBICとWAICを計算できるようにしてある.

訂正(2018-04-02) $4 E_{\varphi(\mu,\lambda)}\left[(\log p_N(x|\mu,\lambda))^2\right]$ の公式の最後の項を修正した(なぜか $\alpha\theta/\lambda=0$ になっていた). WBICを計算するためのコードの方には誤りは無かった.

WBIC

対数分配函数の $\beta$ による微分は次のように書ける:

$$ Z(\beta) = Z^\beta(X_1,\ldots,X_n|\mu_0,\lambda_0,\alpha,\theta), \quad \tilde\varphi(\mu,\lambda) = \varphi(\mu,\lambda|\tilde\mu_0, \tilde\lambda_0, \tilde\alpha, \tilde\theta) $$

とおくと, $Z(0)=1$ でかつ

$$ -\frac{\partial}{\partial\beta}\log Z(\beta)= \int_0^\infty d\lambda\int_{-\infty}^\infty d\mu\; H(\lambda,\mu)\;\tilde\varphi(\mu,\lambda) = E_{\tilde\varphi(\mu,\lambda)}[H(\mu,\lambda)]. $$

ここで

$$ H(w) = - \sum_{i=1}^n \log p(X_i|\mu,\lambda). $$

$\tilde\varphi(\mu,\lambda)$ で省略されているパラメーターの中に逆温度 $\beta$ が含まれていることに注意せよ.

ゆえに,

$$ -\log Z(1) = \int_0^1 d\beta\;E_{\tilde\varphi(\mu,\lambda)}[H(\mu,\lambda)] $$

WBICの理論によれば上の $\log Z(1)$ は $\beta=1/\log n$ における $E_{\tilde\varphi(\mu,\lambda)}[H(\mu,\lambda)]$ の値で近似される. AICとBICの伝統的スケールでWBICは次のように定義される:

$$ \mathrm{WBIC} = 2\left.E_{\tilde\varphi(\mu,\lambda)}[H(\mu,\lambda)]\right|_{\beta = 1/\log n} = 2\left.\sum_{i=1}^n E_{\tilde\varphi(\mu,\lambda)}[-\log p_N(X_i|\mu,\lambda)]\right|_{\beta = 1/\log n}. $$

これは前節の公式を使えば計算できる. WBICは $\beta=1$ での BayesianFreeEnergy の推定値になっている.

WBICについては

http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/wbic2012.html

を見よ. そこでは, このノートにおける $H(w)=H(\mu,\lambda)$ は $E^\beta_w[nL_n(w)]$ と書かれている.

$$ \mathrm{WBIC} = 2\left(\text{逆温度 $\beta=\frac{1}{\log n}$ での事後分布に関する対数尤度函数の $-1$ 倍の期待値}\right). $$

WAIC

AICやBICの伝統的なスケールにおけるWAICは以下のように定義される:

$$ \mathrm{WAIC} = 2(T + \beta V). $$

ここで

$$ \begin{aligned} & T = -\sum_{i=1}^n \log\tilde p^*_N(X_i), \\ & \tilde p^*(x) = E_{\tilde\varphi(\mu,\lambda)}[p_N(X_i|\mu,\lambda)], \\ & V = \sum_{i=1}^n V_{\tilde\varphi(\mu,\lambda)}\left[\log p_N(X_i|\mu,\lambda)\right], \\ & V_{\tilde\varphi(\mu,\lambda)}\left[\log p_N(X_i|\mu,\lambda)\right] = E_{\tilde\varphi(\mu,\lambda)}\left[(\log p_N(X_i|\mu,\lambda))^2\right] - E_{\tilde\varphi(\mu,\lambda)}\left[\log p_N(X_i|\mu,\lambda)\right]^2. \end{aligned} $$

これらは前々節の公式を使えば計算可能である.

WAICについては次のサイトを参照せよ:

http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/waic2011.html

$$ \begin{aligned} \mathrm{WAIC} & = -2 \sum_{i=1}^n\log\left(\text{逆温度 $\beta$ での事後分布に関する $p(X_i|w)$ の期待値}\right) \\ & + 2\beta \sum_{i=1}^n\left(\text{逆温度 $\beta$ での事後分布に関する $\log p(X_i|w)$ の分散}\right). \end{aligned} $$

improved WBIC?

ベイズ情報量規準及びその発展 ~概説編~ 2016年4月19日

には improved WBIC に関する研究が賞をもらったと書いてある:

https://www.albert2005.co.jp/release/archives/201604/19_110050.html

しかし, Googleなどで検索しても論文を見付けることができなかった(2017-11-30).

さらに検索を続けると

http://onlinelibrary.wiley.com/doi/10.1111/rssb.12187/full

のdiscussionに定義式が書いてあるように見えた. その部分を次のセルに引用しておく.

In [1]:
using Base64
open("drton2017_01.jpg") do f
    base64 = base64encode(f)
    display("text/html", """<img src="data:image/jpg;base64,$base64">""")
end

筆者なりに計算した結果は以下の通り. 以下ではAICやBICの伝統的なスケールをWAICやWBICでも採用する. WBICは次の論文の2倍がこのノートでの定義になる(WAICは $2n$ 倍).

WBICの元論文は

[WBIC] http://www.jmlr.org/papers/volume14/watanabe13a/watanabe13a.pdf

である. このノートのWBICの定義はこの論文の定義の2倍である.

WAICについては

[渡辺] 渡辺澄夫『ベイズ統計の理論と方法』

の解説が詳しい. このノートのWAICの定義はこの本の定義の $2n$ 倍である.

[WBIC]のpp.876--877を次のセルに引用しておく.

In [2]:
open("WBIC_paper_01.jpg") do f
    base64 = base64encode(f)
    display("text/html", """<img src="data:image/jpg;base64,$base64">""")
end