比例記号を活用した事後分布の計算テク #統計検定1級

この記事について

ベイズ統計学や統計検定１級で頻出であり，また煩雑になりがちな事後分布の計算を，比例記号を用いて楽に計算しようというものです．有名なテクニックですが，事後分布を手で計算したい人類が少ないのか，ググってもあまり多くの記事がヒットしないので，まとめてみました．

目的

以下の問題を考えます．いわゆる，ベータ二項分布というやつです．

$X \sim Bin(n, p)$ とし，$p$ に事前分布 $Beta(\alpha, \beta)$ を仮定する．ここで、$\alpha, \beta > 0$ は既知とする．
このとき，$p$ の事後分布 $\pi (p\mid x, \alpha, \beta)$ を求めよ．

事前分布と事後分布の関係，$\pi (p\mid x, \alpha, \beta) = f(x\mid p) \cdot \pi(p\mid \alpha, \beta) / f_{\pi}(x\mid \alpha, \beta)$ にすべて代入してゴリゴリ計算しても求めることはできますが，煩雑な計算となってしまいます．特に，統計検定１級では，時間が足りないことも多く，また，計算ミスをしてしまうと雪崩失点してしまう可能性もあります．
そこで，この記事では，これの計算を簡単にする方法を紹介したいと思います．

重要な考え方

分布を表す密度関数には，「核」となる部分（変数が含まれている部分）と，「辻褄合わせ」の部分（正規化定数）があります．例えば，ベータ分布 $Beta(\alpha, \beta)$ の分布は，以下のように表すことができます．

$$f(x) = \frac{1}{B(\alpha, \beta)} x^{\alpha - 1} (1-x)^{\beta -1}, 0 < x < 1$$

ここで，核となる部分は $x^{\alpha - 1} (1-x)^{\beta -1}$ の部分です．この部分がこの分布をベータ分布たらしめています．
一方で，$\frac{1}{B(\alpha, \beta)}$ の部分は正規化定数といい，$\int_{-\infty}^{\infty} f(x) dx = 1$ を満たすために存在している，いわば辻褄合わせの部分です．もっというと，この正規化定数は核の部分を用いて $\frac{1}{\int_{0}^{1} x^{\alpha - 1} (1-x)^{\beta -1} dx}$ と計算することができるので，核の部分だけわかればいつでも導出できます．
なにが言いたいかというと，分布を求める際には，変数に関係ない係数の部分は一旦無視をして，核の部分を求めるだけで十分です．もし正規化定数が必要であれば，核の部分を計算した後に導出すればよいです．この核の部分に注目する考え方が，これから述べるテクニックの根幹となっています．

比例記号の導入

今後のために，ここで比例記号 $\propto$ を導入しておきます．すなわち，$y$ が $x$ に比例するとき，$y \propto x$ と書きます．これがあるとなにが嬉しいかというと，ざっくりいえば，定数倍を無視した式変形を行うことができます．例えば，ベータ分布の例を出すと，

$$\begin{aligned} f(x) &= \frac{1}{B(\alpha, \beta)} x^{\alpha - 1} (1-x)^{\beta -1} \\ &\propto x^{\alpha - 1} (1-x)^{\beta -1} \end{aligned}$$

と書くことができます．この例からわかるように，比例記号 $\propto$ は，前節で述べた核の部分に注目する考え方と非常に相性が良いです．

比例記号を用いた事前分布の計算

準備は整ったので，冒頭の問題を解いていきましょう．今回注目している変数が，$p$ であることに注意して，比例記号 $\propto$ を用いて変形していきます．(スマホ等で，数式の各行とその式番号が被ってしまう場合は，画面を横向きにして試してみてください．)

$$\begin{align} \pi (p \mid x, \alpha, \beta) &= f(x\mid p) \cdot \pi(p\mid \alpha, \beta) / f_{\pi}(x\mid \alpha, \beta)\\ &\propto f(x\mid p) \cdot \pi(p\mid \alpha, \beta)\\ &\propto p^x(1-p)^{n-x} \cdot p^{\alpha - 1} (1-p)^{\beta -1}\\ &= p^{x + \alpha -1}(1-p)^{n-x+\beta -1}\\ \end{align}$$

式変形の解説

順に説明をしていきます．なお，説明中の「無視する」という表現は，「式変形中では，比例記号を用いて表記を省略し，考えないものとする」くらいの意味です．

1. (1) では，事前分布と事後分布の関係を用いました．
2. (2) では，$1/f_{\pi}(x\mid \alpha, \beta)$ を $p$ に依存しない定数として，無視しました．
3. (3) では，それぞれ，$x$の分布（二項分布）と $p$ の事前分布（ベータ分布）を代入しています．ここで，$p$ に依存しない正規化定数は無視しました．
4. 指数法則を用いて整理しました．

