DMPにおけるハイパーパラメータのサンプリングの仕方

ハイパーパラメータの決め方

Collapsed Gibbs samplingなどではパラメータは積分消去されることからハイパーパラメータが通常のパラメータの役割を果たすことが少なくありません。そういうわけで通常だと割と適当に「えいやっ!!」と決めてしまうようなハイパーパラメータをベイジアンな人たちは頑張って決める。LDAのときとかは経験ベイズっぽく最適化したり、DPMのハイパーパラメータのようなものはハイパーパラメータに事前分布(つまり、ハイパーハイパーパラメータが導入されるということである...)をかけて、ハイパーパラメータもサンプリングしてしまうのが普通らしい。どこまで事前分布を置くのが適切なのかは色々あるんだろうけど、とりあえずやり方だけは把握しておくことにする。

"Hyperparameter estimation in Dirichlet process mixture models"を参考にしながらメモる。Antoniak(1974)によると、DPMにおけるクラスタ数kに対して
$p(k | \alpha, n) \propto \alpha^k \frac{\Gamma(\alpha)}{\Gamma(\alpha + n)}$
と書けるらしい。これを使うと、DPMのハイパーパラメータ $\alpha$ の事後分布は
$p(\alpha | k) \propto p(\alpha) p(k | \alpha)$
とできるので、2つの項についてそれぞれ考えていくことにする。どうでもいいですけど $p(k | \alpha)$ ってちょっとすごって思ってしまいました。でも、PRMLのSection 3でも $p(\mbox{Model} | D)$ みたいなものが出てきたことあったなぁというのを思い出すとそんなに特別でもないのかもしれない(ようするによく分かっていない)。

まず、 $\alpha$ のprior $p(\alpha)$ について考えていこう。 $\alpha$ は0以上の値を取る確率変数として考えられるので、Gamma分布を事前分布として置くのが都合がよさそう( $\alpha \sim G(a, b)$ ということ)。ところで、 $p(k | \alpha, n) \propto \alpha^k \frac{\Gamma(\alpha)}{\Gamma(\alpha + n)}$ の $\frac{\Gamma(\alpha)}{\Gamma(\alpha + n)}$ のところは
$\frac{\Gamma(\alpha)}{\Gamma(\alpha + n)} = \frac{(\alpha + n) \beta(\alpha + 1, n)}{\alpha \Gamma(n)}$
と書ける。これはベータ関数とガンマ関数のよく知られた関係式 $\beta(x, y) = \frac{\Gamma(x) \Gamma(y)}{\Gamma(x + y)}$ と $\Gamma(n + 1) = n \Gamma(n)$ の関係式を使うと*1
$\frac{\Gamma(\alpha)}{\Gamma(\alpha + n)} = \frac{(\alpha + n)\Gamma(\alpha + 1)}{\alpha \Gamma(\alpha + 1 + n)} = \frac{(\alpha + n)\Gamma(\alpha + 1) \Gamma(n)}{\alpha \Gamma(n) \Gamma(\alpha + 1 + n)} = \frac{(\alpha + n) \beta(\alpha + 1, n)}{\alpha \Gamma(n)}$
とできるからである。

data augmentation(データ拡大)

最初見たときはこれの何がうれしいのかよく分からないが、ベータ関数を定義にしたがって積分の形で書いてみると
$p(\alpha | k) \propto p(\alpha) \alpha^{k-1} (\alpha + n) \beta(\alpha + 1, n) \propto p(\alpha) \alpha^{k-1} (\alpha + n) \int_0^1 x^\alpha (1 - x)^{n - 1} dx$
とできて、これは確率変数xってものがあったらそれを周辺化したものだと見なすことができる。xを周辺化しないで書くと
$p(\alpha, x | k) \propto p(\alpha) \alpha^{k-1} (\alpha + n) x^\alpha (1 - x)^{n - 1} = \left(\frac{1}{\Gamma(a)} b^a \alpha^{a - 1} \exp{(- b \alpha)}\right) \alpha^{k-1} (\alpha + n) x^\alpha (1 - x)^{n - 1} \propto \alpha^{a + k - 2} (\alpha + n) \exp{(- b \alpha)} x^\alpha$
とできる。また、 $y = \exp{(- b \alpha)} x^\alpha$ とおいて、両辺の対数を取ると $\log{y} = - b \alpha + \alpha \log{x}$ で、もっかい指数に戻すと $y = \exp{(- b \alpha + \alpha \log{x})} = \exp{(- \alpha(b - \log{x}))}$ とできる*2ので
$p(\alpha, x| k) \propto \alpha^{a + k - 2} (\alpha + n) \exp{(- \alpha(b - \log{x}))}$
とできる。

さて、なんでこんなことをやっていたかの種明かしをすると

$p(\alpha | k)$ からサンプリングしたい => しかし、ちょっと見た感じだとこの事後分布からのサンプリングは直接は厳しそう...
Gamma関数をばらして、ベータ関数に持っていき積分の形にすると、追加の変数xが周辺化されたものっぽいということが分かる
xとの同時分布 $p(\alpha, x| k)$ からは(後で見るように)容易にサンプリングできる
こやつを使えば元の欲しかった $p(\alpha | k)$ が手に入りそうだ!!!

ということである。こういうテクニックをデータ拡大アルゴリズム(data augmentation algorithm)と呼ぶらしく「そんなん知るかいな」と思ったが、実はPRMLの251ページにIPアルゴリズム(モンテカルロEMとも呼ばれる?)として紹介されていたorz(そして、この部分昔自分でゼミの時に紹介したことがあったところであったorz)。