『統計学実践ワークブック』の第３１章の勉強メモです。
ベイズ法
ベイズ法によるパラメータ推測
同時確率密度関数で定義されるパラメトリックモデル$f(\boldsymbol{x}|\boldsymbol{\theta})$と$\theta$に関する事前情報を確率分布で表した事前分布$\pi(\theta)$の組をベイズモデルという。ここで、$\boldsymbol{x}=(x_1,x_2,\cdots,x_n)^\top$はパラメトリックモデルから生成されるデータ、$\theta\in\Theta$はモデルを定めるパラメータ、$\Theta$はパラメータ空間とする。
データ$\boldsymbol{x}$を観測したときの$\boldsymbol{\theta}$に関する事後情報を確率分布で表した事後分布$\pi(\boldsymbol{\theta}|\boldsymbol{x})$はベイズの定理によって
\pi(\theta|x)=\frac{f(x|\theta )\pi(\theta)}{\int_\Theta f(x|\theta)\pi(\theta)d\theta}
と表される。この事後分布から最良のパラメータ$\boldsymbol{\theta}$を求めることをベイズ法という。
事後分布$\pi(\boldsymbol{\theta}|\boldsymbol{x})$の分母$f(x)=\int_\Theta f(x|\theta)\pi(\theta)d\theta$は規格化定数で、周辺尤度(marginal likelihood)と呼ばれ、$\boldsymbol{\theta}$によらない。これより、事後分布$\pi(\boldsymbol{\theta}|\boldsymbol{x})$は
\pi(\theta|x)\propto f(x|\theta)\pi(\theta)
と表すことができる。
例１
ロボット掃除機の普及率$p$に関する推測を行うために、20世帯のロボット掃除機の保有の有無を調査したところ、$x$世帯が保有していると回答があった。
$p$の事前分布には、２年前の保有率が5%であったことからベータ分布(p.41)$Beta(2,20)$を改定した。このとき、$p$の事後分布を求める。
一般に、ベータ分布$Beta(a,b)$の確率密度関数は
f(\theta|a,b)=\frac{\theta^{a-1}(1-\theta)^{b-1}}{B(a,b)}
と表せる。ここで、ベータ関数$B(a,b)=\int_0^1 \theta^{a-1}(1-\theta)^{b-1}d\theta$となる。
ベータ分布$Beta(a,b)$は、期待値は$\frac{a}{a+b}$、最頻値は$\frac{a-1}{a+b-2}$をとる。
$p$の事前分布$\pi(p)$は
\pi(p)=\frac{p(1-p)^{19}}{B(2,20)}
パラメトリックモデルは
f(x|p)={}_{20} C_{x}p^x(1-p)^{20-x}
となる。
事後分布は
\begin{multline}
\begin{split}
\pi(p|x)&=\frac{f(x|p)\pi(p)}{\int_0^1f(x|p)\pi(p)dp}=\frac{{}_{20} C_{x}p^x(1-p)^{20-x}\frac{p(1-p)^{19}}{B(2,20)}}{\int_0^1{}_{20} C_{x}p^x(1-p)^{20-x}\frac{p(1-p)^{19}}{B(2,20)}dp}\\
&=\frac{p^{x+1}(1-p)^{39-x}}{\int_0^1p^{x+1}(1-p)^{39-x}dp}=\frac{p^{x+1}(1-p)^{39-x}}{B(x+2,40-x)}\\
&=Beta(x+2,40-x)
\end{split}
\end{multline}
データからパラメータ推測
例１では、事後分布が$Beta(2+x,40-x)$と求められた。
これより、例えば、$x=0$だった場合、事後分布は$Beta(2,40)$となり、普及率$p$の期待値は$\frac{2}{2+40}\fallingdotseq 0.048$、最頻値は$\frac{2-1}{2+40 -2}=0.025 $となる。
$x=5$だった場合、事後分布は$Beta(7,35)$となり、普及率$p$の期待値は$\frac{7}{7+35}\fallingdotseq 0.167$、最頻値は$\frac{7-1}{7+35-2}=0.15$となる。
ベイズ法による点推定
ベイズ法の点推定量としてはよくベイズ推定量とMAP推定量が用いられる。
ベイズ推定量
$\theta$の事後分布に対する期待値
\hat{\theta}_B=\int_{\Theta}\theta\pi(\theta|x)d\theta
を$\theta$のベイズ推定量という。
点推定量$\hat{\theta}$の二乗損失$(\hat{\theta}-\theta)^2$の事後分布に対する期待値
\int_{\Theta}(\hat{\theta}-\theta)^2\pi(\theta|x)d\theta
を最小にする推定量がベイズ推定量となる。
証明
\begin{multline}
\begin{split}
\int_{\Theta}(\hat{\theta}-\theta)^2\pi(\theta|x)d\theta&=\hat{\theta}^2\int_{\Theta}\pi(\theta|x)d\theta-2\hat{\theta}\int_{\Theta}\theta\pi(\theta|x)d\theta+\int_{\Theta}\theta^2\pi(\theta|x)d\theta\\
&=\int_{\Theta}\pi(\theta|x)d\theta\left\{\hat{\theta}^2-2\hat{\theta}\frac{\int_{\Theta}\theta\pi(\theta|x)d\theta}{\int_{\Theta}\pi(\theta|x)d\theta}\right\}+\cdots\\
&=\int_{\Theta}\pi(\theta|x)d\theta\left\{\hat{\theta}-\frac{\int_{\Theta}\theta\pi(\theta|x)d\theta}{\int_{\Theta}\pi(\theta|x)d\theta}\right\}^2+\cdots\\
&=\int_{\Theta}\pi(\theta|x)d\theta\left\{\hat{\theta}-\int_{\Theta}\theta\pi(\theta|x)d\theta\right\}^2+\cdots\\
\end{split}
\end{multline}
ここで、
\begin{multline}
\begin{split}
\int_{\Theta}\pi(\theta|x)d\theta&=\int_{\Theta}\frac{f(x|\theta)\pi(\theta)}{\int_{\Theta}f(x|\theta)\pi(\theta)d\theta}d\theta\\
&=\frac{\int_{\Theta}f(x|\theta)\pi(\theta)d\theta}{\int_{\Theta}f(x|\theta)\pi(\theta)d\theta}=1
\end{split}
\end{multline}
以上より、$\hat{\theta}=\int_{\Theta}\theta\pi(\theta|x)d\theta$で最小値となる。
MAP推定量
$\theta$の事後分布の最頻値
\hat{\theta}_{MAP}=\arg_{\theta\in\Theta}\max \pi(\theta|x)
を最大事後確率推定量という。実用的なベイズモデルでは、周辺尤度を解析的に計算することが困難であることが多いが、MAP推定量は周辺尤度を計算しなくてもよいという利点がある。
また、MAP推定量は
\hat{\theta}_{MAP}=\arg_{\theta\in\Theta}\max\{\log f(x|\theta)+\log\pi(\theta)\}
と表すことができる。
ベイズ法による区間推定
ベイズ法において区間推定する場合は、信用区間が用いられる。$\theta$の$100(1-\alpha)$%信用区間$I=(l,u)$は、事後分布を用いて
\int_l^u\pi(\theta|x)d\theta=1-\alpha
を満たす区間として定義される。
良く取られる信用区間として、最大事後密度区間$I_{HPD}$がある。これは、
\begin{multline}
\begin{split}
I_{HPD}&=\{\theta|\pi(\theta|x)\geq c(\alpha)\}\\
\textrm{s.t.}&\int_{I_{HPD}}\pi(\theta|x)d\theta=1-\alpha
\end{split}
\end{multline}
を満たす区間。
事前分布
ベイズ法において事前分布をどのように選ぶかは重要な問題である。
ベータ二項モデル
事前分布にベータ分布を仮定したモデル。
ベータ二項モデルでは事後分布もベータ分布になる。
事前分布と事後分布が同じ確率分布になるような事前分布を共役事前分布という。
ガンマ・ポアソンモデル
$\boldsymbol{x}$をポアソン分布$Po(\lambda)$からの独立なデータとする。このとき、ポアソン分布のモデル$f(\boldsymbol{x}|\lambda)$は
f(x|\lambda)=\Pi_{i=1}^{n}\frac{\lambda^{x_i}e^{-\lambda}}{x_i!}
となる。ここで、$\lambda$にガンマ分布$Ga(\alpha,1/\beta)$を事前分布と仮定する。
ガンマ分布$Ga(\alpha,1/\beta)$の確率密度関数は
\pi(\lambda)=\frac{\beta^{\alpha}\lambda^{\alpha-1}e^{-\beta\lambda}}{\Gamma(\alpha)}
となる。ここで、ガンマ関数
\Gamma(\alpha)=\int_0^\infty t^{\alpha-1}e^{-t}dt
である。事後分布は
\begin{multline}
\begin{split}
\pi(\lambda|x)&=\frac{\left(\Pi_{i=1}^{n}\frac{\lambda^{x_i}e^{-\lambda}}{x_i!}\right)\frac{\beta^{\alpha}\lambda^{\alpha-1}e^{-\beta\lambda}}{\Gamma(\alpha)}}{\int_0^\infty\left(\Pi_{i=1}^{n}\frac{\lambda^{x_i}e^{-\lambda}}{x_i!}\right)\frac{\beta^{\alpha}\lambda^{\alpha-1}e^{-\beta\lambda}}{\Gamma(\alpha)}d\lambda}\\
&=\frac{\lambda^{\sum_{i=1}^n x_i+\alpha-1}e^{-(\beta+n)\lambda}}{\int_0^\infty\lambda^{\sum_{i=1}^n x_i+\alpha-1}e^{-(\beta+n)\lambda}d\lambda}\\
&=\frac{(\beta+n)^{\sum_{i=1}^n x_i+\alpha}\lambda^{\sum_{i=1}^n x_i+\alpha-1}e^{-(\beta+n)\lambda}}{\Gamma(\sum_{i=1}^n x_i+\alpha)}
\end{split}
\end{multline}
となり、$Ga(\sum_{i=1}^n x_i+\alpha,1/(\beta+n))$である。このモデルは共役事前分布であることがわかる。
ガウス・ガウスモデル
ガウス分布のパラメトリックモデルを
f(x|\theta)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\theta)^2}{2\sigma^2}\right)
とする。事前分布をガウス分布で与える。
\pi(\theta)=\frac{1}{\sqrt{2\pi\sigma_0^2}}\exp\left(-\frac{(\theta-\mu_0)^2}{2\sigma_0^2}\right)
事後分布は
\begin{multline}
\begin{split}
\pi(\theta|x)&=\frac{\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\theta)^2}{2\sigma^2}\right)\frac{1}{\sqrt{2\pi\sigma_0^2}}\exp\left(-\frac{(\theta-\mu_0)^2}{2\sigma_0^2}\right)}{\int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\theta)^2}{2\sigma^2}\right)\frac{1}{\sqrt{2\pi\sigma_0^2}}\exp\left(-\frac{(\theta-\mu_0)^2}{2\sigma_0^2}\right)d\theta}
\end{split}
\end{multline}
例３
ガンマ・ポアソンモデルの事前分布を$Ga(2,1)$とすると、事後分布は
Ga(2+\sum_{i=1}^5x_i,1/(1+6))=Ga(22,1/6)
となる。
ベイズ推定量は期待値となるので、$\hat{\lambda}_B=\frac{22}{6}=\frac{11}{3}$となる。
MAP推定量は最頻値なので、$\hat{\lambda}_B=\frac{22-1}{6}=\frac{7}{2}$
ベイズ予測
観測されたデータ$\boldsymbol{x}=(x_1,x_2,\cdots,x_n)^\top$を用いて、将来観測されるデータ$x_{n+1}$の分布を予測する。$\boldsymbol{x}$、$\boldsymbol{\theta}$を所与としたときの$x_{n+1}$の条件付き確率密度関数を$f(x_{n+1}|x,\theta)$とする。観測データは$\boldsymbol{\theta}$のもとで独立と仮定すると、$f(x_{n+1}|x,\theta)=f(x_{n+1}|\theta)$となる。
ベイズ法では
f(x_{n+1}|x)=\int_{\theta\in\Theta}f(x_{n+1}|x,\theta)\pi(\theta|x)d\theta
を用いて$x_{n+1}$の予測を行う。これをベイズ予測分布という。
マルコフ連鎖モンテカルロ法
求める確率分布をマルコフ連鎖として作成し、乱数の生成により確率分布のサンプリングを求める方法をマルコフ連鎖モンテカルロ法という。
過去の値によらず、現在の値だけで未来の値の確率が決まる確率過程のことをマルコフ過程という。このマルコフ過程のうち、取りうる状態が離散的、特に時間が離散的なものをマルコフ連鎖という。マルコフ連鎖の状態遷移確率は
Pr(X_{n+1}|X_{n},X_{n-1},\cdots,X_0)=Pr(X_{n+1}|X_{n})
となる。
マルコフ連鎖の重要な性質として、定常分布を持つことである。定常分布とは、過程が時間に依存しなくなる状態である。マルコフ連鎖モンテカルロ法はこの定常分布に至る性質も用いて確率分布を求める。
MHアルゴリズム
メトロポリス・ヘイスティング法はマルコフ連鎖の定常分布からマルコフ過程を設計する。
状態が$\boldsymbol{x}$から$\boldsymbol{x}'$に遷移することを考える。状態$\boldsymbol{x}$の状態確率を$P(x)$とすると、定常状態で詳細つり合い条件が成り立つので、
P(x)P(x'|x)=P(x')P(x|x')
となり、
\frac{P(x'|x)}{P(x|x')}=\frac{P(x')}{P(x)}
とできる。
遷移確率を$P(x'|x)=Q(x'|x)A(x',x)$とする。ここで、提案分布$Q(x'|x)$は$\boldsymbol{x}$が与えられたときの状態$\boldsymbol{x}'$を提案する条件付き確率、採択確率$A(x',x)$は$\boldsymbol{x}$が与えられたときの状態$\boldsymbol{x}'$を採択する条件付確率。代入して整理すると、
\frac{A(x',x)}{A(x,x')}=\frac{P(x')}{P(x)}\frac{Q(x|x')}{Q(x'|x)}
となる。この式を満たす採択確率は
A(x',x)=\min\left(1,\frac{P(x')}{P(x)}\frac{Q(x|x')}{Q(x'|x)}\right)
を採用することができる。
ベイズ法では周辺尤度$\int_\Theta f(x|\theta)\pi(\theta)d\theta$の計算困難な場合があるが、MH法では事後分布に依存しているのは採択確率のみであり、
\begin{multline}
\begin{split}
A(x',x)&\to \frac{\pi(\theta^{t+1}|x)}{\pi(\theta^{t}|x)}=\frac{f(x|\theta^{t+1})\pi(\theta^{t+1})}{f(x|\theta^{t})\pi(\theta^{t})}
\end{split}
\end{multline}
とできるので、周辺尤度に依存しないことがわかる。
ギブス・サンプリング
状態$\boldsymbol{x}$から$x_i$だけが変化して状態$\boldsymbol{x}'$になることを考える。このとき、遷移確率は$P(x'|x)=P(x'|x_1,\cdots,x_{i-1},x_{i+1},\cdots,x_n)$とする。詳細つり合い条件について、
\begin{multline}
\begin{split}
P(x)P(x'|x)&=P(x)P(x'|x_1,\cdots,x_{i-1},x_{i+1},\cdots,x_n)\\
&=P(x)\frac{P(x')}{\sum_{x_i}P(x_1,\cdots,x_{i-1},x_i,x_{i+1},\cdots,x_n)}\\
&=P(x')\frac{P(x)}{\sum_{x_i}P(x_1,\cdots,x_{i-1},x_i,x_{i+1},\cdots,x_n)}\\
&=P(x')P(x|x')
\end{split}
\end{multline}
成り立つことがわかる。これより、この遷移確率を用いて、マルコフ過程を設計できることがわかる。
例題
問３１．１
[1]
パラメトリックモデルは二項分布となり、$f(X|\theta)=nC_X\theta^X(1-\theta)^{n-X} $
周辺分布は、$\int_0^1 f(x|\theta)U(0,1)d\theta$
事後分布は
\begin{multline}
\begin{split}
\pi(\theta |3)&=\frac{{}_{10}C_3\theta^3(1-\theta)^{7}U(0,1)}{\int_0^1{}_{10}C_3\theta^3(1-\theta)^{7}U(0,1)d\theta}\\
&=\frac{{}_{10}C_3\theta^3(1-\theta)^{7}}{\int_0^1{}_{10}C_3\theta^3(1-\theta)^{7}d\theta}\\
&=\frac{{}_{10}C_3\theta^3(1-\theta)^{7}}{B(4,8)}=Beta(4,8)
\end{split}
\end{multline}
ベイズ推定量は期待値なので、$\frac{4}{4+8}=\frac{1}{3}$
[2]
事前分布がベータ分布なので、ベータ二項モデルとなる。よって事後分布は$Beta(1+3,4+7)=Beta(4,11)$
MAP推定量は最頻値であるので、$\frac{4-1}{4+11-2}=\frac{3}{13}$
問３１．２
[1]
ガンマ・ポアソン分布なので、事後分布は
Ga(\alpha +\sum_i x_i,\frac{1}{\beta + n})=Ga(17,\frac{1}{8})
となる。確率密度関数は
\pi(\lambda |x)=\frac{\beta^\alpha \lambda^{\alpha-1}e^{-\beta\lambda}}{\Gamma(\alpha)}=\frac{8^{17}\lambda^{16}e^{-8\lambda}}{\Gamma(17)}
となる。
[2]
\begin{multline}
\begin{split}
f(x_6|\lambda)&\equiv f(x_6|x,\lambda)=\frac{f(x_6,x,\lambda)}{f(x,\lambda)}=\frac{f(x_6,x|\lambda)f(\lambda)}{f(x|\lambda)f(\lambda)}\\
&=\frac{f(x_6,x|\lambda)}{f(x|\lambda)}=\frac{\Pi_{i=1}^{6}\frac{\lambda^{x_i}e^{-\lambda}}{x_i!}}{\Pi_{i=1}^{5}\frac{\lambda^{x_i}e^{-\lambda}}{x_i!}}=\frac{\lambda^{x_6}e^{-\lambda}}{x_6!}
\end{split}
\end{multline}
これより、
\begin{multline}
\begin{split}
f(x_6|x)&=\int_0^{\infty}f(x_6|\lambda)\pi(\lambda|x)d\lambda=\int_0^{\infty}\frac{\lambda^{x_6}e^{-\lambda}}{x_6!}\frac{8^{17}\lambda^{16}e^{-8\lambda}}{\Gamma(17)}d\lambda\\
&=\frac{8^{17}\int_0^{\infty}\lambda^{16+x_6}e^{-9\lambda}d\lambda}{x_6!\Gamma(17)}=\frac{8^{17}\Gamma(17+x_6)}{9^{17+x_6}x_6!\Gamma(17)}
\end{split}
\end{multline}
問３１．３
[1]
$\pi(y),\pi(x)$が目標分布になる必要があるので、
C=\min\left\{1,\frac{\frac{1}{8}\phi(y+5)+\frac{3}{4}\phi(y)+\frac{1}{8}\phi(y-5)}{\frac{1}{8}\phi(x^{(t)}+5)+\frac{3}{4}\phi(x^{(t)})+\frac{1}{8}\phi(x^{(t)}-5)}\right\}
[2]
(ア)
$\alpha=1$
理由は、$x$の変動が(エ)の次に緩いため。
(イ)
$\alpha=10$
理由は、$x$の変動が(ア)に次に緩いため。
(ウ)
$\alpha=100$
理由は、$x$の変動が他３つと比べて明らかに少なく、これは採択確率があまりにも低かったと考えられる。採択確率が一番低くなるのは$\alpha$の値が一番大きいものとなる。
(エ)
$\alpha=0.1$
理由は、$x$の変動幅が小さいため。
[3]
初めの方に生成した$x$は初期値の影響を強く受けているため。
問３１．４
[1]
正規分布の条件付き分布はp.44より、
期待値は$\mu_{x|y}=\mu_x+\frac{\sigma_{yx}}{\sigma_y^2}(x_y-\mu_y)=\frac{1}{2}(y^{(t)}-2)$
分散は$\sigma_{x|y}^2=\sigma_x^2(1-\rho^2)=\frac{3}{4}$
となるので、$N(\frac{1}{2}y^{(t)}-1,\frac{3}{4})$
[2]
期待値は$\mu_{y|x}=\mu_y+\frac{\sigma_{xy}}{\sigma_x^2}(x_x-\mu_x)=2+\frac{1}{2}x^{(t)}$
分散は$\sigma_{y|x}^2=\sigma_y^2(1-\rho^2)=\frac{3}{4}$
となるので、$N(\frac{1}{2}x^{(t)}+2,\frac{3}{4})$