上記の式変形から，$\pi (p \mid x, \alpha, \beta) \propto p^{x + \alpha -1}(1-p)^{n-x+\beta -1}$ がわかりました．
ここで，$p^{x + \alpha -1}(1-p)^{n-x+\beta -1}$ は，$Beta(x + \alpha, n-x+\beta)$ の核の部分と一致するので，事後分布は $Beta(x + \alpha, n-x+\beta)$ であるとわかります．なお，正規化定数は，$Beta(\alpha, \beta)$ の正規化定数が $\frac{1}{B(\alpha, \beta)}$ であることから，$Beta(x + \alpha, n-x+\beta)$ の正規化定数は $\frac{1}{B(x + \alpha, n - x + \beta)}$ と簡単に求めることができます（もちろん，積分計算して直接求めてもよいです）．
このように，比例記号を用いて，余計な定数部分の計算を端折ることで，圧倒的に楽に計算することができます．

ちなみに，少し脱線しますが，本問題のように，母数 $p$ が定める確率分布が二項分布であるとき，$p$ の事前分布をベータ分布としてあげれば，事後分布もベータ分布となります．このように，事前分布と事後分布が同じ分布族に入るような事前分布を，共役事前分布 といいます．
共役事前分布であると嬉しいことがたくさんあります．例えば，本問題のように計算が容易であったり，また，この事後分布を次の事前分布として用いたとき，その事後分布もまた，同じ分布族であるので，ベイズの更新を容易に行うことができます．

他の例

有名な共役事前分布の例を用いて，比例記号に慣れていきましょう．

ガンマポアソン分布

次の問題を考えてみましょう．

$X \sim Po(\lambda)$ とし，$\lambda$ に事前分布 $Ga(\alpha, \beta)$ を仮定する．ここで、$\alpha, \beta$ は既知とする．
このとき，$\lambda$ の事後分布 $\pi (\lambda \mid x, \alpha, \beta)$ を求めよ．

ここに，ポアソン分布 $Po(\lambda)$，ガンマ分布 $Ga(\alpha, \beta)$ の分布を示します．

$$\begin{aligned} f(x\mid \lambda) &= \frac{\lambda^x}{x!}e^{-\lambda} \\ f(x\mid \alpha, \beta) &= \frac{1}{\Gamma (\alpha)\beta ^ {\alpha}} x^{\alpha -1} e^{-\frac{x}{\beta}} \end{aligned}$$

これらを用いて，事後分布を計算していきます．

$$\begin{align} \pi (\lambda \mid x, \alpha, \beta) &= f(x\mid \lambda) \cdot \pi(\lambda\mid \alpha, \beta) / f_{\pi}(x\mid \alpha, \beta)\\ &\propto f(x\mid \lambda) \cdot \pi(\lambda\mid \alpha, \beta)\\ &\propto \lambda^x e^{-\lambda} \cdot \lambda ^{\alpha -1} e^{-\frac{\lambda}{\beta}}\\ &= \lambda^{x + \alpha -1} e^{-\left(1 + \frac{1}{\beta}\right)\lambda} \end{align}$$

上記の式変形から，$\pi (\lambda \mid x, \alpha, \beta) \propto \lambda^{x + \alpha -1} e^{-\left(1 + \frac{1}{\beta}\right)\lambda}$ がわかりました．
ここで，$\lambda^{x + \alpha -1} e^{-\left(1 + \frac{1}{\beta}\right)\lambda}$ は，$Ga\left(x + \alpha, \left(1+\frac{1}{\beta}\right)^{-1}\right)$ の核の部分と一致するので，事後分布は $Ga\left(x + \alpha, \left(1+\frac{1}{\beta}\right)^{-1}\right)$ であるとわかります．
ベータ二項分布のときの関係のように，母数 $\lambda$ が定める分布がポアソン分布であるとき，ガンマ分布は共役事前分布になります．

正規分布

もう一問，次の問題を考えてみましょう．みんな大好き正規分布です．

$X_1, X_2, \dots X_n \overset{\text\small\textrm{i.i.d.}}{\sim} N(\mu, \sigma^2)$ とし，$\mu$ に事前分布 $N(\xi, \tau^2)$ を仮定する．ここで，$\sigma, \xi, \tau$ は既知であるとする．
このとき，$\mu$ の事後分布 $\pi (\mu \mid \bm{x}, \sigma, \xi, \tau)$ を求めよ．

ここに，正規分布の分布を示します．

$$f(x\mid \mu, \sigma) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$

これを用いて，事前分布を計算していきます．少し大変です．

