1 序論
Exercise 1.7 \(\alpha, \beta\)は正の実数とする.観測\(x\)は\(\mathcal{E}(\theta)\)に従い,\(\theta\)の事前分布を\(\mathcal{G}(\alpha,\beta)\)とする.
(a). このとき\(\theta\)の事後分布を求めよ.
(b). 事後平均を求めよ.
(c). 事後予測分布を求めよ.
Solution. (a). 事後密度関数は \[\begin{align*} p(\theta|x)&\propto \theta e^{-\theta x} \theta^{\alpha-1}e^{-\beta\theta}\\ &\propto \theta^\alpha e^{-(x+\beta)\theta} \end{align*}\] となる.従って,事後密度関数の形から,事後分布は \(\mathcal{G}(\alpha+1, x+\beta)\).
(b). 練習問題1.5より \[\frac{\alpha+1}{x+\beta}.\]
(c). 正の実数\(x^*\)を将来の観測とすると,事後予測分布の確率密度関数は \[\begin{align*} p(x^*|x)&=\int_0^\infty\theta e^{-\theta x^*}\frac{(x+\beta)^{\alpha+1}\theta^\alpha}{\Gamma(\alpha+1)} e^{-(x+\beta)\theta}\mathrm{d}\theta\\ &=\int_0^\infty\frac{(x+\beta)^{\alpha+1}\theta^{\alpha+1}}{\Gamma(\alpha+1)} e^{-(x+x^*\beta)\theta}\mathrm{d}\theta\\ &=\frac{\Gamma(\alpha+2)}{(x+x^*+\beta)^{\alpha+2}}\frac{(x+\beta)^{\alpha+1}}{\Gamma(\alpha+1)}\\ &=\frac{(\alpha+1)(x+\beta)^{\alpha+1}}{(x+x^*+\beta)^{\alpha+2}}. \end{align*}\]Exercise 1.8 観測\(x\)は\(\mathcal{N}(\theta_1,\theta_2^{-1})\)に従い,\(\theta_1\)の事前分布を\(\mathcal{N}(0,\theta_2^{-1})\)とし,さらに\(\theta_2\)の事前分布を\(\mathcal{E}(1)\)とする.
(a). \(\theta_1\)と\(x\)で条件づけた\(\theta_2\)の条件つき分布を求めよ.
(b). \(\theta_1\)の周辺事後分布を求めよ.
Solution. (a). 尤度は \[p(x|\theta_1,\theta_2)=\sqrt{\frac{\theta_2}{2\pi}}\exp\left(-\frac{(x-\theta_1)^2}{2}\theta_2\right)\] と書ける.したがって事後分布の確率密度関数は \[\begin{align*} p(\theta_1,\theta_2|x) &\propto p(x|\theta_1,\theta_2)p(\theta_1|\theta_2)p(\theta_2)\\ &= \sqrt{\frac{\theta_2}{2\pi}}\exp\left(-\frac{(x-\theta_1)^2}{2}\theta_2\right)\sqrt{\frac{\theta_2}{2\pi}}\exp\left(-\frac{\theta_1^2}{2}\theta_2\right)\exp(-\theta_2)\\ &\propto \theta_2\exp\left(-\left(\frac{(x-\theta_1)^2}{2}+\frac{\theta_1^2}{2}+1\right)\theta_2\right) \end{align*}\] となる.すると\(p(\theta_2|x, \theta_1,)\propto p(\theta_1,\theta_2|x)\)だから\(\theta_1\)と\(x\)で条件づけた\(\theta_2\)の条件つき分布は ガンマ分布とパラメータを見比べると, \[\mathcal{G}_a\left(2,\frac{(x-\theta_1)^2}{2}+\frac{\theta_1^2}{2}+1\right)\] であることがわかる.
(b). ガンマ関数と見比べることで, \[\begin{align*} p(\theta_1|x)&=\int_0^\infty p(\theta_1,\theta_2|x)\mathrm{d}\theta_2\\ &=\left(\frac{(x-\theta_1)^2}{2}+\frac{\theta_1^2}{2}+1\right)^{-2}\int_0^\infty \left(\frac{(x-\theta_1)^2}{2}+\frac{\theta_1^2}{2}+1\right)^2\\ &\quad\theta_2\exp\left(-\left(\frac{(x-\theta_1)^2}{2}+\frac{\theta_1^2}{2}+1\right)\theta_2\right)\mathrm{d}\theta_2\\ &=\left(\frac{(x-\theta_1)^2}{2}+\frac{\theta_1^2}{2}+1\right)^{-2}\Gamma(2) \end{align*}\] となる.だから\(t\)分布のパラメータと見比べれば,これは \[ \mathcal{T}_3\left(x, \frac{2}{3}\left(\frac{\theta_1^2}{2}+1\right)^{-1}\right) \] の確率密度関数である.したがって,周辺事後分布は上記の\(t\)分布である.Remark. ディリクレ分布は和が\(1\)になる正の数列\(\theta_1,\ldots, \theta_K\)に値を取る確率分布である.とくに\(K=2\)としたときがベータ分布である.和をとると\(1\)になる列は多項分布の確率関数になるから,ディリクレ分布は確率分布の確率分布とも捉えられる.うえの練習問題はこの事実を使った.
ベータ分布は\(K=2\)であり,\(\theta_1+\theta_2=1\)となるので,\(\theta_1\)さえわかれば\(\theta_2\)もわかるので確率密度関数を描画するのもたやすい.だから,ベータ分布のパラメータがベータ分布にどのように影響するかをみることも容易である.それに比べると,\(K\ge 3\)のディリクレ分布の形状の把握は難しい.ここでは\(K=3\)の場合を考えてみよう.じつは \[X_1\sim\mathcal{G}(\alpha_1,1),\ldots, X_K\sim\mathcal{G}(\alpha_K,1)\] であり,独立であるとき, \[\theta=(\theta_1,\ldots,\theta_K)=\left(\frac{X_1}{\sum_{k=1}^KX_k},\ldots, \frac{X_K}{\sum_{k=1}^KX_k}\right)\] とすると,\((\theta_1,\ldots,\theta_{K-1})\)は\(\mathcal{D}(\alpha_1,\ldots,\alpha_K)\)に従う.以下では \(K=3\)のときに,乱数\(\theta\)をたくさん生成し,3次元空間に描画した.パラメータを変えて実験してみると良い.set.seed(1234)
N <- 1e3
alpha <- c(1.0, 1.5, 1.8) #Parameter
x <- rgamma(N, alpha[1])
y <- rgamma(N, alpha[2])
z <- rgamma(N, alpha[3])
w <- x + y + z
plot_ly(x=x/w, y=y/w, z=z/w, type="scatter3d", mode="markers", size=0.1)
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
Solution. \(\xi=\sigma^2\)とおく.すると正規分布\(\mathcal{N}(\mu,\xi)\)の確率密度関数は \[\sqrt{\frac{1}{2\pi\xi}}\exp\left(-\frac{1}{2\xi}(x-\mu)^2\right)\] である.だからスコア関数は \[\begin{matrix}S_1\\S_2\end{matrix}:=\begin{pmatrix}\xi^{-1}(x-\mu)\\-\frac{1}{2\xi}+\xi^{-2}\frac{(x-\mu)^2}{2}\end{pmatrix}\] となる.いっぽう,例1.14で紹介したように, \[\mathbb{E}[x-\mu]=\mathbb{E}[(x-\mu)^3]=0,\ \mathbb{E}[(x-\mu)^2]=\xi,\ \mathbb{E}[(x-\mu)^4]=3\xi^2\] となる.したがって,フィッシャー情報行列は \[\mathbb{E}\left[\begin{pmatrix}S_1^2&S_1S_2\\S_2S_1&S_2^2\end{pmatrix}\right]=\begin{pmatrix}\xi^{-1}&0\\0&\xi^{-2}/2\end{pmatrix}\] である.この行列の行列式は\(\xi^{-3}/2\)となるから,パラメータを\((\mu,\xi)\)としたときのジェフリーズ事前分布は \[2^{-1/2}\xi^{-3/2}\] を密度関数として持つ.これでひとつめの直接的な方法でのジェフリーズ事前分布を導出できた.
いっぽう,例1.14では\(\tau=\xi^{-1}\)としてパラメータを\((\mu,\xi)\)としたときのジェフリーズ事前分布が \[\tau^{-1/2}\] を密度関数として持つことがわかっている.また \[(\mu,\tau)\mapsto (\mu,\xi)\] という変数変換によるヤコビ行列式は \[\left|\det\begin{pmatrix}1&0\\0&-\xi^{-2}\end{pmatrix}\right|=\xi^{-2}\] である.よってジェフリーズ事前分布の変数変換による不変性から,パラメータを\((\mu,\xi)\)としたときのジェフリーズ事前分布は \[\xi^{1/2}~\xi^{-2}=\xi^{-3/2}\] となり,たしかにさきほど求めたものと一致する (注: 次のRemark を見よ).Exercise 1.11 \(N\)は正の整数,\(\lambda_1,\lambda_2\)は正の実数とする. \(x_1^1,\ldots, x_N^1|\lambda_1\sim\mathcal{P}(\lambda_1)\), \(x_1^2,\ldots, x_N^2|\lambda_2\sim\mathcal{P}(\lambda_2)\)であり すべて独立とする.二つのモデル \[\mathcal{M}_1:\lambda_1,\lambda_2\sim \mathcal{E}(1), \mathcal{M}_0:\lambda_1=\lambda_2\sim\mathcal{E}(1)\] を考える.二つのモデルへの事前確率は\(0.5\)ずつとする.
(a). モデル\(\mathcal{M}_0\)の事後確率を求めよ.
(b). ベイズ因子\(B_{10}\)を求めよ.
Solution. 例1.15でつかった\(x^N\)や\(x_1^N, x_2^N\)といった記号を再利用したうえで,さらに記号をいくつか用意しよう.各集団の和をそれぞれ \[S_1=\sum_{n=1}^Nx_n^1,\ S_2=\sum_{n=1}^Nx_n^2\] とし,全体の和を \[S=S_1+S_2\] とする.また,各集団の積をそれぞれ \[T_1=x_1^1\cdots x_N^1,\ T_2=x_1^2\cdots x_N^2\] とし,全体の積を \[T=T_1+T_2\] としよう.すると,モデル\(\mathcal{M}_1\)のもとの周辺事後確率密度関数は \[p(x^N|1)=\int_0^\infty p(x_1^N|\lambda_1)p(\lambda_1)\mathrm{d}\lambda_1\times \int_0^\infty p(x_2^N|\lambda_2)p(\lambda_2)\mathrm{d}\lambda_2\] となる.うえの2つの積分の積のうち,最初の一つは \[\int_0^\infty\frac{\lambda_1^{S_1}}{T_1}e^{-N\lambda_1}~e^{-\lambda_1}\mathrm{d}\lambda_1=\frac{\Gamma(1+S_1)}{T_1(N+1)^{1+S_1}}\] と計算できる.2つ目の積についても同じ計算をすれば \[\int_0^\infty\frac{\lambda_2^{S_2}}{T_2}e^{-N\lambda_2}~e^{-\lambda_2}\mathrm{d}\lambda_2=\frac{\Gamma(1+S_2)}{T_2(N+1)^{1+S_2}}\] となることがわかる.よってモデル\(\mathcal{M}_1\)のもとの周辺事後確率密度関数は \[p(x^N|1)=\frac{\Gamma(1+S_1)\Gamma(1+S_2)}{T(N+1)^{2+S}}\] となる.まったく同じ計算で,モデル\(\mathcal{M}_0\)の場合は \[p(x^N|0)=\int_0^\infty\frac{\lambda_1^{S}}{T}e^{-2N\lambda}~e^{-\lambda}\mathrm{d}\lambda=\frac{\Gamma(1+S)}{T(2N+1)^{1+S}}\] を得る.したがって \[\frac{p(x^N|0)}{p(x^N|1)}=\frac{(N+1)^{2+S}}{(2N+1)^{1+S}}\binom{S}{S_1}\] となる.よって
(a). 事後確率は \[\frac{(N+1)^{2+S}}{(2N+1)^{1+S}}\binom{S}{S_1}\left(1+\frac{(N+1)^{2+S}}{(2N+1)^{1+S}}\binom{S}{S_1}\right)^{-1}\].
(b). ベイズ因子は \[\left(\frac{(N+1)^{2+S}}{(2N+1)^{1+S}}\binom{S}{S_1}\right)^{-1}.\]