$$\begin{align} \pi (\mu \mid \bm{x}, \sigma, \xi, \tau) &= f(\bm{x}\mid \mu, \sigma) \cdot \pi(\mu\mid \xi, \tau) / f_{\pi}(\bm{x}\mid \xi, \tau)\\ &\propto f(\bm{x}\mid \mu, \sigma) \cdot \pi(\mu\mid \xi, \tau)\\ &\propto \left(\prod_{i=1}^{n} \exp\left(-\frac{(x_i-\mu)^2}{2\sigma^2}\right)\right)\cdot \exp\left(-\frac{(\mu-\xi)^2}{2\tau^2}\right)\\ &= \exp\left(-\sum_{i=1}^{n}\frac{(x_i-\mu)^2}{2\sigma^2} - \frac{(\mu-\xi)^2}{2\tau^2}\right)\\ &\propto \exp\left(-\sum_{i=1}^{n}\frac{-2x_i\mu + \mu^2}{2\sigma^2} - \frac{-2\xi\mu + \mu^2}{2\tau^2}\right)\\ &= \exp\left(-\left(\frac{n}{2\sigma^2} + \frac{1}{2\tau^2}\right)\mu^2 + \left(\frac{n\bar{x}}{\sigma^2} + \frac{\xi}{\tau^2}\right)\mu \right)\\ &\propto \exp\left(-\frac{\left(\mu - \frac{\frac{n\bar{x}}{\sigma^2} + \frac{\xi}{\tau^2}}{\frac{n}{\sigma^2} + \frac{1}{\tau^2}} \right)^2}{2\cdot\frac{1}{\left(\frac{n}{\sigma^2} + \frac{1}{\tau^2}\right)}} \right) \end{align}$$

式変形の補足

煩雑な変形の部分を補足します．

1. (11) では，$\bm{x}$ が $n$ 個の独立な確率変数からの標本であることに注意して，尤度を計算しました．ここで，$\mu$ に依存しない正規化定数は無視しました．
2. (12) では，指数法則 $e^a\cdot e^b = e^{a + b}$ を用いました．
3. (13) では，指数部を展開しました．ここで，$\mu$ に依存しない部分は無視しました（指数部内の計算であるため，係数ではなく，展開時の項の方に $\mu$ を含まない場合を無視することに注意してください）．
4. (14) では，$\mu$ について整理しました．ここで，$\bar{x}$ は $x$ の平均を表し，$\bar{x} = \frac{1}{n} \sum_{i=0}^{n}x_i$ を用いました．
5. (15) では，$\mu$ について平方完成しました．ここで，$\mu$ に依存しない部分は無視しました（上記4.と同様に，指数部内の計算であるため，係数ではなく，展開時の項の方に $\mu$ を含まない場合を無視することに注意してください）．

上記の式変形から，$\pi (\mu \mid \bm{x}, \sigma, \xi, \tau) \propto \exp\left(-\frac{\left(\mu - \frac{\frac{n\bar{x}}{\sigma^2} + \frac{\xi}{\tau^2}}{\frac{n}{\sigma^2} + \frac{1}{\tau^2}} \right)^2}{2\cdot\frac{1}{\left(\frac{n}{\sigma^2} + \frac{1}{\tau^2}\right)}} \right)$ がわかりました．
ここで，$\exp\left(-\frac{\left(\mu - \frac{\frac{n\bar{x}}{\sigma^2} + \frac{\xi}{\tau^2}}{\frac{n}{\sigma^2} + \frac{1}{\tau^2}} \right)^2}{2\cdot\frac{1}{\left(\frac{n}{\sigma^2} + \frac{1}{\tau^2}\right)}} \right)$ は，正規分布 $N\left(\frac{\frac{n\bar{x}}{\sigma^2} + \frac{\xi}{\tau^2}}{\frac{n}{\sigma^2} + \frac{1}{\tau^2}}, \frac{1}{\left(\frac{n}{\sigma^2} + \frac{1}{\tau^2}\right)}\right)$ の核の部分と一致するので，事後分布は $N\left(\frac{\frac{n\bar{x}}{\sigma^2} + \frac{\xi}{\tau^2}}{\frac{n}{\sigma^2} + \frac{1}{\tau^2}}, \frac{1}{\left(\frac{n}{\sigma^2} + \frac{1}{\tau^2}\right)}\right)$ であることがわかります．
ここまで見てきた例と同様に，母数 $\mu$ が定める分布が正規分布であるとき，正規分布自身は共役事前分布になります．
とても大変な計算でしたが，比例記号を使わない場合は，これに加えて正規化定数の計算も同時に実行しなければならないので，修行のようになってしまいます．比例記号を用いることで，$\mu$ に依存しない部分の計算を大幅に削減できるので，計算スピードの向上と計算ミスの削減に大きく貢献することがわかると思います．

おわりに

比例記号 $\propto$ の難点をあげるとしたら，手書きでうまく書くのが難しいことです．
もし，間違っている箇所があれば，コメント，もしくは，XのDM等でご連絡いただけるとありがたいです．

参考文献

現代数理統計学の基礎，久保川達也．https://www.kyoritsu-pub.co.jp/book/b10003681.html